In this research work, we apply Reduced Order Methods (ROM's) to perform a mode decomposition of the uid ow present in two systems widely used in the bio-chemical industry: stirred tanks and bubble columns. This decomposition allows to identify and classify the uid ow components according to their dynamical relevance. The most important structures identied correspond to large-scale structures with a high degree of organization known as coherent structures. Subsequently, the classication of results lead to a reconstruction of the uid ow, the accuracy and computation time involved in the procedure of ROM's were also determined. These parameters are compared with those corresponding to CFD, this comparison shows that ROM's approach implies a lower time cost with a relatively low error. The ROM's chosen are the most recently used for industrial applications: Proper Orthogonal Decomposition (POD) and Dynamical Mode Decomposition (DMD). The implementation of these methods in the baed stirred tank was approached through two strategies because the simulation of the system (perform in ANSYS/FLUENT) has a xed mesh zone and a sliding mesh zone. The rst strategy is the typical or standard algorithm and the second is a new algorithm not reported in the literature. For the rst strategy or typical approach, the POD was applied in the xed zone and in the rotating zone separately. The results allowed the reconstruction of the mean ow and the trailing vortices generated by the passage of the turbine blades. In this case the velocity elds of the rotating zone were suitably transformed to the rotating reference frame of the turbine before the ROM's procedure. For the second and novel implementation the POD and DMD were applied directly to the entire simulation domain without separating the respective zones. According ndings found, this new proposal also allows the reconstruction of large-scale organized structures (mean ow and trailing vortices in the impeller region). However, for this case it is not possible to perform a physical interpretation of the POD eigenvectors as it is the case for the conventional POD algorithm. The results obtained for the POD and DMD using the proposed approach are equivalent due to the presence of highly periodic uid ow, such as a similarity is mentioned (but not veried) in the literature. Additionally, the DMD algorithm allows the identication of low-frequency, low-energetic structures associated with macro instabilities identied experimentally and numerically in other published works. The identication of macro-instabilities in simulations through DMD has not been reported in the literature.

The numerical data for the stirred tank were generated using the U-RANS turbulence model. The choice of this turbulence model over others available such as Large Eddy Simulation (LES) is due to the objective of generating reconstructed velocity elds to be used in future works for compartmentalization models (CMA) implementations. The CMA models can be fed with mean velocity elds calculated by U-RANS simulations or alternatively reconstructed velocity elds generated by order reduction methods. We also apply POD and DMD to the case of bubble plume in a quasi-2D column. An euler-euler two-uid model was solved for the bubbly ow and Rij-epsilon approach was tested as turbulence model. From this approach, it was possible to identify the dominant organized structures (low-frequency structures) for both phases employing POD and DMD. The simulation for this case was taken from a database generated by David Laupsien (PhD 2017), using the NEPTUNE CFD code. Finally we briey address the basic principles necessary for the coupling of CFD with CMA. In this conceptual framework of mixing we implement in the model the particle transport based on the Monte-Carlo approach simulation. This introduction aims to prepare the way for the future implementation of reduced order models POD + CMA or DMD + CMA.

v Acknowledgements/Agradecimientos Since my childhood, I always dreamed of coming to France to study "those things" that appeared in the books of the great writer Jules Verne. For this reason, my passion for this beautiful country began a long time ago. The day has come when this dream came true and many people supported and inspired me to get here. First of all, I must thank my friend and adviser Alejandro García for his prayers, trust, the testimony of life, his word always full of wisdom, his fatherly love, and unconditional care for me.

Similarly, I need to recognize the immense help that my brother Robert has given me who, with his life example, colossal strength, and determination, has encouraged me to face life. Many thanks also to my friend Eddy Araya for his prayers, unconditional support, and conversations on topics from the Chifrijo recipe to Theology... Thank you Bro for always keeping an eye on me. I deeply appreciate all the patience, help, and kindness provided by my supervisors "los jefes" Jérôme Morchain and Alain Liné.

Without them and their commitment, it would have been impossible to complete my project. My deep gratitude also to César Aceves and Erick Mora. Dear friends, you extended the bridge that allowed me to come here and start this crazy adventure. I am really lucky to have you as friends. I can't stop thanking colleagues who helped me a lot in very important technical topics like CFD and Matlab: Francesco Maluta and Mohammed Arezki Chekroun I will be forever grateful to you. And speaking of co-workers, I have to mention my colleagues from oce 208, the happiest oce in the world. Many thanks to Alex, Eliot, Gaëlle, Maike and Maxence. What a blessing to have met you and to have shared so many special and fun moments. Thanks also to my friend Juan Camilo Acosta, you are a super galleta (in Costa Rican language :-D). You guys are truly PURA VIDA! My thanks also to my EAD10 colleagues: Patrick, Plat, Ryma, Seunghye, Shivesh, Ugo, and Zhou. Thank you for all the beautiful moments we shared. I want to thank the TEC Instituto Tecnológico de Costa Rica, especially the colleagues from the Physics Department and the Scholarship Department. Also a big thank you to the Cooperation and Cultural Action Service for Latin America, especially Mrs. Yadira Vargas for her patience and commitment to my process here in France. Finally, my gratitude to INSA, TBI and the MEGeP doctoral school for receiving me and accepting vi me as a collaborator these three and a half years.

Desde mi infancia siempre soñé con venir a Francia a estudiar aquellas cosas que aparecían en los libros del genial escritor Julio Verne. Por este motivo mi pasión por este bello país comenzó desde hace mucho tiempo. Ha llegado el día en que este sueño se hizo realidad y son muchas las personas quienes me apoyaron e inspiraron para llegar hasta aquí. En primer lugar, me es necesario agradecer a mi amigo y consejero Alejandro García por sus oraciones, conanza, testimonio de vida, su palabra siempre llena de sabiduría, su amor de padre y su cuidado incondicional. De manera similar me es indispensable reconocer la inmensa ayuda que me ha brindado mi hermano Robert quien con su ejemplo, fortaleza colosal y determinación me ha animado a afrontar la vida. Muchísimas gracias también a mi gran amigo Eddy Araya por sus oraciones, apoyo incondicional y conversaciones de temas desde la receta del Chifrijo hasta la Teología...Gracias Bro por estar siempre pendiente de mí. Debo agradecer profundamente toda la paciencia, ayuda y amabilidad brindada por mis supervisores "los jefes" Jérôme Morchain y Alain Liné. Sin ellos y su compromiso hubiese sido imposible poder llegar a un buen término de mi proyecto.

Mi profunda gratitud también para César Aceves y Erick Mora. Amigos queridos ustedes extendieron el puente que me permitió venir aquí y comenzar esta loca aventura. Soy realmente afortunado de tenerles a ustedes como amigos. No puedo dejar de agradecer a colegas que me ayudaron mucho en temas técnicos demasiado importantes como CFD y Matlab: Francesco Maluta y Mohammed Arezki Chekroun les estaré por siempre agradecido a ustedes. Y hablando de compañeros de trabajo como no podría mencionar a mis compañeros de la ocia 208, la ocina más feliz del mundo. Muchas gracias a Alex, Eliot, Gaëlle, Maike y Maxence. Que bendición la mía por haberlos conocido a ustedes y haber compartido tantos momentos especiales y divertidos juntos. Gracias también al parcerito Juan Camilo Acosta, eres súper galleta (en idioma costarricense :-D). ½Ustedes son verdaderamente PURA VIDA! Mi agradecimiento también a los compañeros del equipo EAD10: Patrick, Plat, Ryma, Seunghye, Shivesh, Ugo y Zhou. Gracias por todos los momentos bonitos que compartimos. Quiero agradecer al TEC Instituto Tecnológico de Costa Rica, en especial a los compañeros del Departamento de Física y al Departamento de Becas. También un gran agradecimiento al Servicio de Cooperación y Acción Cultural para América Latina, especialmente a la señora Yadira Vargas por su paciencia y compromiso con mi proceso aquí en Francia. Finalmente, mi gratitud al INSA, TBI y a la escuela doctoral MEGeP por recibirme y aceptarme como colaborador estos tres años y medio. viii ix Résumé Dans ce travail, nous appliquons des méthodes d'ordre réduit (ROM's) pour eectuer une décomposition de mode de l'écoulement uide dans deux systèmes largement utilisés dans l'industrie biochimique : les cuves agitées et les colonnes à bulles. Cette décomposition permet d'identier et de classer les composantes de l'écoulement uide en fonction de leur pertinence dynamique. Les structures les plus importantes identiées correspondent à des structures à grande échelle avec un haut degré d'organisation. La classication des résultats conduit à une reconstruction de l'écoulement uide, la précision et le temps de calcul impliqués dans la procédure des ROM's ont également été déterminés. Ces paramètres sont comparés à ceux correspondant à la CFD ; cette comparaison montre que l'approche des ROM's implique un coût en temps inférieur avec une erreur relativement faible. Les ROM's choisies sont: Proper Orthogonal Decomposition (POD) et Dynamical Mode Decomposition (DMD). La mise en ÷uvre de ces méthodes dans la cuve agitée a été abordée selon deux stratégies car la simulation du système (ANSYS/FLUENT) comporte une zone de maillage xe et une zone glissant. La première stratégie est l'algorithme usuel et la seconde est un nouvel algorithme non rapporté dans la littérature. Pour la première stratégie, la POD a été appliquée dans la zone xe et dans la zone rotative séparément.

Les résultats ont permis de reconstruire l'écoulement moyen et les vortex trainants générés par la rotation de la turbine. Cependant, les vitesses dans la zone tournante doivent être exprimées dans le repère rotatif de la turbine avant la procédure de ROM. Pour la deuxième implémentation, la POD et la DMD ont été appliqués directement à l'ensemble du domaine de simulation sans séparer les zones respectives. Selon les résultats trouvés, cette nouvelle proposition permet également la reconstruction de structures organisées à grande échelle. Cependant, dans ce cas, il n'est pas possible d'eectuer une interprétation physique des vecteurs propres de la POD. Les résultats obtenus pour la POD et la DMD en utilisant l'approche proposée sont équivalents en raison de la présence d'un écoulement uide hautement périodique, comme mentionné dans la littérature. De plus, l'algorithme DMD permet l'identication de structures à basse fréquence et à faible énergie associées à des macro-instabilités identiées expérimentalement et numériquement dans d'autres travaux publiés. L'identication des macro-instabilités dans les simulations par DMD n'a pas été rapportée dans la littérature. Les données numériques pour la cuve agité ont été générées en utilisant le modèle de turbulence U-RANS. Le choix de ce modèle de turbulence par rapport à d'autres modèles disponibles, tels que la simulation des grandes tourbillons (LES), est dû à l'objectif de générer des champs de vitesse reconstruits qui seront utilisés dans des travaux futurs pour la mise en ÷uvre de modèles de compartiment (CMA). Nous appliquons également la POD et la DMD au cas d'un panache de bulles dans une colonne quasi-2D. Un modèle euler-euler à deux uides a été résolu pour l'écoulement de bulles et l'approche Rij-epsilon a été testée comme modèle de turbulence.

A partir de cette approche, il a été possible d'identier les structures organisées dominantes (structures basse fréquence) pour les deux phases en utilisant POD et DMD. La simulation pour ce cas a été prise à partir d'une base de données générée par David Laupsien (PhD 2017), en utilisant le code CFD NEPTUNE. Enn, nous abordons brièvement les principes de base nécessaires au couplage de la CFD avec la CMA. Dans ce cadre conceptuel de mélange, nous implémentons dans le modèle le transport de particules basé sur la simulation par approche Monte-Carlo.

by the red rectangle shown (0 ≤ x ≤ 0.17 cm, 0.6 ≤ z ≤ 0.7 cm). III.39POD and DMD (alternative criteria) spectrum mode for 2.3 minutes. The ranking of the modes is not identical for both methodologies. The energy content for the 13 th and 14 th modes is dierent for each data-driven method. The ranking of the modes is identical for both methodologies. . . . . . . . . .

List of Tables

I.1

Values of the constants for the k -ε model. . . . . . . . . . . . . . . . . .

I.2

Steady state approaches for the impeller rotation and baes modelling. . .

I.3

Unsteady state approaches for the impeller rotation and baes modelling.

F

Matrix of owrates: the owrate from a compartment i to a compartment j is placed on the i th line and the j th column of the matrix

V i or V (i) volume of the i th CMA compartment
General introduction

.

Motivation

The rapid advance of industry and technology requires better control of industrial processes and optimizing their designs to achieve higher eciencies. This situation implies a deeper and more realistic understanding of the physical, biological, and chemical phenomena involved. Experimental and numerical approaches are available to meet this need. In particular, uid dynamics research has beneted immensely from these two strategies. On the experimental side, it is worth mentioning that optical or particle imaging techniques have proven for decades to be eective and practical for the study and analysis of uids for many ow conditions and congurations [START_REF] Tayali | Particle sizing techniques in multiphase ows: A review[END_REF]).

Among the various experimental methodologies available, particle image velocimetry (PIV) has gained interest due to its accuracy, resolution, ease of implementation, and analysis of the data generated [START_REF] Grant | Modern developments in particle image velocimetry[END_REF][START_REF] Adrian | Review of particle image velocimetry research presented at the symposium on optical methods in ow and particle diagnostics[END_REF][START_REF] Adrian | Particle-Imaging Techniques for Experimental Fluid Mechanics[END_REF][START_REF] Adrian | Particle image velocimetry[END_REF][START_REF] Scharnowski | Particle image velocimetry -Classical operating rules from today's perspective[END_REF].

Alternatively, numerical tools provide a means for the study of uid dynamics for any possible condition or conguration. The group of numerical algorithms intended to accomplish such a task is known as Computational Fluid Dynamics or CFD. As its name suggests, CFD uses numerical and computational techniques to solve the equations that satisfy the conservation of mass, momentum and energy. There are a variety of approaches to numerically to solve these nonlinear partial dierential equations. For example, it is possible to obtain a numerical solution over the full range of spatial and temporal scales of turbulence, this strategy is known as Direct Numerical Simulation or (DNS). This approach is computationally demanding and remains limited to academic studies. On the other hand, one can use some methods with lower calculation cost to model certain turbulent scales (Large Eddy Simulation) or all turbulent scales (U-RANS) of the ow uid.

It also is possible to incorporate in the CFD methodology other types of phenomena (biological, electromagnetic, chemical, structural, etc.) that must be coupled with the uid dynamics of the system under study. At present, there is a considerable amount of in-house and commercial CFD codes, which shows its popularity and applicability in dierent areas of the industry [START_REF] Freitas | Perspective: Selected Benchmarks From Commercial CFD Codes[END_REF][START_REF] Keshmiri | Benchmarking of Three Dierent CFD Codes in Simulating Natural, Forced, and Mixed Convection Flows[END_REF].

Implementing an experimental or numerical approach to the study of industrial problems implies the need for high data storage and processing capacity [START_REF] Kozelkov | Investigation of supercomputer capabilities for the scalable numerical simulation of computational uid dynamics problems in industrial applications[END_REF]. This situation is precisely the case when seeking an improved understanding of processes leading to better production eciencies. There are two possibilities to address this subject. The rst is to use high-performance computers whose capacity satises the computational demand. This option is generally not the most feasible because highperformance computing equipment is not always available. Alternatively, a second option is to extract the most relevant numerical information from experimental measurements or computational simulations to perform high-delity simplications of the input data.

The information retained will be that which best describes the hydrodynamic behavior according to the scale of the phenomenon of interest. Mathematical models designed to achieve this objective are called Reduced-Order Methods or simply ROM's (Quarteroni & Rozza, 2014;Rowley & Dawson, 2017;[START_REF] Taira | Modal Analysis of Fluid Flows: An Overview[END_REF]Rozza et al., 2018;[START_REF] Taira | Modal Analysis of Fluid Flows: Applications and Outlook[END_REF]. These methodologies have gained great interest in the industrial sector due to the signicant advantages they oer in terms of practicality and convenience. For example, identifying the most relevant hydrodynamic elements of the system makes it possible to use minimal computational resources. This saving of computational resources can be used for the application or improvement of other models complementary to uid mechanics. This strategy is a way to favor the investigation of very complex industrial processes.

Among the most popular ROM's we nd the Proper Orthogonal Decomposition (POD)

and the Dynamical Mode Decomposition (DMD). The former allows generating an order

.2 Objectives of the research project reduction using the covariance (energy) of the input database as a criterion. The latter leads to an order reduction based on the dierent frequencies contained in the data to be analyzed. As will be seen throughout this document, both ROM's approaches work quite well when used for the study of periodic uid ows. This feature of both methodologies makes them good candidates for coupling with mixing models such as Compartment Model Approach (CMA) for which the use of ROM's results could lead to more accurate calculations.

.

Objectives of the research project

For a particular class of problems involving multiple time scales, it would be advantageous to maintain the high spatio-temporal resolution of CFD while running simulation for several days. This typically includes the case of bioreactors in which the repeated exposure to concentration gradients at the scale of a minute may impact the course of the fermentation process which lasts several days. Among, the solutions prevoiusly examined, the Compartment Model Approach is based on the spatial integration of the time averaged velocity eld. Both the spatial and temporal resolution of the initial CFD simulations are ltered through this process and one might wonder if it is possible to preserve both while reducing the computational cost. Reduced Order Methods are ideal candidates since a velocity eld reconstruction is theoretically feasible from a linear combination of modes.

Because ROMs are data-driven methods, the spatial and temporal resolution of the input data sets are preserved. The main objective of this research is the implementation of POD and DMD order reduction methods to analyse and reconstruct the velocity elds using as input data the outcome of 3D and unsteady CFD numerical simulations. We investigated the case of a baed stirred tank and that of a bubble plume in a quasi 2D column.

To this end, we performed the following actions:

Identication and classication of organized structures (coherent structures) of high and low frequency from the above implementations according to covariance (POD)

and frequency (DMD).

Reconstruction of velocity elds in the analyzed systems using relevant dynamical structures from the POD and DMD implementations.

Evaluation of the computational cost of velocity elds reconstruction and comparison to CFD numerical simulation.

Evaluation of the error associated to the POD and DMD reconstruction taking as reference the velocity elds generated in the CFD numerical simulations.

The secondary objective is to investigate the oportunity to combine Reduced Order Models and Compartiment Model Approach to perform time resolved simulations of chemical and bio-chemical reactors. In order to pave the way for dynamic CMA models that considers deviation from the mean ow (periodic structures and macro-instabilities) we propose:

Implement a methodology for a dynamic CMA model using as input dierent uxes maps of a stationary CMA model.

Evaluate the proposed CMA model for particle transport and chemical reactions.

This is intended to show the feasibility of generating CMA models that are fed by reconstructed ROM velocity elds, such a scenario will be addressed in future work.

.

Contributions of the research project

The work composing this thesis has resulted in the following products: CFD can be used in a wide range of industrial and academic applications such as the aerodynamics of aircraft and vehicles, power plants turbines, turbomachinery, electrical and electronic engineering, chemical engineering, biomedical engineering, environmental engineering, geophysical sciences, and even in food industry. CFD modeling has advantages and disadvantages when compared to experimental results. Among the most relevant advantages can be mentioned:

Once the CFD results have been validated with physical experiments, it can be achieved a substantial reduction of lead times and costs for the study and analysis of new designs.

When DNS simulations are implemented, it is possible to study systems where controlled experiments are dicult or impossible to perform. Such as complex cases in the stationary and non-stationary regimen.

High level of detail of the results if one has a ne enough spatial and temporal resolution. In this case "enough" depends on the phenomenon under study.

The possibility of accessing the pressure eld.

On the other hand, this computational procedure also presents problems such as:

Physical simplications are considered when running simulations that include models such as those used for turbulent ows (e.g. RANS, LES, etc.).

Errors in the results due to the process of numerical discretization and truncation.

Time investment for code implementation (when a simulation software is no available), setup, and run of the computation and simulations.

Data processing can be very arduous and involve a large investment of time and computational resources.

In some cases, the interpretation of the obtained results is dicult.

I.1.1 Stages involved in the implementation of a CFD experiment

The application of the CFD procedure (and therefore the respective codes) implies the following stages: pre-processing stage(geometry design, boundary conditions, model assumptions, and mesh design), discretization and numerical solution of the ow equations, post-processing stage, and nally the storage of the data.

The pre-processing stage: The design of the geometric domain is a fundamental aspect because it inuences the numerical modeling and the discretization to be implemented. Additionally, one must consider spatial features such as possible symmetries (applies for RANS models for example), number of spatial dimensions (1-D,2-D, 3-D) that simplify the numerical experiment without compromising the desired physical validity.

Next, it is necessary to dene the assumptions of the model, such as steady-state regime, incompressible ow, inviscid ow, laminar or turbulent ow. In this part, the boundary conditions must be properly dened to ensure compliance with physical laws, such as the conservation of mass, momentum, and energy. Once the geometry, boundary conditions and model assumptions are dened we can discretize or divide the domain space into unit cells. This step is immensely complex because there are almost innite ways of segmentation for any proposed geometry. Aspects such as the size distribution and quality parameters of the unit cells that make up the mesh further complicate this step.

To exemplify this geometric process of spatial segmentation, let us look at the gure I.1.

The gure shows two dierent ways of geometric partitioning for a 2-D annular domain generated in the ANSYS mesher. In the rst image, one observes a quadrangular mesh with a clear size distribution of the elements. Also in the second image, another type of geometric segmentation consisting of a distribution of triangular elements is proposed.

The renement of the mesh at domain boundaries is related to the implementation of the boundary layer. It is easy to imagine that an almost innite number of mesh congurations can be realized for the indicated domain. According to [START_REF] Versteeg | An introduction to computational uid dynamics : the nite volume method[END_REF], this meshing process can involve about 50% of the total time spent on a CFD implementation. In general, a very ne meshing leads to more accurate results but the size involved may compromise the available computational resources 1 . In the absence of experimental data to validate a numerical simulation, it is imperative to perform a sensitivity study of the results for dierent mesh sizes.

It should not be overlooked that each of the unit cells of the mesh has geometrical properties to be taken into account in the subsequent numerical solution algorithm. Like many geometrical entities, each unit element has vertices or nodes, a center point or cell center, and area (2D) or volume (3D).

Numerical solver and discretization of the physical equations: The cases in which we can nd analytical solutions for the partial dierential equations involved in uid mechanics are very rare. In consequence, it is necessary to get the approximated solutions of the governing equations of uid dynamics. For this purpose, a numerical method called discretization has to be applied to the meshed domain created at the respective stage described above. This procedure consists of the adaptation of the uid equations for the unit cells (or control volumes as we are going to see soon) that make up the simulation grid. The result of this stage is a system of spatially and temporally discretized algebraic equations 2 . Three methods are often used for this purpose: nite element method (FEM), nite dierence method (FDM) and nite volume method (FVM). The FVM has been widely used for CFD because its formulation is based on the conservation of physical quantities Ψ in each of the unit mesh cells or control volumes (see equation I.1). This feature allows a more straightforward interpretation of 1 Obviously, the accuracy of the CFD results depends not only on the number of meshing elements but also (deeply) of the implemented model (e.g. RANS, LES, DNS). 2 The general conservation (transport) equations of mass, momentum, energy are the basic equations to be discretized into a system of algebraic equations. However, depending on the phenomena to be analyzed more equations have to be added to the procedure. the FVM results compared to the FEM and FDM approaches [START_REF] Versteeg | An introduction to computational uid dynamics : the nite volume method[END_REF].

     Rate of change over time of Ψ (in a cell volume)      =         N et rate of increase of Ψ due to convection         +         N et rate of increase of Ψ due to dif f usion         +      N et rate of production of Ψ      (I.1)
After the formulation of the discretized equations, we proceed to their resolution in every node of each mesh unit cell. To obtain the numerical non-linear solutions is necessary an iterative procedure that ensures a numerical convergence according to the numerical residuals that are considered appropriate in each case of study. Naturally, several criteria or conditions facilitate numerical convergence. One of the most notable is known as the Courant-Friedrichs-Lewy condition or simply Courant number. This quantity (valid for time explicit scheme) consists of the ratio between the time interval (time discretization or time step) ∆t and the residence or modelization time (= ∆l

U

) in a nite volume (associated with the size ∆l of a particular unit cell and the corresponding velocity U to be determined).

C = U ∆t ∆l ≤ C max (I.2)
The rst implication of the equation I.2 is that the time step ∆t in a numerical experiment must have a higher bound, otherwise, the simulation will produce incorrect results. Second, when the spacing between grid points is reduced, the upper bound for the time interval is also reduced. This condition could eventually make high computational resources necessary for the proper simulation of turbulent uids. Therefore the CFL condition is an indispensable aspect for the mesh design to be implemented in the CFD implementation.

Once an appropriate CFL number has been dened, the numerical convergence of the simulation will require a certain number of iterations. In this iterative stage, the following guidelines are typically considered:

The changes in the solution variables (e.g. numerical residuals) from one iteration to the next are negligible.

Overall properties conservation are achieved (mass imbalance, energy dissipation).

Monitoring quality key parameters such as torque, velocity proles, pressure protes or others associated with the studied phenomena.

Failure to comply with some of these aspects could imply the redesign of the mesh for the studied domain. In this sense, the implementation of CFD also has an iterative character as it usually happens in the designs or solutions to problems in most engineering elds.

The post-processing stage: The process ends with the presentation and visualization of results. For this purpose, most commercial codes such as ANSYS/FLUENT and COMSOL are provided with user-friendly tools. In the case of in-house codes, software such as Matlab, Python, or Paraview can be used. The most important aspect at this stage is to know the type of format of the results generated in the numerical simulation in order to choose an appropriate visualization tool.

I.1.2 Fluid mechanics equations (turbulence models)

In this section we present the uid ow equations of interest in the present work. Specically, we will provide the basic mathematical relations corresponding to a 3D Newtonian single-phase incompressible uid and to a 3D two-phase uid ow.

I.1.2.1 3D RANS equations and k -ε turbulence model

If one omits gravity forces the governing equations for the rst kind of uid ow are written as follows (Einstein summation convention applies to repeated indices):

∂U i ∂x i = 0 (I.3) ∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂p ∂x i + ν ∂ 2 U i ∂x j ∂x j (I.4)
with:

x i : x,y,z Cartesian directions.

U i : component of the total instantaneous velocity in the direction i.

p : presion.

ρ : density of the uid.

ν : kinematic viscosity of the uid.

