Potential of Lagrangian Analysis Methods in the Study of Chemical Reactors

In most chemical reactors the homogeneous distribution of reactants is a desired prerequisite for guaranteed product quality and purity. However, oftentimes, especially in new reactor concepts this condition may not be met due to evolving or existing structures in the transport patterns. In this article the usage of Lagrangian analysis techniques in order to analyze chemical and biochemical reactors regarding the behavior of Lagrangian tracers is explored. A Lagrangian analysis on two different datasets is performed trying to show the large versatility of the concepts. The Lagrangian analysis techniques allow new insights into the reactor dynamics and could be a valuable tool to identify and characterize areas of different mixing performance.


Introduction to Lagrangian Analysis
Lagrangian analysis is used in a variety of studies concerned with transport phenomena and mixing in unsteady velocity fields.The applications of these analysis techniques reach from physiological studies on the blood flow of the heart to the characterization of ocean currents and atmospheric phenomena [1,2].Interestingly, the recent progress in Lagrangian techniques has had little repercussion in chemical process engineering and to date the authors are only aware of one single published study that tries to apply Lagrangian analysis to a laboratory-scale reactor [3].
The basic concept of Lagrangian analysis was developed when dynamical systems theory was applied to steady incompressible flows.In this time-independent flows, stable and unstable manifolds separate regions in the phase space with different flow and mixing dynamics [1,4].Relevant separatrices in phase space can be analyzed by the rate of separation of particles for infinite times.In the case of a two-dimensional fluid flow the phase space is simply spanned by the v x and v y components of the velocity fields.Structures such as stagnation points (points of zero velocity), which are classified as saddle points or hyperbolic points if approaching and departing streamlines (stable and unstable manifolds) intersect, order the flow field by dividing it into topologically and dynamically different subdomains [1,4].An example of a time independent flow field and the corresponding unstable (vertical) and stable (horizontal) manifold is given in Fig. 1.At the unstable manifold tracers approach each other, while they separate along the stable one.
For unsteady velocity fields v(x) = v(x,t) the classical quantities such as the streamlines of the mean velocity give an erroneous picture of the ordering structures of the transport because the trajectories of tracers _ x ¼ u x; t ð Þ in the fluid flow deviate from the streamlines.Therefore, the concept of stable/unstable manifolds in finite times was developed [1] trying to grasp the relevant ordering and coherent structures that are subject to temporal changes.The most intense of these ordering structures were lately termed Lagrangian coherent structures (LCS) and are believed to play a crucial role for mixing dynamics since they can act as material boundaries where the mixing of different fluid masses or advected species is hindered.This characteristic has led to the application of this concepts to, e.g., oceanography where these coherent structures shape the pattern of a large-scale plankton bloom, see Fig. 2 [5].
Modern theory about these ordering structures in finite times departs from the Cauchy-Green strain tensor, which is defined as the square of the gradient of the flow map F t 0 þt t 0 , also called deformation gradient tensor [1,6]: The flow map F t 0 þt t 0 x 0 ð Þ: ¼ x t 0 þ t; t 0 ; x 0 ð Þis the mapping that takes an initial particle distribution x 0 at time t 0 to its distribution x at a time t 0 + t.The deformation gradient of the flow map F t 0 þt t 0 is thereby the linear map that evolves an initial perturbation x 0 + dx 0 to the perturbation of the trajectory at a later time t 0 + t.For a small initial perturbation dx 0 (t 0 ) the flow map can be approximated with a Taylor expansion around x 0 in order to get the final perturbation dx(t 0 + t).
The gradient of the flow map can, therefore, be approximated at each point by: Here, x ¼ x 1 x 2 and dx ¼ dx 1 dx 2 (cf.[1] for a definition in three dimensions).
The concept can be intuitively understood when one thinks of a small vector that is stretched and rotated in a flow field during a finite time t.If we want to know its elongation or contraction, we can simply advect two tracer particles from the bottom and the tip of this small initial vector to their final locations at time t 0 + t.If we now take the finite difference of the new tracer positions in each dimension, we get the deformation gradient tensor.This tensor is not objective, i.e., its values can change under a rotation or translation.However, its square, the Cauchy-Green strain tensor introduced above, is objective and is therefore computed instead.This tensor evolves the square of the small perturbation to its new value at t 0 + t.If one is just interested in the maximal stretching (contraction) that the fluid parcel starting at its center position x 0 = x(t 0 ) experiences during time interval [t 0 , t 0 + t] then it suffices to know the eigenvalues of the Cauchy-Green strain tensor l i , i = 1,2.The ''rate of stretching (contraction)'' during the finite time-interval of interest is then just the square root of the largest (smallest) eigenvalue: . This value forms the basis to the so-called finite time Lyapunov exponent (FTLE), which can be understood as the exponential stretching rate of two initially close tracer particles in forward or backward finite time: If the FTLE is computed for each tracer in a dense tracer field a map of the stretching intensity for the time interval [t 0 , t 0 + t] results, as shown in Fig. 2a.The time interval suitable to analyze the FTLE is thereby somehow arbitrary and depends on the structures and the time scales that are of interest in each process and flow field.The so-called forward FTLE is obtained by taking t > 0 and integrating the velocity forward in time.On the other hand, the backward FTLE field is obtained by integrating backward in time, i.e., t < 0. The forward (backward) FTLE reveals the strongest stretching rate in forward (backward) time and represent regions of divergence (convergence) [4].Thus, they form a time-dependent version of the stable and the unstable manifolds in classical dynamical systems theory as discussed above.
Traditionally, the ridges of high FTLE values where identified with LCS, the material lines that order the flow.However, recent research has found some inconsistencies with this approach and more intricate methods have been developed that take into account the eigenvectors of the Cauchy-Green strain tensor [1,[6][7][8].From this newer concepts LCS are defined to be the most repelling, attracting or shearing material lines of the tracer field in finite time [1,7,9].The most attracting and repelling material lines in the fluid flow are called hyperbolic LCS.They are calculated from the eigenvector fields x i of the Cauchy-Green strain tensor field as detailed in [7] by solving the differential equation _ r ¼ x i r ð Þ.A sketch of the concept considers the fluid parcels as little spheres that are stretched in one direction of the first eigenvector and squeezed in the perpendicular direction of the second, orthogonal eigenvector during the time of integration t = t 1t 0 [1].
Other recent developments in the theory such as the trajectory encounter volume [10,11] and generalized Lagrangian coherent structures [4] expand the concept of purely advective LCS.Another Lagrangian measure that might be advantageous for the characterization of chemical and biochemical reactors are the synoptic Lagrangian maps (SLM).The SLM approach proposed by [12] allows to estimate exit times or escape times through a meaningful region, which has a great importance in chaotic advection and great potential in chemical engineering where mean residence times are an important characteristic [13].The procedure to obtain SLM consists on labeling tracers according to the domain occupied at an instant of time.For this purpose, the domain is split into meaningful regions and each region is assigned with an identifier.The criteria to divide the region could be, e.g., considering the different dynamical regions separated by ridges on the FTLE field.At some initial time, tracers are labeled according to the region where they depart from.While tracers are advected, they move to the other regions and adopt the label of those specific regions as second labels.Eventually, a measure of the occupancy rates is obtained and measures the gain or loss at each region based on the assigned label of the tracers.Further, one can estimate the mobility patterns of tracers inside the domain by plotting a field of the initial tracer positions and its final position identifier as a color code, or vice versa, the final tracer distribution and the initial identifier as a color code.
In this research article two distinct methods are used in order to extract the coherent and ordering structures in two different fluid flows, both time dependent.The first dataset is an experimentally derived two-dimensional velocity field behind a Taylor bubble.Here, a Matlab toolbox, LCS-tool [6], is used to study the LCS as defined by the newer methods relying on distinct paths in the eigenvector space [1,6].Further, also the FTLE values are calculated.Additional insight is gained by calculation of the forward residence times, which are calculated by adding a small Matlab analysis script after the integration of the particle paths by the toolbox.The second dataset is the velocity field from a numerical simulation of a reactor with structured packings introduced in the next section.Since here the data and the relevant fluid flow is three-dimensional, the relevant LCS is extracted by using the FTLE values from a new Python toolbox LAGAR developed by Garaboa-Paz et al. in the Nonlinear Physics group (https://github.com/DanielGaraboaPaz).Additionally, also synoptic Lagrangian maps and residence times are computed.

Introduction to Taylor Bubbles
The understanding of mass transfer at the level of single bubbles is a necessary ingredient in order to improve the understanding of multiphase reactors with a bubbly flow.In chemical and biochemical reactors, it is recurrently reported that undesired side products might be formed if the reaction kinetics and the transport processes are not set optimally.Therefore, it has been suggested [14] to study the interaction of chemical reactions, mass transfer and mixing in an idealized situation using well-controlled elongated bubbles in round channels, so-called Taylor bubbles.One major advantage of these bubbles is that their rise velocities are independent of their volume and purely controlled by varying the Eo ¨tvo ¨s number , where s is the surface tension, r L and r G are the densities of the liquid and the gaseous phase, g is the gravitational constant and D is the diameter of the round channel.It is known that Taylor bubbles start to rise from a critical Eo ¨tvo ¨s number, Eo crit = 4, onwards.For an aqueous solution, this critical Eo ¨tvo ¨s number corresponds to a minimal inner diameter of the round channel of D = 5.4 mm.For elongated Taylor bubbles at low Eo ¨tvo ¨s numbers, no shape oscillations occur at the bubble interface.The bubble is self-centering within the round channel and during the dissolution process only the length of the bubbles decreases gradually.
In the current results we focus on the fluid dynamics and the mixing in the wake of the Taylor bubble for different Reynolds numbers Re = vD/n by using three different round channel diameters [15].For the analysis of the different wake structures the concepts of Lagrangian analysis introduced above are applied.Even though the Taylor bubble experiment is used as an idealized experiment here, several reactors in industry exist that use the same concept on large scale, e.g., multitubular reactors, monolithic reactors [16,17].

Introduction to New Reactor Design Using Periodic Open Cell Structures
The upcoming needs of the chemical industry and environmental issues such as global warming have accelerated the development of new green technologies more efficient and respectful with the environment.Many contributions in the field of reactors design have been suggested in the past in order to improve the performance of chemical reactors.In this context, structured packings reactors haven been shown as an optimal alternative against unstructured stirred tank reactors (STR) for some applications such as the hydrogenation process in the petrochemical industry as is shown in [18,19] (e.g., monolithic honeycombs, open cell foams and periodic open cellular structures (POCS)).Structured packing reactors overcome the limitations on heat transfer of randomly packed-bed reactors, having a major impact on industrial catalytic processes [20][21][22] in addition to an easi-er installation.Further, the advanced additive construction techniques used for POCS allow to print these structures basically with any desired shape and also with new materials allowing to adhere chemical and biological catalysts such as enzymes directly to the structure.POCS reactors combine beneficial properties such as a high specific surface area and low pressure drop, conferred by its periodical shape, with the general benefits of structured packing reactors (i.e., good heat transfer) [23,24].Furthermore, Inayat et.al. have shown how the strut shape and tortuosity is affecting the pressure drop of periodic open-cell foams [25].
The design of this novel reactors consists on a periodical arrangement of unit cells distributed in the three dimensions of the space.The unit cells are composed by intersecting bars, so-called struts.Despite of their promising features, their usage is not widely extended in the industry because many parameters such as mass transfer, energy transfer or heat transfer are still not well understood as well as the correlations among them.Although the insights on POCS are limited, it is known that an influence of the struts geometry against pressure drop exist as it was shown in the work of [26].In the same way, a correlation between struts geometry and mixing distribution is expected.The mixing distribution, influenced by the geometry of the struts and the fluid flow regime, is a crucial factor in chemical or biochemical reactions that require contact with the catalysts attached to the struts.Tracers inside these compartments are constrained to move through it reacting just with those tracers that lie within this region.Thus, an inhomogeneous distribution on mixing is traduced on an inhomogeneous distribution on the chemical reaction rate, which in turn will reflect on the efficiency and the predictability of a chemical process.
Quantifying mixing inside a POCS reactor and determining the influence of struts on it is a challenging task both experimentally and numerically since it involves tackling with a multiphase flow through an intricate structure.In this study the regions that play a role in separating areas of different mixing behavior in a single-phase flow are obtained by using the Lagrangian coherent patterns analysis as introduced above.The dataset is represented by the 3D velocity field of a single-phase flow through a POCS reactor.The Lagrangian analysis focuses on the transport and mixing patterns within a unit cell.With the aid of the Lagrangian analysis two areas of fundamentally different transport dynamics are revealed: a trapping region over the central strut of the unit cell and a region on the outer part of the unit cell with low residence times.In order to reinforce the information drawn from the LCS analysis further Lagrangian measures have been obtained: residence times of passive tracers through the unit cell and the SLM.An analysis of residence times and SLM reveals consistency with the Lagrangian patterns observed.The implications and the potential applications of this methodology on the predictability of a chemical process are discussed.The experimental setup has been thoroughly described in two recent articles [3,14].The heart of the setup is an interchangeable measuring cell consisting of a vertical round glass channel of 300 mm length and variable diameter (here D = 6 mm, 7 mm, 8 mm were used).The liquid flow rate can be adjusted manually with a control valve producing a countercurrent, allowing to keep the CO 2 bubbles in the field of view of the camera for several seconds until the bubble size is so small that shape oscillations occur.
The 2D velocity field and concentration wake structure were obtained simultaneously by means of particle image velocimetry (PIV) and laser-induced fluorescence (LIF) [3] using tracer particles (microParticles GmbH, PS-FluoRot-3.0,d p = 3.16 mm) and a fluorescein sodium salt solution (Sigma Aldrich , 10 -5 mol L -1 ).For the generation of a light sheet a Nd:YLF laser (Darwin Duo Quantronix) is used.It illuminates a planar area around the hydrodynamically fixed Taylor bubble.The laser light sheet is adjusted to be as thin as possible (< 1 mm) using optics from ILA 5150 GmbH.The emitted light from the fluorescence dye and the tracer particles is recorded with a PCO Dimax HS2 at a rate of 600 fps (dt = 1/600 s = 1667 ms) oriented perpendicular to the laser light sheet.
The mass transfer from the CO 2 bubble into the liquid phase suppresses the fluorescence of the fluorescein and, therefore, indicates the mass transfer based boarder between the liquid bulk phase (bright) and the CO 2 -enriched liquid wake region (dark).A bandpass filter protects the camera from direct laser light and reduces the spurious signals from reflections at the round channel walls.The velocity fields were analyzed from the particle images using the software PIVview 2C 3.63, PIVTec GmbH, with parameters as stated in [3].The resulting accuracy of the PIV was assessed in the case of the PIV measurement of the lowest Reynolds number via the signal to noise ratio (SNR), which has a Gaussian distribution with a mean value of 40.Also, a correlation coefficient of about 0.9 and the peak 1/peak 2 ratio higher than 10 everywhere apart from areas very close to the channel boundaries assure valid measurements.Further, the relative error in the velocity measurement was directly estimated as detailed in [3] and it was found to be lower than (2-3 %).

DNS Simulation
In order to obtain the instant velocity fields of a singlephase flow (water) through a POCS reactor a DNS (direct numerical simulation) was performed using the software Ansys Fluent.The computational domain corresponds to a piece of a POCS reactor geometry used in the Institute of Multiphase flows (TUHH) for experimental researching.The unit cells design consists of a sphere surrounded by 12 crossing bars pointing towards the edges of the unit cell.They are geometrical distributed following a periodical pattern in the three directions of the space.
The subdomain of the full reactor was chosen as a tradeoff between accuracy (DNS instead of a large eddy simulation) and computational effort.The subdomain contains 12 unit cells (2 2 3) and is sketched in Fig. 3a.Since the mean flow direction points upwards in the z-axis direction, the structures in the flow field are expected to be elongated in this direction.Thus, three unit cells are distributed in the z-axis direction while only two unit cells expand into the xand y-directions.The edges of one unit cell have a size of L u-c = 5 mm and the strut diameter size is d s = 0.6 mm.The liquid flow rate is passively adjusted by imposing a pressure drop along the z-axis direction.This pressure drop value was adjusted in order to obtain a transitory fluid flow regime, favorable in terms of mixing and transport.A regime of this type starts to form for a strut Reynolds numbers in the range of 200-300 as was observed in preliminary tests (not shown).
For the determination of the Reynolds number the strut diameter has been selected as a characteristic length.Each strut acts similarly as a cylinder in a transitory regime (starting above Re = 47, based on cylinder diameter) that causes the separation of particle trajectories in its wake [27,28].The strut Reynolds number is defined as where v h i is the average velocity, obtained from the mass flow across the inlet, d s is the strut diameter and n is the kinematic viscosity of water at standard pressure and temperature conditions.At these low Reynolds numbers, the flow is assumed to be incompressible.Further, the contribution of gravity is not considered since we are seeking for a characterization independent to the position of the reactor.Applying the incompressibility condition given to the Navier-Stokes equations, the resulting equation governing the motion of the fluid parcels is: Here, m and r represent the viscosity and density of water at standard pressure and temperature conditions: m water = 0.89 mPa s -1 , r water = 0.997 g cm -3 .
For all six external faces (Γ i ) of the computational domain (Ω), shown in Fig. 3, periodic boundary conditions are implied.A no slip condition is imposed on the struts surface, defined as S. Boundary and initial conditions of the system are mathematically expressed by the following equations: where zjj v mean flow and L Length in z (10) Pressure drop and the initial velocity field have been set to ΔP/L = 16 000 Pa m -1 and v 0 = 0.1 m s -1 .
The DNS was performed with a finite volume method (FVM) based software (Ansys Fluent).Under the FVM approach, the numerical methods are especially sensible to mesh parameters such as skewness and element quality.Skewed elements may result in poor convergence and potential loss of accuracy.For this reason, the mesh was designed ensuring that the skewness values never exceed a reference or target value of Sk ref max ¼ 0:9, which is the default target value suggested by Ansys Fluent.In the same way, the orthogonal quality plays an important role in the reliability of the simulation results.The average orthogonal reference value has been set to 0.5.DNS requires to well resolve all the spatial scales of the fluid flow system.Consequently, the mesh was designed ensuring a maximum face size less than the Kolmogorov scale h k,s .For our settings a Kolmogorov scale of is obtained (Re ~300, d s = 6 Á 10 -4 m, n = 0.89 Á 10 -6 m 2 s -1 ).
The mesh consists of unstructured tetrahedral cells and was designed in the Workbench mesher tool (Ansys).In the region nearby the struts the most relevant gradients in the velocity field are in the direction normal to the fluid flow due to the boundary conditions.The boundary layer has a great importance both in mean flow and in turbulence statistics since integral length scales and power spectra of turbulence will be highly affected by the presence of the boundary layer and for a fairly accurate result it needs to be highly resolved.This is accomplished using so-called prism layers.Here, eight of such prism elements layers were applied in the direction normal to each strut wall.Two different meshes made by tetrahedral elements have been designed to test the mesh sensibility.The difference in average mass flow rate at the inlet for both meshes was below 0.0005 kg s -1 , which referred to maximal deviation of the average mass flow of ~1.2 %.Tab. 1 shows an overview of  the most relevant parameters on the coarse (m1) and fine mesh (m2): number of (N n ), number of elements (N e ), average mesh size (L), average skewness Sk h i, maximum skewness (Sk max ), reference skewness (Sk ref max ), average orthogonal quality ( OQ h i) and orthogonal quality reference OQ h i ref .
The convective term was discretized with a first order upwind method and the gradient terms were computed using the Green-Gauss node-based method [28].The simulations were run in parallel in the cluster of the Institute of Applied Mathematics at the University of Santiago de Compostela in Spain and the computation time ranged from 5 h to two weeks.The total integration time was one month.Convergence condition is roughly met in a rate of 30 iterations per hour.At this stage of the simulation the Reynolds number and the mean velocity still show some oscillations around a mean value v mean 0.4 m s -1 .However, the usage of Lagrangian analysis methods is unaffected by stationarity of the flow field and moreover, in many experimental realizations this stationarity is never fully reached.

Lagrangian Analysis for the POCS Reactor
In order to capture the fluid flow topology inside a unit cell and the dominant structures that order the flow local Lagrangian measures were calculated from the 3D velocity fields drawn from the DNS such as the finite-time Lyapunov exponents (FTLE field), the local residence times distributions and the synoptic Lagrangian map (SLM).The resulting FTLE field is dependent on the integration time t and usually this integration time is adjusted by taking into account the temporal scale of the phenomena under study.The current research lacks information about a specific temporal scale since a distinguished chemical or biochemical process is not considered.Consequently, the integration time has been chosen somewhat arbitrary as the value at which the mean FTLE saturate at regions of interest behind the struts where catalysts are attached.The integration time was set to t = 72 Á 10 -4 s and the time step of the integration to dt = 9 Á 10 -5 s, sufficiently small to assure a good accuracy of the tracer trajectories.At an initial time t 0 particles are equidistant distributed on a grid of 100 100 360 spread out over a unit cell and are advected in forward and backward time t 0 + t and t 0 -t, respectively.The computations have been performed on an Intel Xeon CPU E5-2640 v4 cluster with 20 cores and 96 Gb of RAM.The average computational cost for the numerical setup mentioned is in average 1.5 min per updated particle position field being in total close to 2.5 h of computation to achieve one FTLE field of the targeted integration time t = 72 Á 10 -4 s.

2D Velocity Data from a Taylor Bubble Wake
The wake structure of fluid flow behind the Taylor bubble depends on the Reynolds number and thus, for a determined liquid/gas pair, solely on the round channel diameter.The boarder between the bulk phase (bright) and the different wake structures (dark) can be immediately recognized by the grayscale of the experimental raw data (see Fig. 4 in the background).
From the analyzed instantaneous velocity fields in Fig. 4a it becomes even more obvious that the almost laminar flow regime at low Reynolds number (Re = 36) develops into a pair of counter-rotating vortices (Re = 156), which then becomes more and more disturbed in time for the highest Reynolds number considered (Re = 303).In this study, the different transport patterns associated with these different flow situations are of interest.Therefore, the LCS-toolbox [6] was used to analyze the coherent structures in the flow.Fig. 4b shows the most prominent Lagrangian coherent structures (LCS) for the intermediate Reynolds number (Re = 156).The temporal evolution of a patch of tracers advected with the flow field illustrates the effect of the repelling and the attracting manifolds extracted from the data.Tracer initiated above the very coherent repelling (stable) manifold (red) travel towards the bubble, while those below this manifold are rapidly flushed away by the downward flow.Both, the upward and the downward flowing ones are evolving along the attracting (unstable) manifold (blue).The strongest LCS detected also coincides with the highest FTLE values (only shown for the forward FTLE field in Fig. 4c.The FTLE field was calculated for an integration interval of [t 0 , t 0 + t] where t was set to 120Δt and Δt = 1/600 s was the time between two consecutive images, here t 0 corresponds to the images shown in the background of Fig. 4a, respectively.
The coherency of the LCS for the intermediate Reynolds number can further be analyzed by looking at the mean and the standard deviation from the mean of the FTLE values Fig. 5a,b.To compute a mean FTLE field all FTLE fields for intervals t 0 i ; and the standard deviation of every entry in the field was calculated (N = 1000).While for the two larger Reynolds numbers there is a considerable horseshoe-like FTLE ridge surrounding the counter-rotating vortices there is almost no attenuation in the lowest Reynolds number case.However, the temporal coherency of To analyze the impact of the existence of such a coherent structure onto a chemical process we further additionally analyzed the residence times as also done in [3].The residence times T R x kl ; y kl ; t 0 i À Á for the tracers departing from an equidistant grid (x kl , y kl ) at times t 0 i À Á were calculated by assigning the time when a particle leaves the field of view to its initial grid position.From the different residence time maps a mean residence times map can be calculated in a similar fashion as the mean FTLE: Obviously, the coherent structure has a large impact on the residence times of the tracers in the bubble wake and for the intermediate Reynolds number (Re = 156) the residence times take on the highest values.This fact might be counterintuitive since it tells that even though the mean flow rate is sped up there is more liquid trapped in the reactor for longer times (Fig. 5c, middle panel).However, speeding up the mean flow rate even more results in overall shorter residence times, as intuitively assumed (Fig. 5c, right panel).This simple example highlights the necessity of considering eventual coherent structures forming in a reactor, since partially higher residence times could lead to unwanted effects such as enhanced formation of side products [29].

3D velocity data from DNS of POCS reactor
The fluid flow topology has been identified through a Lagrangian analysis based on the calculation of FTLE, SLM and residence times for passive tracers.Passive tracers are advected in the time-dependent velocity field (Fig. 6a).They are initially positioned equally spaced and labeled with colors in a unit cell as shown in Fig. 6b and advected to the final positions at t = t 0 + t as represented in Fig. 6c.It can be observed that after the integration time t particles have moved and mixed across the domain.
Interestingly there is also a considerable transport in the opposite direction to the mean flow.Additionally, some of the tracers that started at the same height (z-location) are advected downstream rapidly while others are trapped behind the main sphere like strut in the middle of the unit cell.We remind that all the Lagrangian quantities obtained in this study are a function of the initial positions of the tracers such that the FTLE fields, the SLM and the residence times are referred to the positions where particles were initially situated.Generally, the fields (FTLE fields, SLM fields and residence time maps) are in their basic structures independent of the details of the particle positions within the seeded area and quite robust against poor resolution of data in space and time as has been shown previously [30].Increasing the resolution of particle initial locations will result in slightly sharper FTLE, SLM and residence times distribution.However, the integration time t is a crucial parameter and the structures in the fields may change considerably while the temporally most coherent ridges will appear more confined.
The ridges in the forward and backward FTLE fields of one unit cell reveal regions of divergence and convergence, respectively, since the FTLE measures the exponential stretching rate that a fluid parcel experiences.A high value of the forward FTLE magnitude in a point of the domain hints that tracers initially close to each other will be far apart in the future.Regions of maximal divergence (stable manifolds) are observed in the forward FTLE field near the faces of the unit cell (Fig. 7a).These veil-shaped manifolds act as barriers of Lagrangian transport and separate the tracers that merely pass by the unit cell from those that become entrapped in the vertical dynamics behind the center sphere shaped strut.The part of the domain where tracers just pass by rapidly will be loosely termed river shape structure.These rivers and the separating veil give rise to an inhomogeneous mixing distribution.This riverlike behavior also already evidenced in the position of particles after being integrated in forward time as can be appreciated in the upstream direction in Fig. 6c.Above the main sphere strut in the center of the unit cell a cone-like structure develops (better visible in Fig. 8).
The manifolds corresponding to the high FTLE ridges that give rise to these separate dynamical regions evolve with time suffering transformations.The dependency of the manifolds with time is an essential distinctive of the theory and provides a realistic portrait of the mixing distribution in real time-dependent fluid flows.This is demonstrated in Fig. 8, where the forward FTLE field was calculated for six different time series f g and dt = t.It can be observed how the FTLE ridges above and below the central strut are mitigated and experience deformations for the different time intervals of integration.Despite the mitigation, the principal coherent veil structure on both sides of the strut changes little, remains coherent and of similar strength for all time intervals.
To assess the meaning of the veil structure for the mixing dynamics the forward residence times are analyzed.The residence times are defined as in the twodimensional dataset in Sect.3.1, but here the residence times map is three-dimensional like the underlying flow data.
Fig. 9a,c depicts the initial position of the tracers colored according to their forward residence times in the unit cell.Tracers initially positioned over struts undergo much higher residence times while in other regions they remain in the unit cell during very short times.The residence times distribution and the coherent patterns unveiled with the forward FTLE calculation show in different ways how mixing is affected by the struts.The trapping region with veil-shaped structure causes high residence times over the central strut while     The SLM analysis introduced above allows to estimate the exchange of tracers between relevant regions in the domain of interest.Further, they can also serve to determine the exit times or escape times a meaningful region, which has great potential in chemical engineering where mean residence times through different parts of a reactor are an important characteristic [13].The methodology to obtain SLM consists on labeling tracers according to the domain occupied at a distinct instant of time.For this purpose, the domain is split into meaningful regions and each region is assigned with an identifier.In this way, by looking at the tracer identifiers we can see at first glance their initial position and their fate.The criteria used here to divide the region into different parts is based on the different dynamical regions obtained by the FTLE field analysis.
Five meaningful regions have been distinguished as represented in Fig. 10a.The blue region (1) partially contains the cone region, the gray region (2) covers a stable manifold under the strut (acknowledged in Fig. 9b), the pink (3) and red (4) regions are not covering any FTLE ridges but are expected to have a distinct dynamical behavior compared to the dynamics inside the blue and gray region for longer times.It is expected from the analysis of the FTLE fields that particles in this transitional regions between the two strong FTLE structures (cone and stable manifold) are susceptible to both, being trapped in the cone-like structure or being washed away quickly on the outside of the cone.Thus, it is of interest where the particles have come from that enter in these transitory regions.The transparent region (0) corresponds with the rest of the unit cell.The SLM in Fig. 10b depicts the initial position of the tracers colored according to the region reached at t = t 0 + t and the values of the forward FTLE at the same time interval.
Observing the light blue colored SLM region (1) it can be observed how particles placed initially intersecting the FTLE ridges (the veil) are expelled while the others have not left the cone region yet.Particles initially positioned in the red (4), gray (2) and pink (3) region in Fig. 10b are mainly expelled.However, even the tracers in the rest of the unit cell, which are not colored, barely protrude into the inner region due to the fact that the veil-shaped structure detected in the FTLE field isolates them.We suggest that the different regions unveiled by Lagrangian analysis could form an objective basis for the definition of compartments within a reactor for further simulation with the aid of compartment models.Our analysis of a POCS reactor reveals two basic compartments each with similar inner mixing dynamics, i.e., the inner area of the cone, and the outer area (rivers), A third compartment could comprise the transition region (similar to region 3 and 4 of the SLM).

Conclusion
Lagrangian-based analysis as presented in this study allows to dynamically outline a realistic time-dependent scenario of the transport patterns that might develop due to boundary constraints.These patterns influence the mixing and the residence times in different parts of reactors as has been shown for two different datasets.For the two-dimensional dataset from laboratory experiments on the flow behind a Taylor bubble a peak of residence times for intermediate Reynolds numbers and flow rates is observed.The Lagrangian analysis reveals a strong repelling manifold where tracers are trapped.This coherent structure develops due to a counter-rotating vortex pair that loses its temporal coherency as the Reynolds number is further increased as can be seen from the standard deviations of the FTLE fields.The elevated residence times behind the Taylor bubble due to the coherent structure for intermediate Reynolds numbers might affect the efficiency of the reactor and the formation of side products of chemical reactions.
The second dataset is derived from a three-dimensional numerical simulation of a simple POCS reactor.Here, the coherent transport patterns are analyzed using forward and backward FTLE fields, SLM and residence times analysis.All three Lagrangian analysis methods reveal an inhomogeneous mixing distribution for passive tracers, which might have crucial consequences on the yield and selectivity of a chemical or biochemical process.Interestingly, the details of the flow topology analyzed via FTLE and SLM indicate the existence of a coherent pattern over the central strut, explaining long residence times in this region.In other words, particles initially placed in the volume constrained by the FTLE ridges over the central strut get trapped and little exchange of tracers is produced between the trapping region and other parts of the unit cell.The analysis of the ridges on the FTLE field also reflects the presence of a coherent pattern with veil-shaped structure nearby the faces of the unit cell that combines with the lower part of the cone structure.The veil structure also acts as a constraint to the transport of passive tracers.The SLM analysis shows that little exchange of tracers is produced between the inner and outer region of the unit cell consistently with the veilshaped and cone-shaped coherent patterns observed.Therefore, tracers initially positioned on either side of the veil do not mix.This inhomogeneous mixing distribution hints that the efficiency and selectivity of a chemical process will be highly affected resulting in higher interaction rates with catalysts printed onto the structure in the trapping region (cone-shaped structure) and lower ones in the region enclosed by the veil-shaped structure (river-shaped structures).It might be desirable to take these insights into account for the future design of POCS.One possible adjustment might be a certain degree of blockage by slight misalignment of the successive unit cells.
To conclude, a novel set of analysis methods to characterize the transport patterns of reactors is suggested.We believe that the introduced methods are especially useful for the detection of dead zones where fluid does not mix with the surrounding fluid or, vice versa, never meets an immobilized catalyst.In addition to that, sketching the different dynamical regions inside a reactor quantitatively may pave the way to establish problem-oriented design criteria for new reactors and might be useful for an objective declaration of different compartments in numerical compartment simulations.

Figure 1 .
Figure 1.In time independent velocity fields (black arrows) the streamlines (light grey lines), the curves that are tangent to the velocity field vectors at each point, coincide with the trajectories of tracers without inertia (round markers).The dark grey tracers were released at (x,y) = (-10.0,±0.2) mm while the light grey tracers were released at the same time at (x,y) = (-10.0,±0.2) mm.The light grey tracer in the upper right quadrant was released at (x,y) = (-10.0,±2.0) m.The dark vertical line indicates the unstable, attracting manifold, where the tracer of different origin come together, the lighter horizontal line defines the stable, repelling manifold where tracers of similar origin separate strongly.

Figure 2 .
Figure 2. FTLE analysis has been performed on velocity data of ocean currents derived geostrophically from altimetric sea surface height (SSH) maps in order to find the relevant transport pattern that could explain the extreme fast and confined plankton bloom observed in chlorophyll concentration in ocean color data sets from the NASA (http://oceancolor.gsfc.nasa.gov).a) Instantaneous map on 17 February 1999 of added forward and backward FTLE [weeks -1 ] with highlighted jet-like LCS inhibiting meridional transport.(c) SeaWiFS chlorophyll concentration [mg m -3 ] with marked position of the jet like LCS.The image corresponds to an 8-day composite centered around 21 February 1999.(Reprinted from Huhn 2012 [5] with kind permission of Geophysical Research Letters).

Figure 3 .
Figure 3. a) The computational domain mimics a centered part of an experimental POCS reactor.The simulated domain consists of three unit cells in the z-direction and two unit cells in the x-and y-direction.The boundary conditions are periodic in all directions.b) A unit cell and its main dimensions: all the edges have the same size L u-c = 5 mm and the strut diameter is d s = 0.6 mm.c) The black lines enclose the region of the computational domain in which particles are advected in order to calculate the FTLE fields within one unit cell.
the two horseshoe-like FTLE ridges is very different as can be understood by looking at the standard deviation of the mean field, s.For the intermediate Reynolds number (Re = 156) s takes on very small values at the location of the horseshoe-like structure while for the highest Reynolds number (Re = 303) the values of the standard deviation are very high.The difference results from the large temporal fluctuations of the FTLE ridges in case of the highest Reynolds number (Re = 303).

Figure 4 .
Figure 4. a) The Taylor bubble is held in a constant counter flow for three different round channel diameters and, thus, Reynolds numbers, D = 6 mm (Re = 36), D = 7 mm (Re = 156), D = 8 mm (Re = 303).The instantaneous velocity fields are plotted on top of the raw images.Three different wake structures can be observed [3].b) Focus on the region highlighted by the red rectangle in the intermediate Reynolds number wake in a): The red repelling (stable) manifold indicates where tracers separate strongly during successive time steps.The blue attracting (unstable) manifold indicates where the tracers are converging.c) For the intermediate Reynolds number the strongest ridge of the FTLE field and the LCS calculated using the toolbox LCS-tool coincide well.

