Mixing in the NETmix Reactor

NETmix is a static mixing reactor composed of a network of mixing chambers interconnected by channels. The repetitive mixing pattern inside the reactor enables the use of reduced geometries to represent the NETmix network, such as the ExtendedNUB model, used in this work. Mixing in NETmix is based on the impingement of jets, issuing from channels. Inside the chambers, the jets are engulfed by dynamic vortices which can be quantified using Lagrangian techniques. Batch Lagrangian Mixing Simulation (BLMS) is based on successive injections of particles to measure the fraction of the fluids at the outlet of the mixing chambers. The distribution of the outlet fraction of particles indicates that it is possible to have nearly perfect mixing inside the NETmix chambers, depending on the dimensions of the channels and chambers. The NETmix design is here optimized in relation to the chamber diameter to channel width ratio, D / d . Results from BLMS show that best performance in NETmix occurs for 6.65 ≤ D / d ≤ 6.85 .


INTRODUCTION
NETmix is a static mixer patented in 2005 (Lopes et al., 2005). The design of the reactor is based on a network of mixing chambers interconnected by transport channels (Lopes et al., 2005). The mixing chambers are elements with a circular footprint which can be cylinders or spheres, the channels that connect the chambers are cylindrical or prismatic. The centre of these elements lays in a 2D plane.
The fluids enter the NETmix reactor through separated inlets placed at the bottom row of the network and leave through a set of outlets at the top row. Figure 1 shows the NETmix network and flow patterns from a tracer CFD simulation for a chaotic flow regime . The NETmix reactor enables the control of macro and micro-mixing which can be manipulated to improve the performance of the reactor (Laranjeira et al., 2009). Moreover, the periodic nature of the NETmix design enables a direct laboratory to industrial scale by numbering-up the network (Schenk et al., 2004;Saber et al., 2010).
The first NETmix device was designed in 2005 based on a network of spherical chambers connected by cylindrical channels (Laranjeira 2005). In 2008, an industrial prototype was built at Fluidinova S.A. (fluidinova.com), for the continuous large-scale production of synthetic nano-hydroxyapatite (Lopes et al., 2006;Lopes et al., 2007;Silva et al., 2008;Gomes et al., 2009). CFD studies have shown that the mixing mechanisms in NETmix are well described from 2D simulations (Gomes 2011). In this way, other NETmix units were developed, using cylindrical chambers and prismatic channels (Gomes 2011). Nowadays the NETmix network is etched in plates, making it particularly suited for coupling other devices such as heating and radiation plates. The diameter of the chambers and channels was also varied between different NETmix devices which enables to switch from a micro to a meso device. In this way it is possible to have better control of heat and mass transfer (Costa, 2017) which has been widening the applications of this technology. Since 2016, a photocatalytic NETmix was designed (Lima et al., 2016) which is used to perform photocatalysis studies such as kinetic modelling (Marinho et al., 2017;Santos et al., 2019;Santos et al., 2021) and ozonation (Filho et al., 2019). In 2017, a NETmix device coupled to heat exchangers was developed for the continuous production of CO 2 hydrates Lopes et al., 2018;Lopes et al., 2019). The NETmix reactor has also been used for microencapsulation (Moreira et al., 2020) and Pickering emulsions (Ribeiro et al., 2021). Figure 1 shows the mixing patterns inside the NETmix reactor which are repetitive along each row of mixing chambers. In this way, it is possible to study how mixing occurs inside the NETmix chambers, i.e., micromixing, through the simulation of a reduced domain of the network Costa et al., 2017), replacing the full network by periodic boundary conditions (PBCs). Previous reduced physical domains of NETmix include the NETmix Unit Block (NUB) Costa et al., 2017) which will be revised in this work.
The simulation of the NETmix mixing dynamics in each chamber is affected by the boundary conditions. In this way, the reduced geometry needs to ensure that the model does not affect the mixing dynamics. This was thoroughly covered by Torres (2017) from spectral analysis. In this work, the reduced geometric model, ExtendedNUB, will be presented and validated from second-moment statistics. Turbulence intensity is related to the magnitude of the velocity fluctuation of the fluid associated with an enhancement of mixing in laminar chaotic or turbulent flows (Pope 2000).
The main operational parameter in NETmix for single phase flows is the Reynolds number of the channels, Re. The Reynolds number in the NETmix reactor is defined at the inlet channels as Re ρ2dυinj μ , in which d is the width of the channels, υ inj , is the velocity of the fluid at the inlets and ρ and μ are the density and viscosity of the fluid, respectively. This parameter marks the transition from segregated laminar, where fluids flow segregated on their side of the chambers, to chaotic flow regimes (Laranjeira et al., 2011). Chaotic flow regimes occur above a critical Reynolds number of ca. 120 (depending on the network topology), when the impingement point of the inlet jets inside the mixing chambers starts oscillating and onsets a dynamic behaviour of the vortices (Laranjeira et al., 2009;Laranjeira et al., 2011). The engulfment of the vortices that are formed specially at the top and bottom of the chambers is responsible for mixing. In fact, tracer experiments have shown that above a certain row of chambers it is not possible to distinguish the flow patterns in the chambers due to the occurrence of a quasi-homogenization (Laranjeira et al., 2011).
Mixing at larger scales, i.e., macromixing, can be assessed in the NETmix reactor through the visualization of the overall patterns as shown in Figure 1, using tracer simulations which measure the relation between the composition of the fluids at the inlet and outlet. In each NETmix chamber, if the composition of the fluid at the outlet is diluted by a factor of two in relation to the composition of the fluid at the inlet of that chamber, perfect mixing occurs. The overall mixing capacity of the NETmix reactor is then 2 n , in which n represents the number of chambers. The composition of the fluid that leaves the NETmix network at the top row is diluted 2 n times in relation to the composition of the fluid that entered in the NETmix network at the bottom row. A way to evaluate the degree of mixing in the NETmix reactor is to use particles to monitor, at the outlet, the distribution of the fraction of the fluids being mixed. Considering one chamber, if the fraction of particles injected in the left inlet that exit the chamber through the right outlet is 0.5, the dilution factor is 2, and thus, mixing is perfect. In this work, the algorithm Batch Lagrangian Mixing Simulation (BLMS) was developed and implemented to access the mixing of two fluids in each mixing chamber.
The mixing capacity of reactors is highly dependent on the device topology. Since its conception, it was noticeable that the configuration of the network has a major impact on the NETmix performance (Laranjeira 2005;Gomes 2011). One of the first works was focused on the effect of the ratio between the volumes of the chambers and channels on micromixing. Laranjeira et al. (2009) used a network model where the chambers were modelled as perfectly mixing zones and the channels as plug flow segregation zones to study the NETmix network design. Results of the ratio between the volume of the channels and the whole network volume have led to the present NETmix configuration where each chamber is connected to four channels, forming an angle of 45°with the chamber axis, i.e., 90°angle between neighbouring channels. Previous works have also shown that macromixing and micromixing in NETmix depend on the number of chamber rows in the flow direction and on the injection scheme of the reactants (Laranjeira et al., 2009;Laranjeira et al., 2011).
Since the topology of the reactor exerts a major influence in the mixing dynamics of the NETmix chambers, this work is dedicated to the evaluation of mixing inside each chamber and, from this study, improve the design of the reactor. The NETmix topology will be optimised in terms of the ratio chamber diameter/channel width, D/d, to achieve a dilution factor as close as possible to two inside each mixing chamber. For that, BLMS will be applied in the ExtendedNUB model to analyse the influence that the ratio D/d has in mixing.