Equation I.3 denotes the mass conservation and equation I.4 gives the momentum conservation (Newton's Second Law).

The numerical solution of equations I.3 and I.4 for turbulent ows is computationally very expensive because the entire spectrum of Kolmogorov 3 spatial and temporal scales involved (on the order of micro-meters and milli-seconds) must be considered in the corresponding simulation. To tackle this situation there are two options. The rst strategy is called direct numerical solution (DNS) and implies to use supercomputers whose computational capabilities allow the solution of the equations with the spatial and temporal resolutions already indicated. This option is sometimes not very practical, especially in the case of industrial volume domains. The second option is the implementation of turbulence models that make it feasible to obtain approximate numerical solutions. A variety of turbulence models are currently available [START_REF] Tennekes | A rst course in turbulence[END_REF][START_REF] Chassaing | Turbulence en mécanique des uides : analyse du phénomène en vue de sa modélisation à l'usage de l'ingénieur[END_REF][START_REF] Versteeg | An introduction to computational uid dynamics : the nite volume method[END_REF]. Each of these proposals is designed for specic objectives depending on the conditions and characteristics of the uid ow to be studied. For example we can mention the Large Eddy Simulation (LES) and the k -ϵ models. The LES is a very ecient approach, which allows computing the most energetic scales, while modellng only a part of the turbulence scales. The second model, the k -ϵ proposes to model all turbulent scales by introducing two additional variables: turbulent kinetic energy and turbulent dissipation. Below we will provide more details on the implementation of this model.

The most industrially relevant ows there are conned uids in which Reynolds shear stresses are of major relevance. It has been shown that in these conditions solving for mean ow and modeling all turbulent scales leads to results very close to those obtained in physical experiments [START_REF] Versteeg | An introduction to computational uid dynamics : the nite volume method[END_REF]. The simulation model we have just described is known as the RANS model because it is based on an averaged version of the Navier-Stokes ow equations [START_REF] Alfonsi | Reynolds-Averaged NavierStokes Equations for Turbulence Modeling[END_REF]. Consider that the instantaneous local velocity U i (⃗ r, t) can be written as the sum of an average local component Ūi (⃗ r) and a uctuating or random turbulent local component u ′ i (⃗ r, t) as indicated by the equation:

U i (⃗ r, t) = Ūi (⃗ r) + u ′ i (⃗ r, t) (I.5)
One can see that the decomposition equation I.5 is nothing more than a solution strategy for the Navier-Stokes equations. We must therefore include the velocity decomposition in I.3 and I.4 to obtain the corresponding equations to be solved. Such relations are known as Reynolds Averaged Navier-Stokes equations or RANS equations:

∂ Ūi ∂t + Ūj ∂ Ūi ∂x j = - 1 ρ ∂ p ∂x i + ν ∂ 2 Ūi ∂x j ∂x j - ∂u ′ i u ′ j x j (I.6)
The term τ ij = -ρu ′ i u ′ j is known as the Reynolds stress tensor and incorporates turbulent eects in the mean stress.

From a mathematical perspective, equations I.6 do not have a unique solution because there are more unknowns than equations ("non-closed" equations). To solve this inconvenience one proposes additional equations (known as closure models) to equalize the number of equations and unknowns in relation I.6. A typical closure model used in industrial applications for RANS equations is known as the k -ε model standard.

This model proposes to determine the components of the Reynolds tensor τ ij from the resolved mean velocity gradients (Boussinesq hypothesis). Additionally, this proposal ignores small-scale vortices (or eddies) in the motion. This leads to the calculation of a large-scale motion with an eective viscosity ν T or "eddy viscosity" or turbulent, which characterizes the transport and dissipation of energy in the smaller-scale ow. The mathematical relationships corresponding to these quantities are shown below.

τ ij = -ρu ′ i u ′ j = ρν T ∂ Ūi ∂x j + ∂ Ūj ∂x i - 2 3 ρδ ij u ′ i u ′ j (I.7) ν T = C µ k 2 ε (I.8)
where C µ is a dimensionless constant.

From I.7 and I.8 one observes that the k -ε model incorporates the formulation and solution of two additional equations of motion: the k equation of turbulent kinetic energy and the ε equation associated with the turbulent energy dissipation rate:

∂k ∂t + Ūi ∂k ∂x i = ∂ ν + ν T σ k ∂k ∂x i ∂x i -u ′ i u ′ j ∂ Ūi ∂x j -ε (I.9) ∂ε ∂t + Ūi ∂ϵ ∂x i = ∂ ν + ν T σε ∂ε ∂x i ∂x i -C ε1 ε k u ′ i u ′ j ∂ Ūi ∂x j -C ε2 ε 2 k (I.10)
where σ k ,σ ε , C ε1 ,C ε2 are dimensionless constants to be dened.

The values of the constants given in equations 1,2,3 are summarized in the table [START_REF] Launder | Lectures in mathematical models of turbulence[END_REF]): To conclude this section, we will present the model equations corresponding to the 3D two-phase uid ow used in our work. For this case, we will consider the Euler-Euler CFD approach for multiphase ows. In this model, the media involved are considered continuous and mathematically interconnected. Each grid element has a certain amount of gas and liquid. This consideration introduces the concept of phase volume fraction to describe the phase transport α L for liquid and α G for gas. These volume fractions of the liquid and gas phases are taken as functions of space and time, and naturally, their sum is equal to 1. As shown for the RANS equations, conservation laws must be ensured in this case for each phase. These physical requirements become a set of the main equations that must be complemented with additional turbulent and closure equations.

In this conceptual framework the laws of conservation of mass take the following form:

∂α k ρ k ∂t + ∂ α k ρ k U k,i ∂x j = 0 (I.11)
with k = l: liquid phase, or k = g: gas phase.

x i : x,y,z Cartesian directions.

The conservation of the momentum has the following structure:

∂α k ρ k ⃗ U k ∂t + ∇ • α k ρ k ⃗ U k ⃗ U k = -α k ∇ (p) + +α k ρ k ⃗ g + ∇ • (α k (Σ k + Rt k )) (I.12)
where p is the pressure.

g is the gravity acceleration.

Σ k is the molecular stress tensor.

Rt k is the Reynolds stress tensor. and all acting forces in bubble are described by the following expression:

⃗ M k = ⃗ M D k + ⃗ M AM k + ⃗ M L k + ⃗ M T D k (I.13)
The terms of I.13 represent the averaged drag, added mass, lift, and turbulent disspersion forces per unit volume. Now it is convenient to provide the mathematical expressions for each of these terms.

Drag force

⃗ M D g = -⃗ M D l = - 1 8 A i ρ l C D ⃗ U g -⃗ U l ⃗ U g -⃗ U l (I.14)
with A i : interfacial area concentration. C Di : drag coecients for bubbles.

This parameters can be determined experimentally.

Added mass force

⃗ M AM g = -⃗ M AM l = -C lg A 1 + 2α g 1 -α g α g ρ l ∂ ⃗ U g ∂t + ⃗ U g • ∇ ⃗ U g - ∂ ⃗ U l ∂t + ⃗ U l • ∇ ⃗ U l (I.15)
In I.15 C lg A is named added mass coecient and is equal to 0.5 for spherical bubble.

In addition, the factor 1+2αg 1-αg takes into account the eect of the bubbles concentration [START_REF] Zuber | On the dispersed two-phase ow in the laminar ow regime[END_REF][START_REF] Ishii | TWO-FLUID MODEL FOR TWO-PHASE FLOW[END_REF].

Lift force

⃗ M L g = -⃗ M L l = -C L α g ρ l ⃗ U g -⃗ U l ∧ ∇ ∧ ⃗ U l (I.16)
The C L is named the lift coecient and the ∧ symbol stands for vectorial product.

Turbulent dispersion force

According to [START_REF] Lance | PHASE DISTRIBUTION PHENOMENA AND WALL EFFECTS IN BUBBLY TWO-PHASE FLOWS[END_REF] one can express the turbulent dispersion contribution as:

⃗ M T D g = -⃗ M T D l = -C T D ρ l K l ∇α g (I.17)
In this expression K l is the liquid turbulent kinetic energy and C T D is a numerical constant of order 1.

Finally we have to propose an appropriate turbulence model (closure model) for the conservation equations already exposed. We begin this stage by indicating that the proposed denition of dynamic viscosity is taken as a reference the model of Lopez de Bertodano, which consists of a linear superposition according to the relationships:

µ = µ shear + µ BIT (I.18) µ BIT = 0.6α g ρ L D B |U R | (I.19)
The equations I.18 and I.19 represent a bubble induced turbulence approach (BIT).

These expressions are complemented by the closure relation which models the Reynolds tensor Rt ij for the continuous phase. According to this conceptual framework, this tensor has as its predominant component the dissipative turbulence produced by the mean velocity gradient and by the wakes of the bubbles. The corresponding model establishes:

(

1 -α l ) ∂Rt ij ∂t + 3 k=1 U l,k ∂Rt ij ∂x = ∂ ∂x k ρ l ν l + ρ l C s K ϵ Rt ij ∂ ∂x k ((1 -α L ) Rt ij ) + (1 -α l ) (P ij + G ij + Φ ij + ϵ ij ) (I.20)
where C s is a constant, Φ ij corresponds to the pressure uctuations, P ij the production of turbulent kinetic energy, G ij the velocity gradient, and ϵ ij the dissipation of turbulent kinetic energy.

The viscous destruction of the stress tensor can be modeled as follows:

ϵ ij = 2 3 ϵδ ij (I.21)
In addition, the turbulent diusion (of diusive nature) has a very well-known structure in terms of the generalized gradient diusion:

D t ij = ∂ ∂x k C s K ϵ u ′ k u ′ q ∂u ′ i u ′ j ∂x q (I.22)
We also consider that the pressure uctuations aect turbulent structures so that there is a redistribution of energy producing more isotropy. This can be express according to the next equations:

Φ ij = Φ ij,1 + Φ ij,2 + Φ ij,3 (I.23)
where

Φ ij,1 = -C 1 ϵ u ′ i u ′ j K - 2 3 δ ij (I.24) Φ ij,2 = -C 2 P ij - 2 3 P δ ij (I.25) Φ ij,3 = -C 3 G ij - 2 3 Gδ ij (I.26)
with:

C 1 = 1.8

P = 0.5 P kk C 2 = 0.6 G = 0.5 G kk C 3 = 0.55 I.1.2.

Fluid ow components or structures of interest for the present project