Figure 5 .
Figure 5. Focus on the wake regions behind the bubbles for all three Reynolds number cases (Re = 36, Re = 256, Re = 303).a) The mean FTLE fields exhibit strong horseshoe-like ridges for the two highest Reynolds number wakes.However, by looking at the standard deviation of the FTLE fields (b), it can be appreciated that only at the intermediate Reynolds number a temporally coherent pattern exists because only here the temporal fluctuations of the ridge are low.c) The strong coherent structure for the intermediate Reynolds number leads to a peak in mean residence times behind the bubble.

Figure 6 .
Figure 6.The velocity field for three consecutive unit cells is shown color-coded in a).b) Tracer positions at the initial time t = t 0 and c) after being advected for a time t = t 0 + t.It can be immediately observed that some particles starting at the same z-location (i.e.marked yellow) are advected downstream rapidly while others are trapped behind the main strut (sphere in the middle of the unit cell).

Figure 7 .
Figure 7. a) Regions of maximal divergence (forward FTLE ~18-26 s -1 ) are observed close to the faces of the unit cell.b) Geometry of the unit cell for comparison.

Figure 8 .
Figure 8.The forward FTLE field in the central cut plane of the unit cell has been calculated for six different advection intervals.The detailed FTLE values behind the strut change with time while the veil-like region separating the different dynamical regimes (rivers, cones) of the flow changes only little.