NETMIX NETWORK
The NETmix model was studied and validated in previous works numerically and experimentally (Laranjeira et al., 2009;Laranjeira et al., 2011). The NETmix model uses a network composed of chambers interconnected by channels whose arrangement has a crucial impact on mixing.
The NETmix network (Figure 2 left) is composed of cylindrical chambers with a diameter, D, interconnected by parallelepipedic channels with width, d, and length, l, at a 45°a ngle (Lopes et al., 2005;Laranjeira 2005;Gomes 2011). This network can be obtained by the repetition of the NETmix unit cell shown in Figure 2 (right), where ω is the reactor's depth. Experimental and computational results have shown that since mixing mechanisms are mainly 2D in the xOy plane (Gomes 2011), it is possible to accurately represent the mixing dynamics using a 2D model. The mixing patterns from 2D CFD simulations are the same as those observed from tracer experiments, as will be later addressed in this work. The validity of the 2D modelling of the flow dynamics is restricted to the range of depth values tried so far, 2 ≤ ω/d ≤ 3.9, which is close to a similar analysis on the validity of 2D assumption in T-jets reactors (Sultan et al., 2019). For the range of depth values that has been used for the design of NETmix reactors, the 2D physics provide an accurate description of the flow dynamics.

ExtendedNUB Geometry
Previous works have shown that the geometry and the flow in NETmix obey to a repetitive pattern (visible in Figure 1) which enables to simulate only a subdomain of the full network to study the NETmix reactor Costa et al., 2017). Fonte (2013) and Fonte et al. (2013) developed a 3D geometry that consisted of three full chambers interconnected to 6 half chambers, the NETmix Unit Block (NUB) shown in Figure 3 (left). The half chambers in each row of the network are connected to each other by a periodic boundary condition as shown in Figure 3 (left). Torres (2017) studied the possibility of simulating reduced NETmix geometries in 2D NUB models. One of the models, the ExtendedNUB, consisted in using a repetition of two NUBs (Figure 3 right). The ExtendedNUB geometry has two bottom inlets and two top outlets and comprises four full chambers (C2, C4, C6 and C8 in Figure 3 right) and five pairs of half chambers. This network has a total of nine rows of channels issuing from the chambers (nine left channels, l1 to l9 and nine right channels, r1 to r9). This geometry will be validated further in this work.