It is possible to show that any turbulent uid ow (valid for experimental or numerical data) is composed of a set of structures of dierent spatial and temporal scales and energy content. We will see later that this impressive feature will allow the implementation of numerical tools aimed at the identication of the mentioned uid components. Figure I.3 schematizes the typical groups of structures that can be identied in a uid ow according to the scale associated with each group (esa, 2018). The energy ow associated to turbulent ows is also presented in this gure (Richardson's energy cascade). First, the power provided or injected into the system (whether by mechanical, chemical, electrical, etc.) is invested in the creation of large structures known as coherent or organized structures. In this group are located components such as mean ow, large vortices (of the order of the dimensions of the whole domain under consideration, large circulation loops, traveling vortices, etc.). Then the energy is transmitted to the smaller structures (turbulente eddies) and scales that are associated with turbulent phenomena. It is in the latter domain that energy dissipation occurs. In the present research project, we will focus our attention on coherent structures specically on mean ow, traveling vortices, and a type of structure called macro-instabilities or MI's.

The mean ow is the component of the ow that does not change in time but can change in space.

As for the trailing vortices, it is important to mention some characteristics. They are rotating ow structures that also move in space and play an important role in mixing.

They are generally produced in the uid ow generated by rotating turbines (see Figure I.4 (adapted from [START_REF] Lee | Turbulence properties of the impeller stream of a Rushton turbine[END_REF])). The oscillations in the uid produced by Θ (⃗ r, t) = Θ 0 g (⃗ r) e iωt = Θ 0 g (⃗ r) (cos (ωt) + i sin (ωt))

(I.27)
where Θ 0 : is a constant related to the wave amplitude that can be a complex number.

g (⃗ r): denes the spatial part of the wave.

i = √ -1: is the unit imaginary number.

As will be seen in the section I.3.2 the mathematical structure of equation I.27 is not limited only to the case of trailing vortices. In fact, we will see that any uid ow can be represented approximately as a sum of ow components (whether periodic or not) every expression analogous to I.27.

The third type of organized structure that concerns us is the MI's. This kind of ow has been identied in physical experiments [START_REF] Yianneskis | An experimental study of the steady and unsteady ow characteristics of stirred reactors[END_REF][START_REF] Bruha | Experimental Study of Turbulent Macroinstabilities in an Agitated System with Axial High-Speed Impeller and with Radial Baes[END_REF][START_REF] Montes | Velocity eld macroinstabilities in an axially agitated mixing vessel[END_REF][START_REF] Hasal | Macro-instabilities of velocity eld in stirred vessel: detection and analysis[END_REF][START_REF] Micheletti | Precessional ow macro-instabilities in stirred vessels: study of variations in two locations through conditional phase-averaging and crosscorrelation approaches[END_REF][START_REF] Doulgerakis | On the Manifestation and nature of macroinstabilities in stirred vessels[END_REF] and numerical simulations (Hartmann et al., 2004b;[START_REF] Nurtono | Macro-instability characteristic in agitated tank based on ow visualization experiment and large eddy simulation[END_REF] of mixing devices widely used in industry 4 . This type of structures has a spatial conguration of whirlpool type, with frequencies that are sub-harmonic to the frequency associated with the rotating ows of the system (of the order of hundredths of Hz for our cases of interest).

Macro instabilities produce variations in mean ow that can signicantly change ow patterns, thus aecting mixing eciency [START_REF] Ducci | Decomposition of Flow Structures in Stirred Reactors and Implications for Mixing Enhancement[END_REF]. Therefore, critically im- portant processes such as mass and heat transport, local retention of gases or solids, homogenization, mechanical stress distribution [START_REF] Hasal | Force Eects of the Macro-Instability of Flow Pattern on Radial Baes in a Stirred Vessel With Pitched-Blade and Rushton Turbine Impellers[END_REF], and even turbulence levels [START_REF] Nikiforaki | On the origin, frequency and magnitude of macro-instabilities of the ows in stirred vessels[END_REF] can be aected due to the presence of MI's. Due to the aforementioned facts, the analysis and identication of these uid components are of utmost importance for the design of the mixing apparatus.

In the following section, we will address the particular aspects related to the implementation of CFD in the systems of interest for the present study.

I.2 Application of CFD for the study of reactor vessels

The CFD implementation of a specic system has unique particularities that depend on the phenomena involved as well as on the models associated with such phenomena. Due to the above, we will briey describe some aspects associated with the two systems that will be studied in this document: stirred tank reactors with radial ow impellers and bubble column reactors.

One of the most interesting and relevant applications of CFD is related to the numerical study of industrial devices known as reactor vessels (IAEA, 2020). These systems are extremely important for both practical and academic reasons. Reactor vessels are widely used in all types of the chemical industry. In addition, the complexity of the physical, chemical, and biological phenomena involved in these devices make them extremely attractive in the academic eld.

A particular type of reactor vessel that is dynamically activated by a turbine (with one or several impellers) is widely used in industrial chemical processes especially to favor mixing through axial ows. The modeling of these stirred tank elements can be developed following two dierent approaches: steady and unsteady state approaches.

The rst strategy consists of solving all the equations of the system in a steady-state, incorporating the eect of the impellers by imposing boundary conditions between the two zones of the domain. To achieve this there are some models such as the Impeller boundary condition approach (IBC), Source-sink approach (SS), Inner-outer approach (IO), and Multiple reference frame technique (MRF).

The second modeling category, the unsteady state approach, establishes a time-dependent interaction of the impellers with the uid ow of the system. This approach has dierent models such as Impeller modeling within lattice-Boltzmann-LES framework (adaptive force eld technique, AFT), Moving-deforming grid technique, and Sliding Mesh Approach (SM).

For most of these models, the studied domain is divided into two zones, a low inuence It is virtually impossible to apply this method for multi-phase systems.

Source-sink (SS) Modelling of the impeller rotation as a momentum, the baes represent sinks of the momentum. The impeller blades are discretized in number of vertical strips. The full vessel has to be modelled.

Inner-outer (IO) The domain is divided in two ovelapping regions: Inner region (impeller zone),

Outer region (bulk vessel zone). Inner zone simulations are used as conditions to simulated the Outer region.

Multiple reference frame (MRF) The domain is divided in two no-ovelapping regions:

Inner region (impeller zone),

Outer region (bulk vessel zone). Inner zone simulations (rotating frame) are used as conditions to simulated the Outer region. The uid is considered a system of particles obeying the conservation laws.

Each uid particle occupies a corner of the lattice.

As time evolves each particle moves and collides with other particles. Solid walls provide control points. Tangential velocities of the impeller surface matches with the physical impeller rotation.

It is no necessary to mesh impeller and vessel geometries separately.

Implies enormous computational requirements.

Moving-deforming grid

The whole domain has one unique mesh.

Mesh quality control is challenging. Mesh impeller deform the mesh vessel when moving. Mesh motion is determined by time derivative transformations (Eulerian equations to Lagrangian equations).

Highest computation requirements of the unsteady methods.

Sliding Mesh (SM) Full transient simulations are carried out using two cylindrical mesh zones:

rotative zone (impeller rotative speed) and xed zone (inertial frame).

The rotative zone includes the rotating turbine, the xed zone the rest of the vessel volume.

Each zone is meshed separately. Geometry-mesh turbine-impellers can be modeled with a high quality as solid rotative walls.

The ow equations are formally written in the rotative zone.

The two mesh zones are implicity coupled using a sliding-grid algorithm, taking into account the relative motion and the required interpolations are performed.

During the simulation the volume mesh cells does not change over time.

This method can be used as a design tool to assess dierent congurations.

A large computational expense is required.

As we will see in chapter IV our working numerical database for the baed stirred tank was generated using the unsteady-state SM model.

The second type of vessel we are concerned with is known as a bubble column. The relevance of this device lies in the fact that it operates with a two-phase ow that is present The interaction of the rising gas and the surrounding liquid results in a type of ow uid known as bubble swarms or plumes. The type of ow structure generated depends essentially on parameters such as the velocity of the injected gas, the viscosity of the liquid, and the aspect ratio (height H-to-width W ratio). Figure I.8 provides a picture that allows visualizing dierent types of ow according to two dierent gas injection velocities (2.4 mm/s and 21.3 mm/s respectively) reported by [START_REF] Cachaza | Unied study of ow regimes and gas holdup in the presence of positive and negative surfactants in a non-uniformly aerated bubble column[END_REF].

The Euler-Euler approach presented in the section I.1.2.2 is reported in some papers [START_REF] Becker | Gasliquid ow in bubble columns and loop reactors: Part II. Comparison of detailed experiments and ow simulations[END_REF]Peger et al., 1999;[START_REF] Olmos | Numerical simulation of multiphase ow in bubble column reactors. In uence of bubble coalescence and break-up[END_REF][START_REF] Olmos | Numerical description of ow regime transitions in bubble column reactors by a multiple gas phase model[END_REF][START_REF] Fard | Transient large-scale two-phase ow structures in a 3D bubble column reactor[END_REF] Similar to the case of stirred tanks, the implementation of CFD for these devices involves the consideration of additional aspects to those described in section I.1.2.2.

Particularly for the case study of the present investigation, additional approaches included in the NEPTUNE CFD software were considered, which allow modeling the forces of drag, added mass, lift, and turbulent dispersion. Specic details are provided in the section II.2.

I.3 Data Driven Methods and Reduction Order Methods (ROM's) for ow elds

The last few decades have been characterized by unprecedented progress in experimental and numerical simulation techniques, leading to the generation of highly accurate digitized data. The processing requirements of these data often exceed the available computational resources. In this context, the implementation of algorithms or methodologies that allow the extraction, analysis, and classication of information becomes a real necessity. This scenario explains why the attention of a large part of the scientic community has been directed towards this type of numerical methodologies for data processing which helps to exploit the immense potential of large databases [START_REF] Hey | The Fourth Paradigm: Data-Intensive Scientic Discovery[END_REF][START_REF] Montáns | Data-driven modeling and learning in science and engineering[END_REF][START_REF] Mendez | Data-Driven Fluid Mechanics: Combining First Principles and Machine Learning[END_REF][START_REF] Karniadakis | The rise of data-driven modelling[END_REF].

As the family of data-driven methodologies is vast, they can be applied in elds as diverse as fundamental physics, materials science, bio-engineering, chemical industry, petroleum industry, pharmaceutical industry, and, of course, uid mechanics. Our interest is focused on the latter, specically concerning reducing the number of data inputs (velocity elds) to a manageable size while preserving the most relevant information in the data in terms of energy or variance. In purely technical terms, reduction order methods (ROMs) belong to the group of disciplines known as unsupervised Machine Learning [START_REF] Ding | Adaptive dimension reduction for clustering high dimensional data[END_REF][START_REF] Chinesta | Foreword[END_REF][START_REF] Brunton | Machine Learning for Fluid Mechanics[END_REF][START_REF] Delua | Supervised vs. Unsupervised Learning: What's the Dierence? SME, IBM Analytics[END_REF][START_REF] Vinuesa | The Potential of Machine Learning to Enhance Computational Fluid Dynamics[END_REF]. These algorithms allow the identication of patterns present in the databases studied even when this input information is not associated with a particular model. This feature provides data-driven methods with immense versatility.

Within the group of order reduction algorithms, we nd some techniques that have been widely used for CFD and PIV known as Proper Orthogonal Decomposition (POD) and Dynamical Mode Decomposition (DMD). The conceptual basis of these methodologies is the decomposition of the input data in terms of modes. These results are obtained by solving eigenvalue problems through matrix operations such as singular value decomposition (SVD) or QR decomposition.

In the present work, we will use the POD and DMD modal decompositions as order reduction tools for velocity eld data from CFD numerical simulations. Below we will focus on the description of POD and DMD techniques.

I.3.1 The fundamentals of POD I.3.1.1 Principle of the method

The proper orthogonal decomposition (POD) or Kosambi-Karhunen-Loève decomposition [START_REF] Kosambi | Statistics in function space[END_REF][START_REF] Loève | Functions aleatoire de second ordre[END_REF][START_REF] Karhunen | Zur spektral theorie stochastischer prozesse[END_REF][START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows[END_REF][START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF] is a practical procedure for extracting the most energetic components from a set of functions (in the continuous case) or experimental (physical or numerical) data (in the discrete case). The key idea of the technique is to provide a vector space basis with the special condition that the mean square projection of the data set in this vector space is maximal.

The vectors constituting the resulting POD basis become modes that allow reconstructing approximately the information contained in the input database. The supreme advantage of this is that the modes or elements that are statistically most important can be identied and retained. Under these conditions, we can reduce a large number of variables to a much smaller set of spatially uncorrelated variables, while maintaining the spatial and temporal resolution of the original set of variables [START_REF] Liang | PROPER ORTHOGONAL DECOMPOSITION AND ITS APPLICATIONSPART I: THEORY[END_REF][START_REF] Kerschen | The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems: An Overview[END_REF]. When applied to uid ow analysis, the data set can be obtained from a series of experimental or simulated velocity elds.

To begin the mathematical description of the method let us note that the POD reconstruction of the p modes (nite truncations in the modal expansion) is optimal compared to any other possible orthogonal basis having the same dimension p [START_REF] Liang | PROPER ORTHOGONAL DECOMPOSITION AND ITS APPLICATIONSPART I: THEORY[END_REF][START_REF] Aubry | On the hidden beauty of the proper orthogonal decomposition[END_REF]. This surprising feature is because this orthogonal decomposition minimizes the mean square of the associated error resulting from the reconstruction procedure. If the data set consists of a collection of instantaneous three-dimensional velocity elds, the original eld ⃗ U (⃗ r, t) is approximated as a linear combination of p modes, ⃗ ϕ k (⃗ r) weighted by time varying coecients (modal coecients) a k as follows:

⃗ U P OD (⃗ r, t, p) = p k=1 a k ⃗ ϕ k (⃗ r) (I.28)
In summary, POD allows minimizing the amount of information needed to represent statistically dependent data.

Therefore, POD is an extremely useful tool that nds multiple applications in CFD and other areas of computational processing when it is desired to generate a reduced-order description without compromising the most relevant information [START_REF] Lu | Review for order reduction based on proper orthogonal decomposition and outlooks of applications in mechanical systems[END_REF]. Moreover, the characteristic of the POD to generate a reduction that separates the spatial and temporal parts opens the possibility of understanding phenomena of immense complexity in a simpler manner.

As with any approximation method, POD also has some disadvantages. These drawbacks are based on the fact that this procedure is based on a second-order correlation (or energy), so the higher-order correlations are ignored [START_REF] Aubry | On the hidden beauty of the proper orthogonal decomposition[END_REF][START_REF] Taira | Modal Analysis of Fluid Flows: An Overview[END_REF].

For example, [START_REF] Graham | Alternative approaches to the Karhunen-Loève decomposition for model reduction and data analysis[END_REF] points out that the analysis downplays the importance of infrequent events or structures, even though they may be dynamically very important (burst-like phenomena). In the section III.2.2, we will show in detail a case that exemplies this situation. Moreover, [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] pointed out, the variance (or energy) ranking might not be in some cases the most appropriate parameter for the choice of dominant elements in the studied databases. This last aspect is particularly relevant when seeking to identify structures that are monochromatic, i.e. associated with a specic frequency.

The elements required to construct the linear combination of equation I.28 can be determined from two dierent approaches: the direct method and the snapshot method. The choice of the approach to be used depends on the spatial and temporal resolution in the studied database. If one has a spatial resolution that is much higher than the temporal one, the snapshot method is preferable, as it involves lower computational requirements than the direct approach (Sirovich, 1987a;Sirovich, 1987b;Sirovich, 1987c;[START_REF] Smith | Low-Dimensional Modelling of Turbulence Using the Proper Orthogonal Decomposition: A Tutorial[END_REF]. The following section shows the fundamental equations of the orthogonal decomposition.

I.3.1.2 Mathematical details: Direct method and Snapshot method

Let us sample for N instants (each one separated by a time of ∆t) a 3D velocity vector

eld ⃗ U (⃗ r) = (u (⃗ r) , v (⃗ r) , w (⃗ r))
, composed of L points in the space of a dened domain.

Each of the vector elds in the resulting database is called a snapshot since it represents a vector "photograph" of the system under study. The data set can be arranged or structured in a matrix M with 3L lines and N columns (called the snapshot matrix) as follows:

M =                        u (1) 1 • • • u (N ) 1 . . . • • • . . . u (1) L • • • u (N ) L v (1) 1 • • • v (N ) 1 . . . • • • . . . v (1) L • • • v (N ) L w (1) 1 • • • w (N ) 1 . . . • • • . . . w (1) L • • • w (N ) L                        = ⃗ U (1) , ⃗ U (2) , • • • , ⃗ U (N ) (I.29)
The k th column of the relation I.29 represents the velocity eld sampled at the k th instant and ⃗ U (k) is a column vector of 3L lines.

From I.29 the covariance or correlation tensor R can be constructed following the direct or snapshot method as follows:

R =          1 N M M T Direct M ethod 3L × 3L 1 N M T M Snapshot M ethod N × N (I.30)
The correlation tensorR denes the magnitude and manner in which the velocities at dierent points are correlated. Consider the velocity vector ⃗ U (⃗ r i ) at a i th point and the velocity ⃗ U (⃗ r j ) at a j th point in the uid domain. If the velocities ⃗ U (⃗ r i ) and ⃗ U (⃗ r j ) are statistically dependent their correlation will be non-zero, otherwise, the correlation will be zero. In other words,R represents a two-point spatial correlation averaged over time.

If one uses the direct method, the correlation tensor is represented by a 3L × 3L matrix, where L represents the number of data points. This mathematical process causes the computational cost to increase with the square power of the number of points in the domain under study (the size of the tensorR is (3L) 2 ). For this reason, when one has a domain with a large number of points, it is necessary to evaluate whether one has the minimum computational and storage requirements associated with POD matrix operations.

We can also use the snapshot approach in lieu of the direct method. The snapshot method proposed by (Sirovich, 1987a;[START_REF] Smith | Low-Dimensional Modelling of Turbulence Using the Proper Orthogonal Decomposition: A Tutorial[END_REF] is always valid when the ergodic hypothesis is fullled (all accessible microstates in the phase space have the same probability to occur over a long period). As can be deduced from the denition I.29, the implementation of this approach involves the construction of an R tensor with dimensions N × N . In this case the computational cost less prohibitive compared with the Direct method (the size of the tensor R is (N ) 2 ). Due to this advantage, the snapshot method is used in this work.

The implementation of the snapshot method involves solving an eigenvalue problem to determine the lamda eigenvalues and ksi eigenvectors that diagonalize the R correlation tensor:

R⃗ γ k = λ k ⃗ γ k (I.31)
It is pertinent to recall that the correlation tensor R is Hermitian so the rst N (N being the number of instants analyzed) eigenvalues are identical for the direct approach and the snapshot approach. However, the eigenvectors are dierent for each method.

Once one has calculated the eigenvectors ⃗ γ k associated with the snapshot technique it is necessary to determine the vectors ⃗ ϕ k with the appropriate 3L dimensions. The necessary relationship for this procedure is as follows:

⃗ ϕ k = M⃗ γ k 1 λ k , k = 1, ..., N (I.32)
Since the eigenvectors constitute an orthonormal basis, they obey the following relationship:

⃗ ϕ i T ⃗ ϕ j = δ ij , i, j = 1, ..., N (I.33)
After the calculation of the eigenvector, ⃗ ϕ k the modal components a k (t) are obtained as:

a k (t) = (M ) T ⃗ ϕ k , k = 1, ..., N (I.34)
This completes all the elements necessary to reconstruct the velocity eld according to I.28.

At this point, it is important to emphasize one aspect regarding the denition of the correlation tensor R. Whichever the approach one implements (direct or snapshot method), the main objective is to nd a basis of orthonormal ⃗ ϕ k (⃗ r) vectors (3L-I.3 Data Driven Methods and Reduction Order Methods (ROM's) for ow elds dimension)

5 and modal coecients a k (N -dimension) that allow us to construct the vectorial space associated to the matrix M . Note that we are looking for a set of vector elds that minimize the mean squared deviation between the POD approximation and the values contained in M .

I.3.1.3 Examples of POD as a ROM's tool for reactor vessels

As already mentioned in I.2, stirred tank reactors are widely used in chemical and biochemical industrial applications [START_REF] Nienow | Stirring and Stirred-Tank Reactors[END_REF]. The hydrodynamics of such systems has been the subject of numerous experimental [START_REF] Ng | Observations on the Distribution of Energy Dissipation in Stirred Vessels[END_REF][START_REF] Bugay | Experimental analysis of hydrodynamics in axially agitated tank[END_REF][START_REF] Baldi | Determination of Dissipation Rate in Stirred Vessels Through Direct Measurement of Fluctuating Velocity Gradients[END_REF][START_REF] Escudié | Experimental analysis of hydrodynamics in a radially agitated tank[END_REF][START_REF] Ducci | Vortex tracking and mixing enhancement in stirred processes[END_REF] and CFD studies that revealed its complexity (Hartmann et al., 2004a;[START_REF] Delafosse | LES and URANS simulations of hydrodynamics in mixing tank: Comparison to PIV experiments[END_REF][START_REF] Delafosse | Trailing vortices generated by a Rushton turbine: Assessment of URANS and large Eddy simulations[END_REF]. Beside the hydrodynamic aspects, additional models describing the chemical or biochemical reactions have to be implemented in the CFD software in order to couple mixing and reaction and perform predictive simulations [START_REF] Olmos | Eects of bioreactor hydrodynamics on the physiology of Streptomyces[END_REF]. In some applications, it is important to consider the time varying ow eld but the computational cost for solving simultaneously momentum and chemical equations is often prohibitive.

In such situations, it is advantageous to nd a way to access the spatio temporal variation for the velocity eld without solving the Reynolds Averaged Navier-Stokes equations (RANS) throughout the process duration. Since the ow in a stirred tank is periodic and organized, the reduced order methods, such POD or DMD, can be used to achieve this objective.

Due to its elegance and practicality, POD has been used to identify, study and model the dynamics of large-scale average spatial structures of the uctuating velocity eld for turbulent ow (Sirovich, 1987a;Sirovich, 1987b;Sirovich, 1987c;[START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows[END_REF][START_REF] Borée | Extended proper orthogonal decomposition: a tool to analyse correlated events in turbulent ows[END_REF][START_REF] Joshi | Dynamics of Flow Structures and Transport Phenomena, 1. Experimental and Numerical Techniques for Identication and Energy Content of Flow Structures[END_REF][START_REF] Tirunagari | Analysis of Proper Orthogonal Decomposition and Dynamic Mode Decomposition on LES of Subsonic Jets[END_REF][START_REF] Du | POD reduced-order unstructured mesh modeling applied to 2D and 3D uid ow[END_REF][START_REF] El-Adawy | On the Application of Proper Orthogonal Decomposition (POD) for In-Cylinder Flow Analysis[END_REF]. In terms of kinetic energy, these large structures (named coherent structures) dominate in the uctuating ow, it follows that their identication is crucial in the study of the uid dynamics. This feature evidences the relevance and convenience of the POD to identify such organized structures, turning the method into an extremely attractive 5 Typically the 3L-dimensional space refers to the Cartesian coordinates of the points at which the velocity vectors are sampled. However, due to the versatility associated with the POD method, we do not rule out the possibility of its implementation in a space of a dierent nature than the positioning of points in space. In chapter III we will show with numerical arguments the feasibility of this reinterpretation in the denition of the space associated with the ⃗ ϕ k vector basis.

numerical technique [START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF].

The POD methodology has also been shown to be eective for studying instabilities [START_REF] Hasal | Macro-instabilities of velocity eld in stirred vessel: detection and analysis[END_REF][START_REF] Hasal | Force Eects of the Macro-Instability of Flow Pattern on Radial Baes in a Stirred Vessel With Pitched-Blade and Rushton Turbine Impellers[END_REF] and the transition to the turbulent regime in uids [START_REF] Knight | Kolmogorov inertial range for inhomogeneous turbulent ows[END_REF][START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF]. In [START_REF] Arányi | Analysis of dierent POD methods for PIV-measurements in complex unsteady ows[END_REF] the technique was also used to study complex ows. This technique has also proved to be very useful for the analysis of turbulent ow elds in mixing tanks [START_REF] Raju | Reynolds number scaling of ow in a stirred tank with Rushton turbine. Part II Eigen decomposition of uctuation[END_REF] The manner in which energy was transferred between the POD modes was also studied. Also it is reported in the literature the application of the POD on CFD data for the study of turbulent swirling ow in an axisymmetric sudden pipe expansion [START_REF] Howard | Proper Orthogonal Decomposition (POD) analysis of CFD data for ow in an axisymmetric sudden expansion[END_REF].

The order reduction technique used proved to be quite robust in reducing the transient database, showing its potential to provide valuable insights into the ow structure. By using a large array of pressure and velocity data, POD was able to pick out several key ow features, including the movement of vortices, and the structure and period of the precessing ow. This was followed by the reconstruction and assessment of the ow eld in a stirred tank by applying POD to a 2D domain using CFD and PIV data (de Lamotte et al., 2018a; de Lamotte et al., 2018b). As reported in this work, the reduced order method was able to reconstruct the uid ow into dierent structures: mean ow, coherent structures, and turbulent ow. Also [START_REF] Janiga | Large-eddy simulation and 3D proper orthogonal decomposition of the hydrodynamics in a stirred tank[END_REF] applied the POD algorithm for the analysis of coherent structures and macro-instability in a 3D (3-dimensional), 3C

(3-component velocity) LES (Large Eddy Simulation) simulation for an unbaed stirred tank. In this numerical experiment the sliding mesh technique was implemented, thus the simulation domain was divided into a xed part and a rotating part. A similar work was recently published by [START_REF] Mikhaylov | Reconstruction of large-scale ow structures in a stirred tank from limited sensor data[END_REF]. The data was created through a Direct Numerical Simulation (DNS) with a frozen rotor approach and Reynolds numbers of 500, 600 and 700. In this research, POD was useful to reconstruct the temporal evolution of large-scale organized vortical structures behind the blades of a Rushton impeller in a non-baed stirred tank. It was found that the rst two modes dominate the energy spectrum by carrying 90 % of the total kinetic energy. It also is relevant to mention that some higher modes (3 rd , 4 th and 5 th , 6 th ) came in pairs in the energy ranking. These two latest remarkable works show the feasibility of the numerical recipe for order reduction in very complex ow elds generated by numerical simulations. However, some important points were not addressed. For example in [START_REF] Janiga | Large-eddy simulation and 3D proper orthogonal decomposition of the hydrodynamics in a stirred tank[END_REF], the POD decomposition was applied to each part of the mesh separately but a reconstruction of the 3D, 3C velocity eld for the bulk was not performed. In general, it is neither simple nor straightforward to determine the relationship and similarity between the POD method for separate parts of the domain and the POD method corresponding to taking the domain as a whole. This is due to the essentially statistical nature of the method. Additionally, an analysis of the eciency of the POD method in terms of calculation time, global and localized error was not provided. 4 ), using CFD data computed with a sliding mesh approach has not been considered yet.

Among the objectives of our work is to address the indicated aspects not developed in these last two mentioned works.

According to the literature consulted, the assessment of bubbly ows (experimental or numerical) using POD or DMD has not been deeply addressed. In the following, we will mention the most relevant ones, which will serve to give a context of the aspects treated in this type of analysis.

I.3.2 The fundamentals of DMD

Similar to POD, the DMD methodology allows identifying, separating, and classifying the dominant structures of a uid ow. This tool allows decomposing time-resolved data into modes, taking as a criterion the characteristic oscillation frequency f j DM D (or f j ), growth/decay rate σ j DM D (or σ j ) of each DMD mode Φ j . In addition, the dynamic evolution will be given by the time coecient Θ j (t) which allows proposing a reconstruction expression analogous to I.28:

⃗ U DM D (⃗ r, t, p) = p j=1 Θ j (t) Φ j (⃗ r) (I.35)
We will show the relationship between these DMD parameters in the next section.

The DMD procedure consists of the decomposition of a linear best-t operator A that approximates the existing dynamics in the data. The results obtained allow describing the physical mechanisms immersed in the input data and then signicantly reducing the degrees of freedom of dynamic systems. It is important to emphasize that unlike POD modes the DMD modes Φ j (⃗ r) of I.35 are not orthogonal. This characteristic will imply certain conditions in the way one can perform an energy ranking of such modes. We will discuss this topic at the end of this section.

Consider again the arrangement of ) of the equation I.29, but assume that each element ⃗ U (j) represents a vector in the R 3L space. It is proposed to nd an operator A that allows establishing the following linear model:

N chronologically ordered elements ⃗ U (1) , ⃗ U (2) , • • • , ⃗ U (N
⃗ U (k+1) = A ⃗ U (k) (I.36)
where A is an operador in the C 3L×3L space.The transformation in I.36 represents a linear tangent approximation when applied to nonlinear dynamics.

If I.36 is applied to all ⃗ U (j) elements a sequence relationship can be expressed as follows:

⃗ U (2) , ⃗ U (3) , • • • , ⃗ U (N ) = A ⃗ U (1) , ⃗ U (2) , • • • , ⃗ U (N -1) (I.37)
or in a more compact form:

M (N ) (2) = A M (N -1) (1) (I.38)
From the expression, it follows that the eigenvectors and eigenvalues of A allow establishing the approximation of the linear model sought. However, when L is very large the resolution of the corresponding eingenvalue problem implies a huge computational expense. To reduce these requirements, two algorithms known as companion matrix and similar matrix or Singular Value Decomposition procedure (SVD) can be implemented [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF][START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF][START_REF] Taira | Modal Analysis of Fluid Flows: An Overview[END_REF]. In the following, we will give a brief description of each numerical recipe.

I.3.2.1 Mathematical details: Companion matrix and Similar matrix (SVD)

Companion matrix

An alternative to reducing the calculation requirements of I.37 is through matrix equivalence:

M (N ) (2) = A M (N -1) (1) = M (N -1) (1)
S + rẽ N -1

(I.39)
where S is a matrix in the C N -1×N -1 space named companion matrix of A, r is a residual vector in R 3L , and ẽN-1 is an Euclidean vector of size N -1. When N is much smaller than L, the eigenvalue problem associated with S is more convenient than the corresponding one for A.

Clearly, the eigenvalues and eigenvectors of S are an approximation of those corresponding to A. The larger N the better the respective approximation.

We can show that S can be determined by a least-square procedure like qr decomposition of M

(N -1) (1) : M (N ) (2) ≈ M (N -1) (1) S = q r S ⇒ S = (r) -1 q ⋆ M (N ) (2) (I.40)
where q ⋆ is the complex conjugate of q.

By determining the eigenvalues λ j DM D and eigenvectors ⃗ Υ j DM D of S one obtains the following DMD parameters:

ℜ j = M (N -1) (1) ⃗ Υ j DM D (I.41) f j DM D = arg λ j DM D 2π∆t (I.42) σ j DM D = ln λ j DM D ∆t (I.43)
The vectorial quantity ℜ j is named Ritz vector and is associated to a single DMD frequency f j DM D and to a growth rate σ j DM D . To calculate the time coecients Θ j (t) we must rst dene the Gram matrix G of the Ritz vectors:

G = (ℜ) ⋆ ℜ (I.44)
where ℜ = (ℜ 1 , ℜ 2 , ...., ℜ N -1 ) and the (ℜ) ⋆ is the complex conjugate.

Then we can determine the relationship:

Θ = (G) -1 (ℜ) ⋆ M (N ) (1) (I.45)
The expression I.45 results in a matrix (of size N -1 × N ) in which the line j is the time coecient Θ j (t).

Finally we can calculate the DMD vectors dened in equation I.35:

Φ j (⃗ r) = ℜ j Θ j (t = t 0 ) (I.46)
Similar matrix (SVD)

The Similar matrix or simply SVD procedure proposes that the matrix A can be approx- imated by using the SVD decomposition of the M (N -1) (1) matrix as shown below [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF][START_REF] Jovanovi¢ | Sparsitypromoting dynamic mode decomposition[END_REF].

We know that we can express M

(N -1)

(1) according to the SVD decomposition:

M (N -1) (1) = U Σ (V) † (I.47)
One observes that U and V (the symbol † corresponds to the matricial transposition) are unitary, the columns of each of them form a set of orthonormal vectors or basis vectors.

Moreover, Σ is a diagonal matrix whose non-zero elements correspond to the eigenvalues of M

(N -1) (1) .

Now the expression of I.47 can be replaced in I.39 which leads to a new form of S:

S = (U) † M (N ) (2) V (Σ) -1 = (U) † A U (I.48)
We can see that the relation I.48 implies a projection of A into the eigenvector space of

M (N -1) (1)
.

Following the same logic established for the matrix companion method it follows to nd the eigenvalues and eigenvectors of S.

S Y j = m j Y j (I.49)
From the information provided by I.49 we can determine the following parameters:

The new Ritz vectors:

ℜ j = U Y j (I.50)
The frequency and growth rate of each DMD mode.

f j DM D = arg (m j ) 2π∆t (I.51) σ j DM D = ln (|m j |) ∆t (I.52)
Needless to say that both DMD methods are equivalent. However, the SVD algorithm may be more convenient if the matrix (G) -1 in I.45 is close to singular o badly scaled.

Finally, it is necessary to establish a criterion to order the DMD modes according to their energy content. Because the DMD modes are not orthogonal, no single way to perform energy ranking can be determined. Given this situation, the literature proposes two possible criteria [START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF][START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF].

The rst ranking method denes the ranking according to the amplitude of the DMD modes Φ j (⃗ r). However, this criterion does not consider that highly damped DMD modes have a negligible contribution, regardless of their long-term amplitude.

An energetic ranking method that contemplates the dynamic eect of the contribution of each mode is, therefore, more suitable in the presence of modes with a very negative growth rate σ j . [START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF][START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF] proposes to use the following relation:

E j = 1 τ τ 0 Φ j 2 λ j DM D t ∆t = Φ j 2 e 2σ j DM D τ -1 2σ j DM D τ (I.53)
with τ = total analysis duration and ∆t is the time step used in the analysis.

The two criteria mentioned above will be used in the sections III and IV to classify the DMD modes in terms of energy content.

I.3.2.2 Examples of DMD as a ROMs tool for reactor vessels

The eectiveness of this technique has been shown in the study of the dynamic information of ow elds obtained by numerical simulations and by measurements in physical experiments. For example, [START_REF] Semeraro | Analysis of timeresolved PIV measurements of a conned turbulent jet using POD and Koopman modes[END_REF] implemented this numerical tool to analyze 2D PIV data measured in a planar and conned jet with co-ow. The results allowed identifying important dynamical features as the large-scale oscillation of the inner jet and shear-layer oscillations.

In [START_REF] Schmid | Decomposition of timeresolved tomographic PIV[END_REF] DMD made it possible to identify relevant uid ow frequencies from PIV tomographic data performed on a transition water jet. This information led to a low-dimensional description of the ow eld. Some works have implemented DMD to assess the dynamics of very complex ow uids in the human cardiovascular system. For example [START_REF] Habibi | Data-Driven Pulsatile Blood Flow Physics with Dynamic Mode Decomposition[END_REF] study sick patientspecic blood ow generated from CFD. The results showed that the ROMs method represented accurately the analytical solutions for incompressible pulsatile ow in tubes.

In [START_REF] Arzani | Data-driven cardiovascular ow modelling: examples and opportunities[END_REF] a review of examples and applications of ROMs in cardiovascular ow uids was reported. Others autors [START_REF] Yu | Application of Dynamic Mode Decomposition to Study Temporal Flow Behavior in a Saccular Aneurysm[END_REF]) also report the use of DMD to identify the impact of inow conditions on spatio-temporal ow behavior in an cerebral aneurysm ow uid. These researchers perfomed low-order ow reconstruction from PIV data that allowed to understand the complexity of the ow uid investigated.

The article [START_REF] Le | Dynamic modes of inow jet in brain aneurysms[END_REF]) also performed assessment but in this case the DMD reconstruction was compared with in vivo data.

According to the literature consulted, the assessment of bubbly ows (experimental or numerical) using POD or DMD has not been deeply addressed. In the following, we will mention some of the most recent works, which will serve to give a context of the aspects treated in this type of analysis.

In [START_REF] Tabib | Analysis of dominant ow structures and their ow dynamics in chemical process equipment using snapshot proper orthogonal decomposition technique[END_REF] POD was implemented to reveal the dominant ow components and their dynamics in various industrial equipment including a cylindrical bubble column system. The data were generated from PIV measurements and using LES simulations. In the case of the bubble column, the authors implemented a LES simulation to generate the database. The researchers report that it was possible to determine the variation in plume oscillation and ow structures in the vortical region of the bubble column resulting from changes in sparger design. Furthermore, it was reported that a reconstruction of 50 POD modes were necessary to capture 88% of kinetic energy and to observe in a pronounced way the eects of plume oscillations on ow patterns of the liquid phase.

For his part, [START_REF] Kim | Dynamic structures of bubble-driven liquid ows in a cylindrical tank[END_REF] investigated the spatial and temporal structures of turbulent water ows driven by air bubbles in a cylindrical tank. The corresponding data were generated by means of PIV. The dynamics of the ow was analyzed using POD.

This paper showed that the POD technique was able to capture three dominant dynamic structures: recirculating ow, bubble-induced ow, and free surface oscillation for ows with dierent Reynolds numbers (from 8,300 to 21,100). [START_REF] Pang | Experimental investigation on the turbulence channel ow laden with small bubbles by PIV[END_REF]) investigated bubbly ows from data measured with the PIV technique in a prismatic vessel. The authors considered two parts: the pure liquid part and the liquid phase (phase discrimination methods were used to perform this separation procedure). The POD technique was then applied to the liquid phase and to the pure liquid part databases to determine the inuence of bubbles on the turbulence structure. Specifically, POD analyses allowed to establish that the injection of bubbles (present in the liquid phase) increases the energetic contribution of the large-scale turbulence structures.

Likewise bubbly ows was assessed in [START_REF] Aliyu | Characteristics of bubble-induced liquid ows in a rectangular tank[END_REF] through PIV measurements in a square tank. POD analysis on the liquid velocity data led to the identication of dominant structures according to the corresponding energy criteria. From these ndings, it was possible to ascertain that vortices induced in the liquid phase are more evident at the locations away from the stream. Moreover, the POD allowed to establish that strong streamwise ow is associated with the rst POD mode. This mode dominates (over bubble-induced motions) since it had the highest contribution to the turbulent kinetic energy as well as the fewest large-scale vortical structures. Finally, the higher modes showed more vortical structures including smaller scale and lower contribution to the overall kinetic energy.

More recently [START_REF] Laupsien | The organized ow structure of an oscillating bubble plume[END_REF] investigated an oscillating bubble plume created in a quasi-two-dimensional bubble column, three dierent gas ow rates of 50, 100, and 200 l/h were used in the investigation. The POD allowed the description of the overall hydrodynamics of the system. Specically, it was found that the periodic motion structures associated with the large-scale motions (vortices, plume oscillations) were successfully extracted by the orthogonal decomposition.

There seems to be even less published work on implementing DMD in vertical bubbly ows. According to the authors in [START_REF] Klevs | Dynamic mode decomposition of magnetohydrodynamic bubble chain ow in a rectangular vessel[END_REF], their work is the rst to perform an application of DMD to bubble ow with resolved dynamic liquid/gas boundaries. The researchers simulated a system in which horizontal and vertical arrays or "chains" of bubbles ascend through a rectangular liquid vessel. Liquid gallium and argon gas are used as working uids. Additionally, a constant static magnetic horizontal force is applied using permanent magnets. Dierent gas ows with and without magnetic force are simulated.

In this investigation, DMD allows making visible various momentum transfer and bubble interaction mechanisms. Moreover, the mode analysis makes possible the explanation of the observed ow patterns. This work is an excellent example of the applicability of DMD in the study and analysis of highly complex uid ows.

Chapter summary and conclusions

This chapter presents the basics of Computational Fluid Dynamics (CFD): conservation equations, fundamental steps, advantages, and disadvantages of the methodology. Such exposition provides a general overview of the CFD that will allow a better understanding of the content corresponding to the next chapters of this document. One nds that CFD is a widely used tool for the study of ows of academic and industrial interest. The implementation of these algorithms implies the realization of several stages related to:

the phenomenon under study, the creation of the geometry and the discretization of the domain to be analyzed, the formulation and the numerical solution of the ow equations, the evaluation of the numerical residuals and the post-processing of the nal solution of the ow equations.

Next, we present the uid ow equations of interest for this research. In the rst place, the equations for the RANS-k -ε turbulence model are shown. This approach solves the mean ow and models all scales of turbulence. To achieve this, one invokes the Boussinesq hypothesis which proposes that the momentum transfer caused by turbulent eddies can be modeled with an eddy viscosity. This aspect implies the computation of a large-scale motion with an eective turbulent viscosity. Moreover, one incorporates two variables (and therefore two equations) known as turbulent kinetic energy and turbulent dissipation.

One nds that the RANS-k -ε approach is a practical methodology for numerically obtaining the mean ow and ow components associated with organized structures.

Second, the equations for the two-phase uid ow are presented. For this case, the Euler-Euler equations are proposed. In this approach, the media studied is considered continuous and mathematically interconnected. Each element of the grid has a certain amount of liquid (volume fraction of the liquid) and gas (volume fraction of gas). Furthermore, this approach takes into consideration the four forces involved in liquid-gas interaction: the drag force, added mass force, lift force, and turbulent dispersion force.

The turbulent model proposed contemplates the denition of dynamic viscosity as an addition of a term associated with the shear and other related to the bubble-induced turbulence. One also proposes a R ij -ε model that incorporates the anisotropy of turbulence.

Subsequently, one presents some research works that apply the CFD: the RANS-k -ε and the Euler-Euler approaches for ow determination in stirred tanks and bubbly ows.

This section ends with the presentation of the mathematical principles associated with ROMs: POD and DMD. The ROMs allow decomposing of the uid ow into modes.

Each mode has its own amount of energy (POD) or frequency (DMD) and a temporal part and a spatial part. We also provide a large list of research works in which ROMs have been successfully used for the study of uid ow in several systems, including those that will be analyzed in the following sections of this manuscript. The articles cited show the relevance and eectiveness of ROMs for the analysis of complex uid ows.

The following chapter presents the details corresponding to the numerical implementation of the systems to be studied.

Chapter II Brief details of the performed numerical experiments

The following section is devoted to the exposition of the systems and numerical experiments carried out in this project. The working uid was water at room temperature and atmospheric pressure (density ρ = 998.2 kg.m -3 , dynamic viscosity µ = 1.003 × kg.m -1 .s -1 ) is used as working uid.

II.1 Stirred tank and simulation set up

The stirring speed is 150 RPM (revolutions per minute) what corresponds to an angular velocity ω = 5π rad.s -1 or a frequency f = 2.5 Hz, the Reynolds number is 56250 according to the following relation:

Re = ρf D 2 µ (II.1)
The power number N p for the Rushton turbine is calculated using the torque C on the impeller (P = Cω = 6.5 Watts) according to the expression:

N p = P ρf 3 D 5 (II.2)
From equation II.3 one obtains a value N p of 5.4, using this value we determine the volume averaged viscous dissipation of kinetic energy ε according to the following equation:

ε = N p f 3 D 5 V = 4 27π N p f 3 D 2 (II.3)
where V is the volume of the tank.

The resulting value for ε is 0.090 m 2 .s -3 or 90 W.m -3 .

The numerical mesh composed of 1,129,140 cells and 1,184,282 nodes, the models, settings and the numerical simulation procedure are identical to those used in [START_REF] Delafosse | Analyse et étude numérique des eets de mélange dans un bioréacteur[END_REF]. The standard k -ε turbulence model was implemented with a standard wall function.

A symmetry boundary condition is used on the free surface. The domain is divided into The rst step in the mesh design was a projection of all necessary edges on a horizontal plane followed by the creation of as many surfaces as needed to further build the impeller II.1 Stirred tank and simulation set up shape. Each face was meshed and the corresponding volume mesh was obtained using the sweeping method. Due to this strategy, the mesh is made of prismatic cells and it is almost invariant by rotation around the vertical axis. The simulations and data treatment were performed with ANSYS-FLUENT R20 on a parallel computer provided of with 40 processors Intel Xeon(R) E5-2660, 2.60 GHz.

The time step used in the simulation is ∆t CF D ≈ 5 × 10 -4 s corresponding to an angular rotation of 0.5 degrees per time step and a Courant number less than 1 in the whole domain to ensure numerical stability. Third order discretization scheme was implemented. For each time step, 30 iterations of the non-linear solver were necessary to reach a plateau at 10 -4 for the velocity and continuity residuals.

The convergence toward a stationary ow was evaluated through two criteria: the equality between the torque on the turbine and the torque on the external walls of the tank, and the equality between the volume integral of the turbulent dissipation rate and the power input. To conclude with the presentation of the details regarding the simulation of the stirred tank, it is essential to mention that the numerical results obtained by the CFD methodology were compared with the PIV data reported in [START_REF] Escudié | Structure de l'hydrodynamique générée par une turbine de Rushton[END_REF].

II.2 Bubble column and simulation set up

The geometry used is the same as reported in [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF]. The vessel is of the prismatic type with 0.06m deep (th), 0.35m wide (W ) and 2.0 m height (H ). The liquid height is xed to 1.3m corresponding to a ratio of H W > 3.5 (see Figure II.6). The gas injection sparger is located at 1 2 W , the sparger diameter is 20mm. The gas injection ow rate is 50 l h . Water at atmospheric pressure and room temperature was used as work uid.

The numerical simulation was carried out in the CFD code known as NEPTUNE. The mesh is composed entirely of 2,016,000 hexahedron ("bricks").

Details of the force models implemented in this tool are provided below.

The drag force coecient is proposed by [START_REF] Ishii | Drag coecient and relative velocity in bubbly, droplet or particulate ows[END_REF] and is applicable for small bubbles of the spherical type: where Re b is the bubble Reynolds number:

Re b = ρ l d b |U g -U l | µ l (II.5)
Likewise, the added mass force coecient is proposed by [START_REF] Zuber | On the dispersed two-phase ow in the laminar ow regime[END_REF]:

C M V = 1 2 1 + 2α g 1 + 2α g (II.6)
The lift force is proposed by [START_REF] Tomiyama | Transverse migration of single bubbles in simple shear ows[END_REF]. In this approach the lift coecient C L depends on the bubble Eötvös number E 0 .

Finally the turbulent dispersion force is based on the model developed by [START_REF] Laviéville | A Generalized turbulent dispersion model for bubbly ow numerical simulation in NEPTUNE_CFD[END_REF]. According to this approach the coecients of equation I.17 are dened as follows:

C T D = ⟨F D ⟩ τ t lg -1 b + η r 1 + η r + C M V b 2 + η r 1 + η r (II.7)
with the following denitions:

τ t lg = 3 2 C µ K l ϵ l 1 + β V 2 r K l 1 2
(II.8)

C µ = 0.09, β = 2.7, and V r is the relative velocity of bubbles.

⟨F D ⟩ = 1 8 C D 6 d 2 b ⃗ U g -⃗ U l (II.9) η r = τ lg t τ lg F (II.10) τ lg F = ⟨F D ⟩ -1 ρ g ρ l + C V M (II.11) and b = ρ l + ρ l C V M ρ g + ρ l C V M (II.12)
Finally, it is important to mention that the simulation results for this system were compared with PIV data obtained in the window indicated by the red rectangle in the gure II.7 [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF]

Chapter summary and conclusions

This chapter shows the specic details of the numerical experiments carried out. For the case of the stirred tank, the simulation carried out in the doctoral work of [START_REF] Delafosse | Analyse et étude numérique des eets de mélange dans un bioréacteur[END_REF] was taken as a reference. The corresponding geometry consists of a stirred tank provided with four baes, activated by a six-blade Rushton turbine (with a rotational speed of 150 RPM) and water as working uid at standard atmospheric conditions. The simulation (unstable RANS-k -ε model) whose mesh has a million cells implements the sliding mesh approach to model the turbine rotation. The mesh is therefore divided into two zones: a rotating zone that contains the turbine together with a small volume around it, and a xed zone in which the rest of the volume of the stirred tank is located. After 20 s of simulation, the torque on the turbine is essentially equal to the torque on the external walls of the tank. Similarly, for this simulation time, the volume integral of the turbulent dissipation rate is very close to the input power. These ndings conrms the numerical convergence to a stationary ow. In addition, [START_REF] Delafosse | Analyse et étude numérique des eets de mélange dans un bioréacteur[END_REF] compared the numerical simulation results with PIV experimental data [START_REF] Escudié | Structure de l'hydrodynamique générée par une turbine de Rushton[END_REF]. The results of the simulation show a good agreement with the experimental measurements. For its part, the data of the bubble column were taken from the numerical simulation carried out by [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF]. In this case, the mesh is xed and contains around 2 million grid cells. [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF] compared the numerical results with experimental PIV data for a portion of the study domain. This comparison veries that the numerical simulation leads to results quite close to those measured by [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF].

Therefore, the numerical data to be analyzed in the following chapters are reliable, since it comes from properly implemented numerical simulations whose results were veried with experimental data.

The following chapter describes in detail the implementation of the ROMs for the simulated data of the described stirred tank. Subsequently, the results for high and low frequencies are provided.

Chapter III Baed stirred tank case Additionally, it is necessary to check that the number of snapshots and the time window of data acquisition are sucient to capture the structures carrying the highest variance in the ow. The total simulation time must be long enough to collect sucient information of the rotational motion while the sampling frequency must be high enough to prevent the ltering of small-scale uid motion. Attention must be paid regarding the sampling frequency to avoid the collection of snapshots in phase with the impeller 1 In our context, stationary ow condition refers to a ow that is recurrent and to which a well-dened mean can be assigned with respect to time. More precisely this requirement of stationary ow uid does not necessarily imply a ow that is independent of time. . This gure shows that the sampling frequency is high enough to capture the expected periodic behavior while sampling events are not in phase with the rotation motion. Finally, the frequency of data acquisition must meet the Nyquist criterion. A sample rate of f S of 33.3Hz (0.03s) leads to a Nyquist frequency f N of 16.7Hz. This means that frequencies lower than f N will not exhibit the aliasing phenomenon. As will be seen later, the results obtained satisfactorily comply with this upper limit value of frequency. It is indispensable to emphasize that the present analysis will not capture structures whose period is longer than the total sampling time, i.e., 11 seconds. In this sense, we say that we are performing a high-frequency analysis.

Naturally, to capture low-frequency structures, an analysis time much longer than 11 s is required (see section III.2.2 ).

The sampling location in the simulation domain also implies particularities worth mentioning.

A typical POD decomposition requires that velocity vectors and cell location be expressed in the same frame of reference. This is the case when POD is performed on experimental PIV data, since velocities are generally measured in a xed plane and at xed cell locations [START_REF] Liné | On POD analysis of PIV measurements applied to mixing in a stirred vessel with a shear thinning uid[END_REF][START_REF] Rodriguez | Mixing time and kinetic energy measurements in a shaken cylindrical bioreactor[END_REF]. Some particular attention must be paid when dealing with CFD data. In a previous work, Mikhaylov and co-workers used CFD to compute the velocity eld in an unbaed stirred tank. The problem was solved in a rotating reference frame using a single mesh domain [START_REF] Mikhaylov | Reconstruction of large-scale ow structures in a stirred tank from limited sensor data[END_REF]. Here also, velocities are expressed in the rotating frame and calculated at some xed position in that frame. In the present work dealing with a stirred tank, the presence of baes necessitates to split the uid domain into two mesh zones, a xed one and a rotating one, this strategy is described extensively in [START_REF] Joshi | CFD simulation of stirred tanks: Comparison of turbulence models. Part I: Radial ow impellers[END_REF]. In order to satisfy the above-mentioned condition for a typical POD decomposition, velocities must be expressed in the inertial frame for the xed zone and in the rotating frame for the rotating zone. The export of data, in FLUENT, extracts velocities expressed in the inertial frame. Thus, prior to the POD treatment, the data collected in the rotating zone have to be converted into velocities in the rotating frame (see III.3). The algorithm is applied to this transformed velocity eld, afterwards the reconstructed information is expressed in the inertial frame of reference, and it is added to the one coming from the xed zone. Then, the POD treatment was performed separately for each set of velocities. This per zone approach corresponds to the procedure adopted by Janiga when two mesh zones are present in the CFD model [START_REF] Janiga | Large-eddy simulation and 3D proper orthogonal decomposition of the hydrodynamics in a stirred tank[END_REF]. In that frame, the coordinates of the rotating cells are time independent.

On the other hand, a global POD treatment considering the two zones altogether (with the entire velocity eld expressed in the inertial frame) is attempted with a view of performing the 3D reconstruction of the ow eld in the entire uid domain. This global approach has not been yet proposed in the literature for baed stirred tanks and its feasibility remains to be demonstrated. We will show that the global treatment is feasible bearing in mind that its essential purpose is not to perform a hydrodynamic study of the uid ow structure but a reconstruction of the velocity eld through a linear combination of POD modes. Evidence of that will be provided in the section of accuracy and eciency assessment of the ow eld reconstruction. Proposed conguration for the global treatment: To better understand the features of the global approach on POD results, it is of utmost importance to observe that the data extracted from the CFD code are the components of the velocity vectors expressed in the inertial frame (laboratory reference frame). In addition, this information is referenced with respect to a cell index. Thus, each 3D snapshot consists of a matrix whose rows correspond to the index of the meshing cells or Cell Index (CI), and the columns Subsequently, the snapshot matrix M is constructed by placing in chronological order 3.5 degrees decreases as one approaches the shaft. In consequence, the side view of the velocity eld in the sector may not contain enough information close to the shaft resulting in an apparent absence of data. Although, it can be seen, from the top view (III.8 d) that the ow eld reconstruction at every grid point of the computational mesh is actually obtained close to the shaft.

III.1.2.2 POD analysis in the rotating zone:

After having analyzed the xed zone, the results corresponding to the rotating zone are presented. The POD analysis is performed using the velocity vectors belonging to the rotating zone of the mesh using again 386 snapshots. In this zone, the mesh rotates with respect to the inertial reference frame and the velocity elds ⃗ U IF are obtained in dierent locations for each time instant. Consequently, the velocity elds must be expressed in a reference frame in which the velocity measurement coordinates do not change in time.

We recall here that the CFD code export procedure gives access to the Cartesian velocity components in the inertial reference frame. The velocity components in the moving (or rotating) frame RF , ⃗ U RF , were deduced performing the next steps:

⃗ U ′ IF = ⃗ U IF + ⃗ r(t) x ⃗ ω (III.1) ⃗ U RF = M Rotation ⃗ U ′ IF (III.2)
Where ⃗ r(t) is the position vector of any cell expressed in the inertial reference frame IF . Importantly, the vector position of a cell diers from one snapshot to the other due to the rotation of the mesh. This rotation takes place around the vertical axis (z-axis) aligned with the shaft and the position vector is known from the current cell coordinates

⃗ r(t) = ⃗ r(x IF (t), y IF (t), z IF (t)
). M Rotation is the time dependent rotation matrix that is used to express the velocity ⃗ U ′

IF in a frame attached to the turbine (the position of the mesh in the rst snapshot is used as a reference)

2 . Once the entire data set is treated, the snapshot database now contains a time series of the velocity components in the rotating reference frame RF . Since the mesh is subjected to a solid rotation, the location of the cells, in the moving (rotating) frame, is now also time-independent. Once the transformations of equations III.1 and III.2 have been performed, the POD algorithm can be used.

The POD results obtained are displayed in Figure III.9.

According to Figure III.9 a one single mode dominates the ranking of the variance. This dominant mode exhibits a time-independent behavior as indicated in Figure III.9 b. Additionally, the second and third modes both oscillate with a frequency of 10 Hz or a period of 0.1 s. To validate this time behavior the corresponding data was tted using the Matlab curve tting toolbox® and taking as reference the function A sin (ω t + φ) .

The results from Table III.1 conrmed that these modes oscillate with a period of 0.1 s (10 Hz) corresponding to the passage of the baes observed from the moving frame. This makes sense since from the rotating reference frame, four bae periods must be registered 2 It is convenient to indicate that the rotation operation does not produce a change in the elements of the vector eld.

for each complete turn of the turbine. Reconstruction of the velocity eld from the POD results per zone: After applying the POD method to the xed and rotating zones separately, it is now possible to reconstruct the velocity eld (expressed in the inertial reference frame IF ) for the whole domain of the stirred tank. Since the two zones mentioned above were processed independently, the reconstruction of the velocity eld was done by placing the reconstructed vector elds in each region. Special attention must be paid to two important details.

First, the entire velocity eld was built as a superposition of the velocity elds reconstructed in each zone because the POD components involved are synchronized (the same snapshots were used). Second, all the data (xed and rotating zone) must be expressed in the inertial reference frame IF . In the case of the xed zone, non additional procedure is necessary because the POD technique was applied on a xed grid. Therefore, for this area, the reconstruction of the velocity eld is already expressed in the inertial reference frame IF . The foregoing case is not valid for the rotating zone because the CFD data were processed through the transformations of equations III.1 and III.2. Therefore, the reconstructed velocity eld in this region must be processed to express it again in the inertial frame of reference. Specically, the inverse transformations of equations III. POD treatment is applied to these velocity elds to obtain a reconstructed velocity representation. Then the reconstructed information is again transformed back to the inertial frame. Finally, the POD results of the xed zone are added to those of the rotating zone to generate the reconstruction of the entire domain. The described procedure is feasible but long, complex, and cumbersome.

III.1.3 Global POD treatment

In the above treatments, velocities are expressed either in the xed frame or in the moving or rotating frame. Such velocity elds are also referenced by a cell index, which does not change over time, i.e., ⃗ U IF (CI, t), r(t) = CI. The POD method is a statistical method that seeks a correlation between dierent observed values. Also, the POD procedure does not require that all data be recorded in the same location. In the present case, the cell index is independent of time, while the relationship between the cell index and the location of that cell at any time is easily accessible. Consequently, it is postulated that a global POD treatment considering the entire mesh, with all velocities expressed in the inertial frame, might make sense. The collected snapshots were analyzed using the numerical procedure described in the previous sections, producing a decomposition of the timevarying 3D ow eld on an orthogonal basis of 3D vector elds called modes. It is crucial to note that such a POD decomposition is performed in cell index space, which means that the associated POD modes are not spatial modes and their physical interpretation (as descriptors of the ow structure) has to be questioned 3 . From the vector reconstruction point of view, the above methodology means that the velocity is reconstructed for the whole domain under the geometrical conditions (angle of rotation of the rotating zone) corresponding to a given snapshot. Once the reconstruction has been performed under these conditions, the velocity vectors belonging to the rotating zone must be placed in the appropriate location, taking into account the rotation of the mesh. Presumably, this will not aect the velocity reconstruction derived from the POD method, this hypothesis will be conrmed based on the error values calculated at the of this section.

3 This fact is true for the part of a POD mode that belongs to the rotating zone. In contrast, those parts of a POD mode that correspond to the xed zone meet the condition of the typical POD method and therefore have a straightforward physical interpretation. Interestingly, the third mode is a single-mode, suggesting that it is associated with the reconstruction of the mean ow, the time-independent component of the 3D velocity eld. Similar results were obtained using 250, 286, and 386 snapshots, showing statistical convergence of the POD methodology. The use of three dierent time spans will allow validation of the DMD results in section III.2. After this validation, we will be able to compare the POD results with the DMD results. Subsequent discussion will focus on the rst three modes because they contain almost all of the variance in the system. shows a periodic time variation (again with a period of 0.4 s) but its mean value is dierent from zero. Such a characteristic of the mean value together with the variance content of the mode suggests that this third mode will be associated with mean ow reconstruction.

III.1.3.1 Energy analysis of the Global POD treatment spectrum

III.1.3.2 Modes for the Global POD treatment: Modal components and eigenvectors

The relatively small oscillations presented by that component of the POD are related to the interaction between the periodic ow in the turbine region and the other parts of the Moreover, a tting of the modal components of the rst eleven modes was performed, using as reference the function A sin (ω t + φ) with the Matlab curve tting toolbox®.

The intention is to determine if the higher modes with lower variance content have also a periodical temporal behavior. In addition, obtaining continuous functions of the modal components makes possible a reconstruction in time steps shorter than the sampling interval, which is tremendously useful as it will allow the generation of ow elds that have a higher (or customizable) time resolution. Table III.2 summarizes the results of the corresponding mathematical tting. All the modes considered exhibit a highly organized temporal pattern since all the R square correlation coecients are equal to 1. The rst three modes oscillate with a period of 0.4 s (2.5 Hz) corresponding to the impeller rotation speed while the remaining modes turn out to be associated with harmonic frequencies of the main frequency. It also is clear the presence of paired modes with the same frequency and amplitude. Such a modal conguration implies a strong correlation between duets and is usually an indicator of coherent ow structure being present. Furthermore, Figure III.13 presents the relationship of the rst mode with the fourth through tenth modes.

According to the results of Table III.2, each pair of modes correlates with the rst in terms of Lissajous patterns. For each member of each pair, the number of loops is related to the corresponding harmonic mode. So, for example, the fth mode has three loops since this is the third harmonic of the main frequency associated with the rst mode. The same reasoning applies to the sixth, eighth, and tenth modes. The information obtained above concerning the a k mode components provides evidence that the global POD approach produces an unexpected pattern of the temporal behavior of the corresponding modes. Furthermore, a strong correlation is evident between all frequencies and the rotational speed of the impeller. This result suggests that mode components capture information related to mesh rotation rather than ow structure. Since the amplitude factors associated with each mode are expressed as continuous sinusoidal functions, it is possible to perform a continuous or quasi-continuous reconstruction of the 3D velocity eld at any instant and not only at the instants corresponding to the analyzed snapshots. However, it is crucial to note that the modes produced by global POD treatment are related to cell indices, or in other words, they are not spatial modes as it is the case in typical POD analysis:

⃗ U P OD (CI, t) = N k=1 A k sin (ω k + φ k ) ⃗ ϕ k (CI) (III.3)
For the cells belonging to the rotating zone, the relationship between the cell index and the location in the inertial frame is given by a simple rotation around the z-axis. More precisely, the exact location of any cell at any moment can be easily determined by knowing its initial position and the speed of rotation of the rotating zone. The angle of rotation is θ (t) = ω × (t -t 0 ) where t 0 is the instant corresponding to the rst snapshot.

Finally, for those cells that belong to the xed zone, the relationship between cell index and spatial location does not depend on time this implies that the angle of rotation is constant for all snapshots: θ (t) = θ 0 . To sum up, the global procedure to perform the velocity eld reconstruction contains only three steps:

Perform the POD on the entire domain.

Reconstruct the velocity eld in the entire domain using cell index POD modes.

Relocate the velocity vector of the rotating zone at their actual location.

To better identify the ow characteristics derived from the implementation of the previous steps (and the equation III.3) it is necessary to observe or analyze a time instant of the rotating zone. Figure III.14 depicts the reconstructed ow eld close to the blade impeller region. The moment when the blade passes through the middle of a 3.5 degree circular sector is shown in this representation. Three main characteristics of a ow eld generated by a Rushton turbine are visible. First, there is an axial ow that moves towards the blade. This eect is due to the suction produced by the passage of the turbine blade (red boxes # 1). Secondly, a pair of symmetrical vortices develop behind the blade (red boxes # 2). These highly organized structures represent a manifestation of the well-known trailing vortices, which play a very important role in mixing and transport processes. These organized ow structures would not be present if a time average of the instantaneous velocity eld (snapshots) is performed. This result represents an important advantage of using the POD for velocity eld reconstruction in comparison to a mean ow eld description. Third, an axial jet ow is produced at the periphery of the turbine (red box # 3). The jet ow is, on average, slightly deected upwards (despite a small wiggle) in concordance with [START_REF] Escudié | Experimental analysis of hydrodynamics in a radially agitated tank[END_REF][START_REF] Delafosse | Analyse et étude numérique des eets de mélange dans un bioréacteur[END_REF]. This ow feature is present in the third mode and can therefore be related to the mean ow. After presenting the most relevant results generated by the global POD method, it is convenient to provide a general summary of the properties of this new methodology. The results of the global POD treatment are not as intuitive as those obtained in the typical POD method. For example, the eigenvector obtained through the global treatment does not reveal the structure of the ow. Nevertheless, this non-standard treatment oers some striking features worth noting. First, the resulting modes preserve the variance content and allow a direct reconstruction without applying the per-zone method. This fact is an advantage in terms of computational resource savings by avoiding the use the transformations included in Equations III.1 and III.2. Second, the rotational frequency of the mesh is present in the temporal behavior of all modes. This fact does not reveal a periodic uid motion but reects a correlation between velocity data as the mesh occupies the initial location again after every turn. In that sense, this observation is related to the statistical nature of POD. This type of correlation emerges here in the modal components because the velocity eld in every cell is strongly correlated with that obtained one turn later. Higher frequencies can be related to a high angular periodicity of the such. However, as already emphasized, their contributions to the reconstructed velocity eld can be identied. These rst two modes are involved in the description of the trailing vortices.

On the other hand, the third mode is present throughout the tank, including the rotating zone, and all three components (axial, radial, and angular) are present in this mode. Although, the three components are not present everywhere. In the zone swept by the rotating mesh, the third mode contains essentially an axial component above and below the turbine and an intense radial component away from the blades (see Figure III.15). Considering that the corresponding mode component is almost constant, the third mode seems to contain the axial ow corresponding to the pumping ow rate of the turbine and a radial ow rate associated to the jet ow created by the radial impeller.

Note that both of these ow structures are, on average, independent of the actual mesh elements which is equivalent to a data storage of 9.7 GB. Thus, 121 times more storage capacity is required when using CFD data.

III.1.3.4 Accuracy of the ow eld reconstruction using Global POD treatment

Once the results obtained through the global POD method have been presented, it is pertinent to evaluate the accuracy associated with the proposed POD reconstruction.

Equation III.4 provides the quality of the approximation as a function of the number of modes p used for the ow eld reconstruction from global approach. Equation III.4 estimates a maximum average relative error E M AR . First, the average of the maximum 4 GB = Gigabyte.

error between the reconstructed velocity eld and the CFD velocity eld is calculated for every snapshot. Then by dividing this quantity by the tip blade velocity: U tip = 1.18ms -1 , a normalized averaged error is obtained. This quantity is relevant because the POD technique is designed to minimize the average error between the reconstructed data and the experimental data (CFD data in this case).

Equation III.4 was evaluated using 2, 3, 5 and 10 POD modes (time span of 386 snapshots).

The results are provided in the Table III.3 and depicted in Figure III.16. Using the rst and second POD modes, an average error of 62% is obtained. These modes are independent of the third mode whose contribution allows the reconstruction of the mean ow. Therefore, the POD reconstruction associated with the rst two modes does not provide the mean ow information, such a situation implies a considerably high mean error. The inclusion of the third mode in the reconstruction produces a signicantly lower mean error. This result is fully consistent with the analysis already presented.

E M AR = N i=1 max ⃗ U i P OD -⃗ U i CF D N U tip (III.4)
According to the Table III.3, the maximum error becomes approximately 6.0 % when the third mode is included. Adding more modes to the reconstruction slowly reduces the error down to 3.5 %. This information shows that the global treatment allows the reconstruction of the velocity eld with a low average error.

In addition to the error determined by Equation III.4, it is pertinent to quantify the error considering the number of cells in which the error is signicant, that is, the error as a function of the accumulated percentage of cells. This error computation allows a stricter evaluation of the error committed in the reconstructive approximation of the POD. Equation III.5 estimates the local dierence between the global POD reconstruction of the velocity eld and the CFD velocity eld as a fraction of the local CFD velocity magnitude for each cell of the mesh.

E (⃗ r, t, p) = ⃗ U P OD (⃗ r, t, p) -⃗ U CF D (⃗ r, t) ⃗ U CF D (⃗ r, t) (III.5)
In contrast to Equation III.4, Equation III.5 provides the error due to the reconstruction in a localized and instantaneous manner because the respective calculation is performed for each cell. For the present work, i.e., an unsteady representation of the uid ow, the denition of Equation III.5 is more appropriate as it allows a more accurate evaluation of the POD reconstruction.

Figure III

.17 provides a bar chart of the error distribution corresponding to three dierent snapshots (n = 1, 23, 386) using p = 5 and p = 10 POD modes for the POD reconstruction. For the case of 5 POD modes, about 90.5 % of the total number of cells have an error smaller than 3 %. The relative error of the rest of the cells is distributed according to the following categories: 4.5 % of the cells have an error greater than 3 % and less than 5 %. The remaining 5 % has an error greater than 5 % and less than 27 %.

The previous results are improved by using 10 POD modes to perform ow reconstruction:

98.5 % of the total number of cells have an error less than 3 %, 1 % of the cells have an error between 3 % and 5 %. The remaining 0.5 % error is greater than 5 % and less than 27 %. Table III.4 provides the same results in terms of the percentage of cells with an error within a given range.

As a complement to Table III .4,Figure III.18 shows the information of the localized Number of snapshots E (⃗ r, t, p) < 3% 3% < E (⃗ r, t, p) < 5% 5% < E (⃗ r, t, p) p = 5 p = 10 p = 5 p = 10 p = 5 p = 10 and 10 POD modes. For the case of 5 POD modes, 90.5% of the mesh shows an error less than that indicated by the dotted red line (3%). For 10 POD modes, the results are even better since 98.5% of the cells have an error of less than 3%. The local error provided by Equation III.5 also allows the cells associated with a given error to be located in space. Figure III.19 provides an example for the cases evaluated in Table III .4 and Figures III.17 and III.18. The POD reconstruction is inaccurate in a limited number of cells that are located close to the non-conformal interface connecting the two mesh zones. Increasing the number of modes clearly improves the quality of the velocity reconstruction. It is plausible that the inaccuracy in the at bottom of the tank is related to low-frequency macro instabilities that are not present in the reconstruction [START_REF] Yianneskis | An experimental study of the steady and unsteady ow characteristics of stirred reactors[END_REF][START_REF] Nikiforaki | On the origin, frequency and magnitude of macro-instabilities of the ows in stirred vessels[END_REF][START_REF] Micheletti | Precessional ow macro-instabilities in stirred vessels: study of variations in two locations through conditional phase-averaging and crosscorrelation approaches[END_REF][START_REF] Doulgerakis | On the Manifestation and nature of macroinstabilities in stirred vessels[END_REF]. The last part of this chapter will deal with the corresponding low-frequency

analysis.

An exhaustive study of the precision obtained by using the global POD method as a tool to reconstruct the ow uid in the agitated tank was carried out. The results obtained show that the method is feasible and leads to an accurate reconstructive approximation of the uid ow in the studied system.

III.1.3.5 Numerical cost of a spatio-temporal reconstruction using the Global POD technique

The use of the POD methodology also has to consider the eciency of implementation compared to using standard CFD tools. With this idea in mind, it is interesting to propose a quasi-continuous spatio-temporal reconstruction of the 3D unsteady ow eld, i.e., a POD approximation of the ow ied with a higher temporal resolution than the initial CFD data used in the POD analysis (386 snapshots). Based on the previous POD results and after identifying the continuous function associated with the modal components, we reconstructed 714 snapshots (one snapshot every 5.6 × 10 -4 s) of the 3D ow eld to complete 0.4 s, the duration of one rotation. The temporal parameters, i.e. duration and resolution used to perform this POD reconstruction are the same as in the CFD calculation. Obviously, this analysis will allow comparing the time cost of the velocity eld reconstruction and the corresponding computational uid mechanics approach.

Table III.5 summarises the wall clock times required to simulate a complete turn of the entire domain and perform the POD reconstruction by a global treatment using 5 and 10 modes. The rst two lines refer to the results for ve and ten POD modes. For both cases, the setup time (time to read the data and perform the POD decomposition) is the same because the input data and numerical operations are identical. The reconstruction of the 3D ow eld takes longer when considering ten modes than when considering 5 modes, 133.6 s for 10 modes versus 64.5 s for the 5-mode reconstruction. The corresponding time for the CFD simulation is 2.7 × 10 4 s, which corresponds to 58 (50) times the duration of the POD reconstruction using 5 (10) POD modes.

2.7x10 4

N.A.

The savings in calculation time are signicant, especially when the reconstruction process is carried out only once. Comparing the processing times between the two numerical methodologies shows the advantage of using the global POD method for ow reconstruction. Using the reconstructed information to feed other physical and biochemical models is therefore attractive, especially if the latter seek a more accurate description of the phenomena involved. For example, a quasi-continuous 3D POD reconstruction of the velocity eld with minimal computational cost can be interesting for Lagrangian particle tracking and/or compartmentalization of the uid domain.

The following section presents the details corresponding to the implementation of the DMD method for the stirred tank already studied with POD. The DMD constructs a linear dynamic approximation of the uid ow. The results obtained will be compared with those generated by the global POD method. Some works [START_REF] Semeraro | Analysis of timeresolved PIV measurements of a conned turbulent jet using POD and Koopman modes[END_REF][START_REF] Schmid | Decomposition of timeresolved tomographic PIV[END_REF] suggest that POD and DMD methods should be equivalent for a highly periodic system. There are no convincing numerical examples of such verication in the consulted literature. From our results, we will show that both methodologies are indeed equivalent, which represents an interesting benchmark case for the comparison between POD and DMD in a highly complex experimentally veried system. In addition, it will be shown that DMD can identify and capture low-frequency, low-energy structures reported in previous experimental work [START_REF] Yianneskis | An experimental study of the steady and unsteady ow characteristics of stirred reactors[END_REF][START_REF] Nikiforaki | On the origin, frequency and magnitude of macro-instabilities of the ows in stirred vessels[END_REF][START_REF] Micheletti | Precessional ow macro-instabilities in stirred vessels: study of variations in two locations through conditional phase-averaging and crosscorrelation approaches[END_REF][START_REF] Doulgerakis | On the Manifestation and nature of macroinstabilities in stirred vessels[END_REF]. Obviously, POD was not able to capture these peculiar low energy structures. The possible reasons for this nding will also be discussed in the development of the next section III.2.

III.2 Dynamical Mode Decomposition (DMD) III.2.1 High frequency analysis

The exposition and analysis of the DMD (companion matrix algorithm) begin with the high-frequency case. We will present ndings similar to those already discussed in the previous section III.1. In addition, this section provides information on the dynamic stability and statistical convergence linked to the dynamic decomposition method. Arguments supporting the equivalence between the POD and DMD methods are provided.

III.2.1.1 Sampling and data processing

The sampling procedure used for the implementation of the DMD is identical to the one discussed in the subsection IV.1. The information generated by the method will not contain frequencies (ow structures) beyond the cut-o imposed by the Nyquist sampling theorem. As already stipulated the highest frequency captured by DMD analysis is limited by the expression:

1 2 t S
, where t S is the sampling period. Likewise, the approach used for formulating the snapshot matrix M is the same implemented in the Global POD method.

In this way, the application of the DMD method considers the whole simulation domain expressed in the space of cell indexes, considering the three sets of snapshots or durations indicated in IV.1, however the focus will be on the case of 386 snapshots. It is worth mentioning that similar to the case of POD, the literature consulted does not report the implementation of DMD in domains with moving or sliding meshes.

III.2.1.2 Dynamical stability and statistical convergence of the DMD modes

The DMD algorithm oers the possibility to verify the dynamical stability and the statistical convergence of the calculated modes. This information is relevant and practical as it allows to establish whether the method can generate appropriate dynamics associated with the input data. For example, the position in the complex plane of the eigenvalues λ DM D of the companion matrix M makes it possible to identify dynamically decaying, stable (or steady), and dynamically growing structures according to the criteria described below.

When the eigenvalues lie within a unit circle, the corresponding mode structures have a decaying behavior. The steady modes will present eigenvalues lying on the perimeter of the unit circle. Finally, modes whose eigenvalues lie outside the unit circle correspond to growing mode structures. It is worth mentioning that the modes represented by points close to the unit circle have a dynamic attractor-like behavior [START_REF] Mezi¢ | Comparison of systems with complex behavior[END_REF][START_REF] Mezi¢ | Spectral Properties of Dynamical Systems, Model Reduction and Decompositions[END_REF][START_REF] Rowley | Spectral analysis of nonlinear ows[END_REF]. This fact means that the stable modes will evolve to a steady-state even if their conditions are slightly disturbed. Figure III.20 shows the stability DMD results for the three durations or snapshots sets already analyzed with POD. The three eigenvalue spectrums are mainly composed of dynamically stable modes (dynamic attractor). In addition, we nd some λ DM D values within the unit circle that represents dampening structures. The predominant structures must be in the group of stable modes since the dynamics of the system must be represented in the expansion I.35 by dominant terms that do not diminish or grow without limit. We will soon provide numerical arguments to support such a hypothesis. Figure III.20 also provides the maximum real and imaginary values for the studied λ DM D spectra. The presence of maximum real values equal to unity implies zero imaginary parts and, therefore, the existence of non-oscillatory modes whose magnitude is constant, e.g., the mean ow (third POD mode). In contrast, imaginary parts equal or close to unity imply pure oscillatory modes without damping, i.e., structures with a time behavior identical to the rst and second POD modes. These dynamic features are congruent (but not exclusive) with a closed system conguration like the stirred vessel studied in this work.

Another way to visualize the dynamical stability of the spectra consists of separating and analyzing the exponential and oscillatory parts of each j th mode. The rst part (σ j DM D ) (growth rate) gives visibility to the temporal bounding (exponential behavior) of the modes. The second part (ω j DM D ) is related to the oscillatory behavior of every element of the spectrum.

Next, we will analyze the information regarding the exponential behavior of the spectra.

The results corresponding to the oscillatory part will be presented later with the energy 5 analysis of the DMD modes. DM D = 0 s -1 ). Some of the values lie 5 Note that in this case, the term energy refers to the energy content carried by the vector eld and not to the energy of the uid. This result is based on the denition of the snapshot matrix M which does not involve the volumes of each cell that constitute the numerical mesh below the stability line associated with temporally decaying modes. In addition, certain growth rate values are slightly larger than zero. As will be seen later, such DMD modes carry a very small amount of energy. The largest magnitude for each data set shown is of the order of 10 -3 s -1 . This temporal behavior of the spectra means that most of the terms of the DMD expansion I.35 will not decay or increase in an unbounded manner over time. To conclude the stability analysis, it is relevant to note that the dynamic stability condition is obtained for the three sets of snapshots analyzed. Even the increase of 136 snapshots (when going from 250 to 386 snapshots) or 4.1 minutes (10.2 turns) does not modify the stability results in the DMD spectrum. The most important direct implication of this behavior is that the proposed dynamic decomposition is statistically stable or has reached statistical convergence. These statistical conditions imply that parameters such as average and variance of the velocity at any point in the tank are independent of time.

This scenario is favorable if one wishes to use the proposed decomposition (reconstruction)

as input data to feed models such as the compartmentalization model.

III.2.1.3 Energy analysis of the DMD spectrum

Similar to the POD procedure, the DMD methodology allows calculating the specic energy associated with each mode. These calculations lead to a spectrum ranking. It is worth recalling some of the DMD details mentioned in the chapter I. First of all, we have to point out that the DMD modes are not orthogonal, which implies that the energy rankings of the POD and the DMD do not necessarily coincide (except in very particular cases for highly periodic systems). From this feature derives the fact that the total energy of the DMD spectrum is not in general equal to the total energy contained in the vector eld.

Second, there is no single way to determine the energy associated to each DMD mode.

According to the consulted literature, one can calculate the energy contained in each DMD mode by calculating the square of the amplitude of the corresponding DMD vector Φ j

2

. Alternatively, [START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF] suggests that the square of the amplitude of each mode has to be weighted by the corresponding time coecient to correctly calculate the energy of each element of the DMD spectrum

E j = Φ j 2 e 2σ j DM D τ -1 2σ j DM D τ .
Below we present the results based on the two criteria mentioned above. Surprisingly, the energy rankings for each criterion coincide. Unlike the POD algorithm, the dynamic decomposition procedure allows for an analysis of the energy spectrum in terms of three dierent DMD parameters: σ j DM D or σ j (growth rate), f j DM D or f j (oscillation frequency), and mode number. The following discussion will be carried out in the order of appearance of the parameters just cited. Continuing with the energy analysis of the DMD spectrum, we will carry out a discussion similar to that of the growth rate. This time the reference parameter will be the frequency f j DM D associated with each mode. Figure III.24 gives the energy rankings (typical and alternative criteria) according to the frequency f j DM D in Hz corresponding to each DMD mode. There are some features pre- sented in the gure worth highlighting. First, we note that the results are symmetrical about the 0 HZ frequency. This symmetric conguration implies some practical aspects.

First, the existence of "mirror pairs" in frequency allows the mathematical description of traveling waves which are structures expected in the ow eld description (trailing vortices). Furthermore, the expression I.35 has to produce a real vector eld for k = 0, such a result is achieved when the superposition of the DMD modes annihilates or eliminates the imaginary terms of each element. Such a condition is possible when one has conjugated pairs of modes. In addition, the presence of specic frequencies with a predominant energy content is evident. The top three corresponds to 0 and 2.5 Hz. The next ranking positions are modes with frequencies of 7.5, 10, 12.5, and 15 Hz which are harmonics of the 2.5 Hz.

These predominant DMD frequencies coincide with those of the POD method. Thus the predominant DMD and some POD modes oscillate with identical frequencies. Finally, it remains to verify if the decreasing order of energy of the DMD spectrum matches with the corresponding to the POD. Figure III.25 allows us to conclude that, indeed, both spectra coincide when the corresponding information is sorted in decreasing energy order. As in the case of the POD, the dynamic decomposition generates two modes carrying a total of 81 percent and a third mode with a relative contribution of 16 percent. This result suggests that these top three corresponds to the 1 st and 2 nd , and 3 rd POD modes.

To provide more support to the discussion about the energy content of the DMD spectra, we provide in the Appendix A the gures (analogous to gures III.22, III.23, III.24, III.25) corresponding to the sets of 250 and 271 snapshots. A simple inspection of the attached information allows us to conrm that the results are consistent with those exposed for the case of 386 snapshots. interesting facts. The rst three modes will be sucient to reconstruct the uid ow, as they carry almost all the information in the system, and the third mode is related to the mean ow. We also observe the presence of harmonic modes of 2.5 Hz. These elements of the spectrum lead the energy ranking after the rst three modes. The following section III.2.1.4 presents the results associated with the DMD vectors and the uid vector eld reconstruction. The POD and DMD reconstructions are also compared. The ndings obtained will conrm the equivalence between the two data-driven methods. DMD mode number Normalized energy Growth rate σ j (s -1 ) Frequency f j (Hz)

III.2.1.5 Accuracy of the ow eld reconstruction using DMD

To conclude the high-frequency analysis, we will present the accuracy calculations associated with the DMD procedure. In this case, we use 11 DMD modes to include the conjugate pairs of spectral elements 3,61,279 and 335 (see table III.6). Additionally, for this case, we will take into account the instantaneous elds 1, 23, and 300. The idea of this is to evaluate an additional instant to the 386 already considered in the POD. Similar to the POD case discussed in section III.1.3.4, we will review the accuracy according to the maximum average relative error E M AR and the instantaneous localized error.

For the E M AR estimation, one uses the equation III.4. shows the distribution of the localized error for the POD and DMD in four dierent snapshots. We nd that the percentages of cells obtained for each category of the DMD method are essentially identical to those obtained for the POD case. For 5 DMD modes, about 90.5 % of the total number of cells have an error smaller than 3 %. The relative error of the rest of the cells is distributed according to the following categories: 4.5 % of the cells have an error greater than 3 % and less than 5 %. The remaining 5 % has an error greater than 5 % and less than 27 %.

As with the POD, the results improve when using 11 DMD modes to perform the ow reconstruction: 98.5% of the total cells have an error less than 3%, 1% of the cells have an error between 3% and 5%. The remaining 0.5% have an error greater than 5% and less than 27%.

a) Local and instantaneous relative error for three dierent snapshots using 5 and 10 POD modes.

b) Local and instantaneous relative error for three dierent snapshots using 5 and 11 DMD modes. Finally, we show in Figure III.32 the location of the cells with an instantaneous error bigger than 5% for 5 DMD modes and 3% for 11 DMD modes. We nd that the cells with high error percentages are located essentially in the interface zone between the rotating part and the xed part of the mesh. In addition, one observes some elements in the lower part of the domain (see also gure III.30). This fact probably indicates that the mean ow reconstruction needs to be improved. modes and 424 s (1.6 × 10 -2 CFD) for 11 modes. As we can see the eciency values are approximately similar to those associated with the POD. However, it should be noted that the DMD energy ranking is not obtained automatically as in the case of the POD. An additional eort has to be made in DMD to identify the modes that are dominant and stable.

2.7x10 4

N.A.

After having exposed and analyzed the results derived from the POD and DMD methods we can state some facts. The implementation of these data-driven algorithms for the presented CFD simulation is numerically feasible. Moreover, the time and accuracy obtained to validate the practicality of the proposed uid ow reconstruction. For the simulated system, the equivalence of the implemented methods is evident. The DMD did not generate dominant structures dierent from those found by the POD. This nding allows us to conclude that the POD was able to capture the most important dynamic information in the analyzed data set. Such a scenario is not always valid. This fact means that there may be situations in which the DMD captures dominant structures that are not revealed by the POD. An example is some transient type structures or structures that are not representative of the statistics imposed by the POD on a given database [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. In the following section, we will discuss a case that exemplies this situation.

Specically, we will show that the DMD can capture important low-frequency uid ow components that play a relevant role in mixing processes.

III.2.2 Low frequency -low energy analysis

The application of the DMD in this section will be performed using the CFD velocity eld without taking into account the mean ow. This strategy was implemented in the experimental works in which the presence of low-frequency structures or macro-instabilities (MI's) is reported [START_REF] Hasal | Macro-instabilities of velocity eld in stirred vessel: detection and analysis[END_REF][START_REF] Doulgerakis | On the Manifestation and nature of macroinstabilities in stirred vessels[END_REF]. This choice seeks to eliminate the inuence that can produce the components or modes associated with the mean ow.

The subtraction of the mean ow in the database implies the presence of very small values in many mesh cells which favors the use of the SVD-DMD algorithm. This is because the matrix companion algorithm presents certain numerical accuracy problems when one uses very small velocity values 8 . One observes that the dynamical decomposition allows iden- tifying macro-instabilities with a normalized frequency f ' = 0.02 or f M I = 0.05Hz= f ′ N , being N the impeller rotacional frequency 2.5 Hz. The typical snapshot POD is also used for the study of low-frequency structures. We will show that this procedure cannot capture such long-period structures.

III.2.2.1 Sampling and data processing

In order to evaluate the presence of MI's in dierent regions of the domain space, data were collected in the zones shown in the 8 In the companion matrix method, the calculation of the inverses of the matrix is not always appropriate because the existence of matrices that are close to singular or badly scaled, so the results can be inaccurate.

III.2.2.2 Dynamical stability and statistical convergence of the DMD modes

The dynamic stability results will be performed in the region I because it has a larger number of cells. Moreover, some works report the presence of MI's in this region [START_REF] Nikiforaki | On the origin, frequency and magnitude of macro-instabilities of the ows in stirred vessels[END_REF][START_REF] Micheletti | Precessional ow macro-instabilities in stirred vessels: study of variations in two locations through conditional phase-averaging and crosscorrelation approaches[END_REF].The results for regions II and III will be shown in the energy analysis section.

The database for the region I was divided into three subsets with time spans of 1 minute, 1.5 minutes, and 2.3 minutes. Then the SVD-DMD algorithm was applied for each subset of data. This procedure pretends to evaluate the sensitivity of the results for the three dierent time spans [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF].

III.2.2.3 Energy analysis of the DMD spectrum

The evaluation of the energy content was carried out considering the alternative criteria [START_REF] Tissot | Model reduction using Dynamic Mode Decomposition[END_REF] This should not be surprising since the SVD-DMD does not produce orthogonal vectors so the spectral energy could vary according to the time span considered. The three subsets produce a symmetrical spectrum in which the peaks at frequencies 1.5, 2.5, 15 Hz are noticeable. The presence of 1.5 Hz (1 st and 2 nd in the energy ranking) in the energy spectrum is interesting. Such an oscillation is neither a harmonic nor a sub-harmonic of the fundamental frequency of 2.5 Hz. There is also no mention in the literature of such an oscillation in the system studied. For these reasons, we do not consider it correct to associate any physical interpretation to the corresponding mode. Further analysis will be necessary to establish the possible physical meaning of the 1.5 Hz frequency. There is also a possibility that this oscillation is related to the CFD meshing. The 2.5 Hz frequency (3 rd and 4 th in the energy ranking) is associated with the rotation of the sliding mesh. It makes sense that this fundamental frequency remains important even in the upper part of the domain and on large time scales. The frequency of 15 Hz (one-sixth of 2.5 Hz) is related to the passage of the turbine blades. This high frequency is in the 7 th and 8 th positions of the energy ranking There are also important energy values for the frequencies of 5, 7.5, 10, and 12.5 Hz which correspond to harmonics of the fundamental frequency of 2.5 Hz.

What is surprising is the appearance of a peak with a cluster of elements carrying substantially improve the implementation of order reduction models for the stirred tank.

III.2.2.4 Energy analysis of the spectrum: POD vrs DMD

To conclude this section, we will analyze the results obtained using the snapshot POD for the 2.3 minutes database in the region I. Thus, the sampling process is identical to that applied for SVD-DMD.

The objective of this section III.2.2.3 is to determine whether it is possible to identify some low-frequency-low energy structures from the POD. Similar to the case of section III.2.2.3 we will give priority to those frequencies that have already been found through experiments.

Figure III.39 compares the POD spectrum with the corresponding SVD-DMD. We note that the rst four modes carry the same amount of energy for both methods. From the fth mode onwards the energy content is dierent for POD and SVD-DMD. This result is explained by the fact that the POD modes are orthogonal while the DMD modes are not. Special attention should be given to the 13 th and 14 th modes. This group of elements is the ones that oscillate with a frequency of 0.05 Hz according to the DMD energy ranking. The corresponding POD modes carry lower energy than their DMD analogs. Additionally, these POD elements do not constitute a pair (same energy level) so their frequencies cannot be equal.

Figure III.39: POD and DMD (alternative criteria) spectrum mode for 2.3 minutes. The ranking of the modes is not identical for both methodologies. The energy content for the 13 th and 14 th modes is dierent for each data-driven method.

To verify the time behavior of the 13 th and 14 th modes let's take a look at their modal components a j . Figure III.40 provides the a 13 and a 14 .

It is shown that both modes start oscillating with a low frequency at the beginning of the time span. However, after a certain time (approximately 25 s) the modes are attenuated and start oscillating with a frequency much higher than 0.05 Hz or 20 s period.

Thus even though these modes have the same energy ranking as those of the DMD they do not present an organized temporal behavior with low frequency.

Another possibility to address is to search among the remaining POD modes for one or several that have a frequency of 0.05 Hz. This search represents an arduous task since there are approximately 4700 POD modes. However, for verication purposes, we calculated the fast Fourier transform of the rst twenty POD modes. The corresponding results can be found in Appendix A. It is observed that the rst four modes have a dominant frequency of 1.5 and 2.5 Hz respectively. This is in congruence with the fact that the energy spectrum of both methods is the same for those dominant modes. On the other hand, an inspection of the spectral densities of the other modes (from 5 th up) reveals that no mode is related to a frequency of 0.05 Hz. In fact, it is found that most modes have more than one frequency. An exception is the seventh mode POD which has a predominant frequency of 0.05 Hz. However, when reviewing the temporal behavior of the corresponding modal component, it is not found that its oscillations are stable in amplitude or frequency, similar to the 13 th and 14 th modes (see Figure III.41). This is III.2 Dynamical Mode Decomposition (DMD) natural because the POD does not allow the separation of the data in terms of frequencies. Despite the above, there is the possibility of nding a higher mode with a frequency of 0.05 Hz. However, the corresponding energy content would be practically negligible. The above arguments point to the fact that the POD is not as ecient for capturing low frequencies (low energy) as the DMD, at least for the analyzed database. Obtaining MI's is surprising because one uses CFD data generated with a RANS model. This result could be explained by arguing that the implemented mesh was conceived to perform a simulation for an LES turbulence model. This fact implies that it is possible to generate an excellent description of the mean ow even when using a RANS simulation.

In the next section, we will apply the POD and DMD methods to the study of the last system of our work which is the bubble column. The most interesting thing to be addressed corresponds to the study of the gas phase of the stirred tank.

Chapter summary and conclusions

This chapter deals with the details associated with the implementation of ROMs for the simulated data of the stirred tank. It begins with a description of the sampling and data processing strategies. This stage is of vital importance for the adequate application of the ROMs. This is so since the sampling frequency must avoid the record of snapshots that are in phase with the rotation of the turbine. Additionally, data acquisition must comply with the Nyquist criterion, which prevents analysis errors associated with the aliasing phenomenon. The ROMs were applied following two dierent approaches. In the rst approach (typical method) the order reduction tools are applied to the rotating and xed mesh parts separately. The second method called the global method consists of the application of the ROMs to the entire simulation domain without making a distinction between the rotating and xed parts of the mesh. This approach is innovative and is more direct than the typical method. However, the global method does not allow a direct physical interpretation of the spatial part associated with each mode. This fact does not at all prevent the representation of uid ow in terms of POD or DMD modes.

The results obtained from the ROMs implementation allow identifying and classifying organized structures of high frequency and low frequency. Moreover, it was possible for both approaches the reconstruction of the velocity elds in the analyzed system with a good accuracy and time cost compared to the CFD procedure. It is important to mention that 3 POD/DMD modes are enough to capture 99.9% of the system variance.

The POD and DMD results for high frequency are identical which is associated with the highly periodic features of the ow in the stirred tank. Nonetheless, for the case of low frequency, the POD and DMD results are dierent. The rst ROMs method did not allow the identication of low-energy structures known as macro-instabilities. For its part, the DMD did allow the identication of macro-instabilities (low-energy), which is consistent with the properties of this ROMs method of being able to identify ow components associated with specic frequencies.

Chapter IV Bubble plume case IV.1 Sampling and data processing

The sampling of the simulated data consists of a set of 3D velocity elds (snapshots) collected in a plane of symmetry (y = 0.00125m , -0.18 m ≤ x ≤ 0.18 m, 0 m ≤ z ≤ 1.3 m) of the simulated domain. Once a stationary mean ow condition is established (at t = 50.6 s), a total physical time of 5 minutes (which corresponds to 10 periods of oscillation of the bubble plume), every 0.1 s representing around 2996 snapshots was sampled to be used as input data.

The snapshot POD method and the companion matrix DMD method were implemented to the liquid phase eld α L ⃗ U L (⃗ r, t) and the gas phase eld α G ⃗ U G (⃗ r, t) separately. According to the consulted literature the assessment of two-phase ow has not been deeply studied.

IV.2 Proper Orthogonal Decomposition (POD)

IV.2.1 Kinetic energy analysis

In this case, since the mesh is uniform, the eigenvalues POD represent a direct measure of the kinetic energy of the system. Figure IV.1 provides the normalized modal or eigenvalue spectra for the liquid phase and the gas phase. Interestingly, the two spectra do not resemble each other. Regarding the liquid phase, we observe three dominant modes which contribute around 80% of the total kinetic energy of the system. Subsequently we nd a group of three modes (4 th , 5 th , 6 th modes) very close to each other with energies of the order 10 -2 . Due to their closeness, these might be correlated. Finally, we nd that the higher modes begin to cluster so that it is not possible to identify modes that stand out from their neighbors. a) b)

Figure IV.1: a) POD eigenvalue spectrum for the liquid phase: the rst three modes represents around the 80 % of the total kinetic energy of the system. POD eigenvalue spectrum for the gas phase: the rst three modes represents around the 70 % of the total kinetic energy of the system. b) Cumulated normalized energy POD spectrum for both phases show that a lot of modes have to be taken to capture the 100 % of the total system kinetic energy.

The spectrum of the gas phase shows curious characteristics. There is a continuous decrease of the gas phase eigenvalues in the log-log plot provided. One can propose a expression of the form λ i i λ i = e -n i where n is equal to -5 3 . It is interesting to note that this tendency is mentioned in the experimental work carried out by [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF] Finally, it is observed that for both cases about 400 modes are necessary to capture essentially all the kinetic energy of the system.

We will see shortly that these characteristics of the eigenvalue spectra do not necessarily prevent us from identifying relevant organized structures in the liquid and gas phase uid ows.

IV.2.2 Modal components and eigenvectors

Now we beging to analyse the most energetic elements of the POD spectra for both phases.

Liquid phase

The rst and second modes of the liquid phase show an organized periodic time behavior with a period of approximately 24 s. This period value is reasonable when compared with the period of 28.2 s reported by [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF][START_REF] Laupsien | The organized ow structure of an oscillating bubble plume[END_REF] for the plume oscillation period. Additionally, these elements of the spectrum are approximately out of phase by π 2 as can be seen in gure IV.2 b. It is expected that these are related to the generation of vortices that move in space. For its part, the third mode presents small uctuations but is constant. Due to its time behavior and its energy content, this mode must be related to the mean ow. This is in agreement with the POD analysis performed by [START_REF] Laupsien | Hydrodynamics, Mass Transfer and Mixing induced by Bubble Plumes in Viscous Fluids[END_REF] in a smaller region of the ow domain.

Regarding the eigenvectors, we observe a clear spatial organization for the rst 6 modes shown in gure IV.3 and IV.4 . The 1 st mode eigenvector presents two large vortices whose dimensions are of the order of the vessel height. The 2 nd mode eigenvector presents three vortices, this time smaller in size but still comparable to the dimensions of the tank. One also observes that these modes (1 and 2) have similar wavelengths.

These rst two modes have to be related to the plume oscillation and will clearly play a determining role in the mixing process. The 3 rd mode eigenvector reveals two large circulation loops typical of the average ow in this type of device. This spatial conguration conrms that this mode represents the mean ow. Finally, the fourth, fth, and sixth modes also present a certain degree of spatial organization but their structures are smaller than the lower modes (1 and 2). Their wavelengths are half the wavelengths of the modes 1 and 2. Due to this conguration and its temporal behavior, these modes are related to the periodic movement induced by the gas on the liquid phase. Once the eigenvectors and modal components for each phase have been determined, we can continue with the reconstruction of the velocity eld for both the liquid phase and the gas phase. To visually evaluate the results of the POD reconstruction, we present a series of reconstructed instantaneous elds together with the velocity elds corresponding to the simulation. The gures IV.9 and IV.10 provide the vector plots for the instants for both phases. For the instants shown in gure IV.9 we nd that the conguration of vortices for both POD reconstruction (600 modes) and CFD scenarios coincide quite well.

However, it is clearly observed that the reconstructed elds do not faithfully replicate the detailed description of structures present in the CFD. This suggests that the ROM's tool was not able to eciently extract such small structures. However, it is notable that large organized structures are captured by the POD.

For the gas phase, we also observe that the POD and CFD reconstruction agree relatively well. Once again it is observed that the reconstructed elds do not faithfully replicate the detailed description of structures present in the CFD. In addition, a widening of the gas jet is observed. Despite such divergences, it is notable that POD is capable of replicating the overall structure conguration generated by CFD. IV.2.4 Accuracy of the ow eld reconstruction using POD treatment Due to plume oscillations, the velocity value may be zero or very small in some cells for the liquid phase and especially for the gas phase. Such characteristics make it impractical to determine the error from equation III.5 because of the possibility of division by zero.

Therefore we propose the standard calculation of the mean POD error (DMD) according to the following equation IV.1:

E M AR = N i=1   ⃗ U i Reconstruction -⃗ U i CF D N   1 max ⃗ U CF D (IV.1)
where: i corresponds to the i th mesh cell, ⃗ U i Reconstruction -⃗ U i CF D is computed in every mesh cell, N is the total number of the snapshots and max ⃗ U CF D is the biggest velocity magnitude in the system.

To evaluate equation IV.1, the reconstruction of the velocity eld was carried out for each phase, taking 2996 POD modes. Then the corresponding error calculation was performed. As shown in table IV.1, the maximum mean error is similar for both phases.

It is clear that increasing from 600 modes to 1500 modes does not signicantly reduce the error. This is due to the fact that the higher modes in the spectrum are very close to each other. It is noted that the case of the reconstruction with 600 modes gives an error of 3% for the liquid phase and approximately 1% for the gas phase.

IV.3 Dynamical Mode Decomposition (DMD)

As we saw in the POD section, the time behavior for higher modes for the two phases analyzed does not show a high degree of organization. The small uctuations of the mean ow probably prevents capturing organized time behavior for higher modes. This eect could be even more intense for the gas phase for which the mean ow is the highest energy mode. For this reason, we proceed to perform the DMD analysis by subtracting the mean ow from the data set. Once again this strategy will allow us to capture the structures of low frequency and low energy content. In this case we applied the companion matrix DMD method.

IV.3.1 Dynamical stability and statistical convergence of the

DMD modes

As usual, we start with the dynamic stability analysis. As can be seen, essentially all DMD modes for both phases are extremely stable. This result is a good sign since it indicates that the ROM's tool has the possibility of generating a stable linear model of the system.

In gure IV.12 is also shown that the growth rate values for the spectra of both phases are very close to zero. This characteristic of the spectra implies that there will be no burst-like structures. Additionally, makes us choose the typical criterion to determine the energy associated with each mode (growth rates close to zero can produce a numerical divergence for the Tissot criterion, see equation I.53 ). Such results of the energy ranking will be seen below.

IV.3.2 Kinetic energy analysis

In gure IV.13 we see that the two spectra are symmetric in terms of energy ranking.

It is evident that the spectrum is dominated by low frequency structures in the vicinity of 0.043 Hz, which corresponds to a period of approximately 23 s. This period value is very close to that found for the POD case. In this way the results of both tools are coincident. However, in the case of DMD, it is possible to identify several structures with unique periods of low frequency, while with POD we nd a period that corresponds to a superposition of dierent frequencies. In the gas spectrum, modes 59 and 60 is the most energetic and has a period of 0.04339 Hz. It is therefore interesting to observe its spatial structure.

IV.3.3 DMD vectors

The gure IV.14 depicts the real part and imaginary parts of DMD vectors 507 (rst of energy ranking). An organized conguration of vertical structures is clearly observed. The dimensions of these structures are of the order of the size of the vessel. It is interesting to note that the real part of the DMD 507 vector (0.04 Hz) presents a structure very similar to that shown by the second POD eigenvector. Moreover, the corresponding imaginary part shows a spatial conguration that resembles to that of the rst POD eigenvector.

It is important to mention that these spatial characteristics are present in vectors from second to eighth DMD modes in the energy ranking. Finally, gure IV.16 provides the spatial conguration (real and imaginary parts) of the rst two DMD modes according to the corresponding energy ranking (gas phase). One notes that there are some similarities of these vectors to the second, third, and fourth POD the accuracy of the implemented ROM's methods. However, one of the objectives of the order reduction is to determine the possibility of eciently capturing organized ow structures. According to the results obtained for the liquid phase (POD, DMD), less than ten modes are necessary to reconstruct a ow that presents structures associated with large-scale organized motions of the ow uid (see Figures IV.17 and IV.18). In fact, as discussed in section IV.2.2, the second, fourth, fth, and sixth modes are harmonics of the rst POD mode (according to the ratio of the corresponding wavelengths). This implies that for the liquid phase the ROM's used can extract coherent structures associated with dierent scales even taking a relatively small number of modes. In contrast, for the gas phase, it is found that even taking the rst six modes it is not possible to obtain structures that show noticeable changes. This means that in the case of the gas phase it will be necessary to take a larger number of modes compared to the liquid phase to detect structures with dierent wavelengths.

Chapter summary and conclusions

This section describes sampling and data processing. For this case, the data was taken in a plane of symmetry and for a span equivalent to 10 oscillations of the bubble plume.

The implementation of ROMs for the liquid and gas phases was carried out separately, which is an innovative proposal since such a scenario has not been addressed before.

The results obtained show that the ROMs were able to eciently identify large-scale organized structures such as vortices associated with the oscillations of the bubble plume (rst and second POD modes) and the mean ow of the bubble plume (third POD mode) present in the phase liquid, these 3 modes contain 80% of the total kinetic energy of the corresponding phase. The time parts of the rst and second modes have a period of about 24 s whose value is close to the oscillation period of the plume found experimentally.

For the gas phase, ROMs also allowed the identication of large-scale organized structures.

In this case, the rst POD mode represents the mean ow, the second and third modes present an organized conguration with wavelengths similar to the rst and second modes of the liquid phase, respectively. These rst three modes contain about 70% of the total This information can usefully be used to perform the simulation of a (bio)chemical reactor coupling transport and reaction. In the following, the mixing of an inert scalar, the transport of massless particles as well as homogeneous and heterogeneous catalysis numerical experiments will be presented.

V.1.1 Scalar transport

Once built, the matrix allows writing a simple expression for the transport of a scalar in the space of compartments. Indeed, the conservation of c in the compartment i writes

dV (i)c(i) dt = k c(k)F (k, i) -c(i) j F (i, j) (V.1)
The rst term in equation V.1 corresponds to the summation of mass uxes of c originating from the compartments k and entering the compartment i. The second term on the right-hand side is the sum of mass uxes leaving the compartment i. Thus, forming a new matrix M such that the diagonal is lled with the sum of owrates exiting a compartment, is helpful to write the mass balance in a more compact form

     M (i, j) = F (i, j) for i ̸ = j M (i, i) = -j F (i, j) (V.2)
Let C be a N s × N c matrix, where N s is the number of scalars transported by the uid ow. Thus C(m, i) contains the concentration of the scalar m in the compartment i. The mass balance for each scalar over the entire set of compartments can be written in a compact form as

dCV dt = C.M + SV (V.3)
S is a N s × N c matrix of volumetric source terms due to interfacial mass transfer and reaction. V is a N c × N c diagonal matrix made of the volume of each compartment.

It is worth observing that there is no spatial information associated to the compartments. The information regarding the indices of the compartment adjacent to a given compartments is now embedded in the matrix M. Typically, the neighbors of any com- partment, i, are found as the column indices of the non-zero terms on the corresponding line. The knowledge of the compartment numbering procedure allows the creation of visual representation that resembles the physical uid domain.

V.1.2 Particle transport

We shall now consider massless particles and their movement in the space of compartments, one takes as reference the approach reported in [START_REF] Delafosse | EulerLagrange approach to model heterogeneities in stirred tank bioreactors Comparison to experimental ow characterization and particle tracking[END_REF]. Any particle can be indexed by the compartment it belongs to at time t. Let, idxloc(p) be the index of the compartment in which particle p resides. The displacement of particles actually consists in updating that table of particle location index idxloc.

The probability that a particle present in the compartment i, will leave the compartment during the time interval ∆t is

     P out = 1 -e -∆t/τ (i) P out ≈ ∆t/τ (i) if ∆t << τ (i) (V.4)
τ (i) is the mean residence time in the compartment i dened as

τ (i) = V (i) j F (i, j) (V.5)
Assuming that the particle under consideration is to leave the compartment i, the probability to reach the compartment j writes

P j i = F (i, j) F (i, j) (V.6)
A given compartment has a maximum of six neighbors, so for any compartment i one will actually compute a maximum of six probabilities. These are stored in a N c × 6 table, P (i, k). In that table, some elements may be zero if the number of neighbors is actually less than 6. We also dene another N c × 6 connectivity table, idxngb, which contains the compartement indices of the said six neighbors.

In order to select the destination compartment from a random number trail, it appears interesting to build a new N c × 6 matrix CP of the cumulative probability as

CP(i, m) = m k=1 P (i, k) (V.7)
By construction, CP is a stepwise function that exhibits at most 6 dierent values in [0, 1].

It is shown in gure V.2 how the destination compartment is identied using a random number sampled in a uniform distribution. The y-axis presents the cumulated probabilities to exit a compartement. Two examples are illustrated : the blue one corresponds to a compartment with ve neighbors, the orange one corresponds to a compartment with 4 neighbors. Accordingly, the [0, 1] interval on the y-axis is split into 5 and 4 subintervals, marked by blue and orange dots respectively.

The procedure start with the trail of a random number β ∈ [0, 1]. The destination compartment is found through the identication of the interval β belongs to.

CP (i, k) < β < CP (i, k + 1) (V.8)
The procedure rst identies an integer between 1 and 6 and then the connectivity table idxngb provides the corresponding compartment index. From there the table of particle location index is updated. The following algorithm describes the step by step procedure used to update the particle location index using two random number trials for each particle.

Algorithm 1 Calculate the new location index of a particle for p = 1 to N p do trial a random number α in the uniform distribution U(0, 1) if α < 1 -e -∆t/τ (i) then the particle leaves the compartment trial a random number β in the uniform distribution U(0, 1) idxloc(p) ← idxngb(i, 1) (assuming β < CP (1))

for k = 1 to 5 do if CP (i, k) < β < CP (i, k + 1) then idxloc(p) ← idxngb(i, k + 1)
end if end for end if end for

V.1.3 Data processing and preliminary calculations

The CFD data correspond to those computed in the previous chapter, i.e. in a 70 L reactor lled with water and equiped with a Rushton turbine rotating at 150 rpm. The CFD data set stored over one rotation of the impeller are processed so as to produce 14 matrices F corresponding to the instantaneous owrates for 14 dierent positions of the impeller.

Hence, an identical number of matrices M, P and CSP are build and stored. In contrast, the connectivity table idxngb is unique since the indices of neighbors are independent of the actual owrates between compartments. Considering instantaneous owrates and consequently the uctuations with respect to the time average value is interesting here because it naturally produces some kind of dispersion in the particle pathways. In contrast to the strategy adopted by [START_REF] Delafosse | EulerLagrange approach to model heterogeneities in stirred tank bioreactors Comparison to experimental ow characterization and particle tracking[END_REF], we do not refer to a turbulent transition matrix to introduce some stochasticity in the particle displacement algorithm. Indeed, the turbulent motions act on the particle trajectory over a length scale that is much smaller than the compartment characteristic size. Turbulent velocity uctuations most probably cannot be considered as the source of particle dispersion at the compartment scale. Our approach do not require a model for particle dispersion (whatever the physical justication supporting this model). Indeed, looping over a set of instantaneous owrates allows simulating the blade passages in a more realistic way and produces some changes in the transition matrix (probabilities P k ) from one instant to another.

Transition Matrix: As far as the particle motion is concerned, the prerequisite is that the statistics of particle displacement match the expected values computed from the known owrates between compartments. In order to check this point, 10 5 particles were randomly distributed in the compartments and mixed during 30 seconds while the transition matrices P were periodically updated every ∆t u = 1/(2.5 * 14) ≈ 0.0286 seconds.

During that process, all exit events were recorded and cumulated in a N c × 6 table. This allows to compute the eective transition probabilities n(i, k)/ n(i, k) and compare them with the expected values. Alltough the transition matrices P i (k) change in time, the accumulation of all exiting events in a single table, without refering to the current transition matrix, provides an information about the transport of particle by the mean ow. Hence, the expected values were computed from equation V.9, using the temporal mean of the 14 matrices F. About 9 Million exit events were collected.

P i (k)expected = F (i, j) F (i, j) (V.
V.2 Results

V.2.1 Mixing of a inert scalar

We rst study the mixing of an inert scalar. In this context, one usefully denes

The mean volumetric concentration of a scalar c as

⟨c(t)⟩ = 1 V Nc i=1 c(i, t)V (i) (V.10)
The variance of the concentration distribution

σ c (t) = 1 N c Nc i=1 (c(i, t) -⟨c⟩) 2 (V.11)
The segregation index, or the normalized variance Mixing is often described as a process that reduces the variances of a distribution.

I c = σ c (t) σ c (0) (V.12)

V.2.2 Mixing of inert particles

In this part, the mixing of massless particles is simulated. In this context, we dene

The particle density (in number), i.e the number of particles, n p (i, t), in the compartment of volume V (i)

ϕ p (i, t) = n p (i, t)/V (i) (V.13)
The mean particle density

⟨ϕ p (t)⟩ = 1 V Nc i=1 n p (i, t) = 1 V Nc i=1 ϕ p (i, t)V (i) (V.14)
the variance of the particle density and the denition of the segregation index are unchanged. This because, the particle density in number is analogous to a concentration : a number of molecules/particle per unit volume.

All particles are initially placed in a given compartment and then displaced following the algorithm previously described. The same experiment is repeated for dierent numbers of particles N p ∈ {0.5, 1, 5, 10} × 10 5 . The number of particles per compartment at each time step is computed and stored.

The local particle density in ve randomly selected compartment are presented on gure V.5. The broken line aspect is due to the fact that results are collected every second.

V.2 Results Particle homogenization takes place between 15 s and 20 s whereas the scalar mixing time computed by Delafosse using CFD (LES) was found to be as 14 s. This slight overestimation is a consequence of the severe reduction in ow description (500 compartments against 1,2 M cells). It therefore appears that the proposed strategy provides reasonably accurate results for a minor computational cost, on condition that the unsteadiness of the ow is considered. However, it appears that the particle density uctuations are all the more pronounced than the number of particle is small. This point will be discussed in the next paragraph.

V.2.3 Cross-validation of scalar and particle mixing

We observe in gure V.6 that the reduction of variance is very similar for scalar mixing and particle mixing. however, as explained before, the residual variance of the particle density is dependent on the total number of particles. In order to provide a more detailled view on the dynamics of mixing, the normalized concentration and normalized particle density are plotted against time for ve randomly selected compartments, in Figure V.7. A very good agreement is found between the local densities computed with an Eulerian approach and a stochastic method for particles. This indicates that the particle displacement algorithm is not only good to describe the overall dynamics of mixing, it is also capable, if the number of particles is large enough to describe V.2 Results the local evolution in particle density throughout the mixing process. V.2.4 Understanding noise on the number particle density Despite the mixing process has reached the point where the mean number of particles per unit volume is constant, the number of particle per compartiment still uctuates from one iteration to another due to the stoachastic nature of the particle displacement. Moreover, the magnitude of this uctuation seems to be related to the compartment itself and the total number of particles, as illustrated in Figure V.5.

Let us dene the discrete random variable X i,n (t) such that X i,n (t) = 1 if a particle is in the volume V (i) and X i,n (t) = 0 if a particle is not in the volume V (i). This variable follows a law of Bernoulli with parameter u i such that      P(X i,n (t) = 0) = 1 -u i P(X i,n (t) = 1) = u i (V.15)

A property of the law of Bernoulli is that the expectation of the associated variable is equal to the parameter u i and the variance is equal to u i (1 -u i ).

The expectation of X i,n (t) is the number of events true divided by the total number of events. In our particular case, the total number of events is limited to the total number of particles tracked. The event true reads the particle is in the compartment i .

V.2 Results since V (i) << V in the present case, the expression may be simplied into

∆n p ⟨n p (i)⟩ = N p V (i) V -1/2 (V.24)
As the total number of particles goes to innity, the relative change in particle number goes to zero. However, this convergence is slower for smaller compartments. In other words, once the stationary condition is reached, the relative change in the local number of particles is larger for small compartments and it is inversely proportional to the square root of the total number of particles.

In order to check the consistency of the algorithm implementation, the number of particles in each compartment is stored every seconds form t = 20 to t = 30 seconds, that is to say once mixing is achieved. The standard deviation in each compartment is quantied and normalized by the expected number of particles, readily accessible from equation (V.18). Results are presented in Figure V.8. It is clear that the standard deviation decreases as the number of particle increases : groups of black circles are much closer to the origin than green circles are. One can also observe several horizontal lines of circles among the same color which correspond to compartments with similar volumes. Each of these clouds is more or less centered on the parity curve. Averaging over compartments of the same size would cleary improve the results. The relatively large spreading is actually due to the fact that we only used 10 particles counts to perform the statistics. Anyway, the trend issuing from the particle mixing experiment is in very good agreement with the analytical prediction given by equation (V.24).

V.2.5 Coupling with the reaction

In this part, we perform the comparison between to options to calculate the source term in a scalar transport equation. For illustration purposes, a rst order kinetic reaction is The scalar concentration in the reactor is initially homogeneous as well as the particle density in number. In order to perform the comparison, two scalars are transported by the ow. However, the source term for the rst one obeys the homogeneous catalysis formulation whilst the source term for the second scalar is computed through the heterogeneous approach (as the sum of particle contributions). Obviously, in this example, both approaches should give the same results, i.e. C(t) = C 0 (1 -e -kt ). Indeed, gure V.9 shows the evolution of the concentration of a scalar in ve compartments (randomly selected among 500) using either the homogenous or the heterogeneous catalysis approach. The continuous line is actually a superposition of ve concentration proles using the homogeneous catalysis approach. The ve proles Dashed lines stand for the normalized particle density which remain constant over time.

the calculation of the local source term is not strictly identical in every compartment due to the noise in the particle number. As a results, the open circle do not exactly overlap at each instant. However, in the present case, the total number of particles (around 10 5 ) is large enough to ensure that the deviation in terms of particle number is negligible.

Hence the local source terms computed using the particle based approach are practically identical.

V.3 Extension to biological populations

From a fundamental point of view, the combination of a compartiment model approach and the representation of the biological population as an ensemble of particles permits to overcome the usual limitations encoutered when the number of properties required to describe the state of the bioparticle is large. Indeed, the set of equation describing the changes in dissolved species concentrations contains an integral term refering to the number density function. One of the possible strategy in that case is to solve a population balance equation.

∂n(x, ξ, t) ∂t + ∂ ẋn(x, ξ, t) ∂x + ∂ ξn(x, ξ, t) ∂ξ -Γ(ξ)n(x, ξ, t) . . .

• • • = Γ(ξ ′ )K(ξ, x ′ i )n(x, ξ ′ , t)dξ ′ (V.27)
ẋ it the velocity of particles, ξ describes the velocity in the space of particle properties, Γ(ξ) is the division frequency, K(ξ, ξ ′ ) the redistribution kernel dening the probability to form a particle with properties ξ from the division of a particle with properties ξ ′ .

The equation for the dissolved scalars consummed by the biological particles is generally expressed as ∂C(x, t) ∂t + ∂ ẋC(x, t) ∂x = q(C, ξ)m(x, ξ, t)dξ (V.28) q(C, ξ) is the vector of specic consumption rates per unit mass of bio particle, m(x, ξ, t)

is the mass of bioparticles in the state ξ. The latter can be accessed from the size and density of particles which can actually be part of the vector of properties.

The last equation refers to the dynamics of cell properties: ξ = f (q, ξ)

(V.29)
We showed in the previous section that we could compute the source term in the scalar equation as a sum over the particles present in a compartment at a given time.

q(C, ξ)m(x, ξ, t)dξ ≈ np(i,t) p=1 q(C, ξ p )m p x i V (i)

np(i,t) p=1 m p (V.30)
In other words, the knowledge of the location of particles is sucient to solve the conservation equation for the scalar. Moreover, the evolution of particle properties must be updated using equation V.29 which poses no real diculty. Finally, the remaining task lies in the denition of the division kernel as well as the algorithm to set the properties of the daughter cell. Here again, a constant number Monte-Carlo strategy can be used.

Chapter summary and conclusions

This section proposes a CMA approach that incorporates time into the modeling process. This approach consists in switching (computing) periodically the CMA uxes matrix. In this way, the ows between compartments will change over time due to the rotation of the turbine. In this framework, the time uctuations of the ux values naturally produce some kind of dispersion in the particle pathways. The methodology proposes to model scalar and particle transport phenomena. For the case of scalar transport, a mass conservation equation is proposed in which the incoming and outgoing ows of each compartment depend on time. For the case of particle transport, a stochastic approach is used. In this methodology each particle is referenced by the compartment it belongs to, at time t. Two stages are considered that allow modeling for each snapshot the transport of the particles in the studied tank. In the rst stage, it is evaluated if a particle is going to leave its residence compartment. Subsequently, if a particle leaves its compartment of residence, one determines to which neighboring compartment it will migrate. The phenomena of scalar and particle transport (mixing), and coupling with a chemical reaction (homogeneous and heterogeneous catalysis) were modeled. The corresponding CMA was based on CFD data from the stirred tank simulation described in section II.1 for ve hundred compartments and one full turn of the Rushton turbine. For the case of scalar transport and particle transport, the results show that the temporal evolution of the concentration per compartment is consistent, since a state of homogenization is reached in a time very close to that calculated by means of a LES simulation [START_REF] Delafosse | Analyse et étude numérique des eets de mélange dans un bioréacteur[END_REF]. Finally, the results for the homogeneous and heterogeneous catalysis reactions coincide quite well, which shows the correct modeling of the transport phenomena involved. These ndings show that the proposed dynamic CMA approach is promising and opens one way for the generation of models that integrate phenomena associated with the biological population in stirred tanks.

Conclusions and perspectives

POD and DMD were applied to U-RANS CFD data describing the turbulent ow in a baed stirred tank.

Energy spectrum and ow structure analysis can be performed in a classical manner in each zone of the mesh (xed or rotating). As expected the main ow structure and the trailing vortices could be identied from the data analysis.

An original methodology is proposed to perform the spatio-temporal reconstruction of the velocity eld in the entire volume. In principle, this objective is out of reach due to the fact that the cell location is time dependent. We showed that, although the physical meaning of modes is lost, the ROM's methods (POD and DMD) can serve the objective of reconstructing the ow eld.

Both POD and DMD are adapted to this reconstruction as proved by an error analysis between the initial (CFD elds) and reconstructed ow elds. As expected, the accuracy increases with the number of modes. In the very case of an unsteady turbulent ow simulated using the U-RANS approach, the rst 15 modes are sucient to recover the ow eld with an accuracy lower than 3%.

The 3D unsteady ow eld is reconstructed at a very moderate cost compared to that of the CFD solver. Obviously, the generation of CFD data remains time consumming but the investment is highly interesting if the ow is periodic.

The POD results are identical to those obtained by the DMD for the high-frequency analysis. Whereas POD cannot identify low-energy structures, DMD is able to identify low energy (low frequency) structures (macroinstabilities) because such structures are not associated with high energy contents. Tracking such macro-instabiities in a set of numerical data coming from U-RANS CFD simulation represented a considerable investment and we showed that the correct strategy is to sample over a long period of time and treat the data with the DMD method. Also, the observation of macro-instabilities in our U-RANS simulation is clearly related to the renement of the mesh and the small time step.

The ROM's methods (POD and DMD) were also applied to CFD data describing the turbulent ow for a bubble plume in a quasi-2D column. In this case, a xed uniform mesh was used in the numerical simulation so that the POD and DMD energy spectra can be directly associated with the kinetic energy of the system. Both POD and DMD were applied for the velocity elds of the liquid and gas phases weighted by their respective volume fractions. The POD implementation for the liquid phase reveals that the rst 3 modes are sucient to capture 80% of the kinetic energy of the system. The rst two POD modes are associated with the oscillation of the plume and the third mode represents the mean ow. The modal components a k of the rst and second POD modes show a periodic behavior with an oscillation period of 24 s which is close to the 28.2 s oscillation period of the plume found experimentally. Likewise, the respective eigenvectors present a spatial organization whose wavelengths are comparable to the dimensions of the tank. It also is observed that the higher POD modes are the harmonics of the rst mode. Regarding the gas phase, the POD energy spectrum indicates that the rst three modes represent approximately 70% of the corresponding kinetic energy. Interestingly the slope of this spectrum in the log-log plot showed a slope of -5 3 such a result may give clues about important uid ow phenomena. In addition, the rst six modes presented a slightly organized behavior in time, suggesting an oscillation period of approximately 20 s. Concerning the spatial part of the POD modes, it was found that the rst mode is associated with the mean ow of the gas, this is the gas jet. The second mode shows an organized structure with a wavelength similar to that of the rst mode of the liquid phase. The same can be said if the third eigenvector of the gas phase is compared to the second eigenvector of the liquid phase.

Similar results were found for DMD in both phases when the data was analysed without the mean ow. For this case, very stable structures are observed with a maximum at the frequency of 0.043 Hz or approximately 23 s. This nding is in agreement with the POD results.

This work opens interesting questions related to the POD and DMD analysis of the gas phase.

Finally, a strategy to perform the transient simulations of reacting unsteady ows for mixing tanks in real-time is proposed. It consists simply in switching periodically the matrix of ow rates in a Compartment Model Approach.

As far as the application to stirred tank is targetted, several point would deserve further attention. The full integration of the velocity elds reconstruction results with mixing models such as the CMA is envisioned for future work. To achieve this goal it should be foreseen that the reconstructed velocity eld data are rewritten according to the formats used by the CMA tools available in the research group. Indeed, for the moment, matrices of owrates are build directly from CFD data and not from the ROM's reconstruction of velocity elds. The benet of using a reconstructed ow eld is that it contains the contribution of macro-instabilities, especially if the reconstruction is performed with DMD.

Once the conditions are in place to generate reduced order models: CMA + POD or CMA + DMD, the impact on mixing performance due to the use of reconstructed velocity elds can be veried.

The main diculty encountered in this work, among several others, was related to the sliding mesh. The fact that the input data are generated at nodes whose location is time dependent creates serious diculties and also forbid the standard interpretation of modes as velocity vector elds. One option to circumvent this limitation is to interpolate the instantaneous velocity elds on a xed grid when sampling it. Whether this can be done within the CFD code or after exporting the data has to be examined.

Saving the data on a xed grid would also greatly facilitate the construction of the compartments and the calculation of the matrix of owrates. Quite obviously, this improvement would allow to analyse time resolved input data issued from Large Eddy Simulation. This is actually ongoing work but the number of snapshot to consider is enormous due to the much wider range of scales that are produced in such highly resolved simulation.

The ow structure analysis of LES input data is certainly interesting but the use of

  . . . . . . . . . 49 III.1 Sampling data sets collected for apply POD. Each rectangle represents a snapshot of registered velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 III.2 Schematic representation of the data sampling frequency. The continuous plot is given by y = sin(5πt). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 III.3 2D schematic representation for the velocity vector expressed in the inertial and rotating frames. The inertial frame is equivalent to the laboratory frame of reference. The inertial and rotating frame share the same origin. The rotating frame rotates around the z -axis with angular velocity ω. In that frame, the coordinates of the rotating cells are time independent. . . . . . . . . . . . . . . 54 III.4 Schematic representation of the typical and global treatments. . . . . . . . . . . 55 III.5 Schematic representation for the n th snapshot of the 3D velocity eld ⃗ U -column array representation for the n th snapshot of the 3D velocity eld ⃗ U (n) IF (CI). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 III.7 Schematic representation for the snapshot matrix M . The superscript notation (N ) refers to the N (th) snapshot. . . . . . . . . . . . . . . . . . . . . . . . . . 57 III.8 POD analysis of the ow in the xed zone of the mesh: a) mode spectrum, b) mode component associated with the top three ranking modes, circles: rst mode, x's second mode and rhomboids third mode (each data set is normalized by the variance corresponding to each mode.), c) First eigen-vector visualized in a vertical XZ plane, d) Top view of the rst eigen-vector. The main feature of the mean ow, i.e. radial jet ow, re-circulation loops and coherent vortices behind the baes are clearly visible (red boxes 1, 2 and 3). . . . . . . . . . . . 58 LIST OF FIGURES xix III.9 POD results for the velocity elds in rotating zone measured in the moving frame. a) Mode spectrum: one mode dominates the variance ranking; b) mode components associated to the three modes, circles: rst mode, x's second mode and rhomboids third mode (each data set is normalized by the variance corresponding to each mode). The modal component of the rst ranking mode is time independent and the second and third modes both oscillate with a frequency of 10 Hz. These oscillations are related to the passage of the baes observed from the moving frame. c) The rst eigenvector in the plane Y Z three degrees behind one of the turbine blades reveal the presence of the trailing vortices (red box 1), d) The rst eigenvector in the plane Y Z six degrees behind one of the turbine blades reveals the radial displacement of the trailing vortices (red box 2). . . . . 61 III.10a) POD eigenvalue spectrum for the rst twenty modes and three dierent time spans: the rst three modes represents almost the 100 % of the total system variance. It is also evident a correlation between pairs in the upper modes until the eleventh one. b) Cumulated variance POD spectrum for the time span of 386 snapshots. The cumulated variance until the third mode already represents the 99.9 % of the total system variance. . . . . . . . . . . . . . . . . . . . . . . . 64 III.11Graphical representation of the rst three modes identied with the global treatment. For clarity, results are illustrated during 2 seconds only. a) Normalized time variation of the rst and second mode. The circles and the crosses correspond to the modal component a 1 and a 2 respectively. The dotted and continuous lines correspond to their tting by continuous sinusoidal functions. b) The circular conguration reveals a phase shift of π 2 between the rst and second modes. c) Normalized time variation of the third mode. The solid circle corresponds to the modal component a 3 . The dotted line corresponds to the sinusoidal tting of the third mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 III.12Probability Density Functions for a 1 (t), a 2 (t) a 3 (t). . . . . . . . . . . . . . . . 66 III.13Lissajous patterns obtained from the temporal variation of the rst to tenth modes. 67 xx LIST OF FIGURES III.14a) Reconstruction of the velocity using the rst three modes. From the left to the right side the most important features of the eld are marked by numbered dot rectangles: Upward and downward axial pumping ow (red boxes #1), vortices generated by the passage of the blade (red boxes #2), and radial ow out of the impeller-sweeping zone (red box #3), b) Velocity eld from the CFD data in the same snapshot of the POD reconstruction. . . . . . . . . . . . . . . . . . . . . 70 III.15Flow reconstruction using the third mode. a) In an angular sector of 3.5 degrees visualized in a vertical XZ plane. The circulation loops are marked with red boxes 1 and 2. b) Coherent structures made of cylindrical eddies rotating clockwise are clearly observed behind the blades reveals the presence of clockwise vortices next to the baes as marked by red boxes 3 and 4. c) and d) shows a zoom of two of the coherent structures behind the baes. . . . . . . . . . . . . . . . . . . . . 72 III.16Maximum average relative error for 2, 3, 5 and 10 POD modes. . . . . . . . . . 74 III.17Bar charts of local and instantaneous relative error for three dierent snapshots using 5 and 10 POD modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 III.18Cumulative percentage of the number of cells as a function of the relative error calculated by equation III.5, using either 5 and 10 POD modes. An increase in the number of POD modes leads to an upward shift of the curve. The vertical red line indicates 3% relative error. . . . . . . . . . . . . . . . . . . . . . . . . 77 III.19a) Localized error bigger than 5 % for the 1 th , 23 rd and 386 th snapshots from a 5 POD modes reconstruction. b) Localized error bigger than 3 % for the 1 th , 23 rd and 386 th snapshots from a 10 POD modes reconstruction. . . . . . . . . . . . 78 III.20Summary of the dynamical stability results for the spectra corresponding to 250, 271 and 386 snapshots. Obviously, the parameter λ DM D should not be confused with the parameter λ P OD . The former is a complex number containing the DMD frequency and the DMD growth rate, the latter is a measure of the covariance calculated in the POD method. . . . . . . . . . . 82 III.21Summary of the growth stability results for the spectra corresponding to 250, 271 and 386 snapshots. . . . . . . . . . . . . . . . . . . . . . . . . . . 83 LIST OF FIGURES xxi III.22Spectral ranking results for σ j DM D corresponding to the DMD modes for the 386 snapshots case typical ranking criteria: a) Modes with a very negative sigma value carry very small amounts of energy compared to modes whose sigma value is approximately zero. , b) Approach in the neighborhood of σ j DM D = 0s -1 . Both positive and negative sigma values have a very small energy content compared to modes whose sigma values are approximately zero. . . . . . . . . . . . . . . 85 III.23Spectral ranking results for σ j DM D corresponding to the DMD modes for the 386 snapshots case Tissot ranking criteria: a) Modes with a very negative sigma value carry very small amounts of energy compared to modes whose sigma value is approximately zero. , b) Approach in the neighborhood of σ j DM D = 0s -1 . Both positive and negative sigma values have a very small energy content compared to modes whose sigma values are approximately zero. . . . . . . . . . . . . . . 86 III.24Frequency spectral ranking results corresponding to the DMD modes for the 386 snapshots case: a) Typical ranking criteria, b) Alternative ranking criteria. The black boxes to the left of 0 Hz represent the "mirror" pairs of 2.5, 7.5, 10.0, 12.5, and 15.0 Hz. The horizontal dashed line represents a normalized energy value of 10 -6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 III.25POD and DMD (alternative criteria) spectrum mode for the 386 snapshots case. The ranking of the modes is identical for both methodologies. . . . . . . . . . 89 III.26Graphical representation of the top three modes (temporal part) identied with the DMD method. For clarity, results are illustrated during 2 seconds only. Real (red rhomboids) and imaginary (blue circles) parts of the 83 rd ,112 th , and 113 rd DMD modal coecients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 III.27Comparison between the POD reconstruction and the corresponding DMD reconstruction for the mean ow using: a) 3 rd POD eigenvector eld (6.5 degrees), b) 83 rd DMD vector eld (6.5 degrees), c) 3 rd POD eigenvector eld, d) 83 rd DMD vector eld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 III.28Reconstruction of the velocity using the real parts of the 83,112, and 113 DMD modes. Important features (also obtained by the POD method) are identied: Upward and downward axial pumping ow, trailing vortices generated by the passage of the blade, and the radial ow out the impeller-sweeping zone, b) Velocity eld from the CFD data in the same snapshot of the DMD reconstruction. 92 III.29Maximum average relative error for POD and DMD for a time span of 386 snapshots or 11 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 III.30Bar charts of local and instantaneous relative error for four dierent snapshots using POD and DMD procedures. . . . . . . . . . . . . . . . . . . . . . . . . 94 III.31Cumulative error for the POD and DMD in four dierent snapshots. . . . . . . 95 III.32a) Localized error bigger than 5 % for the 1 th , 23 rd and 300 th snapshots from a 5 DMD modes reconstruction. b) Localized error bigger than 3 % for the 1 th , 23 rd and 300 th snapshots from a 11 DMD modes reconstruction. . . . . . . . . . . . 96 III.33Sampling regions for the low-frequency analysis. . . . . . . . . . . . . . . . . . 99 III.34Sampling time for the low-frequency analysis. . . . . . . . . . . . . . . . . . . 99 III.35DMD dynamical stability plots for the subsets: 1 min, 1.5 mins and 2.3 mins corresponding to the region I. . . . . . . . . . . . . . . . . . . . . . . . . . . 100 III.36Summary of the growth stability results for the three time spans corresponding to the region I. In the blue dotted line rectangle on nds the smallest growth rate values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 III.37Energy spectrum vrs frequency for the subsets: 1 min, 1.5 mins and 2.3 mins corresponding to the region I. . . . . . . . . . . . . . . . . . . . . . . . . . . 102 III.38Energy spectrum vrs frequency for the regions I,II and III using a time span of 2.3 mins (around 4700 snapshots). . . . . . . . . . . . . . . . . . . . . . . . 103

  components for the 13 th and 14 th POD mode. . . . . . . . . . . . . . . 105 III.41Modal component for the 7 th POD mode. . . . . . . . . . . . . . . . . . . . . 105IV.1 a) POD eigenvalue spectrum for the liquid phase: the rst three modes represents around the 80 % of the total kinetic energy of the system. POD eigenvalue spectrum for the gas phase: the rst three modes represents around the 70 % of the total kinetic energy of the system. b) Cumulated normalized energy POD spectrum for both phases show that a lot of modes have to be taken to capture the 100 % of the total system kinetic energy. . . . . . . . . . . . . . . . . . . . 108 xxiii IV.2 Normalized mode components associated to the rst three POD modes for the liquid phase. b) Relation between the rst and the second POD components modes, the circular conguration evidences la periodic behavior between these parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 IV.3 Eigenvectors associated to the rst three POD modes for the liquid phase. The red loops represent the sense of circulation for the organized structures in each eigenvector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 IV.4 Eigenvectors associated to the fourth, fth and sixth POD modes for the liquid phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 IV.5 Normalized mode components associated to the rst three POD modes for the gas phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 IV.6 Normalized mode components associated to the fourth, fth and sixth POD modes for the gas phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 IV.7 Eigenvectors associated to the rst three POD modes for the gas phase. . . . . . 113 IV.8 Eigenvectors associated to the fourth, fth and sixth POD modes for the gas phase.114 IV.9 a) POD reconstruction using 600 POD modes for the liquid phase in the snapshots 1 th , 49 th , 87 th , and 133 rd . b) CFD velocity eld for the liquid phase in the snapshots 1 th , 49 th , 87 th , and 133 rd . . . . . . . . . . . . . . . . . . . . . . . . 115 IV.10a) POD reconstruction using 600 POD modes for the gas phase in the snapshots 1 th , 49 th , 87 th , and 133 rd . b) CFD velocity eld for the gas phase in the snapshots 1 th , 49 th , 87 th , and 133 rd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 IV.11a) Summary of the DMD dynamical stability plots for the subsets: 1000,1500, 2996 snapshots for the liquid phase. b) Summary of the DMD dynamical stability plots for the subsets: 1000,1500, 2996 snapshots for the gas phase . . . . . . . . 119 IV.12a) Summary of the growth stability results for the subsets: 1000,1500, 2996 snapshots for the liquid phase. b) Summary of the growth stability results for the subsets: 1000,1500, 2996 snapshots for the gas phase . . . . . . . . . . . . . 120 IV.13a) Frequency spectral ranking results corresponding to the DMD modes for the liquid phase using 1000,1500 and 2996 snapshots. b) Frequency spectral ranking results corresponding to the DMD modes for the gas phase using 1000,1500 and IV.14DMD vector 507: real and imaginary parts for the liquid phase. The red loops represent the sense of circulation for the organized structures in the DMD vector. 122 IV.15DMD vector 563: real and imaginary parts for the liquid phase. . . . . . . . . . 122 IV.16DMD vectors 607, 605: real and imaginary parts for the gas phase. . . . . . . . 123 IV.17POD reconstruction of the velocity eld for the liquid phase using the rst three modes in the snapshots 1, 49, 87, and 133. Reconstruction with 3 POD modes allows obtaining coherent structures whose conguration is roughly similar to the CFD results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 IV.18POD reconstruction of the velocity eld for the liquid phase using the rst six modes in the snapshots 1, 49, 87, and 133. Even taking only 6 POD modes the ow reconstruction is relatively good when compared to the CFD results shown earlier in this section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 V.1 Exit owrates between adjacent compartments and construction of the matrix of owrates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 V.2 Identication of the destination compartment from the know cumulated probability CP and the random number β. For any compartment, neighboring compartments are indexed form 1 to 6 (x-axis) . . . . . . . . . . . . 131 V.3 Parity plot of the eective transition probabilities vs. expected probabilities. About 9 Million exit events were collected. . . . . . . . . . . . . . . . 133 V.4 Scalar mixing validation : Mean scalar concentration (left) and normalized concentration proles in ve randomly seleted compartments. . . . . . . . . 134 V.5 Local particle density n p (i, t)/V (i) in ve randomly selected compartments. a. 5 10 4 particles, b. 10 5 particles, c. 5 10 5 particles, d. 10 6 particles. . . . . 135 V.6 Segregation index for scalar mixing and particle mixing. The residual variance evolves with the inverse of the total number of particles. . . . . . . . 136 V.7 Normalized scalar and particle densities in ve randomly selected compartments. Continuous line : scalar mixing, symbols : particle mixing (10 6 particles), 500 compartments. . . . . . . . . . . . . . . . . . . . . . . . . . 137 V.8 Relative standard deviation in local particle number for all compartments. The parity curve represents equation (V.24). x-axis: measured values, yaxis : expected values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 xxv V.9 First order reaction calculation in a 500 compartments, 70 L reactor. In ve randomly selected compartments, the concentration proles computed using particle information (open circles) matches the usual homogeneous approach (continuous line). Dashed lines stand for the normalized particle density which remain constant over time. . . . . . . . . . . . . . . . . . . . A.1 Spectral ranking results for σ j corresponding to the DMD modes for the 250 snapshots case typical ranking criteria: a) Modes with a very negative sigma value carry very small amounts of energy compared to modes whose sigma value is approximately zero, b) Approach in the neighborhood of σ j = 0s -1 . Both positive and negative sigma values have a very small energy content compared to modes whose sigma values are approximately zero. . . . . . . . . . . . . . . A.2 Spectral ranking results for σ j corresponding to the DMD modes for the 250 snapshots case alternative ranking criteria: a) Modes with a very negative sigma value carry very small amounts of energy compared to modes whose sigma value is approximately zero, b) Approach in the neighborhood of σ j = 0s -1 . Both positive and negative sigma values have a very small energy content compared to modes whose sigma values are approximately zero. . . . . . . . . . . . . . . A.3 Frequency spectral ranking results corresponding to the DMD modes for the 250 snapshots case: a) Typical ranking criteria, b) Alternative ranking criteria. . . . A.4 POD and DMD (alternative criteria) spectrum mode for the 250 snapshots case. The ranking of the modes is identical for both methodologies. . . . . . . . . . A.5 Spectral ranking results for σ j corresponding to the DMD modes for the 271 snapshots case alternative ranking criteria: a) Modes with a very negative sigma value carry very small amounts of energy compared to modes whose sigma value is approximately zero, b) Approach in the neighborhood of σ j = 0s -1 . Both positive and negative sigma values have a very small energy content compared to modes whose sigma values are approximately zero. . . . . . . . . . . . . . . A.6 Frequency spectral ranking results corresponding to the DMD modes for the 271 snapshots case: a) Typical ranking criteria, b) Alternative ranking criteria. . . . A.7 POD and DMD (alternative criteria) spectrum mode for the 271 snapshots case.

