Fully Resolved Dynamics of Mixing in Confined Impinging Jets Reactors
Full text
Fully Resolved Dynamics of Mixing in Confined Impinging Jets Reactors Master Dissertation in Computational Mechanics by Hélder Manuel Marques Salvador Supervisors Ricardo Jorge Nogueira dos Santos Co-Supervisor José Carlos Brito Lopes Laboratory of Separation and Reaction Engineering - LSRE Laboratório Associado LSRE/LCM Departamento de Engenharia Química Faculdade de Engenharia Universidade do Porto September 2015
Acknowledgments I wish to thank my supervisor, Professor Ricardo Santos, for his guidance, helpful technical discussions, continual availability during the development of my work and for encouraging me to carry out this master thesis at LSRE - Laboratory of Separation and Reaction Engineering, in collaboration with the Computational Mechanics Department. I am indebted to my co-supervisor, Professor José Carlos Lopes, for his expert support, for his trust on my work and his enthusiasm for this project. I also acknowledge Professor Madalena Dias, director of the LSRE - Laboratory of Separation and Reaction Engineering, for her availability and for providing me all the necessary facilities to perform my work. I would also like to express my appreciation for the knowledge all the professors have transmitted throughout these two years. I would like to thank specially to Professor José Sá for his availability and for his efforts in opening the Computational Mechanics master course. During my time at FEUP the help and support of my friends and the collaboration with my classmates have been indispensable, for which I thank them dearly. In particular, I owe a special thanks to my colleagues Mohamed Lotfi, Behzad Farahani and José Berardo. My gratitude also goes to the colleagues in the LSRE, particularly the Mixing Group, which made this work possible through their help and support, namely Carlos Teixeira, Gabriela Ruphuy, Joana Pereira, Marcelo Costa, Mónica Silva, Nelson Gonçalves and Rómulo Oliveira. I also have to acknowledge my brothers in life: Ivo Cardoso, João Faustino, Luís Branco, Marco Pereira and Paulo Correia. I will see you next summer. Finally, I wish to convey my deepest gratitude to my family, particularly my parents, for their unconditional support.
To my parents, family and friends
“Big whorls have little whorls that feed on their velocity, and little whorls have lesser whorls and so on to viscosity” Lewis F. Richardson “When the present determines the future, but the approximate present does not approximately determine the future” Edward Norton Lorenz
Abstract Confined Impingement Jets mixers (CIJ) have attracted wide interest in the past years because of the efficient micromixing and possible application to nanoparticles production. The Reaction Injection Molding (RIM) technology employs CIJ due to the fast reaction between the two monomers used to create plastic components. In the present work, the flow scales were calculated for Reynolds numbers 150, 200, 250 and 300 until the smallest hydrodynamic scale, the Kolmogorov scale. From the numerical analyses of the flow field, the oscillatory behaviour of CIJ was studied and the energy containing scales identified. Additionally, the turbulent kinetic energy and the turbulent energy dissipation were calculated for all the numerical studies. A multiphase VOF model is used for the study of 3D mixing inside this reactor with an edge length equal to the Kolmogorov scale, in order to evaluate the mixing related mechanisms and the adequacy of the numerical grid for a non-diffusive model. Particle Image Velocimetry (PIV) measurements were carried out in order to compare the flow field behaviour of the front plane with the numerical analyses performed. Additionally, two methods previously developed for the CIJ jets impingement point, the elastic analogue model and the jets pressure model, were tested for a wide range of Reynolds numbers in order to assess their accuracy in predicting the impingement point position under different flow conditions. The flow field was also visualized with the Planar Laser Induced Fluorescence (PLIF) experimental technique using a high speed camera to capture all the large flow scales. An image processing algorithm was developed to compute the jets impingement point position and the jets angle. The PLIF technique was performed for several Reynolds numbers and the results obtained from the image processing and by the differential pressure transducer were processed using Fast Fourier Transforms (FFT). The frequency of the jets oscillations was then compared with results reported in the literature. In addition, the PLIF technique was used to study the impingement mixing of fluids with different rheological properties. The experimental measurements were subsequently compared with the results provided by the two RIM control models for fluids with different viscosities in each injector, in order to validate them for a wide range of jets’ momentum ratios.
Resumo A mistura por jatos opostos (CIJ) tem atraído grande interesse nos últimos anos devido à eficiente micromistura e à potencial aplicação na produção de nanopartículas. A máquina de moldagem por injeção com reação química (RIM), a qual é utilizada para produzir componentes de plástico, usa a mistura por jatos opostos devido à reação rápida que ocorre entre os dois monómeros. No presente trabalho, as escalas do escoamento foram calculadas para números de Reynolds de 150, 200, 250 e 300, até à menor escala hidrodinâmica, a escala de Kolmogorov. A partir da análise numérica do campo de velocidades, foi estudado o comportamento oscilatório dos jatos opostos e a energia das escalas identificadas. Além disso, a energia cinética turbulenta e a dissipação de energia turbulenta foram determinadas para todos os estudos numéricos. O modelo multifásico VOF é usado para o estudo 3D do interior deste reator com uma discretização espacial igual à escala de Kolmogorov, a fim de identificar os mecanismos de mistura e de avaliar a escala de discretização espacial acoplada com um modelo nãodifusivo. A técnica de Velocimetria por Imagem de Partículas (PIV) foi realizada com o intuito de comparar o comportamento do campo de velocidades no plano dos eixos da câmara de mistura e injetores com as análises numéricas realizadas. Além disso, dois métodos desenvolvidos previamente para o controlo do processo RIM, o modelo elástico e o modelo de pressão, foram testados para vários números de Reynolds no sentido de avaliar a sua adequação para a previsão da posição do ponto de impacto dos jatos em condições de escoamento diversas. O campo de velocidades foi igualmente visualizado por via da técnica experimental de Fluorescência Induzida por Laser (PLIF). Com recurso a uma câmara de alta velocidade, que permite capturar instantaneamente imagens da câmara de mistura, foram obtidas as grandes escalas do escoamento. Foi ainda construído um algoritmo para avaliar a posição do ponto de estagnação e o ângulo formado pelos jatos. As medições experimentais com PLIF foram realizadas para uma vasta gama de números de Reynolds e os resultados obtidos quer pelo algoritmo, quer pelo transdutor diferencial de pressão, foram processados por transformadas rápidas de Fourier (FFT). A frequência de oscilação de cada variável foi em seguida comparada com resultados reportados na literatura. Finalmente, a técnica PLIF foi aplicada no estudo da injeção de fluidos com diferentes viscosidades. Os resultados experimentais obtidos foram comparados com os resultados dos modelos elástico e de pressão, de modo a verificar-se a sua validade nestas condições de operação.
Contents 1 Introduction ................................................................................. 1 1.1 Relevance and Motivation ............................................................ 1 1.2 Mixing Process, an Introduction for the Case Study .............................. 2 1.3 Reaction Injection Moulding in Confined Impingement Jets .................... 3 1.4 Previous work .......................................................................... 7 1.5 Thesis objective and layout .......................................................... 9 2 State of art of RIM and CIJ………………………………………………………………………………..11 2.1 Introduction ........................................................................... 11 2.2 Mixing Study ........................................................................... 12 2.2.1 Characterization of the mixtures ............................................. 12 2.2.2 Residence time distribution RTD .............................................. 13 2.2.3 Measurement of the degree of mixing ....................................... 14 2.2.4 Mixing mechanism ............................................................... 17 2.2.5 Mixing and hydrodynamics simulation........................................ 20 2.3 Mixing in Confined Impinging Jets ................................................. 28 2.3.1 Hydrodynamics in CIJ - Experimental Studies ............................... 28 2.3.2 Hydrodynamics in CIJ - Numerical Studies .................................. 36 2.3.3 Mixing and Chemical Reaction in CIJ ......................................... 39 3 CFD Model and Simulation of CIJ ..................................................... 45
Table of Contents 3 𝑢𝑖𝑛𝑗 Average velocity at the injectors [m/s] V Volume [𝑚3] 𝑉𝑑𝑖𝑠𝑠𝑖𝑝𝑎𝑡𝑖𝑜𝑛 Volume where dynamic impingement mixing occurs [𝑚3] W Dimensionless weight x Position vector [m] X Initial position vector [m] 𝑥𝑖 Component of x in i-direction X𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒 Mass fraction of glycerine in the glycerol solution 𝑥𝐼𝑃 Impingement point position [m] Greek letters 𝛼𝑖 Volume fraction of phase i 𝛼 Spatial volume fraction average ∆𝑡 Integration time step [s] ∆𝑡𝑏𝑒𝑎𝑚𝑠 Interval between consecutive laser flashes in the PIV experiments [s] ∆𝑥 Numerical grid element length [m] ∆𝑥𝑖𝑛𝑡,𝑟𝑒𝑔𝑖𝑜𝑛 Length of the interrogation region in PIV experiments [m] ∆𝑝 Differential pressure [Pa] ∆𝐸𝑝 Storage potential energy [J] 𝛾 Strain rate [𝑠−1] 𝜀 Turbulent energy dissipation rate [𝑚2/𝑠3] 𝜙 Frequency [𝐻𝑧] Κ Turbulent kinetic energy [J/kg] 𝜆 Wave length [𝑚] 𝜆𝐵 Batchelor Scale[𝑚] 𝜆𝑘 Kolmogorov Scale [𝑚] 𝜉 Impingement point position [𝑚] 𝜌 Fluid density[𝑘𝑔/𝑚3] 𝜌𝑝 Particles density[𝑘𝑔/𝑚3] 𝜌 Average fluid density weighted by the flow rate [𝑘𝑔/𝑚3] 𝜎 Standard deviation 𝜎2 Variance 𝜎𝛼 2 Spatial variance of the volume fraction 𝜇 Viscosity [𝑃𝑎. 𝑠]
Table of Contents 4 𝜇 Average fluid viscosity weighted by the flow rate [𝑃𝑎. 𝑠] 𝜏 Passage time in the mixing chamber [𝑠] 𝜏𝑖𝑖 Normal Stress [𝑃𝑎] 𝜏𝑖𝑗 Shear stress tensor [𝑃𝑎] 𝜙𝐾 Kinetic energy rate ratio 𝜙𝑀 momentum rate ratio 𝜙𝐹𝑅 Flow rate ratio Ω Volume or area domain ∂Ω Domain boundary surface or perimeter 𝜔 Non-dimensional vorticity 𝜒∗ Normalized impingement point position Mathematical operations ⟨ ⟩ Spatial Average 𝜕𝑖(∙) Differential operator 𝐞 𝒊 Normalized vector x Cross Product ∙ Dyadic product 𝛻(∙) Gradient operator 𝛻 𝑥(∙) Differentiation with respect of 𝑥 𝛻2(∙) Laplacian operator (∙)𝑇 Transpose (∙) Average Indices ∗ Dimensionless variable ‘ Fluctuating term 2𝐷 For 2D of the chamber 3𝐷 For 3D of the chamber A Chemical specie A B Chemical specie B C Chemical specie C i Index of the counter inj At injector
Table of Contents 5 j Index of the counter k Index of the counter max Maximum of a variable min Minimum of a variable R Chemical specie R S Chemical specie S x x-component of the variable y y-component of the variable z z-component of the variable 1 Left injector 2 Right injector 0..10 Probe points Abbreviations 2D Two-dimensional 3D Three-dimensional CFD Computational Fluid Dynamics CIJ Confined Impinging Jets CPU Central Processing Unit DIP MODA Digital Image Processing Multiobjective Optimization Deterministic Algorithm GPU Graphics Processing Unit DFT Direct Fourier Transform DNS Direct Numerical Simulation FFT Fast Fourier Transform LDA Laser Doppler Anemometry LES Large Eddy Simulation LIF Laser Induced Fluorescence NAJ Narrow Axisymmetric Jets NS Navier-Stokes PDF Probability Distribution Function PIV Particle Image Velocimetry PLIF Planar Laser Induced Fluorescence RAM Random Access Memory RANS Reynolds Averaged Navier-Stokes RIM Reaction Injection Moulding
Table of Contents 6 RTD Residence Time Distribution TKE Turbulent Kinetic Energy VOF Volume-Of-Fluid GUI Grafical User Interface
List of Figures Figure 1.1 - Scheme of a RIM machine adapted from Macosko (1989) .................... 5 Figure 1.2 – RIM machine commercialized by KraussMaffei ................................ 5 Figure 3.1 - Schematic representation of the CIJ mixing chamber and related physical dimensions. Retrieved from Fonte (2012) ................................................... 47 Figure 3.2 – Hexahedral mesh with used in the simulations (left). Detail of the triangular cells used to couple both meshes avoiding distorted elements (right) ..... 49 Figure 3.3 - Perturbation time function on the right injector affecting the normal velocity to the injector's cross section ....................................................... 55 Figure 3.4 - Velocity maps for Reynolds 150 at time / = ............................. 60 Figure 3.5 - Velocity maps for Reynolds 150 at time / = ............................. 60 Figure 3.6 - Velocity maps for Reynolds 150 at time / = ............................. 60 Figure 3.7 - Velocity maps for Reynolds 200 at time / = ............................. 61 Figure 3.8 - Velocity maps for Reynolds 200 at time t/τ=2 ............................... 61 Figure 3.9 - Velocity maps for Reynolds 200 at time / = ............................. 61 Figure 3.10 - Velocity maps for Reynolds 250 at time / = ........................... 62 Figure 3.11 - Velocity maps for Reynolds 250 at time / = ........................... 62 Figure 3.12 - Velocity maps for Reynolds 250 at time / = ........................... 62 Figure 3.13 - Velocity maps for Reynolds 300 at time / = ........................... 63 Figure 3.14 - Velocity maps for Reynolds 300 at time / = ........................... 63 Figure 3.15 - Velocity maps for Reynolds 300 at time / = ........................... 63 Figure 3.16 - Vorticity map for Reynolds 150 at time / = ............................ 66
Figure 3.17 - Vorticity map for Reynolds 150 at time / = ............................ 66 Figure 3.18 - Vorticity map for Reynolds 150 at time / = ............................ 66 Figure 3.19 - Vorticity map for Reynolds 200 at time / = ............................ 67 Figure 3.20 - Vorticity map for Reynolds 200 at time / = ............................ 67 Figure 3.21 - Vorticity map for Reynolds 200 at time / = ............................ 67 Figure 3.22 - Vorticity map for Reynolds 250 at time / = ............................ 68 Figure 3.23 – Vorticity map for Reynolds 250 at time / = ............................ 68 Figure 3.24 - Vorticity map for Reynolds 250 at time / = ............................ 68 Figure 3.25 - Vorticity map for Reynolds 300 at time / = ............................ 69 Figure 3.26 - Vorticity map for Reynolds 300 at time / = ............................ 69 Figure 3.27 - Vorticity map for Reynolds 300 at time / = ............................ 69 Figure 3.28 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 150 ................................................................................. 72 Figure 3.29 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 200 ................................................................................. 72 Figure 3.30 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 250 ................................................................................. 73 Figure 3.31 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 300 ................................................................................. 73 Figure 3.32 - Schematic illustration of the energy cascade. It shows the amount of turbulent energy E associated with different wave number k ............................ 75 Figure 3.33 - Turbulent kinetic energy along the chamber axis .......................... 76 Figure 3.34 - Cascade of transport and dissipation of turbulent energy ................ 77 Figure 3.35 - Lines aligned with origin of the chamber .................................... 78
Figure 3.36 - Impingement point position, X-velocity along the chamber axis ......... 79 Figure 3.37 - Impingement point position, Y-velocity along the injector axis ......... 79 Figure 3.38 - Impingement point position, Z-velocity along the injector axis ......... 80 Figure 3.39 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 150 ..................................................................................... 83 Figure 3.40 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 200 ..................................................................................... 84 Figure 3.41 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 250 ..................................................................................... 85 Figure 3.42 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 300 ..................................................................................... 86 Figure 3.43 - Algorithm to obtain fluid flow variables and ........................... 87 Figure 3.44 -Convergence due averaging time steps for and ......................... 88 Figure 3.45 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 150 ........................................................................ 91 Figure 3.46 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 150 ........................................................................ 92 Figure 3.47 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 200 ........................................................................ 93 Figure 3.48 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 200 ........................................................................ 94
Figure 3.49 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 250 ........................................................................ 95 Figure 3.50 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 250 ........................................................................ 96 Figure 3.51 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 300 ........................................................................ 97 Figure 3.52 - Variation of along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of for the top part of the chamber (bottom). Reynolds 300 ........................................................................ 98 Figure 3.53 - Variation of and along the numerical experiments. Observation energy containing inside the impingement structure (left). Normalization of this energy with energy containing in the central area (right) ............................................. 100 Figure 3.54Power spectra form all the Reynolds number performed [150, 200, 250, 300] for the probe ........................................................................ 104 Figure 3.55Power spectra form all the Reynolds number performed [150, 200, 250, 300] for the probe ........................................................................ 105 Figure 3.56 - Phase maps for Reynolds 150 at time / = ............................ 112 Figure 3.57 - Phase maps for Reynolds 150 at time / = ............................ 112 Figure 3.58 - Phase maps for Reynolds 150 at time / = ............................ 112 Figure 4.1 – Schematic of the transparent mixing chamber and the mold ............ 117 Figure 4.2 – Representation of the RIM equipment: a) hydraulic loop, b) transparent mixing chamber and mold ................................................................... 118 Figure 4.3 – Variation of the glycerine-water solution with temperature ............ 119 Figure 4.4 – PIV experimental setup ........................................................ 121
Figure 4.5 – Schematic of the PIV experimental setup .................................. 122 Figure 4.6 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 45 ................................................................................ 127 Figure 4.7 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 45 ................................................................................ 128 Figure 4.8 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 45 ................................................................................ 129 Figure 4.9 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratio = . at Reynolds 45 ......... 130 Figure 4.10 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 111............................................................................... 131 Figure 4.11 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 111............................................................................... 132 Figure 4.12 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 111............................................................................... 133 Figure 4.13 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . ] at Reynolds 111 ................................................................................................... 134 Figure 4.14 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 159............................................................................... 135 Figure 4.15 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 159............................................................................... 136
Figure 4.16 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 159............................................................................... 137 Figure 4.17 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . ] at Reynolds 159 .................................................................................. 138 Figure 4.18 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 220............................................................................... 139 Figure 4.19 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 220............................................................................... 140 Figure 4.20 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 220............................................................................... 141 Figure 4.21 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 220............................................................................... 142 Figure 4.22 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 220............................................................................... 143 Figure 4.23 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 313............................................................................... 144 Figure 4.24 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 313............................................................................... 145 Figure 4.25 – Velocity field vector maps and non-dimensional velocity fluctuations contours for jets’ kinetic energy rate ratios = [. , . , . , . ] at Reynolds 313............................................................................... 146
List of Tables Table 1.1 - Annual revenue in billions of euros of the three major chemical industries ........................................................................................ ..7 Table 3.1 - Computational cost of each simulation performed in HPC environment .. 57 Table 3.2 - Computational cost of the VOF methodology comparing with previous simulation for the same Reynolds number ................................................ 110
1 1 Introduction 1.1 Relevance and Motivation In the second half of the 20th century, plastic became a most universally-used and multipurpose material in the global economy. From the beginning of production in 50’s it was boosted with an increasing rate of almost 9% per year from 1950 to 2012 (PlasticsEurope, 2015). Nowadays plastic is globally used in components which require low weight, excellent finishing and low cost. Depending on the conditions necessary to produce a plastic component, two major groups emerge to classify those plastics, thermoplastic and thermoset. Thermosoftening (thermoplastic) plastics is a material which requires a specific temperature to be mouldable, weakening the intermolecular bound forces, fulfilling the mould and obtaining their shape upon cooled. Thermosettting plastics is a petrochemical material which requires a chemical reaction to create permanent intermolecular bounds fusing inside the mould and therefore obtaining their shape. A wide range of industrial processes to produce plastic components exist in the industry in order to adapt for the materials used due to chemical transformation, number of parts to be produced, size and parameters which constrains the mould and process. One of this process is Reaction Injection Moulding (RIM) that uses confined impingement jets mixers
1 Introduction 2 (CIJ), which are the topic of the present thesis. This method have a set of advantages being used particularly to produce component for the aeronautic and automotive industries. In Section 1.2 will be introduced the process of mixing for the particular case study, followed by a description of the RIM process in the Section 1.3 where will be introduced the process of producing plastic components in CIJ reactors. In the end, in the Section 1.5 will be presented an overview of the thesis layout and objectives. 1.2 Mixing Process, an Introduction for the Case Study Mixing is a critical stage in any chemical industry and it is used to combine two or more material (fluids or solids) into a homogenous solution. A detailed study of the mixing process is necessary due to limitations such as space, material and operation cost which strongly contribute for the performance of the industrial process. In the early ages, mixing was performed using a continuous flow stirred-tank reactor (CSTR), based on a tank containing the solution to be homogenized and an impeller to perform the mixing. This methodology is still used nowadays for some specific process, but due to energy costs and poor distribution of residence time, it starts to be converted into other mixing techniques more environmentally friendly, with lower production cost and better control of the mixing process. The work presented in this thesis introduces the study of the flow in Confined Impingement Jets (CJI) where two fluids are injected in a confined chamber though opposed jets. In the impingement zone where the two jets contact a rapid deceleration occurs, dissipating energy radially across the chamber. In this methodology two different geometries emerge: T-Jets with a rectangular section for the injectors and chamber; CIJ with a cylindrical cross section of the injectors and chamber. Although both geometries used the same technique to mix the fluids, different flow structures emerge along the flow regimes studied in the literature. This geometry is quite simple to manufacture, with no moving components, but a very low roughness in the chamber and injector walls is required. This mixing technology is very advantageous for rapid chemical processes where small micromixing scales have to be reached in short times, for example in Reaction Injection Moulding (RIM), and in applications such as the production of nanoparticles via wet chemistry routes.
1 Introduction 3 1.3 Reaction Injection Moulding in Confined Impingement Jets The Reaction Injection Moulding process (RIM) was initially developed in Bayer AG in 1964 in order to produce plastic components initially from polyurethanes, but due to the increasing demand of plastic production was adopted to produce polyureas, polyamides, polyesters and polyepoxides (Fonte, 2012). This method is suited to produce large plastic parts and this is why nowadays is used in the aeronautical and automotive industries. In this process two fluid streams enter in the mixing chamber as opposed jets and mixing occurs rapidly, meaning low residence times order of 1/10 s. In an industrial machine fluids exit the injectors with a velocity of 10 to 150 m/s but due to high viscosity rate of the polymers, up to 1 Pa.s, the Reynolds number are low from 100 to 500. Due to this factor the process demands powerful pumps to produce pressures in order of 100 to 200 bar and therefore equipment to support such pressures. The industrial mixing chambers range from 3 to 15 mm of diameter, with a height five times the diameter and injectors diameter in range 1 to 3 mm (Oertel, 1985; Santos, 2003). Although such pressures are required in the fluids injection, the mixture enters in the mould at low pressure, requiring low cost moulds with low weight (Teixeira, 2000). Adding to this factor the large flow rates, it means the small machines can produce large plastic parts. Inside the mixing chamber the instability of the jets creates mechanisms of mixing due to the fast deceleration of the jets towards the impingement point. Although high velocities are produced, the fluid motion is dominated by viscous forces relatively to inertia, meaning no turbulent flow is observed. Due to the low residence times, there is no meaningful polymerization of the fluids inside the mixing chamber, and also due to the high viscosity no considerable molecular diffusion occurs. The mixture mechanisms consist of the stretching of fluid from the jets, followed by bending and folding, which creates a lamellar structure. It is the thickness of these lamellar structures that contributes to mixing and it is considered that both fluids are mixed when their laminas reach a thickness that is sufficient to allow for a full polymerization to occur inside the mould. The mixture is the most important factor of this machine and due to this factor a precise control of the flow rate is needed to deliver the fluids in stoichiometric proportions (Santos, 2003) and perfectly mixed. This process is highly adaptive since it can also include a vast number of reinforcements, fibbers and reagents in order to add stiffness to the component. Each combination of reinforcement fibbers and reagents species introduced in the mould changes the
1 Introduction 4 mechanical response to applied forces, and therefore the part produced can be optimized regarding its weight and strength. Due to this feature the author believes in the future this process will be more widely used, due to the demand in the industry, particularly in automotive and aeronautical, by trustful, low weight parts with good mechanical proprieties that enable substituting expensive materials such as aluminium, steel and others. Figure 1.1 shows a graphical representation of a RIM machine. In this scheme two major circuits are observable, an injection and a recirculation circuit. A circuit represented by the storage tanks (7), heat exchanger (5) and recirculation pumps (4) recirculate the fluid for the thermal stabilization of the fluids. The injection circuit is presented next 1. Monomers and polymers, in this example Polyol and Isocyanate are storaged in separated storage tanks (7). 2. High pressure cylinders (6) are used to deliver fluid into the mixing chamber (1). In this process a good control of the flow rates guaranties control of the process and deliver the fluids with accurate stoichiometric ratios 3. The reactants are injected into the mixing chamber where they are mixed until the mould (2) is filed. 4. In the end of the process a piston is shot into the mixing chamber to remove any remaining material. 5. Curing of the polymer in the mould, which depending on the mould size and formulation used, could be 1 to 30s. The time of this process is the limiting step of the entire RIM process (Santos, 2003) 6. Demoulding Figure 1.2 shows a RIM machine commercialized by KraussMaffei.
1 Introduction 5 Figure 1.1 - Scheme of a RIM machine adapted from Macosko (1989) Figure 1.2 – RIM machine commercialized by KraussMaffei
1 Introduction 6 According the advantages and disadvantages, it is possible to say Advantages o Suitable for low quantities. If the cost of production overweights the cost of investment this methodology starts to become inefficient. o Low tooling investment, due to low pressure injection it requires lower mould strength and as a consequence this process can be used for prototyping. o The RIM process is not dependent of the parts dimensions or weight, being capable to handle complex geometries with thin walls, producing components with excellent finishing allowing high quality painting. o By adopting different reaction species it is possible to choose specific mechanical proprieties Disadvantages o RIM requires excellent control of the environmental conditions and operational variables due to high sensitivity of the process o RIM require an excellent know-how to modify/adjust the process variables to different reactors in order to obtain a proper mixture o Knowledge in mould designing to archive fast and proper demolding with a good finishing quality Financial data of the materials for RIM is shown in Table 1, from the annual revenue of the top chemical industries: Bayer AG, Huntsman, BASF (BASF, 2014; Bayer, 2014; Huntsman, 2014). It is separated the global revenue of the companies versus the revenue for polyurethanes, principal RIM material, respectively for the years 2014 and 2013. It is conclude the market for polyurethanes represents an important share of the companies revenue. An average growing rate of polyurethanes market of 4% demonstrates the economic importance of this market in the production of plastic components. The report of Huntsman it was converted from Dollars to Euros by an average exchange in the two years of this study of 0.79 €/$.
1 Introduction 7 Table 1 - Annual revenue in billions of euros of the three major chemical industries Global revenue polyurethanes Growing rate of polyurethanes (%) 2014 2013 2014 2013 BASF 74.326 73.973 6.135 5.708 6.96 Huntsaman 9.15 8.75 3.76 3.92 1.35 Bayer AG 42.239 40.157 6.285 6.054 3.67 The increasing demand for plastic components allied with the introduction of new capabilities in RIM technology and the poor knowledge of the full resolved fluid flow motion inside the mixing chamber, in order to control the process, constitute the driving forces for the present thesis that pursues a research line that was initiated by Teixeira (2000), and continued with Santos (2003), Fonte (2012) and Gomes (2015). 1.4 Previous work This work is the follow up of the studies performed by Teixeira (2000), Santos (2003), Fonte (2012) and Gomes(2015) in their PhD thesis. Therefore, before presenting the thesis objective and layout is important to situate the starting point. Teixeira (2000) performed an experimental and computational studies of the dynamics of the flow field using for the experimental side Laser Doppler Anemometry (LDA) comparing with a 2D model representative of the mixing chamber simulated using the CFD code FIDAP. In the experimental work, Teixeira (2000) measured the two velocity components normal to the chamber axis for a Reynolds range from 50 to 600. From the LDA experiments the flow variables were computed, namely the Reynolds stresses, turbulence intensity and the distribution of flow field for several planes inside the chamber. From these studies was related the flow field inside the chamber with the Reynolds number in the injector showing the existence of a critical Reynolds for where exist a full segregation of the jets inside the chamber. Also it was observed for the first time the existence of characteristic frequencies of the jets oscillation.
1 Introduction 8 Teixeira (2000) also performed Direct Numerical Simulation (DNS) of the flow field using a 2D computational geometry. This simulation was performed for the same Reynolds range as the experiments, being thoroughly studied the numerical parameters that affect the numerical results. It was proven by comparing with other published works the effectiveness of the numerical simulation to predict the actual behaviour of the flow field. Santos (2003) performed experiments with LDA for Reynolds numbers ranging from 250 to 600 in which was studied the dynamic behaviour of the flow field and the oscillation frequencies. Using a 2D geometry the flow was modelled using the commercial CFD code FIDAP to correctly replicate the dynamic behaviour of the flow field providing answers for the jets oscillation mechanisms. This was also done using experimental PIV technique. Santos (2003) work is based in experimental work with LDA and PIV to study the oscillation frequencies and the oscillations mechanism of the flow field. 2D numerical simulation were performed to accurately reproduce the flow field dynamics and the mass transfer mechanism. LDA was used to study the flow for ranges of Re from 260 to 600, the oscillation frequencies were studied for velocity components with the same direction of the injection in order to compare the results with the 2D numerical simulations, a study made near the impingement region to access the nature of the flow structures is performed. From PIV experimental work the flow field maps were created for Reynolds ranged from 100 to 500 helping to clarify and validate the hydrodynamics of the flow field when compared with numerical data. From experimental work operation parameters such as Reynolds and Froude number and jet’s momentum ratio as study helping to provide data for better control of the process. Using 2D CFD simulations with a commercial software Fluent, Santos (2003) studied the mass transfer mechanisms in the RIM process, where the numerical and physical mechanisms was presented. Also a chemical reaction was perform to better understanding the role played by the flow field in the RIM process. Fonte (2012) based on 3D numerical simulation made a Lagrangian particle tracking to study the distribution of striation thickness along the chamber using the commercial software Fluent. Also, using numerical simulation it was performed a study using a multiphase model, Volume of Fluid (VOF), coupled with an adaptive mesh refinement based on the gradient of the chemical concentrations for a 2D representation of the mixer chamber with the objective to model the smaller mixing length scales, proven the
1 Introduction 9 applicability of this technique to flow visualization although with extreme computational cost. In the experimental side PLIF technique was used to qualify the mixing patters of the entire flow field. PLIF experiment results, validate an elastic model, mathematically deduced by Fonte (2012), to predict the position of the impingement point for a different set of flow conditions. It proves the validity of the model handling different flow rates and for different circular injector’s cross sections. 1.5 Thesis objective and layout The present thesis is divided into two main objectives. The first objective is to study the fluid flow and mixing mechanisms in the RIM mixing chamber, using more than one technique. The second objective is to study the jets interaction under balanced and unbalanced flow conditions to evaluate more than one mathematical equation used to developed a robust control for RIM processes. This thesis is divided into six chapters. The first chapter corresponds to the introduction of the theme and the relevant motivation to study it. The second chapter presents the state of the art of mixing and flow dynamics applied to RIM. The rest of the chapters are divided into different techniques for flow analyses. The main objectives are to investigate and to contribute to widen the knowledge of the governing mixing mechanisms and process control of the RIM machines. In Chapter 3 DNS simulation is performed for Reynolds 150, 200, 250 and 300 in order to study the energy containing flow scales where the energy cascades are studied to compare the vortices length in relation to the chamber dimensions. Additionally turbulent quantities such as the turbulent kinetic energy and the turbulent energy dissipation rate are presented for all the Reynolds numbers studied. VOF model is implemented for Reynolds 150 in order to study the advection mechanisms without the effect of diffusion and not considering surface tension. In Chapter 4, PIV measurements in a plane passing thought the injector axis are presented in order to perform the study of balanced and unbalanced jets using two mathematical models to predict the impingement point position. The results of the mathematical models to control the mixing process will be compared with the experimental results
2 State of art of RIM and CIJ 16 Where A, B, R and S are chemical species and 𝑘1, 𝑘2 and 𝑘3 are chemical reaction kinetic constants. Consider the 1st order chemical reaction: the reactant A will then create R after the time 𝑘1 it is complete. The kinetic constant 𝑘1 access the time required for all A is transform to R, and if this time it is not archived the reaction is incomplete. In a case of a 2nd order reaction two different reactions are parallel created and in case of 𝑘1 is fast enough and B are completely converted to R so, no subproduct S will be created. The full conversion of product A into R for the 1st order reaction or the non-existence the reactant S in the 2nd order chemical reaction is used to measure the micromixing efficiency. The main conclusions presented by Danckwerts (1952, 1958) and by Weinstein and Adler (1967) show that 1sr order reactions are almost full dependent of the macromixing rather than micromixing wile 2nd order reaction are greatly affected by both micro and macromixing. In a case of a 1st order reaction segregation can cause the increase of the kinetic reaction constant which can be overpass by increasing the reactor size, but, in case of 2nd reactions it cannot be possible. This reaction is nowadays extremely important to access the distribution of consecutive reactions and due their sensibility to both micro and macromixing to evaluate the mixing performance. Kusch et al. (1989) and Mahajan and Kirwan (1996) used 2nd order reaction to study a wide range of reactors particularly confined impingement jets reactors. Triplet consecutive reactions (Bourne & Yu, 1994; El-Hamouz & Mann, 1998) allow wide flexibility to characterize the mixing performance, when compared with 2nd and 1st order consecutive reactions. 𝐴+𝐵 𝑘1 𝑅 (2.10) 𝑅+𝐵 𝑘2 𝑆 (2.11) 𝑆+𝐵 𝑘3 𝑇 (2.12) In this case, perfect mixing is obtained by full conversion of A and B into R. Complete segregation favour creation of T while the creation of S is due complete segregation of A and B from the initial conditions.
2 State of art of RIM and CIJ 17 Parallel competitive reaction with different constant kinetic reaction coefficient one fast and other slower is study by Bourne, Kut, Lenzner, and Maire (1990) and J. Villermaux, Falk, Fournier, and Detrez (1992) 𝐴+𝐵 𝑘1 𝑅 (2.13) 𝐴+𝐶 𝑘2 𝑆 (2.14) In the fast reaction is occur if a past homogenization of the reactant is present, in the other and, the slower reaction helps to compare the mixing efficiency. 2.2.4 Mixing mechanism By Danckwerts (1952) mixing process can be divided in two main space scales, micromixing and macromixing which is directly related to the scale where the mixing occur. Micromixing is processed in to microscopic level and macromixing is due the effect of large scale eddies in the fluid flow system. In any fluid flow system exist two mixing mechanisms identified by Bourne (1997) relating how they occur. If two reactants are mixing in a molecular level the mixture is called diffusive mixing in opposite convective mixing occur due the effect of convection in spreading the fluid along the reactor. In the smaller scales, microscale, the fluids are homogenized by diffusive mixing and convective mixing at the lower scales. In bigger scales, macroscales, the fluid is homogenized only due convective mixing but remains full heterogeneous in the lower scales Although mixing mechanism is not fully understood, Beek (1959) and Jacques Villermaux and David (1983) proposed as a ground foundation three stages in where the mixing is created. This three stages is them directly related for this particular process, RIM. Fluid is grossly distributed along the reactor with heterogynous mixing in the macro scale: In a RIM two fluids are injected by opposed jets where a instability region created by the large scale eddies spread the fluid material along the chamber. Reduction of the mixing scales: After big eddies distribute the fluids, smaller hydrodynamic scales will homogenized the fluid in a macro and micro scales along is going thought the reactor. In this stage the surface area will exponential increase by stretching increasing and thickness reduction.
2 State of art of RIM and CIJ 18 Mixture is produced by diffusion: Due the low residence times, molecular diffusion is not expected to occur significantly inside the reactor. Molecular diffusion will only play a significant role by reaction inside the mould. 2.2.4.1 Lamellar model In a laminar systems the striation thickness of the fluid is increased along the chamber due viscous forces creating inside the chamber, by bending and folding, lamellar structures of fluid with low thickness and big surface area (Ottino et al., 1979; Ranz, 1979). In turbulent flow systems this elongation of the fluid streams could be break due inertial forces overlapping viscous forces. In the end, inside the mould molecular diffusion enters in place, described by the Fick’s law. In fluids as the Schmidt number is big the diffusion is slow and controlled by interfacial area. Ranz (1979) used a Lagrangian reference frame fixed on a laminae to study the variation of this laminae with the striation thickness, interfacial area and the velocity gradient. A reducing of the striation thickness, consequently and increasing of the surface area augments the direction of diffusion, correspondent to the normal direction with the surface area, decreases improving molecular diffusion and therefore the mixing quality. Ottino et al. (1979) derived and equation, for Newtonian incompressible fluids, to obtained the interfacial area generation as a function of viscous dissipation. Ottino et al. (1979) also presented a set of tools to analyse lamellar structures in a microscale environment without surface tension effects and with instabilities at microscale. Chella and Ottino (1984); Ottino et al. (1979) applied the lamellar model present previously to find concentrations field in a system with multicomponent diffusion and n-linear chemical reactions. Chella and Ottino (1984) performed a parametric study for a reactor with an arbitrary flow field and segregated feeding with the lamellar model, where the parameters to study was the initial reactant segregation, the ratio of the reactants diffusivities, the first and second Damkõhler number correspondent to the relation between the mixing and reaction times for the first and the relation of mixing and diffusional time for the second number. The reactive schemes used was a second order bimolecular reactions, parallel competitive reactions and consecutive competitive reactions. Fields and Ottino (1987) used lamellar model to predict the effect of striation thickness distribution on the polymerisation. A competitive-parallel reaction scheme was used to evaluate the formation of soft and hard segments were was tested two distributions of
2 State of art of RIM and CIJ 19 striation thickness: uniform and one obtain by Kolodziej (1980). The author reported the best result was obtained by the uniform distribution rather than the Kolodziej (1980) distribution. Fields and Ottino (1987) used a limped competitive-parallel reaction scheme to study the effect of stretching of two fluids, combining the polymerization model with the lamellar model. The authors study the stretching effect where the flow is driven by shear or elongations. The authors conclude the stretching have a strong effect on the reaction system. As a conclusion the authors observed the elongation flow, when compared with shear flow, have a higher mixing efficiency Due high viscosities of the flow stream the main mixing mechanism present is the reduction of the thickness to a point where diffusion can enter to a level where polymerization can be archived. 2.2.4.2 Engulfment Deformation Diffusion Model Baldyga and Bourne (1989a, 1989b) proposed a mixing model where the mixing is produced in three steeps: Engulfment caused by vorticity followed by elongation due shear stresses and in the end diffusion. Consider two elongate different fluids in layers A and B with scales equal to the Kolmogorov scale 𝜆𝑘, Baldyga and Bourne (1984c) proposed the velocity of the reduction of the lamina thickness is a function of the energy dissipation rate and viscosity of the fluids. To model the chemical species transport is used the diffusionconvection-reaction equation. The engulfment is caused by a vortices with lamina thickness of order 𝜆𝑘 which are initially twisted in vortices with 𝜆𝑘 of diameter, but increasing until twice the size when the two fluids are incorporated. 2.2.4.3 The E Model Baldyga and Bourne (1989a, 1989b) introduced the E model assuming Schmidt number lower than 4000 and more than two engulfment are needed for complete mixing. The author proved experimental and theoretical that deformation and diffusion it is not a limiting steps of the mixing process. Experiments revel the E and EDD model can be used for calculating the product distribution for the consecutive competitive reaction schemes with no differences, although in a numerical component
2 State of art of RIM and CIJ 20 2.2.5 Mixing and hydrodynamics simulation In a numerical model the only way to study the fluid structure and the reaction mechanics is to describe the domain, in three-dimension form and discretize into small control volumes, where in each one the Navier-Stokes time dependent equation is calculated 𝜌𝐷𝑣 𝐷𝑡 =−∇𝑝+𝜇∇2𝑣 +𝜌𝑔 (2.15) And the continuity equation ∇.𝑣 =0 (2.16) Coupled with advection-diffusion reaction equation 𝜕𝐶 𝜕𝑡+𝑣 −∇ 𝐶=𝐷𝑚∇2𝐶+𝑟 (2.17) Where 𝑣 is the three dimensional flow field 𝑣 =𝑢1𝐞 𝟏+𝑢2𝐞 𝟐+𝑢3𝐞 𝟑, 𝜌 the fluid density, 𝑝 is the pressure, 𝑔 is the gravity vector, C is the concentration of a chemical species, 𝐷𝑚the molecular diffusivity and 𝑟 it is the net rate of creation or destruction of a specie by chemical reaction. The solving of the Navier-Stokes, in a three dimension time dependent resolution for a certain domain and boundary conditions is called Direct Numeric Simulation DNS (Baldyga & Bourne, 1999). Some authors take the liberty of calling DNS of a simulation, both 2D or 3D flow field where all the scales are calculated until the lower scales, Kolmogorov scaled. Although this method can obtain all the scales present inside a reactor, the used of a fine computation grid to discretize the domain make almost impractically, due time limitations, any study until this scale. The grid points required to perform this study is proposition of the ratio off all the scales and the Reynolds number, given by 𝑁𝑥𝑦𝑧 ≈(𝐿 𝜆𝑘)3∝(𝑅𝑒)9 4 (2.18) 𝐿 is The largest length scale of the problem (Baldyga & Bourne, 1999; Givi, 1989) and 𝜆𝑘 is the smallest scale in a turbulent flow, in which the viscosity overtakes the problem and turbulent kinetic energy is dissipated into heat and is function of kinematic viscosity 𝜈 and the average rate of dissipation of turbulent kinetic energy per unit mass 𝜀
2 State of art of RIM and CIJ 21 𝜆𝑘=(𝜈3 𝜀)1 4 (2.19) Due to high computation coast of a DNS, nowadays only was possible to study small and moderate Reynolds number, and for higher Re common is the industry and atmosphere the computation cost required will largely exceed the capability of the most powerful computer currently available If now consider the computational coast required to perform a DNS of the flow field but now considering also the mixing effect. It is know the scale of mixing in turbulent flows is the Batchelor scale 𝜆𝑁=𝜆𝑘 √𝑆𝑐 (2.20) And Sc is the Schmidt number defined as 𝑆𝑐=𝜇 𝜌𝐷𝑚 (2.21) Where 𝐷𝑚is the molecular diffusivity of a chemical species, 𝜌 is the density and 𝜇 is the viscosity. If is considered the mixture is liquid-liquid then the Schmidt number is order of 103 and the mixing scale is one order of magnitude smaller than the smallest hydrodynamic scale, Kolmogorov scale. Coupling this parameters in a 3D simulation considering mixing the grid scales required to perform this study is 𝑁𝑥𝑦𝑧 ∝(𝑅𝑒)9 4 𝑆𝑐3 2 (2.22) For 𝑆𝑐>1 it will be necessary a gird refining. In case of liquids, where the Schmidt is ranged of 103 then to perform a DNS to a mixing problem the grid nodes will increase in a factor of 3 ×104 Due to high computation coast of a DNS other methods was developed to study the flow and mixing mechanics up to fully turbulent flows. This methodologies used is the Reynolds average Navier –Stokes equation (RANS) now Large Eddy Simulation LES is highly used in turbulent flows to perform computation until low hydrodynamic scales
2 State of art of RIM and CIJ 22 2.2.5.1 RANS Simulation In the study of Reynolds-Average Navier-Stokes equation provide a methodology is compute turbulent flows in a certain range for the majority engineering applications. In this method the smaller scales are not computed and the flow is average in each time required. This provided a good stability in the progress of the time integration whit minimum computational cost. In this method all the variables, such as velocity field 𝑣𝑖, are divided on time average term, 𝑣𝑖 , and the fluctuation term 𝑣𝑖′ 𝑣𝑖=𝑣𝑖 +𝑣𝑖 (2.23) In this method the governing equation is solved, by replacing each variable as their average plus fluctuations in the Stokes equation. A new apparent stress emerge, known as Reynolds stress, by adding a second order tensor in the equation, which from that is possible deduce various models and adding new variables to access different engineering problems, resulting in: 𝜌(𝜕𝑣𝑖 𝜕𝑡 +𝑣𝑗 𝜕𝑣𝑖 𝜕𝑥𝑗)=−𝜕𝑝 𝜕𝑥𝑖+𝜇𝜕2𝑣𝑖 𝜕𝑥𝑗2−𝜌𝜕(𝑣𝑖′𝑣𝑗′) 𝜕𝑥𝑗+𝜌𝑔 𝑖 (2.24) Was refer early the term 𝑣𝑖𝑣𝑗 is referred to incompressible flow and named as the Reynolds stress tensor 𝜏𝑖𝑗 =−𝜌𝑣𝑖′𝑣𝑗′ (2.25) RANS methods can be one or two equation (Pope, 2000). Prantl’s, Baldwin and SparlatAllmaras are examples of one equation RANS models, while for two equation exist the Realisable, K-omega and the K-Epsilon, which is the most used method in RANS nowadays (Baldyga & Bourne, 1999; Givi, 1989; P. A. W. Libby, F.A., 1994). A insight for each RANS method is given by (Pope, 2000; Speziale, 1991). The same methodology used in RANS can be applied to the advection-diffusion-reaction equation, where the species concentration is given by 𝐶, and their time averaging 𝐶 plus the fluctuation 𝑐′ resulting in 𝜕𝐶 𝜕𝑡+𝑣𝑖 𝜕𝐶 𝜕𝑥𝑖=𝜕 𝜕𝑥𝑖[−𝑣𝑖′𝑐′ +𝐷𝑚𝜕𝐶 𝜕𝑥𝑖]+𝑟 (2.26)
2 State of art of RIM and CIJ 23 Similarly with the RANS equation it was adding a second term 𝑣𝑖′𝑐′ which is related with the turbulent diffusivity 𝐷𝑇 𝑣𝑖′𝑐′ =𝐷𝑇𝜕𝐶 𝜕𝑥𝑖 (2.27) The modelling of the Diffusivity term is mostly perform using the 𝑘−𝜀 model (Baldyga & Bourne, 1999), but others methods exist in the literature. A second order reaction is also divided into averaging term plus fluctuation of the concentration for a given species 𝑟1 =𝑘1(𝐶𝐴 𝐶𝐵 +𝑐𝐴′𝑐𝐵′ ) (2.28) 𝑟2 =𝑘2(𝐶𝑅 𝐶𝐵 +𝑐𝑅′𝑐𝐵′ ) (2.29) Where 𝑐𝑖′𝑐𝑗′ is the unclosed term and 𝑘𝑖,𝑗 is the kinetic constant for each reaction. Toor (1969), Vassilatos and Toor (1965) and Kosály (1987) presented a method to closer the second order reaction and validated with experimental data. Bourne and Toor (1977), Brodkey and Lewalle (1985), Dutta and Tarbell (1989), Li and Toor (1986) and Wang and Tarbell (1993), proposed models more closely from the reaction reality and therefore more used nowadays. A study of one of this method is presented which was proposed by Bourne and Toor (1977) and Brodkey and Lewalle (1985), respectively. 𝑐𝐴 ′𝑐𝐵′ =𝐼𝑆 𝐶𝐴0 𝐶𝐵0 (2.30) 𝑐𝐴 ′𝑐𝐵′ =𝐼𝑆 𝐶𝐴0 𝐶𝐵0𝐶𝑅 𝐶𝐴 (2.31) Where 𝐼𝑆 is the intensity of segregation and 𝐶𝐴0 and 𝐶𝐵0 are the mean concentration of specie A and B, respectively. Other method is used to calculate the unclosed terms. The Probability Density Function PDF methods is particularly used in combustion (Givi, 1989; Pope, 1985, 1994) but have also proved their value in chemical reaction engineering (R. O. Fox, 1992, 1998; R. O. G. Fox, M.R., 1993; Tsai & Fox, 1994). The main difference of this method from the other method presented is there is no transport of the variables, but by using a subgrid-scale a PDF of this variable is transported (Tsai and Fox (1994). This approach was employed by Li and Toor (1986) for a consecutive reaction CFD based.
2 State of art of RIM and CIJ 24 Initially the flow is calculated normally using a generic grid and a method suited for the problem in study, where information of the fields (velocity, pressure, temp and others in interest) are extract to a subgrid domain. From this subdomain the chemical reaction will be transported using PDF for all the subdomain in study similar to a post processing stage. This method requires further development (Baldyga & Bourne, 1999). 2.2.5.2 Large Eddy Simulation Large Eddy Simulation (LES) it is become a very important tool to simulate turbulent flows. In Kolmogorov (1991) theory indicates the large eddies of the flow field are, for confined flows, function of the geometry, but by reducing the scales their effect start to be reduced until, in the lower scale, the kinetic energy is dissipated due viscosity. This feature allows to explicitly solve large eddies and afterward calculate the small eddies in a subgrid scale. This method is formulated in three major steps (Ghosal & Moin, 1995) The first step consists in filtering small spatial scales less than the computation grid. The result equation is used to describe the temporal and special evolution of the large eddies with the effect of the subgrid scale stress tensor from the unresolved small scales into the resolved scales. The second step is substitute the subgrid stress tensor by a model which can be function of the resolved scales and tuning parameters. The last step will performed numerical simulation of the large scale field into a grid small enough to resolved the smallest scales resulting from the large eddies, but much larger than the lower theoretical Kolmogorov scale. This method is not such computational cost of a DNS, where all the scales of the flow are calculated and it is possible to visualize but also do not provided only information of the flow average behaviour as the RANS method, but instead can give a clear picture of the important scales, where the large scales are accurate simulated while also present spectral characteristics of the flow field presented in the smallest scales. For this reason this method is used in wide engineering problems where a picture of the scales from 6𝜆𝑘 up to 60 𝜆𝑘 is needed. In terms of computational cost this method will requires less effort when compared with DNS but still much larger when compared with RANS. Due to stability issues this method required an excellent initial solution therefore many author’s start initially to calculate the large scales using a RANS model and import as initial solution for LES. The same procedure can be by done by resolving in a large grid the Navier-Stokes equation. For mixing proposes the LES can be used with the advection-diffusion-reaction with a low pass filter resulting in a unclosed problem. In the same way as the eddy viscosity is is
2 State of art of RIM and CIJ 25 treated in LES it can be used the same way for the diffusivity eddy and therefore solving the advection-diffusion problem (Givi, 1989). For a reaction problem since the reaction scale, correspondent to the Batchelor scale presented previously, produce small scales than the hydrodynamic, it is possible to couple in a sub-grid domain others methods such as PDF or RANS (Givi, 1989) (deBruynKops & Riley, 2001) to solve this issue 2.2.5.3 DNS of turbulent flows Due to high computational cost off a DNS for turbulent flow this method is still very limited for small domains our low Re as reported previously. This method allow understand and determine the flow structures for all the mixing scales, from largest to the lowest, reproducing faithfully the flow dynamics by directly solving the Navier-Stokes equations for a given domain, without any simplification or model being used and therefore is the only way validate numerically the mixing models used nowadays. (Givi & Mcmurtry, 1988) study a second order bimolecular reaction (Toor, 1969) (Kosály, 1987) revealing when compared with experimental data an advantage of knowing the instantaneous values for the chemical in space and time domain. Other advantages of a DNS into study the hydrodynamics of a chemical reaction is the parametric study of a particular reaction, in where the factors which affect the mixing can be study singular, or parallel with complex flow dynamics. A particular case is present in Chakrabarti, et al. (1995) where was observed the flow dynamics of a cubical domain with periodic boundary conditions for study the fluid flow inside a Plug Flow Reactor PFR, (Gerlinger, Schneider, Falk, and Bockhorn (2000), a 2D DNS simulation using a pseudospectral method to for Re equal to 2000 with different flow structures, where 2 and 3 vortices and fully turbulent flow existent in the domain. The replication of exact conditions for the flow structures, creation of vortices, structured flow, and the position distribution and dimension of such structures creates a big advantage in using numerical methods to study such flows, where due to unpredictable behaviour of such flows their replication in laboratory will be extremely difficult being basically a lucky experiment. 2.2.5.4 DNS of chaotic flows In mixing process a chaotic flow it is a flow which under a tracer experiment and subjected to flow field will produce complex and unpredictable fractal structures with a time dependent exponential growing. In fact the term chaotic flow is often used to indicate a flow, below the turbulent flow regime and where in mixing process can be an alternative
2 State of art of RIM and CIJ 32 representation of the flow structures formed in this type of mixers. They calculated the intensity of segregation for different Reynolds numbers and found that the intensity of segregation declines abruptly for Reynolds numbers between 100 and 300. From Re = 300 to Re = 500 the intensity of segregation continues to decrease before reaching a bottom. Unger and Muzzio (1999) also compared the mixing performance of two impingement jets geometries: directly opposed jets and the jets pointing to the close end of the chamber with 8º deviation from the vertical axis. They concluded that up to Re = 300 the intensity of segregation was lower for the chamber with opposed jets, whereas for higher Reynolds numbers the intensity of segregation is slightly lower in the asymmetric jet geometry. In Unger and Muzzio (1999) work, the region of sharp transition was offset to higher Reynolds numbers, most probably due to the differences between the geometry studied and typical RIM mixing chambers, namely the ten times larger dimensions. The flow visualisation studies based on LIF provided a better understanding of the phenomenon underlying the sharp increase of mixing efficiency at Reynolds numbers between 100 and 200: the flow evolves from a steady symmetric state with high segregation between both streams to complex flow patterns that generate mixing under laminar regimes. 2.3.1.2 Measurement of Hydrodynamic Variables For an accurate characterisation of the flow field, quantitification of hydrodynamic variables is required, which is not attainable with the visualisation techniques previously described. Moreover, with a quantitative description of the flow field, namely the velocities field, it is possible to gain further insight into the effects of relevant parameters on impingement mixing, such as the Reynolds number. The following techniques, which will be considered in the next sections, have been applied in experimental hydrodynamic studies in impinging jets reactors: Laser Doppler Anemometry, LDA Johnson (2000a, 2000d); Johnson et al. (1996); R. J. Santos, A. M. Teixeira, M. R. P. F. N. Costa, and J. C. B. Lopes (2002); R. J. Santos, Teixeira, and Lopes (2005); R. J. N. d. Santos (2003); A. M. Teixeira (2000); André M. Teixeira, Santos, Costa, and Lopes (2005); Wood et al. (1991) Particle Tracking Velocimetry, PTVZhao and Brodkey (1998a, 1998d) ,Particle Image Velocimetry, PIV R. Santos, A. Teixeira, M. Costa, and J. Lopes (2002); Ricardo J. Santos, Erkoc, Dias, Teixeira, and Lopes (2008); Unger et al. (1998)
2 State of art of RIM and CIJ 33 A. Laser Doppler Anemometry (LDA) LDA is used to measure the fluid velocity at specific single points of the mixing chamber. The spatial distribution of the average velocities and of its statistics, such as the standard deviation, is then assembled from multiple measurements at different locations. Despite allowing high data rates and extensive time series, is not possible to obtain the instantaneous velocity fields or the circulation paths of the fluid inside the chamber under dynamic flow fields with LDA. Johnson et al. (1996) were not able to characterize the flow field with LIF for a Reynolds number of 300 since the jets were the only observable structure. As a result, they employed LDA to assess the influence of the geometric parameters at this regime. By measuring the vertical velocity at sucessive planes along the chamber axis by LDA and subsequently calculating the averaged velocity component in the chamber axis, Johnson et al. (1996) obtained the average recirculation patterns within the mixing chamber. R. J. Santos et al. (2002); A. M. Teixeira (2000) reported the detailed flow field characterisation inside a RIM machine mixing chamber by means of LDA measurements and quantified the effect of the Reynolds number on the flow field. From the distribution of the average velocity components, measured in normal planes to the mixing chamber axis within a thin measurement grid, it was found that the averaged flow field exhibited no variation with the Reynolds number. R. J. Santos et al. (2002); A. M. Teixeira (2000) also computed the distribution of higher statistical moments, namely the pdf of the lateral velocities and the Reynolds stresses, from the LDA measurements. They found that the flow fluctuations intensify rapidly with the Reynolds number up to Re = 100 but afterwards display a slight evolution only. Additionally, regions of higher velocity fluctuations were seen to occur near the jets impingement region and to increasingly spread throughout the chamber for higher Reynolds numbers. Wood et al. (1991) performed LDA measurements to characterise the oscillatory behaviour of the impinging jets in a RIM machine mixing chamber. The axial and lateral velocity components were measured at 4 jets diameters downstream of the impingement point at a data rate of 300 Hz and the Strouhal number was determined at Re = 125 for two different fluids with viscosity of 27 and 55 mPa.s. Wood et al. (1991) showed that, for the same Reynolds number, the Strouhal number is constant, irrespective of the velocity of the jets.
2 State of art of RIM and CIJ 34 Johnson and Wood (2000) studied the typical frequencies in opposed jets mixing chambers through LDA measurements. Two mixing chambers geometries were considered, one with a circular cross-section and the second with a parallelepiped shape, for which three injectors diameters were also tested. Johnson and Wood (2000) computed the Strouhal number for the Reynolds number at which the flow field becomes time dependent, that is Re ≈ 100, and for the limit Reynolds number where oscillations with distinct frequency still persist, namely for 130 < Re <220. It was not feasible to calculate the Strouhal number for higher Reynolds numbers as the oscillations became markedly random and without a defined frequency. Johnson and Wood (2000) noticed that the Strouhal number increases with the Reynolds number from the lowest values where flow field oscillations occur to the highest values where the oscillations still have a defined frequency. In addition, it was concluded that the Strouhal number is affected both by the ratio of the injectors and the chamber diameter. Johnson (2000a, 2000d) studied the effect of the jets momentum ratio with LDA and found that, for momentum ratios different from one, the impingement point moves towards the nozzle where the jet has less momentum and the flow field also displays more stability. R. J. Santos (2003) determined the Strouhal number from LDA measurements for Reynolds numbers ranging from 250 to 600. In addition to a LDA probe with a high area of signal capture that achieves high data rates, data processing was employed to attain defined frequencies at these high Reynolds numbers regimes. In general, researchers have faced difficulties in finding the typical frequency above the critical Reynolds number, due to the intensification of the random behaviour of the jets oscillations (Johnson & Wood, 2000) However, the technique used by R. J. Santos (2003) made possible to establish that for 250≤Re≤600 the jets oscillate with frequencies around typical values which scale with the Reynolds number, while the Strouhal number is kept almost constant, meaning the main characteristics of the flow field are retained. Accordingly, R. J. Santos (2003) pointed out that the dependence of the Strouhal number on the Reynolds number only occurs at the first transition regimes to dynamic states. R. J. Santos (2003) also showed that the dependence of the Strouhal number on the injectors and chamber diameters ratio stems from the flow structures that are formed for the different geometries, for which specific jet oscillations frequencies develop. Additionally, R. J. Santos (2003) studied the influence of the jets momentum ratio on the flow field dynamic behaviour and concluded that if the
2 State of art of RIM and CIJ 35 momentum ratio is different from 1 the jets do not oscillate, regardless of the flow regime being above the critical Reynolds number. B. Particle Tracking Velocimetry (PTV) PTV is a techique used to compute the velocity, in a Langragian approach, by measuring the location of particles in three-dimensional space as they move under the influence of the flow field. Zhao and Brodkey (1998a, 1998d) obtained the averaged flow field in an opposed jets precipitator chamber from PTV measurements with water at Re = 200 and Re = 4000. The average velocity field determined by for the unsteady flow in the chamber at Re = 200 was similar to the steady state reported by other authors such as Wood et al. (1991). However, in the velocity maps where time average was done over smaller periods the flow field dependence with time was undoubtedly detected as the flow structures present were not the same as the ones found in the long time averaged sequences. As for the studies at Re = 4000, it can be seen that both the flow patterns and the overall circulation patterns change entirely. The unsteady nature of the flow was evidenced by Zhao and Brodkey (1998d) C. Particle Image Velocimetry (PIV) PIV is an Eulerian method that measures the distribution of the instantaneous velocities of the flow field in a 2D plane, allowing to represent the flow structures without the influenced of mixing effects like in the case of tracer techniques. Unger et al. (1998) characterised the flow field in an opposed jets chamber typical of precipitators with PIV and glycerine aqueous solutions for Re≤80 only, as the flow starts to evolve into unsteady flow for higher Reynolds numbers. From the velocity vectors maps of the flow field, it can be noticed that the jets begin to impinge when the Reynolds number exceeds 10 and recirculation zones are also formed for those regimes.
2 State of art of RIM and CIJ 36 2.3.2 Hydrodynamics in CIJ - Numerical Studies Computational fluid dynamics numerical codes have the capability to full characterize a reactor functioning, ranging from the hydrodynamic study to more complex analyses such as reaction studies, multiphase studies and others. CFD is used in RIM machines to major quantify the mixing mechanism from the study of reaction and hydrodynamic as also the reaction parameters such as RTD, pressure drop and other crucial in any project. In the beginning of the development of RIM CFD analyses of a 3D domain of the mixing chamber was impossible due computational limitations and this is why several simplifications and approaches was perform to approximate the experiment to the numeric work. One big advantage of CFD simulations is the knowledge of all the variables in any single point of the computational domain, making the possibility to study the reactor in a macro and microscopic scale, depending of the discretization used. In this subsection will be presented the studies performed for a steady state flow regime and for a time-dependent flow regime 2.3.2.1 Mixing Models employed in CIJ The models used to access the mixing performance of a CIJ is a lamellar model, since the flow inside the mixing chamber is driven by viscosity, therefore a way to perform the mixing quantification is by calculate the stretching of particles and their striation thickness L. J. Lee et al. (1980) develop an expression to calculate the striation thickness base on a stretching function, where this function correlate the particles stretching versus the viscous energy dissipation. The validation of the equation was perform for a range of Reynolds number and the geometric parameters of the chamber. 𝑠=(𝑛𝑑𝐷3 1+𝑟𝑠−1)1 4𝑅𝑒−3 4 (2.33) C.L. Tucker and N.P. Suh (1980) formulate the final particle striation thickness is related to the size of the smallest eddies and therefore function of the Reynolds number for a fix geometry 𝐼𝑠∝𝑅𝑒−3/4 (2.34) And where the striation thickness can be calculated by
2 State of art of RIM and CIJ 37 𝑠∝∝ 𝑑 𝑅𝑒3 4 (2.35) Baldyga and Bourne (1983) developed a model to predict the striation thickness for T-Jet’s mixer based on the theory of turbulent diffusion and validate the model with experimental work performed by Kolodziej, Macosko, and Ranz (1982). The model reproduced a good agreement for striation thickness less than 10𝜇𝑚, value for which a good mixing is achieve for polymerization can occurs Fonte (2012) used a lagrangian approach so study the stretching of particles in a 3D domain and consequently the striation thickness by assuming a conservation of a particle volume. It conclude for the Reynolds in study the CIJ produce a striation thickness small enough to ensure a controlled polymerization in the mould. 2.3.2.2 Steady CFD Simulations In a steady state regime a 3D mixing chamber used as precipitator simulation was perform by Zhao and Brodkey (1998a, 1998d), where by presenting the flow vector map for a plane aligned with the chamber axis a recirculation vortex was observed for the top and bottom of the chamber although the main flow it is conducted through the chamber axis. For the same precipitator described early Unger et al. (1998) the hydrodynamic behavior of the flow field where observed a transition of regime located in Re equal to 80 achieve a self-sustainable chaotic flow regime. Using the methodology of the shear rate the author observed a linear relationship between the Reynolds and the strain rate, even in the zones where high shear stress occur, at the impingement point and outlet. A study of the jets imbalance with the intensity of segregation shows a slight change in the position of the stagnation point in 1º can reduce the intensity of segregation in 80% by increasing the recirculation zones. 2.3.2.3 Unsteady CFD Simulations Dynamics CFD analyses is used in RIM mixing chamber to study the flow mechanism involved in the mass transfer, simply from advection terms with only resolved NavierStokes equation, or when the advection terms are couple with reaction equation, permits a full study, comparing with experimental data the perform and optimization of the mixing chambers. Due the unstable nature of the flow and the necessity to quantify the flow,
2 State of art of RIM and CIJ 38 steady CFD simulations didn’t provide enough tools needed to provide a time-dependent flow characterization. Wood et al. (1991) perform a 3D CFD simulation for a transient flow and compared with visualization studies and LDA described in the previous subsections. The results of the simulation prove, when comparing with LDA and visualization studies a good agreement. The axial velocity proved consistent with LDA and the flow pattern when compared with visualization studies shows similarity. For the Reynolds in study the oscillatory movement of the jets was study using a frequency spectra and also the Strouhal number was study for different geometries dimensions. A. M. Teixeira (2000) perform a unsteady 2D CFD simulation of CIJ’s for Reynolds range from 100 to 500. The main conclusion was the self-sustainable chaotic regime was archive for Re > 250 and, for Re < 250 the flow tends to be segregated. The study of the flow perturbations bellow the critical Reynolds shows stability even when exited, but with the increasing of Reynolds above the critical value shows the flows will never tend to be segregated. Johnson (2000a, 2000d) study the effect of the jests momentum imbalance in the flow field. He demonstrate a slight change of the jets momentum will cause a deviation of the stagnation point from the chamber axis. This work proves the conservation of jets momentum with different geometries dimensions. Testing different injector diameters but with the flow rate compensated to maintain the same momentum is proven to maintain the impingent point in the center. R. J. Santos (2003) presented numerical simulation in CIJ to study the effect of unequal momentum ratio in the mixing performance where ti conclude by study the mass transfer efficiencies using a reaction model it is possible the system is not sensible to slightly changes in the jets momentum ratio. Bierdel and Piesche (2001) study for a 3D mixing chamber in Reynolds equal to 50, which is know it produces a segregated state. By inducing oscillatory perturbation in the jets velocity unsteady flow field is created for a Reynolds bellow the critical. The mixing quality study was perform using particle tracking, where the author conclude the mixing efficiency increases with the Reynolds number and with the frequency of the timedependent oscillatory perturbation. In the end this results was compared with visualization and tracer experiments.
2 State of art of RIM and CIJ 39 2.3.3 Mixing and Chemical Reaction in CIJ Until now the focus was present the work perform in RIM and their mixing chamber CIJ’s, where initially the hydrodynamic behaviour of the flow was study using experimental and numerical tools. It is know the hydrodynamics greatly affect the mixing degree and therefore their study is critical to ensure a good mixing of the fluid flow stream and consequent a good polymerization inside de mould. In this subsection will be addressed the mixing study’s in CIJ, studied by experimental and numerical tools. 2.3.3.1 Experimental Studies A. Tracer Experiments Tracer experiment is an old technique to quantify the mixing performance of any reactor used nowadays. This technique is very used due their simplicity and reduced cost when compared with another’s. S.C. Malguarnera and N.P. Suh (1977); Salvatore C. Malguarnera and Nam P. Suh (1977) used glycerin as a fluid and acetic acid solution as a tracer in one of the fluid streams where a several samples was perform and the standard concentration was compute and used as mixing quality. The main objective is to perform a dimensional analyses of CIJ understanding the main factor which can affect the mixing mechanism. The author identifies some key aspects: The Reynolds number calculate at one injector was ranged from 50 to 500 where the author conclude o The standard deviation of the tracer concentration as 2.5 times lower comparing Re=100 with Re=50 o For Reynolds lower than 50 no mixing occur o The moment ration between both streams should be the equal for effective mixing, and calculated as 𝑅𝑀=𝜌1𝑣𝑖𝑛𝑗1 2𝑑1 2 𝜌2𝑣𝑖𝑛𝑗2 2𝑑2 2 (2.36)
2 State of art of RIM and CIJ 40 Where 𝜌 is the fluid density, v the mean injector velocity and d the injector diameter. The subscript 1 and 2 corresponded to the left and right injector respectively According the geometrical characteristics of the mixing chamber the author concluded o There is no effect of the runner length in the mixing quality o The chamber diameter don’t affect directly the mixing performance, although the volume of the chamber implies: for bigger chamber have difficulty of expelling the air inside the chamber and for smaller volumes a tight control of the flow is required. o The distance of the injector to the top of the chamber cause no impact on the mixing quality C.L. Tucker and N.P. Suh (1980) measure the mixing quality by quantifying the spatial variance of the absorbance of the tracer on the outlet. From this method it is observed an exponential decreasing of the absorbance with the Reynolds number by a factor -9/4 and regarding the study of the moment ratio it was obtained no changes on the mixing efficiency up to values of 2.5. The author’s conclude the decreasing of the scale of segregation would be related with the decrease of the smallest scales on the flow field. From this work some important notes should be said. According the mixing efficiency with the moment rate ratio present a opposition with the previous work form S.C. Malguarnera and N.P. Suh (1977); Salvatore C. Malguarnera and Nam P. Suh (1977). The second note is regarding the decreasing of the scale of segregation with the flow scales by a factor of -9/4. It is know the smaller cases in turbulent flow Tennekes and Lumley (1997) varies as: 𝜆𝑘∝𝑅𝑒−3 4 (2.37) In fact this statement cannot be very accurate since the flow regime is far from turbulent and therefore, in a three dimensional geometry, it cannot relate the flow scales with the smaller scales. In fact (R. J. Santos, 2003) present a power spectra where prove the flow scales did not have such big variation by the Reynolds is increasing. According the study of the moment ratio with the mixing quality it can be said the method did not guaranty the separation of the influence in the practical effect in the mixing mechanism inside the chamber. Face of what is described and because the work present others error which cannot be denied and therefore the conclusions about the influence of the moment rate ratio with the Reynolds number should be carefully attended.
2 State of art of RIM and CIJ 41 B. Adiabatic Temperature Rise One process to measure the mixing uses the measurement of the rising temperature in the experiments with the exothermic reactions in case of a mixing constrain reaction. Lee et all (1980) study the adiabatic temperature rise of the polymerization reaction in where various Reynolds number were studied whit some geometrical parameter. The author conclude according the effect of the Reynolds number it was observed a increase of the temperature rise up to Reynolds 200 although the bigger change had occur between Reynolds 100 and 160. According the second part of the study, the effect of geometrical parameters namely the ratio of the injectors’ diameters, the camber diameter and the jets alignment and direction. It was observed for all the geometrical parameter present the same critical Reynolds number in which the adiabatic temperature rise is constant. It was also verified that directly opposed jets present lower efficiently. This study is is very important since present some conclusion where at date it was new, but since the experiment was performed for jets momentum ratio different than one then the is necessary carful attention when compared with the results of this thesis since ratios different than one complete changes the reactor hydrodynamics and therefore the mixing performance. Sebastian and Boukobbal (1986) studied the stoichiometric reaction between polyol and isocyanate the effect on the momentum ratio in the adiabatic temperature. For the same stoichiometric ratio, where for equal injector diameter the moment ratio is 2.42, it was perform a changing in the moment ratio by changing the injectors diameters. It was observed the rise of the adiabatic temperature is was up to Reynolds 200 where for each Reynolds studies the maximum temperature rise was for equal momentum ratios Harris, Anderson, and Shannon (1992) a study to quantify the effect on the Reynolds number into the adiabatic temperature of a polymerization reaction. The authors observed a increase in the adiabatic temperature until Reynolds 200, proving the mixing performance is only improve until that Reynolds reinforcement the idea of concentering the study of the mixture efficiency by studding the mixing mechanism and reaction until that particular Reynolds.
3 CFD Model and Simulation of CIJ 48 steady state simulation is performed since it is known from previous works that a fully segregated flow without any chaotic flow motion is obtained in these conditions. Therefore, for the steady state flow regime the symmetry plane is activated and set with zero shear forces in each direction allowing free flow motion in one half of the mixing chamber. Conversely, when a transient simulation is performed the symmetry plane is removed, thus allowing flow motion in the whole chamber domain. 3.2.2 Domain Discretization The physical domain is discretized into small control volumes where the NS equation will be solved. The construction of the geometrical domain and its discretization to create a simulation grid were performed using the DesignModeler and the Meshing software included in the ANSYS® v15 suite. The entire domain was discretized into small sections in order to manually ensure mesh control, both in cell size and meshing method, generating 2.3 million uniformly sized hexahedral finite volumes with the maximum edge length of ∆𝑥=100 𝜇𝑚. The mesh generator automatically adapts the mesh to discretize curved surfaces as well as maps the face cell between each section so that a perfect cell match is attained. A few triangular cells are used in order to match the injector hexahedral cells with the chamber hexahedral cells. Figure 3.2 shows a plane along the chamber and the injectors axes of the 3D hexahedral mesh in a cut passed in the injector’s axis. It can be observed the presence of matching quadrilateral cells along the domain, which means a good information transfer between each cell since there are not present discontinuities the mesh domain. Triangular cells is used to match the cell region form the injectors with the cell region form the chamber due to the fact a no using of this element’s will force high distorted elements in the region which is worse than have non align face cells with the primary direction of the flow. This effect is clearly visualize in the velocity contour where small numeric diffusion of the velocity is observed, particularly align with the injector’s walls. Since a parabolic profile exist in the injector, meaning high velocity in the center of the injector’s, then the consequence for the oscillatory movement of the jets is primarily due the high velocities, correspondent in the injector axis, rather than low velocities, and therefore the using of this cells will not have a big impact in the end results
3 CFD Model and Simulation of CIJ 49 Figure 3.2 – Hexahedral mesh with 𝟏𝟎𝟎 𝝁𝒎 used in the simulations (left). Detail of the triangular cells used to couple both meshes avoiding distorted elements (right) 3.2.2.1 Grid size The choice of the maximum edge length is chosen according to the Kolmogorov length scale which corresponds to the lowest scale in a turbulent flow and varies with Reynolds number as 𝑅𝑒−3/4. Despite the Reynolds in study is laminar and far beyond the transition region choosing the lower scale as the Kolmogorov ensures all the hydrodynamics length scales are being simulated without a need of a sub-grid turbulence modelling as the one used in LES. Since the Navier-Stokes equation are directly simulated for all the possible length scales existent for a particular hydrodynamic simulation this is referred in the literature as Direct Numerical Simulations (DNS). The Kolmogorov length scale 𝜆𝐾 in isotropic and homogeneous turbulent flow can be estimated as 𝜆𝐾=(𝜐3 𝜀)−14 (3.1) Where 𝜐 is the kinematic viscosity and 𝜀 the turbulent energy dissipation rate which can be calculated by integrate over a boundary surface 𝜕Ω neglecting the contribution of the friction in the walls, resulting in:
3 CFD Model and Simulation of CIJ 50 𝜀≈∯ |𝒖|2 2𝜌 (𝐮.𝐧) 𝑑𝑆 𝑑Ω (3.2) Where the vector 𝐧 represent the normal vector’s to the geometry axis. Consider the chamber volume where dissipation occur as 𝑉𝑑𝑖𝑠𝑠𝑖𝑝𝑎𝑡𝑖𝑜𝑛=(𝜋4𝐷2)(𝑛𝐷) (3.3) Where 𝑛𝐷 correspond the distance between the top of the chamber until the point where dynamic mixing occur, which can take value as 𝑛𝐷≈2𝐷 to 3𝐷 (kolodziej et all.,1982). A new variable can be assumed as the average density weighted by the volumetric flow rate Q, defined as 𝜌=𝑄1𝜌1+𝑄2𝜌2 𝑄1+𝑄2 (3.4) Where the indices 1 and 2 correspond to the left and right injector. This weighted density can be also extended for the average viscosity, 𝜇, defined as: 𝜇=𝑄1𝜇1+𝑄2𝜇2 𝑄1+𝑄2 (3.5) The variation of kinetic energy can be neglected since due to the areas ratio, the average velocity in the outlet is about 5% of the average velocity in inlet. Consider the flow in the injector is laminar resulting in a parabolic velocity profile, and the introduced new variables, the turbulent energy dissipation rate can be now expressed as: 𝜀=∫[2 𝑢𝑖𝑛𝑗,1(𝑑1 2)2−𝑟2 (𝑑1 2)2]𝜌1 𝜋 𝑟 𝑑𝑟 𝑑1 2 0+ ∫[2 𝑢𝑖𝑛𝑗,2(𝑑2 2)2−𝑟2 (𝑑2 2)2]𝜌2 𝜋 𝑟 𝑑𝑟 𝑑2 2 0 𝜌 𝑉𝑑𝑖𝑠𝑠𝑖𝑝𝑎𝑡𝑖𝑜𝑛 (3.6) Where 𝑟∈[0,𝑑2] is the radial position at the injector’s circular cross section. This equation can be introduced in the Kolmogorov length scale formulation, resulting in: 𝜆𝐾=(𝐷3𝑛𝑑 1+1 ∅𝐾)14𝑅𝑒−34 (3.7)
3 CFD Model and Simulation of CIJ 51 Where the kinetic energy rate ratio can be calculated as: ∅𝐾=𝜌1𝑑12𝑢𝑖𝑛𝑗,1 2 𝜌2𝑑22𝑢𝑖𝑛𝑗,2 2 (3.8) For a isotropic and homogeneous turbulent flow, meaning there is no mean shear due rotational our boundary effect or there is no flow gradients, the Kolmogorov scale correspond to the lower hydrodynamic scale in the flow and if the grid spacing is lower than this value for Navier-Stokes can be applied means a full DNS with all the flow scales is numerically calculated. Libby (1996) say acceptable accuracy in DNS is obtain if grid spacing is between four to six times greater than the Kolmogorov scale. In the present work and for the largest Reynolds in study, 300, it is obtain for n=3 Kolmogorov scales a value closer to 95 𝜇𝑚 which is a bit lower than the computational grid used but inside the Libby theory range. In this particular case only the lower scales will not be accurate calculated but, as it will be seen forward this scales represent to the dissipation range which don’t affect he major effect in study which occur for lower space scales. For the lower Reynolds in study, [150, 200 and 250] it produces scale values bigger than the computational mesh meaning in this case all the scales it will accurate calculated. The mesh size chosen it guaranties a full representation of the hydrodynamic phenomena for all the Reynolds number in study. 3.2.3 Governing equations Mathematical equations which describe the flow motion was originally introduced by Euler (18th century) described the conservation of mass (continuity) and balance of momentum and energy for incompressible fluid (or compressible if the divergence of the flow field is zero) without viscosity. Latter Navier-Stokes equation (19th century) described the rate of change of momentum at each point in a viscous fluid with constant density 𝜌 and constant kinematic viscosity 𝜈 𝜕𝐮 𝜕𝑡+𝐮∙∇𝐮=−𝛻𝑝 𝜌+𝜇∇2𝐮 (3.9) With ∇∙𝐮=0, representing the conservation of mass equation evaluated for the boundary of the problem consider incompressible fluid. 𝐮(𝐗,t) is the fluid velocity field in a Lagrangian reference frame. 𝐏(𝐗,t) is the pressure field resulting from the preservation of incompressibility of the fluid. This equation is similar to the Newton’s law in which force is equal to mass times acceleration. In the left hand side of the equation represent the
3 CFD Model and Simulation of CIJ 52 acceleration of the fluids particles and on the right hand side is the sum of the forces per unit mass in a control volume: the pressure force and the viscous forces arising from momentum diffusion though molecular collision (Ecke, 2005). The nonlinear term 𝐮∙∇𝐮 describes the adjective transport of the fluid momentum which produce a particular phenomenon. The Navier Stokes equation due the non-linearity of their equation are very sensible to the initial conditions which can translate into complete different solutions by s slight change in the initial conditions. This effect was visualized by Edward Lorenz (1972) in atmospheric flows introducing a new field: the Chaos theory. Changing the flow regime can produce smooth laminar flow to more complex flow motions quantified by their length or time scale. Reynolds study this concept in a pipe flow where visualize the concept of the non-linear term of the equation and by introducing a new flow parameter, the Reynolds number, could according to this an-dimensional number predict the dominance of the non-linear term over the viscous term, which is the case of turbulent flow motion. The Reynolds number is calculated as 𝑅𝑒=𝑖𝑛𝑒𝑟𝑡𝑖𝑎𝑙 𝑓𝑜𝑟𝑐𝑒𝑠 𝑣𝑖𝑠𝑐𝑜𝑢𝑠𝑓𝑜𝑟𝑐𝑒𝑠=𝜌𝑣𝐿 𝜇 (3.10) For CIJ’s is defined as 𝑅𝑒=𝜌𝑣𝑖𝑛𝑗𝑑 𝜇 (3.11) Where v is the average velocity of the fluid, L the characteristic linear dimension and 𝜌 and 𝜇 the fluid density and dynamic viscosity, respectively. For CIJ’s the characteristic dimension is the internal diameter of the injector and the velocity is the average velocity for that injector. In the case in study the flow is far from turbulent and situate in the laminar region meaning in the fluid flow the viscous forces over dominates the inertial forces produced by the non-linear term. The work perform so far clearly illustrate that fact Santos (2003) and therefore as it can be seen from previous studies (Teixeira, 2000; Santos, 2003; Fonte, 2012) after a certain Reynolds, critical Reynolds number, exist in time and spatial domain a big variation of the flow field inside the mixing chamber which since the flow is driven by viscosity is not caused by turbulence effects but only caused by the highly chaotic nature of the flow.
3 CFD Model and Simulation of CIJ 53 3.2.4 Boundary conditions and initial conditions 3.2.4.1 Injectors The fluid enter in the chamber in a normal direction of the chamber axis in the injector’s section. The injectors have a small length which will going to allow the introduction of a steady full developed, laminar flow with a parabolic velocity profile, Couette flow, in a cross section of the injector. The equation to describe the parabolic profile is 𝐮=−2𝑢𝑖𝑛𝑗((𝑑2)2−𝑟2 (𝑑2)2) 𝐧 (3.12) Where 𝑢𝑖𝑛𝑗 stands for the mean velocity in the cross section, 𝐧 the normal direction of the flow and 𝑟 is defined as 𝑟=𝑟(𝑥,𝑧)=√𝑥2+𝑧2 (3.13) Additionally the velocity gradient is assumed as null in each injector resulting in introducing zero pressure to not include addition flow forces which will going to affect the flow behaviour. 3.2.4.2 Walls For all the wall of the mixing chamber the velocity component are set as zero producing a non-slip wall. In this problem the walls are consider without roughness 3.2.4.3 Outlet The outlet boundary condition is used assuming parallel flow in the boundary with velocity gradient null resulting in 𝑣𝑥,𝑧=0 (3.14) 𝜏𝑦𝑦=−𝑝0 (3.15) Where 𝜏𝑦𝑦 is the normal total shear stress in the direction of the chamber axis and 𝑝0 the pressure at outlet. For a Newtonian fluid 𝜏𝑦𝑦 is defined as 𝜏𝑦𝑦=2𝜇𝑑𝑣𝑦 𝑑𝑦−𝑝 (3.16)
3 CFD Model and Simulation of CIJ 54 Where 𝜇 is the fluid viscosity and 𝑝 the pressure. In the outlet 𝑝=𝑝0 and 𝜏𝑦𝑦=−𝑝0 resulting in 𝑑𝑣𝑦 𝑑𝑦=0 (3.17) 3.2.4.4 Steady state Steady state solution is archive in order to obtain an initial solution for transient flow regime. This simulation is obtain by setting the time dependent term of the NS equation as zero which is allowed iterate until the residuals set minimum as 10−6 in the continuity and momentum equation. To ensure a fast process and since in a steady state flow regime both fluids streams remains segregated, as previous studies demonstration, then a symmetry plane, demonstrated in Figure 3.1, is set as “on” with no shear stress between the flow and each face of the symmetry plane 3.2.4.5 Transient state The previous steady state simulation provide an initial for the transient nature of the flow, by introducing in the Stokes equation the time derivative term. The symmetry plane used in the previous run was disable allowing passage of the bow between one side of the chamber to another side. To start a transient simulation a non-symmetric perturbation is introduced in one injector while the opposite injector remained the same flow rate. Such perturbation is directly introduced in the velocity profile, using a User Defined Function (UDF) increasing the injector’s velocity. In the right injector the velocity profile was replaced by 𝐮=−2𝑢𝑖𝑛𝑗((𝑑2)2−𝑟2 (𝑑2)2)𝑓(𝑡/𝜏) 𝐧 (3.18) Where 𝑓(𝑡/𝜏) is the perturbation time function, depending on the flow time 𝑡 and the chamber mean residence time, 𝜏, defined as: 𝜏= 𝐷2𝐻 2𝑢𝑖𝑛𝑗𝑑2 (3.19) The methodology of implementing a perturbation time function was previously performed by Teixeira (2000), Santos (2003) and Fonte (2012) where it was study in order to ensure a
3 CFD Model and Simulation of CIJ 55 introduction of a discontinuity in the flow field. Therefore a smooth time function with a continuous first order derivative was chosen 𝑓(𝑡𝜏)=1+𝑎2(1−𝐻(𝑡 𝑡𝑎𝑢−𝑏) ) (1−cos(2 𝜋 𝑏𝑡 𝜏)) (3.20) Where 𝐻(𝑡 𝑡𝑎𝑢−𝑏) represent the Heaviside function. The parameter a is defined as the maximum amplitude of the perturbation, set as a=0.1. Parameter b represent the time, according the chamber mean residence time in where the perturbation is applied, set as b=0.121. Figure 3.3 represent graphically the perturbation function where it can be observed the perturbation have very small effect on the solution since 𝑣𝑖𝑛𝑗 𝑝𝑒𝑟𝑡𝑢𝑟𝑏 is set as maximum correspondent to 10% 𝑣𝑖𝑛𝑗. After 0.1 𝜏 no further perturbation is introduced on the domain and the velocity parabolic profile is set as equal in both fluids streams. Santos evaluate the influence of the perturbation on the flow dynamics and observed there is no influence of the initial perturbation after a transitory state, where the system evolves into chaotic or into segregated. From the results it was observed for all the Reynolds in study the transition region always lies prior than 1𝜏, therefore all the numeric results presented is was perform form for 𝜏=2 and 𝜏=3. Figure 3.3 - Perturbation time function on the right injector affecting the normal velocity to the injector's cross section Courant number emerges of solving partial differential equations using finite differences as approximation. In CFD time dependent problem the advance is time is major dependent of
3 CFD Model and Simulation of CIJ 56 the advective term, which dominates the solution. Therefore is necessary a progressed time less than a certain time for explicit time-marching method. Courant-Friedrichs-Lewy (CFL) is defined as 𝐶=𝑢 ∆𝑡 ∆𝑥≤𝐶𝑚𝑎𝑥 (3.21) Where 𝑢 in the magnitude of the velocity, correspondent in this case of 𝑢𝑚𝑎𝑥 of the parabolic profile, ∆𝑡 the incremental time step and ∆𝑡 the length of the interval, in this case correspond to the mesh size. If explicit (time-marching) solution is used then 𝐶𝑚𝑎𝑥<1 to ensure numerical stability, while implicit method usually don’t required so lower values of 𝐶𝑚𝑎𝑥. To guaranty uniformity in the results, the time step marching is calculates based on 𝐶𝑚𝑎𝑥=0.5 for each Reynolds in study 3.2.5 Discretization and Numerical methods The numerical implementation of the mathematical techniques used for temporal and spatial discretization of the flow equations along with the gradients discretization into the discrete domain is presented in Appendix A. 3.2.5.1 Steady and Unsteady Numerical Methods In a steady state flow regime the pressure-velocity coupling is used SIMPLEC since low gradient pressure field exists in the mixing chamber resulting in fast convergence and accurate results particularly in low-Reynolds flows. For the momentum discretization a third-Order MUSCL scheme is used because the flow behaviour naturally exhibit discontinuities in the flow field, particularly in areas of interest, therefore the using of a high accurate numerical scheme is required. The areas of interest also required higher interpolation of the pressure and therefore a second order discretization is used. In the end the gradients are evaluated using Least Squares Cell based scheme since is the best suited method to be used in structure polyhedral grids. When transient state is used it was chosen to have a second order implicit formulation to avoid temporal discontinuities. 3.2.6 Computational performance The DNS simulation was performed in a HPC environment using Dell PowerEdge R420 chassis connected with 40 Gbpds InfiniBand network, each node have 2 Intel E5-2450 CPUs, 68GB of DDR3 RAM clocked at 1600MHz and SATA HDD disk with 500GB 7.2K. For each Reynolds in study two full nodes was used where the computation cost of each one present in Table 3.1
3 CFD Model and Simulation of CIJ 57 Table 3.1 - Computational cost of each simulation performed in HPC environment 3.3 Results of Numerical Simulation In this section will be presented and discuss the results provided by the numerical simulation for Reynolds 150, 200, 250 and 300. The choice of the Reynolds lies on previous numerical and experimental results, presented in chapter 2, which clear conclude the mixing reach their maximum efficiency at Reynolds 300. Since the flow is driven by viscosity over inertia, the shear force of fluids laminae will stretch the particle of fluid into even smaller thickness allowing diffusion and reaction in the mould enter in place and create a uniform plastic component. The aim of this chapter is then identify the hydrodynamic scales which leads to the striation of the particles of fluid, provided emphasis in the location where the hydrodynamic scales are created and dissipated allowing in the future a tighter and efficient control of the process with the lower production cost possible. Additionally this works present a new alternative to study mixing in a laminar chaotic environment leading to a new world of possibilities in field of mixing. 3.3.1 Flow visualization Results from transient simulation for each Reynolds in study is presented in terms of velocity and vorticity maps to allowed better understanding of the flow dynamic behaviour. As introduced early the chaotic flow motion is initialized using a perturbation, increasing the velocity in one injector in a specific time. Previous study Santos (2003) prove the flow is independent form this initial perturbation after one residence time. Therefore for all the simulations performed will be study ranging from 𝜏=1 to 𝜏=3. For better visualization of the flow, animations is provided in a DVD, were is observable the velocity maps, 2D vorticity maps and 3D vorticity maps, allowing a better visualization of the impinging jets iteration with the direction of the flow rotation 3.3.1.1 Instantaneous flow field velocity maps From previous works it is know the flow from each injectors feeding streams impinging impinge align with the injector axis where the flow rapidly decelerate spreading radial the flow in a form of ellipsoid pancake. Oscillation of this impinging surface, observe upstream Reynolds Hours CPU h 150 521 16675 200 595 19057 250 670 21440 300 744 23822
3 CFD Model and Simulation of CIJ 64 3.3.1.2 Instantaneous flow field Vorticity map In this section will be presented 3D vorticity maps for each Reynolds in study. The movies correspond to this figures are presented in (FOLDER). Each figure represent a representation of the impinging surface by setting a constant value of vorticity enough to allowed clear visualization of the oscillatory motion patterns. Additionally in this figure is presented the vorticity normal to the front plane calculated as 𝜔 =∇×𝑣 =(𝜕 𝜕𝑥,𝜕 𝜕𝑦,𝜕 𝜕𝑧)×(𝑣𝑥,𝑣𝑦,0)=(𝜕𝑣𝑦 𝜕𝑥−𝜕𝑣𝑥 𝜕𝑦)𝑧 (3.23) The two dimensional voticity is defined as the cross product of the gradient ∇ with the flow field 𝑣 . It represent the rotation of the flow parallel to the 𝑧 direction. When coupled with the study of the strain rate It is a value tool to access phenomena of stretching of fluid particles since represents the areas where increasing both interfacial area (area between two fluid flow streams for mixing) and local gradients in fluid proprieties (driving potential for mixing across this interface). For this case the impingement surface is created with two colours. Orange is associated with positives values of vorticity while the blue colour represents negative vorticity values. By observing the normal impingement surface oscillation is observed when the angle of this surface go in clockwise then the region downstream the impingement point is filled with positive values and otherwise, when the impingement surface go anticlockwise then the same region is filled with negative values. The oscillatory behaviour of this surface will then create regions with a certain thickness and different rotations which will improving the stretching of the fluid particles. The thickness of this regions and the intensity of this rotation will have a huge influence in the mixing performance, and it was observed the thickness reduces with the increasing of the Reynolds number as opposite the intensity of the rotation will increase with the increasing of the Reynolds. It was observed some disconnections of small structures of high vorticities from the impingement surface with this small structures moving downstream to the chamber outlet. Visualization studies suggests with the increasing of the Reynolds number will increase the number of vorticity regions loosen and their intensity and size. The loosen of this regions will then transform the flow behaviour along the chamber.
3 CFD Model and Simulation of CIJ 65 By observing the oscillatory behaviour the impingement surface it is clear the high angle of the vorticity will cause the breaking of the jets movement. This can be cause since there is two jets with different direction but the same rotation direction which can cause an intensification of the jets bending by repelling each other. When occur then a large number of smaller structures with high intensity vorticity values is created in the impingement region. The small regions with high intensity but with different directions will couple together and create complex flow patterns. This small structures will then flow downstream of the chamber, significantly impacting the mixing performance. This small structures when loosen from the main structure and hit the jets can cause a local deformation of the jets with severe impact in the impingement point position and/or jets angle and oscillation. If breaking occur can cause the rupture of the fluids lamina Fonte (2012) which can cause the stop of the stretching of the fluids particle. Additionally the smalls regions with high intensity will travel along the chamber and perturb the flow behaviour. It is know this small regions usually carry along their path zones with very concentrate fluids phase, which in terms of mixing will reduced the performance. Alongside is it unknown the effect of this small vortex regions when travel along the chambers, but it is know this small structures are the main feature in reactors type (4 inlets), in terms of mixing performance, where the main feature is the confining the mixing in the top of the reactor while the phases are fully dissipated by small and high intensity vorticities regions It is important focus all of this visualization studies were performed for the same fluid with same momentum in each side. Different types of fluids in each side, especially if the viscosity and density differs allied with different stoichiometric ratios will produce a far more complex flow pattern which need be carefully study and could differ from those study in this section
3 CFD Model and Simulation of CIJ 66 Figure 3.16 - Vorticity map for Reynolds 150 at time 𝒕/𝝉=𝟏 Figure 3.17 - Vorticity map for Reynolds 150 at time 𝒕/𝝉=𝟐 Figure 3.18 - Vorticity map for Reynolds 150 at time 𝒕/𝝉=𝟑
3 CFD Model and Simulation of CIJ 67 Figure 3.19 - Vorticity map for Reynolds 200 at time 𝒕/𝝉=𝟏 Figure 3.20 - Vorticity map for Reynolds 200 at time 𝒕/𝝉=𝟐 Figure 3.21 - Vorticity map for Reynolds 200 at time 𝒕/𝝉=𝟑
3 CFD Model and Simulation of CIJ 68 Figure 3.22 - Vorticity map for Reynolds 250 at time 𝒕/𝝉=𝟏 Figure 3.23 – Vorticity map for Reynolds 250 at time 𝒕/𝝉=𝟐 Figure 3.24 - Vorticity map for Reynolds 250 at time 𝒕/𝝉=𝟑
3 CFD Model and Simulation of CIJ 69 Figure 3.25 - Vorticity map for Reynolds 300 at time 𝒕/𝝉=𝟏 Figure 3.26 - Vorticity map for Reynolds 300 at time 𝒕/𝝉=𝟐 Figure 3.27 - Vorticity map for Reynolds 300 at time 𝒕/𝝉=𝟑
3 CFD Model and Simulation of CIJ 70 3.3.2 Spatial evolution of the flow In the present case the flow is mainly driven by viscosity forces rather than inertial forces, which indicates the flow are laminar. The non-linearity of the NS equation produces an uncertainty flow behaviour, which can cause the chaotic nature or turbulent nature. Chaotic is defined by a set of non-linear equation which with a slightly change in the initial condition will greatly affect the flow behaviour in the next time steep. A succession of sight differences in each time step would greatly differ the results in the end of the simulation. Exist in the literature tree main approaches to measure chaotic systems (Argyris, Faust, & Haase, 1993; Froehling, Crutchfield, Farmer, Packard, & Shaw, 1981; Lester, Smith, Metcalfe, & Rudman; Y. C. Li, 2013; Narasimha, 1986; Ngan & Vanneste, 2011; Research, 2003; Sommerer, Ott, & Tél, 1997; Southerland, Frederiksen, Dahm, & Dowling, 1994; Stroock et al., 2002; Voth, Saint, Dobler, & Gollub, 2003; Wong, 2010): Stochastic approach – This is a probabilistic method to characterize to the behaviour of the flow field. Exist two main methods. The first is the Lyapunov exponents which quantifies the rate of separation of an infinitesimal close trajectories. The second is the correlation length and time scales which measure the dimensionality of the space and time occupied by a set of random points. Roughly speaking the Lyapunov exponents is the same as the inverse of the correlation length and time. This methods are generic and don’t provide relevant informations. Chaotic advection approach – This method focus in specific features of the flow. The main methods used to calculate are spatial distributions of Lyapunov exponents and stable and unstable manifolds. Those methods are incredible hard to implement and don’t provide such important quantification tools and this is why are often used as diagnostic tool. Eigenfunction approach – This method solve infinite-dimensional eigenvalue problem for advection-diffusion problems, by calculate the Eigenmodes of the data series and only focus in the dominant Eigenmodes. This method is very hard to implement and focus in special characteristics of the problem. All of the methods used to quantify and qualify the results from chaotic systems don’t provided enough practical information to access to main problem, thus is required to use other methods to study spatially to access the flow hydrodynamic temporal and spatial scales. Turbulent methods in other hand ensure quality of the result this relevant and practical information therefore will be used in the present work.
3 CFD Model and Simulation of CIJ 71 Along the domain a set of point, align with the axis of the chamber is positioned to record each time step, the velocity components of the flow in each point. The points was chosen based on the following equation. Data series provide the simulation is then treated where is only consider 1 to 3 𝜏 data, since the first 𝜏 is used to create the oscillatory behaviour. 0,𝑑,𝑑+2−6𝐷,...,𝑑+22𝐷 (3.24) In order to quantify the flow, turbulent methods of investigating the flow behaviour was used. Initially is assumed the flow can be characterize by their mean velocity and the fluctuations 𝑣∗=𝑣+𝑣′ (3.25) Where 𝑣 is the average velocity, 𝑣(𝑡) is the instantaneous velocity and 𝑣∗ is the velocity which contain the velocity variation along the time series. The oscillatory time dependent velocity is shown in Appendix B. It is observed in this appending an irregular behaviour with different magnitude and frequency passing through the point monitor. Such region is referred to as an eddy or turbulent structure which simply pass by this monitor. The size of the region and its representative velocity are namely the length scale and the velocity scale. The Appendix B contain the time series used in the analysis of this problem. It is observable a vanishing of the velocity magnitude and a slightly increase of the length scale as the point is progress towards the chamber outlet. Figure 3.28 to Figure 3.31 show the each velocity component against the position of the point normalized by the chamber diameter. Being 𝑣𝑦 the velocity of the injector and 𝑣𝑥 the velocity parallel with the chamber axis and 𝑣𝑧 the velocity normal to boot of them. As expected the 𝑣𝑦 take in the point zero, the velocity close to the injector axis velocity, with fluctuations corresponding to the side of the impingement point position with a fast vanishing of those fluctuations in the downstream points. The impingement will then convert velocity in y direction to velocity in z direction thus is observable by the bigger z velocity in the impingement point. According to the x velocity is observable a higher value in point 1 due the spreading of material from the impingement point. Generally speaking is possible to conclude the following. Fluctuations in y-velocity is mainly driven by the oscillation in the impingement point position and his velocity almost don’t lose kinetic energy from the injection to the chamber to impingement. The magnitude of the xvelocity when compared with z-velocity, prove the spreading of fluid is mainly parallel with the chamber axis rather than opposite direction.
3 CFD Model and Simulation of CIJ 72 Figure 3.28 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 150 Figure 3.29 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 200
3 CFD Model and Simulation of CIJ 73 Figure 3.30 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 250 Figure 3.31 - Mean and fluctuations velocities components toward the chamber outlet, for Reynolds 300
3 CFD Model and Simulation of CIJ 80 Figure 3.38 - Impingement point position, Z-velocity along the injector axis The maximum edge size for all the simulation are 100 𝜇𝑚. Consider all the values bigger than the resolution of the mesh. It is concluded, face of the error in resolution and averaging, the impingement point position, in other words, the location where the velocity change in direction is in [0, 0, 0] m, which means the impingement point position are centred with the origin. This prove will be used when was develop a deterministic algorithm to obtain the stagnation point position and the jets angle. Important fact to refer is the highest rate of change in the velocity, in order words, the higher partial derivative, is in all representations around the stagnation point position, being for higher Reynolds number a slight increase in this value. Higher derivatives are in y and x direction, namely the injector axis direction and the chamber axis direction. This conclusion will later held to understand the rate of dissipation of turbulence 3.3.3.2 Velocity The average velocity component is show in Figure 3.39 to Figure 3.42 for each Reynolds in study. In the left side is shown the top plane contour of the velocity and on the right side is show the velocity along a line aligned with the injector axis and which cross the points used to record the velocity (Equation 3.24).
3 CFD Model and Simulation of CIJ 81 For all the Reynolds in study is observable the same behaviour: In the plot correspondent to the x velocity the point 0 have an inverse direction of the velocity, since the position of the stagnation point, obtained in the simulation is around 50 𝜇𝑚 from the origin tending downstream. Although there is not resolution to confirm this value numerically this produce a shift in the sign. The adjacent point p1 have the maximum x-velocity since all the volume contained energy provided by the injectors are directionally downstream. The following point have consequently a reduced velocity since there is a spreading of material in more than one direction thus is clearly observable an cone with high velocity, being this cone formatted due the oscillatory behaviour of the impingement surface. Downstream of the jets a negative sigh of the x-velocity is observable since it correspond to a low speed recirculation zone. Along the chamber a vanishing of the high central located velocity into a more spread distribution, being in the end almost uniform x-velocity in the chamber outlet meaning in the total height of the chamber a vanishing of high velocity and high intensity fluctuations into a uniformly distributed velocity profile translating no more mixing advection mechanics close to the outlet of the chamber. The y-velocity profile presents in the injector axis, p0, the maximum velocity magnitude. It is clear from this graphic the reduction of the velocity since enter in the chamber don’t lose its kinetic energy, meaning, the spreading of the fluid due impingement is mainly influence by the injector’s velocity. In impingement the velocity quickly decays in very low distance, producing a fast deceleration of the fluid which is the mainly cause of the spreading and oscillatory motion of the impingement surface. After the impingement point and progress towards the chamber outlet there is low value of velocity, slightly vanishing until the chamber outlet. The z-velocity in the Reynolds in study present to have a random behaviour without any common reference. Generally speaking the magnitude of this velocity don’t present a relevant factor in the advection mixing since the value is low. Is observable an homogeneous and uniform velocity in the chamber outlet for all the Reynolds. The low value and random pattern could be due averaging, since velocity in this component is mainly driven by the location of the impingement point and the oscillatory motion of the impingement surface. All the comments was perform for all the Reynolds in study in a general sense. An important note, which was referred previously the averaging of the velocity component was necessary to compute the rate of creation and dissipation of turbulent kinetic
3 CFD Model and Simulation of CIJ 82 energy, thus the stabilization value of this variables in the whole was consider to limiting number of time steeps for averaging due space limitations. The values and contours present was averaged from the instantaneous velocity field considering this limiting factor, which could cause a lack of resolution in some locations or, in case of low values and shifting velocity directions each time step.
3 CFD Model and Simulation of CIJ 83 A. Reynolds 150 Figure 3.39 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 150
3 CFD Model and Simulation of CIJ 84 B. Reynolds 200 Figure 3.40 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 200
3 CFD Model and Simulation of CIJ 85 C. Reynolds 250 Figure 3.41 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 250
3 CFD Model and Simulation of CIJ 86 D. Reynolds 300 Figure 3.42 - Contour of the velocity components in the front plane (left). Velocity components along the probe points in a line aligned with the injectors axis (right). Reynolds 300
3 CFD Model and Simulation of CIJ 87 3.3.3.3 Rate of creation and dissipation of turbulent kinetic energy The rate of creation and dissipation or turbulence was calculated based on the following scheme: Figure 3.43 - Algorithm to obtain fluid flow variables 𝒌 and 𝜺 Initially was read a specific number of data files, each one correspondent to a certain time step. The choice of the time step was used as quasy-spaced time step along 1𝜏 to 3𝜏 of data files, with a random choice of the time step according the mutation rate, in other words was chosen data files correspondent even distributions which after wards will be randomly Read data files Create UDS with 𝑣𝑥,𝑣𝑦𝑣𝑧 Write UDS Read 𝑣𝑥,𝑣𝑦𝑣𝑧 Averaging 𝑣𝑥,𝑣𝑦𝑣𝑧 Write 𝑣𝑥 ,𝑣𝑦 ,𝑣𝑧 Read 𝑣𝑥 ,𝑣𝑦 ,𝑣𝑧 as UDS Using UDF calculate 𝑠,𝑘,𝜀 as UDS Write UDS Read 𝜀,𝑠,𝑘 Averaging 𝜀,𝑠,𝑘 Write 𝜀,𝑠,𝑘 Read 𝜀,𝑠,𝑘 as UDS Post processing UDS
3 CFD Model and Simulation of CIJ 88 chosen data file in a certain range (more and less) the specific time step. This method help to prevent stroboscopic effects of the fluid in the results, particularly in the averaging, if it happens. Additionally after all the process being repeated help to analyse the stability of the method and the accuracy of the result. The previous scheme is set with colours to represent each step used, being the orange colour associated with fluent process and grey as of the side created code in c to run in windows OS. First using fluent the selected time steps are read and the velocities are exported using a UDS (used-defined scalar). After a set number of time steps are read the UDS’s are written in a text file. A code is created in c language to read the text file, containing all the time steps for each variable and calculate the averaging of each cell in the domain. The averaged values are then written into a text file to be used in the third phase. A UDF (user-defined function) is created to calculate the rate of creation of turbulent kinetic energy and the spatial derivatives which is necessary to calculate the rate of dissipation of turbulence. Using the same process all of this this variables are averaged. In the end all the variables are read and set as UDS and all the essential post-processing is perform. Since exist close to three million cells in the domain and each one individual cell and their correspondent variable are save in the hard drive, an increasing of the time steps used in the calculation will greatly affect the available space needed for the calculations. Since this physical constrain then was necessary to proceed to a study of the strict necessary number of time steps used in the calculations. Because of this restriction was perform a test to evaluate the essential number of time steps Figure 3.44 -Convergence due averaging time steps for 𝒌 and 𝜺
3 CFD Model and Simulation of CIJ 89 In Figure 3.44 is presented the convergence history for 𝜀,𝑘. Both of them were calculated based the volume integral of this scalar in the whole domain. Those convergence histories is presented for one set of averaging, spite the fact that more than one averaging was perform and the results seem to agree with the current presentation. Due lack of space it was impossible run for more than 300 averaged results. It is observed for an averaged time steps of 100 and 150 the values seem to have a big difference from averaged values. Although this is true the relative difference of the error is not so high. For a trustful results was chosen that an average of 300 time steps guaranties in both high number of time steps to averaged, low relative error from the previous averages and with reasonable space to perform the calculations. Recall the physical meaning of the creation and dissipation of turbulence. 𝑘, creaction of turbulent kinetic energy is the variable which count the fluctuation of velocity of the flow field, meaning an local increase of those fluctuations will then concentrate the production of turbulence in that area. In opposite 𝜀 stands for the dissipation of turbulence, which is calculated by the partial derivatives of the flow field, meaning high local derivatives will then create a local area where occur turbulent dissipation. The results shows contours of both variables and also a plot of those variables in lines parallel with the injector axis along the chamber crossing the points chosen to record the velocity each time step. The importance of this work lies in identification of the zones exhibiting local production and dissipation of turbulence which physically represent the creation of vortices and their transformation into small vortices by turbulent cascade, which are the main phenomena of advection mixing in turbulent flows. This technique is widely used in full turbulent flow systems although for the first time will be applied into a laminar chaotic internal flow. In some publish works some authors prefer to use turbulence modelling for modelling a laminar chaotic system (Icardi et al., 2011). Those modelling are by using RANS, LES to model the flow field. RANS consider a time space averaging without scales resolving, as opposite LES uses a sub-grid to compute scales up to a certain level, and therefore only compute larger and medium scales. It is true both techniques perform better than DNS, when only computational cost is considered, although they are not considered the quality of the results. All the results presented with this techniques only show the effect of the larges scales modelling of the problem to identify factors such as impingement point position and local velocities. But when medium and lower scales are modelled the author don’t agree with the using of such tools. In part it is know that turbulence model are statistically created for
3 CFD Model and Simulation of CIJ 96 Figure 3.50 - Variation of 𝒌 along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of 𝒌 for the top part of the chamber (bottom). Reynolds 250
3 CFD Model and Simulation of CIJ 97 D. Reynolds 300 Figure 3.51 - Variation of 𝜺 along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of 𝜺 for the top part of the chamber (bottom). Reynolds 300
3 CFD Model and Simulation of CIJ 98 Figure 3.52 - Variation of 𝒌 along the chamber for lines parallel with the injector axis and passing in the probe point (top), Contour of 𝒌 for the top part of the chamber (bottom). Reynolds 300
3 CFD Model and Simulation of CIJ 99 For all the Reynolds number in study it is clear the results tends to be symmetric regarding the centre of the chamber. Non perfect symmetry is only due lack of more averaging time steps. Spite this fact it is only present the important variables for a averaging of 300 time steps, others averaging such as 50, 100 and 150 was visualized, and qualitatively have the same shape although with slightly different values. In the present work don’t focus in the value itself but in the shape of interesting variables. According the creation of turbulent energy it is clear that a cone without turbulent emerges. This is due the jet´s don’t carry turbulent energy and its high kinetic energy don’t allowed any creation in this region. Is the perimeter of the jets with low velocity which in contact with the flow create this conical structure. All the turbulent energy is created in the impingement point and the dissipated trough the pathline of the flow due the oscillatory mechanics of the impingement surface (translation of the impingement point and rotation of the surface). From the contour of this variable it is clear the main creation of turbulent energy in the front plane fulfil in height approximately six times the injector diameter and in width 4 times the injector diameter. In the top plane similar relation emerge, begin the height almost fulfil the entire chamber diameter and the width have the same relation, 4 times the injector diameter. This relations are in the same range for all the Reynolds in study. It is observed for Reynolds 150 the production of turbulent energy have a large planar structure for production of such energy. An increasing of the Reynolds number will then tend to reduce the size of such structure into a format of a boll. It is observable than for high Reynolds number this high energy region will tend to locate the production into a single point, where opposite for low Reynolds number, the pancake form will help to distribute the region where this energy is created almost fulfil the top of the entire chamber. The dissipation of turbulent energy create some structures which go from the high intensity dissipation into the camber walls by following the path of the low velocity of jet’s, particularly in the top of the jets. Recalling the mathematical expression of this variable, it is directly correlated by the high velocity gradient. Therefore is normal to have high dissipation where each component of the velocity fast changes, which is normal close the jet’s and in the impingement point. Regarding about the location where the energy is dissipated it is clear almost the same relationship between height and width for the top and front surface. Even that it is observable more complexes patterns than 𝑘.
3 CFD Model and Simulation of CIJ 100 A careful attention should be made when is comparing the contour for 𝑘 and 𝜀. For Reynolds 150 the contour of 𝜀 is cropped to the half oh the maximum value while for 𝑘 the contour almost present all the values. For Reynolds 200 booth contours are cropped for almost 2/3 of the maximum value for each variable. For the other Reynolds the contour of 𝜀 is almost cropped to ½ of the maximum value while 𝑘 is only 1/3. To better analyse this effect a comparison between 𝑘,𝜀 is presented by comparing the difference between the maximum value for the point 0 with the maximum value of the point 5, which closely correspond the end of the impingement surface for the average Reynolds in study. Defining ∆𝜀=𝜀𝑝0−𝜀𝑝5 (3.34) Similar for 𝑘 ∆𝑘=𝑘𝑝0−𝑘𝑝5 (3.35) And ∆𝑝/𝑝0=∆𝜀/𝜀𝑝0 (3.36) ∆𝑝/𝑝0=∆𝑘/𝑘𝑝0 (3.37) Based on the previous definitions is possible to identify the behaviour of those variables 𝜀,𝑘 by representing into a plot which is in Figure 3.53. Figure 3.53 - Variation of 𝒌 and 𝜺 along the numerical experiments. Observation energy containing inside the impingement structure (left). Normalization of this energy with energy containing in the central area (right)
3 CFD Model and Simulation of CIJ 101 According to Figure 3.53 it is possible to observe a constant growing by an increase of the Reynolds number. What is surprising and was observed in k obtaining from the points is a almost stagnation of the flow behaviour from 200 to 250 in both 𝑘 and 𝜀. It is unknown what mechanism produce such stabilization of the flow behaviour but this means in terms of both variables an increasing of the Reynolds number don’t produce such advantage in the increasing of the hydrodynamic mixing in the range 200 to 250. An increasing is observable in the range 150 to 200 and from 250 to 300, being this range where the main growing is visualized with more than 250% of increasing of the value for each variable in study. About the Figure 3.53 if it obtain the average between 𝑘 and 𝜀 for the increasing Reynolds number is observable a small increase of both relationship variables when the Reynolds number is increased. After Reynolds 200 is observed a fluctuation of this relationship for both 𝑘 and 𝜀 with a small difference between each one. For Reynolds 150 the difference is much higher, where the area for the dissipation overhauls the area of the production of turbulent energy Knowing the length scale relationship for isotropic and homogenous turbulence is defined as 𝑙~𝑘23 𝜀 (3.38) Then for an increasing in Reynolds is observable a decreasing of the length scale for the range 200 to 300, while in observable an increasing of the length scale from 150 to 200. This effect indicate than Reynolds 250 makes the transition from an incremental length scale to a reduction of length scale which allied with a fast growing rate of those variable have as consequence the reduction of the eddies sizes and the amount of eddies created which facilitate the transfer of energy by energy cascade and the mixing overall performance.
3 CFD Model and Simulation of CIJ 102 3.3.4 Spectral analysis Time series provided by the sampling points usually don’t provide direct information about the flow behaviour, then is necessary to use other tools to quantify the flow. As previous turbulent methods was used to quantify the flow field using spectral analysis methods to transform the data points into the energy-containing-frequency domain. In Appendix C is presented all the power spectra for all the points and velocity components for each Reynolds in study. Small conclusions was presented in Appendix B was the influence of the Reynolds number in the flow regime, the non-random velocity oscillation, an increase of the frequency of such oscillations with the increasing of the Reynolds number and the vanishing of such velocities and oscillation when the flow go downstream the mixing chamber. The transformation from time into frequency domain was done using the Discrete Fourier transform DFT, defined as 𝐺=(𝜙= 𝑘 𝑁∆𝑡)=∑𝑣𝑦(𝑖∆𝑡)𝑒−𝑗 2 𝜋 𝑖𝑘/𝑁, 𝑁−1 𝑖=1 𝑘=0,1,2,..,𝑛−1 (3.39) Where 𝑣𝑦(𝑖∆𝑡) correspond to the velocity time series, ∆𝑡 is the time interval, 𝜙 is the frequency, 𝑁 is the total number of points existent in the time series and 𝑒−𝑗 2 𝜋 𝑖𝑘/𝑁 is the complex number where 𝑗 is the imaginary part. Fast Fourrier Transformers, FFT, substitutes DFT, in terms of implementation since improve the computational effort to produce such calculations. FFT provides the complex number series, which to translate into power spectrum is necessary to calculate the modulus |𝐺(𝜙)|.The power spectra presented in this chapter and in the appendix (xxx) is normalized by the outlet velocity for each simulation performed. The outlet velocity is calculated based on the conservation off mass in the control volume 𝑚𝑖𝑛𝑗1+𝑚𝑖𝑛𝑗1=𝑚𝑜𝑢𝑡 (3.40) Where 𝑖𝑛𝑗 represent in injector and the indices 1 and 2 represent the left and right injector. Since the is an incompressible fluid with the same characteristics in each boundary and each injector have the same velocity then the previous equation can be simplified into 2𝑄𝑖𝑛𝑗=𝑄𝑜𝑢𝑡 (3.41)
3 CFD Model and Simulation of CIJ 103 2𝜋𝑑2 4∙𝑣𝑖𝑛𝑗=𝜋𝐷2 4∙𝑣𝑜𝑢𝑡 (3.42) 𝑑2 𝐷2∙𝑣𝑖𝑛𝑗=𝑣𝑜𝑢𝑡 (3.43) Yielding the normalization of the frequency series to obtaining the power spectrum 𝐺′(𝜙) is performed by 𝐺′(𝜙)=|𝐺(𝜙)| 𝑣𝑜𝑢𝑡 (3.44) 3.3.4.1 Power spectra Using MATLAB an automated code was created to read the data sampling and compute FFT and in the end plotting the respective data. From this results will be investigated the results provided by the point 0 and point 1, correspondent to the point align with the injector axis and the first point after this following downstream the chamber. Although all the results is presented in the Appendix C, only go to be presented this two points for they-velocity component just to latter being possible to compared with experimental data. By turbulent theory it is know the power spectra have three range of scales: Energy containing scales: is where most of the energy is containing into large length scales associated to large eddies Inertial range: is where the energy is transfer. In this range two main slopes can emerge if the case of an injection of energy is observed. -5/3 slope is where an inversion of the energy cascade, where smaller scales feed up the larger scales. -3 slope represent the direct entrosphy cascade where the direction of the scales is positive, meaning the larger scales are dissipated into smaller ones and the process continue until the dissipation range. Is the injection where the shift between these two slopes emerge Dissipation range is where the scales reach the “final destination” and are transform into heat This knowledge found the ground base to access the problem in hand, where will be discussed
3 CFD Model and Simulation of CIJ 104 Figure 3.54Power spectra form all the Reynolds number performed [150, 200, 250, 300] for the probe 𝒑𝟎
3 CFD Model and Simulation of CIJ 105 Figure 3.55Power spectra form all the Reynolds number performed [150, 200, 250, 300] for the probe 𝒑𝟏
3 CFD Model and Simulation of CIJ 112 Figure 3.56 - Phase maps for Reynolds 150 at time 𝒕/𝝉=𝟏 Figure 3.57 - Phase maps for Reynolds 150 at time 𝒕/𝝉=𝟐 Figure 3.58 - Phase maps for Reynolds 150 at time 𝒕/𝝉=𝟑
3 CFD Model and Simulation of CIJ 113 3.5 Conclusions From the hydrodynamic behaviour is concluded the following DNS simulation was perform for Reynolds number 150, 200, 250, 300 since it is necessary to evaluate the hydrodynamic behaviour to create and develop new technologies to improve the mixing dynamics. The breaking of the oscillatory movement of the impingement surface due jets bending or high angles, cause complex behaviour in the impingement region with biggest changes in the flow field, particularly in the vorticity. The impingement point and a impingement surface have a 3D motion but the main variations are in the front plane. Turbulent quantitative methods was successful applied to study laminar chaotic flows with emphases in the creation and dissipation of turbulent energy and correlation with the geometrical location of such variables and the relation between both for each Reynolds in study. Additionally was successful study the correlation between the eddies size with the geometrical characteristics of the chamber. The author recommends the development of new technologies for low Reynolds number since it reduces the equipment costs due high pressures and tight control of the flow dynamic. New works tend to progress in study the flow using mechanisms to induce oscillation in the jets. The author believe the usage of such technology in low Reynolds number will greatly affect positively the flow dynamics and consequently improve mixing performance. From the mixing dynamics it is concluded the following: For the first time was perform a representation of the mixing dynamics in a 3D environment using VOF for interface tracking and without molecular diffusion. Mesh resolution was not enough to calculate the smaller mixing scales and therefore the results appear to have a big numerical diffusion thus only was possible to verify the larger vortices. The author recommends in the future, when computer make accessible a creation of the same method VOF for the 3D domain but with interface taking based on the gradients of the phases (Fonte, 2012) Note: this was test but due the large necessity to refinement is practice impossible nowadays, with the resources available to tackle this problem.
114 4 Development of a Real-time Control System for the CIJ Flow Field with PIV 4.1 Introduction Over the past years, researchers have used the flow visualization technique based on the Particle Image Velocimetry (PIV) to study the flow field behaviour in several industrial and academic applications and the same holds true for studies in Reaction Injection Moulding. The PIV technique was used in the present work to compute the 2D velocity field from the displacement in 2 consecutive images of a small point which flows with the current. The method will be explained later in section 4.4 and the essential hardware and software employed will be described in subsection 4.4.1. This method has proved to be extremely useful to compare the experimental results with the CFD data. Although the hydrodynamics of the CIJ flow field was already investigated in previous works (Santos, 2003; Fonte, 2012; Gomes, 2015), it was necessary to replicate some parts of the work performed in order to ensure the quality of the results. Based upon the results of this hydrodynamic study, it was possible to analyse and validate other tools that can be used to establish operational and design parameters of the CIJ mixer. First, the elastic model proposed by (Bird et al., 2002; white, 2006) and further developed by (Fonte, 2012)was considered. To allow validation of the model, when possible similar conditions to those of (Fonte, 2012) work were employed. In addition, a Reynolds number dependence study was
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 115 performed. Next, the pressure model proposed by (Gomes, 2015) was examined and, for the first time, it was validated against experimental data. This chapter is structured as follows. Section 1.2 defines the hydrodynamic variables that are relevant to describe the operating conditions in the system. Section 1.3 describes the main characteristics of the experimental facility. Section 1.4 explains the PIV technique and includes a description of both the hardware and the software used. Section 1.5 concerns the study of the turbulence intensity within the mixing chamber. Section 1.6 gives an overview of the elastic model, presenting its validation and an extension of the model for different Reynolds numbers. Section 1.7 introduces the pressure model and discusses its ability to reproduce the experimental measurements. Section 1.8 addresses the main conclusions that can be drawn from this chapter. 4.2 Hydrodynamic variables The hydrodynamic variables introduced (Malguarnera and Suh, 1977; Macosko, 1989) to carry out the dimensional analysis of impingement mixing in RIM machines will be adopted hereafter to describe the designated operating conditions. Initially, the Reynolds number is used to describe the flow motion. In impinging jets mixers, the Reynolds number is defined at the injectors as 𝑅𝑒=𝜌𝑣𝑖𝑛𝑗𝑑 𝜇 (4.1) where 𝜌 is the fluid density, 𝑣𝑖𝑛𝑗 is the velocity at the injectors, 𝑑 is the injectors diameter and 𝜇 is the fluid viscosity. The fluid properties are highly dependent on the operating temperature, therefore, in addition to a careful data acquisition during the experiment, a good knowledge of the fluid rheology is necessary. Additionally, several parameters can be employed to express the jets’ flow unbalancing in impinging jets mixers (Macosko, 1989). The jets’ momentum rate ratio, 𝜙𝑀, is defined as 𝜙𝑀=𝜌1𝑑12𝑢𝑖𝑛𝑗,1 2 𝜌2𝑑22𝑢𝑖𝑛𝑗,2 2 (4.2) Where the indices 1 and 2 represent the left and the right injector, respectively. The kinetic energy rate ratio between the two jets, 𝜙𝐾, is defined as
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 116 𝜙𝐾=𝜌1𝑑12𝑢𝑖𝑛𝑗,1 3 𝜌2𝑑22𝑢𝑖𝑛𝑗,2 3 (4.3) As for the jets’ mass flow rate ratio, it is defined as 𝜙𝐹𝑅=𝜌1𝑑12𝑢𝑖𝑛𝑗,1 𝜌2𝑑22𝑢𝑖𝑛𝑗,2 (4.4) All of these definitions will be mentioned later when the analysis of the flow regime is coupled with the identification of the impingement point by means of two new methods that were developed in previous works (Fonte, 2012; Santos, 2003). The first method to be considered is an elastic analogue model. Then, at the end of this chapter, a pressure model will be discussed. These two models are extremely important in order to perform real-time control of the system as their application makes it possible to select operational and geometrical conditions that set the impingement point at the centre of the chamber, promoting maximum oscillatory motion of the impingement surface and thus enhanced mixing. 4.3 Domain 4.3.1 Physical characterization of the mixing chamber To allow visualization of the flow behaviour, the mixing chamber was machined from a Plexiglas block and the internal walls were polished to remove defects caused by the manufacturing process. Figure 4.1 shows a schematic image of the chamber. The injectors are perfectly aligned with each other and with the chamber axis. The chamber diameter, D, is 10 mm and the injectors diameter, d, is 1.5 mm. The injectors have a length of 6D to allow a fully developed parabolic flow and are positioned 5 mm below the top of the chamber. The chamber height, H, is 50 mm. The fluid enters in the chamber through the opposed injectors and is mixed by impingement prior to the disposal in the mold. The liquid is then discarded from the mold to the exterior and it will go on to be reused during the experiments. The fluid physical properties do not suffer significant variations during the recirculation since the tanks, the chamber, the mold and the storage vessel were clean and dry before the experiments took place.
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 117 Figure 4.1 – Schematic of the transparent mixing chamber and the mold 4.3.2 Schematic representation of the RIM machine The Reaction Injection Moulding (RIM) machine where the experiments were carried out comprises the equipment depicted in Figure 4.2. Two containers (1) feed the working fluid to each injector whereas the additional container (6) is used to collect the fluid discharged from the mould. A hydraulic loop feeds each container (1) with the working fluid that is collected in the disposal container (6). The two containers (1) are air pressurized in order to guarantee a smooth injection of the fluid without pulsation, which occurs if the fluid injection is done with pumps, causing pressure spikes. This ensures total control of the RIM process and allows making direct comparisons between the results of the PIV technique and those of the CFD simulation.
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 118 Figure 4.2 – Representation of the RIM equipment: a) hydraulic loop, b) transparent mixing chamber and mold 4.3.3 Physical characterization of the working fluid A tight control of the working fluid composition is required since the process variables, namely the Reynolds number and the jets’ kinetic energy rate ratio, which is a function of both the density and the viscosity, are highly dependent on the fluid properties. A further issue that arises is that those properties change drastically with the operating temperature. Therefore, it is necessary to determine the density and the viscosity of the fluid as a function of temperature in the designated operating conditions, which range from 20º to 30º. The liquid used in the experiments is a glycerin (99%) solution in water. The mass fraction of the glycerine solution, defined as 𝜒𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒=𝑀𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒 𝑀𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 (4.5) where 𝑀𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒 Is the mass of glycerol in the solution and 𝑀𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 is the total mass of the solution, was set at 73%. The density of the solution is determined by
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 119 𝜌=(1−𝜒𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒)𝜌𝑤𝑎𝑡𝑒𝑟+𝜒𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒𝜌𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒 (4.6) where 𝜌𝑤𝑎𝑡𝑒𝑟 and 𝜌𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒 (kg/m3) are the water and the glycerine densities, respectively. The variation of the water and the glycerine densities with temperature 𝜃(𝐾) is given by (Santos, 2003) 𝜌𝑤𝑎𝑡𝑒𝑟=1064.6−0.23𝜃 (4.7) 𝜌𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒=1440.5−0.61𝜃 (4.8) In addition, the viscosity of the solution is computed by the following expression 𝜇=((1−𝜒𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒)(𝐴𝑤𝑎𝑡𝑒𝑟𝑒𝐵𝑤𝑎𝑡𝑒𝑟 𝜃)+𝜒𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒(𝐴𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒𝑒𝐵𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒 𝜃)𝛼)1𝛼 (4.9) where 𝜇 is the viscosity in 𝑚𝑃𝑎.𝑠, 𝐴𝑤𝑎𝑡𝑒𝑟=5.17𝑥10−4 𝑚𝑃𝑎.𝑠, 𝐵𝑤𝑎𝑡𝑒𝑟=2.22𝑥103 𝐾, 𝐴𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒=7.06𝑥10−9 𝑚𝑃𝑎.𝑠, 𝐵𝑔𝑙𝑖𝑐𝑒𝑟𝑖𝑛𝑒=57.64𝑥103𝐾 and 𝛼=−0.32. This theoretical curve was compared with viscosity measurements made with a rheometer (Paar Physica UDS 200) with a plane (MK24) and plate in the temperature range of 20º to 30º. Figure 4.3 demonstrates the dependence of the working fluid viscosity with temperature. Figure 4.3 – Variation of the glycerine-water solution with temperature 0.01 0.012 0.014 0.016 0.018 0.02 18 20 22 24 26 28 30 32 34 Viscosity (Pa.s) Temperature (oC)
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 120 4.3.4 Recording hardware The flow is controlled with two needle valves, one in each injector line, and adjusted through two ROTA-Yokogama Coriolis mass flow meters (model RCCCS33 M01A1SH) with a precision of ±0.05𝑔/𝑠 in the range of 1.2g/s to 13.2 g/s. For the lowest flow rate (0.5 g/s) this corresponds to a relative error of 10%. Additionally, the mass flow meters enable the measurement of the fluid temperature and density for a single experiment, thus allowing a more precise control of the experimental variables. For each run, the pressure is recorded by a calibrated differential pressure transducer (Validyne P305D) fitted at the start of the injectors. The signal is processed through a data acquisition system (GW Instruments MacADIOS II) connected to an external board (GW Instruments MacADIOS ABO) with a computer using a program written in LABVIEW. The same recording system is employed for the mass flow readers. The data acquired is then processed in MATLAB to quantify the results obtained. 4.4 PIV technique In the late 90’s, developments in the acquisition of high-quality digital images allowed new techniques for flow visualization to emerge, one of them being Particle Image Velocimetry (PIV) (Westerweel, 2008; Westerweel, Elsinga, and Adrian ,2013). The PIV technique requires a sequence of snapshots of the flow field in order to correlate seeding points in consecutive images and create a displacement vector for the domain under study. These points are created experimentally by seeding the fluid with particles that reflect the light. By adjusting a beam planar laser it is possible to obtain the displacement field in any 2D plane in the mixing chamber. 3D PIV is also possible with the implementation of an additional camera, but was not performed in the present work. A charge-coupled device (CCD) camera with a high transfer rate allows taking shots at the plane illuminated by the beam laser. From two consecutive frames the displacement field is computed using a cross-correlation method that consists of a digital image processing technique to identify high gradients in the image created by the seeding particles. This is performed for each division of the target area on the chamber, known as interrogation spot. The capture time is chosen according to the size of the interrogation spot and the velocity field passed by. It should be assured that the particle moves less than a fourth of the interrogation spot size between two frames. Due to the low acquisition rate of the CCD cameras relative to the instantaneous velocity this method can only be applied when low
4 Development of a Real-time Control System for the CIJ Flow Field with PIV 121 velocity fields are expected. In order to overcome this limitation, a laser beam is pulsed close to the camera shutter is almost closing and the second beam is pulsed as soon as the camera shutter opens again. This reduces the time between the capture of two consecutive pictures drastically and therefore boosts the technique capabilities. At the experimental facilities available in the LSRE laboratory, a TSI PIV system is used to perform 2D analyses of the flow field. To carry out the PIV measurements it is necessary to employ a pulsed Yag laser, a CCD camera and a synchroniser, as well as a computer to control each component, acquire the images and post process them with in-build software to yield the instantaneous velocity vector maps. The image acquisition was done in an Intel Xeon CPU at 2.33 GHz with 2GB of RAM. Additionally, a PCI card was utilised to transfer the images from the camera to the computer prior to their analysis. A photograph and a schematic of the apparatus are shown, respectively, in Figure 4.4 and in Figure 4.5. Figure 4.4 – PIV experimental setup
5 Real time visualization of the flow field using PLIV 224 Starting by analysing the frequency recorded by the differential pressure transducer. For the case of segregated flow dynamic is observable some noise in the range of 60/70 Hz, which could indicate some electromagnetic interference which could be introduced by external equipment such as computer monitor or the high speed camera. This hypotheses is give more weight to validate since when the self-sustainable chaotic regime is obtained this value vanish. When the Reynolds is increase it would be expected the maximum frequency of the system increase as documented by Santos (2003) but this value tend to fluctuate far away from the expected frequencies. For the stagnation point position it is observable similar values such in the documented range until Reynolds 123. Far from this Reynolds this variable don’t present physical reliability to predict the main frequencies in the flow field. For the jets angle same analyses was conducted. It is observable a very good agreement was obtain until Reynolds 273 which indicates the driving force of the flow field is directly related with the jets angle rather than the impingement point position. For Reynolds above this value there is not observable any correlation with documented expectable frequency. Above Reynolds 273 and due the combination of small angle of the impingement jets, associated with succession of breaking of the impingement surface and fast displacement of the impingement point would not allowed record the frequencies existent in the impingement surface An inverse energy cascade, represented by the slope -5/3 is observable for all the variables in all the Reynolds in study, meaning it was capture the regime where injection of energy is obtain until the large energy containing scale, which is only due the strong dynamics of the impingement surface.
5 Real time visualization of the flow field using PLIV 225 Figure 5.51 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 26 Figure 5.52 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 65 Figure 5.53 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 83
5 Real time visualization of the flow field using PLIV 226 Figure 5.54 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 90 Figure 5.55 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 100 Figure 5.56 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 123
5 Real time visualization of the flow field using PLIV 227 Figure 5.57 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 134 Figure 5.58 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 171 Figure 5.59 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 205
5 Real time visualization of the flow field using PLIV 228 Figure 5.60 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 270 Figure 5.61 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 273 Figure 5.62 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 275
5 Real time visualization of the flow field using PLIV 229 Figure 5.63 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 344 Figure 5.64 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 421 Figure 5.65 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 568
5 Real time visualization of the flow field using PLIV 230 Figure 5.66 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 723 Figure 5.67 – Power spectra for Jet’s angle (left), impingement point position (middle) and pressure transducer (right) for Reynolds 857 5.4 Parametric analyses of the flow field using two different viscosity fluids Maximising the flow dynamics requires the impingement point position lie in the middle of the mixing chamber allowing the bending and stretching of the fluid laminas occupied all the chamber diameter. Ensuring the flow dynamics occur in the centre of the chamber requires thigh control of the flow rate in order to ensure not only this condition but others more problematic as it is different stoichiometries and different fluid proprieties. This is only possible with advance techniques to control the mixing process and expensive equipment not just for operating in high pressure but also so ensure low operation errors. In the present work the flow parameters are study using two different viscosity fluids each one entering in one side of the chamber in order to validate some mathematical models existent in the literature The importance of the present chapter is not only to validate some control models but to access their efficacy in a wide range of possible operation conditions, which is known to be one of main difficulty to control different fluids proprieties, which could be used to design and construct a multipurpose RIM machine. In the present section will be study two main mathematical models. The first model study is the elastic analogue model proposed by Fonte (2012) which due their versatility will be
5 Real time visualization of the flow field using PLIV 231 tested against experimental data. This step is new in literature since there is no prove of existence of this test for this model or similar models for RIM machine. Additionally the pressure model developed by Gomes (2015) which is the main driven force in the RIMcop technology into create a multipurpose machine which can allowed different wide range of fluid proprieties. This section is separated in two main parts. The first is presented the experimental PLIF setup and additional with the experimental conditions which going to be introduced all the tools and equipment used for this testing and also the conditions for each individual experience. In the end the results will be presented with addition of a brief description of the models and the sequence performed in the experiments. 5.4.1 Experimental setup 5.4.1.1 The RIM machine In the present subsection the RIM machine used in experiments is equal to the RIM present in section 5.3.1.1 in both mixing chamber and additional components, but with a small difference in the operation, in which in the present case the two fluids are storage in different reservoir tanks in order to feed different injectors 5.4.1.2 PLIF Setup In the PLIF setup big changes were performed, now combining the CCD camera used for PIV and a generator of laser beams. The transparent mixing chamber is illuminated with a laser sheet cur the chamber axially through injectors, as previous experiments, with a thin slice of 1mm. In the present work the laser source used is a Nd:YAG double-head pulsed laser from Litrom Lasers (model NanoL50-100), which emits 4ns pulsed beams with maximum energy of 400 mJ/pulse at a maximum rate of 100 Hz and correspondent wavelength of 532 nm. The transport and conversion of the laser beam into a laser sheet is the same as previously reported. Due existence of errors in the hardware, particularly in the computer memory only was possible to obtain 50 pairs of laser-synchronized frames of the illuminated flow. The record of the flow field was performed using a CCD camera (model TSI 630157) with the ability of capturing 11192 x 1600 pixels with 16 bit of greyscale at a maximum acquisition rate of 30
5 Real time visualization of the flow field using PLIV 232 Hz. The transport of the images into a computer was performed through a frame grabber (model TSI 600067). A filter was adapted into the camera lens as explain in chapter xxx to avoid reflections. The frequency of the be used in the experiments 𝑓𝑙𝑎𝑠𝑒𝑟 as defined as function of the Reynolds number (Re) and the maximum allowed frequency of the PLIF system 𝑓𝑙𝑎𝑠𝑒𝑟,𝑚𝑎𝑥, defined as 𝑓𝑙𝑎𝑠𝑒𝑟(𝑅𝐸)=𝑅𝑒 600𝑓𝑙𝑎𝑠𝑒𝑟,𝑚𝑎𝑥 (5.16) The default exposure time of 405 𝜇𝑠 is it sufficient low to ensure a freezing flow into images. The TSI software (Insight 3G version 9.0) is responsible input the experimental parameters into a synchronizer (model TSI 610035) which is used to coordinate the image shutter with the input laser beams. The image acquisition was done using a Dell Precision 690 PC with a dual core Intel Xenon CPU at 2.33 GHz and 2 GB of RAM 5.4.2 Experimental conditions In the present experiments it was set during each Reynolds in study, a constant flow rate for the right side injector while small increments of 0.1 g/s was introduced in the left side of the injector in order to dislocate the impingement point position from the left side of the chamber until the right side on the chamber. 5.4.2.1 Rheology of the solution The main objective presented in this chapter is study the parameters which directly affected the flow impingement point. Therefore it was chosen as working fluids, fluids with complete different rheology. In the left side of the chamber fluid close to 20 cPA of viscosity was used, and for the other side of the chamber a fluid with three times more viscosity was chosen in order to access the efficacy of such control models. The fluids rheology was perform for a temperature range between 20 and 34 ºC. Figure 5.68 shows the reduction of viscosity with temperature for each working fluid.
5 Real time visualization of the flow field using PLIV 233 Figure 5.68 -Rheology of both working fluids 5.4.2.2 Calibration Since the control equipment, and data acquisition system was the same used in the previous PLIF experiment, such as the Coriolis mass flow meter and the differential pressure transducer, the calibration and assessment of the error was already performed in the previous part of the present chapter. This results will then be used to evaluate the location of the impingement point in function of the control parameters of the RIM system. In the fluid with higher viscosity is added a dye marker of rhodamine 6G at the concentration of 0.4 g/l. No need for calibration of the concentration of dye marker in the solution since will be only interested for the impingement point position. 5.4.3 Results The images obtained from the PLIF technique are averaged for each experiment. Using the average image will be presented and visualization study is performed. Additionally and using the DIP MODA program the impingement point position is obtain and the results are used to compare with both mathematical models in study. Before initialize the comment of the results is necessary to define some parameters which going to be used for this validation. This parameters was already explained in previous chapter and therefore will going only be presented the mathematical equations Defining the momentum rate ratio as 𝜙𝑀=𝜌1𝑑1 2𝑢𝑖𝑛𝑗,1 2 𝜌2𝑑2 2𝑢𝑖𝑛𝑗,2 2 (5.17) 0,01 0,02 0,03 0,04 0,05 0,06 20 22 24 26 28 30 32 34 Viscosity (Pa.s) Temperature (oC) 60 Cp 20 Cp
5 Real time visualization of the flow field using PLIV 240 Figure 5.76Momentum rate ratio for Reynolds 118 (left) and Reynolds 138 (right) Figure 5.77Momentum rate ratio for Reynolds 162
5 Real time visualization of the flow field using PLIV 241 5.4.3.3 Pressure model The pressure model presented in chapter 4 will be tested for a case of two different fluid viscosities. The mathematical model present a decreasing of the pressure when the stagnation point is increase, being this increase be exponential for low values of normalized pressure. The experimental work present a impingement point position situated in high values of the normalized pressure, which is perform for the more viscosity fluid. As a conclusion the mathematical model which describes the pressure variation don’t present any match with the experiments for fluids with different viscosities, and therefore the author recommends in new works the study of this methodology of control the process by developing a new fitting mathematical model which can be used for a wide range of operational conditions such as different diameters, different fluid proprieties and which can be used for a variety of Reynolds number. Figure 5.78Pressure model validation for Reynolds 72 (left) and Reynolds 81 (right)
5 Real time visualization of the flow field using PLIV 242 Figure 5.79Pressure model validation for Reynolds 118 (left) and Reynolds 138 (right) Figure 5.80Pressure model validation for Reynolds 162
5 Real time visualization of the flow field using PLIV 243 5.5 Conclusions The main presented in the present chapter is divided into two main sections. The first section analyse the results provided by PLIF experiments using a high speed camera to record the main flow scales. The second work tends to analyse the two main mathematical models to predict the impingement point position by knowing two different parameters, namely the mass flow rate and the Reynolds number in each injector for the elastic analogue model, and the differential pressure between both injectors for the pressure jets model. According the first part of the work it was observed the FFT of the pressure transducer didn’t obtain any physical meaning of the frequencies of the main flow scale and therefore it was necessary to go to previous works where this frequency was already obtain for a set of Reynolds number. It became clear with this work the main flow scale correspond to the impingement surface rotation rather than the impingement point position, which could be due the magnitude of both variables or the tracking of the flow, while the impingement point position is 1D compared with the 2D for the jets angle. Instantaneous flow field visualization identified three main flow dynamics: segregated flow, periodic segregated flow and self-sustainable chaotic regime. The present work identify for the self-sustainable chaotic regime three main different flow conditions: the periodic, flow frequency random oscillation and high frequency random oscillation. Using fluids with different viscosities entering the mixing chamber by different jets two mathematical models was used to access the viability of such techniques in order to be used as control model for a wide variety of operating conditions such as different stoichiometry, different fluids rheology’s and geometrical conditions. Therefore the elastic analogue model and the pressure model was tested. It was conclude the elastic analogue model suited perfectly in the results provided by the experimental data. In the other hand the jets pressure model, which is the main concept of the RIMcop technology, didn’t fit in the results since completely different curves was obtain for each case which proves the complexity and difficulty of modelling the impingement point position with this technology and therefore improvement is necessary to be made.
244 6 Final remarks 6.1 Conclusions The present work is divided into three main chapters according to the type and methodology to access different aspects of the flow field in RIM mixing head. First results chapter access using numerical tools the flow field in a RIM machine which was study in order to quantify and qualify the scales of the flow. In the second results chapter experimental technique, PIV, was perform in a front plane of the mixing chamber. This allowed to quantify and qualify the flow field behaviour and the access the effectiveness of two main mathematical equations models to predict the impingement point position base on a various flow parameters. In the last results chapter the experimental technique, PLIF, was performed for the same front plane of the mixing chamber. By using an algorithm, DIP MODA, the impingement point position and the jets angle was calculated for the domain of each experiment. Additionally using the same technique it was access the effectiveness of the same mathematical equations models, using different fluid viscosities entering opposed in the mixing chamber. In the present chapter will be resumed the main conclusions for each topic studied, for the first time in the literature, presented by the present work. The main topics are the flow hydrodynamic behaviour and scale containing eddies, frequency of the larger eddies, the mixture using a non-diffusivity model (VOF) and using an large size grid and in the end
6 Final remarks 245 mixing control using two different mathematical models the elastic analoge model and the pressure model. 6.2 DIP MODA for tracking impingement point and jets angle Digital Image Processing Multiobjective Optimization Deterministic Algorithm (DIP MODA) was created in order to obtain the jets impingement point and the jets angle. This code is divided into two parts, where the first part is the image processing of the images obtained from the PLIF experiments for a defined set images. It was used a high speed camera which allowed to track, in real time, the flow behaviour. A multiobjective optimization algorithm allowed to track the impingement point position, by searching in a line central to the jets centre, the value of the interface, combined by others variables using a bell shape weighted method. The same idea was performed for the jets angle, where different variables used and the interface was track and the angle calculated using a linear regression. This code prove to be very effective by tested against a benchmark test and afterwards applied to the real problem in study. 6.3 Hydrodynamic behaviour 6.3.1 Visualization studies of the hydrodynamic flow behaviour Numerical studies was performed for Reynolds 150, 200, 250 and 300. All the flow hydrodynamic scales was calculated since the flow was solved until the Kolmogorov scale. The visualization studies shows with an increasing of the Reynolds number an increase of the oscillation of the impingement surface due impingement jets. Also with this increasing it is observable a breaking of the normal oscillatory motion of the impingement surface which cause a creation of small vortices in the impingement region rather than vortices flowing along the chamber height. This indicate a large energy dissipation in this small region and with an increasing of this region the size of this dissipation tends to small due increasing of the creation of this small vortices. Generally speaking it is observable a large oscillation, in terms on angle in the front plane rather than the top plane. Additionally is observable the large size eddies is majority contain in the front plane since the top plane the eddies are containing in the chamber diameter while the front plane allowed bigger area available to jets motion.
6 Final remarks 246 6.3.2 Qualitative and quantitative analyses of 𝒌 and 𝜺 The turbulent kinetic energy and the turbot energy dissipation was calculated for each numerical studies performed. It is visualized for each variable the location is in the impingement area since it is here exist the biggest fluctuation velocity field and the higher velocity gradients of the flow field. In term of the vertical dimension of the higher concentration of this variables it is observed a small growing in terms of the magnitude from 150 to 200 while in the range 200 to 250 remain almost constant. The biggest rate of growing is observed in range 250 to 300. By normalizing this variables by the value in the central point it is observed a fast growing rate for 𝒌 with a high decreasing of 𝜺 for the range 150 to 200. Generally speaking is observed a opposed variation of this variables with a slight increasing of the Reynolds number 6.3.3 Scales containing eddies Using Fast Fourier Transformed was used to calculate the scale containing eddies which was compared with the chamber dimensions. Two main cascades was obtain: the inverse energy cascade with slope -5/3 which means the larger low scales are feed with low dimension scales and the direct entrosphy cascade which means the flow flows in a positive direction where the larger scales give energy to the lower scales. In the dissipation range all the flow scales are dissipated into heat in the size equal as the kolmogoroff scale. In the inertial subrange the inversion of cascades from a slope -5/3 into -3 is a point where the energy is injected which correspond a size equal the injector’s diameter. The larger energy containing scale have an eddy size between the chamber diameter and the height of the chamber size. 6.3.4 Frequency analyses Using the DIP MODA the jets angle and the impingement point position was obtain for each set of capture images referred to and experiment. Additionally was also recorder the value differential pressure between each injector. FFT was perform for each case and the results compared with documented work. It was proven the frequency of the jets angle matches the documented work frequency until Reynolds close to 300. After this value it was not observed any significant match, since complex phenomena’s in impingement jet iteration.
6 Final remarks 247 All the others variables such as the impingement point position and the differential pressure didn’t obtained any relevant data 6.4 Mixture using a non-diffusivity model (VOF) In the present work Volume of Fluid (VOF) model was used to access the study of mixing mechanism using a non-diffusivity model. Due the high computational effort of such method it was choose to setup a mesh length equal to the Kolmogorov rather than the bachelor scale. In the present work it was clear after the large scale vortices created due impingement the reduction of scale was more than the mesh size and therefore numerical diffusion dominate the entire domain. An adaptive grid gradient based was performed to effectiveness obtain smaller scales but the tremendous computational effort was necessary to abort such technique due computational limitations. 6.5 RIM process control 6.5.1 Elastic analogue model The elastic analogue model was verified for a variety of Reynolds number with equal fluids rheology’s, with 20 mPa.s in each injector using PIV and with different fluid rheology’s with 20 and 60 mPa.s of viscosity respectively in each injector using in this case PLIF technique. The elastic analogue model was proven to successful determine the impingement point position for a wide fluid proprieties and Reynolds number. Previous works was already validate for different injectors diameters which concludes this model to be one of the most well adapted model in the literature 6.5.2 Pressure jets model The pressure jets model was tested in the same conditions as the elastic analogue model and using the same methodologies. The present work concludes this model, which only model half of the chamber, didn’t perform very well with the same fluid theologies in each injector but it was acceptable for the impingement point close the walls and the centre of the chamber. For the test with different fluid proprieties the model didn’t pedict and even behave opposite as the experimental work. The author then conclude further investigation in the development of such model in order to be used with RIMcop technology for a wide range of control conditions.
6 Final remarks 248 6.6 Future work recommendations The present work is divided into two main research. For the first research related to study the hydrodynamics mixing mechanism for the fluid flow system. For this case some suggestions for future work are presented. Create a mechanical/electromechanical system to oscillate the flow rate in the jets. New works in the field of combustion uses this methodology to improve the mixing efficiency. In the present work this methodology can improve the formation of large eddies, which going to break into small eddies for a longer distance inside the mixing chamber, improving the overall mixing performance particularly at low Reynolds number. Further investigate, using experimental and numerical techniques, the mixing between fluids com different rheology’s, particularly with high viscosity ratios which is usually found in RIM process. Implement a non-diffusive dynamic adaptive grid to model the flow field inside a 3D model of the mixing chamber in order to obtain the full hydrodynamic behaviour, and the full mixing behaviour by model the flow field until the bachelor scale. The second research theme it was the validation of some mathematical models to be used as control the mixing process. For this case some suggestions for future work is presented Improvement of the control of the mixing process using the RIMcop technology for a wide variety of situations requires a development of a new mathematical model which can be used for more than one control variable, such as Reynolds number, impingement point position, jets ratio and others. Development of a multi-variable control using the elastic analogue model to predict the position of the impingement point, which can be coupled with the differential pressure of the impingement jets and can be validated by a real time image processing technique with the same concept used to calculate the impingement point position using the DIP MODA algorithm
249 References Adrover, A., Cerbelli, S., & Giona, M. (2002). A Stectral Approach to Reaction/Diffusion Kinetics in Chaotic Flows. Computers and Chemical Engineering, 26, 125-139. Angst, W., Bourne, J. R., & Sharma, R. N. (1982). Mixing and Fast Chemical Reaction - V. Influence of Diffusion Within the Reactor Zone on Selectivity. Chemical Engineering Science, 37(8), 1259-1264. Baldyga, J., & Bourne, J. R. (1983). Distribution of Striation Thickness from Impingement Mixers in Reaction Injection Molding. Polymer Engineering and Science, 23(10), 556-559. Baldyga, J., & Bourne, J. R. (1984a). A Fluid Mechanical Approach to Turbulent Mixing and Chemical Reaction Part Iii Computational and Experimental Results for the New Micromixing Model. Chemical Engineering Communications, 28(4-6), 259-281. doi: Doi 10.1080/00986448408940137 Baldyga, J., & Bourne, J. R. (1984c). A Fluid Mechanical Approach to Turbulent Mixing and Chemical Reaction PartII Microming in the Light of Turbulence Theory. Chemical Engineering Communications, 28, 243-258. Baldyga, J., & Bourne, J. R. (1984f). Mixing and Fast Chemical Reaction-VIII. Initial Deformation of Material Elements in Isotropic, Homogeneous Turbulence. Chemical Engineering Science, 39(2), 329-334.
References 256 Sebastian, D. H., & Boukobbal, S. (1986). Mixhead Parameters Governing Impingement Mixing Effectiveness for Polyurethane Reactive Injection Molding Processes. Polymer Process Engineering, 4(1), 53-70. Speziale, C. G. (1991). Analytical Methods for the Development of Reynolds-Stress Closures in Turbulence. Annual Reviews of Fluid Mechanics, 23, 107-157. Teixeira, A. M. (2000). Escoamento na cabeça de mistura de uma máquina de RIM Caracterização experimental por LDA e simulação dinâmica por CFD. (PhD Thesis), FEUP - Faculdade de Engenharia da Universidade do Porto, Porto. Teixeira, A. M., Santos, R. J., Costa, M. R. P. F. N., & Lopes, J. C. B. (2005). Hydrodynamics of the mixing head in RIM: LDA flow-field characterization. AIChE Journal, 51(6), 1608-1619. doi: 10.1002/aic.10454 Tennekes, H., & Lumley, J. L. (1997). A First Course in Turbulence. Cambridge: MIT Press. Toor, H. L. (1969). Turbulent mixing of two species with and without chemical reactions. Industrial & Engineering Chemistry Fundamentals, 8(4), 655-659. Tosun, G. (1987). A Study of Micromixing in Tee Mixers. Industrial & Engineering Chemistry Research, 26, 1184-1193. Tropea, C., Yarin, A. L., & Foss, J. F. (2007). Springer Handbook of Experimental Fluid Mechanics Springer. Tsai, K., & Fox, R. O. (1994). PDF Simulation of a Turbulent Series-Parallel Reaction in an Axismmetric Reactor. Chemical Engineering Science, 49(24B), 5141-5158. Tucker, C. L., & Suh, N. P. (1980). Mixing for reaction injection molding. I. Impingement mixing of fiber suspensions. Polymer Engineering and Science, 20(13), 887-898. doi: 10.1002/pen.760201308 Tucker, C. L., & Suh, N. P. (1980). Mixing for Reaction Injection Molding. I. Impingement Mixing of Liquids. Polymer Engineering and Science, 20(13), 875-886. doi: 10.1002/pen.760201307 Unger, D. R., & Muzzio, F. J. (1999). Laser-induced Fluorescence Technique for the Quantification of Mixing in Impinging Jets. AIChE Journal, 45(12), 2477-2486. doi: 10.1002/aic.690451203 Unger, D. R., Muzzio, F. J., & Brodkey, R. S. (1998). Experimental and Numerical Characterization of Viscous Flow and Mixing in an Impinging Jet Contactor. The Canadian Journal of Chemical Engineering, 76, 546-555. Vassilatos, G., & Toor, H. L. (1965). Second-Order Chemical Reactions in a Nonhomogeneous Turbulent Fluid. AIChE Journal, 11(4), 666-673. doi: 10.1002/aic.690110419
References 257 Villermaux, J., & David, R. (1983). RECENT ADVANCES IN THE UNDERSTANDING OF MICROMIXING PHENOMENA IN STIRRED REACTORS. Chemical Engineering Communications, 21(1-3), 105-122. doi: 10.1080/00986448308940280 Villermaux, J., Falk, L., Fournier, M.-C., & Detrez, C. (1992). Use of Parallel Competing Reactions to Characterize Miromixing Efficiency. Paper presented at the AIChE Symposium Series. Wang, D. M., & Tarbell, J. M. (1993). Closure Models for Turbulent Reacting Flows with a Nonhomogeneous Concentration Field. Chemical Engineering Science, 48(23), 3907-3920. Wehrmeyer, J. A., Cheng, Z., Mosbacher, D. M., Pitz, R. W., & Osborne, R. (2002). Opposed Jet Flames of Lean or Rich Premixed Propane-Air Reactants versus Hot Products. Combustion and Flame, 128, 232-241. Weinstein, H., & Adler, R. J. (1967). Micromixing Effects in Continuous Chemical Reactors. 22, 65-75. Wood, P., Hrymak, A. N., Yeo, R., Johnson, D. A., & Tyagi, A. (1991). Experimental and Computational Studies of the Fluid Mechanics in an Opposed Jet Mixing Head. Physics of Fluids A, 3(5), 1362-1368. Zalc, J. M., & Muzzio, F. J. (1999). Parallel-Competitive Reactions in a Two-Dimensional Chaotic Flow. Chemical Engineering Science, 54, 1053-1069. Zhao, Y., & Brodkey, R. S. (1998a). Averaged and Time-resolved Full-field (threedimensional), Measurements of Unsteady Opposed Jets. The Canadian Journal of Chemical Engineering, 76, 536-545. Zhao, Y., & Brodkey, R. S. (1998d). Particle Paths in Three-dimensional Flow Fields as a Means of Study: Opposing Jet Mixing System. Powder Technology, 100(2-3), 161-165. doi: 10.1016/S0032-5910(98)00136-3 Zwietering, N. (1959). The Degree of Mixing in Continuous Flow Systems. Chemical Engineering Science, 11(1), 1-15.
258 A Mathematical discretization techniques The present appendix present the mathematical techniques used in the special and temporal discretization to evaluate the dependent and independent variables. Also the gradient discretization is presented. The intent of this chapter is provide to the reader background of the mathematical discretization implementation used in the CFD program. A.1 Space discretization A.1.1 MUSCL scheme The MUSCL scheme uses the local cell averages values of the quantities in transport at time 𝑡 to evaluate the value of such quantities in cell edges at 𝑡+∆𝑡/2 and afterward to predict values at 𝑡+∆𝑡. This process is done in four steps I. Calculating the slopes: Local cell averages, for the desired quantity, can be represented as 𝑊(𝑥,𝑡)=𝑊𝑖+𝑥−𝑥𝑖 ∆𝑥 𝛿𝑊𝑖 (A.1) Where 𝑥𝑖−12<𝑥<𝑥𝑖−12 , 𝑥𝑖≡(𝑥𝑖+12+𝑥𝑖−12)/2 , ∆𝑥≡𝑥𝑖+12+𝑥𝑖−12 and 𝛿𝑊𝑖 are the reconstructed slope calculated as
Mathematical discretization techniques 259 𝛿𝑊𝑖=𝑎𝑣𝑒 (𝑊𝑖−𝑊𝑖−1,𝑊𝑖+1−𝑊𝑖) (A.2) Where 𝑎𝑣𝑒 is the desired averaging function, where 𝑎𝑣𝑒(𝑎,𝑏)=(𝑎+𝑏)/2 lead to the central differences scheme and 𝑎𝑣𝑒(𝑎,𝑏)=0 leads to the zero slope II. Predict solution ate half step: Using the sloples previously calculated the solution will be progressed half time-step ∆𝑡/2, where the predict values in the cell can be calculated as 𝑉𝑗=𝑉𝑗−∆𝑡 2∆𝑥𝐴𝑝(𝑉𝑗)𝛿𝑉𝑗 (A.3) III. Predict values at cell edges: Using the predict solution, the values at cell edges will be calculated assuming a linear solution 𝑊𝑖−12 +=𝑊 𝑖− 𝛿𝑊𝑖/2 (A.4) 𝑊𝑖+12 −=𝑊 𝑖+ 𝛿𝑊𝑖/2 (A.5) IV. Calculate values at cell edges: Used Reimann solver to update desired variables at time 𝑡 𝑉𝑖𝑛+1=𝑉𝑖𝑛−∆𝑡 ∆𝑥(𝐹𝑖+12−𝐹𝑖−12) (A.6) Where 𝐹𝑖±1/2 correspond the numerical fluxes calculated form values predicted in III 𝐹𝑖−12≡𝐹(𝑊𝑖+12 −,𝑊𝑖−12 +) (A.7) A.1.2 Coupled scheme Coupled scheme is an full implicit algorithm which solved the momentum equation and the pressure-based continuity equation in simultaneous. This process uses the implicit discretization of the pressure gradients using the momentum equation and an implicit discretization for the mass flux. The momentum equation is in the form 𝑎𝑝𝑢=∑𝑎𝑛𝑏𝑢𝑛𝑏 𝑛𝑏 +∑𝑝𝑓𝐴⋅î+𝑆 𝑓 (A.8) Where the pressure values at the faces is interpolated using the momentum equation
Mathematical discretization techniques 260 𝑝𝑓=𝑃𝑐0 𝑎𝑝,𝑐0+𝑃1 𝑎𝑝,𝑐1 1 𝑎𝑝,𝑐0+1 𝑎𝑝,𝑐1 (A.9) Resulting in a smooth variation of the pressure. This equation is also known as standard pressure scheme. From the momentum equation the pressure gradient for component 𝑘 is calculated as ∑𝑝𝑓𝐴𝑘= 𝑓∑𝑎𝑢𝑘𝑝𝑗 𝑓 (A.10) Where 𝑎𝑢𝑘𝑝𝑗 is the coefficient provided from the Gauss divergence theorem and the coefficient of the pressure interpolation scheme. The momentum equation for a given cell 𝑖 and for a component 𝑢𝑘 is defined as ∑𝑎𝑖𝑗 𝑢𝑘𝑢𝑙𝑢𝑘𝑗 𝑗+∑𝑎𝑖𝑗 𝑢𝑘𝑝𝑗 𝑗=𝑏𝑖𝑢𝑘 (A.11) The value of the velocity in a cell face isn’t linear average but instead a weighted method is used based on the coefficients of the momentum equation, therefore the flux over a face cell can be calculated as 𝐽𝑓 =𝜌𝑓𝑎𝑝,𝑐0𝑣𝑛,𝑐0+𝑎𝑝,𝑐1𝑣𝑛,1 𝑎𝑝,𝑐0+𝑎𝑝,𝑐1+𝑑𝑓((𝑝𝑐0+(∆𝑝)𝑐0⋅𝑟0 )−(𝑝𝑐1+(∆𝑝)𝑐1⋅𝑟1 ))𝐽𝑓 +𝑑𝑓(𝑝𝑐0−𝑝𝑐1) (A.12) Where 𝑝𝑐0,𝑝𝑐1 the pressure at two cells on the same side of the face and 𝑣𝑛,𝑐0,𝑣𝑛,𝑐1 are the normal velocity ate the same location. 𝐽𝑓 is the influence of the velocity in this cells. 𝑑𝑓is function of average momentum equation coefficients. Substituting the flux formulation in the momentum equation results in ∑∑𝑎𝑖𝑗𝑝𝑢𝑘𝑢𝑘𝑗 𝑗𝑘 +∑𝑎𝑖𝑗 𝑝𝑝 𝑝𝑗 𝑗=𝑏𝑖𝑝 (A.13) This equation is then transform into a matricial form in order to be calculated
Mathematical discretization techniques 261 ∑[𝐴]𝑖𝑗𝑋𝑗 𝑗=𝐵𝑖 (A.14) Where each matrix have the following form 𝐴𝑖𝑗= [ 𝑎𝑖𝑗 𝑝𝑝 𝑎𝑖𝑗 𝑝𝑢 𝑎𝑖𝑗 𝑝𝑣 𝑎𝑖𝑗 𝑝𝑤 𝑎𝑖𝑗 𝑢𝑝 𝑎𝑖𝑗 𝑢𝑢 𝑎𝑖𝑗 𝑢𝑣 𝑎𝑖𝑗 𝑢𝑤 𝑎𝑖𝑗 𝑣𝑝 𝑎𝑖𝑗 𝑣𝑢 𝑎𝑖𝑗 𝑣𝑣 𝑎𝑖𝑗 𝑣𝑤 𝑎𝑖𝑗 𝑤𝑝 𝑎𝑖𝑗 𝑤𝑢 𝑎𝑖𝑗 𝑤𝑣 𝑎𝑖𝑗 𝑤𝑤 ] (A.15) 𝑋𝑗=[𝑝𝑖´ 𝑢𝑖´ 𝑣𝑖´ 𝑤𝑖´] (A.16) 𝑋𝑗= [ −𝑟𝑖𝑝 −𝑟𝑖𝑢 −𝑟𝑖𝑢 −𝑟𝑖𝑤 ] (A.17) A.1.3 SIMPLE-C scheme The SIMPLE-C algorithm is an adaptation of the SIMPLE scheme for velocity-pressure coupling. This method, introduced by Pantankar(), initially start to predict the velocity components and the pressure (𝑢,𝑣,𝑤,𝑝). After the guess is performed momentum equations is used to calculate the velocity component’s (𝑢∗,𝑣∗,𝑤∗) based on the momentum equation which afterwards is used to calculate the pressure 𝑝′. 𝑎𝑒𝑢𝑒′=∑𝑎𝑛𝑏 𝑛𝑏 𝑢𝑛𝑏 ′+∆𝑦(𝑝𝑝′−𝑝𝑒′) (A.18) 𝑎𝑒𝑣𝑒′=∑𝑎𝑛𝑏 𝑛𝑏 𝑣𝑛𝑏 ′+∆𝑥(𝑝𝑠′−𝑝𝑝′) (A.19) 𝑎𝑒𝑤𝑒′=∑𝑎𝑛𝑏 𝑛𝑏 𝑤𝑛𝑏 ′+∆𝑥(𝑝𝑤 ′−𝑝𝑝′) (A.20) Where 𝑒,𝑠,𝑤,𝑝 correspond to the east, south, direction perpendicular to e,s plane and the position of the point respectively. Approximation of the velocity correction are made ignoring the first term of the momentum equation ∑𝑎𝑛𝑏 𝑛𝑏 𝑢𝑛𝑏 ′, ∑𝑎𝑛𝑏 𝑛𝑏 𝑣𝑛𝑏 ′, ∑𝑎𝑛𝑏 𝑛𝑏 𝑤𝑛𝑏 ′. The
Mathematical discretization techniques 262 correct velocities are substituting the continuity equation yielding a discrete pressure correction equation 𝑎𝑝𝑝𝑝′=∑𝑎𝑛𝑏 𝑛𝑏 𝑝𝑛𝑏+𝑏 (A.21) After is obtain the corrected velocity and pressure 𝑝′,𝑢′,𝑣′,𝑤′, the variables are updated using: 𝑢=𝑢∗+𝛼𝑝𝑢′ 𝑣=𝑣∗+𝛼𝑝𝑣′ 𝑤=𝑤∗+𝛼𝑝𝑤′ 𝑝=𝑝∗+𝛼𝑝𝑝′ (A.22) Where 𝛼𝑝 is the under relaxation factor being necessary since the nonlinear nature of the equations. Whit all the variables corrected the face flux cell is then calculated using 𝐽𝑓=𝐽𝑓∗+𝑑𝑓(𝑝𝑐0 ′−𝑝𝑐1 ′) (A.23) Where the coefficient 𝑑𝑓 is defined as 𝑑𝑓=(𝑎𝑝−∑ 𝑎𝑛𝑏 𝑛𝑏 ) (A.24) A.2 Time discretization Finite–difference method are a numerical method used to solve a differential equation based on the approximation into a finite space/time. Initial a function 𝜑=𝑓(𝑥) is pointed approximated using Taylor series expansion into a polynomial form in the space discretization 𝑓(𝑥,𝑡) =𝑓(𝑎)|𝑡=1+𝑓′(𝑎) 1! (𝑥−𝑎)|𝑡=𝑖+𝑓(2)(𝑎) 2! (𝑥−𝑎)2|𝑡=𝑖+⋯+𝑓(𝑛)(𝑎) 𝑛! (𝑥−𝑎)𝑛|𝑡=𝑖 (A.25) Where 𝑛 is the order of the interpolation and 𝑎 the point where the function is discretize. Solving 𝑓(𝑥,𝑡), time dependent, can be done using three methodologies, full implicit, full explicit and Crank-Nicolson scheme (central differences). 𝑥𝑡+1−𝑥𝑡 ∆𝑡 =𝑘 . 𝑓(𝑥𝑡+1) + (1−𝑘) .𝑓(𝑥𝑡) (A.26) Where 𝑘=0 for full explicit scheme, 𝑘=1 for full implicit scheme and 𝑘=0.5 for central difference scheme. The function 𝑓(𝑥,𝑡) will be then integrated over a discrete domain
Mathematical discretization techniques 263 ∫ 𝑓(𝑥)𝑑𝑡=[𝑘.𝑓𝑥𝑡+∆𝑡+(1−𝑘).𝑓𝑥𝑡]∆𝑡 𝑡+∆𝑡 𝑡 (A.27) Explicit scheme the function is evaluated at the current time (𝑡+1) provided by previous values 𝑡 and added the current increment value ∆𝑡 𝑓(𝑥𝑡), resulting in 𝑥𝑡+1=𝑥𝑡+∆𝑡 𝑓(𝑥𝑡) (A.28) This method is fast since only is necessary to add value, but conditionally stable, since the time step marching is limited by the Courant–Friedrichs–Lewy condition. Implicit scheme evaluate the current time 𝑥𝑡+1 based on the previous values 𝑥𝑡 and adding values correspondent to the current time values ∆𝑡 𝑓(𝑥𝑡+1),. This methodology implies solving a linear system of equations resulting in a computational effort to invert the constant’s matrix. Although this factor this method is unconditionally stable for any increment of time 𝑥𝑡+1=𝑥𝑡+∆𝑡 𝑓(𝑥𝑡+1) (A.29) This formulation is referred for a first order implicit, but can be extender for higher orders 𝑓(𝑥𝑡)=∑(𝑡𝑘)(−1)𝑡−𝑘𝑓(𝑥+𝑘) 𝑟 𝐾=0 (A.30) Where (𝑡𝑘) is the binomial coefficient wheke the row of pascal triangle provide the coefficient for 𝑘. A.3 Gradient discretization A.3.1 Least Squares cell-Based Gradient evaluation This method assume a linear variation of a propriety 𝜙 between a cell 𝑐0 and adjacent cells 𝐶𝑖 with the distance between the adjacent cell’s center 𝛿𝑟𝑖. The cell centroid evaluation can be expressed as (∆𝜙)𝑐0∙∆𝑟𝑖=(𝜙𝑐𝑖−𝜙𝑐0) (A.31) Where the vectorial distance between the two center cells can be calculated as ∙∆𝑟𝑖=∆𝑥𝑖 𝐢+∆𝑦𝑖 𝐣+∆𝑧𝑖 𝐤 (A.32)
Mathematical discretization techniques 264 Where the change of proprieties between the two center cells can be estimated as ∆𝑥𝑖 𝜕𝜙 𝜕𝑥|0+∆𝑦𝑖 𝜕𝜙 𝜕𝑦|0+∆𝑧𝑖 𝜕𝜙 𝜕𝑧|0=𝜙𝑖−𝜙0 (A.33) Which can be transform into an algebraic system of equations 𝐌 𝐝=∆𝜙 (A.34) Being each one variable transform into matricial form ∆𝜙=[𝜙1−𝜙0 𝜙2−𝜙0 ⋮ 𝜙𝑖−𝜙0] 𝒅= [ 𝜕𝜙 𝜕𝑥|0 𝜕𝜙 𝜕𝑦|0 𝜕𝜙 𝜕𝑧|0 ] 𝐌=[∆𝑥1∆𝑦1∆𝑧1 ∆𝑥3∆𝑦3∆𝑧3 ⋮ ⋮ ⋮ ∆𝑥𝑖∆𝑦3∆𝑧3] (A.35) Since the system is over-determinated, a weighted method, based on the Gram-Schmidt is used, where the decomposition of the matrix’s yields a weight matrix for each cell, and therefore three individual weights, for each vector in the Cartesian reference frame will be calculated in order to produce values at each face of the cell 𝐶0 , (𝜙𝑥)𝑐0=∑𝑊𝑥𝑖0∙(𝜙𝑐𝑖−𝜙𝑐0) 𝑛 𝑖=1 (A.36) (𝜙𝑦)𝑐0=∑𝑊𝑦𝑖0∙(𝜙𝑐𝑖−𝜙𝑐0) 𝑛 𝑖=1 (A.37) (𝜙𝑧)𝑐0=∑𝑊𝑧𝑖0∙(𝜙𝑐𝑖−𝜙𝑐0) 𝑛 𝑖=1 (A.38)
265 B Spatial Evolution of the flow along the relevant time domain in the probe points B.1 Reynolds 150 Figure B.1 – Temporal evolution for vx, vy and vz at Reynolds 150 in p0
B Spatial evolution of the flow along the relevant time domain in probe points 272 Figure B.14 – Temporal evolution for vx, vy and vz at Reynolds 200 in p2 Figure B.15 – Temporal evolution for vx, vy and vz at Reynolds 200 in p3
B Spatial evolution of the flow along the relevant time domain in probe points 273 Figure B.16 – Temporal evolution for vx, vy and vz at Reynolds 200 in p4 Figure B.17 – Temporal evolution for vx, vy and vz at Reynolds 200 in p5
B Spatial evolution of the flow along the relevant time domain in probe points 274 Figure B.18 – Temporal evolution for vx, vy and vz at Reynolds 200 in p6 Figure B.19 – Temporal evolution for vx, vy and vz at Reynolds 200 in p7
B Spatial evolution of the flow along the relevant time domain in probe points 275 Figure B.20 – Temporal evolution for vx, vy and vz at Reynolds 200 in p8 Figure B.21 – Temporal evolution for vx, vy and vz at Reynolds 200 in p9
B Spatial evolution of the flow along the relevant time domain in probe points 276 Figure B.22 – Temporal evolution for vx, vy and vz at Reynolds 200 in p10 B.3 Reynolds 250 Figure B.23 – Temporal evolution for vx, vy and vz at Reynolds 250 in p0
B Spatial evolution of the flow along the relevant time domain in probe points 277 Figure B.24 – Temporal evolution for vx, vy and vz at Reynolds 250 in p1 Figure B.25 – Temporal evolution for vx, vy and vz at Reynolds 250 in p2
B Spatial evolution of the flow along the relevant time domain in probe points 278 Figure B.26 – Temporal evolution for vx, vy and vz at Reynolds 250 in p3 Figure B.27 – Temporal evolution for vx, vy and vz at Reynolds 250 in p4
B Spatial evolution of the flow along the relevant time domain in probe points 279 Figure B.28 – Temporal evolution for vx, vy and vz at Reynolds 250 in p5 Figure B.29 – Temporal evolution for vx, vy and vz at Reynolds 250 in p6
B Spatial evolution of the flow along the relevant time domain in probe points 280 Figure B.30 – Temporal evolution for vx, vy and vz at Reynolds 250 in p7 Figure B.31 – Temporal evolution for vx, vy and vz at Reynolds 250 in p8
B Spatial evolution of the flow along the relevant time domain in probe points 281 Figure B.32 – Temporal evolution for vx, vy and vz at Reynolds 250 in p9 Figure B.33 – Temporal evolution for vx, vy and vz at Reynolds 250 in p10
268 C Power spectra for velocity components in probe points
C Power spectra for velocity components in probe points 269 C.1 Reynolds 150 Figure C.1 – Power Spectra for vx, vy and vz at Reynolds 150 in p0
C Power spectra for velocity components in probe points 270 Figure C.2 – Power Spectra for vx, vy and vz at Reynolds 150 in p1
C Power spectra for velocity components in probe points 271 Figure C.3 – Power Spectra for vx, vy and vz at Reynolds 150 in p2
C Power spectra for velocity components in probe points 272 Figure C.4 – Power Spectra for vx, vy and vz at Reynolds 150 in p3
C Power spectra for velocity components in probe points 273 Figure C.5 – Power Spectra for vx, vy and vz at Reynolds 150 in p4
C Power spectra for velocity components in probe points 274 Figure C.6 – Power Spectra for vx, vy and vz at Reynolds 150 in p5
C Power spectra for velocity components in probe points 275 Figure C.7 – Power Spectra for vx, vy and vz at Reynolds 150 in p6
276 Figure C.8 – Power Spectra for vx, vy and vz at Reynolds 150 in p7
C Power spectra for velocity components in probe points 277 Figure C.9 – Power Spectra for vx, vy and vz at Reynolds 150 in p8
C Power spectra for velocity components in probe points 284 Figure C.16 – Power Spectra for vx, vy and vz at Reynolds 200 in p4
C Power spectra for velocity components in probe points 285 Figure C.17 – Power Spectra for vx, vy and vz at Reynolds 200 in p5
C Power spectra for velocity components in probe points 286 Figure C.18 – Power Spectra for vx, vy and vz at Reynolds 200 in p6
C Power spectra for velocity components in probe points 287 Figure C.19 – Power Spectra for vx, vy and vz at Reynolds 200 in p7
C Power spectra for velocity components in probe points 288 Figure C.20 – Power Spectra for vx, vy and vz at Reynolds 200 in p8
C Power spectra for velocity components in probe points 289 Figure C.21 – Power Spectra for vx, vy and vz at Reynolds 200 in p9
C Power spectra for velocity components in probe points 290 Figure C.22 – Power Spectra for vx, vy and vz at Reynolds 200 in p10
C Power spectra for velocity components in probe points 291 C.3 Reynolds 250 Figure C.23 – Power Spectra for vx, vy and vz at Reynolds 250 in p0
C Power spectra for velocity components in probe points 292 Figure C.24 – Power Spectra for vx, vy and vz at Reynolds 250 in p1
C Power spectra for velocity components in probe points 293 Figure C.25 – Power Spectra for vx, vy and vz at Reynolds 250 in p2
C Power spectra for velocity components in probe points 300 Figure C.32 – Power Spectra for vx, vy and vz at Reynolds 250 in p9
C Power spectra for velocity components in probe points 301 Figure C.33 – Power Spectra for vx, vy and vz at Reynolds 250 in p10
C Power spectra for velocity components in probe points 302 C.4 Reynolds 300 Figure C.34 – Power Spectra for vx, vy and vz at Reynolds 300 in p0
C Power spectra for velocity components in probe points 303 Figure C.35 – Power Spectra for vx, vy and vz at Reynolds 300 in p1
C Power spectra for velocity components in probe points 304 Figure C.36 – Power Spectra for vx, vy and vz at Reynolds 300 in p2
C Power spectra for velocity components in probe points 305 Figure C.37 – Power Spectra for vx, vy and vz at Reynolds 300 in p3
C Power spectra for velocity components in probe points 306 Figure C.38 – Power Spectra for vx, vy and vz at Reynolds 300 in p4
C Power spectra for velocity components in probe points 307 Figure C.39 – Power Spectra for vx, vy and vz at Reynolds 300 in p5
C Power spectra for velocity components in probe points 308 Figure C.40 – Power Spectra for vx, vy and vz at Reynolds 300 in p6
C Power spectra for velocity components in probe points 309 Figure C.41 – Power Spectra for vx, vy and vz at Reynolds 300 in p7
D Average flow field visualization for different viscosity ratios 316 𝜙𝐾= 1.9473 𝜙𝐾= 2.4825 𝜙𝐾= 3.3024 𝜙𝐾= 3.6071 Figure D.3-Average flow field visualization, using PLIF, for Reynolds 72
D Average flow field visualization for different viscosity ratios 317 𝜙𝐾= 4.3171 𝜙𝐾= 6.5673 𝜙𝐾= 7.1909 Figure D.4-Average flow field visualization, using PLIF, for Reynolds 72
D Average flow field visualization for different viscosity ratios 318 D.2 Reynolds 81 𝜙𝐾= 0.0185 𝜙𝐾= 0.0438 𝜙𝐾= 0.1048 𝜙𝐾= 0.1422 Figure D.5-Average flow field visualization, using PLIF, for Reynolds 81
D Average flow field visualization for different viscosity ratios 319 𝜙𝐾= 0.3358 𝜙𝐾= 0.4181 𝜙𝐾= 0.5954 𝜙𝐾= 0.9136 Figure D.6-Average flow field visualization, using PLIF, for Reynolds 81
D Average flow field visualization for different viscosity ratios 320 𝜙𝐾= 1.0786 𝜙𝐾= 1.3867 𝜙𝐾= 1.7419 𝜙𝐾= 2.3349 Figure D.7-Average flow field visualization, using PLIF, for Reynolds 81
D Average flow field visualization for different viscosity ratios 321 𝜙𝐾= 2.7144 𝜙𝐾= 3.0808 𝜙𝐾=10.0679 Figure D.8-Average flow field visualization, using PLIF, for Reynolds 81
D Average flow field visualization for different viscosity ratios 322 D.3 Reynolds 118 𝜙𝐾= 0.0097 𝜙𝐾= 0.0277 𝜙𝐾= 0.0572 𝜙𝐾= 0.1052 Figure D.9-Average flow field visualization, using PLIF, for Reynolds 118
D Average flow field visualization for different viscosity ratios 323 𝜙𝐾= 0.2027 𝜙𝐾= 0.2502 𝜙𝐾= 0.2830 𝜙𝐾= 0.4255 Figure D.10-Average flow field visualization, using PLIF, for Reynolds 118
D Average flow field visualization for different viscosity ratios 324 𝜙𝐾= 0.6754 𝜙𝐾= 0.6785 𝜙𝐾= 0.8940 𝜙𝐾= 1.3266 Figure D.11-Average flow field visualization, using PLIF, for Reynolds 118
D Average flow field visualization for different viscosity ratios 325 𝜙𝐾= 1.7358 𝜙𝐾= 2.0207 𝜙𝐾= 1.9620 𝜙𝐾= 2.4638 Figure D.12-Average flow field visualization, using PLIF, for Reynolds 118
D Average flow field visualization for different viscosity ratios 332 𝜙𝐾= 0.5343 𝜙𝐾= 0.6434 𝜙𝐾= 0.6906 𝜙𝐾= 0.8226 Figure D.19-Average flow field visualization, using PLIF, for Reynolds 162
D Average flow field visualization for different viscosity ratios 333 𝜙𝐾= 1.0854 𝜙𝐾= 1.1841 𝜙𝐾= 1.13912 𝜙𝐾= 1.6631 Figure D.20-Average flow field visualization, using PLIF, for Reynolds 162
D Average flow field visualization for different viscosity ratios 334 𝜙𝐾= 1.8038 𝜙𝐾= 1.7639 Figure D.21-Average flow field visualization, using PLIF, for Reynolds 162