CFD Model and Boundary Conditions
The 2D domain was discretized with a symmetrical mesh in relation to the x and y axis, using ANSYS ICEM. A quarter of chamber is shown in Figure 4. The cell elements have an average edge length, Δx 1.7 × 10 −5 m.
The commercial finite volume CFD software ANSYS Fluent ™ was used to solve the mass and momentum conservation equations for the ExtendedNUB model, respectively defined as where υ is the velocity vector, p is the pressure, ρ and μ are the density and viscosity of the fluid, respectively. The SIMPLEC scheme was used with PRESTO! for the pressure-velocity coupling and Third-Order MUSCL for the spatial discretization. The convergence criterion for the maximum residuals of continuity and momentum is 10 −5 .
At the inlets, a constant velocity is imposed, υ inj , which is obtained from the Reynolds number, according to where D h 2d is the hydraulic diameter for the 2D geometry, ρ 998.2 kg/m 3 and μ 1.0 mPa.s. At the outlets, a constant and uniform pressure profile is used. As extensively analysed by Torres (2017), it is possible to use reduced geometries with only one column to represent the NETmix network if periodic boundaries are used at the half chambers, as shown in Figure 3. A periodic boundary condition is then set at each pair of half lateral chambers. Velocity components x and y are set to zero at the walls of NETmix.
To perform the transient CFD simulation, the initial solution was obtained by setting all variables to zero.