II. 1

 1 Abbreviations

  Figure I.1: a) Quadrangular segmentation of an annular domain. b)Triangular segmentation of an annular domain.

Figure I. 2

 2 Figure I.2 shows the most important steps involved in the implementation of a CFD experiment. The dotted line shows a possible route to follow when key criteria such as maximum numerical residuals and conservation laws are not fullled in the results obtained.

Figure

  Figure I.2: c of a CFD experiment.

Figure I. 3 :

 3 Figure I.3: Schematic view of Richardson's energy cascade in turbulence.

Figure I. 4 :

 4 Figure I.4: Schematic representation of the trailing vortices in the blade of a Rushtone turbine. The structures start on the internal side of the blade and afterward travel outward.

  Figure I.5 shows the typical conguration of an open agitation tank activated by a Rushton turbine. The device shown has a at bottom and also four vertical elements called baes whose function is to promote mixing and at the same time avoid parabolic velocity proles.

Figure I. 5 :

 5 Figure I.5: Typical conguration of a baed stirred tank activated by a Rushton turbine.

  zone and a high inuence zone of the impellers (see Figure I.6).

Figure I. 6 :

 6 Figure I.6: Conguration of two mesh zones for a baed stirred tank activated by a Rushton turbine. The rotative zone contains the turbine, the xed zone the rest of the domain.

  in numerous industrial chemical processes. The typical conguration of this type of device is very simple as shown in the Figure I.7. The vessel has a cylindrical or prismatic tank containing a liquid material, typically water. In the lower part of the vessel, there are orices through which gas is injected at a given rate and velocity.

