1-1 API Gravity (o) v/s sulfur content (%) of major crude oils in 2005, ..... and x=16D downstream of the center of the mixing joint (left to right) 110. 4-30 Yd ..... The effects of the free- ... They neglected the mass transfer between the two phases.
Numerical Simulations of Supercritical Water-Hydrocarbon Mixing in a 3-D Cylindrical Tee Mixer by
Ashwin Raghavan B.Tech., Mechanical Engineering, Indian Institute of Technology Bombay (2011) Submitted to the Department of Mechanical Engineering in partial fulfillment of the requirements for the degree of Master of Science in Mechanical Engineering at the MASSACHUSETTS INSTITUTE OF TECHNOLOGY Feb 2014 c Massachusetts Institute of Technology 2014. All rights reserved.
Author . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Department of Mechanical Engineering Nov 27th, 2013 Certified by . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ahmed F. Ghoniem Professor, Department of Mechanical Engineering Thesis Supervisor Accepted by . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . David E. Hardt Professor, Department of Mechanical Engineering Graduate Officer, Department of Mechanical Engineering
2
Numerical Simulations of Supercritical Water-Hydrocarbon Mixing in a 3-D Cylindrical Tee Mixer by Ashwin Raghavan Submitted to the Department of Mechanical Engineering on Nov 27th, 2013, in partial fulfillment of the requirements for the degree of Master of Science in Mechanical Engineering
Abstract Supercritical water upgrading and desulfurization (SCWUDS) is a new concept in the oil refining industry wherein, crude oil is mixed with supercritical water in a reactor leading to chemical breakdown of the sulfur containing compounds (desulfurization) and cracking of long chain hydrocarbons to shorter chain compounds closer to commercial fuel components (upgrading). The focus of the present work is the development of a numerical tool to investigate the mixing of water and hydrocarbons under supercritical fully-miscible conditions (water and hydrocarbon forming a single phase) in a realistic 3-D cylindrical tee mixer geometry so as to develop an understanding of the effects of geometry, flow rates and fluid properties on the mixing dynamics which in turn will influence the rate of thermal cracking reactions of hydrocarbons and organosulfur compounds as well as the final product distributions. This work includes a consistent treatment of near-critical thermodynamics and transport property variations of real fluid mixtures. A Peng-Robinson EoS with simple van der Waals mixing rules is employed to model the near-critical thermodynamic behavior, with the mixture binary interaction parameter obtained from a Predictive Peng-Robinson approach using a group contribution method (PPR78). A 2nd order accurate finite-volume methodology is used for the numerical solution of the conservation equations. The developed numerical tool was used to investigate the mixing of supercritical water and a model hydrocarbon (n-decane) in a small-scale cylindrical tee mixer (pipe ID = 2.4mm) under fully miscible conditions. For a Reynolds number at the water inlet of 500 and a ∆T between the two streams of 100K (Tw,in = 800K, Td,in = 700K), the flow downstream of the mixing joint remained laminar. Most of the mixing and heat transport occurs due to the circulating action of a counter-rotating vortex pair (CVP) in the hydrocarbon jet formed due to the reorientation of the vorticity preexisting in the hydrocarbon stream flowing through the vertical pipe. This CVP gets progressively weaker as it is advected downstream, due to vorticity diffusion and species and heat transport is dominated by molecular diffusion over small length scales in the far downstream region. Consequently, the mixing rate plateaus in the far 3
downstream region of the tee mixer. Near-critical property variations were found to have a negligible impact on the flow field and mixing behavior under these conditions. However, for a 300K temperature difference between the two streams (Tw,in = 1000K, Td,in = 700K), the water-HC shear layer becomes unstable and rolls up downstream of 5 diameter lengths from the mixing joint. The onset of instability in the shear layer also triggers the stretching and breakdown of the hydrocarbon jet CVP leading to a significant enhancement in mixing manifested as a jump in the mixing rate and a thickening of the mean mixing layer. However, water n-decane mixing under identical inlet conditions but with constant physical properties, showed a stable shear layer with the flow reaching steady state. For a large ∆T between the streams of 300K, the strong density increase (due to cooling of the water component) and the strong viscosity decrease (due to heating of the n-decane component) leads to a local increase in the Re within the mixing layer, resulting in the instability of the shear layer. SCW n-decane mixing with ∆T = 100K was also simulated for increasingly higher Reynolds numbers up to the transition to turbulence. When the Reynolds number at the water inlet is increased to 700, the shear layer between the water and n-decane streams is found to become unstable near x = 6D downstream of the mixing joint followed by the subsequent roll up of the shear layer. The local increase of Re within the mixing layer due to mechanisms similar to the ∆T = 300K case was found to be the cause of the shear layer instability. At Rew,in = 800 the unsteady small scale flow structures in the mixing layer and the consequent flow field fluctuations due to them are much stronger. The stretching and breakdown of the CVP in this case, is accompanied by stronger streamwise vorticity enhancement resulting in much faster mixing compared to the case of Rew,in = 700. Key words: crude oil upgrading, desulfurization, supercritical water, supercritical fluids, mixing, CFD, mixing tee, intermediate Reynolds number, laminar, shear layer instability, transition. Thesis Supervisor: Ahmed F. Ghoniem Title: Professor, Department of Mechanical Engineering
4
Acknowledgments I would like to sincerely thank my research advisor, Prof. Ahmed Ghoniem. His timely advice and guidance in matters of not only research, but also academic life in general have been priceless. I am really grateful for his efforts in helping me channelize my thought, interests and enthusiasm towards concentrated research effort. His constant encouragement and motivation have helped me push myself to broaden and deepen my knowledge and understanding. I would like to express my greatest gratitude to our industry collaborator and sponsor, Saudi Aramco without whose funding and support, this work would not have been possible. It has indeed been a wonderful experience working with them. I would also like to thank Prof. William Green, Prof. Michael Timko, Dr. Guang Wu and all the past and present members of the Supercritical Water Desulfurization and Upgrading team at MIT, as also our counterparts at Saudi Aramco for their constant feedback, encouragement, support and fruitful collaboration. I would also like to express my heartfelt gratitude to Dr. Jose-Sierra Pallares for his invaluable insights and assistance in the thermodynamic modeling of non-ideal mixtures and the state-of-the-art equations of state. The completion of this work would truly not have been possible without him. Special thanks to Dr. Santosh Shanbhogue for his efforts in maintaining our computing cluster and helping us out with all our computing problems, big and small. I am also grateful to Kushal, Gaurav, Konstantina, Christos and all other members of the Reacting Gas Dynamics Laboratory. It has been a pleasure to work with them all and I have learned a lot from each one in the process. I would also like to mention Pavitra, Kashi, Abhishek, Sudeep, Carl, Sameer, Ujwal, Suvinay, Tapovan, Arun, Akash, Surekha, Abaya, Snegdha, Shreya, Shobhna, Jhanvi, Kat and the endless list of my friends who have made my life at MIT worth cherishing forever. My thanks are also due to Lorraine and Leslie for helping me out with all kinds of administrative issues. The acknowledgment cannot be complete without profusely thanking my parents 5
Shri. T.V. Raghavan and Smt. Renuka Raghavan, who taught me the value of hard work, perseverance and learning and showered upon me their constant blessings and support for all these years. Most of all, I would like to thank my elder brother Karthikeyan Raghavan for his perennial support and wise guidance. He has truly been a great mentor and friend throughout my life. I definitely owe my good understanding of fundamentals in mathematics and science to his clear and lucid teaching during my early days as a student. As such, he undoubtedly deserves a share of the credit in all my current and future research endeavors. Last but not the least, I thank The Brahman (God), the supreme cosmic energy which is the basis of this vast, incredible universe; in understanding the mysteries of which, we dedicate our lives. This work has been sponsored by Saudi Aramco Contract No. 6600023444.
6
Contents Abstract
4
Acknowledgments
6
1 Introduction
23
1.1
Role of water in desulfurization and upgrading . . . . . . . . . . . . .
26
1.2
Dynamics of mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
1.2.1
Reactor geometry and flow Reynolds number . . . . . . . . . .
27
1.2.2
Near-critical thermodynamic and transport property variations
29
1.2.3
Water-hydrocarbon phase equilibrium . . . . . . . . . . . . . .
31
Thesis objectives and scope . . . . . . . . . . . . . . . . . . . . . . .
32
1.3
2 Near-critical thermodynamics and transport property variations 2.1
2.2
2.3
35
Near-critical thermodynamics . . . . . . . . . . . . . . . . . . . . . .
37
2.1.1
Peng-Robinson Equation of State for real fluid mixtures . . . .
39
Near-critical transport properties . . . . . . . . . . . . . . . . . . . .
49
2.2.1
Viscosity and Thermal Conductivity . . . . . . . . . . . . . .
49
2.2.2
Mass Diffusivity . . . . . . . . . . . . . . . . . . . . . . . . . .
52
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
3 Problem Formulation and Methodology
57
3.1
Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
3.2
Geometry and Boundary Conditions . . . . . . . . . . . . . . . . . .
60
3.3
Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
7
3.4
3.5
Numerical methodology . . . . . . . . . . . . . . . . . . . . . . . . .
62
3.4.1
Discretization Schemes . . . . . . . . . . . . . . . . . . . . . .
63
3.4.2
Discretized Equations . . . . . . . . . . . . . . . . . . . . . . .
63
3.4.3
Solution Algorithm . . . . . . . . . . . . . . . . . . . . . . . .
65
3.4.4
Numerical Stability . . . . . . . . . . . . . . . . . . . . . . . .
68
3.4.5
Linear System Solvers . . . . . . . . . . . . . . . . . . . . . .
69
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
4 Mixing at intermediate Re: Flow dynamics and impact of temperature difference
71
4.1
Validation: Grid convergence tests . . . . . . . . . . . . . . . . . . . .
73
4.2
Case I: Water n-decane mixing with small temperature difference (Rew,in = 500, Tw,in = 800K, Td,in = 700K) . . . . . . . . . . . . . . .
4.3
75
Case II: Water n-decane mixing with large temperature difference (Rew,in = 500, Tw,in = 1000K, Td,in = 700K) . . . . . . . . . . . . . .
92
4.4
Quantification of mixing . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.5
Impact of near-critical property variations: Cause of shear layer instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.6
Dynamics of flow in the cylindrical tee mixer . . . . . . . . . . . . . . 111
4.7
Vorticity dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.7.1
4.8
Species and thermal transport enhancement due to fluctuations 118
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
5 Mixing at intermediate Re: Impact of Re
127
5.1
Validation: Grid convergence tests . . . . . . . . . . . . . . . . . . . . 127
5.2
Mixing Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.3
Quantification of mixing . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.4
Impact of near-critical property variations: Cause of shear layer instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
5.5
Vorticity dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
5.6
Species and thermal transport enhancement due to fluctuations . . . 154 8
5.7
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
6 Summary and Future Work
161
6.1
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
6.2
Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
9
10
List of Figures 1-1 API Gravity (o ) v/s sulfur content (%) of major crude oils in 2005, bubble sizes proportional to 2005 production volumes, [55] . . . . . .
24
1-2 Distribution of proved world crude reserves (1991, 2001 and 2011) [38]
25
1-3 Schematic of the cylindrical tee reactor geometry used at Saudi Aramco, Dhahran, Saudi Arabia . . . . . . . . . . . . . . . . . . . . . . . . . .
28
1-4 Rectangular opposed-flow tee micromixer [9] . . . . . . . . . . . . . .
29
2-1 Water thermodynamic properties: Comparison between ideal-gas equation predictions and NIST data at P = 25MPa (a) ρ (in kg/m3 ) and (b) Cp (in kJ/kg − K) . . . . . . . . . . . . . . . . . . . . . . . . . .
36
2-2 Comparison of density (in kg/m3 ) predictions by the Peng-Robinson EoS and the Ideal Gas EoS with NIST data at P = 25MPa (a) water and (b) n-decane . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
2-3 Comparison of density predictions by the Volume-Translated PengRobinson EoS (VT-PR EoS), the simple Peng-Robinson EoS and the Ideal Gas EoS with NIST data at P = 25 MPa (a) water and (b) n-decane 45 2-4 Comparison of constant pressure specific heat (Cp in kJ/kg − K) predictions by the Peng-Robinson EoS and the Ideal Gas EoS with NIST data at P = 25 MPa (a) water and (b) n-decane . . . . . . . . . . . .
48
2-5 Comparison of dynamic viscosity (µ in P a−s) predictions using Chung’s generalized correlations [6] with NIST data at P=25MPa (a) water and (b) n-decane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
3-1 Cylindrical tee mixer geometry and boundary conditions . . . . . . .
61
11
3-2 Cylindrical tee reactor mesh . . . . . . . . . . . . . . . . . . . . . . .
62
4-1 ρ (in kg/m3 ) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
4-2 µ (in Pa-s) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
4-3 u (in m/s) v/s z (in m) for simulated Case I for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . . . .
76
4-4 umean (in m/s) v/s z (in m) for simulated Case II for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . .
77
4-5 Tmean (in K) v/s z (in m) for simulated Case II for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
4-6 Yd,mean v/s z (in m) for simulated Case II for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . . . .
79
4-7 Contour plots on the centerplane (y=0 plane) at steady-state for simulated Case I of (top) n-decane mass fraction (Yd ) field and (bottom) Temperature (T) field; The white vertical lines indicate the positions x=2D, x=4D, x=8D, x=10D, x=12D, x=14D downstream of the center of the mixing joint from left to right
. . . . . . . . . . . . . . . . . .
81
4-8 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case I of (a) n-decane mass fraction (Yd ) field and (b) temperature (T in K) field at different downstream locations (x=2D, x=4D, x=6D, x=8D) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
82
4-9 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case I of (a) velocity magnitude (Umag in m/s) (b) density (ρ in kg/m3 ) (c) dynamic viscosity (µ in Pa-s) and (d) Reynolds number (Re = ρUmag D/µ) at different downstream locations (x=2D, x=4D, x=6D, x=8D) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
4-10 (top) Mass fraction of n-decane (Yd ) contours ( middle) Temperature (T in K) contours and (bottom) Streamwise vorticity (ωx in s−1 ) contours for simulated Case I (Rew,in = 500, Tw,in = 800K, Td,in = 700K) at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
84
4-11 Contour plots on the centerplane (y=0 plane) at steady-state for the initially isothermal simulation: (top) Temperature (T) field (in K) ∂T (middle) Rate of heating/cooling due to heat diffusion ρC1 p ∂x∂ j λ ∂x j (in K/s) (bottom) Rate of heating/cooling due to species diffusion P k ∂hk (in K/s); The white along partial enthalpy gradient ρC1 p k ρDk ∂Y ∂xj ∂xj vertical lines indicate the positions x=2D, x=4D, x=6D downstream of the center of the mixing joint from left to right . . . . . . . . . . .
88
4-12 partial enthalpy of water, hw and n-decane, hd (in kJ/kg) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) . . . . .
89
4-13 Contours of the rate of fluid heating/cooling (in K/s) due to (a) heat ∂T diffusion ( ρC1 p ∂x∂ j λ ∂x ) and (b) species diffusion along partial enj P k ∂hk thalpy gradient ( ρC1 p k ρDk ∂Y ) at different downstream locations ∂xj ∂xj (x=2D, x=4D, x=8D, x=16D, x=22D from left to right) . . . . . . . 90 ∂T 4-14 rate of fluid heating/cooling (in K/s) due to (a) heat diffusion ( ρC1 p ∂x∂ j λ ∂x ) j P k ∂hk and (b) species diffusion along partial enthalpy gradient ( ρC1 p k ρDk ∂Y ) ∂xj ∂xj along the vertical centerline at x = 4D for simulated Case I
. . . . .
91
4-15 Lewis number, Le = α/D at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D) . . . . . . . . . . . . . . . . . . . . . .
91
4-16 Fourier transform of the Yd temporal variation at x = 6D, y = 0, z = 0 93 13
4-17 Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field at (a) t = 2.0s (2.5 tf low−through ) (b) t = 2.4s (3 tf low−through ) (c) t = 2.8s (3.5 tf low−through ) (d) t = 3.2s (4 tf low−through ) and (e) Mean field; The white vertical lines indicate the positions x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right) . . . . . . . . . . . . . . . . . . . . . . . .
95
4-18 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case II of (a) mean n-decane mass fraction (Yd ) field and (b) mean temperature (T in K) field at different downstream locations (x=2D, x=4D, x=6D, x=8D) . . . . . . . . . . . . . . . . . . . . . .
96
4-19 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case II of (a) mean velocity magnitude (Umean,mag in m/s) (b) mean density (ρmean in kg/m3 ) (c) mean dynamic viscosity (µmean in Pa-s) and (d) mean Reynolds number (Remean = ρmean Umean,mag D/µmean ) at different downstream locations (x=2D, x=4D, x=6D, x=8D) . . .
97
4-20 Contours of Yd for simulated Case II at different downstream crosssections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) at (a) t = 2.4s (3 tf low−through ) (b) t = 2.8s (3.5 tf low−through ) (c) t = 3.2s (4 tf low−through ) and (d) Mean . . . . . . . . . . . . . . . . . . . . . . . .
98
4-21 Streamwise vorticity (ωx in s−1 ) contours for simulated Case II at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) at (a) t = 2.4s (3 tf low−through ) (b) t = 2.8s (3.5 tf low−through ) (c) t = 3.2s (4 tf low−through ) and (d) Mean . . . . . . . .
99
4-22 Contour plots on the centerplane (y=0 plane) of the streamwise velocity (u in m/s) field at (a) t = 2.0s (2.5 tf low−through ) (b) t = 2.4s (3 tf low−through ) (c) t = 2.8s (3.5 tf low−through ) (d) t = 3.2s (4 tf low−through ) and (e) Mean field; The white vertical lines indicate the positions X=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right) . . . . . . . . . . . . . 101 14
4-23 Contour plots on the centerplane (y=0 plane) of the vertical velocity (w in m/s) field at (a) t = 2.0s (2.5 tf low−through ) (b) t = 2.4s (3 tf low−through ) (c) t = 2.8s (3.5 tf low−through ) (d) t = 3.2s (4 tf low−through ) and (e) Mean field; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right) . . . . . . . . . . . . . 102 4-24 φ v/s x/D: Variation along the length of the tee of the (a) the species mixing quality (φspecies ) for simulated Cases I-IV and (b) the thermal mixing quality (φthermal ) for simulated Cases I and II . . . . . . . . . 104 4-25 Streamwise vorticity (ωx in s−1 ) contours for simulated Case I (top) and Case III (bottom) at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D) . . . . . . . . . . . . . . . . . . . . . . 106 4-26 n-decane mass fraction (Yd ) contours for simulated Case I (top) and Case III (bottom) at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D) . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4-27 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases I and III of (a) density (ρ in kg/m3 ) and (b) dynamic viscosity (µ in Pa-s) at x=5D . . . . . . . . . . . . . . . . . . . . . . 108 4-28 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases I and III of the local Reynolds number at x=5D . . . 109 4-29 Yd contours on the centerplane for simulated Case IV; The white vertical lines indicate the positions x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right) 110 4-30 Yd contours for simulated Case IV at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) . . . . . . . . 111 4-31 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases II and IV of (a) density (ρ in kg/m3 ) and (b) dynamic viscosity (µ in Pa-s) at x=5D . . . . . . . . . . . . . . . . . . . . . . 112 4-32 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases II and IV of the local Reynolds number at x=5D . . 113 15
4-33 Streamwise vorticity (ωx in s−1 ) contours for simulated Case I (Rew,in = 500, Tw,in = 800K, Td,in = 700K) at different downstream cross-sections (x=-0.5D, x=-0.25D, x=0, x=0.25D, x=0.5D, x=1D, x=2D) . . . . . 115 4-34 Transverse velocity vectors (velocity component along constant X planes, √ utrans = v 2 + w2 ) for simulated Case I (Rew,in = 500, Tw,in = 800K, Td,in = 700K) at two downstream cross-sections (x=0.25D, x=2D) Note: The length of the vectors is not representative of the magnitude in the figures, the color of the vectors represents the magnitude . . . . . . . 116 4-35 (left) Contour plots on the centerplane (y=0 plane) of the velocity 1/2
magnitude (Umag = (u2 + v 2 + w2 )
) in the joint region; (right) n-
decane mass fraction profile along the vertical line at x = 2D (the position indicated by the white line in the left figure) . . . . . . . . . 117 4-36 Contours of the mean field for simulated Case II at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) of (a) ω˙ st,x in s−2 (b) ω˙ bg,x in s−2 (c) ω˙ v,x in s−2 and (d) ω˙ d,x in s−2 119 4-37 Profiles along the vertical line at y=0.0004 (y=D/6) for simulated Case II of (a) ω˙ st,x,mean (in s−2 ) and (b) ω˙ bg,x,mean (in s−2 ) at different downstream locations (x=4D, x=5D, x=6D, x=8D) . . . . . . . . . . . . . 120 4-38 Profiles along the vertical line at y=0.0004 (y=D/6) for simulated Case II of ωx,mean (in s−1 ) at different downstream locations (x=4D, x=5D, x=6D, x=8D) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 0 1/2 4-39 RMS fluctuation of Yd (Y¯d 2 ) contours for simulated Case II at dif-
ferent downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4-40 Contours for simulated Case II at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=16D from left to right) of (a) 0 wmean Yd,mean /uref (b) w0¯Yd /uref ; uref = 0.14 m/s . . . . . . . . . . . 123 4-41 Contours for simulated Case II at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=16D from left to right) of (a) 0 vmean Yd,mean /uref (b) v 0¯Yd /uref ; uref = 0.14 m/s . . . . . . . . . . . . 124 16
5-1 u (in m/s) v/s z (in m) for simulated Case IV for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5-2 Yd v/s z (in m) for simulated Case IV for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5-3 T (in K) v/s z (in m) for simulated Case IV for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5-4 Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field at steady-state for (a) Case I and (b) Case II; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (from left to right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5-5 Contour plots of Yd on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=16D and x=22D from left to right) at steady-state for (a) Case I and (b) Case II . . . . . . . . . . . . . . 133 5-6 Contour plots of T (in K) on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=16D and x=22D from left to right) at steady-state for (a) Case I and (b) Case II . . . . . . . . 134 5-7 Contour plots of the streamwise vorticity, ωx (in s−1 ) on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=16D and x=22D from left to right) at steady-state for (a) Case I and (b) Case II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5-8 Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field for simulated Case III at (a) t = 2.0s (4 tf low−through ) (b) t = 2.2s (4.4 tf low−through ) (c) t = 2.5s (5 tf low−through ) and (d) Mean field; The white vertical lines indicate the positions x=6D, x=8D, x=10D, x=12D, x=14D, x=16D, x=18D and x=20D downstream of the center of the mixing joint (left to right) . . . . . . . . . . . . . . . . . . . . . 137 17
5-9 Fourier transform of the Yd temporal variation at x = 8D, y = 0, z = 0 for Case III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 5-10 Contours of Yd for simulated Case III at different downstream crosssections (x=2D, x=4D, x=8D, x=12D, x=16D from left to right) at (a) t = 2.0s (b) t = 2.2s (c) t = 2.5s and (d) Mean . . . . . . . . . . 139 5-11 Contours of ωx (in s−1 ) for simulated Case III at different downstream cross-sections (x=2D, x=4D, x=8D, x=12D, x=16D from left to right) at (a) t = 2.0s (b) t = 2.2s (c) t = 2.5s and (d) Mean . . . . . . . . . 140 0
5-12 Contour plots on the centerplane (y=0 plane) for Case III of u /U , where U is the reference velcoity used for normalization taken to be the water inlet average velocity; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D, x=14D and x=16D downstream of the center of the mixing joint (left to right) . . 141 5-13 Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field for simulated Case IV at (a) t = 1.8s (4 tf low−through ) (b) t = 2.05s (≈ 4.5 tf low−through ) (c) t = 2.25s (5 tf low−through ) and (d) Mean field; The white vertical lines indicate the positions x=4D, x=6D, x=8D, x=10D, x=12D, x=14D, x=16D and x=18D downstream of the center of the mixing joint (left to right) . . . . . . . . . . . . . . . . . 142 5-14 Contours of Yd for simulated Case IV at different downstream crosssections (x=4D, x=6D, x=8D, x=12D, x=16D from left to right) at (a) t = 1.8s (b) t = 2.05s (c) t = 2.25s and (d) Mean . . . . . . . . . 143 5-15 Contours of ωx (in s−1 ) for simulated Case IV at different downstream cross-sections (x=4D, x=6D, x=8D, x=12D, x=16D from left to right) at (a) t = 1.8s (b) t = 2.05s (c) t = 2.25s and (d) Mean
. . . . . . . 144
5-16 φ v/s x/D: Variation along the length of the tee of the (a) the species mixing quality (φspecies ) and (b) the thermal mixing quality (φthermal ) for simulated Cases I-IV . . . . . . . . . . . . . . . . . . . . . . . . . 146 18
5-17 ρ (in kg/m3 ) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 5-18 µ (in Pa-s) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 5-19 Yd contours on the centerplane for simulated Case V; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=14D downstream of the center of the mixing joint (left to right) 149 5-20 Yd contours for simulated Case V at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D from left to right) . . . . . . . 150 5-21 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases III and V of (a) density (ρ in kg/m3 ) and (b) dynamic viscosity (µ in Pa-s) at x=6D . . . . . . . . . . . . . . . . . . . . . . 152 5-22 Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases III and V of the local Reynolds number at x=6D . . 153 5-23 Contours of the mean field for simulated Case IV at different downstream cross-sections (x=4D, x=6D, x=8D, x=10D, x=16D from left to right) of (a) ω˙ st,y in s−2 (b) ω˙ bg,y in s−2 (c) ω˙ v,y in s−2 and (d) ω˙ d,y in s−2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 0 1/2 5-24 RMS fluctuation of Yd (Y¯d 2 ) contours for simulated Case IV at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 5-25 Contours for simulated Case IV at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=10D from left to right) of 0 (a) wmean Yd,mean /uref (b) w0¯Yd /uref ; uref = 0.13 m/s . . . . . . . . . 157 5-26 Contours for simulated Case IV at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=10D from left to right) of 0 (a) vmean Yd,mean /uref (b) v 0¯Yd /uref ; uref = 0.13 m/s . . . . . . . . . . 158
19
20
List of Tables 2.1
Molecular weight, critical properties and acentric factor for water and n-decane [32] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2
Group interaction parameters for the PPR78 model, Akl = Alk and Bkl = Blk (in MPa) [14][39]
2.3
. . . . . . . . . . . . . . . . . . . . . . .
51
Function coefficients used for generalized thermal conductivity correlations from Chung et al.[6] . . . . . . . . . . . . . . . . . . . . . . .
2.7
49
Function coefficients used for generalized viscosity correlations from Chung et al.[6] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6
44
Acentric factor, dipole moment and association factor for water and n-decane [32] [8] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5
42
Molar critical volume and volume translation parameter for water and n-decane Mathias et al. [27] . . . . . . . . . . . . . . . . . . . . . . .
2.4
40
54
Lennard-Jones force constants for water and n-decane Liu et al. (1998) [24] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
4.1
Inlet flow conditions for the cases discussed in this chapter . . . . . .
72
4.2
Properties of water and n-decane at 25MPa and extreme temperatures of the simulations I and II . . . . . . . . . . . . . . . . . . . . . . . .
74
5.1
Inlet flow conditions for the cases simulated in this study . . . . . . . 128
5.2
Properties of water and n-decane at 25MPa and extreme temperatures of the simulations in this study . . . . . . . . . . . . . . . . . . . . . 147
5.3
Inlet flow conditions for the comparison test Case V simulated in this study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 21
22
Chapter 1 Introduction The crude oil obtained, on an average, from oil wells worldwide continues to become sourer (higher sulfur mass fraction) and heavier (higher heavy hydrocarbon fraction). Fig.1-1 shows the API gravity (in o , see Eq.1.1) v/s sulfur content (%) scatter for the most widely produced crudes in the world in 2005 illustrating how the Arab, Ural and Iran crudes which constitute the bulk of world oil production are of the middle to heavy and sour grade. Moreover, close to 75% of proved oil reserves in the world exist in the Middle-East, Ural region and South America (mainly Venezuela) most of which is heavy, sour crude oil (see Fig.1-2). On the other hand, regulatory standards on fuels keep getting tighter with the maximum allowable level of sulfur set by the US EPA as 15 ppm for diesel [51] and 30 ppm for gasoline [50]. Since a majority of the crude oil demand comprises of light (for gasoline) and middle (for diesel and jet fuel) distillates, the increasing heavy fraction of the crude oil presents a serious demand-supply problem. All this, coupled with the perpetual effort of oil companies to improve their margins necessitates the development of new efficient cost-effective technologies for the desulfurization (removal of sulfur compounds) and upgrading (increasing the light fraction) of crude oil feedstock. AP Igravity =
141.5 − 131.5 ρoil /ρwater
(1.1)
One such new concept in the oil industry is supercritical water upgrading and 23
Figure 1-1: API Gravity (o ) v/s sulfur content (%) of major crude oils in 2005, bubble sizes proportional to 2005 production volumes, [55]
desulfurization(SCWUDS) of crude oil wherein, crude oil is mixed with water at supercritical pressure and temperature (Pcritical = 22.1M P a, Tcritical = 647K), causing the sulfur compounds to breakdown, releasing hydrogen sulfide (H2 S) in the process, while also resulting in the cracking of long chain heavy hydrocarbons into smaller chains (which are components of conventional fuels). This technology offers economic benefits over the conventional desulfurization methodology in which transition metal catalysts and hydrogen gas are used (hydro-desulfurization). In order to develop SCWUDS into a viable alternative to hydro-desulfurization followed by upgrading in the oil industry, one first needs to understand clearly the physics and chemistry governing the process and the interactions of the two therein. The thermal breakdown reactions of organic sulfur compounds and heavy hydrocarbons are strongly depen24
Figure 1-2: Distribution of proved world crude reserves (1991, 2001 and 2011) [38]
dent on the local temperature and concentration of the respective species. Therefore, understanding and characterizing the mixing dynamics in a realistic SCWUDS reactor geometry is crucial to be able to estimate conversion rates and product distributions of the SCWUDS process. The focus of this work and ongoing work is to develop a computational tool capable of simulating mixing of supercritical water and hydrocarbons at a range of temperature and pressure conditions pertinent to the SCWUDS process so as to develop insight into the physical mechanisms governing the process. Consequently, we would be in a position to tweak and optimize the design of the reactor and/or the process as a whole to make it more effective and efficient. 25
1.1
Role of water in desulfurization and upgrading
When crude oil is heated up to sufficiently high temperatures (typically > 400o C), the sulfur containing compounds and long hydrocarbon chains breakdown due to thermal cracking reactions. This typically results in the formation of shorter hydrocarbon chains, increasing the light fraction. Unfortunately, the cracking reactions also give rise to coke precursor radicals like olefinic radicals which can further undergo cyclization and condensation reactions to form coke (carbonaceous agglomerates insoluble even in toluene). Coke formation not only eats into potential light hydrocarbon yields but also leads to carbonaceous deposits on the reactor walls thus clogging up the reactor. The presence of water is known to suppress coke formation by inhibiting the cyclization and condensation reactions through the solvation and dispersion of coke precursors. The solvation and dispersion effect of water prevents coke precursors from coming together in the medium, thus suppressing coke formation pathways. Water may also suppress coke formation by capping the coke precursor radicals by acting as a hydrogen donor. However, this role of water as a hydrogen donor under supercritical conditions is disputed in the literature. The role of water in suppressing coke formation in the supercritical water upgrading of vacuum residue and bitumen is well illustrated in the works of Cheng et al. [5] and Vilcaez et al. [53] respectively. Thus water transports heat to the hydrocarbons in a controlled fashion at the same time inhibiting coke formation pathways. The downside to the presence of water is the reduction in the rate of hydrocarbon cracking reactions due to dilution of the species. Whether water plays an active chemical role in the primary breakdown processes of large hydrocarbons and sulfur compounds is still uncertain and a topic of continuing research.
1.2
Dynamics of mixing
Since the rates of the different cracking reactions and coke formation reactions are strongly dependent on the local temperature and species concentrations, it is im26
portant to understand the transport of heat and water to the hydrocarbons through mixing. The mixing of supercritical water with hydrocarbons is influenced by a variety of factors which are listed below: • Reactor geometry and flow Reynolds number • Near-critical thermodynamic and transport property variations • Water-hydrocarbon phase equilibrium
1.2.1
Reactor geometry and flow Reynolds number
The present study was focused on a cylindrical tee reactor geometry shown in Fig.13. Flow in a cylindrical tee mixer exhibits a variety of flow features likely to be encountered in any general reactor. Also, this reactor geometry is currently being used by our collaborators at Saudi Aramco for their tests on SCWDS and upgrading (small scale test reactor, tube ID = 2.4 mm). As such, it makes sense for us to begin our modeling efforts with this reactor geometry and study the flow structures and their effect on the mixing dynamics. The bulk flow Reynolds number (Re) which depends on the flow rates, densities and viscosities of both the water and hydrocarbon streams will determine whether the flow downstream of the mixing joint will remain laminar or undergo transition to turbulence. Dreher et al. [9] performed CFD simulations of mixing of two streams in a rectangular cross-section opposed-flow tee micromixer with channel width of 300µm (see Fig.1-4) for Re up to 1000. They were able to observe the transition of the flow regime from straight laminar flow for Re < 10 to steady engulfment flow for 10 < Re < 240 (the two fluid streams swap to the opposite side) to periodically fluctuating flow for 240 < Re < 500 and finally to chaotic flow for Re < 500 (turbulence). Hoffmann et al. [12] experimentally characterized the flow and mixing in a rectangular opposedflow tee micromixer and also observed the flow transition from the straight laminar to engulfment flow regime. Correia et al. [7] numerically investigated the mixing in opposed-flow micromixers with a tortuous downstream flow channel and found an unsymmetrical downstream channel configuration results in slightly better mixing than 27
Inlet 2 Cold pressurized (25MPa) crude oil feed Inlet 1 Hot SC water (25MPa)
3D
3D
22D
Outlet
Figure 1-3: Schematic of the cylindrical tee reactor geometry used at Saudi Aramco, Dhahran, Saudi Arabia
a symmetric configuration. Though these studies focus on the laminar to transition flow regimes which are of interest to us in the preliminary small-scale tests of SCWDS, they study mixing of two streams of the same fluid in rectangular opposed-flow tees. Mixing of hot and cold water in cylindrical mixing tees (similar geometry to ours) at high Re has been extensively studied in the nuclear industry community using the Large Eddy Simulation (LES) and Unsteady Reynolds Averaged Navier-Stokes (URANS) methodologies. LESs were performed by Kuczaj et al. [20, 19], Jayaraju et al. [15], Westin et al. [54] and Odemark et al. [33] for very high Re (∼150,000). They were all able to obtain satisfactory agreement with the experimental results of Andersson et al. [2] using appropriate mesh resolutions, inlet boundary conditions and sub-grid scale models. Frank et al. [11] and Merzari et al. [28] performed URANS 28
of the same flow but were unable to capture the transient behavior as well as the LES calculations. The focus of all these studies was mainly to validate LES and URANS as viable predictive techniques for the prediction of thermal fluctuations in tee junctions in the nuclear industry. As such, they did not investigate the laminar to turbulent transition in cylindrical tee mixers.
Figure 1-4: Rectangular opposed-flow tee micromixer [9]
1.2.2
Near-critical thermodynamic and transport property variations
The peculiar variations of thermodynamic and transport properties near the critical point can be expected to affect the flow and mixing dynamics. Chapter 2 includes a detailed discussion of real fluid behavior near the fluid critical point and the near29
critical variation of thermodynamic (density, specific heat) and transport (viscosity, thermal conductivity, mass diffusivity) properties for both pure real fluids and nonideal mixtures of real fluids. Cubic equations of state which can capture real-fluid thermodynamics have been formulated by Peng at al.[37] (Peng-Robinson EoS) and Soave et al.[47] (Redlich-Kwong-Soave EoS). These equations of state are the most widely used in numerical simulations of real-fluid flows due to their simplicity of implementation and effectiveness in capturing near-critical variations in fluid densities and specific heats (especially for hydrocarbons). The variations of viscosity and thermal conductivity of fluids and fluid mixtures near the critical point can be modeled using the correlations developed by Chung et al.[6] for a large number of substances and a wide range of temperature-pressure conditions. Liu et al.[23][24] and Silva et al.[46] have developed expressions for the variation of mass diffusivities with temperature and pressure of substances in dense non-ideal mixtures based on the hard-sphere Lennard-Jones model. Miller et al. performed Direct Numerical Simulations (DNS) of a supercritical heptane-nitrogen temporal mixing layer (both two and three dimensional) [29] using a Peng-Robinson EoS to model the near-critical thermodynamic behavior and studied the effect of density stratification across the shear layer on the transition to turbulence. Okong’o et al. performed a similar DNS of a supercritical liquid oxygen/hydrogen three dimensional mixing layer [34] and observed a similar effect of the density stratification causing a suppression of the transition to turbulence. They investigated this effect by looking at the vorticity budget. The effects of the freestream density ratio on the evolution of incompressible, high Reynolds and Froude number, confined mixing layer was investigated numerically by Soteriou et al. [48]. They found that a non-unity density ratio alters the flow characteristics significantly, influencing the entrainment patterns by means of baroclinic generation of vorticity. Zong et al. [58] have performed the LES of a two-dimensional cryogenic supercritical nitrogen jet at high Re and studied the spatial and temporal evolution of the jet for different ambient pressures which lead to varying amounts of density stratification. Their findings seem consistent with previous results of Miller et al. [29] and Okong’o 30
et al. [34] regarding the effect of density stratification leading to the suppression of shear layer transition to turbulence. A similar study of the dynamics of cryogenic nitrogen jets at supercritical pressures was performed by Kim et al. [17] using a RANS framework with a k- turbulence model and a presumed probability-density function for the conserved scalars. Park [36] performed a series of RANS and LES simulations of a cryogenic liquid nitrogen jet with a variety of turbulence models and three different equations of state (ideal gas, P-R and SRK) and found that the choice of an appropriate real-fluid EoS is more important in capturing the flow and mixing dynamics than the choice of turbulence models. Schmitt et al. [44] analyzed the behavior of a nitrogen coaxial shear jet under supercritical pressures through LES and experiments and investigated the effect of externally imposed acoustic perturbations on the jet mixing efficiency. Narayanan et al. [31] simulated a supercritical water oxidation (SCWO) hydrothermal flame of methanol using RANS and an eddy dissipation combustion model with single-step kinetics. They were successful in obtaining a fair agreement with their experiments for the flame position but over-predicted the flame temperature. SierraPallares et al. [45] have also simulated hydrothermal flames in SCWO of methanol using an eddy dissipation concept along with a micro-mixing model and were able to predict flame structure and temperature to a fair degree of accuracy. Kim et al. [18] validated the use of a real-fluid flamelet model in the simulation of gaseous hydrogen/cryogenic oxygen coaxial jet flames under supercritical pressures. All of the aforementioned numerical studies provide good insight into the flow and mixing behavior of near-critical fluids and were very helpful in formulating the methodology of the study of 3-D SCW-HC mixing in a tee reactor performed in the present work.
1.2.3
Water-hydrocarbon phase equilibrium
The effect of water-hydrocarbon phase equilibrium on the rate of mixing of the two fluids is captured very well by Dabiri et al. [8] in fundamental numerical studies of the mixing process of a cold hydrocarbon droplet in a vast reservoir of heated supercritical water. Their studies of water-toluene and water-decane mixing highlight 31
the importance of the upper-critical solution temperature (UCST) of the two fluids in determining the persistence of the phase interface which greatly retards the mixing rate. They also observed large abrupt increases in mixing rate for the water-toluene system as the water temperature is increased from below the UCST to above the UCST due to the flat nature of the water-toluene phase equilibrium curve. Guang et al. [56] extended this numerical framework to ternary water-HC-HC systems of water-toluene-decane and water-toluene-tetralin. They demonstrated an interesting fractionation effect (selective buildup of one of the hydrocarbons in the water phase) when the water temperature is between the UCSTs of the two different water-HC binary systems (UCST of water-toluene system is 308o C while UCST of water-decane system is 359o C). Erriguible et al. [10] simulated jet breakup in a supercritical antisolvent process (SAS) where a is solution injected into CO2 under conditions below the mixture critical point (liquid and gaseous phases coexist) using a Volume of Fluid (VoF) one-fluid method. They neglected the mass transfer between the two phases across the interface and were able to correlate the jet breakup length to the Weber number of the liquid phase. However, numerical studies of multiphase flow at nearcritical conditions in general 3-D geometries are not reported in literature. This thesis focuses on the mixing of water and hydrocarbon under fully miscible conditions (above the UCST). Development of a robust numerical method to study multiphase mixing of near-critical fluids in generic 3-D geometries using sophisticated interface-tracking methods (like VoF and level-set) coupled with mass and heat transfer across the phase interface is the focus of ongoing and future work in this research project.
1.3
Thesis objectives and scope
The work described in this thesis is the first major step in our efforts to ultimately simulate the reactive mixing process of supercritical water, multiple hydrocarbons and organosulfur compounds in the SCWUDS of crude oil. The focus of the present work was to study the mixing of supercritical water and a model hydrocarbon (chosen to be n-decane in this study) in a 3-D cylindrical tee reactor geometry under fully miscible 32
conditions (cold hydrocarbon temperature above UCST of water-decane system) at flow Re up to transition to turbulence. The following was accomplished under the scope of this work: • Development of a robust numerical code to perform 3-D near-critical fluid mixing simulations using existing open-source CFD libraries (OpenFOAM [35]) with a consistent treatment of near-critical fluid thermodynamics and transport property variations (completely new additions to the libraries) • Characterization of the mixing and flow structures in a cylindrical tee reactor at intermediate Reynolds numbers under supercritical fully-miscible conditions • Investigation of the impact of the temperature difference between the streams on the flow dynamics and mixing behavior under fully miscible conditions • Investigation of the impact of the Reynolds number on the flow dynamics and mixing behavior under fully miscible conditions up to the transition to turbulence The thesis is organized as follows. Chapter 2 includes a detailed discussion of real-fluid thermodynamics and transport property variations near the critical point. This is followed by a detailed discussion of the problem formulation and numerical methodology employed including the governing equations, geometry, meshing, boundary conditions, flow conditions and numerical methods in Chapter 3. Chapter 4 presents and discusses results of the simulation of water and n-decane mixing at a water inlet Re of 500 under fully miscible conditions and discusses the flow features and mixing dynamics in a cylindrical tee reactor under laminar flow conditions. The impact of the temperature difference between the streams on the stability of the water-HC shear layer and consequently, the mixing dynamics is also discussed in detail in this chapter. Chapter 5 presents mixing results for different inlet Re values in the range of 500-800 and evaluates the impact of Re on the flow and mixing dynamics. The Re at which the transition to turbulence occurs is determined and the physical mechanisms leading to the shear layer instability are explained. Finally, Chapter 6 33
summarizes the work done and discusses the future goals and plans for the research project.
34
Chapter 2 Near-critical thermodynamics and transport property variations
While simulating the flow and mixing of fluids, it is crucial to be able to model the thermodynamic (density, specific heat, enthalpy) and transport (viscosity, thermal conductivity, mass diffusivity) properties at the relevant temperature and pressure conditions since they directly affect the different heat and mass transport mechanisms. The distinctive feature of flows of near-critical fluids is the peculiar variation of these fundamental properties near the critical point as well as the non-linear interactions between the different components in a fluid mixture. Fig.2-1 clearly illustrates the inability of simple equations of state like the ideal-gas EoS to capture these variations in thermodynamic fluid properties (density, specific heat). Also, simplified calculations of transport properties based on the kinetic theory of gases is inadequate. Hence, there is a need to understand the near-critical behavior of pure fluids and their mixtures and select appropriate equations of state and transport models to be used in the numerical simulations of the problem at hand. This chapter gives a brief overview of popular methods of modeling near-critical fluid behavior and discusses the details of the specific equation of state and transport models selected for our particular application. 35
(a) ρ
(b) Cp
Figure 2-1: Water thermodynamic properties: Comparison between ideal-gas equation predictions and NIST data at P = 25MPa (a) ρ (in kg/m3 ) and (b) Cp (in kJ/kg − K) 36
2.1
Near-critical thermodynamics
As a fluid is pressurized at a given temperature, the distinction between the liquid and vapor phases disappears at a particular pressure. This is the critical pressure of the fluid. Similarly, the temperature above which there is no distinction between the liquid and vapor phases is the critical temperature of the fluid. So if either the pressure or the temperature is above the critical value, the fluid is said to be in a supercritical state and does not undergo a phase transition. However, even at supercritical pressures the properties of the fluid change rapidly in the vicinity of the critical temperature (rapid drop in density, huge spike in specific heat). Simplified equations of state like the Ideal Gas EoS cannot capture these variations in thermodynamic properties giving rise to the need for state equations that model more complex behavior of fluid molecules in order to predict the correct property variation trends. The earliest such equation of state that tried to predict the thermodynamic properties of a fluid over a wide range of temperatures and pressures was developed by van der Waals (1873) [52]. He included the effects of intermolecular forces and the physical volume occupied by the fluid molecules in developing a cubic equation of state given by Eq.2.1.
P =
Ru T a − 2 V −b V
(2.1)
where P is the pressure, T is the temperature, V is the molar specific volume and Ru is the universal gas constant (8.314 J/mol-K). a and b are constants for a given substance with a representing intermolecular interactions and b representing the physical volume excluded by molecules of the fluid (which are not negligible at such high pressures). Though this equation correctly represents qualitative trends in density, it is quantitatively inaccurate with the PVT predictions getting considerably worse towards higher densities of the fluid. Moreover, it does not yield accurate estimates of derived thermodynamics properties like enthalpy and entropy. Also, since the EoS involves only two constants for a given substance (both of which are 37
determined to closely satisfy the stability criteria at the critical point), a universal critical compressibility factor of Zc = 0.375 (see Eq. 2.2) is predicted for all fluids.
Zc =
Pc Vc Ru Tc
(2.2)
where, Pc , Tc and Vc are the critical pressure, critical temperature and critical specific molar volume respectively. In order to overcome these drawbacks of the van der Waals EoS, many researchers have proposed complicated equations of state which are typically formulated in a virial format with exponential terms for the high-density regime. These include the state equations proposed by Benedict et al. [3], Starling and Han [49], Martin and Stanford [25] and Lee and Kesler [22]. These equations of state are cumbersome to use in a CFD calculation and work well only for a specific set of compounds for which they were developed. A good critical review of these methods can be found in the dissertation of Kutney [21]. Improved easy-to-use cubic equations of state were later formulated by Redlich and Kwong [40] (Eq.2.3), Soave [47] (Eq.2.4) and Peng and Robinson [37] (Eq.2.6). Redlich-Kwong EoS (1949) [40]: P =
Ru T a − 1/2 V − b T V (V + b)
(2.3)
Redlich-Kwong-Soave EoS (1972) [47]: P =
Ru T a − V − b V (V + b)
(2.4)
where, a is not a constant but a temperature dependent function (see Eq.2.5) with the acentric factor (ω) as a parameter. The acentric factor for a substance represents how much its behavior deviates from that of a monoatomic gas. a = ac α(ω, Tr )
(2.5)
where ac is the intermolecular force constant at the critical point and Tr is the reduced 38
temperature (T /Tc ). Peng-Robinson EoS (1976) [37]: P =
a Ru T − V − b V (V + b) + b(V − b)
(2.6)
with a similar temperature dependence for a as in the RKS EoS. In this study, the Peng-Robinson EoS was chosen for its simplicity and particular effectiveness in dealing with the thermodynamic behavior of hydrocarbons and water. Specifically, the extension of the PR EoS to real-fluid mixtures is employed as is explained in detail in the following section.
2.1.1
Peng-Robinson Equation of State for real fluid mixtures
The cubic Peng-Robinson EoS for real fluid mixtures used in conjunction with standard van der Waals mixing rules is given by: P =
Ru T am − 2 V − bm V + 2bm V − b2m
(2.7)
am =
XX
√ Xi Xj (1 − kij ) ai aj
(2.8)
i
j
bm =
X
X i bi
(2.9)
i
where, Xi is the mole fraction of specie i and the parameters ai and bi for each species in the mixture is calculated using its critical temperature (Tci ), critical pressure (Pci ) and acentric factor (ωi ): " r !#2 Ru2 Tci2 T ai = 0.457236 1 + 0.3746 + 1.5423ωi − 0.2699ωi2 1 − Pci Tci 39
(2.10)
substance water n-decane
Mw (kg/kmol) 18.0 142.285
Pc (M P a) Tc (K) 22.06 647 2.123 618.5
Vc (m3 /kmol) ω 0.0571 0.344 0.6031 0.484
Table 2.1: Molecular weight, critical properties and acentric factor for water and n-decane [32]
bi = 0.0777961
Ru Tci Pci
(2.11)
The molecular weight (Mw ), critical temperature (Tc ), critical pressure (Pc ), critical volume (Vc ) and acentric factor (ω) for water and n-decane are given in table 2.1. The molecular weights are required to determine the effective molecular weight of the mixture, which in turn, is required to calculate the mixture density (ρ) from the mixture specific volume (V) using ρ = Mw /V . In Eq.2.8, kij is a binary interaction parameter for the species pair i and j. In this study, the binary interaction parameter used for the water-n-decane system is a temperature dependent quadratic polynomial function as given in Eq.2.12 below, applicable to the range of T from 700K to 1000K at P=25MPa, which encompass the thermodynamic conditions investigated in this study. k12 (T ) = 1.14558E(−06)T 2 − 2.80487E(−03)T + 1.49148
(2.12)
The Predictive Peng-Robinson or PPR78 approach [14] with a group contribution method to capture inter-molecular interactions between the components was employed to determine the T dependent function for the binary interaction parameter. In the PPR78 approach, the same cubic PR EoS is employed with the same linear mixing rule as in Eq.2.9 for the mixture excluded volume bm . However, the mixing rule for the parameter am is formulated on a stronger physical basis as given below in Eq.2.13: " am = b m ∗
X i
ai g E Xi − res bi C
# (2.13)
E where, gres is the residual part of the molar excess Gibbs energy when the pressure
40
goes to infinity and is given by Eq.2.14 below: E 1 gres = C 2
P P i
j
Xi Xj bi bj Eij (T ) bm
(2.14)
where, the constant C in Eq.2.13 and Eq.2.14 above is given by: √ √ 2 C= ln 1 + 2 ≈ 0.6232 2
(2.15)
The interaction parameters Eij (T ) in Eq.2.14 are estimated using a group contribution method (GCM) based on the works of Kehiaian et al. [16] and Abdoul et al. [1], the equation for which is given below in Eq.2.16:
Eij (T ) = −
Ng Ng 1 XX
2
(αik − αjk ) (αil − αjl ) Akl
k=1 l=1
298.15 T
Bkl −1 Akl
(2.16)
where, Ng is the number of different functional groups defined by the method. αik is the fraction of molecule i occupied by group k, that is: αik =
occurrence of group k in molecule i total number of groups in molecule i
(2.17)
The constants Akl = Alk and Bkl = Blk were determined by Jaubert et al. in their previous studies by minimizing the deviations between calculated and experimental vapor-liquid equilibrium (VLE) data for numerous binary systems. The constants Akl and Bkl for the functional groups present in the molecules in this study (−CH3 , −CH2 andH2 O are given (in MPa) in Table2.2. Also, Akk = Bkk = 0. The capability of the PPR78 model in predicting mixture properties and phase equilibria accurately for systems/conditions outside the date-set used for fitting the above parameters has been demonstrated extensively, for example in [14], [39]. Due to its predictive capability, this model was selected for the water-HC system in the present study for which fitting data under the T,P conditions of interest was not available in the literature. The advantage of this formulation is that the binary interaction parameter (kij ) of Eq.2.8 can be easily recast as a function of temperature using Eq.2.13 41
CH3 (group 1) CH2 (group 2) H2 O (group 16) CH2 (group 2)
CH3 (group 1) 0 A12 = 74.81 B12 = 165.7 A1−16 = 3557 B1−16 = 11195
CH2 (group 2) 0
H2 O (group 16) -
A2−16 = 4324 B2−16 = 12126
0
Table 2.2: Group interaction parameters for the PPR78 model, Akl = Alk and Bkl = Blk (in MPa) [14][39]
to Eq.2.16 as below: − 12 kij (T ) =
PNg PNg k=1
l=1
298.15 T
(αik − αjk ) (αil − αjl ) Akl √ ai (T )aj (T ) 2 bi bj
√
√
Bkl −1 Akl
−
ai (T ) bi
−
aj (T ) bj
(2.18) The above equation was used to determine the quadratic variation of kij with T given in Eq.2.12. Fig.2-2 shows the variation of density for water and n-decane with temperature at P = 25MPa (close to the operating pressure of SCWUDS reactor). It compares the density predictions of the PR EoS and the Ideal Gas EoS with experimental data values obtained from NIST [32]. One can see that the PR EoS predicts the density variation near the critical temperature and above very well while the Ideal Gas EoS fails to do so. However, the error in density prediction by the PR EoS keeps increasing towards the high-density (liquid-like) region. This can be corrected using the volume translation approach of Mathias et al. [27] to correct for the molar specific volumes predicted by the PR EoS as below: Vmcorr
= Vm + tm + fc
V2 δ=− m Ru T
∂P ∂V
0.41 0.41 + δ
(2.19)
(2.20) T
fc = Vcm − (3.946bm + tm ) 42
(2.21)
2
(a) water
(b) n-decane
Figure 2-2: Comparison of density (in kg/m3 ) predictions by the Peng-Robinson EoS and the Ideal Gas EoS with NIST data at P = 25MPa (a) water and (b) n-decane 43
substance water n-decane
Vc (cm3 /mol) t(cm3 /mol) 57.1 −3.40 603.1 0.90
Table 2.3: Molar critical volume and volume translation parameter for water and n-decane Mathias et al. [27]
where, the constant 0.41 was determined by Mathias et al. [27] by regressing data for many substances and the volume translation parameter (tm ) and molar critical volume (Vcm ) for the mixture can be determined using simple mixing rules given below: tm =
X
Xi ti
(2.22)
Xi Vci
(2.23)
i
Vcm =
X i
where, the critical molar volumes (Vci ) and volume translation parameter (ti ) for water and n-decane are given in Table 2.3. The volume translation parameter for n-decane was not provided in the work of Mathias et al. [27]. The value for n-decane was roughly tuned to obtain the best fit for predicted densities with NIST data in the range of 300 K to 650 K at P = 25MPa. Fig.2-3 compares the density variation with temperature (at P=25MPa) predicted by the PR EoS, the volume translated PR EoS and the Ideal Gas EoS along with the experimental data values from NIST for both water and n-decane. The improvement in liquid-like density predictions is evident from these plots. The VT correction was, however, not employed in the current study because in the temperature range under investigation (700K to 1000K) the densities are predicted accurately without it.
Simulating the thermal transport processes
also requires the determination of derived thermodynamic properties like the specific enthalpy of the mixture (hm ), constant pressure specific heat of the mixture (Cp,m ) and partial molar enthalpies of the component species (hi ). The functional form of the mixture specific internal energy (um ) is first determined using the departure function methodology and the rest of the properties can be derived from it as shown 44
(a)
(b)
Figure 2-3: Comparison of density predictions by the Volume-Translated PengRobinson EoS (VT-PR EoS), the simple Peng-Robinson EoS and the Ideal Gas EoS with NIST data at P = 25 MPa (a) water and (b) n-decane 45
in Eq.2.24 to Eq.2.30: um = um,IG + ∆u
(2.24)
where, um,IG is the specific internal energy of the mixture at the same temperature but P = 0 or V = ∞ (where the mixture would behave as an ideal gas mixture) and ∆u is the internal energy departure function which can be derived by integrating the fundamental relation of thermodynamics for the system from V = Vm to V = ∞ and using the Maxwell’s Relations: ZVm ∂P T ∆u = − P dV ∂T V
(2.25)
∞
The term
∂P ∂T V
in Eq. 2.25 above can be determined from the PR EoS and inte-
grating gives us the final form of the internal energy functional: √ # " V + 1 − 2 bm 1 dam √ um = um,IG + √ am − T ln dT V 2 2bm V + 1 + 2 bm
(2.26)
The specific enthalpy functional is easily obtained from Eq. 2.26 keeping in mind that h = u + P V : √ # " V + 1 − 2 bm 1 da m √ hm = hIG am − T ln (2.27) m + P V − Ru T + √ dT V 2 2bm V + 1 + 2 bm Differentiating the specific enthalpy with respect to T at a constant P yields the specific heat capacity at constant pressure (Cp,m ) as below: Cp,m =
" √ # V + 1 + 2 bm d2 am T IG √ ln = Cp,m − Ru + √ 2 2bm V + 1 − 2 bm dT 2 P 2 Ru dam 1 −Ru T 2am (V + bm ) −T − / + V − bm dT V ∗ V ∗2 (V − bm )2
∂hm ∂T
V ∗ = V 2 + 2bm V − b2m
(2.28)
(2.29)
Differentiating the specific enthalpy with respect to the mole fraction of specie i gives 46
the functional form for the partial molar enthalpy of specie i: hi =
∂hm ∂Xi
"
P,T
√ # " V m + 1 − 2 bm ∂am 1 √ ln + √ ∂Xi T,P 2 2bm Vm + 1 + 2 bm
#
Vi − Vm /bm bi V∗ P,Xi # ∂ 2 am 1 ∂am −T − am − T bi ∂Xi ∂T bm ∂T
= hIG i + P Vi − Ru T + am − T
"
dam dT
(2.30)
∂am ∂Xi
=2 T,P
X
√ xj (1 − kij ) ai aj
(2.31)
j
where, V¯i is the partial molar volume of species i given by Eq.2.32 below: Vi =
∂Vm ∂Xi
T,P
=
−1
∂P ∂Vm
Ru T 2am (Vm − bm ) Ru T + bi − 2 bi + (Vm − bm ) V ∗2 (Vm − bm )
∂am ∂Xi
V∗
T,Xi
(2.32)
∂P ∂Vm
= T,Xi
−Ru T 2am (Vm − bm ) 2 + V ∗2 (Vm − bm )
(2.33)
IG The ideal gas mixture properties at temperature T like Cp,m and hIG m are obtained by
simple mole fraction weighted average of the ideal gas properties of the constituent IG and hIG species (Cp,i i ). These in turn, are calculated using the temperature dependent
polynomials for ideal gas Cp and h variations from Yaws Handbook [57]. Fig. 2-4 shows the variation of Cp with temperature (at P = 25MPa) for water and n-decane predicted using the PR EoS and using the Yaws polynomials only along with the experimental data values from NIST. It is seen that the cubic PR EoS satisfactorily predicts the Cp values while the temperature dependent polynomials alone fail to do so.
47
T
(a) water
(b) n-decane
Figure 2-4: Comparison of constant pressure specific heat (Cp in kJ/kg − K) predictions by the Peng-Robinson EoS and the Ideal Gas EoS with NIST data at P = 25 MPa (a) water and (b) n-decane 48
substance water n-decane
ω 0.344 0.484
Π(Debye) κ 1.8 0.1 0.07 0.0
Table 2.4: Acentric factor, dipole moment and association factor for water and ndecane [32] [8]
2.2 2.2.1
Near-critical transport properties Viscosity and Thermal Conductivity
The calculation of viscosities and thermal conductivities was done using the correlations developed by Chung et al. [6]. They extended the Chapman-Enskog theory for dilute gas viscosities [4] by correcting for high pressure dense conditions using functions involving the acentric factor (ω), reduced dipole moment (Πr ) and association factor (κ) of the substance as parameters. The coefficients of these functions were determined by them by fitting experimental data for a vast number of pure substances. The ω, Π and κ values for water and n-decane are given in Table 2.4. The fluid viscosity (µ) is calculated as follows: µ = µk + µp
(2.34)
where the viscosities are in units of Poise µk = µ0
µp = 36.344 ×
1 + A6 Y G2
1/2 −6 (M Tc ) 10 A7 Y 2 G2 2/3 Vc
µ0 = 4.0785 ×
(2.35)
A9 A10 exp A8 + ∗ + ∗2 T T
1/2 −5 (M Tc ) 10 Fc 2/3 Vc Ω∗
(2.36)
(2.37)
where, µ0 is the dilute gas viscosity in Poise according to the Chapman-Enskog theory, M is the molecular weight of the substance in g/mol and Vc is the molar critical volume of the substance in cm3 /mol. Ω∗ is the reduced collision integral given by Eq.2.38 49
and Fc is a linear function in ω, Π4r and κ given by Eq.2.40. Ω∗ =
E 1.1615 C ∗B ∗W + + GT sin ST − H + exp (DT ∗ ) exp (F T ∗ ) T ∗B
(2.38)
where, A = 1.16145, B = 0.14874, C = 0.52487, D = 0.77320, E = 2.16178, F = 2.43787, G = −6.435 × 10−4 , H = 7.27371, S = 18.0323 and W = −0.76830. T ∗ = 1.2593
T Tc
(2.39)
Fc = 1 − 0.2756ω + 0.059035Π4r + κ
(2.40)
Y = ρVc /6
(2.41)
G1 =
1 − 0.5Y (1 − Y )3
1 1 − exp (−A4 Y ) + A2 G1 exp (A5 Y ) + A3 G1 G2 = A1 Y A1 A4 + A2 + A3
(2.42)
(2.43)
where, A1 to A10 are linear functions in ω, Π4r and κ given by: Ai = a0 (i) + a1 (i)ω + a2 (i)Π4r + a3 (i)κ
(2.44)
The polynomial coefficients a0 (i), a1 (i), a2 (i) and a3 (i) are given in Table 2.5. The reduced dipole moment (Πr ) can be calculated from the dipole moment of the substance using Eq. 2.45 below: Πr = 131.3
Π (Vc Tc )1/2
(2.45)
where, the dipole moment (Π) is in units of Debye, the critical temperature (Tc ) is in K and the molar critical volume is in cm3 /mol. A similar procedure is used to
50
i 1 2 3 4 5 6 7 8 9 10
a0 (i) 6.32402 0.12102 × 10−2 5.28346 6.62263 19.74540 −1.89992 24.27450 0.79716 −0.23816 0.68629 × 10−1
a1 (i) 50.41190 −0.11536 × 10−2 254.20900 38.09570 7.63034 −12.53670 3.44945 1.11764 0.67695 × 10−1 0.34793
a2 (i) −51.68010 −0.62571 × 10−2 −168.48100 −8.46414 −14.35440 4.98529 −11.29130 0.12348 × 10−1 −0.81630 0.59256
a3 (i) 1189.02000 0.37283 × 10−1 3898.27000 31.41780 31.52670 −18.15070 69.34660 −4.11661 4.02528 −0.72663
Table 2.5: Function coefficients used for generalized viscosity correlations from Chung et al.[6]
determine the thermal conductivity (λ) of near-critical fluids as follows: λ = λk + λp
(2.46)
where the thermal conductivities are in units of cal/(cmsK) λk = λ0
λp = 3.039 ×
1 + B6 Y H2
(Tc /M ) 10−4 2/3 Vc
1/2
B7 Y 2 H2 (T /Tc )1/2
λ0 = 7.452 (µ0 /M ) Ψ
Ψ = 1 + α1
0.215 + 0.2829α1 − 1.061α2 + 0.2667α3 0.6366 + α2 α3 + 1.061α1 α2
α1 =
(2.47)
Cv 3 − Ru 2
(2.48)
(2.49)
(2.50)
(2.51)
where, Cv is the ideal gas heat capacity at constant volume. α2 = 0.7862 − 0.7109ω + 1.3168ω 2 51
(2.52)
α3 = 2.0 + 10.5
T Tc
2
1 − exp (−B4 Y ) 1 H2 = B1 + B2 G1 exp (B5 Y ) + B3 G1 Y B1 B4 + B2 + B3
(2.53)
(2.54)
where, B1 to B7 are linear functions in ω, Π4r and κ given by: Bi = b0 (i) + b1 (i)ω + b2 (i)Π4r + b3 (i)κ
(2.55)
The polynomial coefficients b0 (i), b1 (i), b2 (i) and b3 (i) are given in Table 2.6. The viscosity and thermal conductivity of the water-n-decane mixture is calculated by a simple mass-fraction weighted average of the properties of the individual species at the particular temperature and pressure. This methodology to determine transport properties of a real-fluid mixture has also been employed before in the work of Narayanan et al. [31]. Chung et al. [6] in their work, do propose some complicated mixing rules to determine the effective properties of the mixture like Vcm , Tcm , ωm , Mm , Π4rm and κm which can then be used in Eq.2.34 to Eq.2.54 above to calculate the viscosity and thermal conductivity of the mixture. However, the application of these mixing rules yields unphysical negative values for the viscosity and thermal conductivity at certain values of the mixture mass-fractions. This causes artificial numerical instabilities in the simulation. Hence, the current study employs simple mass-fraction weighting for the calculation of mixture viscosity and thermal conductivity. Fig.2-5 shows the variation of viscosity with temperature (at P=25MPa) for water and ndecane. The predictions match closely with NIST data with small inaccuracies near room temperature conditions.
2.2.2
Mass Diffusivity
The mass diffusivity for the water-n-decane binary mixture is calculated using the Tracer Liu-Silva-Macedo (TLSM) model formulated by Liu et al. [23] [24] [46] as sug52
(a) water
(b) n-decane
Figure 2-5: Comparison of dynamic viscosity (µ in P a − s) predictions using Chung’s generalized correlations [6] with NIST data at P=25MPa (a) water and (b) n-decane 53
i b0 (i) b1 (i) b2 (i) b3 (i) 1 2.41657 0.74824 −0.91858 121.72100 2 −0.50924 −1.50936 −49.99120 69.98340 3 6.61069 5.62073 64.75990 27.03890 4 14.54250 −8.91387 −5.63794 74.34350 5 0.79274 0.82019 −0.69369 6.31734 6 −5.86340 12.80050 9.58926 −65.52920 7 81.17100 114.15800 −60.84100 466.77500 Table 2.6: Function coefficients used for generalized thermal conductivity correlations from Chung et al.[6]
gested by Kutney in his comprehensive review of near-critical fluid mass diffusivities [21]. The model is based on the hard-sphere diffusivity theory. The diffusivity of species 1 in species 2 (D12 ) is given by Eq. 2.56 below: T LSM D12
669.1M12 = ρNA σT2 LSM
r
Ru T 0.75ρNA σT3 LSM T12LSM exp − − 0.27862 M12 1.2588M12 − ρNA σT3 LSM kT (2.56)
σT2 LSM =
2 21/3 σ12 1/3 p 1 + 1.2 kT /T12LSM
(2.57)
T12LSM x1 x2 = 1 2 k k
(2.58)
σ12 = x1 σ1 + x2 σ2
(2.59)
M12 = x1 M1 + x2 M2
(2.60)
where, NA = 6.023×1023 /mol is the Avogadro constant and k = 1.3806×10−23 m2 kg/s2 is the Boltzmann constant. In Eq. 2.56, D12 is in units of cm2 /s, ρ is in units of gm/cm3 , M is in units of gm/mol, σ is in units of ˚ Aand /k is in units of K. The Lennard-Jones force constants for water and n-decane σ and /k were taken from Liu et al. (1998) [24] and are shown in Table 2.7. They may also be calculated from the 54
substance water n-decane
˚) σ(A 1.53091 6.71395
/k(K) 3788.51 434.86
Table 2.7: Lennard-Jones force constants for water and n-decane Liu et al. (1998) [24] critical properties of the substance using Eq. 2.61 and Eq. 2.62 below. 1/3
σi = 0.809Vci
(2.61)
where, σ is in ˚ A and Vc is in cm3 /mol i Tci = k 1.2593
2.3
(2.62)
Summary
The ability to capture the variations of thermodynamic and transport properties with temperature and pressure is crucial while simulating the heat and mass transfer processes of near-critical fluids. In this chapter, the details of the consistent treatment of near-critical thermodynamics of real fluid mixtures using the Peng-Robinson cubic equation of state were discussed. The use of the PPR78 approach with a group contribution method for inter-molecular interactions for the purpose of determining the thermodynamic parameters of the non-ideal mixture was also discussed. The calculation of viscosities and thermal conductivities are done using the generalized correlations of Chung et al. [6]. Mass diffusivity of one component in the other in the water-n-decane binary mixture is calculated based on the Tracer Liu-Silva-Macedo (TLSM) model developed by Liu et al. [23] [24] [46].
55
56
Chapter 3 Problem Formulation and Methodology 3.1
Governing equations
The conservation equations used to simulate the mixing of supercritical water and hydrocarbon (here, n-decane) are given below in Eq.3.1 to Eq.3.6, in the most general form in order to accommodate variable density, compressibility and variable transport property effects. Since the equations are solved numerically using a finite-volume (FV) formulation, they are expressed in the strong conservative form. The assumptions employed in each of the equations are also explained below. Continuity Equation: ∂ρ ∂ρuj + =0 ∂t ∂xj
(3.1)
where, ρ is the bulk mixture density and u = ui is the velocity vector. The Newton’s indicial notation has been used in all the equations. Momentum Equation (Navier-Stokes): ∂ρui ∂ρuj ui ∂P ∂τij + =− + + ρgi ∂t ∂xj ∂xi ∂xj
(3.2)
where, P is the pressure field, τ = τij is the stress tensor and g = gi is the gravitational 57
acceleration vector. Here, Stokes’ formulation of the stress tensor in terms of the symmetric part of the velocity gradient tensor (the strain tensor Sij ) is used to close the stress terms in the momentum conservation equation (Eq. 3.2) and is given by:
1 τij = 2µ Sij − Sii δij 3
=µ
∂ui ∂uj + ∂xj ∂xi
2 ∂ui − µδij 3 ∂xi
(3.3)
where, µ is the viscosity of the mixture (calculated using the methodology outlined in Sec.2.2.1) and δij is the Kronecker delta tensor. Stokes’ assumption of zero bulk viscosity is inherent in the above formulation of the viscous stress tensor. Species Transport Equation: ∂ρYk ∂ρuj Yk ∂ + = ∂t ∂xj ∂xj
∂Yk ρDk ∂xj
(3.4)
where, the subscript k denotes specie k, Yk is the mass fraction of specie k in the mixture and Dk is the mass diffusivity of specie k in the bulk mixture. In general, for a mixture with n species, one needs to solve (n−1) species conservation equations like Eq. 3.4 above and the mass fraction of the nth specie is obtained from the constraint P k Yk = 1. Since this study involves only binary mixtures of water and n-decane, only 1 species transport equation is solved. In Eq. 3.4, the species diffusive flux vector has been modeled using the Fick’s Law of Diffusion and Soret and Dufour diffusive effects are neglected. The mass diffusivity of specie k in the bulk mixture (Dk ) can be obtained from its binary diffusivities in each of the other components (Dkl ) using Eq. 3.5 below: 1 − Xk k6=l Xl /Dkl
Dk = P
(3.5)
where, Xk is the mole fraction of specie k in the mixture and the binary diffusivities (Dkl ) are determined using the TLSM model as explained in Sec.2.2.2. For a binary mixture, the diffusion coefficient of both species in the bulk mixture simply reduces to their binary diffusion co-efficient.
58
Enthalpy Transport Equation: X ∂ρh ∂ρuj h DP ∂ ∂T ∂ ∂Yk ∂ui + = + λ + hk ρDk + τij ∂t ∂xj Dt ∂xj ∂xj ∂xj ∂xj ∂xj |{z} | {z } |k | {z } {z } A B
C
(3.6)
D
where, h is the specific enthalpy of the mixture, T is the temperature, λ is the mixture thermal conductivity and hk is the partial enthalpy of specie k. Term B in the energy equation is the heat diffusion term modeled using the Fourier Law of Heat Diffusion. Term A represents the work done by the pressure forces on the fluid. The magnitude of this pressure work term with respect to the dominant enthalpy advection term in the energy equation can be estimated using the non-dimensional Eckert number defined below: Ec =
u2 ∆h
(3.7)
where, dh is a characteristic enthalpy difference, which could be taken to be the enthalpy difference between the water and hydrocarbon streams. At the low flow Mach numbers encountered in this study, the Eckert number ∼ 10−7 and so, the pressure work term may be safely neglected. Using a similar scaling argument, we can also neglect the viscous dissipation term in the energy equation (term D) since the Brinkmann number (see Eq. 3.8) Br ∼ 10−7 . µu2 Br = λ∆T
(3.8)
where, ∆T is a characteristic temperature difference in the problem (the temperature difference between the water and hydrocarbon streams). Term C is the transport of enthalpy due to species diffusion. Since the water and hydrocarbon have markedly different specific enthalpies at these conditions, this transport mechanism may be significant and hence must be considered to make accurate thermal transport predictions. Since the enthalpy of a real-fluid mixture depends on the temperature, pressure and species mass fractions, the correct formulation of the enthalpy differential must
59
include all of these contributions as shown in the equation below:
dh = Cp dT + v − T
∂v ∂T
dP + P
X
hk dYk
(3.9)
k
The enthalpy change due to pressure variation in Eq. 3.9 above can be neglected since the Eckert number is very small. Recasting the temperature gradient in the heat diffusion term in Eq.3.6 using Eq.3.9, the enthalpy transport equation can be reformulated including only the significant terms as: Modified Enthalpy Transport Equation: ∂ρh ∂ρuj h ∂ + = ∂t ∂xj ∂xj
λ ∂h Cp ∂xj
X ∂ λ ∂Yk + hk ρDk − ∂x Cp ∂xj j k
(3.10)
The above formulation of the enthalpy equation contains an implicit treatment of both the enthalpy advection and heat diffusion terms which has numerical advantages as explained in Section 3.4. The system of equations is closed using the equation of state for real-fluid mixtures to determine ρ and derived thermodynamic properties like Cp , h as functions of P, T and the mixture composition. Since, the enthalpy field is calculated by solving the enthalpy transport equation, an iterative procedure to determine the T field from the enthalpy, pressure and composition fields has to be employed because of the non-linearity of the enthalpy constitutive relationship for near-critical fluid mixtures. Typically, 3-4 iterations are sufficient to determine the temperature with acceptable accuracy (error ∼ 10−6 ). The viscosity and thermal conductivity are calculated using the procedure explained in Sec.2.2.1.
3.2
Geometry and Boundary Conditions
The mixer geometry is a cylindrical tee as shown in Fig. 3-1. The geometry corresponds to a small-scale test reactor with cylindrical tube inner diameter (d) of 2.4 mm. The two inlet branches of the tee meet at an angle of 90o with a long downstream section of length 22D. The inlet sections are each of length 3D from the inlet 60
Figure 3-1: Cylindrical tee mixer geometry and boundary conditions
to the center of the mixing joint. Hot supercritical water enters from the horizontal inlet while the colder hydrocarbon (n-decane) stream enters from the vertical inlet. Since the Reynolds number at both inlets is in the laminar regime for all of the simulated cases, a parabolic inlet profile velocity boundary conditions is imposed at both the water and hydrocarbon inlet with specified average bulk velocities for each. No slip boundary condition is applied at all walls. At the outlet, a zero gradient BC with flux correction is applied (to ensure that the total mass flux entering the reactor is same as the total mass flux exiting). A fixed value temperature BC is applied at both inlets with values corresponding to the specific inlet conditions of the simulation. All the walls are adiabatic and the outlet is far enough downstream to assume that the temperature gradient at the outlet is zero. A zero gradient (normal to the face) boundary condition is imposed for the pressure at both inlets as also on all walls. This boundary condition is not an intuitive one but can be easily derived from the momentum equation (in the boundary normal direction) considering that all other terms are zero. At the outlet, a fixed back-pressure of 25 MPa (the nominal reactor operating pressure) is specified. Since the flow Mach number is low, one does not have to worry about pressure waves being reflected back into the domain from the 61
Figure 3-2: Cylindrical tee reactor mesh outlet face.
3.3
Meshing
The computational domain is meshed using the meshing software GAMBIT by ANSYS Inc. (see Fig. 3-2). The entire domain has to be smartly decomposed into many blocks and each block meshed separately so as to obtain the best quality mesh possible (minimize non-orthogonality of mesh elements). All mesh elements are hexahedra. The adequacy of the spatial resolution is confirmed a posteriori by grid independence checks as shown in Sec.4.1 and Sec.5.1.
3.4
Numerical methodology
The conservation equations in Section 3.1 are discretized using the finite-volume methodology (FVM). The computational code was developed using the set of opensource CFD libraries, OpenFOAM by OpenCFD Ltd.[35]. The underlying use of the mid-point rule to approximate the integrals of fluxes over the control volume surfaces 62
and the integrals of quantities over cell volumes limits the order of accuracy of the numerical method to 2nd order at best. The details of the numerical methodology employed are presented in the following sub-sections 3.4.1 to 3.4.5.
3.4.1
Discretization Schemes
A 2nd order accurate central difference scheme with a correction for mesh nonorthogonality is used for the gradient of velocity, temperature, pressure and mass fractions at the face centers. A simple bi-linear interpolation is used to determine the mass flux per unit area at the face centers. The velocity, temperature and mass fractions are interpolated to the cell face centers using a 2nd order accurate linear interpolation scheme with a limiter, which adds upwind scheme based numerical diffusion in regions of steep gradients to smooth out unphysical oscillations. The transport properties and pressure are interpolated to the cell face centers using a 2nd order accurate linear interpolation scheme. Time marching is done using a 2nd order accurate BDF2 (Backward difference formula) semi-implicit scheme.
3.4.2
Discretized Equations
The discretized forms of the continuity, momentum and enthalpy equations are shown in Eq. 3.11, Eq. 3.12 and Eq. 3.13 respectively. Discretized form of the Continuity Equation: h in+1 X V n+1 n n−1 (3ρ − 4ρ + ρ ) + Af ρf (ui,f )nf =0 3δt f
(3.11)
where, the superscripts n+1, n, n-1 denote the new time-step level, old time-step level and old old time-step level respectively. Discretized form of the Momentum Equations: [Ani ] un+1 = [Hin ] un+1 − ∆i P n+1 + [Sin ] i i 63
(3.12)
n+1 ui is the ui solution vector of length N (where N is the total number of grid cells). Three sets of N equations like the one shown above in Eq. 3.12 need to be solved; one for each velocity component. [Ani ] is the diagonal part of the co-efficient matrix of the ui equation and includes contributions of all terms involving un+1 i,m ; the subscript m referring to the cell number m. [Hin ] is the off-diagonal part of the co-efficient matrix including contributions from all terms involving un+1 i,nb ; where the subscript nb refers to all neighbor cells of cell m. The co-efficients in both these NXN matrices are calculated using variable values from the previous time-step level; hence the superscript n. ∆i is the discrete gradient operator. [Sin ] is the equation source vector (RHS vector) and includes contributions of all terms involving uni,m and uni,nb (the explicitly treated terms) as well as the source terms of the momentum equation (like gravity force term). The advective transport terms in the momentum conservation equations are non-linear and are linearized while writing the discrete form of the equations by treating the mass flux terms explicitly using their values from the nth time-step level. Therefore, the discretized equations above are not fully implicit in ui and hence, the time-stepping scheme has been referred to as semiimplicit. Discretized form of the Enthalpy Equation: [B n ] hn+1 = [G n ] hn+1 − [Qn ]
(3.13)
where, the matrices [B n ] and [G n ] are analogous to [Ani ] and [Hin ] in the discretized momentum equations. [Qn ] is the enthalpy equation source vector and includes the heat transport due to species diffusion and the part of the temporal derivative of enthalpy involving previous time-step/iteration values. All the advective and diffusive terms are treated implicitly in the h equation and their contributions appear in [Bin ] and [Gin ]. The discrete species conservation equation also has a similar form as Eq. 3.13:
64
Discretized form of the Species Conservation Equation: [Ckn ] Ykn+1 = [Fkn ] Ykn+1 − [Rnk ]
(3.14)
It is clear from Eq. 3.12 that the pressure and velocity fields at the new time-step level are coupled. In order to decouple the equations, an operator-splitting approach is adopted based on the Pressure Implicit Splitting of Operators (PISO) algorithm proposed by Issa [13]. In this approach, an intermediate velocity field u∗i is predicted by solving the set of discretized momentum equations by omitting the contribution of the pressure term. This predicted velocity field in general, will not satisfy mass conservation. Therefore, the new time-step pressure field pn+1 is calculated such that use of this field to correct the velocities will cause both mass and momentum conservation equations to be satisfied. A discrete equation to determine the correct pressure field can be formulated by interpolating the velocities in Eq. 3.12 to the face centers and substituting these in the discrete continuity equation (Eq. 3.11) giving Eq. 3.15 below: Discrete Pressure Equation: " # h i X n+1 X V ∂P −1 Af ρn+1 Ani,f Af ρn+1 u∗i,f nf + (3ρn+1 −4ρn +ρn−1 )+ =0 f f 3δt ∂xi nf f f (3.15) n −1 where, Ai,f is the inverse of the diagonal matrix [Ani ] in Eq. 3.12 interpolated to the face center of face f. This technique of interpolating the pressure correction term to the face centers is called the Rhie-Chow interpolation proposed by Rhie et al. [41] and prevents odd-even decoupling of the pressure and velocity fields on collocated grids.
3.4.3
Solution Algorithm
The solution algorithm is constructed using the above discrete momentum, species, temperature and pressure equations and is explained below in a step-wise fashion. Step 1 (Velocity Predictor): 65
Predict the velocity-field by solving the discrete momentum equation by omitting the pressure contribution and using mass flux values (ρf (ui,f )nnf ) from the previous time-step (n): [u∗i ] = ([Ani ] + [Hin ])−1 ([Sin ])
(3.16)
Step 2 (Species Mass Fractions Predictor): Predict the species mass fraction fields Yk∗ by solving the discrete species conservation equations using mass flux values (ρf (ui,f )nnf ) from the previous time-step (n): [Yk∗ ] = ([Ckn ] + [Fkn ])−1 ([Rnk ])
(3.17)
Step 3 (Enthalpy Predictor): Predict the enthalpy field h∗ by solving the discrete enthalpy equation using mass flux values (ρf (ui,f )nnf ) from the previous time-step (n) and predicted species mass fractions (Yk∗ ): [h∗ ] = ([B n ] + [G n ])−1 ([Qn ])
(3.18)
Step 4 (Temperature, Density and Thermodynamic Properties Correction): Update the temperature field to T ∗ using the enthalpy field calculated in Step 3 using an iterative procedure. Correct the density field and derived thermodynamic properties (like Cp ) using the EoS and P n , T ∗ and Yk∗ . Step 5 (Pressure and Velocity Correction): Solve the discrete pressure equation to obtain the corrected estimate of the pressure field (P ∗ ): " ∗ n+1 # h i X X ∂P V −1 =0 (3ρ∗ − 4ρn + ρn−1 ) + Af ρ∗f u∗i,f nf + Af ρ∗f Ani,f 3δt ∂xi nf f f (3.19) Here, the estimate of the new density field ρ∗ is used instead of the actual new density field ρn+1 . This is the variable density, low Mach number formulation which yields sufficient accuracy and numerical stability for the low Mach number flows simulated in this work. For higher Mach number (highly compressible) flows, this term can be
66
easily corrected by accounting for the density change with pressure as below: ρn+1 = ρ∗ + ψ ∗ (P ∗ − P n ) ∂ρ ψ= ∂P T,Yk
(3.20)
This correction to the variable density term though small at low Mach numbers, improves the numerical stability of the pressure equation and hence, has been employed. In Eq. 3.19, the pressure gradient at the cell face centers has to be corrected for the non-orthogonality of the cells. This is done using an explicit correction term in the numerical scheme for the gradient (deferred correction approach). As such, Eq. 3.19 needs to be solved iteratively (typically only 2-3 times) in order to get an accurate estimate of the pressure field. The velocity field is corrected using P ∗ to give u∗∗ i . Since ui at cell m depends on the ui values at neighboring cells, the corrected velocity field now does not satisfy the momentum equation. Hence, a second corrector step for pressure and velocity must be executed to obtain p∗∗ and u∗∗∗ . Issa [13] showed in his work that 2 corrector steps are sufficient to ensure 2nd order accuracy for the velocity field. Since, the FVM formulation itself is only 2nd order accurate at best, further corrector steps will yield no benefit. Step 6 (Species Mass Fractions Corrector): Correct the species mass fraction fields to obtain the new time-step value Ykn+1 by solving the discrete species conservation equations corrected mass flux values (ρ∗f u∗∗∗ i,f nf ): n+1 (3.21) Yk = ([Ck∗∗∗ ] + [Fk∗∗∗ ])−1 ([R∗∗∗ k ]) Step 7 (Enthalpy Corrector): Correct the enthalpy field to obtain the new time-step value hn+1 by solving the discrete enthalpy equation using corrected mass flux values (ρ∗f u∗∗∗ i,f nf ) and new time-step species mass fractions (Ykn+1 ):
hn+1 = ([B ∗∗∗ ] + [G ∗∗∗ ])−1 ([Q∗∗∗ ]) 67
(3.22)
Step 8 (Second Temperature, Density and Thermodynamic Properties Correction): Obtain the new time-step values of the temperature field, density field and derived thermodynamic properties using the EoS and P ∗∗ , hn+1 and Ykn+1 . Step 9: Repeat Step 5 to obtain the new time-step values of the pressure and velocity fields (P n+1 and un+1 ). i Step 10 (Transport Properties Correction): Correct transport properties using Ykn+1 , Tkn+1 and Pkn+1 .
3.4.4
Numerical Stability
The implicit formulation of the species, enthalpy and pressure equations ensures unconditional numerical stability for these equations. If the diffusion term in the enthalpy equation were treated explicitly, a highly restrictive stability condition on the time-step size would have had to be followed: δt ≤
(δx)2 2α
(3.23)
where α is the thermal diffusivity. This means that if the mesh size is refined by a factor of 2, the time-step has to be decreased by a factor of 4. This is the reason why the enthalpy equation was formulated in an implicit fashion as in Eq.3.10 . As explained previously, the momentum equations are actually semi-implicit with the mass-flux terms taken from the previous time-step and hence the time-step must satisfy the Courant-Lewy-Friedrichs (CFL) stability condition: δt ≤
1 ux δx
+
uy δy
+
uz δz
(3.24)
The time-step size is adjusted dynamically at each time-step to ensure that the above condition is satisfied at every cell in the domain. 68
3.4.5
Linear System Solvers
Solution of each of the discretized equations involves solving a large linear system of equations in the respective variable. Appropriate choice of the linear system solver for each equation set is crucial for computational speed. For the momentum, species and enthalpy equations, a Preconditioned Bi-Conjugate Gradient (PBiCG) method with a diagonal-incomplete LU preconditioner for the matrix is used. For the pressure equation, a Conjugate Gradient (CG) method with a diagonal-incomplete Cholesky preconditioner is used. Since the size of the equation set is very large (N ≈ 1million), the equation set is distributed to 48 processor cores (on the Pharos computing cluster) using the domain decomposition method and solved in parallel with OpenMPI interprocessor communication.
3.5
Summary
In this chapter, the fundamental conservation equations governing the mixing of supercritical water and hydrocarbons are presented and the various simplifying assumptions made are discussed (Sec.3.1). In Sec.3.2 the cylindrical tee-reactor geometry is described and the boundary conditions imposed are discussed. Sec.3.3 focuses on computational mesh generation. This is followed by a detailed discussion of the finitevolume methodology, discretization schemes, discrete forms of conservation equations, the solution algorithm, linear system solvers and numerical stability in Sec.3.4.
69
70
Chapter 4 Mixing at intermediate Re: Flow dynamics and impact of temperature difference In this chapter, the mixing of supercritical water and n-decane was simulated at a Reynolds number of 500 at the water inlet for two different inlet temperature conditions as shown in Table 4.1. In both cases, the mass flow rates through the water and n-decane inlets is the same (a design condition found to be favorable to the process from the chemical kinetics point of view). The nominal pressure in the tee mixer was ≈ 25MPa (ensured by fixing the back-pressure at the outlet) and the n-decane inlet temperature was 700 K, which is sufficiently above the Upper Critical Solution Temperature (UCST) of the water n-decane system (≈ 632 K) thereby ensuring complete miscibility of the two components. The water inlet temperature in Case I was 800 K and that in Case II was 1000 K. This will help us examine the impact of the different degrees of variation of properties like density, viscosity with temperature on the flow and mixing dynamics in a cylindrical tee mixer. In order to understand the differential impact that the near-critical conditions have on the flow field and the dynamics of mixing, we also perform a set of simulations (Case III and Case IV) with the same inlet conditions as in Case I and Case II, but with the physical properties of the pure fluid components held constant at their inlet values. The thermal trans71
Case
Rew,in
Red,in
I II III IV
500 500 500 500
204 239 204 239
m ˙ w,in (kg/hr) 0.11 0.13 0.11 0.13
m ˙ d,in (kg/hr) 0.11 0.13 0.11 0.13
uw,in (m/s) 0.0832 0.1416 0.0832 0.1416
ud,in (m/s) 0.0144 0.0168 0.0144 0.0168
Tw,in (K) 800 1000 -
Td,in (K) 700 700 -
Property variations yes yes no no
Table 4.1: Inlet flow conditions for the cases discussed in this chapter port is not solved for in these cases. The effect of gravity forces is not included in these simulations so as to allow us to study the impact of geometry, Reynolds number and property variation effects on the flow field and mixing dynamics. The impact of gravity forces, which could play a significant role considering that the Froude number, Fr =
u2 gD
∼ 1, will be the focus of a future study.
The physical properties of the mixture like the density and viscosity can be expected to have a significant impact on the dynamics of the flow. Therefore, it is important to note at this stage, the variation of the mixture density and viscosity with both temperature and composition. Fig.5-17 and Fig.5-18 show the variation of the mixture density (ρ in kg/m3 ) and dynamic viscosity (µ in P a − s) with temperature for the range of temperatures pertinent to this study, for different mixture compositions from pure water to pure n-decane. Table 4.2 shows the ρ and µ for the water and n-decane pure components at the extreme temperatures within the tee mixer for the cases I and II. It can be seen that the mixture density and viscosity are strong functions of the temperature under these conditions. The density of the water component increases by around 100% as it gets cooled from 1000K to 700K within the tee mixer domain in Case II while the density of the n-decane component decreases by around 30% as it gets heated from 700K to 1000K. The viscosity of water does not change appreciably in this temperature range. However, the viscosity of the n-decane component decreases by more than 100% as it gets heated from 700K to 1000K. These variations in ρ and µ influence the local Reynolds number which in turn can be expected to have a significant bearing on the local flow dynamics and stability. These variations in fluid physical properties within the flow domain will be much smaller in Case I with only a 100K temperature range. The comparison of these 72
Figure 4-1: ρ (in kg/m3 ) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) two cases (Case I and Case II) will help us evaluate the impact of varying degrees of property changes on the flow dynamics and mixing in the tee mixer.
4.1
Validation: Grid convergence tests
In order to verify the adequacy of the chosen spatial resolution and validate the numerics of the code, a grid independence study was performed for both Case I and Case II. Fig.4-3 shows the streamwise velocity profile along the vertical direction in Case I at three different axial locations (x=4D, x=6D and x=8D) for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm). A mesh resolution of 0.10 mm seems to be sufficient to obtain grid-independence. Fig.4-4 shows the mean streamwise velocity profile along the vertical direction in Case II at three different axial locations (x=4D, x=6D and x=8D) for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm). The averages are taken over 4 tf low−through time periods after the initial transient of 2 tf low−through time periods. Grid-independence is attained for a mesh 73
Figure 4-2: µ (in Pa-s) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane)
Case I I I I II II II II
Specie water water n-decane n-decane water water n-decane n-decane
ρ(kg/m3 ) 125.62 84.22 486.27 424.45 125.62 57.86 486.27 329.94
T (K) 700 800 700 800 700 1000 700 1000
µ(P a − s) 3.1773e − 05 3.3626e − 05 8.2378e − 05 5.8599e − 05 3.1773e − 05 3.9332e − 05 8.2378e − 05 4.1035e − 05
Table 4.2: Properties of water and n-decane at 25MPa and extreme temperatures of the simulations I and II
74
resolution of 0.06 mm in this case. The velocity field thus seems to be well resolved. Fig.4-5 and Fig.4-6 show the mean temperature and n-decane mass fraction profile respectively along the vertical direction in Case II at three different axial locations (x=4D, x=6D and x=8D) for the three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm). These figures indicate that the mixing layer is well resolved in terms of the scalar fields. The mixing layer thickness from these plots can be seen to be around 0.8 mm around x = 6D. We will later see that in Case II, we observe that the shear layer becomes unstable. The wavelengths of the unstable perturbations or waves on the shear layer are multiples of the mixing layer thickness as stated in literature, for example [30]. Since the grid size of 0.06 mm gives at least 10 cells within the mixing layer, we can conclude that the spatial resolution is sufficient to capture the scales of instability in the shear layer. This gives us reasonable confidence in the consistency of the numerical methods and code. Unfortunately, there is no experimental data available in the literature for mixing of near-critical fluids in a cylindrical tee reactor geometry.
4.2
Case I: Water n-decane mixing with small temperature difference (Rew,in = 500, Tw,in = 800K, Td,in = 700K)
Under these flow conditions, the flow remains laminar and steady state is reached. Fig.4-7 shows the contour plots of the n-decane mass fraction and temperature fields on the centerplane (y=0 plane) at steady-state. The variations in these scalars can be seen more clearly (especially in the thin mixing layer), in the plots of Yd and T along the vertical centerline at different downstream locations (x=2D, x=4D, x=6D, x=8D) in Fig.4-8. It is clear from these plots, that the mixing layer is slowly advected downwards as it travels downstream. Based on rough mixing layer thickness estimates from these plots, we can see that the mixing layer thickness increases from around 0.4 mm at x = 2D to around 0.6 mm at x = 8D. The growth of the mixing layer under 75
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 4-3: u (in m/s) v/s z (in m) for simulated Case I for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D 76
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 4-4: umean (in m/s) v/s z (in m) for simulated Case II for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D 77
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 4-5: Tmean (in K) v/s z (in m) for simulated Case II for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D 78
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 4-6: Yd,mean v/s z (in m) for simulated Case II for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D 79
these laminar conditions is only due to the diffusion of the scalars in the transverse direction across the mixing layer. While the Yd and T remains fairly unchanged in the bulk of the water stream except near the thin mixing layer, the T increases and the Yd decreases by a significant amount in the n-decane stream close to the top of the pipe. This is due to the advective mixing brought about by the action of a counter-rotating vortex pair (CVP) in the body of the n-decane jet and a secondary flow of water around the HC jet as explained in detail in Sec.4.6. Fig.4-9(a) shows the profile of the velocity magnitude along the vertical centerline at different downstream locations. It shows a 3 times increase in the velocity of the n-decane stream from x=2D to x=8D. The peak velocity in the water stream remains fairly constant between these locations however, the thickness of the water stream as seen on the centerplane is continually reducing. These geometry dependent flow dynamic effects are explained in detail in Sec.4.6. The changes in local composition and temperature are reflected in changes in the local density and viscosity as shown in Fig.4-9(b)&(c) respectively. The reduction in density of the fluid mixture in the n-decane rich stream at the top of the pipe due to mixing with water is seen to be stronger than the reduction in viscosity due to the fact that the local density is a function of the local mole fractions of water and n-decane. Since water has a much smaller molecular weight than n-decane, a given mass fraction of water will contribute more moles to the mixture than the same mass fraction of n-decane. As a result, the mixture density is always strongly biased towards the water density. The net effect of the flow acceleration of the n-decane rich stream and the density, viscosity variations on the local flow Reynolds number is shown in Fig.4-9(d). Here, the pipe diameter is used as the length scale for the Re definition. The Re in the n-decane rich stream continually increases as we go downstream, while that in the water rich stream does not change much. This seems to be mainly due to the flow acceleration rather than the variations in physical properties. The local Re does not go above 1360 and it seems that the flow is stable to perturbations at these values of the local Re. Fig.4-10 shows contour plots of Yd (top), T (middle) and the streamwise component of the vorticity vector, ωx (bottom) on cross-sections of the cylindrical tee mixer 80
Figure 4-7: Contour plots on the centerplane (y=0 plane) at steady-state for simulated Case I of (top) n-decane mass fraction (Yd ) field and (bottom) Temperature (T) field; The white vertical lines indicate the positions x=2D, x=4D, x=8D, x=10D, x=12D, x=14D downstream of the center of the mixing joint from left to right (planes of constant X) at different locations downstream of the mixing joint (x=2D, x=4D, x=8D, x=16D, x=22D). It can be seen that mixing predominantly occurs due to the large scale transport of species caused by the action of a counter-rotating vortex pair (CVP) in the hydrocarbon jet coming in from the top (labeled (A) in Fig.4-10 (bottom)) and the flow of water around the jet towards the top (labeled (B) in Fig.4-10 (bottom)). Far downstream, as the strength of the CVP diminishes due to diffusion and the upward flow of water around the HC jet becomes weak, mixing is dominated by diffusion of species over small length scales and hence, is extremely slow. The main role of the supercritical water stream in this process being to heat up the HC stream in a controlled fashion (along with providing a flux of water into the HC stream), it is important to understand the physical mechanisms contributing to the local heating/cooling of the fluid in the tee mixer. This physical insight is difficult to obtain by looking at the enthalpy transport equation, wince the local enthalpy of the fluid mixture is a complex function of both the local temperature and the local 81
(a) Yd
(b) T (in K)
Figure 4-8: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case I of (a) n-decane mass fraction (Yd ) field and (b) temperature (T in K) field at different downstream locations (x=2D, x=4D, x=6D, x=8D)
82
(a) Umag (in m/s)
(b) ρ (in kg/m3 )
(c) µ (in Pa-s)
(d) Re = ρUmag D/µ
Figure 4-9: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case I of (a) velocity magnitude (Umag in m/s) (b) density (ρ in kg/m3 ) (c) dynamic viscosity (µ in Pa-s) and (d) Reynolds number (Re = ρUmag D/µ) at different downstream locations (x=2D, x=4D, x=6D, x=8D)
83
2D
4D
8D
16D
22D
2D
4D
8D
16D
22D
2D
4D
8D
16D
22D
Figure 4-10: (top) Mass fraction of n-decane (Yd ) contours ( middle) Temperature (T in K) contours and (bottom) Streamwise vorticity (ωx in s−1 ) contours for simulated Case I (Rew,in = 500, Tw,in = 800K, Td,in = 700K) at different downstream crosssections (x=2D, x=4D, x=8D, x=16D, x=22D)
84
composition. As such, an increase in enthalpy of a fluid elements as it moves in the tee mixer cannot be attributed directly to a corresponding increase in its temperature. It is therefore necessary to analyze the terms in a transport equation formulated in terms of the temperature to gain insight into the local heating/cooling of the fluid. A temperature equation can be obtained by substituting for dh from Eq.3.9 in Eq.3.6: X X DT ∂ ∂T ∂ ∂Yk ∂Yk ∂Yk ρCp = λ + ρDk hk − + ρuj hk ρ Dt ∂xj ∂xj ∂x ∂x ∂t ∂xj j j k k Substituting
∂ ∂xj
(4.1)
h i h i ∂Yk ∂Yk k + ρu ρDk ∂Y for ρ j ∂xj from the species conservation equation ∂xj ∂t
Eq.3.4 gives: 1 ∂ 1 X DT ∂T ∂Yk ∂hk = + λ ρDk Dt ρCp ∂xj ∂xj ρCp k ∂xj ∂xj | {z } | {z } A
(4.2)
B
In Eq.4.2 above, term A is the rate of increase/decrease of T (heating/cooling) of a moving fluid element due to heat diffusion (energy transfer through molecular collisions). Term B represents the rate of heating/cooling of the fluid element due to the diffusion of species along a gradient of partial enthalpy. As molecules of a particular species (say, specie k) diffuse along a gradient in partial enthalpy, the difference in partial enthalpy is taken from/released to the fluid element at that location leading to the cooling/heating of the fluid element represented by term B in the T equation, Eq.4.2. This gradient in the partial enthalpy, in general, can arise due to a gradient in mixture temperature, pressure or composition. The effect of variation of species partial enthalpies with pressure can be safely neglected due to the very low Mach numbers (≈ 0.0001) and consequently, low P differences (≈ 1-2 Pa) in this study. This heating/cooling mechanism produces an interesting effect for mixing of nonideal mixtures such as the water n-decane system under isothermal conditions. When ideal-gases mix under initially isothermal conditions, the flow remains isothermal (in the absence of chemical reactions) since term (A) (heat diffusion) vanishes and term (B) in Eq.4.2 also vanishes since the partial enthalpies of components in ideal mixtures are only functions of temperature and as such, there are no partial enthalpy
85
gradients in the flow domain. However, for mixing of non-ideal fluids as is the case with the water n-decane system in this study, the partial enthalpies of the species (hk ) are functions of not only temperature, but also mixture composition. As a result, the heating/cooling term (B) of Eq.4.2 is non-zero even for initially isothermal flows. Consequently, the flow becomes non-isothermal. This local cooling/heating of an initially isothermal flow is well illustrated in Fig.4-11(top) which shows the T contours on the tee centerplane at steady-state for a simulation with the same flow conditions as Case I, but with both streams entering at the same Tinlet of 750K. In the mixing joint region, the T drops to as low as 737 K in the mixing layer, where the diffusion of species is strongest. The thickness of the fluid layer affected by the cooling continues to grow downstream where it is also assisted now by heat diffusion which starts to grow. Fig.4-11(bottom) shows P k ∂hk contours of term (B) ( ρC1 p k ρDk ∂Y ) on the tee centerplane. It is clear that this ∂xj ∂xj term is negative in the mixing layer and is responsible for the T reduction of the fluid below 750K. Because the near-critical mixture of water and n-decane is a non-ideal mixture, the partial enthalpy (hk ) of both species in the mixture is dependent on not only the mixture temperature, but also the mixture composition. This fact is clearly illustrated in Fig.4-12 which shows curves for the variation of the partial enthalpy of both water and n-decane with T for different mixture compositions. We can see that the partial enthalpy of water in the mixture increases with higher n-decane mass fraction of the mixture. Similarly, the partial enthalpy of n-decane in the mixture increases with higher water mass fraction of the mixture. Water molecules have a greater affinity for other molecules than for n-decane molecules. Consequently, ndecane molecules taking the place of water molecules in a mixture (increase of Yd ) leads to more repulsive interactions and hence a greater partial enthalpy for the water component. This is the reason why we observe the trend of increase in hwater with increasing Yd . The same explanation holds for the trend in hn−decane . The species (water and n-decane) diffuse from regions in which their concentration is higher to regions in which their concentration is lower (in accordance with the 2nd Law of Thermodynamics). Thus, they are diffusing from towards regions in which 86
their partial enthalpy in the mixture is higher. That is, the gradient in Yk and hk are in opposite directions. Consequently, the term B of Eq.4.2 is negative (or zero) everywhere in the domain as observed.
This local cooling of fluid due to species diffusion along a positive partial enthalpy gradient also affects the heating of the HC stream (when Tw,in > Td,in ) as in Case I. Fig.4-13 shows term A and Term B of Eq.4.2 on cross-sections of the tee at different downstream locations for Case I. We can see that the local rate of cooling due to species diffusion is finite. However, the heat diffusion term is strongly dominant. In Fig.4-14 we can see that the rate of cooling due to heat diffusion is more than 3 times that of the rate of cooling due to species diffusion at most points in the mixing layer along the vertical centerline at x=4D. Thus, this cooling effect only slightly hinders the process of heating the n-decane stream which is a working objective of the mixing process in the tee mixer. Also, we can see from Fig.4-13 and Fig.4-14 that the term B is very small (almost zero) in the water-rich part of the mixing layer even though significant species diffusion occurs in this region. This is due to the fact that the partial enthalpy curves for both water and n-decane are very close together in the low Yd range (The curves corresponding to Yd = 0 and Yd = 0.2 in Fig.4-12 almost overlap). As a result, the heating/cooling of fluid due to species diffusion is small in the water-rich region.
The transport of heat between the two streams is seen to be quicker than the transport of species due to molecular heat diffusion being 2-3 times faster than molecular mass diffusion in regions of high n-decane concentration. This is well illustrated in Fig.4-15 which shows contours of the Lewis number, the ratio of the thermal diffusivity to the mass diffusivity (Le = α/D) over the cross-section of the tee at different downstream locations. The Lewis number ≈ 2 − 3 in regions rich in n-decane and ≈ 1 in water-rich regions. This difference in the rates of thermal and species mixing is more noticeable at far downstream locations where transport is mainly diffusion controlled. 87
T (in K)
1 ∂ ρCp ∂xj
1 ρCp
P
k
∂T λ ∂x j
(in K/s)
k ∂hk ρDk ∂Y ∂xj ∂xj (in K/s)
Figure 4-11: Contour plots on the centerplane (y=0 plane) at steady-state for the initially isothermal simulation: (top) Temperature (T) field (in K) (middle) Rate of 1 ∂ ∂T heating/cooling due to heat diffusion ρCp ∂xj λ ∂xj (in K/s) (bottom) Rate of heatP k ∂hk ing/cooling due to species diffusion along partial enthalpy gradient ρC1 p k ρDk ∂Y ∂xj ∂xj (in K/s); The white vertical lines indicate the positions x=2D, x=4D, x=6D downstream of the center of the mixing joint from left to right
88
Figure 4-12: partial enthalpy of water, hw and n-decane, hd (in kJ/kg) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane)
89
(a)
(b)
1 ∂ ρCp ∂xj
1 ρCp
P
k
∂T λ ∂x (in K/s) j
k ∂hk ρDk ∂Y ∂xj ∂xj (in K/s)
Figure 4-13: Contours of the rate of fluid heating/cooling (in K/s) due to (a) heat 1 ∂ ∂T diffusion ( ρCp ∂xj λ ∂xj ) and (b) species diffusion along partial enthalpy gradient P k ∂hk ( ρC1 p k ρDk ∂Y ) at different downstream locations (x=2D, x=4D, x=8D, x=16D, ∂xj ∂xj x=22D from left to right)
90
Figure 4-14: rate of fluid heating/cooling (in K/s) due to (a) heat diffu∂T sion ( ρC1 p ∂x∂ j λ ∂x ) and (b) species diffusion along partial enthalpy gradient j P k ∂hk ( ρC1 p k ρDk ∂Y ) along the vertical centerline at x = 4D for simulated Case I ∂xj ∂xj
2D
4D
8D
16D
22D
Figure 4-15: Lewis number, Le = α/D at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D)
91
4.3
Case II: Water n-decane mixing with large temperature difference (Rew,in = 500, Tw,in = 1000K, Td,in = 700K)
In Sec.4.2, we saw that the shear layer between the water and n-decane streams is stable and the flow remains laminar for a temperature difference of 100K between the two streams. When, the water enters at a temperature of 1000K (300K greater than the n-decane inlet), the water n-decane shear layer becomes unstable and the flow in the cylindrical tee mixer is unsteady. Fig.4-17 shows the n-decane mass fraction contours on the centerplane of the tee mixer at different time instants (t=2.0s, t=2.4s, t=2.8s, t=3.2s roughly corresponding to 2.5, 3, 3.5 and 4 tf low−through respectively) and the mean field. The flow was allowed to develop for 2 tf low−through to wipe out the initial conditions in the tee mixer and thereafter, averaging of field variables was performed over the next 2 tf low−through . Statistically stationary state is achieved after about 2 to 2.5 tf low−through . Since the statistics of the field variables do not change with time after about 2 to 2.5 tf low−through , we suspect some periodicity in the temporal variations of the fields. Fig.4-16 shows the Fourier transform of the n-decane mass fraction field at x = 6D, y = 0, z = 0 just after the shear layer destabilizes. We can see a cluster of frequency modes around 50 Hz which dominate the flow field fluctuations. The plots of Fig.4-17 show the location of onset of instability to be close to 5 diameters downstream of the mixing joint center with the shear layer developing some waviness at this location. Further downstream, the shear layer starts to roll up leading to faster mixing through large scale advective transport of species across the two streams. This effect is manifested in the thickening of the mean mixing layer beyond x=5D as seen in the mean Yd plot in Fig.4-17(e). The sudden acceleration in the growth of the mixing layer is more clearly visible in the Yd and T plots at different downstream locations in Fig.4-18. Fig.4-19(a) shows the profile of the mean velocity magnitude along the vertical centerline at different downstream locations. Similar to Case I, the acceleration of the n-decane rich stream 92
Figure 4-16: Fourier transform of the Yd temporal variation at x = 6D, y = 0, z = 0
93
from x=2D to x=8D is evident. The peak velocity in the water stream remains fairly constant between these locations however, the thickness of the water stream as seen on the centerplane is continually reducing. These geometry dependent flow dynamic effects are similar to those in Case I and are explained in detail in Sec.4.6. The changes in local composition and temperature are reflected in changes in the local density and viscosity as shown in Fig.4-19(b)&(c) respectively. These variations in local fluid properties are stronger than in Case I due to the larger temperature difference between the two streams. The net effect of the flow acceleration of the n-decane rich stream and the density, viscosity variations on the local flow Reynolds number is shown in Fig.4-19(d). Here, we notice a stark difference to Case I. We see that the local flow Re in the shear mixing layer increases from x = 2D to x = 8D, so that the peak Re in the shear layer is above 1500 at x = 8D. This local Re increase in the shear layer will be shown to be mainly due to the reduction in the viscosity of the fluid mixture due to the heating up of the n-decane component in Sec.4.5. At these higher values of local Re, the flow becomes unstable to perturbations leading to the growth of instabilities and waves followed by the rollup of the shear layer and finally, the collapse of these structures into small scale turbulence. Fig.4-20 and Fig.4-21 show the Yd and ωx contours respectively, on cross-sections of the tee at different downstream locations at three different time instants and also the mean field. Up to x = 4D, fluid mixing over the cross-section is predominantly caused by the circulating action of the CVP in the body of the HC jet and the secondary flow of water around the HC jet, similar to Case I. The destabilization of the shear layer and its subsequent rollup leads to intense stretching and breakdown of the streamwise vorticity accompanied by a strong enhancement of the same. This dynamics of the vorticity field is explained in detail in Sec.4.7. The enhanced vorticity stretches the material line between the two fluids, leading to enhanced mixing as is evident in Fig.4-20(a)to(d). Fig.4-22 and Fig.4-23 show the streamwise (X) and cross-stream (Z) velocity components on the tee centerplane at different time instants (t=2.0s, t=2.4s, t=2.8s, t=3.2s roughly corresponding to 2.5, 3, 3.5 and 4 tf low−through respectively) and the 94
(a) t = 2.0s
(b) t = 2.4s
(c) t = 2.8s
(d) t = 3.2s
(e) Mean
Figure 4-17: Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field at (a) t = 2.0s (2.5 tf low−through ) (b) t = 2.4s (3 tf low−through ) (c) t = 2.8s (3.5 tf low−through ) (d) t = 3.2s (4 tf low−through ) and (e) Mean field; The white vertical lines indicate the positions x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right)
95
(a) Yd,mean
(b) Tmean (in K)
Figure 4-18: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case II of (a) mean n-decane mass fraction (Yd ) field and (b) mean temperature (T in K) field at different downstream locations (x=2D, x=4D, x=6D, x=8D) 96
(a) Umean,mag (in m/s)
(b) ρmean (in kg/m3 )
(c) µmean (in Pa-s)
(d) Remean = ρmean Umean,mag D/µmean
Figure 4-19: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Case II of (a) mean velocity magnitude (Umean,mag in m/s) (b) mean density (ρmean in kg/m3 ) (c) mean dynamic viscosity (µmean in Pa-s) and (d) mean Reynolds number (Remean = ρmean Umean,mag D/µmean ) at different downstream locations (x=2D, x=4D, x=6D, x=8D)
97
(a) t = 2.4s
(b) t = 2.8s
(c) t = 3.2s
(d) Mean
Figure 4-20: Contours of Yd for simulated Case II at different downstream crosssections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) at (a) t = 2.4s (3 tf low−through ) (b) t = 2.8s (3.5 tf low−through ) (c) t = 3.2s (4 tf low−through ) and (d) Mean
98
(a) t = 2.4s
(b) t = 2.8s
(c) t = 3.2s
(d) Mean
Figure 4-21: Streamwise vorticity (ωx in s−1 ) contours for simulated Case II at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) at (a) t = 2.4s (3 tf low−through ) (b) t = 2.8s (3.5 tf low−through ) (c) t = 3.2s (4 tf low−through ) and (d) Mean
99
mean field. As the stream of water enters the joint region, the jet of HC entering from the top creates an effective area reduction to the flow of water, thus accelerating the fluid in the water stream casuing the peak fluid velocity in the stream to reach as high as 2.5 times the average value at the water inlet. This effect is clearly observable in the mean streamwise velocity plot in Fig.4-22(e). In Fig.4-23, we see that this accelerated water stream has a vertically upward velocity. This is due to the motion of the water from the stream at the bottom, around the HC jet, towards the space vacated near the top of the pipe cross-section by the downward motion of the HC jet. These flow effects are expalined in detail in Sec.4.6. At around x=5D, the water-HC shear layer becomes unstable and develops the first signs of waviness as can be seen in the instantaneous streamwise velocity plots of Fig.4-22(a)to(d). Thereafter, the shear layer undergoes rollup and vortices are continuously shed from it. The waves on the shear layer induce vertical velocity fluctuations in the fluid near the shear layer of magnitude comparable to the vertical velocity of the HC jet in the vertical pipe. This is due to the fact that the induced vertical velocity fluctuations are components of the streamwise velocity of the water stream at the shear layer which is much larger than the velocity of the HC stream as can be seen from the inlet flow conditions in Table 4.1. The eddies on the shear layer and those shed from it, also cause the velocity profile in the downstream section of the tee to become more fuller than in the laminar Case I.
4.4
Quantification of mixing
In order to make a quantitative comparative assessment of the mixing performance at different flow Reynolds numbers, we define a parameter called the ’mixing quality’ (φ) for both the species and thermal transport which represents the degree of mixing achieved at a particular cross-section downstream of the mixing joint (see Eq. 5.1). This definition is based on the ratio of the mass flux weighted mean square deviation of the mass fraction/temperature from the well mixed value over the entire crosssection to that over the inlet faces. A mixing quality of 1 indicates completely mixed 100
(a) t = 2.0s
(b) t = 2.4s
(c) t = 2.8s
(d) t = 3.2s
(e) Mean
Figure 4-22: Contour plots on the centerplane (y=0 plane) of the streamwise velocity (u in m/s) field at (a) t = 2.0s (2.5 tf low−through ) (b) t = 2.4s (3 tf low−through ) (c) t = 2.8s (3.5 tf low−through ) (d) t = 3.2s (4 tf low−through ) and (e) Mean field; The white vertical lines indicate the positions X=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right)
101
(a) t = 2.0s
(b) t = 2.4s
(c) t = 2.8s
(d) t = 3.2s
(e) Mean
Figure 4-23: Contour plots on the centerplane (y=0 plane) of the vertical velocity (w in m/s) field at (a) t = 2.0s (2.5 tf low−through ) (b) t = 2.4s (3 tf low−through ) (c) t = 2.8s (3.5 tf low−through ) (d) t = 3.2s (4 tf low−through ) and (e) Mean field; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right)
102
fluid and a mixing quality of 0 indicates completely unmixed fluid.
φspecies,x
v R u u S ρux (Yd − Yd,f ullymixed )2 .dAs Ainlets = 1 − tR x ρux (Yd − Yd,f ullymixed )2 .dAs Asx Sinlets
(4.3)
v R u u S ρux (T − Tf ullymixed )2 .dAs Ainlets = 1 − tR x ρux (T − Tf ullymixed )2 .dAs Asx Sinlets
(4.4)
φthermal,x
where, Yd,f ullymixed = 0.5. Tf ullymixed is 737.68 K in Case I and 826.85 K in Case II. The well-mixed temperature was calculated by performing an energy balance on the tee reactor. The well-mixed temperature is lower than the mean of the two inlet temperatures because of the positive heat of mixing of the water n-decane non-ideal mixture at these conditions. This means that the enthalpy of the water n-decane mixture is greater than the mass weighted sum of the enthalpies of the individual components at the same temperature and pressure. This is possibly due to the repulsive forces between the water and HC molecules. As a result, the water n-decane mixing process is endothermic. Fig.5-16(a) shows the variation of φspecies for all four simulated cases I-IV. Fig.5-16(b) shows the variation of φthermal for the simulated cases I and II. It is evident from the plots that the rate of mixing is very similar in Cases I, III and IV since the physical mixing mechanism in all three cases is the same, dominated by the action of the CVP in the HC jet and molecular diffusion far downstream (where the CVP is very weak). The rate of mixing in Case II also follows a similar rate until x = 5D, when the flow in the tee is steady. In contrast, downstream of x = 5D, the onset of waves and rollup of the mixing layer in Case II leads to intense streamwise vortical structures over the tee cross-section as shown in Sec.4.7. This leads to a significant enhancement in the rate of mixing after x = 5D as seen in Fig.5-16(a) and Fig.5-16(b). Far downstream (beyond x=16D), as these vortical structures start to weaken due to vorticity diffusion, the rate of mixing starts to plateau.
103
(a) φspecies
(b) φthermal
Figure 4-24: φ v/s x/D: Variation along the length of the tee of the (a) the species mixing quality (φspecies ) for simulated Cases I-IV and (b) the thermal mixing quality (φthermal ) for simulated Cases I and II 104
4.5
Impact of near-critical property variations: Cause of shear layer instability
In order to evaluate the impact of the variation of fluid properties of real-fluid mixtures, for a 100K temperature difference between the streams, on the flow and mixing dynamics in a tee mixer, we performed a simulation of mixing of water and n-decane with the same inlet conditions as in Case I above, but with the fluid properties of each pure component held fixed at the inlet values. The corresponding inlet conditions and fluid properties are given under Case III in Table 4.1. Fig.4-25 shows the contour plots of the streamwise vorticity on cross-sections of the tee at different downstream locations for both Case I (top) and Case III (bottom) and Fig.4-26 shows the contour plots of the n-decane mass fraction on cross-sections of the tee at different downstream locations for the same. The two sets of plots do not show any noticeable differences indicating that the density and transport property variations for a 100K temperature difference between the two streams does not significantly impact the flow field and mixing. Table 4.2 shows the density and viscosity of the pure components, water and n-decane at the two extreme temperatures in the domain at 25MPa (the nominal pressure in the tee mixer). Fig.4-27(a)&(b) show the variation of the mixture density and dynamic viscosity along the vertical centerline for the two cases I and III at x = 5D. The temperature variations in Case I lead to an increase in mixture density in the n-decane rich stream and within the mixing layer compared to Case III. The density of the water component of the mixture increases as it cools down (to temperatures lower than the inlet value of 800K) while that of the n-decane component decreases as it heats up (to temperatures higher than the inlet value of 700K). However, due to the mixture density being determined by the local mole fractions of the components, the effect of temperature on the water component density tends to dominate, resulting in the observed increase in mixture density due to T variations throughout the domain (except in the bulk of the water stream) compared to Case III. This increase in density due to T variations is mostly less than 10% of the local mixture density. Fig.4-27(b) shows that the 105
2D
4D
8D
16D
22D
2D
4D
8D
16D
22D
Figure 4-25: Streamwise vorticity (ωx in s−1 ) contours for simulated Case I (top) and Case III (bottom) at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D)
mixture viscosity is lower in the n-decane rich stream and the mixing layer in Case I compared to Case III. This effect can be understood by recalling (from Fig.5-18 and Table 4.2) that the viscosity of the water component decreases with decreasing temperature and the viscosity of the n-decane component decreases with increasing temperature at these conditions. Therefore, both the heating up of n-decane and the cooling of water leads to a decrease in mixture viscosity compared to Case III. The reduction in mixture viscosity due to T variations is below 10% of the local mixture viscosity. The net effect of the density and viscosity changes due to T variation on the local Reynolds number can be observed in Fig.4-28. The Re within the shear layer in Case I is greater than in Case III (both increase in ρ and decrease in µ lead to an increase in Re). However, the increase in Re is small and the shear layer remains stable at these values of the local Re. In order to evaluate the impact of the variation of the fluid mixture properties, for the larger ∆T between the streams (300K), on the flow and mixing dynamics in 106
2D
4D
8D
16D
22D
2D
4D
8D
16D
22D
Figure 4-26: n-decane mass fraction (Yd ) contours for simulated Case I (top) and Case III (bottom) at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D)
the tee mixer, we performed a simulation of mixing of water and n-decane with the same inlet conditions as in Case II above, but with the fluid properties of each pure component held fixed at the inlet values. The corresponding inlet conditions and fluid properties are given under Case IV in Table 4.1. As seen in Fig.5-13, the shear layer remains stable and the flow remains steady in this case. As a result, mixing occurs mainly due to the circulating action of the CVP in the HC jet and the flow around the jet from the water stream at the bottom to the top of the pipe, very similar to Case I. This can be seen in Fig.5-14 which shows contour plots of Yd on cross-sections of the tee. Clearly, for a ∆T of 300K between the two streams, the changes in local density and viscosity due to temperature variations are large enough to influence the stability characteristics of the flow. Fig.4-31(a)&(b) show the variation of the mixture density and dynamic viscosity along the vertical centerline for the two cases II and IV at x = 5D, the location at which the shear layer is observed to become unstable in Case II. The temperature 107
(a) ρ (in kg/m3 )
(b) µ (in Pa-s)
Figure 4-27: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases I and III of (a) density (ρ in kg/m3 ) and (b) dynamic viscosity (µ in Pa-s) at x=5D
108
Figure 4-28: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases I and III of the local Reynolds number at x=5D
109
Figure 4-29: Yd contours on the centerplane for simulated Case IV; The white vertical lines indicate the positions x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (left to right)
variations in Case II lead to an increase in mixture density in the n-decane rich stream and within the mixing layer compared to Case IV due to reasons similar to Case I. Since the ∆T between the streams is 3 times that in Case I, the water within the mixing layer is cooled to a greater degree from its inlet temperature value. This results in density increases as much as 20% of the local density within the mixing layer as compared to Case IV. Fig.5-18 and Table 4.2 indicate that the viscosity of the n-decane component can reduces more than 100% as it is heated from 700K (ndecane inlet T in case II) to 1000K (water inlet T in Case II). The variation of water viscosity with temperature is very small in comparison. In accordance with this, the local viscosity within the mixing layer is lower by as much as 50% in Case II compared to Case IV due to the heat up of the n-decane component. This large increase in mixture density and decrease in mixture viscosity within the mixing layer due to the cooling of water and the heating of n-decane respectively, leads to an increase of the local flow Reynolds number as shown in Fig.4-32 with the peak value of Re reaching close to 1500 in case II leading to the shear layer instability. In contrast, the local Re remains below 1300 within the mixing layer in Case IV just as in cases I and III and the shear layer remains stable. Thus, though the inlet Re number in Cases I and II are the same, larger local variations in density and viscosity in Case II due to the larger ∆T cause the local Re to increase to values in the unstable flow regime. 110
2D
4D
6D
8D
16D
Figure 4-30: Yd contours for simulated Case IV at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right)
4.6
Dynamics of flow in the cylindrical tee mixer
In Sec.4.2, we saw that mixing in the cylindrical tee mixer under laminar conditions is predominantly caused due to the advective action of a CVP in the HC jet and the flow of the water stream around this jet towards the top of the pipe. The formation of a CVP in transverse jets has been studied by Marzouk et al. [26]. Schlegel et al. [43] investigated the CVP formed in transverse jets and the contributions of the wall boundary layer to its formation. Schlegel et al. [42] also studied the impact of the CVP formed in high Re reactive transverse jets on the flame front. All these studies illustrate the strong impact of the CVP on the flow dynamics of a jet in a cross-flow (the HC stream entering the joint region from the vertical pipe bears resemblances to this flow configuration). In this section, we investigate the origin of the CVP in the HC jet and other flow features in the mixing tee and study the dynamics of the same. Fig.4-33 shows contour plots of the streamwise vorticity component (ωx ) on crosssections of the tee at different locations withing the mixing joint and slightly downstream of it (x=-0.5D, x=-0.25D, x=0, x=0.25D, x=0.5D, x=D, x=2D). Fig.4-33(a) shows a snapshot of the water stream just before it comes into contact with HC in the tee joint. Fig.4-33 (a),(c) and (d) are snapshots of the streamwise vorticity field within the mixing joint, where the two streams first come into contact. Fig.4-33 (e) 111
(a) ρ (in kg/m3 )
(b) µ (in Pa-s)
Figure 4-31: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases II and IV of (a) density (ρ in kg/m3 ) and (b) dynamic viscosity (µ in Pa-s) at x=5D 112
Figure 4-32: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases II and IV of the local Reynolds number at x=5D
113
and (f) are in the downstream pipe. These plots illustrate the process of inception of the CVP within the bulk of the HC jet entering from the top. In Fig.4-33(b), we can see the vorticity in a layer of the HC jet in the vertical pipe (anticlockwise (red) on the right and clockwise (blue) on the left). This vorticity (labeled (C)) was created due to the action of the vertical pipe wall on the hydrocarbon stream. However, the flow itself is straight down the vertical pipe without any circulation. As the HC jet flows into the cross-stream of water, a part of this vorticity gets oriented in the streamwise direction (X direction). Due to the absence of the viscous stress from side walls after the HC jet has plunged out of the vertical pipe, this streamwise vorticity leads to circulation of the fluid (anticlockwise (red) on the right and clockwise (blue) on the left) and the inception of the CVP (labeled (D) in Fig.4-33(e)). The strength of the CVP continually decreases as it is advected downstream due to vorticity diffusion and is less than 10% of its peak value beyond x=10D. Fig.4-34 shows the transverse velocity vectors (the velocity vector component along the constant X plane) colored by their magnitude on constant X planes at two X locations. At x=-0.5D, the water stream just comes in contact with the HC stream entering the mixing joint from the vertical pipe. In Fig.4-34(a) we can see the HC jet entering from the top and pushing the stream of water downwards. This downward motion of the water leads to the creation of the vorticity marked (E) in Fig.4-33(b) and also causes the water stream to accelerate in the longitudinal direction (X direction) in order to prevent the accumulation of mass near the bottom. This longitudinal acceleration of the water stream close to the mixing joint region is clearly visible in the contour plot of the streamwise velocity component in Fig.4-35. A secondary flow of HC, around the main HC jet, into the water stream is also observed in the mixing joint region (x=-0.5D to x=0.5D) in Fig.4-34(a). This secondary HC flow is responsible for the small pockets of well-mixed fluid at the peripheries of the HC jet, labeled (G) in Fig.4-10(top). Beyond x=0.5D (the end of the mixing joint), as the HC jet falls further downwards into the water stream, it vacates space near the top of the pipe cross-section. This causes fluid from the bulk of the water stream to rush upwards around the HC jet 114
(a) -0.5D
(b) -0.25D
(c) 0D
(d) 0.25D
(e) 0.5D
(f) 1D
(g) 2D
Figure 4-33: Streamwise vorticity (ωx in s−1 ) contours for simulated Case I (Rew,in = 500, Tw,in = 800K, Td,in = 700K) at different downstream cross-sections (x=-0.5D, x=-0.25D, x=0, x=0.25D, x=0.5D, x=1D, x=2D)
115
(a) 0.25D
(b) 2D
Figure 4-34: √ Transverse velocity vectors (velocity component along constant X planes, utrans = v 2 + w2 ) for simulated Case I (Rew,in = 500, Tw,in = 800K, Td,in = 700K) at two downstream cross-sections (x=0.25D, x=2D) Note: The length of the vectors is not representative of the magnitude in the figures, the color of the vectors represents the magnitude
into the space at the top and can be observed in Fig.4-34(b). This flow of water to the top part of the tee is further aided by the circulating action of the CVP. The addition of mass to the n-decane rich stream at the top causes the acceleration of this stream as observed in Fig.4-9(a) and Fig.4-19(a) in spite of the slow growth of the shear layer.
4.7
Vorticity dynamics
The streamwise vorticity, ωx plays a significant role in mixing as explained in Sec.4.6 and hence, we need to look at the evolution of this component of vorticity and the physical mechanisms affecting it. Eq.5.3 below is the vorticity transport equation 116
Figure 4-35: (left) Contour plots on the centerplane (y=0 plane) of the velocity mag1/2 nitude (Umag = (u2 + v 2 + w2 ) ) in the joint region; (right) n-decane mass fraction profile along the vertical line at x = 2D (the position indicated by the white line in the left figure) obtained by applying the curl operator to the momentum transport equation: ∇.τ 1 ∂ω + (u.∇) ω = (ω.∇) u − ω (∇.u) + 2 (∇ρ × ∇P ) + ∇ × | {z } | {z } ρ ∂t ρ | {z } | {z } A B C
(4.5)
D
where, term (A) is the enhancement of vorticity due to the stretching and tilting of vortex lines (ω˙ st ), term (B) is the vorticity enhancement due to the dilatation of fluid elements (ω˙ d ), term (C) is the baroclinic vorticity generation (ω˙ bg ) due to the presence of a density gradient perpendicular to a gradient in pressure and term (D) represents the viscous term (ω˙ v ) which includes in it, the vorticity diffusion (µ∇2 (ω)). Fig.4-36 shows the mean field of the above-mentioned streamwise vorticity source terms on different downstream cross-sections for simulated Case II. The formation of the CVP in the HC jet occurs in a similar fashion to Case I. Thereafter, till x=5D, the vorticity diffusion dominates as seen in Fig.4-36(c). This leads to a decrease in the strength of the CVP from the joint region to about x = 5D. This CVP decay process is clearly visible in Fig.4-21. The destabilization of the shear layer near x=5D triggers 117
stretching and breaking of the CVP into smaller coherent structures. This stretching and breakdown of the vortices can be clearly observed in Fig.4-21. There is also a significant enhancement of these streamwise vortices due to the stretching of fluid elements (ω˙ st,x ). Fig.4-38 shows the ωx profile along the z-direction at y=0.0004 and different downstream locations (x=4D, x=5D, x=6D, x=8D). The sudden increase in the magnitude of ωx near z = 0.0001 from x=5D to x=8D can be seen in this plot. This streamwise vorticity enhancement is due to two effects: (a) the strain induced in the streamwise direction due to the waviness and rollup of the shear layer resulting in vorticity enhancement due to fluid stretching as seen in Fig.4-37(a) and (b) baroclinic generation of vorticity resulting from the average effect over time of strong instantaneous local density and pressure gradients due to the unsteady motions of the shear layer, as seen in Fig.4-37(b). The enhanced vorticity stretches the material line between the two fluids, leading to enhanced mixing beyond x=5D as is evident in Fig.4-20(a)to(d). As these vortices travel further downstream (beyond x=10D), the viscous term starts to dominate again and their strength continually reduces.
4.7.1
Species and thermal transport enhancement due to fluctuations
The advective fluxes of the conserved scalars like Yd and h can be broken down into contributions due to the mean field and the fluctuations as below, in Eq.5.4: u¯i φ = u¯i φ¯ + u0i¯φ0 |{z} |{z} A
(4.6)
B
where, φ denotes a conserved scalar (like Yd or h), (¯) denotes mean value and (’) denotes the fluctuation. Term (A) represents the advective flux contribution due to the mean field; Term (B) represents the advective flux contribution due to the fluctuations in the velocity and scalar fields. In Sec.4.7, we saw that the rollup of the shear layer triggers intense stretching and breakdown of the streamwise vorticity into smaller and stronger vortical struc118
(a) stretching and tilting
(b) baroclinic
(c) viscous diffusion
(d) dilatation
Figure 4-36: Contours of the mean field for simulated Case II at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right) of (a) ω˙ st,x in s−2 (b) ω˙ bg,x in s−2 (c) ω˙ v,x in s−2 and (d) ω˙ d,x in s−2
119
(a) ω˙ st,x,mean (in s−2 ), Case II
(b) ω˙ bg,x,mean (in s−2 ), Case II
Figure 4-37: Profiles along the vertical line at y=0.0004 (y=D/6) for simulated Case II of (a) ω˙ st,x,mean (in s−2 ) and (b) ω˙ bg,x,mean (in s−2 ) at different downstream locations (x=4D, x=5D, x=6D, x=8D) 120
Figure 4-38: Profiles along the vertical line at y=0.0004 (y=D/6) for simulated Case II of ωx,mean (in s−1 ) at different downstream locations (x=4D, x=5D, x=6D, x=8D)
121
2D
4D
6D
8D
16D
0 1/2 Figure 4-39: RMS fluctuation of Yd (Y¯d 2 ) contours for simulated Case II at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right)
tures. This fluctuating streamwise vorticity and the resultant fluid circulation produces strong fluctuations in the scalar variables like Yd over the cross-section of the tee. This can be seen in Fig.4-39 which shows the RMS fluctuations of Yd on different cross-sections of the tee for Case II. Fig.4-40 shows the advective flux contributions due to w on the tee cross-sections. Fig.4-41 shows the advective flux contributions due to v on the tee cross-sections. The fluctuations in the velocity components are normalized with uref = 0.14 m/s, the water inlet average velocity. It is clear that the contribution due to the fluctuations contribute significantly to the transport of species over the cross-section of the tee. The resultant enhancement in mixing can be visualized in Fig.4-20 which shows the instantaneous and mean Yd field on different cross-sections downstream of the mixing joint.
4.8
Summary
In this chapter, simulations of mixing of supercritical water and a model hydrocarbon (n-decane), under fully-miscible conditions, in a cylindrical tee mixer geometry for a water inlet Reynolds number of 500, for two different inlet temperatures of the water stream of 800K (Case I) and 1000K (Case II) were discussed. The mass flow rate of the n-decane stream was equal to that of the water stream and its inlet temperature was 700K in both cases. 122
(a) wmean Yd,mean /uref
(b) w0¯Yd /uref 0
Figure 4-40: Contours for simulated Case II at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=16D from left to right) of (a) 0 wmean Yd,mean /uref (b) w0¯Yd /uref ; uref = 0.14 m/s
123
(a) vmean Yd,mean /uref
(b) v 0¯Yd /uref 0
Figure 4-41: Contours for simulated Case II at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=16D from left to right) of (a) vmean Yd,mean /uref 0 (b) v 0¯Yd /uref ; uref = 0.14 m/s
124
It was found that the flow downstream of the mixing joint remained laminar in case of the 100K temperature difference between the streams (Case I). Most of the mixing and heat transport occurs due to the circulating action of a counter-rotating vortex pair (CVP) in the hydrocarbon jet formed due to the reorientation of the vorticity pre-existing in the hydrocarbon stream flowing through the vertical pipe. This CVP gets progressively weaker as it flows downstream due to vorticity diffusion and species and heat transport is dominated by molecular diffusion over small length scales. Consequently, the mixing rate plateaus in the far downstream region of the tee mixer. Comparison with mixing of water and n-decane at the same inlet conditions but without variations of the physical properties with temperature (Case II) suggests that the near-critical property variations, for a 100K temperature difference between the streams, have a negligible impact on the flow field and mixing. For a 300K temperature difference between the two streams (Case II), the waterHC shear layer becomes unstable and starts to form waves near x=5D downstream of the mixing joint center. Further downstream, the shear layer rolls up and vortices are shed from it. The onset of instability in the shear layer also triggers the stretching and breakdown of the CVP in the body of the hydrocarbon jet. These smaller, stronger vortical structures cause faster advective transport of species and enthalpy over the cross-section of the tee leading to a significant enhancement in mixing. This manifests as a jump in the mixing rate beyond the location of onset of instability and a thickening of the mean mixing layer. However, water n-decane mixing under identical inlet conditions but with constant physical properties (Case IV), shows a stable shear layer with the flow reaching steady state. In Case II, the temperature of the water component changes strongly (from 1000K to 700K) within the mixing layer and this cooling of the water component leads to an increase in the density of the fluid mixture within the mixing layer compared to Case IV. The cooling of water and heating up of n-decane within the mixing layer also leads to significant reduction in the mixture viscosity (as much as 50%) within the layer compared to Case IV. Both, the density increase and viscosity decrease lead to an increase in the local flow Reynolds number within the mixing layer with the peak 125
Re value reaching close to 1500. The shear layer becomes unstable to perturbations in the flow at these higher local Re. This results in the growth of perturbations in the shear layer and the consequent waviness and rollup of the layer followed by the collapse of these structures into small-scale turbulence.
126
Chapter 5 Mixing at intermediate Re: Impact of Re In this chapter, the mixing of supercritical water and n-decane is simulated for four different values of the Reynolds number at the water inlet (Rew,in = 500, 600, 700 and 800) as shown in Table 5.1. In all cases, the mass flow rates through the water and n-decane inlets is the same (a design condition found to be favorable to the process from the chemical kinetics point of view). The nominal pressure in the tee mixer was ≈ 25MPa (ensured by fixing the back-pressure at the outlet) and the n-decane inlet temperature was 700 K, which is sufficiently above the Upper Critical Solution Temperature (UCST) of the water n-decane system (≈ 632 K) thereby ensuring complete miscibility of the two components. The water inlet temperature was 800 K. The effect of gravity forces is not included in these simulations so as to allow us to study the impact of geometry, Reynolds number and property variation effects on the flow field and mixing dynamics.
5.1
Validation: Grid convergence tests
In order to verify the adequacy of the chosen spatial resolution and verify the numerics of the code, a grid independence study was performed for the largest Re case, Case IV (Rew,in = 800). Fig.5-1 shows the mean streamwise velocity profile along the vertical 127
Case
Rew,in
Red,in
I II III IV
500 600 700 800
204 245 286 327
m ˙ w,in (kg/hr) 0.11 0.14 0.16 0.18
m ˙ d,in (kg/hr) 0.11 0.14 0.16 0.18
uw,in (m/s) 0.0832 0.0998 0.1164 0.1331
ud,in (m/s) 0.0144 0.0173 0.0202 0.0230
Tw,in (K) 800 800 800 800
Td,in (K) 700 700 700 700
Property variations yes yes yes yes
Table 5.1: Inlet flow conditions for the cases simulated in this study direction for Case IV (the highest Re simulated) at three different axial locations (x=4D, x=6D and x=8D) for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm). The averages are taken over 4 tf low−through time periods after the initial transient of 2 tf low−through time periods. Grid-independence is attained for a mesh resolution of 0.06 mm in this case. The velocity field thus seems to be well resolved. Fig.5-2 and Fig.5-3 show the mean n-decane mass fraction and mean temperature profiles respectively, at three different axial locations (x=4D, x=6D and x=8D) for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm). These figures indicate that the mixing layer is well resolved in terms of the scalar fields. The mixing layer thickness from these plots can be seen to be around 0.5 mm around x = 6D. We will later see that in Case IV, we observe that the shear layer becomes unstable. The wavelengths of the unstable perturbations or waves on the shear layer are multiples of the mixing layer thickness as stated in literature, for example [30]. Since the grid size of 0.06 mm gives around 10 cells within the mixing layer, we can conclude that the spatial resolution is sufficient to capture the scales of instability in the shear layer. This gives us reasonable confidence in the consistency of the numerical methods and code. Unfortunately, there is no experimental data available in the literature for mixing of near-critical fluids in a cylindrical tee reactor geometry.
5.2
Mixing Results
Fig.5-4 shows the n-decane mass fraction (Yd ) on the tee centerplane for Case I and Case II. At these flow Re, the water-HC shear layer remains stable and the flow in the tee is laminar with a steady-state being reached. Fig.5-5 and Fig.5-6 show 128
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 5-1: u (in m/s) v/s z (in m) for simulated Case IV for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D 129
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 5-2: Yd v/s z (in m) for simulated Case IV for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D
130
(a) x = 4D
(b) x = 6D
(c) x = 8D
Figure 5-3: T (in K) v/s z (in m) for simulated Case IV for three different mesh resolutions (0.10 mm, 0.08 mm and 0.06 mm) at (a) x = 4D (b) x = 6D and (c) x = 8D 131
contour plots of the Yd and T fields respectively on different cross-sections of the tee downstream of the mixing joint, illustrating more clearly, the features of the multidimensional mixing process in the tee mixer. Fig.5-7 shows the streamwise vorticity field (ωx ) on these cross-sections for both cases. It is clear that most of the mixing over large length scales is due to the circulating action of a counter-rotating vortex pair (CVP) within the body of the HC jet entering through the vertical pipe. This CVP is formed due to the re-orientation of the pre-existing vorticity in the HC jet flowing through the vertical pipe as it enters the mixing joint and travels downstream. The details of the process of inception of the CVP are presented in Sec.4.6. As the HC jet moves further downwards in the downstream pipe, it vacates space near the top of the pipe. This causes a motion of water from the stream at the bottom, around the periphery of the HC jet, towards the space at the top. This secondary flow around the HC jet is another physical mechanism leading to large scale advective mixing. Fig.5-7 also shows that the strength of the streamwise vorticity and the resultant circulating motions is higher in Case II compared to Case I. This is due to the fact that, the pre-existing vorticity in the HC jet in the vertical pipe is larger in Case II due to the larger flow velocity in this section. The re-orientation of this stronger vorticity leads to a stronger CVP in the higher Re case. The enhanced advective mixing due to the stronger streamwise vorticity can be clearly seen in the comparative plots of the Yd and T fields in Fig.5-5 and Fig.5-6. The strength of these vortices continually decreases as they travel downstream due to vorticity diffusion by viscous action, falling to less than 10% of the value just after the mixing joint by around x = 12D. Further downstream, mixing is mainly controlled by molecular diffusion. When, Rew,in = 700, the shear layer becomes unstable, with waves forming on the surface followed by the roll up of the shear layer as can be seen in the instantaneous centerplane Yd plots in Fig.5-8(a)to(c). The formation of these smaller scale flow structures on the shear mixing layer leads to faster advective transport of the scalars resulting in enhanced mixing and heat transport. This is manifested in a thickening of the mean mixing layer beyond x = 8D as can be seen in the mean Yd field plot in Fig.5-8(d). In the unsteady cases, the flow is allowed to develop for 2tf low−through to 132
(a) Case I, Rew,in = 500
(b) Case II, Rew,in = 600
Figure 5-4: Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field at steady-state for (a) Case I and (b) Case II; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=16D downstream of the center of the mixing joint (from left to right)
(a) Case I, Rew,in = 500
(b) Case II, Rew,in = 600
Figure 5-5: Contour plots of Yd on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=16D and x=22D from left to right) at steady-state for (a) Case I and (b) Case II
133
(a) Case I, Rew,in = 500
(b) Case II, Rew,in = 600
Figure 5-6: Contour plots of T (in K) on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=16D and x=22D from left to right) at steady-state for (a) Case I and (b) Case II
wipe out the initial transients and thereafter, averaging of field variables is performed over the next 4tf low−through . Statistically stationary state is typically reached within 2 − 2.5tf low−through . Since the statistics of the field variables do not change with time after about 2 to 2.5 tf low−through , we suspect some periodicity in the temporal variations of the fields. In Fig.5-8(b) a clear spatial periodicity is observed in the roll up structures on the water-HC shear layer with a wavelength close to 0.5D. Fig.5-9 shows the Fourier transform of the n-decane mass fraction field at x = 8D, y = 0, z = 0 (a location after the shear layer destabilizes). We can see a cluster of frequency modes around in the 0-10 Hz and 90-100 Hz bands which dominate the flow field fluctuations. In addition, a group of frequency modes in the 40-60 Hz range is also observed but of lesser amplitude. It is not possible to identify from these plots, an exact location at which the onset of instability occurs in this case, since the waviness in the shear layer is intermittent in time, with the shear layer being smooth at certain time instants for the same downstream location (even locations as far as x 134
(a) Case I, Rew,in = 500
(b) Case II, Rew,in = 600
Figure 5-7: Contour plots of the streamwise vorticity, ωx (in s−1 ) on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=16D and x=22D from left to right) at steady-state for (a) Case I and (b) Case II
135
= 10D). The effective location of onset of instability is identified to be around x = 6D 0
from the normalized root-mean square (RMS) velocity fluctuation (u /U ) plot on the centerplane in Fig.5-12. Fig.5-10 shows the instantaneous and mean Yd contour plots on cross-sections of the tee at different downstream locations (x=2D, x=4D, x=8D, x=12D and x=16D). Fig.5-11 shows the instantaneous and mean ωx field on the same cross-sections. Analogous to cases I and II, the flow is steady just downstream of the mixing joint (up to about x = 6D) and mixing predominantly occurs due to the circulating action of the CVP in the HC jet and the secondary flow of water around the jet. This action becomes progressively weaker as the CVP strength decreases due to vorticity diffusion as it travels downstream. Near x=8D, the unstable motions of the shear layer trigger the intense stretching and breakdown of the streamwise vorticity accompanied by an enhancement of the same as seen in Fig.5-11. The development of these stronger smaller scale vortical structures promotes faster advective transport of the scalars over the cross-section of the tee. The enhanced vorticity stretches the material surface between the two fluid streams leading to an enhancement in mixing as seen in Fig.5-10. At Rew,in = 800 (Case IV), the water-HC shear layer destabilizes earlier, around x = 5D. The unsteadiness in the shear layer is not intermittent as in Case III. The waves on the shear layer and roll up structures are larger than in Case III (Rew,in = 700) as can be seen in the instantaneous Yd plots on the centerplane in Fig.5-13(a)to(c). The small scale flow structures formed on the shear layer lead to fast advective transport of the scalars resulting in enhanced mixing. This is manifested in the thickening of the mean mixing layer beyond x = 6D and continuing far downstream as seen in the mean Yd plot on the centerplane in Fig.5-13(d). This effect is much more pronounced than in Case III illustrating the efficacy of the small scales in enhancing the mixing process. Fig.5-14 shows the instantaneous and mean Yd contour plots on cross-sections of the tee at different downstream locations (x=4D, x=6D, x=8D, x=12D and x=16D). Fig.5-15 shows the instantaneous and mean ωx field on the same cross-sections. Near x=6D, the destabilization of the shear layer triggers the intense stretching and breakdown of the streamwise vorticity accompanied by an 136
(a) t = 2.0s
(b) t = 2.2s
(c) t = 2.5s
(d) Mean
Figure 5-8: Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field for simulated Case III at (a) t = 2.0s (4 tf low−through ) (b) t = 2.2s (4.4 tf low−through ) (c) t = 2.5s (5 tf low−through ) and (d) Mean field; The white vertical lines indicate the positions x=6D, x=8D, x=10D, x=12D, x=14D, x=16D, x=18D and x=20D downstream of the center of the mixing joint (left to right)
137
Figure 5-9: Fourier transform of the Yd temporal variation at x = 8D, y = 0, z = 0 for Case III
138
(a) t = 2.0s
(b) t = 2.2s
(c) t = 2.5s
(d) Mean
Figure 5-10: Contours of Yd for simulated Case III at different downstream crosssections (x=2D, x=4D, x=8D, x=12D, x=16D from left to right) at (a) t = 2.0s (b) t = 2.2s (c) t = 2.5s and (d) Mean
139
(a) t = 2.0s
(b) t = 2.2s
(c) t = 2.5s
(d) Mean
Figure 5-11: Contours of ωx (in s−1 ) for simulated Case III at different downstream cross-sections (x=2D, x=4D, x=8D, x=12D, x=16D from left to right) at (a) t = 2.0s (b) t = 2.2s (c) t = 2.5s and (d) Mean
140
0
u /U
0
Figure 5-12: Contour plots on the centerplane (y=0 plane) for Case III of u /U , where U is the reference velcoity used for normalization taken to be the water inlet average velocity; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D, x=14D and x=16D downstream of the center of the mixing joint (left to right) enhancement of the same as seen in Fig.5-15. The development of these smaller scale vortical structures promotes faster advective transport of the scalars over the crosssection of the tee. The enhanced vorticity stretches the material surface between the two fluid streams leading to an enhancement in mixing as seen in Fig.5-14. This stretching and breakdown of the eddies into smaller structures is more vigorous as compared to Case III and is a classic feature of a turbulent flow field.
5.3
Quantification of mixing
In order to make a quantitative comparative assessment of the mixing performance at different flow Reynolds numbers, we define a parameter called the ’mixing quality’ (φ) for both the species and thermal transport which represents the degree of mixing achieved at a particular cross-section downstream of the mixing joint (see Eq. 5.1). This definition is based on the ratio of the mass flux weighted mean square deviation of the mass fraction/temperature from the well mixed value over the entire crosssection to that over the inlet faces. A mixing quality of 1 indicates completely mixed fluid and a mixing quality of 0 indicates completely unmixed fluid.
φspecies,x
v R u u S ρux (Yd − Yd,f ullymixed )2 .dAs Ainlets = 1 − tR x ρux (Yd − Yd,f ullymixed )2 .dAs Asx Sinlets 141
(5.1)
(a) t = 1.8s
(b) t = 2.05s
(c) t = 2.25s
(d) Mean
Figure 5-13: Contour plots on the centerplane (y=0 plane) of n-decane mass fraction (Yd ) field for simulated Case IV at (a) t = 1.8s (4 tf low−through ) (b) t = 2.05s (≈ 4.5 tf low−through ) (c) t = 2.25s (5 tf low−through ) and (d) Mean field; The white vertical lines indicate the positions x=4D, x=6D, x=8D, x=10D, x=12D, x=14D, x=16D and x=18D downstream of the center of the mixing joint (left to right)
142
(a) t = 1.8s
(b) t = 2.05s
(c) t = 2.25s
(d) Mean
Figure 5-14: Contours of Yd for simulated Case IV at different downstream crosssections (x=4D, x=6D, x=8D, x=12D, x=16D from left to right) at (a) t = 1.8s (b) t = 2.05s (c) t = 2.25s and (d) Mean
143
(a) t = 1.8s
(b) t = 2.05s
(c) t = 2.25s
(d) Mean
Figure 5-15: Contours of ωx (in s−1 ) for simulated Case IV at different downstream cross-sections (x=4D, x=6D, x=8D, x=12D, x=16D from left to right) at (a) t = 1.8s (b) t = 2.05s (c) t = 2.25s and (d) Mean
144
φthermal,x
v R u u S ρux (T − Tf ullymixed )2 .dAs Ainlets = 1 − tR x ρux (T − Tf ullymixed )2 .dAs Asx Sinlets
(5.2)
where, Yd,f ullymixed = 0.5. Tf ullymixed is 737.68 K. The well-mixed temperature was calculated by performing an energy balance on the tee reactor. The well-mixed temperature is lower than the mean of the two inlet temperatures because of the positive heat of mixing of the water n-decane non-ideal mixture at these conditions. This means that the enthalpy of the water n-decane mixture is greater than the mass weighted sum of the enthalpies of the individual components at the same temperature and pressure. This is possibly due to the repulsive forces between the water and HC molecules. As a result, the water n-decane mixing process is endothermic. Fig.5-16(a) shows the variation of φspecies for all four simulated cases I-IV. Fig.5-16(b) shows the variation of φthermal for the simulated cases I-IV. It is evident from the plots that the rate of mixing is very similar in Cases I and II since the dominant physical mixing mechanism in both cases is the same, dominated by the action of the CVP in the HC jet and molecular diffusion far downstream (where the CVP is very weak). The rate of mixing in Case III also follows a similar rate until x = 6D, when the flow in the tee is steady. In contrast, downstream of x = 6D, the onset of waves and rollup of the mixing layer in Case III leads to the formation of intense streamwise vortical structures over the tee cross-section as shown in Sec.5.5. This leads to a significant enhancement in the rate of mixing after x = 6D as seen in Fig.5-16(a) and Fig.5-16(b). Similarly, in Case IV the steep increase in mixing rate when the mixing layer becomes unstable near x = 5D is evident in Fig.5-16(a) and Fig.5-16(b). The streamwise vorticity that stretches the material line between the two streams and enhances mixing is much stronger in Case IV compared to Case III, due to the higher flow Re. Hence, the mixing rate after mixing layer transition is also much greater for Case IV. Far downstream (beyond x=16D), as these vortical structures start to weaken due to vorticity diffusion, the rate of mixing starts to plateau.
145
(a) φspecies
(b) φthermal
Figure 5-16: φ v/s x/D: Variation along the length of the tee of the (a) the species mixing quality (φspecies ) and (b) the thermal mixing quality (φthermal ) for simulated Cases I-IV 146
Case I to IV I to IV I to IV I to IV
Specie water water n-decane n-decane
T (K) 700 800 700 800
ρ(kg/m3 ) 125.62 84.22 486.27 424.45
µ(P a − s) 3.1773e − 05 3.3626e − 05 8.2378e − 05 5.8599e − 05
Table 5.2: Properties of water and n-decane at 25MPa and extreme temperatures of the simulations in this study
5.4
Impact of near-critical property variations: Cause of shear layer instability
It is important to note at this stage, the variation of the mixture density and viscosity with both temperature and composition. Fig.5-17 and Fig.5-18 show the variation of the mixture density (ρ in kg/m3 ) and dynamic viscosity (µ in P a − s) with temperature for the range of temperatures pertinent to this study, for different mixture compositions from pure water to pure n-decane. Table 5.2 shows the ρ and µ for the water and n-decane pure components at the extreme temperatures within the tee mixer for the cases I and II. It can be seen that the mixture density and viscosity are strong functions of the temperature under these conditions. The density of the water component can increase by around 50% as it gets cooled from 800K to 700K within the tee mixer domain while the density of the n-decane component can decrease by around 10% as it gets heated from 700K to 800K. The viscosity of water does not change appreciably in this temperature range. However, the viscosity of the n-decane component can reduce by about 30% as it gets heated from 700K to 800K. These variations in ρ and µ will influence the local Reynolds number which in turn can be expected to have a significant bearing on the local flow dynamics and stability. We intend to evaluate the impact of these property variations (with temperature) within the mixing tee on the flow field and mixing dynamics. In order to evaluate the impact of near-critical property variations in the tee mixer, an additional test simulation was performed (Case V) with inlet conditions as specified in Table 5.3. In Case V, the inlet conditions are identical to Case III (Rew,in = 700) but with the physical properties (ρ, µ, λ, D) held fixed at the inlet 147
Figure 5-17: ρ (in kg/m3 ) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane)
Figure 5-18: µ (in Pa-s) v/s T (in K) at P = 25 MPa for different mixture compositions of water and n-decane from Yd = 0 (pure water) to Yd = 1 (pure n-decane) 148
Case
Rew,in
Red,in
V
700
286
m ˙ w,in (kg/hr) 0.16
m ˙ d,in (kg/hr) 0.16
uw,in ud,in (m/s) (m/s) 0.1164 0.0202
Tw,in (K) −
Td,in (K) −
Property variations no
Table 5.3: Inlet flow conditions for the comparison test Case V simulated in this study
Figure 5-19: Yd contours on the centerplane for simulated Case V; The white vertical lines indicate the positions x=2D, x=4D, x=6D, x=8D, x=10D, x=12D and x=14D downstream of the center of the mixing joint (left to right) values. Though the variation of physical properties with temperature are not allowed in this simulations, the local variation in the mixture properties due to the mixing of water and n-decane (since water and n-decane have significantly different physical properties) are accounted for. Comparison of this case with Case III will tell us how the near-critical variations of mixture properties affect the flow and mixing. As seen in Fig.5-13, the shear layer remains stable and the flow remains steady in Case IV. As a result, mixing occurs mainly due to the circulating action of the CVP in the HC jet and the flow around the jet from the water stream at the bottom to the top of the pipe, very similar to the laminar cases I and II. This can be seen in Fig.5-20 which shows contour plots of Yd on cross-sections of the tee. Clearly, for a ∆T of 100K between the two streams, the changes in local density and viscosity due to temperature variations influence the stability characteristics of the flow. Fig.5-21(a)&(b) show the variation of the mixture density and dynamic viscosity along the vertical centerline for the two cases III and V at x = 6D, the location at which the shear layer is observed to become unstable in Case III.. The temperature variations in Case III lead to an increase in mixture density in the n-decane rich stream and within the mixing layer compared to Case V. The density of the water 149
2D
4D
8D
16D
22D
Figure 5-20: Yd contours for simulated Case V at different downstream cross-sections (x=2D, x=4D, x=8D, x=16D, x=22D from left to right) component of the mixture increases as it cools down (to temperatures lower than the inlet value of 800K) while that of the n-decane component decreases as it heats up (to temperatures higher than the inlet value of 700K). However, due to the mixture density being determined by the local mole fractions of the components, the effect of temperature on the water component density tends to dominate, resulting in the observed increase in mixture density due to T variations throughout the domain (except in the bulk of the water stream) compared to Case V. This increase in density due to T variations is mostly around 10% of the local mixture density. Fig.5-21(b) shows that the mixture viscosity is lower in the n-decane rich stream and the mixing layer in Case III compared to Case V. This effect can be understood by recalling (from Fig.5-18 and Table 5.2) that the viscosity of the water component decreases with decreasing temperature and the viscosity of the n-decane component decreases with increasing temperature at these conditions. Therefore, both the heating of ndecane and the cooling of water leads to a decrease in mixture viscosity compared to Case V. The reduction in mixture viscosity due to T variations is close to 10% of the local mixture viscosity. The net effect of the density and viscosity changes due to T variation on the local Reynolds number can be observed in Fig.5-22. The Re within the shear layer in Case III is greater than in Case V (both increase in ρ and decrease in µ lead to an increase in Re). This increase in the local Re within the shear layer makes it unstable to the growth of perturbations resulting in the shear 150
layer instability and roll up observed in Case III.
5.5
Vorticity dynamics
The streamwise vorticity, ωx plays a significant role in mixing as explained in Sec.4.6 and hence, we need to look at the evolution of this component of vorticity and the physical mechanisms affecting it. Eq.5.3 below is the vorticity transport equation obtained by applying the curl operator to the momentum transport equation: ∇.τ ∂ω 1 + (u.∇) ω = (ω.∇) u − ω (∇.u) + 2 (∇ρ × ∇P ) + ∇ × | {z } | {z } ρ ∂t ρ {z } | | {z } A B C
(5.3)
D
where, term (A) is the enhancement of vorticity due to the stretching and tilting of vortex lines (ω˙ st ), term (B) is the vorticity enhancement due to the dilatation of fluid elements (ω˙ d ), term (C) is the baroclinic vorticity generation (ω˙ bg ) due to the presence of a density gradient perpendicular to a gradient in pressure and term (D) represents the viscous term (ω˙ v ) which includes in it, the vorticity diffusion (µ∇2 (ω)). In order to understand the impact of the shear layer instability on the streamwise vorticity field, which plays a crucial role in advective transport, we look at the contributions to the enhancement/production of streamwise vorticity due to ω˙ st,y , ω˙ bg,y , ω˙ v,y and ω˙ d,y in Fig.5-23. We look at Case IV in particular, but the same general phenomena are observed in Case III as well. The formation of the CVP in the HC jet occurs in a similar fashion to Case I. Thereafter, till x=5D, the vorticity diffusion dominates leading to a reduction in the strength of the CVP as seen in Fig.5-23(c) and Fig.5-15. The destabilization of the shear layer near x=5D triggers stretching and breaking of the CVP into smaller coherent structures with strong vortical structures spreading over most of the tee cross-section. This stretching and breakdown of the vortices can be clearly observed in Fig.5-15. There is also a significant enhancement of these streamwise vortices due to the stretching of fluid elements. This is due to the strain induced in the streamwise direction due to the waviness and rollup of the shear layer. Downstream of x=5D, the unsteady motion of the shear layer also leads to 151
(a) ρ (in kg/m3 )
(b) µ (in Pa-s)
Figure 5-21: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases III and V of (a) density (ρ in kg/m3 ) and (b) dynamic viscosity (µ in Pa-s) at x=6D 152
Figure 5-22: Profiles along the vertical centerline (y=0 plane) at steady-state for simulated Cases III and V of the local Reynolds number at x=6D
153
strong density and pressure gradients giving rise to a significant baroclinic generation of vorticity as seen in Fig.5-23(b). As these vortices travel further downstream, the viscous diffusion starts to dominate and their strength continually reduces. Even so, the streamwise vorticity is strong even at far downstream locations (as high as 100s−1 at x=16D) in this case.
5.6
Species and thermal transport enhancement due to fluctuations
The advective fluxes of the conserved scalars like Yd and h can be broken down into contributions due to the mean field and the fluctuations as below, in Eq.5.4: u¯i φ = u¯i φ¯ + u0i¯φ0 |{z} |{z} A
(5.4)
B
where, φ denotes a conserved scalar (like Yd or h), (¯) denotes mean value and (’) denotes the fluctuation. Term (A) represents the advective flux contribution due to the mean field; Term (B) represents the advective flux contribution due to the fluctuations in the velocity and scalar fields. In Sec.5.5, we saw that the rollup of the shear layer triggers intense stretching and breakdown of the streamwise vorticity into smaller and stronger vortical structures. This fluctuating streamwise vorticity and the resultant fluid circulation produces strong fluctuations in the scalar variables like Yd over the cross-section of the tee. This can be seen in Fig.5-24 which shows the RMS fluctuations of Yd on different cross-sections of the tee for Case IV. Fig.5-25 shows the advective flux contributions due to w on the tee cross-sections. Fig.5-26 shows the advective flux contributions due to v on the tee cross-sections. The fluctuations in the velocity components are normalized with uref = 0.13 m/s, the water inlet average velocity. It is clear that the contribution due to the fluctuations contribute significantly to the transport of species over the cross-section of the tee. The resultant enhancement in mixing can be 154
(a) stretching and tilting term
(b) baroclinic term
(c) viscous term
(d) dilatation term
Figure 5-23: Contours of the mean field for simulated Case IV at different downstream cross-sections (x=4D, x=6D, x=8D, x=10D, x=16D from left to right) of (a) ω˙ st,y in s−2 (b) ω˙ bg,y in s−2 (c) ω˙ v,y in s−2 and (d) ω˙ d,y in s−2
155
2D
4D
6D
8D
16D
0 1/2 Figure 5-24: RMS fluctuation of Yd (Y¯d 2 ) contours for simulated Case IV at different downstream cross-sections (x=2D, x=4D, x=6D, x=8D, x=16D from left to right)
visualized in Fig.5-14 which shows the instantaneous and mean Yd field on different cross-sections downstream of the mixing joint.
5.7
Summary
In this chapter, simulations of mixing of supercritical water and a model hydrocarbon (n-decane), under fully-miscible conditions, in a cylindrical tee mixer geometry for four values of the water inlet Reynolds number in the range of 500 to 800 were discussed. The mass flow rate of the n-decane stream was equal to that of the water stream. The inlet temperatures of water and n-decane were 800K and 700K respectively. It was found that the flow downstream of the mixing joint remained laminar in Cases I and II (Rew,in = 500 and Rew,in = 600 respectively). Most of the mixing and heat transport occurs due to the circulating action of a counter-rotating vortex pair (CVP) in the hydrocarbon jet formed due to the reorientation of the vorticity pre-existing in the hydrocarbon stream flowing through the vertical pipe. This CVP gets progressively weaker as it flows downstream due to vorticity diffusion and species and heat transport is dominated by molecular diffusion over small length scales. Consequently, the mixing rate plateaus in the far downstream region of the tee mixer. When the Reynolds number at the water inlet is increased to 700 and the n-decane 156
(a) wmean Yd,mean /uref
(b) w0¯Yd /uref 0
Figure 5-25: Contours for simulated Case IV at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=10D from left to right) of (a) 0 wmean Yd,mean /uref (b) w0¯Yd /uref ; uref = 0.13 m/s
157
(a) vmean Yd,mean /uref
(b) v 0¯Yd /uref 0
Figure 5-26: Contours for simulated Case IV at cross-sections downstream of location of onset of instability (x=6D, x=8D, x=10D from left to right) of (a) vmean Yd,mean /uref 0 (b) v 0¯Yd /uref ; uref = 0.13 m/s
158
inlet Re also increased correspondingly to maintain the same mass flow rate (Case III), the shear layer between the water and n-decane streams is found to destabilize near x = 6D downstream of the mixing joint. Further downstream, the shear layer rolls up and vortices are shed from it. The onset of instability in the shear layer also triggers the stretching and breakdown of the CVP in the body of the hydrocarbon jet. These smaller, stronger vortical structures cause faster advective transport of species and enthalpy over the cross-section of the tee leading to a significant enhancement in mixing. This manifests as a jump in the mixing rate beyond the location of onset of instability and a thickening of the mean mixing layer. At Rew,in = 800 (Case IV) the unsteady small scale flow structures in the mixing layer and the consequent flow field fluctuations due to them are much stronger. The stretching and breakdown of the CVP in this case, is accompanied by stronger streamwise vorticity enhancement resulting in much faster mixing compared to Case III. However, water n-decane mixing under identical inlet conditions to Case III but with constant physical properties (Case V), shows a stable shear layer with the flow reaching steady state. In Case III, the temperature of the water component changes from 800K to 700K within the mixing layer and this cooling of the water component leads to an increase in the density of the fluid mixture within the mixing layer compared to Case V. The cooling of water and heating up of n-decane within the mixing layer also leads to a reduction in the mixture viscosity within the layer compared to Case IV. Both, the density increase and viscosity decrease lead to an increase in the local flow Reynolds number within the mixing layer. The shear layer becomes unstable to perturbations in the flow at these higher local Re. This results in the growth of perturbations in the shear layer and the consequent waviness and rollup of the layer followed by the collapse of these structures into small-scale turbulence. Thus, the near-critical variations in mixture properties have a significant impact on the flow field and mixing dynamics in a tee mixer by influencing the local flow stability.
159
160
Chapter 6 Summary and Future Work 6.1
Summary
A robust numerical tool for simulating mixing of water and hydrocarbons under near-critical fully miscible conditions in a complicated 3-D mixer geometry was developed using the open-source CFD libraries of OpenFOAM. A consistent treatment of near-critical fluid thermodynamics and transport property variations is employed as explained in Chapter 2. The cubic Peng-Robinson equation of state [37] with van der Waals mixing rules is used to model the thermodynamic behavior of waterhydrocarbon mixtures at conditions near the critical point. The Predictive Peng Robinson (PPR78) approach [14] with a group contribution method (GCM) for modeling inter-molecular interactions is used to determine the binary interaction parameter for the water-hydrocarbon pair as a function of temperature. Viscosity and thermal conductivity calculations are performed using the generalized correlations of Chung et al. [6]. Mass diffusivities are calculated using the Tracer Liu-Silva-Macedo (TLSM) model by Liu at al. [23] [24] [46]. A 2nd order accurate finite-volume methodology is used for the numerical solution of the conservation equations as detailed in Chapter 3. An operator-splitting approach based on the PISO algorithm of Issa [13] is employed to handle the pressure-velocity coupling along with semi-implicit numerical schemes to ensure the stability of the solution. The developed numerical tool was used to investigate the mixing of supercritical 161
water and a model hydrocarbon (n-decane) in a small-scale cylindrical tee mixer (pipe ID = 2.4mm) under fully miscible conditions (T¿Upper Critical Solution Temperature of water n-decane system). First, in Chapter 4, mixing simulations were presented for a water inlet Reynolds number of 500 (the lowest Re simulated in this study), for two different inlet temperatures of the water stream of 800K and 1000K. The mass flow rate of the n-decane stream was equal to that of the water stream and its inlet temperature was 700K in both cases. This being the lowest temperature in the domain, was carefully chosen to be above the UCST of the water n-decane system (632K). It was found that the flow downstream of the mixing joint remained laminar in case of the 100K temperature difference between the streams. Most of the mixing and heat transport occurs due to the circulating action of a counter-rotating vortex pair (CVP) in the hydrocarbon jet formed due to the reorientation of the vorticity pre-existing in the hydrocarbon stream flowing through the vertical pipe. This CVP gets progressively weaker as it is advected downstream, due to vorticity diffusion and species and heat transport is dominated by molecular diffusion over small length scales in the far downstream region. Consequently, the mixing rate plateaus in the far downstream region of the tee mixer. Comparison with mixing of water and n-decane at the same inlet conditions but without variations of the physical properties with temperature suggests that for a 100K temperature difference between the streams, near-critical property variations have a negligible impact on the flow field and mixing behavior. For a 300K temperature difference between the two streams, the water-HC shear layer becomes unstable and begins to display sustained waviness near x=5D downstream of the mixing joint center. Further downstream, the shear layer rolls up and vortices are shed from it. The onset of instability in the shear layer also triggers the stretching and breakdown of the CVP in the body of the hydrocarbon jet. These smaller, more intense vortical structures cause faster advective transport of species and enthalpy over the cross-section of the tee leading to a significant enhancement in mixing. This manifests as a jump in the mixing rate at the location of onset of instability and a thickening of the mean mixing layer. However, water n-decane mix162
ing under identical inlet conditions but with constant physical properties, showed a stable shear layer with the flow reaching steady state. In the case of SCW n-decane mixing with a ∆T between the streams of 300K, the temperature of the water component changes strongly (from 1000K to 700K) within the mixing layer and this cooling of water leads to an increase in the density of the fluid mixture within the mixing layer compared to the simulated case where the properties of the components are held constant. The cooling of water and heating of n-decane within the mixing layer also leads to significant reduction in the mixture viscosity (as much as 50%) within the layer compared to the constant properties case. Both, the density increase and the viscosity decrease lead to an increase in the local flow Reynolds number within the mixing layer with the peak Re value reaching close to 1500. The shear layer becomes unstable to perturbations in the flow at these higher local Re. This results in the growth of perturbations in the shear layer and the consequent waviness and rollup of the layer followed by the collapse of these structures into small-scale turbulence. In Chapter 5, simulations of mixing of supercritical water and n-decane, were presented for four values of the water inlet Reynolds number in the range of 500 to 800. The mass flow rate of the n-decane stream was equal to that of the water stream. The inlet temperatures of water and n-decane were 800K and 700K respectively. It was found that the flow downstream of the mixing joint remained laminar for Rew,in = 500 and Rew,in = 600. For Rew,in = 600, the flow and mixing behavior were similar to the case of Rew,in = 500 with a slight enhancement in mixing rate due to the increase in the strength of the CVP with Re. When the Reynolds number at the water inlet is increased to 700 and the n-decane inlet Re also increased correspondingly to maintain the same mass flow rate, the shear layer between the water and n-decane streams is found to destabilize near x = 6D downstream of the mixing joint followed by subsequent shear layer roll up and vortex shedding. The onset of instability in the shear layer also triggers the stretching and breakdown of the HC jet CVP into smaller, stronger streamwise vortices which significantly enhance mixing by stretching the material surface between the two streams. This manifests as a jump in the mixing rate at the location of onset of instability and a thickening of the mean mixing layer 163
similar to the case of large ∆T in Chapter 4. At Rew,in = 800 the unsteady small scale flow structures in the mixing layer and the consequent flow field fluctuations due to them are much stronger. The stretching and breakdown of the CVP in this case, is accompanied by stronger streamwise vorticity enhancement resulting in much faster mixing compared to the case of Rew,in = 700. However, water n-decane mixing under identical inlet conditions to the Rew,in = 700 case but with constant physical properties, showed a stable shear layer with the flow reaching steady state. In SCW n-decane mixing with Rew,in = 700, the temperature of the water component changes from 800K to 700K within the mixing layer and this cooling of water leads to an increase in the density of the fluid mixture within the mixing layer compared to the constant properties case. The cooling of water and heating of n-decane within the mixing layer also leads to a reduction in the mixture viscosity within the layer compared to the constant properties case. Both, the density increase and the viscosity decrease lead to an increase in the local flow Reynolds number within the mixing layer. The shear layer becomes unstable to perturbations in the flow at these higher local Re. This results in the growth of perturbations in the shear layer and the consequent waviness and rollup of the layer followed by the collapse of these structures into small-scale turbulence. Thus, the near-critical variations in mixture properties have a significant impact on the flow field and mixing dynamics in a tee mixer by influencing the local flow stability.
6.2
Future Work
The work presented in this thesis is the first important step in an endeavor to simulate the reactive mixing of water and crude oil in complicated 3-D reactor geometries. Though the present work deals with some aspects of this complicated problem, it still has a number of shortcomings which will be the focus of future work in this project. The numerical tool developed in this study does not have the ability to track the phase interface between immiscible phases. In a realistic SCWDS and upgrading process, the local environment in the reactor may not allow complete miscibility of 164
water and hydrocarbon phases. As such, the ability to track the phase interface and in general, multiple phase interfaces is essential. This can be done using sophisticated interface tracking methods like the Level-Set method or the Volume of Fluid method. In addition, phase equilibrium calculations need to be performed dynamically to determine the equilibrium compositions of the phases near at the interface. Also, a consistent treatment of the heat and species transport across the interface is required to simulate transport and mixing in this multiphase scenario. This development of the numerical tool is crucial in our efforts to realistically simulate the mixing of water and hydrocarbons in a range of reactor conditions. Computational cost constraints currently limit our ability to simulate mixing at Reynolds numbers higher than 1000. Since it has been found desirable to operate the reactor under turbulent conditions in this study, we need to develop a framework to simulate high Reynolds number flows with reasonable computational effort. For this, turbulence modeling using a Large Eddy Simulation methodology must be implemented. This will involve the development of turbulent mixing models which incorporate the key physics of near-critical fluids. Inclusion of chemical reactions in the numerical tool is also an essential future step. Simplified reaction mechanisms for thermal cracking reactions of hydrocarbons and sulfur compounds as well as cyclization and condensation reactions leading to coke formation need to be considered in order to investigate the coupling between transport and chemical reaction. This will help in obtaining realistic predictions of conversion rates and product distributions which is the ultimate goal of our collaborative research effort.
165
166
Bibliography [1] W. Abdoul, E. Rauzy, and A. P´eneloux. Group-contribution equation of state for correlating and predicting thermodynamic properties of weakly polar and non-associating mixtures: Binary and multicomponent systems. Fluid Phase Equilibria, 68:47–102, 1991. [2] U. Andersson, J. Westin, and J. Eriksson. Thermal Mixing in a T-junction. Technical report, Vattenfall Research and Development AB, 2006. [3] Manson Benedict, George B. Webb, and Louis C. Rubin. An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures I: Methane, Ethane, Propane and n-Butane. The Journal of Chemical Physics, 8:334, 1940. [4] Sydney Chapman and Thomas George Cowling. The mathematical theory of nonuniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases. Cambridge university press, 1991. [5] Zhen-Min Cheng, Yong Ding, Li-Qun Zhao, Pei-Qing Yuan, and Wei-Kang Yuan. Effects of supercritical water in vacuum residue upgrading. Energy & Fuels, 23(6):3178–3183, 2009. [6] Ting Horng Chung, Mohammad Ajlan, Lloyd L. Lee, and Kenneth E. Starling. Generalized multiparameter correlation for nonpolar and polar fluid transport properties. Industrial & engineering chemistry research, 27(4):671–679, 1988. [7] A. Correia, A. Afonso, A. Cavadas, M.S.N. Oliveira, M.A. Alves, and F.T. Pinho. Numerical investigation of laminar flow in symmetric and asymmetric complex Tchannels for micro-mixing. In III Conferˆencia Nacional em Mecˆanica de Fluidos, Termodinˆamica e Energia (MEFTE-BRAGANC ¸ A 09), 2009. [8] Sadegh Dabiri, Guang Wu, Michael T. Timko, and Ahmed F. Ghoniem. Mixing of single-component hydrocarbon droplets and water at supercritical or nearcritical conditions. The Journal of Supercritical Fluids, 2012. [9] Simon Dreher, Norbert Kockmann, and Peter Woias. Characterization of laminar transient flow regimes and mixing in T-shaped micromixers. Heat Transfer Engineering, 30(1-2):91–100, 2009. 167
[10] Arnaud Erriguible, St´ephane Vincent, and Pascale Subra-Paternault. Numerical investigations of liquid jet breakup in pressurized carbon dioxide: Conditions of two-phase flow in Supercritical Antisolvent Process. The Journal of Supercritical Fluids, 63:16–24, 2012. [11] T.H. Frank, C. Lifante, H-M. Prasser, and F. Menter. Simulation of turbulent and thermal mixing in T-junctions using URANS and scale-resolving turbulence models in ANSYS CFX. Nuclear Engineering and Design, 240(9):2313–2328, 2010. [12] Marko Hoffmann, Michael Schl¨ uter, and Norbert R¨abiger. Experimental investigation of liquid–liquid mixing in T-shaped micro-mixers using µ-LIF and µ-PIV. Chemical engineering science, 61(9):2968–2976, 2006. [13] Raad I. Issa. Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational physics, 62(1):40–65, 1986. [14] Jean-No¨el Jaubert, Romain Privat, and Fabrice Mutelet. Predicting the phase equilibria of synthetic petroleum fluids with the PPR78 approach. AIChE Journal, 56(12):3225–3235, 2010. [15] S.T. Jayaraju, E.M.J. Komen, and E. Baglietto. Suitability of wall-functions in Large-Eddy Simulation for thermal fatigue in a T-junction. Nuclear Engineering and Design, 240(10):2544–2554, 2010. [16] H.V. Kehiaian, K. Sosnkowska-Kehiaian, and R. Hryniewicz. Enthalpy of mixing of ethers with hydrocarbons at 25o C and its analysis in terms of molecular surface interactions. J. Chim. phys, 68:922–934, 1971. [17] Taehoon Kim, Yongmo Kim, and Seong-Ku Kim. Numerical study of cryogenic liquid nitrogen jets at supercritical pressures. The Journal of Supercritical Fluids, 56(2):152–163, 2011. [18] Taehoon Kim, Yongmo Kim, and Seong-Ku Kim. Real-fluid flamelet modeling for gaseous hydrogen/cryogenic liquid oxygen jet flames at supercritical pressure. The Journal of Supercritical Fluids, 58(2):254–262, 2011. [19] A.K. Kuczaj and E.M.J. Komen. An assessment of Large-Eddy Simulation toward thermal fatigue prediction. Nuclear technology, 170(1):2–15, 2010. [20] A.K. Kuczaj, E.M.J. Komen, and M.S. Loginov. Large-Eddy Simulation study of turbulent mixing in a T-junction. Nuclear Engineering and Design, 240(9):2116– 2122, 2010. [21] Michael Charles Kutney. Thermodynamic and transport property modeling in super critical water. PhD thesis, Massachusetts Institute of Technology, 2005. 168
[22] Byung Ik Lee and Michael G. Kesler. A generalized thermodynamic correlation based on three-parameter corresponding states. AIChE Journal, 21(3):510–527, 1975. [23] Hongqin Liu, Carlos M. Silva, and Eug´enia A. Macedo. New equations for tracer diffusion coefficients of solutes in supercritical and liquid solvents based on the Lennard-Jones fluid model. Industrial & engineering chemistry research, 36(1):246–252, 1997. [24] Hongqin Liu, Carlos M. Silva, and Eug´enia A. Macedo. Unified approach to the self-diffusion coefficients of dense fluids over wide ranges of temperature and pressure, hard-sphere, square-well, lennard-jones and real substances. Chemical engineering science, 53(13):2403–2422, 1998. [25] J.J. Martin and T.G. Stanford. Development of high precision equations of state for wide ranges of density utilizing a minimum of input information: example argon. AIChE Symp. Ser., 70:1–13, 1974. [26] Youssef M. Marzouk and Ahmed F. Ghoniem. Vorticity structure and evolution in a transverse jet. Journal of Fluid Mechanics, 575:267, 2007. [27] Paul M. Mathias, Tarik Naheiri, and Edwin M. Oh. A density correction for the Peng Robinson equation of state. Fluid Phase Equilibria, 47(1):77–87, 1989. [28] Elia Merzari, Azizul Khakim, Hisashi Ninokata, and Emilio Baglietto. Unsteady Reynolds-averaged Navier-Stokes: toward accurate prediction of turbulent mixing phenomena. International Journal of Process Systems Engineering, 1(1):100– 123, 2009. [29] Richard S. Miller, Kenneth G. Harstad, and Josette Bellan. Direct numerical simulations of supercritical fluid mixing layers applied to heptane-nitrogen. Journal of Fluid Mechanics, 436(6):1–39, 2001. [30] R.D. Moser and M.M. Rogers. Mixing transition and the cascade to small scales in a plane mixing layer. Physics of Fluids A: Fluid Dynamics, 3:1128, 1991. [31] C. Narayanan, C. Frouzakis, K. Boulouchos, K. Pr´ıkopsk` y, B. Wellig, and P. Rudolf von Rohr. Numerical modelling of a supercritical water oxidation reactor containing a hydrothermal flame. The Journal of Supercritical Fluids, 46(2):149–155, 2008. [32] National Institute of Standards and Technology. NIST Chemistry WebBook, 2013. [33] Ylva Odemark, Torbj¨orn M. Green, Kristian Angele, Johan Westin, Farid Alavyoon, and Staffan Lundstr¨om. High-cycle thermal fatigue in mixing tees: New Large-Eddy Simulations validated against new data obtained by PIV in the Vattenfall experiment. ASME, 2009. 169
[34] Nora Okong’o, Kenneth Harstad, and Josette Bellan. Direct numerical simulations of O/H temporal mixing layers under supercritical conditions. AIAA journal, 40(5):914–926, 2002. [35] OpenCFD Ltd. The open source CFD toolbox, 2013. [36] Tae Seon Park. LES and RANS simulations of cryogenic liquid nitrogen jets. The Journal of Supercritical Fluids, 2012. [37] Peng, Ding-Yu and Robinson, Donald B. A new two-constant equation of state. Industrial & Engineering Chemistry Fundamentals, 15(1):59–64, 1976. [38] British Petroleum. BP Statistical Review of World Energy, 2012. [39] Junwei Qian, Romain Privat, and Jean-Noeal Jaubert. Predicting the phase equilibria, critical phenomena and mixing enthalpies of binary aqueous systems containing alkanes, cycloalkanes, aromatics, alkenes and gases (N2, CO2, H2S, H2) with the PPR78 equation of state. Industrial & Engineering Chemistry Research, 2013. [40] Otto Redlich and J.N.S. Kwong. On the thermodynamics of solutions V: An Equation of State, Fugacities of gaseous solutions. Chemical Reviews, 44(1):233– 244, 1949. [41] C.M. Rhie and W.L. Chow. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal, 21(11):1525–1532, 1983. [42] Fabrice Schlegel and Ahmed F. Ghoniem. Simulation of a high Reynolds number reactive transverse jet and the formation of a triple flame. Combustion and Flame, 2013. [43] Fabrice Schlegel, Daehyun Wee, Youssef M. Marzouk, and Ahmed F. Ghoniem. Contributions of the wall boundary layer to the formation of the counter-rotating vortex pair in transverse jets. Journal of Fluid Mechanics, 676:461, 2011. [44] T. Schmitt, J. Rodriguez, I.A. Leyva, and S. Candel. Experiments and numerical simulation of mixing under supercritical conditions. Physics of Fluids, 24:055104, 2012. [45] Jos´e Sierra-Pallares, Mar´ıa Teresa Parra-Santos, Juan Garc´ıa-Serna, Francisco Castro, and Mar´ıa Jos´e Cocero. Numerical modelling of hydrothermal flames: Micromixing effects over turbulent reaction rates. The Journal of Supercritical Fluids, 50(2):146–154, 2009. [46] Carlos M. Silva, Hongqin Liu, and Eug´enia A. Macedo. Models for self-diffusion coefficients of dense fluids, including hydrogen-bonding substances. Chemical engineering science, 53(13):2423–2429, 1998. 170
[47] Giorgio Soave. Equilibrium constants from a modified redlich-kwong equation of state. Chemical Engineering Science, 27(6):1197–1203, 1972. [48] Marios C. Soteriou and Ahmed F. Ghoniem. Effects of the free-stream density ratio on free and forced spatially developing shear layers. Physics of Fluids, 7:2036, 1995. [49] K.E. Starling and M.S. Han. Thermo data refined for LPG-14 mixtures. Hydrocarbon processing, 51(5), 1972. [50] US EPA. Code of Federal Regulations, Schedule 40, 80H, 2013. [51] US EPA. Code of Federal Regulations, Schedule 40, 80I, 2013. [52] J.D. van der Waals. Over de continuiteit van den gas-en vloeistoftoestand. PhD thesis, Leiden, Holland, 1873. [53] Javier Vilc´aez, Masaru Watanabe, Noriaki Watanabe, Atsushi Kishita, and Tadafumi Adschiri. Hydrothermal extractive upgrading of bitumen without coke formation. Fuel, 2012. [54] Johan Westin, C. Mannetje, F. Alavyoon, P. Veber, L. Andersson, U. Andersson, J. Eriksson, M. Henriksson, and C. Andersson. High-cycle thermal fatigue in mixing tees; Large-Eddy Simulations compared to a new validation experiment. ICONE16, Orlando, 2008. [55] David Wood. Consequences of a heavier and sourer barrel. Petroleum Review, pages 30–32, 2007. [56] Guang Wu, Sadegh Dabiri, Michael T. Timko, and Ahmed F. Ghoniem. Fractionation of multi-component hydrocarbon droplets in water at supercritical or near-critical conditions. The Journal of Supercritical Fluids, 2012. [57] Carl L. Yaws. Yaws handbook of thermodynamic properties for hydrocarbons and chemicals. Knovel, 2009. [58] Nan Zong, Hua Meng, Shih-Yang Hsieh, and Vigor Yang. A numerical study of cryogenic fluid injection and mixing under supercritical conditions. Physics of fluids, 16:4248, 2004.
171