2D ExtendedNUB Geometry Validation
To validate the 2D ExtendedNUB model, the dynamic flow field was simulated for a geometry with a chamber diameter, D 6.65 × 10 −3 m, channel length, l 2.0 × 10 −3 m and channel width, d 1.0 × 10 −3 m. Several Reynolds numbers were used from 100 to 340. The time step was calculated according to the Courant number, C υ inj Δt Δx 1. As an example, for Re 300, Δt 1.14 × 10 −4 s .
The residence time in one NETmix chamber is given by τ c πD 2 8υinjd 0.12s for the geometry dimensions above mentioned. The residence time of the entire ExtendedNUB is then τ 1.1s. In CFD simulations, the oscillating behaviour referred in the Introduction and shown in Figure 1 takes about 9 τ (10 s  flowtime) to be reached as will be shown further in this work. The simulations need, then, to be run for more than 10 s in order to generate dynamic data that represents the characteristic oscillating behaviour of the NETmix reactor. In this case, the simulation was run during 34 s flowtime and the results were recorded for analysis only from 17 to 34 s. As an example, the contour plot of the z-component of the vorticity for the 2D ExtendedNUB simulation with Re 300 is shown in Figure 5 (left), where the oscillating behaviour inside the NETmix chambers along with the vortex formation on the top and bottom of each chamber is observed. The vortices engulf the jets issuing from the channels and make the angle of contact between the two jets oscillate periodically from left to right. The overall mixing pattern is consistent with the one shown in the tracer simulation of Figure 1, which is also observed from experiments (Laranjeira et al., 2011). In contrast, if a Reynolds number below the critical value is chosen, the oscillating behaviour inside the chambers does not occur and there is no formation of vortices that dynamically evolve in the chamber, as shown in Figure 5 (right) for Re 100. When boundary conditions are imposed to a reduced model, it is necessary to ensure that the flow dynamics in the mixing chambers are not affected. As described in CFD Model and Boundary Conditions, at the inlet, there is a constant velocity and at the outlet, the constant and uniform pressure profile imposes that the velocity vectors are normal to the outlet surface, i.e., parallel outlet flow, which causes damping of the flow dynamics. In this way, turbulence intensity, I υ x , studies were performed in the centre of all full chambers C2, C4, C6 and C8, in the direction of the main flow oscillations.
The turbulence intensity is defined by where υ inj is the surface velocity at the inlet, v xi is the x component of the velocity at several flowtimes, υ x is the average of υ xi and t i is the flowtime.
The onset of chaotic mixing mechanisms in NETmix occurs above a critical Reynolds (ca. 120) where the steady flow evolves to fully developed laminar chaotic until Re 300. For Re > 300, the head loss increases without significant improvement in the mixing performance . To perform the turbulence intensity study, different Reynolds numbers were chosen from 140 (just above critical) to 300, by varying the inlet velocity. The time-step size was also recalculated to maintain the Courant number at 1.
The inlet velocity and time-step size of each case are: υ inj 7.02 × 10 −2 m/s, Δt 2.44 × 10 −4 s for Re 140; υ inj 9.03 × 10 −2 m/s, Δt 1.90 × 10 −4 s for Re 180; υ in j 1.10 × 10 −1 m/s, Δt 1.55 × 10 −4 s for Re 220; υ inj 1.31 × 10 −2 m/s, Δt 1.31 × 10 −4 s for Re 260; and υ inj 1.51 × 10 −2 m/s, Δt 1.14 × 10 −4 s for Re 300. Figure 6 shows that the turbulent intensity is larger than zero for all Reynolds numbers. The flow is then above the critical Re which onsets dynamic flows. For lower Re, between 140 and 220, the turbulence intensity is nearly constant. As Re increases, the turbulence intensity is higher, as expected. At Re 260, an increase is observed for chamber C8. At Re 300 it is visible that the turbulence intensity is higher in chamber C6 and then decreases in chamber C8. This is explained by the effect of the outlet boundary condition, which dampens the oscillations. In this way, since this effect is not noticeable in chamber C6, the ExtendedNUB model can be used to represent the NETmix reactor for Re 300 if the studies are performed in chamber C6.

NETMIX NETWORK TOPOLOGY OPTIMISATION
In the NETmix reactor, the degree of mixing that occurs in one chamber can be analysed through the outlet distribution of the fluids being mixed, i.e., the fraction of fluid injected at the left inlet that goes out through the right outlet and vice-versa. In this work, Lagrangian methods will be used to track the injected fluids from the inlet to the outlet of one chamber. This algorithm will then be used to optimise the NETmix network topology in terms of the ratio chamber diameter/channel width, D/d, which is expected to have a major influence on the mixing capacity of the NETmix reactor.