Figure I. 7 :

 7 Figure I.7: Typical conguration of a bubble column vessel.

Figure I. 8 :

 8 Figure I.8: Typical conguration of a bubble column vessel.

(

  [START_REF] Mikhaylov | Reconstruction of large-scale ow structures in a stirred tank from limited sensor data[END_REF] performed the reconstruction of coherent structures just in the zone around the impeller blades at low Reynolds number through the POD treatment of the velocity eld expressed in a rotating frame. Finally, in both investigations the presence of baes in the simulated stirred tanks was not considered. The incorporation of these elements is of great importance due to their contribution to the mean ow (circulation loops) and their decisive role in mixing phenomena. To sum up, running time-resolved CFD (Large Eddy Simulation or U-RANS) and applying Reduced Order Methods such as POD (Proper Orthogonal Decomposition or Karhunen-Loève decomposition) is appealing because it generates a time-varying velocity elds at a moderate expense while preserving the spatial resolution[START_REF] Du | POD reduced-order unstructured mesh modeling applied to 2D and 3D uid ow[END_REF] de Lamotte et al., 2018b). This opens the route to both rened and cost-eective description of the unsteady velocity elds. From this brief review of the existing literature on the subject, it appears that the POD analysis of the ow eld in a baed stirred tank, in the fully turbulent regime (at a Reynolds number > 10

(

  [START_REF] Kou | An improved criterion to select dominant modes from dynamic mode decomposition[END_REF] used DMD to predict the ow regime of ow past a cylinder at the transient state of Re = 60. Additionally, this work studied the transient ow in a NACA0012 airfoil (Reynolds numbers of 60 and 3.0 x 10 6 ).[START_REF] Clainche | Higher order dynamic mode decomposition to identify and extrapolate ow patterns[END_REF]) simulated a three-dimensional ow around a circular cylinder (Reynolds number = 220). The resulting data allowed the implementation of a DMD reduction order model. This research proposes using this ROM to nd ow patterns in complex ows, highlighting the advantages of reducing the computational cost of numerical simulations or the necessary amount of data collected in experiments.Likewise, in (de Lamotte et al., 2018a; de Lamotte et al., 2018b) DMD allowed extracting the spatial structures according to their temporal behavior (frequency) from PIV data and CFD simulations in a stirred tank with baes. In these outstanding works, a 2D plane analysis led to compare the numerical simulation results with those of the PIV experiments. In[START_REF] Weheliye | Mode decomposition and Lagrangian structures of the ow dynamics in orbitally shaken bioreactors[END_REF] PIV data was used as input to the DMD technique to assess the ow dynamics in an orbital shaken bioreactor of cylindrical geometry and at bottom. The results allowed to assess the ow dynamics in dierent operating conditions of the vessel. Also[START_REF] Wu | Proper orthogonal decomposition and dynamic mode decomposition of jet in channel crossow[END_REF] used this reduced order methodology to generate a reduced model of a low-moment laminar jet discharged in a laminar channel cross-ow through a circular orice (Reynolds number = 3333). More recently,[START_REF] Jin | PIV experimental study on ow structure and dynamics of square stirred tank using modal decomposition[END_REF] investigated the ow structure in a square stirred tank without baes and with a Rushton impeller (RT) using (PIV) technique. In this publication, the ROM's technique was practical for analyzing velocity elds, ow structure, and dynamic information out of the impeller area. It was possible to identify and extract large-scale ow structures such as the mean ow and wake structure generated by the periodic jet in blade passage.

  Our numerical data were generated from an unsteady RANS CFD simulation in a 70 L stirred tank with baes developed by[START_REF] Delafosse | Analyse et étude numérique des eets de mélange dans un bioréacteur[END_REF]. The geometry and conguration of the at-bottom stirred tank equipped with a Rushton turbine and four baes are presented in gure II.1.

Figure

  Figure II.1: Schematic representation of the baed stirred tank simulated.

  two zones: the xed zone contains the walls, baes, the major part of the shaft, and the volume outside the rotating zone RZ, the latter is a cylindrical domain, which contains the impeller and a small portion of the shaft, as depicted in Figure II.1. A structured mesh made of hexahedrons was built in the RZ zone (see Figure II.2).

Figure II. 2 :

 2 Figure II.2: Mesh FLUENT views of the most important parts of the simulated tank: a) Top surface of the tank, b) External walls of the tank, c) Shaft and Rushton turbine, d) Shaft and rotating zone of the simulation domain.

