Selecting for lactic acid producing and utilising bacteria in anaerobic enrichment cultures

Potential of Lagrangian Analysis Methods in the Study of Chemical Reactors Cristina Garcia Llamas 1 , Claas Spille 2 , Sven Kastens 2 , Daniel Garaboa Paz 3 , Michael Schlu ¨ter 2 , and Alexandra von Kameke 2, * DOI: 10.1002/cite.201900147 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. 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. Keywords: Reactor dynamics, Coherent structures, Lagrangian analysis Received: September 15, 2019; revised: December 06, 2019; accepted: January 31, 2020 1 Introduction 1.1 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 vxðÞ¼ax y  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¼ux;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]. www.cit-journal.com ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Chem. Ing. Tech. 2020,92, No. 5, 540–553 – 1 Cristina Garcia Llamas Eindhoven University of Technology (TU/e), Multi-Scale Modelling of Multi-phase Flows, Het Kranenveld 14, 5612 AZ Eindhoven, Netherlands. 2 Claas Spille, Sven Kastens, Prof. Dr.-Ing. Michael Schlu ¨ter, Dr. Alexandra von Kameke [email protected] Hamburg University of Technology (TUHH), Institute of Multiphase Flows, Eissendorfer Straße 38, 21073 Hamburg, Germany. 3 Dr. Daniel Garaboa Paz University of Santiago de Compostela, Group of Nonlinear Physics, Xose Marı ´a Suarez Nun ˜ez s/n, Santiago de Compostela, 15782, Spain. 540 Research Article Chemie Ingenieur Technik 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 Ft0þt t0, also called deformation gradient tensor [1, 6]: Ct0þt t0x0 ðÞ¼Ft0þt t0x0 ðÞ  TFt0þt t0x0 ðÞ (1) The flow map Ft0þt t0x0 ðÞ:¼xt0þt;t0;x0 ðÞis the mapping that takes an initial particle distribution x 0 at time t 0 to its distribution xat a time t 0 +t. The deformation gradient of the flow map Ft0þt t0is 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). dxt0þtðÞ»Ft0þt t0dx0t0 ðÞ (2) The gradient of the flow map can, therefore, be approximated at each point by: Ft0þt t0x0 ðÞ» x1t0þt;t0;x0þdx1 ðÞ x1t0þt;t0;x0 ðÞ dx1 jj x1t0þt;t0;x0þdx2 ðÞ x1t0þt;t0;x0 ðÞ dx2 jj x2t0þt;t0;x0þdx1 ðÞ x2t0þt;t0;x0 ðÞ dx1 jj x2t0þt;t0;x0þdx2 ðÞ x2t0þt;t0;x0 ðÞ dx2 jj 0 B B B B B B @ 1 C C C C C C A (3) Here, x¼x1 x2  and dx¼dx1 dx2  (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: Si¼ffiffiffiffi li p. 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: L–x0;t0;tðÞ¼ 1 tln ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l1t0–t;t0;x0 ðÞ p (4) 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 CauchyGreen strain tensor [1, 6–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¼xirðÞ. A sketch of the concept considers the fluid parcels as little spheres that are stretched in one direction of Chem. Ing. Tech. 2020,92, No. 5, 540–553 ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.cit-journal.com 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. Research Article 541 Chemie Ingenieur Technik the first eigenvector and squeezed in the perpendicular direction of the second, orthogonal eigenvector during the time of integration t=t 1 –t 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. www.cit-journal.com ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Chem. Ing. Tech. 2020,92, No. 5, 540–553 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). 542 Research Article Chemie Ingenieur Technik 1.2 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 EoD¼rLrG ðÞgD2 s, where s is the surface tension, r L and r G are the densities of the liquid and the gaseous phase, gis the gravitational constant and Dis the diameter of the round channel. It is known that Taylor bubbles start to rise from a critical Eo¨tvo¨s number, Eocrit = 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/nby 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]. 1.3 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–22] in addition to an easier 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. Chem. Ing. Tech. 2020,92, No. 5, 540–553 ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.cit-journal.com Research Article 543 Chemie Ingenieur Technik 2 Materials and Methods 2.1 2D Velocity Data from Taylor Bubble Wake 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-FluoRot3.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 %). 2.2 3D Velocity Data from DNS of POCS Reactor 2.2.1 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 Re ¼v hi ds n(5) where v hiis the average velocity, obtained from the mass flow across the inlet, d s is the strut diameter and nis 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: r¶v ¶tþvv  ¼pþmvþvðÞ T  (6) Here, mand rrepresent 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 (G i ) of the computational domain (W), 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: vt;xþLiei ðÞ¼vt;xðÞx˛Gi;i¼1;2;3;t‡0 (7) www.cit-journal.com ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Chem. Ing. Tech. 2020,92, No. 5, 540–553 544 Research Article Chemie Ingenieur Technik pt;xþLiei ðÞ¼pt;xðÞx˛Gi;i¼1;2;3;t‡0 (8) pt;xðÞz¼DP Lx˛W;t‡0(9) where zjj vmean flow and L”Length in z(10) vt;xðÞ¼0x˛S;t‡0 (11) v0;xðÞ¼v0x˛W(12) Pressure drop and the initial velocity field have been set to DP/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 Skref 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 hk;s¼n3ds v3  0:25 ¼ds Re0:75  ~105m (13) 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 Chem. Ing. Tech. 2020,92, No. 5, 540–553 ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.cit-journal.com 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 xand 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. Table 1. Most relevant parameters of the coarse (m1) and fine (m2) meshes used in the computation. N n N e L hi [m] Sk hi (Sk max )Skref max OQ hi OQ hi ref m1 17 973 516 75 249 216 2.5 10 –5 0.2 0.85 0.9 0.79 0.5 m2 19 799 701 83 381 513 2.3 10 –5 0.2 0.85 0.9 0.79 0.5 Research Article 545 Chemie Ingenieur Technik the most relevant parameters on the coarse (m1) and fine mesh (m2): number of nodes (N n ), number of elements (N e ), average mesh size (L), average skewness Sk hi , maximum skewness (Sk max ), reference skewness (Skref max), average orthogonal quality ( OQ hi ) and orthogonal quality reference OQ hi 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. 2.2.2 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 +tand 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. 3 Results and Discussion 3.1 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 twas set to 120Dtand Dt= 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 t0i;t0iþt  with t0i˛N120Dt;1120 Dt½were averaged FTLE hi ¼1 NX i FTLEi(14) and the standard deviation s¼1 N1ðÞ X i FTLEiFTLEi hiðÞ 2(15) 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 www.cit-journal.com ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Chem. Ing. Tech. 2020,92, No. 5, 540–553 546 Research Article Chemie Ingenieur Technik the two horseshoe-like FTLE ridges is very different as can be understood by looking at the standard deviation of the mean FTLE field, s. For the intermediate Reynolds number (Re = 156) stakes 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). 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 TR xkl;ykl;t0i  for the tracers departing from an equidistant grid (x kl ,y kl ) at times t0i 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 Chem. Ing. Tech. 2020,92, No. 5, 540–553 ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.cit-journal.com 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. Research Article 547 Chemie Ingenieur Technik a mean residence times map can be calculated in a similar fashion as the mean FTLE: TR hi ¼1 NX i TRxkl;ykl;t0i  (16) 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]. 3.2 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 +tas represented in Fig. 6c. It can be observed that after the integration time tparticles 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 tis 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 www.cit-journal.com ª2020 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim Chem. Ing. Tech. 2020,92, No. 5, 540–553 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. 548 Research Article Chemie Ingenieur Technik