Batch Lagrangian Mixing Simulation (BLMS)
Batch Lagrangian Mixing Simulations (BLMS) was based on the computational implementation of a previous algorithm, Lagrangian Mixing Simulation (LMS) presented by Matos et al. (2018). In LMS, a continuous injection of particles is used to track the filling front of fluids being mixed, according to in which X i,j is the matrix of the positions of the particles, Δt is the time step (which is equal to the time step size used for the CFD simulation that provided the flow-field) and υ is a matrix with the velocities of the particles, i represents the time (iteration) and j the particle. The time difference between the velocity fields used for the tracking is 10Δt. In Eq. 5, the time interval is then set accordingly. The position matrix of the particles, X i,j , is growing in LMS since new particles are inserted when two consecutive particles become apart more than 1.0 × 10 −4 m. This step of inserting new particles is used to keep the accuracy of the fluid front tracking since the main objective of LMS is to determine all convective mixing scales.
The degree of mixing in the NETmix reactor can be assessed from macromixing in a unitary cell where the dynamic flow is fully developed. Since this can be done without simulating the full complexity of the interface between the fluids being mixed, independent discrete particles are used. In BLMS the number of particles is constant at each injection, which enables to follow the path of each particle from the inlet to the outlet and associate it to a fraction of each fluid being injected.
A first batch injection is set at flowtime t 0 . Each injection is performed during ni iterations until all the particles have left the mixing chamber. The process is then restarted and a new injection is set at t 0 + 10Δt, in which Δt is the time-step size of the CFD simulation that provided the velocity-fields. In this way, a series of injections are performed with interval 10Δt until the total number of injections, nt, is reached. The injection procedure is represented in Figure 7, in which two series of np particles (the number of particles will be discussed further in this section) are injected, a red series at the right inlet and a blue series at the left inlet. The example of Figure 7 shows BLMS applied at chamber C6 of the 2D ExtendedNUB model with D/d 6.65 for Re 300.
In Figure 7A-D, the velocity field vectors show oscillating behaviour inside the mixing chamber. The oscillations are also visible from the path of the particles, which follows the inlet streams ( Figures 7B,C) and then are engulfed by the dynamic vortices that onset in NETmix above the critical Re ( Figure 7D). In Figure 7D it is possible to see the particles leaving the mixing chamber, which enables to register the number of blue and red particles that exit the mixing chamber through the left or right outlet. The complete BLMS algorithm was developed in MATLAB and considers the following steps: Bach Lagrangian Mixing Simulation (BLMS) Algorithm 1. Specify t 0 , np, Δt, ni and nt; 2. Load the first velocity-field file, i 1, for the first injection, n 1, at t 0 + (n − 1)10Δt; 3. Position the particles over the inlet; 4. Calculate the velocities of the particles with MATLAB function interp2(); 5. Recalculate the position of the particles with the LMS equation (Eq. 5); 6. Register the number of blue and red particles that exit the mixing chamber on the left or right outlet, n l b , n r b , n l r and n r r ; 7. Load the next velocity field, i+1, at t 0 + i10Δt; 8. Repeat steps 4 to 6 until i ni; 9. Load the first velocity-field file, i 1, for the next injection, n n + 1, at t 0 + (n − 1)10Δt; 10. Repeat steps 3 to 8 until n nt. The respective block diagram is shown in Figure 8.
In BLMS, the position of the particles at the inlets was set according to the flow-field, the zones with higher velocity had a larger density of particles. i.e., the distance between the particles is smaller. The velocity profile, υ x i , at the inlet is close to a parabolic profile. The distance between two consecutive particles, d x i , was than calculated in order to have d x i (υ x i − υ x i −1 ) constant. A schematic representation is in Figure 9 and the algorithm for the calculation of the distance between the particles, dx i is presented next.
Particles Distribution Over the Inlet Algorithm 1. Select the coordinates x and y for the first and last particles (P1 and Pn), according to the coordinates of the inlet channel; 2. Draw a line between P1 and Pn; 3. Build vectors x0 and υ0 with x-coordinate and x-velocity of the mesh cells that are crossed by the line in 2;  Frontiers in Chemical Engineering | www.frontiersin.org November 2021 | Volume 3 | Article 771476 7 4. Start with dx i 2.5 × 10 −5 m ; 5. Calculate the area below υ x curve: A i υ x i +υ x i−1 2 dx i ; 6. Calculate the difference between consecutive areas: ΔA |A i − A i−1 |; 7. Reduce dx i by a factor of 10 −7 m and repeat steps 5 and 6 until ΔA ≤ 10 −8 m 2 /s. 8. Repeat steps 4 to 7 until the coordinates of Pn are reached. 9. The number of intervals, i, that are needed to reach Pn is used to calculate the number of particles of each series, np i+1. Figure 10 shows the beginning of the particle tracking after one injection where it is visible that in the tip of the profile the density of particles is larger.