Figure

  Figure II.3 shows that the torque on the turbine and the turbine shaft is practically identical to the torque generated on the external walls of the tank

  Figure II.5 (Delafosse, 2008) gives the experimental and numerical (RANS-k -ε) vertical proles of the mean velocity for three dierent positions close to the Rushton turbine. It is observed that the PIV experimental data agree quite well with the results of the numerical simulation.

Figure II. 5 :

 5 Figure II.5: Experimental and numerical vertical proles (RANS-k -ε) of the three components of the mean velocity in the turbine jet.

Figure II. 6 :

 6 Figure II.6: Typical conguration of a bubble column vessel.

  . A good match between the local vertical velocity values of the PIV and those predicted by the R ij +BIT (Reynolds ij model + Bubble Turbulence Induced) model is conrmed (see section I.1.2.2).

Figure

  Figure II.7: Comparison between PIV values and simulated values for the window indicated by the red rectangle shown (0 ≤ x ≤ 0.17 cm, 0.6 ≤ z ≤ 0.7 cm).

III. 1

 1 Proper Orthogonal Decomposition (POD) III.1.1 Sampling and data processing After reaching the stationary ow condition 1 (139 complete revolutions), the velocity eld data was collected (see chapter IV).The sampling of the simulated data consists of a set of 3D velocity elds (snapshots) collected in the entire uid domain (1,129,140 cells). A total physical time of 11 s, representing around 28 complete rounds of the turbine and requiring 26 days of computational calculation, was spanned. During this time lapse, 386 snapshots were taken every 53 computational time steps. Consequently, the time interval between two samples is ∆t st = 0.03s (∆t st = 53∆t CF D ). As will be shown later, the application of POD to 250, 271, and 386 snapshots will allow us to assess the sensitivity of the results. Figure III.1 provides an overview of the sampling data chosen in the present study.