Figure 9 .
Figure 9. a) Initial positions of tracers at t = t 0 in the center cut plane of the unit cell, color-coded with their residence times.b) Forward FTLE field for same cut plane and time interval.c) Three-dimensional distribution of residence times in the unit cell.The cut plane is indicated by dashed lines.The ratio of the highest values of the residence times to the lowest ones is t res,long /t res,short > 7. The distribution of residence times is very inhomogeneous and largely influenced by the struts.
times distribution are due to the existence of underlying transport patterns.

Figure 10 .
Figure 10.a) Five SLM regions (0-4) are defined in the unit cell based on the knowledge about the ridges in the FTLE field.b) Here the tracers are plotted at their initial positions but color-coded with their final position after a time t.It is easily seen that most of the tracers started in the light blue region have not left it after a time t = t 0 + t.Additionally, also the forward FTLE values in the indicated cut planes (1-4) are plotted to show the close correspondence of the ridges with the boundaries of the light blue SLM structure.The strut is colored in deep blue.
Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/cite.201900147 by Technical University Eindhoven, Wiley Online Library on [06/03/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2 Materials and Methods 2.1 2D Velocity Data from Taylor Bubble Wake

Table 1 .
Most relevant parameters of the coarse (m1) and fine (m2) meshes used in the computation.