Diameter Ratio Study
The ratio chamber diameter/channel width, D/d, is the critical dimension for mixing in NETmix, when the depth is in the range of depth values where the 2D assumption is valid. BLMS was applied in the 2D ExtendedNUB model for chamber diameters in the range 6.25 × 10 −3 m ≥ D ≥ 7.25 × 10 −3 m; the channel length, l 2.0 × 10 −3 m and the channel width, d 1.0 × 10 −3 m were kept constant. For each D/d, a CFD simulation, similar to the one shown in Figure 5 (left), was performed for Re 300 to provide the velocity field for the injection of particles with BLMS. For each D/d the velocity field data of chamber C6 was used to perform nt 5000 injections, each one during 0.58 s flowtime (ni 510 iterations). Note that, since each NETmix chamber has a residence time, τ 0.12s (for D/d 6.65), the particles are tracked for ca. 3 τ. The first injection was made at t 0 17s flowtime. The number of particles for each series depends on the ratio D/d and varied from 25 to 35.
To evaluate the particles distribution at the outlet, the exit conditional probabilities were computed. P l b is the probability of a blue particle exiting the mixing chamber on the left outlet, P r b is the probability of a blue particle exiting the mixing chamber on the right outlet, P l r is the probability of a red particle exiting the mixing chamber on the left outlet and P r r is the probability of a red particle exiting the mixing chamber on the right outlet. These probabilities were calculated according to in which n l b is the number of blue particles that exit on the left outlet, n r b is the number of blue particles that exit on the right outlet, n l r is the number of red particles that exit on the left outlet and n r r is the number of red particles that exit on the right outlet (see point six of the BLMS algorithm).
The results presented in Figure 11 show that the relation between macromixing in the NETmix chambers and the ratio D/ d is represented by a butterfly shaped graphic centred at D/d 6.75. For all ratios, the plot is symmetrical in relation to the probability line, P 0.5. For both the red and blue particles, the ratios 6.65 and 6.85 are the ones that have the best mixing since the exit conditional probability is very close to 0.5 which indicates that the dilution factor of two is attained. Looking to all the ratios, 6.75 is the one that has a probability closest to 0.5; however, this only happens for the blue particles. As the ratios move from the middle value 6.75 to the edges, 6.25 and 7.25, the probability moves away from 0.5 and the mixing rate decreases. D/d 6.65 is the ratio which has the overall average, i.e., the average for both, the blue and the red sets of particles, closest to 0.5, so this ratio is considered the optimal geometry for unitary cells.