Figure

  Figure III.1: Sampling data sets collected for apply POD. Each rectangle represents a snapshot of registered velocity.

Figure

  Figure III.2: Schematic representation of the data sampling frequency. The continuous plot is given by y = sin(5πt).

Figure

  Figure III.3: 2D schematic representation for the velocity vector expressed in the inertial and rotating frames. The inertial frame is equivalent to the laboratory frame of reference. The inertial and rotating frame share the same origin. The rotating frame rotates around the z -axis with angular velocity ω. In that frame, the coordinates of the rotating cells are time independent.

Figure

  Figure III.4 depicts schematically the steps involved in the typical and global treatments.

Figure

  Figure III.4: Schematic representation of the typical and global treatments.

Figure III. 5 :

 5 Figure III.5: Schematic representation for the n th snapshot of the 3D velocity eld ⃗ U (n) IF (CI).

Figure

  Figure III.6: One-column array representation for the n th snapshot of the 3D velocity eld ⃗ U (n) IF (CI).

  Figure III.9: POD results for the velocity elds in rotating zone measured in the moving frame. a) Mode spectrum: one mode dominates the variance ranking; b) mode components associated to the three modes, circles: rst mode, x's second mode and rhomboids third mode (each data set is normalized by the variance corresponding to each mode). The modal component of the rst ranking mode is time independent and the second and third modes both oscillate with a frequency of 10 Hz. These oscillations are related to the passage of the baes observed from the moving frame. c) The rst eigenvector in the plane Y Z three degrees behind one of the turbine blades reveal the presence of the trailing vortices (red box 1), d) The rst eigenvector in the plane Y Z six degrees behind one of the turbine blades reveals the radial displacement of the trailing vortices (red box 2).

Figure

  Figure III.10 shows the mode spectrum of the POD obtained by the global treatment. The rst three modes have 99.99 % of the total system variance (Figure III.10 b). The rst and second modes share the same variance content, while the third mode has slightly lower variance. In addition, from modes 4 to 11, pairs of modes of equal variance are observed. In general, one or more pairs of modes with the same variance content are typical of coherent structures propagating in space (The phase dierence between the time coecients and the spatially shifted conguration between modes leads to the movement of the structures). Hence, this decomposition evidences a strong spatial organization of the ow in the form of multiple coherent structures with decreasing variance content.