TRACER SIMULATION
In the previous section it was concluded that the ratio D/d 6.65 is the one that has better mixing. The ExtendedNUB geometry with this ratio is then going to be used to perform a tracer simulation. The tracer simulation will be used to assess the performance of the BLMS method to characterize mixing and validate the ExtendedNUB model for the simulation of mixing in NETmix.
The commercial finite volume CFD software ANSYS Fluent ™ was used for the tracer simulation, considering the models described in CFD Model and Boundary Conditions, with the additional species transport equation in which Y i is the mass fraction of each species, i, and J i is the diffusion flux vector, defined as in which D i 10 −9 m 2 /s is the mass diffusion coefficient of species i in the mixture. The spatial discretization of the geometry is also described in CFD Model and Boundary Conditions. The tracer simulation was performed for Re 300 and the parameters of 2D ExtendedNUB Geometry Validation were maintained. Two fluids, Fluid1 and Fluid2, with the same properties of water were used. The reactor was initially filled with Fluid1. Fluid2 was injected through the right inlet when the dynamic simulation was started. Figure 12 shows the evolution of the molar fraction of Fluid2. Figure 12A shows the injection of Fluid2 at the right inlet. In Figure 12B, Fluid2 is only present in the right side of the reactor since the oscillating behaviour of the jets formed in the NETmix channels is still being developed. In less than 2 s, mixing starts to occur and Fluid2 appears in the left side of the reactor ( Figure 12C). The mixing process continues and in less than 6 s, the average of the molar fraction of Fluid2 in the top chambers, C6 and C8 is already close to 0.5, as shown in Figure12F. Figure 13 (left) shows the history of the molar fraction of Fluid2 in the centre of chamber C6 where from approximately 6 s flowtime the molar fraction of Fluid2 tends to 0.5. The first 5 s correspond to the injection of Fluid2, which starts to occupy the reactor already filled with Fluid1. Although it takes just 6 s to have nearly perfect mixing inside the NETmix reactor, the oscillatory behaviour takes about 10 s to completely develop, as shown in Figure 13 (right) where only from 10 s, the amplitude and frequency of the flow oscillations become constant.
In order to understand the mixing process of Fluid1 and Fluid2 inside the reactor, the average of the molar fractions between 10 and 17 s flowtime was evaluated in the centre of each channel at the right, r, and at the left, l (see Figure 3 right). The result is shown in Figure 14 where the fifth row of chambers is close to 0.5. From the fifth row of chambers, the average of the molar fractions continues to approach 0.5, reaching exactly these values in row 8 which are maintained in row 9.
This result shows that the ExtendedNUB geometry is capable of representing the behaviour of a larger NETmix network since perfect mixing is reached in the top chambers. Moreover, nearly perfect mixing is obtained from the fifth row, i.e., before chamber C6. The mixing patterns obtained from tracer simulations accurately described those obtained experimentally in geometries with depths in the range: 2 ≤ ω/d ≤ 3.9. This validates the 2D ExtendedNUB for simulation of flow and mixing in NETmix reactors in the range of depth values that are most typical in industrial applications. Moreover, this study was performed for single-phase mixing without heat exchange. For multiphase flows as well as reactive systems, especially if they include temperature variations, mixing features might change.
In the centre of chamber C6, the average of the molar fraction is nearly the same as the one reported by the BLMS study in 2D ExtendedNUB Geometry Validation. Table 1 shows the comparison of the values for Fluid2 which was injected in the right side, represented in red, with the red particles that were injected in BLMS, P l r and P r r . For Fluid1, which was injected in the left side, represented in blue, the comparison is performed with the blue particles of BLMS, P l b and P r b . The similarity of the values between BLMS and the tracer simulation shows that the injection of particles is a good method to evaluate the performance of this type of reactor.

CONCLUSION
For the first time, the degree of mixing of the NETmix chambers was assessed and the ratio chamber diameter/channel width, D/d, was studied. Mixing simulations show that it is possible to have nearly perfect mixing in the NETmix reactor for geometries in which the unit cell has 6.65 ≥ D/d ≥ 6.85. The design of the unit cells of the NETmix network has then a major impact in the performance of the reactor. To characterize macromixing and evaluate the degree of mixing of the NETmix chambers, Batch Lagrangian Mixing Simulations (BLMS) were performed. The results have shown that for 6.65 ≥ D/d ≥ 6.85 the average of the exit conditional probability is close to 0.5. This indicates that a dilution factor of two is obtained between the inlet and outlet of the NETmix chambers, under fully developed flow conditions above the critical Reynolds number. A tracer simulation was performed for the best ratio, D/d 6.65, which shows that nearly perfect mixing is obtained within five rows of chambers, and perfect mixing is obtained with eight or more rows. A new reduced geometric model to represent the NETmix reactor, ExtendedNUB, was presented and validated. It was shown that the ExtendedNUB model simulates the hydrodynamics in the NETmix network without being affected by the inlet and outlet boundary conditions imposed on CFD.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.