Figure

  Figure III.11 depicts the time variation of the normalized amplitude factors associated with the rst and second modes over 28 impeller turns (time span 386 snapshots). The temporal organization of these modes is clear; their behavior is oscillatory, with an identical period (0.4 s equal to the rotation period of the turbine) and amplitude. In addition, the phase shift between these modes is π 2 as can be seen from the circular conguration of Figure III.11 b. Similarly, the amplitude factor of the third mode, presented in III.11 c,

  Figure III.10: a) POD eigenvalue spectrum for the rst twenty modes and three dierent time spans: the rst three modes represents almost the 100 % of the total system variance. It is also evident a correlation between pairs in the upper modes until the eleventh one. b) Cumulated variance POD spectrum for the time span of 386 snapshots. The cumulated variance until the third mode already represents the 99.9 % of the total system variance.

  Figure III.12: Probability Density Functions for a 1 (t), a 2 (t) a 3 (t).

  Figure III.14 b shows the velocity eld obtained from the CFD simulation at the same instant of the POD reconstruction (Figure III.14 a). A simple inspection of both gures shows a very close resemblance between the reconstructed eld and the eld generated by the simulator.

  Figure III.14: a) Reconstruction of the velocity using the rst three modes. From the left to the right side the most important features of the eld are marked by numbered dot rectangles: Upward and downward axial pumping ow (red boxes #1), vortices generated by the passage of the blade (red boxes #2), and radial ow out of the impeller-sweeping zone (red box #3), b) Velocity eld from the CFD data in the same snapshot of the POD reconstruction.

  Figure III.15: Flow reconstruction using the third mode. a) In an angular sector of 3.5 degrees visualized in a vertical XZ plane. The circulation loops are marked with red boxes 1 and 2. b) Coherent structures made of cylindrical eddies rotating clockwise are clearly observed behind the blades reveals the presence of clockwise vortices next to the baes as marked by red boxes 3 and 4. c) and d) shows a zoom of two of the coherent structures behind the baes.

Figure

  Figure III.16: Maximum average relative error for 2, 3, 5 and 10 POD modes.

Figure

  Figure III.17: Bar charts of local and instantaneous relative error for three dierent snapshots using 5 and 10 POD modes.

  cumulated percentage of the number of cells for the reconstructions of 5

Figure

  Figure III.18: Cumulative percentage of the number of cells as a function of the relative error calculated by equation III.5, using either 5 and 10 POD modes. An increase in the number of POD modes leads to an upward shift of the curve. The vertical red line indicates 3% relative error.

  Figure III.19: a) Localized error bigger than 5 % for the 1 th , 23 rd and 386 th snapshots from a 5 POD modes reconstruction. b) Localized error bigger than 3 % for the 1 th , 23 rd and 386 th snapshots from a 10 POD modes reconstruction.

  Figure III.20: Summary of the dynamical stability results for the spectra corresponding to 250, 271 and 386 snapshots. Obviously, the parameter λ DM D should not be confused with the parameter λ P OD . The former is a complex number containing the DMD frequency and the DMD growth rate, the latter is a measure of the covariance calculated in the POD method.

Figure

  Figure III.21 provides the values of σ j DM D for the three cases analyzed. Most point values are close to the horizontal stability line (σ j

  Figure III.21: Summary of the growth stability results for the spectra corresponding to 250, 271 and 386 snapshots.

Figure

  Figure III.22 consolidates the energy information of the spectrum as a function of the growth rate for the 386 snapshots data set. The normalized energy values presented are those corresponding to the typical ranking criteria. Figure III.22 a reveals that modeswith very negative growth rate values carry a very small quantity of energy on the order of 10 -7 or less. This behavior implies that these elements of the spectrum quickly vanish and contribute slightly to the total energy content of the vector eld. Furthermore, one observes very small negative (to the left of zero) and very small positive (to the right of

  Figure III.24: Frequency spectral ranking results corresponding to the DMD modes for the 386 snapshots case: a) Typical ranking criteria, b) Alternative ranking criteria. The black boxes to the left of 0 Hz represent the "mirror" pairs of 2.5, 7.5, 10.0, 12.5, and 15.0 Hz. The horizontal dashed line represents a normalized energy value of 10 -6 .

  representation of the reconstructed ow using the three rst DMD modes We will start the section by discussing the time parameters corresponding to the top three modes DMD 112, 113, and 83 (see table III.6). Figure III.26 provides a graphical representation of these leading DMD modes for 2 seconds. It is evident the oscillatory behavior presented in the real parts (red rhomboids) and the imaginary parts (blue circles) of the modal coecients 112 and 113. Additionally, the real and imaginary part of each mode presents a phase of π 2 (see the Figure A.8 in the appendix A). These kinds of temporal structures are associated with the representation/reconstruction of traveling waves. Regarding mode 83 one observes that the real part is constant and non-zero while the imaginary part is zero. Therefore mode 83 does not oscillate and has a constant behavior in time. Similar results to those of modes 112, 113, and 83 were already reported in an experimental work in which the DMD method was also used in a stirred tank (de Lamotte et al., 2018a).

Figure

  Figure III.26: Graphical representation of the top three modes (temporal part) identied with the DMD method. For clarity, results are illustrated during 2 seconds only. Real (red rhomboids) and imaginary (blue circles) parts of the 83 rd ,112 th , and 113 rd DMD modal coecients.

Figure

  Figure III.27 depicts a comparison (planes XZ and XY) between the POD reconstruction (3 rd POD mode) and the corresponding DMD reconstruction (83 rd DMD mode) of the mean ow. A simple inspection of the information shown allows the equivalence between the two modes to be conrmed. All structural features obtained by POD reconstruction

  Figure III.29 provides the results obtained for both POD and DMD. As can be seen, the resulting match is perfect for the E M AR .

Figure

  Figure III.29: Maximum average relative error for POD and DMD for a time span of 386 snapshots or 11 s.

Figure

  Figure III.30: Bar charts of local and instantaneous relative error for four dierent snapshots using POD and DMD procedures.

Figure

  Figure III.31 depicts the results of the cumulative percentage of the number of cells versus the instantaneous localized error for the POD and the DMD. The error distribution for both methods is essentially identical: a 5 DMD modes reconstruction leads to 90.5% of the mesh with an error less than that indicated by the dotted red line (3%). By using 11 DMD modes, the results improve as 98.5% of the cells have an error of less than 3%.

Figure

  Figure III.31: Cumulative error for the POD and DMD in four dierent snapshots.

  Figure III.33. It is observed that zones close to the turbine (regions I and II) are taken into account in the analysis. This is important to highlight because the experimental works report the presence of IM's just in the upper zone or region I and/or below the turbine. Regarding time sampling, a time span of 2.3 mins with a time step of 0.03 s was analyzed, which represents approximately 4670 snapshots (see Figure III.34). This time step is used because it is intended in the future to use the low-frequency results together with the high-frequency results obtained in the previous section. The sampling process started at 32 seconds of simulation after the turbine had rotated 80 times (see the graph of physical convergence of the numerical results in the chapter II).

Figure

  Figure III.33: Sampling regions for the low-frequency analysis.

Figure

  Figure III.34: Sampling time for the low-frequency analysis.

Figure

  Figure III.35 shows the real and imaginary values for the studied λ DM D spectra. Oneobserves that all values lie on the unit circle of the complex plane for the three subsets processed. As explained in section III.2.1.2 this implies dynamic stability for each resulting DMD spectra. This nding means good numerical stability in the database.

Figure

  Figure III.35: DMD dynamical stability plots for the subsets: 1 min, 1.5 mins and 2.3 mins corresponding to the region I.

  exposed in section III.2.1.3. This choice is based on the fact that such a criterion considers the energy contribution over time of each mode. Two dierent scenarios are analyzed to visualize the energy spectrum according to the number of modes and the three regions of the Figure III.33. In the rst scenario, the modal spectra for

Figure

  Figure III.36: Summary of the growth stability results for the three time spans corresponding to the region I. In the blue dotted line rectangle on nds the smallest growth rate values.

Figure

  Figure III.37 summarizes the corresponding results. We note several interesting aspects.First, a change in the time span results in a change in the total energy of the spectrum.

Figure

  Figure III.38: Energy spectrum vrs frequency for the regions I,II and III using a time span of 2.3 mins (around 4700 snapshots).

Figure

  Figure III.40: Modal components for the 13 th and 14 th POD mode.

Figure

  Figure III.41: Modal component for the 7 th POD mode.

Figure

  Figure IV.2 provides the modal components for the liquid phase associated to the rst POD modes.

  Figure IV.2: Normalized mode components associated to the rst three POD modes for the liquid phase. b) Relation between the rst and the second POD components modes, the circular conguration evidences la periodic behavior between these parameters.

Figure IV. 3 :

 3 Figure IV.3: Eigenvectors associated to the rst three POD modes for the liquid phase. The red loops represent the sense of circulation for the organized structures in each eigenvector.

Figure IV. 5 :

 5 Figure IV.5: Normalized mode components associated to the rst three POD modes for the gas phase.

Figure IV. 6 :

 6 Figure IV.6: Normalized mode components associated to the fourth, fth and sixth POD modes for the gas phase.

Figure IV. 7 :

 7 Figure IV.7: Eigenvectors associated to the rst three POD modes for the gas phase.

Figure

  Figure IV.8: Eigenvectors associated to the fourth, fth and sixth POD modes for the gas phase.

  Figure IV.9: a) POD reconstruction using 600 POD modes for the liquid phase in the snapshots 1 th , 49 th , 87 th , and 133 rd . b) CFD velocity eld for the liquid phase in the snapshots 1 th , 49 th , 87 th , and 133 rd .

  Figure IV.11 provides the DMD dynamical stability plots for the subsets: 1000,1500, 2996 snapshots for the two phases.

  Figure IV.11: a) Summary of the DMD dynamical stability plots for the subsets: 1000,1500, 2996 snapshots for the liquid phase. b) Summary of the DMD dynamical stability plots for the subsets: 1000,1500, 2996 snapshots for the gas phase

  Figure IV.12: a) Summary of the growth stability results for the subsets: 1000,1500, 2996 snapshots for the liquid phase. b) Summary of the growth stability results for the subsets: 1000,1500, 2996 snapshots for the gas phase

Figure

  Figure IV.14: DMD vector 507: real and imaginary parts for the liquid phase. The red loops represent the sense of circulation for the organized structures in the DMD vector.

Figure

  Figure IV.15: DMD vector 563: real and imaginary parts for the liquid phase.

Figure

  Figure IV.17: POD reconstruction of the velocity eld for the liquid phase using the rst three modes in the snapshots 1, 49, 87, and 133. Reconstruction with 3 POD modes allows obtaining coherent structures whose conguration is roughly similar to the CFD results.

Figure

  Figure IV.18: POD reconstruction of the velocity eld for the liquid phase using the rst six modes in the snapshots 1, 49, 87, and 133. Even taking only 6 POD modes the ow reconstruction is relatively good when compared to the CFD results shown earlier in this section.

Figure V. 1 :

 1 Figure V.1: Exit owrates between adjacent compartments and construction of the matrix of owrates

Figure V. 2 :

 2 Figure V.2: Identication of the destination compartment from the know cumulated probability CP and the random number β. For any compartment, neighboring compartments are indexed form 1 to 6 (x-axis)

  Figure V.3 indicates that the algorithm is correctly implemented as the eective transition probabilities (based on the number of displacements) closely match the expected values based on the mean owrates.

Figure V. 4 :

 4 Figure V.4: Scalar mixing validation : Mean scalar concentration (left) and normalized concentration proles in ve randomly seleted compartments.

Figure V. 5 :

 5 Figure V.5: Local particle density n p (i, t)/V (i) in ve randomly selected compartments. a. 5 10 4 particles, b. 10 5 particles, c. 5 10 5 particles, d. 10 6 particles.

Figure V. 6 :

 6 Figure V.6: Segregation index for scalar mixing and particle mixing. The residual variance evolves with the inverse of the total number of particles.

Figure V. 7 :

 7 Figure V.7: Normalized scalar and particle densities in ve randomly selected compartments. Continuous line : scalar mixing, symbols : particle mixing (10 6 particles), 500 compartments.

  Figure V.8: Relative standard deviation in local particle number for all compartments. The parity curve represents equation (V.24). x-axis: measured values, y-axis : expected values.

  Figure V.9: First order reaction calculation in a 500 compartments, 70 L reactor. In ve randomly selected compartments, the concentration proles computed using particle information (open circles) matches the usual homogeneous approach (continuous line).

  Figure A.2: Spectral ranking results for σ j corresponding to the DMD modes for the 250 snapshots case alternative ranking criteria: a) Modes with a very negative sigma value carry very small amounts of energy compared to modes whose sigma value is approximately zero, b) Approach in the neighborhood of σ j = 0s -1 . Both positive and negative sigma values have a very small energy content compared to modes whose sigma values are approximately zero.

Figure A. 4 :

 4 Figure A.4: POD and DMD (alternative criteria) spectrum mode for the 250 snapshots case. The ranking of the modes is identical for both methodologies.

Figure

  Figure A.7: POD and DMD (alternative criteria) spectrum mode for the 271 snapshots case. The ranking of the modes is identical for both methodologies.

Figure A. 8 :

 8 Figure A.8: Imaginary part vrs Real part for the DMD coecients 112 and 113 (386 snapshots case). The circular conguration implies that the imaginary and real parts are shifted by π 2 .

Figure A. 9 :

 9 Figure A.9: FFT for the 1 st POD mode

Figure A. 11 :

 11 Figure A.11: FFT for the 3 rd POD mode

Figure A. 13 :

 13 Figure A.13: FFT for the 5 th POD mode

Figure A. 15 :

 15 Figure A.15: FFT for the 7 th POD mode

Figure A. 17 :

 17 Figure A.17: FFT for the 9 th POD mode

Figure

  Figure A.19: FFT for the 11 th POD mode

Figure A. 21 :

 21 Figure A.21: FFT for the 13 th POD mode

Figure A. 23 :

 23 Figure A.23: FFT for the 15 th POD mode

Figure A. 25 :

 25 Figure A.25: FFT for the 17 th POD mode

Figure A. 27 :

 27 Figure A.27: FFT for the 19 th POD mode

TABLE I .

 I 1: Values of the constants for the k -ε model.

	C ε1 C ε2	C µ	σ k	σ ε
	1.44	1.92	0.09	1.00	1.30
	I.1.2.2 3D two-phase uid ow: Bubbly ows

Table I

 I

.2 and table I.3 provide a basic description of each modeling approach

[START_REF] Joshi | CFD simulation of stirred tanks: Comparison of turbulence models. Part I: Radial ow impellers[END_REF]

.

TABLE I

 I

	Steady state approaches Details	The boundary conditions	(velocities, kinetic and dissipation turbulence)	are measured experimentally. Deeply dependence on	the experimental data	availability.
	Model			Impeller boundary condition	(IBC)

.2: Steady state approaches for the impeller rotation and baes modelling.

TABLE I

 I

	Unsteady state approaches	Ideal for computing parallelisation.
	Model	Lattice-Boltzmann	LES

.3: Unsteady state approaches for the impeller rotation and baes modelling.

  ;[START_REF] Gabelle | Experimental determination of the shear rate in a stirred tank with a non-newtonian uid: Carbopol[END_REF] de Lamotte et al., 2018b;[START_REF] Janiga | Large-eddy simulation and 3D proper orthogonal decomposition of the hydrodynamics in a stirred tank[END_REF][START_REF] Fernandes Del Pozo | Hydrodynamic analysis of an axial impeller in a non-Newtonian uid through particle image velocimetry[END_REF][START_REF] Mikhaylov | Reconstruction of large-scale ow structures in a stirred tank from limited sensor data[END_REF]. For example[START_REF] Raju | Reynolds number scaling of ow in a stirred tank with Rushton turbine. Part II Eigen decomposition of uctuation[END_REF], implemented the POD algorithm to study the structures associated to the periodic uctuations about the mean ow induced by the impeller in a stirred tank. In this work, experimental PIV data were used as input to the method for dierent diameters of a Rushton-type turbine, and Reynolds numbers between 4,000 and 80,000. Moreover[START_REF] Liné | On POD analysis of PIV measurements applied to mixing in a stirred vessel with a shear thinning uid[END_REF] used experimental PIV data from

a mixing tank to study the presence of coherent structures and the turbulent phenomena associated. In addition, the kinetic energy was successfully obtained from the POD methodology. Likewise

[START_REF] Gabelle | Experimental determination of the shear rate in a stirred tank with a non-newtonian uid: Carbopol[END_REF][START_REF] Gabelle | Kinetic Energy Transfer between First Proper Orthogonal Decomposition Modes in a Mixing Tank[END_REF][START_REF] Fernandes Del Pozo | Hydrodynamic analysis of an axial impeller in a non-Newtonian uid through particle image velocimetry[END_REF]

used the decomposition method in 2D PIV data to reconstruct organized motion induced by impeller blades without performing angle-averaged sampling.

TABLE II .

 II 1: Dimensions (in meters) of the simulated stirred tank.

	C	D	D RZ	H	h z	l B	t B	l p	t p	w
	0.1500	0.1500	0.3000	0.4500	0.0600	0.0450	0.0050	0.0375	0.0020	0.0300

TABLE III .

 III 1: Sinusoidal tting of the mode components, a k (t) associated to the second and third POD modes.

	Modal	A k	ω k	f k	φ k	R square	f (2.5Hz)
	coecient	Amplitude	Angular frequency	Frequency	Phase
		(m/s)	(rad/s)	(Hz)	(rad)
	a 1 (t)	8.320	62.83	10	21.81	1	4
	a 2 (t)	8.310	62.83	10	23.83	1	4

TABLE III .

 III 2: Sinusoidal tting of the mode components, a k (t) associated to the rst eleven POD modes.

	Modal	A k	ω k	f k	φ k	R square	f (2.5Hz)
	coecient	Amplitude	Angular frequency	Frequency	Phase
		(m/s)	(rad/s)	(Hz)	(rad)
	a 1 (t)	454.50	15.71	2.500	17.21	1	1
	a 2 (t)	454.60	15.71	2.500	15.65	1	1
	a 3 (t)	6.510	15.71	2.500	15.64	1	1
	a 4 (t)	6.840	47.12	7.500	7.840	1	3
	a 5 (t)	6.840	47.12	7.500	12.55	1	3
	a 6 (t)	5.150	62.83	10.00	21.91	1	4
	a 7 (t)	5.150	62.83	10.00	26.62	1	4
	a 8 (t)	3.470	78.54	12.50	5.970	1	5
	a 9 (t)	3.470	78.54	12.50	4.390	1	5
	a 10 (t)	1.530	94.25	15.00	21.61	1	6
	a 11 (t)	1.530	94.25	15.00	23.17	1	6

Figure III.13: Lissajous patterns obtained from the temporal variation of the rst to tenth modes.

TABLE III

 III

.3: Normalized Averaged Error as a function of the number of modes used for the reconstruction.

TABLE III

 III

	Procedure	Number of POD modes	Set up time (s)	Reconstruction time (s)	Total time (s)	t P OD Reconst t CF D
	POD Reconstruction	5 10	400	64.5 133.6	464.5 533.6	1.7x10 -2 2.0x10 -2
	CFD		N.A.

.5: Comparative summary of the time duration for the CFD and POD reconstruction procedure. (N.A. = Not Applicable)

TABLE III .

 III 6: Summary of DMD parameters for the top seven modes.

TABLE III

 III

	.7: Comparative summary of the time duration for the CFD and DMD recon-
	struction procedure. (N.A. = Not Applicable)
	Procedure	Number of DMD modes	Set up time (s)	Reconstruction time (s)	Total time (s)	t DM D Reconst t CF D
	DMD Reconstruction	5 11	377	26 46	404 424	1.5x10 -2 1.6x10 -2
	CFD		N.A.

TABLE IV .

 IV 1: Maximum mean error for the POD reconstruction: liquid and gas phases.

	Phase	Number of modes for the reconstruction	Mean Error	Phase	Number of modes for the reconstruction	Mean Error
		600	0.03		600	0.008
	LIQUID	1500	0.02	GAS	1500	0.003
		2996	8×10 -11		2996	8.4×10 -12

TABLE IV .

 IV 2: Maximum mean error for the DMD reconstruction: liquid and gas phases.

	Phase	Number of modes for the reconstruction	Mean Error	Phase	Number of modes for the reconstruction	Mean Error
		600	0.05		600	0.02
	LIQUID	1500	0.02	GAS	1500	0.01
		2996	1×10 -11		2996	3×10 -13

This spectrum is constituted by the smallest scales of the turbulent motion of the uid. In these dimensions of space and time the viscous and inertial eects are of the same magnitude.

We will give more details of these devices in section I.2. For now, our interest is to provide the most relevant characteristics of MI's.

summarizes the information of the dominant DMD modes in the spectrum. All growth rates provided are close to zero. As obtained for the POD case, the rst three modes contribute almost 100 percent of the total energy. The rst and second modes have a frequency of 2.5 Hz. The third mode of ranking has an energy content of 16 percent and a frequency of 0 Hz. The information from the DMD analysis suggests some6 The table shows only one conjugate pair of modes3,61,279,335.

The DMD modes have a real and an imaginary part then the storage space is bigger than the POD modes.

III.2.1.6 Numerical cost of a spatio-temporal reconstruction using the DMD technique

To conclude the exposition of the high-frequency analysis, we present the times involved for the DMD reconstruction of the uid ow. The table III.7 compares the time resource invested for the implementation of each procedure. We need about 377 s for the setup of the information to be used in the dynamical decomposition. This time is slightly less than the corresponding time for the POD. This result is because in the latter the determination of the coecients a j (projection of each CFD snapshot along with each POD mode) consumes an important amount of time.

The reconstruction and implementation of the DMD involve a total investment of 404 s (1.5 × 10 -2 CFD) for 5 frequencies very close to zero, this means the presence of a family of low-frequency modes. We will show more details about this special group shortly. and its "companions" are less and less important. Also, the high frequencies become more and more important for the zones near the turbine.

The conguration of the spectra is similar to that of Figure III.37. The frequencies of 1.5, 2.5, 5, 7.5, 10, and 12.5 Hz continue to show peaks. Additionally, a peak at 3.7

Hz is observed. As with the 1.5 Hz frequency, there is not enough information to assign a physical interpretation to this oscillation.

The ndings of this section allow us to argue that the SVD-DMD procedure is able to identify a uid structure with a frequency of 0.05 Hz (20 s period) that represents the 50 th sub-harmonic of the fundamental frequency of 2.5 Hz. This sub-harmonic is reported in some experimental works [START_REF] Nikiforaki | On the origin, frequency and magnitude of macro-instabilities of the ows in stirred vessels[END_REF][START_REF] Micheletti | Precessional ow macro-instabilities in stirred vessels: study of variations in two locations through conditional phase-averaging and crosscorrelation approaches[END_REF]) and represents a macro-instability whose inclusion in the uid reconstruction would eigenvectors of the gas phase. However the similarity this time is less than the case for the liquid phase. The error results for this method are quite similar to those of the POD. Again, the error associated with 1500 modes is not much smaller than that corresponding to 600 modes.

Moreover, the error committed for both phases is essentially similar.

So far we have shown that the implementation of the DMD for the bubby ow analyzed is feasible. The results show the presence of stable low-frequency modes with an organized spatial conguration for the liquid and gas phases. The ndings also show that the two phases have essentially identical frequencies suggesting coupling for the ows of both phases. We were also able to identify some similarities in the POD eigenvectors and those corresponding to the DMD methodology. These resemblances are to be expected due to the periodicity of the analyzed ow.

Our incursion into the DMD analysis of this system has reached this point. A deeper analysis of the DMD results for this type of system will be devoted to future works.

Remarks on the ROM's results for the liquid and gas phases

The calculation of the errors as presented in this section responds to the need to evaluate kinetic energy of the corresponding phase. A period of about 20 s is associated with the second and third modes. Findings for DMD are similar to those for POD for both phases.

The results obtained from ROMs for the analyzed two-phase ow system are promising, this opens a door for future research work on this type of system.

The next section presents the last part of this document. In such a chapter a novel CMA approach is proposed which includes time evolution. As will be seen, this proposal allows modeling scalar and particle transport phenomena. The proposed dynamic CMA model can be generated from velocity elds that come from CFD or a ROM reconstruction.

Chapter V

Coupling CFD for CMA and particle tracking

V.1 from CFD to CMA

In this part, we briey present the tool developed by [START_REF] Pigou | Modélisation du comportement cinétique, des phénomènes de mélange et de transfert locaux, et des eets d'hétérogénéité de population dans les fermenteurs industriels[END_REF] and Bastien POLIZZI (post-doctoral fellow, IDEX program BIREM) to convert the results coming form a CFD analysis into a matrix of owrates between a set of compartments. The approach is very similar to that developed by [START_REF] Delvigne | Structured mixing model for stirred bioreactors: An extension to the stochastic approach[END_REF] and later adapted to compartment models generated from CFD data by [START_REF] Delafosse | Development of a compartment model based on CFD simulations for description of mixing in bioreactors[END_REF][START_REF] Delafosse | CFD-based compartment model for description of mixing in bioreactors[END_REF].

A compartment is a volume of uid, much larger than the computational cell size used in the CFD. The collection of compartments occupies the same volume as the total volume simulated in the CFD. The compartments are numbered and the owrate leaving or entering a given compartment is extracted from the integration of the velocity eld around its contour. As depicted in gure V.1, the ow through a given face between two adjacent compartments is characterized by a pair of values which dene the entering and outing owrates with respect to the said compartments. This tool can handle single and two phase ows in cylindrical geometries. Besides the owrates between compartments, an integration of some quantities is performed so as to provide the mean value of the gas hold-up or the kinetic energy dissipation within each compartment.

From these values, a square matrix N c × N c is built as follows: the owrate from a compartment i to a compartment j is placed on the i th line and the j th column of the said matrix. N c refers to the number of compartments.

Let n p (i, t) be the number of particles found in compartment i at time t.

We further assume that there is no correlation between the particle location and this allows to express the mean number of particles that one expects to nd in a compartment as

By denition of the particle density, the following property holds in well mixed reactor

which leads to the expectation of X i,n (t)

The variance of the number of particle in a compartment expresses as

which leads after calculations to

Thus the relative deviation to the expected number of particle in a given compartment writes

Concerning the bubble plume in a quasi 2D column, a short term perspective would be to further analyse the dierences and similarities in the energy spectra of the two phases.

The spectra mode are markedly dierent in the two phases for the low order modes. This might contains some information regarding the energy transfer between phases. Since the input data are well resolved spatially and temporally, one can expect to nd the evidence of physical phenomena in the entire range of modes. It is expected for example that the slopes of the spectra is related to the turbulence and possibly helps to identify a relationship between gas and liquid turbulence.