NUMERICAL STUDY OF VARIABLE DENSITY FLOW INTERACTION WITH A SHOCK WAVE By Yifeng Tian A DISSERTATION Michigan State University in partial fulfillment of the requirements Submitted to for the degree of Mechanical Engineering – Doctor of Philosophy 2019 ABSTRACT NUMERICAL STUDY OF VARIABLE DENSITY FLOW INTERACTION WITH A SHOCK WAVE By Yifeng Tian Fundamental understanding and modeling of multi-fluid miscible Shock-Turbulence Interaction (STI) and the corresponding post-shock turbulence are critically important to many different ap- plications, such as supersonic combustion, nuclear fusion, and astrophysics. This thesis presents a comprehensive study of the multi-fluid Shock-Turbulence flow using accurate flow-resolving, shock- capturing, and shock-resolving simulations. The objective is to develop a better understanding of underlying mechanisms of the variable density fluid effects on shock-turbulence interactions, post- shock turbulence evolution and mixing in high speed flows. Theoretical and numerical analyses of data confirm that all turbulence scales as well as the STI are well captured by the computational method, which is based on a high order hybrid monotonicity preserving-compact finite difference scheme. Linear Interaction Approximation (LIA) convergence tests are conducted to show that shock-capturing numerical simulations exhibit similar converging trend to LIA predictions as more demanding shock-resolving Direct Numerical Simulations (DNS) method. The effects of density variations on STI are studied first by comparing the “Eulerian” results obtained from both Eulerian (grid) and Lagrangian (particle) data for a multi-fluid mixture with the corresponding single-fluid results. The comparison shows that the turbulence amplification by the normal shock wave is much higher and the reduction in turbulence length scales is more significant when strong density variations are present. Turbulent mixing enhancement by the shock is also increased and stronger mixing asymmetry in the post-shock region is observed when there are significant density variations. The dominating mechanisms behind STI influence on post-shock turbulence and mixing are identified by analyzing the transport equations for the Reynolds stresses, vorticity, normalized mass flux, and density specific volume covariance. Statistical analyses of the velocity gradient tensor (VGT) show that the density variations also significantly change the turbulence structure and flow topology. Compared to the single-fluid case, the correlation between rotation and strain is found to be weaker in the multi-fluid case, which is shown to be the result of complex role density plays when the flow passes through the shock. Furthermore, a stronger symmetrization of the joint PDF of second and third invariants of the anisotropic velocity gradient tensor, as well as the PDF of the vortex stretching contribution to the enstrophy equation, are observed in the multi-fluid case. Lagrangian dynamics of the VGT and its invariants are studied by considering particle residence times in different flow regions and the conditional mean rate of change vectors. The pressure Hessian contributions to the VGT invariants transport equations are shown to be strongly affected by the shock wave and local density, making them critically important to the flow dynamics and turbulence structure. Lagrangian statistics of non-inertial particles in post-shock turbulence show that the single-particle dispersion rate and anisotropy can be correlated to the development of post-shock Reynolds stress. Lagrangian integral time scales of fluid particles with different densities are obtained from velocity autocorrelation functions. Particle pair dispersion generally follows the temporal scaling for isotropic turbulence. In this thesis, the propagation of shock waves in non-uniform density media is also studied. The 1D numerical results are shown to agree well with theoretical solutions obtained from classical Chisnell-Whitham theory for weak shocks and linearly varying density fields. For more significantly varying density profiles, the numerical results deviate from the theoretical solutions and exhibit additional long-wavelength oscillations, which are shown to be related to the re-reflected waves. A simplified model (named CWRW) that includes the effects of re-reflected waves is proposed to compensate for the CW model definitions. The results obtained by the proposed CWRW model show significant improvement over the classical CW model. The three-dimensional (3D) effects concerning the shock propagation in 3D non-uniform density media are studied by considering the wavenumber ratio of the streamwise and spanwise length scales in the density field. The asymptotic limits for wavelength ratio approaching zero and infinity are investigated and used to provide a bound for shock propagation in 3D non-uniform media. Copyright by YIFENG TIAN 2019 ACKNOWLEDGEMENTS First and foremost, I would like to express my deepest gratitude to my advisor Dr. Farhad Jaberi and my co-advisor Dr. Daniel Livescu for their invaluable guidance and support through my Ph.D. study. This work would not be possible without their patience and consistent encouragement. It has been a truly wonderful experience to work with Dr. Jaberi and Dr. Livescu. I would also like to thank my Ph.D. committee members Dr. Ahmed Naguib and Dr. Mohsen Zayernouri for their constructive suggestions and comments on my research. The courses they provided have been very helpful and prepared me with fundamental knowledge of compressible flow and numerical methods. Additionally, I would like to acknowledge Dr. Zhaorui Li for his guidance and suggestions on the code development. Financial support for this work was provided by Los Alamos National Laboratory, under Grant No. 319838. Computational resources were provided by the High Performance Computing Center at Michigan State University and Texas Advanced Computing Center (TACC) at The University of Texas at Austin. I would also like to thank my former colleague Dr. AbdoulAhad Validi and Dr. Husam Abdulrahman who helped me start my research at the CFDLab. They have given me many helpful pieces of advice on research over the past years. I would like to thank all of my beloved friends for their emotional support and involvement in various academic and extracurricular activities with me. Last but not least, I would like to thank my express my deepest appreciation and gratitude to my parents and my girlfriend. This dissertation would not have been possible without their warm love, continued patience, and endless support. v TABLE OF CONTENTS LIST OF TABLES . LIST OF FIGURES . x KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . xix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . CHAPTER 1 GENERAL BEHAVIOR OF VARIABLE DENSITY TURBULENCE IN- . . . . . . . . Introduction . TERACTION WITH A NORMAL SHOCK WAVE . . . . . . . . . . . . . 1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Governing Equations and Numerical Methodology . . . . . . . . . . . . . . . . . . 1.3 Simulations Setup and Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Results and Discussions . 1 1 4 7 9 1.4.1 Accuracy of Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4.2 Convergence to LIA for Single-fluid Case . . . . . . . . . . . . . . . . . . 14 1.4.3 Effects of Density Variations on STI . . . . . . . . . . . . . . . . . . . . . 19 1.4.4 Reynolds Stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.5 Vorticity Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.4.6 Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.4.7 Normalized Mass Flux and Density Specific Volume Covariance . . . . . . 44 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 1.5 Conclusions . . . . . . . . . . . . . CHAPTER 2 DENSITY EFFECTS ON THE POST-SHOCK TURBULENCE STRUC- . . . . . Introduction . 2.1 . 2.2 Numerics and Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . TURE AND LAGRANGIAN FLOW DYNAMICS . . . . . . . . . . . . . . 53 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 . 55 2.2.1 Governing Equations and Numerical approach . . . . . . . . . . . . . . . . 55 2.2.2 Numerical Setup of Eulerian Simulation . . . . . . . . . . . . . . . . . . . 56 2.2.3 Numerical Detials of the Lagrangian Analyses . . . . . . . . . . . . . . . . 56 2.2.4 Grid and Statistical Convergence . . . . . . . . . . . . . . . . . . . . . . . 58 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.3 Results and Discussions . . 2.3.1 Density Effects on the Post-shock Turbulence and Its Evolution Away . From the Shock . 2.3.1.1 2.3.1.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Turbulence State Immediately After the Shock . . . . . . . . . . 62 Evolution of Turbulence State Away from the Shock Wave . . . . 67 2.3.2 Topological Analysis of the Post-shock Turbulence . . . . . . . . . . . . . 74 2.3.3 Lagrangian Dynamics of VGT . . . . . . . . . . . . . . . . . . . . . . . . 81 2.3.4 Lagrangian Particle Dispersion . . . . . . . . . . . . . . . . . . . . . . . . 94 Single Particle Statistics . . . . . . . . . . . . . . . . . . . . . . 94 Particle Pair Statistics . . . . . . . . . . . . . . . . . . . . . . . 98 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 . . CHAPTER 3 EFFECTS OF DENSITY VARIATIONS ON A NORMAL SHOCK WAVE . 106 2.3.4.1 2.3.4.2 2.4 Conclusions . . . . vi . . Introduction . . 3.1 3.2 Numerical Methodology . . . . . . 3.4.1 3.4.2 3.2.1 3.2.2 3.2.3 3.3 Theoretical Analysis . 3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 . 109 Shock-capturing Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 109 Shock-resolving Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation Setup . . 110 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 . . 111 . 118 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Shock Wave Propagation in Media with Linearly Varying Density . . . . . 120 Shock Wave Propagation in Media with Fluctuating Density . . . . . . . . 122 3.4.2.1 3.3.1 Linear Interaction Approximation . . . . . . . . . . . . . . . . . . . . . 3.3.2 Chisnell-Whitham Theory . . . . . . . . . . . . . . . . . . . . . . . . . The Similarity of Density Fluctuations from Entropy and Com- position Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . 122 3.4.2.2 Shock Propagation in Sinusoidal Density Media . . . . . . . . . 124 3.4.2.3 CW Solution with Correction for Re-reflected Waves . . . . . . . 128 3.4.2.4 Shock Propagation in Random Density Media . . . . . . . . . . . 132 Shock Wave Propagation in Three-dimensional Non-uniform Density Field 133 3.4.3.1 . 135 . 139 3.4.3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 3D Sinusoidal Density Field . . . . . . . . . . . . . . . . . . . 3D Random Density Field . . . . . . . . . . . . . . . . . . . . . . 3.4.3 . . . . . . . . . . . . BIBLIOGRAPHY . . . vii LIST OF TABLES Table 1.1: Details of the simulations used in the grid convergence tests. . . . . . . . . . . . 10 Table 1.2: Cases considered for LIA convergence tests. . . . . . . . . . . . . . . . . . . . . 17 viii LIST OF FIGURES Figure 1.1: Scalar structure in multi-fluid turbulence interacting with a Mach 2 shock identified by the iso-surface of heavy fluid mole fraction and colored with instantaneous density fluctuations. The black plane located in the middle of the domain represents the instantaneous shock surface. The iso-surfaces correspond to mole fraction values of (a) 0.05 (b) 0.95 (c) 0.25 and (d) 0.75. . . 6 Figure 1.2: Instantaneous contours of 3D pressure fluctuations and shock surface resulting from isotropic turbulence interacting with a Mach 2 shock. (a) Iso-surface of heavy fluid mole fraction (φ = 0.25) colored by instantaneous pressure fluctuations. Instantaneous shock surface colored by the ratio of local pressure jump across the shock for (b) multi-fluid A (see definition in section 1.4.3) and (c) single-fluid cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Figure 1.3: Results of multi-fluid grid convergence tests at Reλ = 45 and Mt = 0.09. (a) Turbulent kinetic energy, (b) Kolmogorov length scale η, (c) streamwise Taylor micro scale λ1 and (d) transverse vorticity variance ω(cid:48) 2ω(cid:48) 2. The region . . . . . . . . . . . of unsteady wrinkled shock movement is marked in grey. . 11 Figure 1.4: Results of multi-fluid grid convergence tests at Reλ = 45 and Mt = 0.09. (a) mass fraction variance φ(cid:48)φ(cid:48) and (b) Batchelor scale λB. . . . . . . . . . . . . . 11 Figure 1.5: Relation between numerical shock width, δn, and grid size in the shock region. 12 Figure 1.6: Comparison of Reynolds stress obtained from single-fluid simulations with LIA results. Mt ≈ 0.24 except for Reλ = 45. (a) Streamwise Reynolds stress and (b) transverse Reynolds stress. . . . . . . . . . . . . . . . . . . . . . . . . 18 Figure 1.7: Comparison of Reynolds stress obtained from multi-fluid, single-fluid simu- lation with LIA results. Reλ ≈ 30. (a) Streamwise Reynolds stress and (b) transverse Reynolds stress. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Figure 1.8: 2D contours of instantaneous scalar and density fields for (a,b) multi-fluid A simulation and (c,d) single-fluid simulation. The 2D contours are taken at the same time step, in planes x-z. . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 1.9: Plots of (a) turbulent kinetic energy, (b) Kolmogorov length scale, (c) Taylor micro scale and (d) transverse vorticity variance for multi-fluid A (red) and single-fluid (blue) simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . 22 ix Figure 1.10: Plots of normalized (a) scalar variance, (b) Batchelor scale for multi-fluid A (red) and single-fluid (blue) simulations. . . . . . . . . . . . . . . . . . . . . . 22 Figure 1.11: Conditional expectation of scalar as a function of density at various stream- wise location for : (a) multi-fluid A and (b) single-fluid cases. The conditional mean is calculated using variables in a spanwise plane at certain streamwise locations and then averaged over time. Peak TKE positions are different between multi-fluid(k0x ≈ 2.0) and single-fluid(k0x ≈ π) cases . . . . . . . . . 23 Figure 1.12: Conditional expectation of pressure as a function of density at various stream- wise locations for : (a) multi-fluid A close to shock, (b) multi-fluid A at far field and (c) single-fluid cases. . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Figure 1.13: Conditional expectation of temperature as a function of density at various streamwise location for : (a) multi-fluid A case and (b) single-fluid case. . . . . 25 Figure 1.14: Conditional expectation of TKE as a function of density at various streamwise location for : (a) multi-fluid A and (b) single-fluid case. . . . . . . . . . . . . . 26 Figure 1.15: Conditional expectation of enstrophy as a function of density at various streamwise location for : (a) multi-fluid A and (b) single-fluid case. . . . . . . . 27 Figure 1.16: Contributions of different terms in the R11 transport equations for multi-fluid A simulation. (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. . . . . . . . . . . . 30 Figure 1.17: Contributions of different terms in the R11 transport equations for multi- fluid A simulation averaged over 5 different realizations (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 1.18: Contributions of different terms in the R22 transport equations for multi- fluid A simulation averaged over 5 different realizations. (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 1.19: Normalized contributions of different terms in the R11 transport equations averaged over 5 different realizations for (a) multi-fluid A and (b) single-fluid simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 . . . . Figure 1.20: Normalized contributions of different terms in the R22 transport equations averaged over 5 different realizations for (a) multi-fluid A and (b) single-fluid simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 . . . . x Figure 1.21: Two dimensional spectra of transverse turbulent kinetic energy. The compar- ison are made at right before shock wave and k0x ≈ 4.0. . . . . . . . . . . . . 34 Figure 1.22: Streamwise and transverse vorticity variance for multi-fluid and single-fluid case 35 Figure 1.23: Contributions of different terms in the transport equation for (a) ω(cid:48) 1ω(cid:48) 1 and . . . . . . . . . . . . . . . . . . . . . . . 36 (b) ω(cid:48) 2ω(cid:48) 2 for multi-fluid simulation. Figure 1.24: Contributions of different terms in the transport equation for (a) ω(cid:48) 1ω(cid:48) 1 and . . . . . . . . . . . . . . . . . . . . . . (b) ω(cid:48) 2ω(cid:48) 2 for single-fluid simulation. . 36 Figure 1.25: Vortex-stretching term in the enstrophy equation, normalized by ω2ω2 and turbulence time scale TKE/ . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Figure 1.27: Contributions of different terms in the scalar variance(cid:103)φ Figure 1.28: Contributions of different terms in the(cid:103)φ Figure 1.26: Plots of scalar Taylor micro scale for multi-fluid (red) and single-fluid (blue). (cid:48)(cid:48)2 transport equations (a) Plots of all the contributing terms in the for multi-fluid simulation. transport equation and (b) balance of LHS and RHS of the transport equation. (cid:48)(cid:48)2 transport equations for (a) multi- . . . . . . . . . . . . . . . . . . . . . . . 42 fluid and (b) single-fluid simulations. . 40 . 41 Figure 1.29: Comparison of scalar dissipation and advection terms obtained from multi- fluid and single-fluid simulations. . . . . . . . . . . . . . . . . . . . . . . . . 43 Figure 1.30: 2D spectra of scalar. The comparison are made at various locations: (a) right before shock wave and k0x ≈ 4.0 and (b) right before shock wave and peak TKE position. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 . . . . Figure 1.31: PDFs of multi-fluid A case: (a) scalar (mole fraction) and (b) density fluc- tuation at various locations. The PDFs are calculated at specified y-z planes and then averaged over time. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Figure 1.32: PDFs of multi-fluid B case: (a) scalar (mass fraction) and (b) density fluc- tuation at various locations. The PDFs are calculated at specified y-z planes and then averaged over time. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Figure 1.33: PDFs of single-fluid case: (a) passive scalar and (b) density fluctuation at various locations. The PDFs are calculated at specified y-z planes and then averaged over time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 . xi Figure 1.34: Plots of (a) normalized mass flux a1 and (b) density specific volume covari- ance − (cid:104)v ρ(cid:105) for multi-fluid and single-fluid simulations. . . . . . . . . . . . . Figure 1.35: Contributions of different terms in the normalized mass flux a1 transport equations for multi-fluid simulation. . . . . . . . . . . . . . . . . . . . . . . . 48 . 48 Figure 1.36: Contributions of different terms in the density specific volume covariance b transport equations for multi-fluid simulation. . . . . . . . . . . . . . . . . . . 49 Figure 2.1: Instantaneous contours of vorticity and shock surface in isotropic turbulence interacting with a Mach 2 shock. (a) Vortex structure, identified by the Q criteria (Q = 2 (cid:104)Qw(cid:105)) and colored by the mole fraction of the heavy fluid. Fluid particles are initialized as a sheet that spans over the homogeneous directions at a given post-shock streamwise position and allowed to develop with the flow. (b) Visualized particle sheet, convected and distorted by the post-shock turbulence. The instantaneous shock surface is colored by the shock intensity across the shock for (c) single-fluid and (d) multi-fluid cases. . . 56 Figure 2.2: Results of multi-fluid grid convergence tests at Reλ = 45 and Mt = 0.09. Streamwise development of (a) velocity dissipation rate ε and (b) mass frac- tion dissipation rate εφ is shown. The region of unsteady shock move- ment is marked in grey. Grid numbers for Grid 1 to 5 are 256×256×1024, 384×384×1024, 384×384×1536, 512×512×1536, 512×512×2048. . . . . . . . 60 Figure 2.3: The statistical convergence for (a) (DQ/Dt)/(cid:104)Qw(cid:105)3/2 and (DR/Dt)/(cid:104)Qw(cid:105)2 and (b) their standard deviations conditioned at point (3.0,3.0) in the (Q, R) phase plane for multi-fluid case. The number of bins are 30 × 30 (solid), 40 × 40 (dashed) and 60 × 60 (dotted). . . . . . . . . . . . . . . . . . . . . . . 61 Figure 2.4: The statistical convergence for (a) (DQ/Dt)/(cid:104)Qw(cid:105)3/2 and (DR/Dt)/(cid:104)Qw(cid:105)2 and (b) their standard deviations conditioned at point (3.0,3.0) in the (Q, R) phase plane for single-fluid case. The number of bins are 30 × 30 (solid), 40 × 40 (dashed) and 60 × 60 (dotted). . . . . . . . . . . . . . . . . . . . . . . 61 Figure 2.5: Comparison of the PDFs of the normalized post-shock velocity derivatives with a Gaussian distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Figure 2.6: PDF of the strain-enstrophy angle Ψ in radian for post-shock turbulence. . . . . 64 Figure 2.7: Conditional expectation of the magnitude of (a) rotation tensor and (b) strain rate tensor as a function of density before and after the shock. . . . . . . . . . . 64 Figure 2.8: (a) Demonstration of the baroclinic torque generation across the shock wave and (b) PDF of the orientation between the vorticity vector and density gra- dient in y-z direction immediately after the shock wave. . . . . . . . . . . . . . 66 xii Figure 2.9: Vortex structures captured using the Q-criterion, colored by density, for multi- fluid (a) pre-shock turbulence and (b) post-shock turbulence. . . . . . . . . . . 67 Figure 2.10: Development of (a) TKE and dissipation, (b) pressure variance, (c) vortex stretching, and (d) anisotropy (b11) of Reynolds stress and vorticity. . . . . . . . 69 Figure 2.11: Development of (a) skewness, and (b) flatness, of the streamwise and trans- verse components of velocity derivatives. . . . . . . . . . . . . . . . . . . . . . 71 Figure 2.12: PDF of the density gradient at different streamwise locations for: (a) single- fluid case and (b) multi-fluid case. . . . . . . . . . . . . . . . . . . . . . . . . 71 Figure 2.13: PDFs of the cosine angle between eigenvectors of the strain rate tensor and streamwise (x) axis for regions with (a) dρ/dx > 0 and (b) dρ/dx < 0 and for multi-fluid (solid lines) and single-fluid (dashed lines) cases. . . . . . . . . . 73 Figure 2.14: Iso-contour lines of joint PDFs of normalized second and third invariants of the anisotropic part of the velocity gradient tensor, (Q∗,R∗), for (a) pre- shock, (b) single-fluid post-shock turbulence, and (c) multi-fluid post-shock turbulence. The lateral lines denote the locus of zero discriminant. . . . . . . . 76 Figure 2.15: Iso-contour lines of post-shock (k0x ≈ 0.44) joint PDF of second and third invariants of the anisotropic part of the velocity gradient tensor, (Q∗,R∗), (a) regions with high density values, in regions with different densities. ρ > (ρ + 90%ρ(cid:48) rms), (b) regions with density around the post-shock mean value, and (c) regions with low density values, ρ < (ρ − 90%ρ(cid:48) rms). . . . . . . 77 Figure 2.16: Color illustration of the flow topology for the multi-fluid STI. The flow topology is represented by the quadrants (denoted as Qi) of the joint PDF of (Q∗,R∗). (a) 2D color contours in the x-z plane at y = 3.14 (half y-domain). The shock wave is located in the middle of the domain at k0x ≈ 0.0. The ratio of the fluid volume in different quadrants in the pre- shock region is Q1 : Q2 : Q3 : Q4=26.7%:38.7%:7.8%:26.8%. The 2D color contours in the homogeneous y-z plane at streamwise locations of: (b) k0x ≈ 0.2 (28.7%:34.4%:14.3%:22.6%), (c) k0x ≈ 2.0, peak TKE location (26.7%:36.9%:11.2%:25.2%) (d) k0x ≈ 4.0 (26.3%:37.9%:9.3%:26.2%). . . . . 78 Figure 2.17: Iso-contour lines of joint PDF of (-Q∗ s , Σ∗) for (a) isotropic box turbulence and (b,c) single-fluid and multi-fluid turbulence at post-shock position of k0x ≈ 0.44. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 . . . . . Figure 2.18: Iso-contour lines of joint PDF of (-Q∗ s , Σ∗) for different quadrants right after the shock wave. (a) Q2, (b) Q1 ,(c) Q3 and (d) Q4. . . . . . . . . . . . . . . . . 81 xiii Figure 2.19: Percentage of fluid particles that stay in each quadrant following particles initialized uniformly in (a) isotropic turbulence and (b,c) single-fluid and multi-fluid turbulence at post-shock position of k0x = 0.44. . . . . . . . . . . . 83 Figure 2.20: Visualization of turbulence structure using iso-surfaces of Q∗ colored by R∗ in the multi-fluid post-shock turbulence. (a) Q1 and Q2 (vorticity-dominated re- gion), and (b) Q3 and Q4 (strain-dominated region). The black lines represent the instantaneous streamlines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Figure 2.21: Visualization of the temporal development (left to right) of the turbulence structure using iso-surfaces of Q∗ colored by density for the multi-fluid post- shock turbulence. These structures are captured immediately after the shock wave. (a) vorticity-dominated structure, and (b) strain-dominated structure. . . 84 Figure 2.22: Contributions to the vortex stretching rate from particles starting in each quadrant. The particles are initialized uniformly at the post-shock position k0x ≈ 0.44 and traced downstream till the vorticity returns to an isotropic state. (a) single-fluid and (b) multi-fluid cases. . . . . . . . . . . . . . . . . . . 85 Figure 2.23: PDFs of (a) (DQ/Dt)/(cid:104)Qw(cid:105)3/2 and (b) (DR/Dt)/(cid:104)Qw(cid:105)2 for fluid particles with different densities at streamwise location of k0x ≈ 0.5. . . . . . . . . . . . 87 Figure 2.24: Conditional mean rate of change vectors of (DQ/Dt/(cid:104)Qw(cid:105)3/2,DR/Dt/(cid:104)Qw(cid:105)2) in the (Q, R) plane for (a) isotropic turbulence, (b) single-fluid post-shock tur- bulence, and (c) multi-fluid post-shock turbulence at streamwise location of k0x ≈ 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 . . . . . . . Figure 2.25: Conditional mean vectors in the (Q, R) invariants plane for (a) light fluid, (b) medium density fluid and (c) heavy fluid at streamwise location of k0x ≈ 0.5. Figure 2.26: Contributions to the transport equations of the VGT invariants by different terms for isotropic turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term, and (d) viscous term. . . . . . . . . 90 . 89 Figure 2.27: Contributions to the transport equations of the VGT invariants by different terms for single-fluid post-shock turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term, and (d) viscous term. . 92 Figure 2.28: Contributions to the dynamics of the VGT invariants by different terms for multi-fluid post-shock turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term, and (d) viscous term. . . . . . . 93 Figure 2.29: Contributions from pressure Hessian to the dynamics of the VGT invariants in (a) light fluid region, (a) medium density fluid region and (c) heavy fluid region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 xiv Figure 2.30: Single particle dispersion rate normalized using t2 for both single and multi- fluid post-shock turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Figure 2.31: The temporal development of the skewness µ3 of the streamwise particle displacement PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Figure 2.32: Lagrangian velocity autocorrelation function for particles staring immediately after the shock wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Figure 2.33: Lagrangian velocity autocorrelation function for particles staring immediately after the shock wave in the multi-fluid STI simulations. The heavy fluid particles (dashed line) and light fluid particles (dotted line) are shown for comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 . . . . . Figure 2.34: Evolution of the pair particle dispersion rate normalized by t2 for both single- fluid (black) and multi-fluid simulations (red). The direction of the initial separation is in the spanwise direction with lyz = l0. The relative particle dispersion in three direction are separated into streamwise direction (dotted lines), longitudinal direction (solid line) and transverse direction (dashed lines). The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Figure 2.35: Evolution of the pair particle dispersion rate normalized by t2 for both single- fluid (black) and multi-fluid simulations (red). The relative particle dispersion in three direction are separated into streamwise direction (dotted lines), lon- gitudinal direction (solid line) and transverse direction (dashed lines). The direction of the initial separation is in the spanwise direction with lx = l0. The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 . . . Figure 2.36: Evolution of the pair particle velocity correlation for both single-fluid (black) and multi-fluid simulations (red). Velocities in the streamwise direction (solid lines), longitudinal direction (dashed line) and transverse direction (dotted line) are compared separately. The direction of the initial separation is in the spanwise direction with lyz = l0. The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. . . . . . . . . . . . . . . 104 Figure 2.37: Evolution of the pair particle velocity correlation for both single-fluid (black) and multi-fluid simulations (red). The direction of the initial separation is in the spanwise direction with lx = l0. Velocities in the streamwise direction (solid lines) and transverse directions (dotted and dashed lines) are compared separately. The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 xv Figure 3.1: Shock front in wrinkled shock regime (a,b) and broken shock regime (c). Shock intensity is represented by density jump ∆ρ. Shock intensity ratio is calculated as ∆ρ/∆ρ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 . Figure 3.2: (a) Demonstration of numerical setup of the 1D simulation for both shock- (b) Space-time diagram of the capturing and shock-resolving simulations. 1D shock propgation (thick solid line) in non-uniform density media. In the meantime, a reflected wave/forward-going wave will be generated (thin solid line), which will then interact with the post-shock non-uniform density media at the particle path (dashed line). As a result, a re-reflected wave/backward- going wave will be generated (dotted line), which will catch up with the shock wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Figure 3.3: Shock wave propagation in the medium with linearly varying density for density ratios of 8 and 0.125. (a) M0 = 1.1 and (b) M0 = 4.0. . . . . . . . . . . 122 Figure 3.4: Evolution of shock Mach number in linearly varying density field for ρ1/ρ0 = 8 and 0.125. Initial shock Mach number M0 is 4.0. . . . . . . . . . . . . . . . . 123 Figure 3.5: Plots for flow variables in streamwise directions for a Mach 3 shock propa- gating in a sinusoidal shaped density field. (a) Pressure, (b) density and (c) temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 . . . . Figure 3.6: Shock propagation in a fluctuating density field. (a,b) M0 = 2.0, density ratio 2.0, and ks = 1 and ks = 4 and (c) M0 = 1.1, density ratio 2.0, and ks = 4. . . . 126 Figure 3.7: Comparison of the wavelength ratios between numerical results and theoret- ical results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Figure 3.8: Comparison of the shock position fluctuations between numerical results and theoretical results. Results from shock-resolving simulations are maked using symbols for Ms = 1.1 (black circle), Ms = 1.3 (black cross), Ms = 1.5 (black star), Ms = 2.0 (red circle), Ms = 2.5 (red corss), Ms = 3.0 (red star). . . . . . 128 Figure 3.9: Comparison of the relative magnitude zb/z f of the re-reflected waves.Results from shock-resolving simulations are maked using symbols for At = 0.1 (black circle), At = 0.2 (black cross), At = 0.333 (black star), At = 0.4 (red circle), At = 0.5 (red corss), At = 0.6 (red star), At = 0.9 (blue circle). . . . . . 130 Figure 3.10: Demonstration of the catchup of the re-reflected waves. . . . . . . . . . . . . . 132 Figure 3.11: Shock propagation in sinusoidally varying density field. The temporal de- velopment of shock position is compared for the shock-resolving simulation, CW theory and the proposed CWRW theory. . . . . . . . . . . . . . . . . . . . 133 xvi Figure 3.12: Shock propagation in random density media. Different realizations of the random density field are represented using different colors.(a) ks = 4.0 and (b) ks = 8.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 . . . . . Figure 3.13: Shock propagation in 3D sinusoidal shaped density media. The temporal development of shock position on certain point on the shock surface is plotted. kx = 4 for all the cases and ky = kz unless the values are indicated by the legend.136 Figure 3.14: The dependence of mean shock speed on wavenumber ratio for shock propa- gation in 3D non-unform density media with kx = 4. . . . . . . . . . . . . . . . 137 Figure 3.15: The dependence of mean shock speed on wavenumber ratio for shock propa- gation in 3D non-uniform density media with kx = 4. . . . . . . . . . . . . . . 139 Figure 3.16: The dependence of L2 norm of the shock fluctuations on wavenumber ratio for shock propagation in 3D non-unform density media with kx = 4. . . . . . . 140 Figure 3.17: The mean propagation of shock wave in 3D random density media. Red lines represent points on the shock surface. The blue and black curves represent the bounds at the two asymptotic limit. . . . . . . . . . . . . . . . . . . . . . . 141 xvii KEY TO SYMBOLS AND ABBREVIATIONS Reλ Taylor Reynolds number Ms Flow Mach number Mt Turbulent Mach number M0 reference Mach number Pr Prandtl number Sc Schmidt number At Atwood number MWi molecular weight of ith specie µi chemical potential of ith specie δ shock thickness δn numerical shock thickness η Kolmogorov length scale ρ density ui flow velocity p pressure s entropy E total energy R specific gas constant T temperatures Y mass fraction of heavy fluid τi j, σi j viscous stress tensor thermal diffusion vector qi qdi enthalpy diffusion vector Ji mass fraction flux vector µ dynamic viscosity xviii Cp specific heat at constant pressure Cv specific heat at constant volume γ specific heat ratio k0 peak wavenumber of the turbulence energy spectrum λ Taylor micro scale λφ scalar Taylor micro scale λB Batchelor scale ε dissipation tp pass-over time ωi vorticity Rαα Normal Reynolds stress ai normalized mass flux b density specific volume covariance i particle position x+ i particle displacement relative to the mean flow y+ i particle velocity u+ Ai j velocity gradient tensor Si j strain rate tensor Wi j rotation tensor θ dilatation P first invariant of the velocity gradient tensor Q second invariant of the velocity gradient tensor R third invariant of the velocity gradient tensor Qw second invariant of the rotation tensor Qs second invariant of the strain rate tensor Ψ strain-enstrophy angle Σ vortex stretching rate xix Hp i j pressure hessian tensor Hb i j baroclinic tensor Ti j viscous contribution tensor to VGT dynamics l0 initial particle pair separation distance lx initial particle pair separation in streamwise direction lyz initial particle pair separation in spanwise directions z wave intensity calculated by the pressure ratio across the wave xs shock position DNS Direct Numerical Simulation RANS Reynolds Averaged Navier-Stokes LIA Linear Interaction Approximation STI Shock-Turbulence Interaction TKE turbulent kinetic energy RMI Richtmyer-Meshkov Instability MP Monotonicity Preserving VD variable density RHS right hand side LHS left hand side IT isotropic turbulence R-H Rankin-Hugoniot ICF Inertial Confinement Fusion VGT velocity gradient tensor CMT conditional mean trajectory SFS stable-focus/stretching UFC unstable-focus/contracting SN/S/S stable-node/saddle/saddle UN/S/S unstable-node/saddle/saddle CW Chisnell-Whitham xx CHAPTER 1 GENERAL BEHAVIOR OF VARIABLE DENSITY TURBULENCE INTERACTION WITH A NORMAL SHOCK WAVE 1.1 Introduction The interaction between a normal shock wave and isotropic turbulence is an important funda- mental problem that has been extensively studied. Understanding the physics behind this problem is beneficial to many applications such as hypersonic combustion, inertial confinement fusion and astrophysics. However, the existence of a wide range of length/time scales and other complicating effects in flows involving both turbulence and shock waves have posed serious challenges on the study of these flows. Ribner (1954) proposed a theoretical model for the description of Shock-Turbulence Interaction (STI), in which Eulerian Rankine-Hugoniot equations were solved. In this theory, referred to as Linear Interaction Approximation (LIA), the pre-shock turbulence was assumed to be made of small amplitude disturbances, so nonlinear and viscous effects were excluded from the analysis. This allowed the changes of different modes across the shock to be separately investigated. The amplification of turbulent kinetic energy (TKE) and the reduction of turbulence length scales by the shock are successfully predicted by LIA and they are treated as the upper limits for the interaction with nonlinear and viscous effects. Early Direct Numerical Simulation (DNS) of shock-isotropic turbulence by Lee et al. (1993) (flow Mach number Ms (cid:54) 1.2) showed reasonably good agreement with LIA for the vorticity amplification when the turbulent Mach number, Mt, was kept small. However, due to the low resolution of the simulations, quantities which peak behind the shock were not compared to LIA. Jamme et al. (2002) and Mahesh et al. (1997) have also used DNS to study the effects of different upstream flow parameters on STI. Lee et al. (1997), Larsson and Lele (2009), Larsson et al. (2013) have conducted turbulence-resolving simulations for different flow Mach numbers using a high order shock-capturing scheme and found good agreement with 1 LIA results for some of the statistics, like turbulent kinetic energy and vorticity variance, but a poor match of individual Reynolds stress components with LIA was observed. More recently, Ryu and Livescu (2014) have conducted shock-resolving DNS for a wide range of parameters and have shown that DNS results converge to LIA when the ratio of shock thickness, δ, to Kolmogorov length scale, η, becomes small, which was achieved by decreasing Mt. This work has established the reliability of LIA as a prediction tool for low Mt STI, even at low upstream Reynolds numbers. LIA and shock-capturing simulation data have been used for a detailed study of the turbulent energy flux in the post-shock region by Quadros et al. (2016). Using LIA to alleviate the need to resolve the shock allows studying of post-shock turbulence at arbitrarily high shock Mach numbers (Livescu and Ryu, 2016). However, Shock-LIA studies (following the notation from Livescu and Ryu, 2016) do not address the decay away from the shock or situations where the nonlinear and/or viscous effects can have an effect on the interaction. On the other hand, shock-resolving DNS studies are limited in the range of Reynolds and shock Mach numbers achievable by today’s computers. Turbulence-resolving shock-capturing simulations could extend the parameter range and flow complexity (Jammalamadaka et al., 2014) achievable by shock-resolving DNS, provided that the reliability of the tool can be established for such problems. This is also addressed in this paper. The early theoretical study of small amplitude turbulent fluctuations conducted by Kovasznay (1953) showed the existence of three different components in compressible turbulence: vorticity, acoustic and entropic. The effect of upstream entropy fluctuations on the shock-turbulence inter- action was further studied by Mahesh et al. (1997) and Jamme et al. (2002). Their results show that a negative correlation between the upstream velocity and temperature will cause a higher am- plification of turbulent kinetic energy and vorticity variance. The presence of acoustic fluctuations in the upstream flow is found to cause less amplification of turbulent kinetic energy but more amplification of transverse vorticity variance (Jamme et al., 2002; Mahesh et al., 1995). Hannappel and Friedrich (1995) observed in their study that the turbulent fluctuations of the transverse vor- ticity increase more and those of streamwise vorticity increase less due to the interaction with the 2 shock wave when the inflow compressibility increases. Shock waves also enhance mixing (Menon, 1989; Kim et al., 2003), a well-known effect used for increasing the fuel-air mixing in practical combustion/propulsion systems (e.g. Budzinski et al., 1992; Yang et al., 1993). In the above STI studies, the pre-shock turbulence is single-fluid, so the density variations are due to either thermodynamic or acoustic fluctuations. In order to understand the process of hypersonic combustion where strong density variations exist in the fuel-air mixture, the multi-fluid variable density effects (i.e. due to compositional variations) should be taken into consideration. In a hypersonic reacting flow, the species mixing enhancement is affected by the density variation and expected to be different from that of a passive scalar. For example, strong mixing asymmetry exists in variable density turbulence (but not observed in single-fluid flow), as predicted by Livescu and Ristorcelli (2008) and Livescu et al. (2010). Such mixing asymmetry is further amplified after passing through the shock wave as shown in figure 1.1. For a system like this, where nonlinear effects become dominant due to the strong density variations, LIA becomes ineffective, leaving high-order numerical simulation as the feasible approach. Li and Jaberi (2012) have recently developed a new shock-capturing finite difference method based on a Monotonicity-Preserving (MP) (Suresh and Huynh, 1997) scheme and have tested the method for different supersonic flows. This work is based on their method, extended to multi-component flows. Our previous study has shown general characteristics of multi-fluid STI (Tian et al., 2017d). The goal of current study is to develop a better understanding of the density effects on the shock-turbulence interaction and scalar mixing. The configuration addressed here, where the shock is stationary and the turbulence is fed through the inlet of an open-ended domain, is also relevant to the Richtmyer-Meshkov Instability (RMI) as it represents a statistically-stationary version of this problem. The classical RMI problem describes the transient linear and the subsequent nonlinear growth of interface and mixing layer (see Zabusky, 1999; Brouillette, 2002) , while the current study resembles the re-shock problem where variable density effects are important for the shock-driven mixing (e.g. Lombardini et al., 2011). In the case where the upstream density field is isotropic, converged statistical data can be obtained with significantly reduced computational cost. 3 The rest of the paper is organized as follows. The governing equations are presented in section 1.2 and the simulations parameters in section 1.3. In order to establish the accuracy of current simulations, the idea of scale separation between Kolmogorov length scale and shock thickness (Ryu and Livescu, 2014) (in shock-capturing simulations it becomes numerical shock thickness δn) is adopted to limit the viscous effects during the interaction, and thus determine whether the interaction is correctly resolved. Combining grid convergence tests with theoretical analysis, we prove in subsection 1.4.1 of the main results section, section 1.4, that the statistics are independent of the mesh and the simulations are reliable. LIA convergence tests are conducted in subsection 1.4.2 to show that LIA limits can be approximated using shock-capturing simulations when certain criteria are satisfied. The effects of density variations on STI are investigated by comparing current results with that of a single-fluid simulation in subsection 1.4.3. For detailed analysis of mechanisms and physics behind this problem, the structure and budgets of important turbulence quantities are examined in subsections 1.4.4 and 1.4.5. The simulation data are then used to address mixing and modeling in variable density turbulence in subsection 1.4.6 and 1.4.7. Finally, section 1.5 presents the conclusions of the study. 1.2 Governing Equations and Numerical Methodology In this work, the following dimensionless compressible Navier-Stokes equations for continuity, momentum, energy and mass fraction, in a binary fluid system following the ideal gas equation of state, are solved numerically: ∂U ∂t + ∂(F − Fv) ∂x1 + ∂(G − Gv) ∂x2 + ∂(H − H v) ∂x3 = 0, p = ρRT γM2 0 (1.1) (1.2) The solution vector U = {ρ, ρu1, ρu2, ρu3, E, ρY} contains the primary variables: density, ρ, velocity components, ui, total energy, E, and mass fraction of the heavy fluid, Y. t, x1, x2, x3 are time 4 and the three coordinate axes. In equation 1.2, R and T denote the specific gas constant and tem- perature. M0 is the reference Mach number, which is generated from the non-dimensionalization. The inviscid flux vectors (F,G,H) and the viscous flux vectors (Fv,Gv,H v) are defined as (cid:110) (cid:110) (cid:110) F = G = H = and ρu1, ρu2 ρu2, ρu1u2, ρu2 ρu3, ρu1u3, ρu2u3, ρu2 1 + p, ρu1u2, ρu1u3,(E + p)u1, ρYu1 2 + p, ρu2u3,(E + p)u2, ρYu2 3 + p,(E + p)u3, ρYu3 (cid:111) (cid:111) (cid:111) , , (1.3a) (1.3b) (1.3c) (1.4a) (1.4b) (1.4c) Fv = {0, τ11, τ21, τ31, uiτi1 − q1 − qd1,−J1} , Gv = {0, τ12, τ22, τ32, uiτi2 − q2 − qd2,−J2} , H v = {0, τ13, τ23, τ33, uiτi3 − q3 − qd3,−J3} The viscous stress tensor τi j, thermal diffusion vector qi, enthalpy diffusion vector qdi and mass fraction flux vector Ji are calculated by the following equations: (cid:18) ∂ui ∂xj τi j = µ Re0 + ∂uj ∂xi − 2 3 ∂uk ∂xk δi j (cid:19) (1.5) (1.6) (1.7) (1.8) qi = −µCp (γ − 1)M2 0 Re0Pr ∂T ∂xi qdi = −µT (γ − 1)M2 0 Re0Sc ∂Cp ∂xi Ji = −µ Re0Sc ∂Y ∂xi 5 Figure 1.1: Scalar structure in multi-fluid turbulence interacting with a Mach 2 shock identified by the iso-surface of heavy fluid mole fraction and colored with instantaneous density fluctuations. The black plane located in the middle of the domain represents the instantaneous shock surface. The iso-surfaces correspond to mole fraction values of (a) 0.05 (b) 0.95 (c) 0.25 and (d) 0.75. where Re0 is reference Reynolds number. The ratio of specific heats is γ = Cp/Cv = 1.4. The fluid viscosity is assumed to be dependent on dimensionless temperature as µ = T0.75 for all species in the flow. The Prandtl and Schmidt numbers are Pr = Sc = 0.75. This set of non-dimensional fully compressible equations are solved numerically in the conservative form. The inviscid fluxes are computed by the fifth-order MP scheme as described in (Li and Jaberi, 2012) with the exception that, here, the mass fraction equation is also computed in the conservative form using the MP scheme. The fluxes are divided into forward and backward parts via Lax-Friedrichs flux splitting method after projecting them into the characteristics field. Both parts are then reconstructed using the MP scheme and combined before projecting back to the physical space. The same tolerance value (10−10) for the MP scheme as in Li and Jaberi (2012) is used for this study. The viscous and diffusive fluxes are calculated by the sixth-order compact scheme (Lele, 1992) coupled with the fifth-order one-sided compact boundary scheme of Cook and Riley (1996). Time advancement is achieved via the classical 3rd-order Runge-Kutta scheme. 6 Figure 1.2: Instantaneous contours of 3D pressure fluctuations and shock surface resulting from isotropic turbulence interacting with a Mach 2 shock. (a) Iso-surface of heavy fluid mole fraction (φ = 0.25) colored by instantaneous pressure fluctuations. Instantaneous shock surface colored by the ratio of local pressure jump across the shock for (b) multi-fluid A (see definition in section 1.4.3) and (c) single-fluid cases. 1.3 Simulations Setup and Parameters Figure 1.2(a) shows the three-dimensional (3D) iso-surface of heavy fluid mole fraction from the base multi-fluid STI simulation, which is colored by local pressure fluctuations. This high- lights the setup of the current simulations: the dimensions of the computational domain are (L1, L2, L3)=(4π,2π,2π); the streamwise direction is denoted by x (or x1) and transverse directions are denoted by y and z (or x2 and x3); the normal shock is nearly stationary and is initialized in the middle of the domain consistent with the laminar Rankine-Hugoniot relations. Periodic boundary conditions in transverse directions are implemented and a buffer layer is set at the end of the domain from x = 4π to x = 6π. Turbulence is assumed to be homogeneous in transverse directions and isotropic before the shock wave. Averages are computed over homogeneous directions to obtain statistics of the flow. Reynolds averages are denoted by an overbar, f , while Favre averages are denoted by a tilde, (cid:101)f ; the corresponding fluctuations around these averages are denoted by f (cid:48) and f (cid:48)(cid:48). To provide inflow conditions, isotropic turbulence fields are superposed on a uniform mean 7 1/2/(cid:113) i flow with Mach number of 2.0 and then advected through the inlet using Taylor’s hypothesis. At relatively low turbulent Mach numbers Mt = u(cid:48) iu(cid:48) γRT, the Taylor’s hypothesis is a good approximation for upstream isotropic turbulence. Moreover, the turbulence is allowed to develop before reaching the shock wave so as to achieve a realistic state. We use the same method as that of Ristorcelli and Blaisdell (1997), and generate the turbulence by a separate temporal simulation of decaying isotropic turbulence. The velocities for the isotropic box simulations are initialized randomly following a 3D Gaussian spectral density function. The mean velocity is set to zero. The peak value of the energy spectrum is k0 = 4.0. The temperature field is uniform initially and the initial pressure is calculated by solving the Poisson equation. The simulation is then conducted until the skewness of velocity derivative reaches around −0.5 and the flatness reaches around 4.0. The range of Mt values of the final fully developed turbulence is 0.03-0.38 and that iu(cid:48) u(cid:48) i/3λ/µ, is 9-45. The Taylor micro of Reynolds number based on Taylor micro scale, Reλ = ρ scale is computed from viscosity, µ, and turbulent energy dissipation, , as λ = (15µu(cid:48)2 rms/)0.5 and u(cid:48) iu(cid:48) rms = (u(cid:48) i/3)0.5. These isotropic turbulence boxes are taken from the same temporally decaying isotropic turbulence simulations at different times. (cid:113) The variable density (multi-fluid) effects arise from variations in composition field, by correlat- ing the density to an isotropic scalar field (heavy fluid mole fraction or mass fraction). The scalar field was generated as a random field following a Gaussian spectrum with a peak at ks = 8.0 and has double-delta probability density function (PDF) distribution, so that the scalar value is either 0.0 or 1.0. The scalar field was then smoothed out by solving the pure diffusion equation, until the scalar gradients became well resolved on the mesh. The resulting scalar field was allowed to decay in fully-developed isotropic turbulence simulations for about one eddy turnover time as a passive scalar (denoted as φ). The density field is then calculated by making X = φ for case A (where X is the mole fraction of the heavy fluid) or Y = φ for case B and using the ideal gas equation of 0/RT. The ratio of the molar masses of the two fluids state in its multicomponent form, ρ = pγM2 is 1.78, resulting in an Atwood number, At = (MW2 − MW1)/(MW2 + MW1), of 0.28. This value of the Atwood number was chosen such that the variable density effects are non-negligible, yet 8 the interaction with the shock wave is still in the wrinkled-shock regime. Higher At values cause significant distortions of the shock front and even breaking of the shock, depending on the speed of sound in the light fluid. A ‘buffer’ layer is set at the end of the domain from x = 4π to x = 6π to prevent reflections of waves from the outflow boundary. In this region, the flow variables are smoothly dissipated to a laminar solution through a damping function, effectively controlling the unphysical oscillations around the outflow (Larsson and Lele, 2009). A constant back pressure calculated from Rankine- Hugoniot equations is applied as the outflow boundary condition. This outflow boundary condition has been shown to cause small movement of the shock wave (Larsson and Lele, 2009). However, the Mt value is low in the final simulations, and the shock drifting speed is calculated to be less than 0.2% of the mean flow speed, so we are confident that the shock movement does not affect the statistics presented. 1.4 Results and Discussions In this section, the variable density (VD) effects on the STI and the mechanisms behind the interaction are discussed in details. However, before discussing these effects, the accuracy of the simulations is established by grid convergence tests and analysis of various flow variables. Convergence to LIA is then studied for shock-capturing simulations and the criteria for convergence are identified. The comparison between multi-fluid and single-fluid cases is then made for various turbulent statistics. A series of in-depth analyses of turbulence budgets for important turbulence and mixing quantities like turbulent kinetic energy and variance of mole fraction are presented for better understanding of VD effects on STI. To eliminate the statistical variability, most results are space averaged over homogeneous directions and time averaged for more than two pass-over times after the flow reaches a statistically steady state, which is achieved after the acoustic wave propagate from the outlet to the shock wave (Larsson et al., 2013). The pass-over time is calculated as tp = 2π/u1,u + 2π/u1,d, where u1,u and u1,d are the pre-shock and post-shock mean streamwise velocities. Instantaneous results are 9 Mesh info 1024 × 2562 1024 × 3842 1536 × 3842 1536 × 5122 2048 × 5122 Grid-1 Grid-2 Grid-3 Grid-4 Grid-5 ∆xd/∆y Reλ Mt 0.09 0.09 0.09 0.09 0.09 2.0 1.3 2.0 1.5 2.0 45 45 45 45 45 (kmaxη)min 0.96 1.45 1.45 1.93 1.93 η/δn 0.93 0.93 1.40 1.40 1.86 Table 1.1: Details of the simulations used in the grid convergence tests. also presented in the form of constant property surfaces and energy spectra when needed. For all the results presented in this section, the coordinate is shifted such that the shock is approximately located at k0x = 0.0. 1.4.1 Accuracy of Numerical Results In the shock-turbulence interaction problem, the grid resolutions used for both turbulence and shock wave are critical to the accuracy of results. Here, the inflow turbulence is generated separately using 2563 grid points and has a kmaxη value around 2.3, where kmax is the maximum turbulence wavenumber and η is the Kolmogorov length scale. kmaxη is usually taken to be greater than 1.5 to resolve small-scale turbulence. It is relatively easy to construct a mesh to fully resolve turbulence in the pre-shock region. However, the post-shock small-scale turbulence is expected to decrease in size, especially in the shock normal direction. Therefore, to ensure that the smallest post-shock turbulence scales are well resolved, a detailed grid convergence test was conducted using 5 different meshes, as shown in table 1.1. The calculation of kmaxη is based on the largest grid size among all three directions to give safer and more conservative estimates. In figure 1.3, several large and small scales turbulent statistics obtained from the convergence test are shown. The region of unsteady wrinkled shock movement is marked as a grey area in this and following figures. Due to the unsteady shock movement and wrinkled shock surface, the averaged shock thickness is much larger than the instantaneous numerical shock thickness. As the mesh gets refined from Grid-1, Grid-3 to Grid-5 with the same compression ratio, all turbulent 10 Figure 1.3: Results of multi-fluid grid convergence tests at Reλ = 45 and Mt = 0.09. (a) Turbulent kinetic energy, (b) Kolmogorov length scale η, (c) streamwise Taylor micro scale λ1 and (d) transverse vorticity variance ω(cid:48) 2ω(cid:48) 2. The region of unsteady wrinkled shock movement is marked in grey. Figure 1.4: Results of multi-fluid grid convergence tests at Reλ = 45 and Mt = 0.09. (a) mass fraction variance φ(cid:48)φ(cid:48) and (b) Batchelor scale λB. 11 k0x-505101520TKE00.050.10.150.2Grid-1Grid-2Grid-3Grid-4Grid-5k0x-505101520η00.0050.010.0150.02Grid-1Grid-2Grid-3Grid-4Grid-5k0x-505101520λ100.050.10.150.20.250.3Grid-1Grid-2Grid-3Grid-4Grid-5k0x-505101520ω′2ω′20102030405060Grid-1Grid-2Grid-3Grid-4Grid-5(a)(b)(c)(d)k0x-505101520φ′φ′00.010.020.030.040.05Grid-1Grid-2Grid-3Grid-4Grid-5k0x-505101520λB00.020.040.060.080.10.12Grid-1Grid-2Grid-3Grid-4Grid-5(a)(b) Figure 1.5: Relation between numerical shock width, δn, and grid size in the shock region. statistics converge, proving the sufficiency of grid resolution. Results from meshes with different compression ratios, i.e. Grid-5 and Grid-4, are compared to make sure that the grid is fine enough to resolve the decreased length scale in x-direction. Statistics that characterize the behavior of mass fraction are also examined. Figure 1.4 shows the mass fraction variance and the Batchelor scale, representing large scale and small scale statistics of mixing, respectively. Again, these statistics are grid converged and the errors are small when using Grid-5. The convergence rate of the Batchelor scale is slower than for all the other statistics, a consequence of VD effects. Larsson et al. (2013) have investigated how the large scales are affected by the finite com- putational box size. Their results show that the effects of box size are minimal for the range of parameters analyzed. For current study, the same peak wavenumber k0 = 4 for the energy spectra and the same box size are used, so the effects from finite computational box size should also be negligible. Another important issue in STI simulations is whether the interaction between shock wave and turbulence is well captured. When a shock capturing scheme is used, a scale separation between shock width and turbulence length scales is desirable. As suggested by Ryu and Livescu (2014), 12 ∆x00.0050.010.0150.020.025δn00.0050.010.0150.020.0250.030.0350.04 the ratio of Kolmogorov length scale to numerical shock width, δn, is considered to be an important parameter to assess this separation, where δn is calculated as (u1,u − u1,d)/|∂u1/∂x1|max, and |∂u1/∂x1|max denotes the maximum magnitude of streamwise velocity gradient. In the coarsest simulation using Grid-1, the ratio of pre-shock Kolmogorov length scale to numerical shock width is estimated to have a value of about 0.93. However, since the Kolmogorov length scale is calculated using average dissipation, locally, turbulence eddies can become commensurate with the numerical shock thickness. Indeed, using instantaneous dissipation values to estimate the minimum turbulence length scale, yields a ratio with the numerical shock thickness that is much smaller. In order to achieve a better scale separation, the numerical shock width needs to be reduced to a much smaller value than that corresponding to the Grid-1 simulation. The dependence of δn on the grid size is investigated by conducting a series of tests using a range of grid sizes around the shock region. The results are shown in figure 1.5, which indicates that the numerical shock width is linearly correlated with the grid size, as expected. Therefore, by refining the mesh in x-direction, the scale separation can be controlled. Results from Grid-2 and Grid-3 show good resolution of post-shock Kolmogorov length scale with a kmaxη value around 1.45, but considerably different post-shock prediction of ω(cid:48) 2ω(cid:48) 2. This difference can be attributed to the effects of scale separation ratio (0.93 for Grid-2 and 1.40 for Grid-3) as shown in table 1.1. Further increasing the scale separation ratio has a small effect on the resolution of post-shock statistics, e.g. from Grid-4 (1.4) to Grid-5 (1.86). To confirm this, we have also examined another set of simulations with lower Reλ but wider range of scale separation ratios, and noticed that having a larger than 1.4 scale separation ratio has minimal effect on the statistics presented here. These results are not shown here for brevity. Based on these results, we can safely say that for the final simulation with Grid-5, which has a scale separation ratio of 1.86, the scale separation ratio is large enough to correctly resolve STI for the current shock-capturing method. 13 1.4.2 Convergence to LIA for Single-fluid Case In the fully resolved STI simulations of Ryu and Livescu (2014), it is shown that as the ratio of shock thickness, δ, to Kolmogorov length scale, η, becomes small, the DNS results converge to the LIA solutions. This suggests that the scale separation between δ and η can be a criterion for controlling the viscous effects on the interaction. In Livescu and Ryu (2016), high Reynolds number / high shock Mach number post-shock turbulence data are generated using the LIA procedure, which resolved the issue of high or prohibitive computational cost in fully resolved DNS. For shock-capturing simulation, such convergence has not been investigated in previous studies. In DNS, the scale separation is controlled by the ratio δ/η , which can be calculated as δ/η = 7.69Mt/(Re0.5 (Ms −1)) and it was varied in Ryu and Livescu (2014) by changing Mt, which also minimized the non- linear effects. In shock-capturing simulations, δn/η, Mt, Reλ, and Ms are, in general, independent parameters; however, the same general principle can be applied to study the convergence of the numerical results to LIA. λ Unlike DNS, where an overlap between shock width and turbulence scales does represent a physical problem, albeit different than the LIA limit, shock-capturing simulations need to have separation of scales to ensure that the numerical shock-capturing algorithm does not alter the physics of the problem. For current simulations, the numerical shock thickness depends on the mesh size in streamwise direction instead of any turbulent statistics, so it can be independently varied without changing Mt and Reλ. The shock thickness for different shock-capturing schemes is affected by the amount of numerical dissipation introduced around shock. It is shown in Li and Jaberi (2012) that the MP5 scheme has less numerical dissipation than the WENO scheme, so the numerical shock thickness is smaller for MP5. Generally, the shock jump is represented by 2 to 3 grid points and the numerical shock thickness calculated using (u1,u − u1,d)/|∂u1/∂x1|max is around 1.6∆x. By changing ∆x, the scale separation ratio can be independently varied. Figures 1.3 and 1.4 (figures correspond to multi-fluid cases, but the ideas are the same) show that for a ratio η/δn greater than 1.40, there is no significant error in the post-shock statistics. Compared to shock-resolving simulations, this feature of shock-capturing simulations relaxes the requirement of 14 fully resolved shock profile and makes the simulations cheaper without diminishing the accuracy of turbulent statistics. Larsson (2010) took a different approach towards this problem. A theoretical model of the error introduced by the shock-capturing scheme was used to predict the post-shock turbulence statistics. The error of post-shock Reynolds stresses was found to scale with grid spacing to the second order. Using this model, a critical value of k0∆x can be calculated to control the error of Reynolds stresses for a certain shock-capturing scheme. Compared to this method, the scale separation criterion naturally takes into account the artificial dissipation added in the shock region by using numerical shock thickness as a shock length scale, which makes it a simpler criterion. Typically, a scale separation ratio η/δn ≥ 1.4 is recommended for current shock- capturing simulations. We also recognize that, even though the scale separation ratio provides guidance in conducting shock-capturing turbulence resolving simulations, the specific value for the scale separation criterion is not universal. It depends on the shock-capturing method and the way the numerical shock is quantified. While δn/η is still an important parameter to assess the accuracy of the results, it is unlikely that the results in the shock region have the same physical meaning as those from shock-resolving DNS. Therefore, to study the convergence to LIA in shock-capturing simulations, Mt and Reλ need to be considered separately, instead of being one single parameter as in shock-resolving DNS. Single-fluid simulations are conducted covering a wide range of Mt and Reλ parameter space to study the convergence to LIA. Reλ immediately upstream of the shock varies between 9 and 45 and Mt right before the shock ranges from 0.03 to 0.38. More details of the LIA convergence tests can be found in table 1.2. The minimum scale separation in these tests is δn/η = 1.4 to ensure accuracy. 2u(cid:48) 1 and u(cid:48) 1u(cid:48) In figure 1.6, the streamwise and transverse Reynolds stress components (u(cid:48) 2) obtained from the single-fluid simulations are compared with the LIA solution. All the plots are normalized by the values immediately upstream of the shock wave. The effects of Reλ are first considered and the Mt values are kept around 0.24 (except for Reλ = 45 case, for which Mt = 0.09). The results for u(cid:48) 1u(cid:48) 1 are shown on figure 1.6(a). As Reλ increases, the post-shock peak values of 15 1 converge to LIA prediction in some region behind the shock. For u(cid:48) 1u(cid:48) u(cid:48) 2u(cid:48) 2, the post-shock value reaches the peak right after the shock wave, and then it monotonically decreases. The peak right after the shock can be affected by shock thickness, shock corrugation and unsteady shock movement, so it is not a good representation of the shock amplification of turbulence. As suggested by Ryu 2 at corresponding u(cid:48) and Livescu (2014), the values of u(cid:48) 2u(cid:48) 1u(cid:48) 1 peak positions are examined instead 2u(cid:48) to evaluate the convergence. The trend of converging to LIA prediction for u(cid:48) 2 is shown in figure 1.6(b), and it qualitatively matches with DNS results. But compared to u(cid:48) 1u(cid:48) 1, the convergence of u(cid:48) 2u(cid:48) 2 is slower. This can be explained by the fact that Mt values in these simulations are not low enough for u(cid:48) 2u(cid:48) 2 to approximate the LIA solution. As Mt increases, the shock wave becomes more wrinkled and it strongly affects the convergence of u(cid:48) 2u(cid:48) 2, since this quantity is more sensitive to shock wrinkling (Lee et al., 1997). To further examine the role of the shock wrinkling, the effects of Mt on the convergence to LIA predictions are examined by setting the Reλ immediately before the shock wave around 30 and varying Mt. In figure 1.7, the normalized plots of Reynolds stresses are compared and it can be seen that the amplification ratios of u(cid:48) 1u(cid:48) 1 for all the cases are very close (but different than the LIA limit) and those of u(cid:48) 2u(cid:48) 2 converge to LIA as Mt decreases. This indicates that for current shock capturing 2u(cid:48) simulations, the nonlinear effects induced by high Mt are more important for u(cid:48) 1u(cid:48) 2 than for u(cid:48) 1. The results are also consistent with our previous statement that u(cid:48) 2u(cid:48) 2 is more sensitive to Mt. When 2 exhibits the same converging trend, but u(cid:48) 2u(cid:48) comparing to DNS results, u(cid:48) 1u(cid:48) 1 does not converge to the LIA limit at Reλ = 30. It seems the LIA limit of u(cid:48) 1u(cid:48) 1 cannot be approximated by decreasing Mt at low Reλ and Reλ needs to be large enough for convergence to be achieved. This is confirmed by another test with low Mt (0.09) and higher Reλ (45) as shown in figure 1.6. For Reλ greater than 45 and Mt lower than 0.1, the match to LIA results is fairly good, suggesting a minimum value of Reλ around 45 needed for the streamwise Reynolds stress component to converge. In summary, the results show that shock-capturing simulations exhibit similar converging trend to LIA solutions as shock-resolving DNS. LIA solutions for individual Reynolds stresses can be approached provided that there is some separation of scale between turbulence scales and shock 16 LIA convergence case Reλ Mt 0.09 0.24 0.24 0.24 0.03 0.09 0.17 0.24 0.38 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8 Case 9 Case 1(final case) 45 30 16 9 30 30 30 30 30 Table 1.2: Cases considered for LIA convergence tests. 1), and Mt is low enough (for u(cid:48) 1u(cid:48) width, Reλ is large enough (for u(cid:48) 2u(cid:48) 2). On the other hand, the approach to the LIA prediction is different for the streamwise and spanwise components of the Reynolds stress tensor. For the transverse component, the converging trend agrees very well with the DNS results as we vary either Mt or Reλ. However, the streamwise component does not converge to the LIA prediction at low Reynolds numbers by varying Mt, which is different from the DNS study by Ryu and Livescu (2014). This is attributed to the fundamental differences between shock-capturing turbulence-resolving simulations and shock-resolving DNS. A further test then shows that the LIA limit of u(cid:48) 1u(cid:48) 1 can be approached at both large Reλ (45) and low Mt (0.09), which represents a regime that viscous and nonlinear effects are not important for the prediction of u(cid:48) 1u(cid:48) 1. Based on these findings, it would be safer to extend the single-fluid STI to multi-fluid STI at higher Reλ. The requirement on Reλ makes it relatively more expensive for shock-capturing simulation to approach LIA than full DNS. However, shock-capturing simulations can still achieve a large scale separation between Kolmogorov length scale and physical shock thickness using coarser grids than shock-resolving DNS. In practice, using coarser grid for simulating the shock should dominate over the requirement of increased grid resolution for the turbulence to be able to capture a higher Reynolds number flow. Therefore, shock-capturing simulations should be an effective tool for studying STI so long as the results have been carefully verified as shown above. 17 Figure 1.6: Comparison of Reynolds stress obtained from single-fluid simulations with LIA results. Mt ≈ 0.24 except for Reλ = 45. (a) Streamwise Reynolds stress and (b) transverse Reynolds stress. Figure 1.7: Comparison of Reynolds stress obtained from multi-fluid, single-fluid simulation with LIA results. Reλ ≈ 30. (a) Streamwise Reynolds stress and (b) transverse Reynolds stress. 18 k0x05101520u′1u′100.511.522.53LIAsingle-fluid:Reλ=45single-fluid:Reλ=30single-fluid:Reλ=16single-fluid:Reλ=9k0x05101520u′2u′200.511.522.53LIAsingle-fluid:Reλ=45single-fluid:Reλ=30single-fluid:Reλ=16single-fluid:Reλ=9(a)(b)k0x05101520u′1u′100.511.522.5LIAsingle-fluid:Mt=0.03single-fluid:Mt=0.09single-fluid:Mt=0.167single-fluid:Mt=0.237single-fluid:Mt=0.38k0x05101520u′2u′200.511.522.5LIAsingle-fluid:Mt=0.03single-fluid:Mt=0.09single-fluid:Mt=0.167single-fluid:Mt=0.237single-fluid:Mt=0.38(a)(b) 1.4.3 Effects of Density Variations on STI In this section, the effects of density variations on STI are examined in detail by comparing the results obtained from variable density or ‘multi-fluid’ cases with those obtained from a reference ‘single-fluid’ simulation. As demonstrated in the previous section, a minimum value Reλ ≈ 45 is needed for streamwise Reynolds stress to converge to LIA prediction. At lower Reλ, the convergence trend of u(cid:48) 1u(cid:48) 1 is different from fully resolved DNS studies (Ryu and Livescu, 2014). Therefore, the value of Reλ = 45 was chosen for the results presented in the rest of the paper. These results were obtained using Grid-5 from table 1.1 and are grid converged, as shown in section 1.4.1. For the multi-fluid case, two different ways are used to generate density variations at inflow: (a) the heavy fluid mole fraction is correlated with the random scalar used for initialization and (b) the heavy fluid mass fraction is correlated with the random scalar. In the following discussions, these two cases will be referred to as the multi-fluid A and multi-fluid B cases. The random scalar fields used to generate density in both multi-fluid cases are the same and have symmetric PDFs. Evidently, this initialization of density field will make the density PDF of multi-fluid A case symmetrical and that of multi-fluid B case asymmetrical. Below, the multi-fluid A case will be treated as the default case, because it is a better representation of most mixing processes, in which the initial volumes occupied by the fluids in the mixture are specified. Moreover, the reactants in a stoichiometric combustion process are commonly introduced based on their mole fractions, so a specified mole fraction field at the inlet is preferred. The multi-fluid B case is also of interest to us, as it features a different inflow density structure (for the same density ratio) compared to the default case. The single-fluid reference simulation was conducted using the same inflow conditions for turbulent variables except density. In this reference case, the mass fraction of the heavy fluid is set to 1.0. At the same time, a passive scalar equation, which is the same as the mass fraction equation in the multi-fluid case, is solved for comparison. This case is referred to as just the single-fluid case and used as a reference to study the effects of VD on STI. For all cases, the turbulence is allowed to adjust itself to the scalar field in the pre-shock region before interacting with the normal shock. In figures 1.2(b) and (c), the instantaneous shock fronts colored by the ratio of local pressure 19 Figure 1.8: 2D contours of instantaneous scalar and density fields for (a,b) multi-fluid A simulation and (c,d) single-fluid simulation. The 2D contours are taken at the same time step, in planes x-z. jump as obtained from multi-fluid A and single-fluid cases are compared. As expected, the shock wave of multi-fluid case becomes more wrinkled and has stronger pressure jump variation. Figure 1.8 shows the instantaneous 2D contours of density and scalar for multi-fluid A (heavy fluid mole fraction) and single-fluid (passive scalar) cases. After scaling the contours to the same range, figures 1.8(a) and (c) show that the density variation is much stronger in the multi-fluid case due to the compositional change. In contrast, the variation of density in the reference (single-fluid) case results from thermodynamic fluctuation, which makes it very small for the parameters considered here. When examining the passive scalar fields and mass fraction fields (figure 1.8(b) and (d)) for each case, we observe that the scalar fields have similar structure and magnitude before the shock wave as they are generated similarly at inflow. After passing through the shock wave, the scalar fields and mixing are different. Evidently, when the density variation is significant in pre-shock turbulence, the mixing enhancement by shock is stronger. Figures 1.9-1.15 show several important flow statistics obtained from the multi-fluid and single- 20 fluid simulations. These statistics are gathered by averaging over homogeneous directions and time, so that they only depend on the streamwise direction. Averaged flow statistics are compared in figure 1.9 for single-fluid and multi-fluid A cases. Before the shock wave, all cases yield the same results. This observation not only confirms that the inflow conditions are somewhat similar in these cases, but also implies that for current simulation, the effect of density variations on turbulence is small in the pre-shock region. When comparing the multi-fluid turbulent kinetic energy and vorticity variance with the corresponding single-fluid values in figures 1.9(a) and (d), it is noted that the amplification in these turbulent statistics is much more significant in the multi-fluid cases. Furthermore, the multi-fluid turbulent kinetic energy reaches a peak around k0x ≈ 2.0, which is closer to the shock than k0x ≈ π for single-fluid case. Figures 1.9(b) and (c) show the comparison for streamwise turbulence Taylor micro length scale, λ1, and Kolmogorov length scale, η. The reduction in turbulence length scales across the shock wave is evident in these figures; the multi-fluid cases show more reduction than the single-fluid case. Note that the changes in turbulence statistics in multi-fluid cases are expected to depend on the scalar structure and the Atwood number (e.g. Lombardini et al., 2011); these are not discussed in this paper. In figure 1.10, statistics related to the scalar field (heavy fluid mass fraction for multi-fluid case A and passive scalar for the single-fluid case) and mixing are compared. Both scalar variance φ(cid:48)φ(cid:48) and Batchelor scale λB are shown. λB is calculated based on the scalar dissipation and is the representative of the smallest scales in the scalar field. The scalar statistics are normalized by the values immediately before the shock wave. After passing through the shock wave, the faster decay of scalar variance for the multi-fluid case indicates stronger shock enhancement of scalar mixing. The Batchelor scale, however, shows a more complex behavior. Unlike the Kolmogorov length scale, the same reduction ratio of Batchelor scale across the shock wave is observed between the multi-fluid and single-fluid cases. After passing through the shock wave, the Batchelor scales of multi-fluid cases exhibit a transient process of decreasing before returning to the pre-shock value, during which an even smaller structure of the scalar field is generated. In the single-fluid case, 21 Figure 1.9: Plots of (a) turbulent kinetic energy, (b) Kolmogorov length scale, (c) Taylor micro scale and (d) transverse vorticity variance for multi-fluid A (red) and single-fluid (blue) simulations. Figure 1.10: Plots of normalized (a) scalar variance, (b) Batchelor scale for multi-fluid A (red) and single-fluid (blue) simulations. 22 k0x051015TKE00.050.10.150.2multi-fluidsingle-fluidk0x051015η00.0050.010.0150.020.025multi-fluidsingle-fluidk0x051015λ100.050.10.150.20.250.3multi-fluidsingle-fluidk0x051015ω′2ω′20102030405060multi-fluidsingle-fluid(a)(b)(c)(d)k0x051015φ′φ′00.20.40.60.81multi-fluidsingle-fluidk0x051015λB00.20.40.60.811.2multi-fluidsingle-fluid(a)(b) Figure 1.11: Conditional expectation of scalar as a function of density at various streamwise location for : (a) multi-fluid A and (b) single-fluid cases. The conditional mean is calculated using variables in a spanwise plane at certain streamwise locations and then averaged over time. Peak TKE positions are different between multi-fluid(k0x ≈ 2.0) and single-fluid(k0x ≈ π) cases however, the Batchelor scale increases monotonically back to its pre-shock value. We also note that after k0x ≈ 10.0, the multi-fluid λB values are larger than the single-fluid values as the faster mixing immediately after the shock smooths out the small scalar scales. The different behaviors of scalar variance and Batchelor scale in the post-shock region can be explained by using the scalar transport equation. This will be discussed in Section 4.5. We also note that results from multi-fluid case B show almost the same behavior as those from multi-fluid case A for the above statistics, despite their completely different density PDFs. It was shown in figure 1.1 that the turbulence structure is affected very differently by the shock in regions with different density values. To further understand this behavior, conditional expectations of several turbulence quantities, conditioned on the density, are calculated and examined. Figures 1.11 - 1.15 show these conditional means at various locations for the multi-fluid A and single-fluid cases. These statistics are calculated in selected spanwise planes (y-z planes) and then averaged in time. The conditional mean of scalar, shown in figure 1.11(a), confirms the correlation of density field to the scalar field (heavy fluid mole fraction in this case). The dependence of mole fraction on density is almost linear on average, which suggests that the density fluctuations are mainly caused by compositional variations, instead of thermodynamics fluctuations. After passing through the 23 ρ−ρ-0.6-0.4-0.200.20.40.6φ0.10.20.30.40.50.60.70.80.9rightbeforeshockpeakTKEk0x≈4.0k0x≈9.24ρ−ρ-0.1-0.0500.050.1φ00.20.40.60.81rightbeforeshockpeakTKEk0x≈4.0k0x≈9.24(b)(a) Figure 1.12: Conditional expectation of pressure as a function of density at various streamwise locations for : (a) multi-fluid A close to shock, (b) multi-fluid A at far field and (c) single-fluid cases. shock wave, a decrease in the slope is observed due to enhancement of density variation across the shock wave, which translates into an increased extent of the x axis values. This results from a larger local variation of the shock strength, with the subsequent amplification of the post-shock values in the multi-fluid case. For comparison, figure 1.11(b) shows that the correlation between the passive scalar and density is very low for the single-fluid case, which is expected. Figure 1.12 shows the development of the correlation between pressure and density fluctuations as the flow passes through the shock and moves away from the shock. Little correlation between pressure and density is observed before the shock, as expected. The increase in the correlation between pressure and density can be mostly attributed to the acoustic field. Due to the fast evolving nature of the acoustic field in STI (as evidenced in the rapid decay predicted by LIA), it is worth comparing the pressure-density correlation close to the shock, where the transient effects dominate, with the far field behavior. Figure 1.12(a), shows the development of the conditional mean of pressure close to the shock wave. Right after the shock, at k0x ≈ 0.64, there exists a strong positive correlation, which is the direct effect of shock wave on turbulence. During the post-shock flow/turbulence development, RMS value of pressure fluctuations, p(cid:48) rms, decreases rapidly, indicated by the decreasing range of pressure fluctuations. This process corresponds to the conversion of potential energy to turbulent kinetic energy. This is confirmed by the fact that p(cid:48) rms remains almost unchanged after reaching the peak TKE position. In the far field (as shown in figure 1.12b), the pressure fluctuations become less correlated with density fluctuations in the mixed fluid regions. The loss of pressure-density correlation is hypothesized to be caused by the nonlinearity 24 ρ−ρ-0.500.51P−P-2024rightbeforeshockk0x≈0.64k0x≈0.825peakTKEρ−ρ-0.6-0.4-0.200.20.40.6P−P-1-0.500.511.5rightbeforeshockpeakTKEk0x≈5.0k0x≈9.24ρ−ρ-0.1-0.0500.050.1P−P-1-0.500.51rightbeforeshockpeakTKEk0x≈5.0k0x≈9.24(a)(b)(c) Figure 1.13: Conditional expectation of temperature as a function of density at various streamwise location for : (a) multi-fluid A case and (b) single-fluid case. introduced by the strong density fluctuations and the coupling between different modes. Figure 1.12(c) shows the density-pressure correlation at various positions for the single-fluid case. In contrast to the multi-fluid case, the small perturbation assumption is satisfied, so the mode decomposition by Kovasznay (1953) can be applied. Based on the way the inflow turbulence is generated (Ristorcelli and Blaisdell, 1997), the pre-shock turbulence is dominated by vorticity and acoustic modes. In this case, the pressure and density fluctuations should follow the isentropic process: γ(γ − 1)ρ(cid:48)/ρ = (γ − 1)p(cid:48)/p (Kovasznay, 1953). For current single-fluid case results, the density and pressure fluctuations indeed follow the isentropic process scaling after being normalized by local mean, in agreement with previous work (Lee et al., 1993; Mahesh et al., 1997). After passing through the shock wave, the isentropic scaling is no longer satisfied due to the generation of entropy mode, which is featured by a negative correlation between density and temperature fluctuations (figure 1.13b). This should be contrasted to the study by Mahesh et al. (1997) whose focus is on the interaction between turbulent boundary layer and shock wave. In such case, the inflow fluctuations satisfy the weak form of Morkovin’s hypothesis (Morkovin, 1962), instead of following the isentropic relations. For the current case, despite the generation of entropy mode, the pressure fluctuations cannot be neglected even at k0x = 9.24 and, therefore, Morkovin’s hypothesis is not satisfied throughout the domain. 25 ρ−ρ-0.8-0.6-0.4-0.200.20.40.60.8T−T-0.04-0.03-0.02-0.0100.010.020.030.040.05rightbeforeshockpeakTKEk0x=5.0k0x=9.24ρ−ρ-0.1-0.0500.050.1T−T-0.015-0.01-0.00500.0050.010.0150.02rightbeforeshockpeakTKEk0x=5.0k0x=9.24(a)(b) Figure 1.14: Conditional expectation of TKE as a function of density at various streamwise location for : (a) multi-fluid A and (b) single-fluid case. For a complete discussion on thermodynamic properties, the conditional expectation of temper- ature for both multi-fluid A case and single-fluid case are shown in figure 1.13. For the multi-fluid case, the direct effect of upstream Mach number fluctuations on temperature is strong, resulting in a positive correlation between density and temperature. The generated entropy mode as predicted by LIA in the single-fluid case seems to be overshadowed by this effect even further downstream. For the single-fluid case, the scaling between density and temperature fluctuations satisfies that of an isentropic process before the shock wave. Through the interaction with the shock wave, negative correlation is established as previously mentioned. However, the correlation is not fully linear, which can be attributed to the presence of the acoustic mode and the coupling between this and the entropy mode. In figure 1.14(a), the conditional expectation of the turbulent kinetic energy, TKE, is shown. We note that TKE has a preferential distribution in the relatively high or low density regions in the post-shock regions. One possible explanation is that in high and low density regions, the local sound speed has different values from that of average sound speed, so that the local shock velocity, u1,s becomes nonzero (in the reference frame of the laminar shock wave) and changes significantly. The local movement of the shock surface then further changes the post-shock velocity, u1,d and makes it deviate from the averaged post-shock velocity u1,d. The deviation from u1,d or u(cid:48) i in low 26 ρ−ρ-0.6-0.4-0.200.20.40.6TKE00.10.20.30.40.50.60.7rightbeforeshockpeakTKEk0x≈4.0k0x≈9.24ρ−ρ-0.1-0.0500.050.1TKE0.060.080.10.120.140.160.18rightbeforeshockpeakTKEk0x≈4.0k0x≈9.24(a)(b) Figure 1.15: Conditional expectation of enstrophy as a function of density at various streamwise location for : (a) multi-fluid A and (b) single-fluid case. and high density regions is much larger in magnitude than that in regions with moderate density, which results in larger TKE in the high and low density regions. We have computed the conditional average of u(cid:48) i on density. Results agree very well with our explanation for conditional TKE. We also note that TKE is larger in the light fluid regions compared to heavy fluid regions. This is due to the low inertia of the light fluid, which responds faster to the changes in the local strain and, thus, accelerates faster (Livescu and Ristorcelli, 2009; Livescu et al., 2010) . This explanation is also applicable to figure 1.14(b) for single-fluid case, which shows a preferential distribution of TKE in the lower density fluid region before the shock wave. After passing through the shock wave, a stronger amplification in the high and low density regions is noted but this is relatively weak, which agrees with the previous arguments. The correlation between vorticity and density has a completely different behavior than TKE. Evidently, figure 1.15(a) shows that more vorticity is generated in the mixed fluid regions. Before reaching the shock wave, the mixing process is relatively slow, so that there are still large regions with pure or partially mixed fluids and only narrow regions with fully mixed fluids. In these regions, the density gradients remain large. Through the interaction with the shock wave, the large density gradients lead to the generation of vorticity through the baroclinic torque (see Eq. 1.10). For the single-fluid case (figure 1.15b), the distribution of vorticity is not affected much by the shock 27 ρ−ρ-0.8-0.6-0.4-0.200.20.40.60.8√ωiωi2345678910rightbeforeshockpeakTKEk0x=5.0k0x=9.24ρ−ρ-0.1-0.0500.050.1√ωiωi234567rightbeforeshockpeakTKEk0x=5.0k0x=9.24(a)(b) because of the absence of large density variations. The discussion above using conditional expectations emphasizes the VD effects on the modifi- cation of the turbulence structure. These effects reveal the underlying physics of multi-fluid STI. It’s worth mentioning that the characteristics of conditional expectations for multi-fluid B case (not shown) are similar to those of the multi-fluid A case. 1.4.4 Reynolds Stresses αu(cid:48)(cid:48) The modification of Favre-averaged normal Reynolds stress componentsu(cid:48)(cid:48) α or Rαα (no summa- tion over α) across the shock wave is an important feature of STI. As mentioned in previous section, the pre-shock behavior of Reynolds stresses is very similar for multi-fluid and single-fluid cases and is dominated by viscous dissipation. The variable density effects on the viscous dissipation before interacting with the shock are small at the Atwood number considered here. After passing through the shock wave, the main contributing mechanisms to the transient development of Reynolds stresses are the pressure-velocity correlation and pressure-strain correlation (Lee et al., 1997; Jamme et al., 2002; Larsson et al., 2013). This corresponds to the decay of the acoustic mode generated by STI. However, when the variable density effect is introduced, the dominating terms in the Reynolds stress budgets are unknown, especially after passing through the shock wave. To investigate this, the budgets of the Favre-averaged normal Reynolds stresses need to be considered. Before further discussing the turbulence budgets, we need to differentiate the modification of turbulence by the shock from the evolution away from the shock. For shock-capturing simulations, when there is an overlap in scale, turbulence and TKE budgets in the shock region cannot be correctly captured. However, when there is a scale separation between the shock and the turbulence, as is the case here, the numerical viscous effects become, in principle, negligible during the interaction. Thus, the interaction becomes physical; however, it also becomes uninteresting from the point of view of the budget equations. Nevertheless, the evolution away from the shock is accurately described by the current approach and, for the VD case, has unknown evolution mechanisms that can be elucidated by the budget equations. Thus, the following discussions are focused on the evolution away from 28 ∂(ρ(cid:101)uk Rαα)/2 ∂xk = −ρRαk ∂xk ∂(cid:101)uα (cid:124)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) +p(cid:48) ∂u(cid:48) α (cid:18) ∂p (cid:19) ∂xα term II (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) − ∂ταk ∂xk term I −u(cid:48)(cid:48) ∂xα term V α (p(cid:48)u(cid:48)(cid:48) α) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) − ∂ ∂xα term III (ρu(cid:48)(cid:48) αu(cid:48)(cid:48) αu(cid:48)(cid:48) k/2) − ∂ ∂xk −τ(cid:48) αk (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) ∂u(cid:48) α ∂xk term IV αku(cid:48)(cid:48) (τ(cid:48) α) ∂ ∂xk term VII + term VI the shock, with the immediate post-shock value treated as boundary value. Livescu et al. (2009) derived the Reynolds stress budgets for compressible turbulence as following: (1.9) The left hand side of Eq. 1.9 contains terms describing the advection of Reynolds stresses by mean velocity, and it is balanced by summation of all terms on the right hand side, which includes: (I) production; (II) turbulent pressure-strain; (III) velocity-pressure transport; (IV) viscous dissipation; (V) production due to mean pressure gradient and stress gradient; (VI) turbulent transport, and (VII) viscous diffusion, respectively. Figures 1.16- 1.20 show the contributions of individual terms in the Rαα equation and the summation of all the right hand side terms, denoted by RHS. As mentioned before, all the statistics are time- and space-averaged to achieve converged results. The streamwise Reynolds stress R11 budgets are shown in figure 1.16 (the results of multi-fluid A and multi-fluid B cases are close, so only multi-fluid A results are shown). We note that in the post-shock region, there exists a strong transient process that is dominated by the combination of pressure-strain (II) and pressure-velocity (III) correlations. These correlations are related to pressure and quickly decrease, which is hypothesized to correspond to the fast decay of the acoustic field generated by STI. All the other terms remain small compared to these three terms in the entire post-shock region. Figure 1.16(b) gives the balance of R11 transport terms. The good balance of advection and RHS in both shock zone and post-shock zone confirms the accuracy of the simulation. After the initial strong transient process, the pressure-strain(II) and pressure-velocity (III) terms exhibit long-lasting small magnitude fluctuations and the dissipation (IV) becomes the leading order term. These fluctuations are partly due to the limited size of the inflow turbulence box and can be smoothed by repeating the simulations using additional realizations of the isotropic turbulence box and then averaging the results. In order to remove the effects from limited turbulence box size, 29 Figure 1.16: Contributions of different terms in the R11 transport equations for multi-fluid A simulation. (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. we have conducted 5 more new simulations using different realizations of statistically identical isotropic inflow turbulence and then averaged the budget terms for all these realizations in addition to time. Other methods have also been used to reduce the statistical variability, such as the blending method in Larsson and Lele (2009) or an auxiliary forced IT simulation in Ryu and Livescu (2014). It can be seen in figure 1.17 that the secondary fluctuations are greatly smoothed. All following budget plots in this subsection are averaged over five realizations. In Figure 1.18, the R22 budgets are presented. For R22, the pressure-velocity(III) term vanishes due to averaging in homogeneous directions, but the pressure-strain remains dominant in the transient region. When comparing this term for R11 and R22 budgets, we note a similar behavior, with a quick return to zero while oscillating around x-axis, but opposite sign. After k0x ≈ 2.0, the dissipation starts to become the leading contributing term in the R22 budget just like the R11 budget. The good balance of R22 budget is shown in figure 1.18(b). To further illustrate the effects of multi-fluid density variations on Reynolds stresses, the contributing terms are normalized by the corresponding Reynolds stress values at the same location, as shown in figure 1.19. We note that both multi-fluid and single-fluid values are now of similar magnitude, in contrast to the large differences in their absolute values. The single-fluid budgets agree qualitatively well with those from previous STI studies (Lee et al., 1997; Jamme et al., 2002; 30 k0x-202468-0.500.511.5Pressure-StrainPressure-VelocityDissipationProductionProduction (by mean pressure and stress gradient)Turbulent DiffustionViscous Diffusionk0x-202468-0.500.511.5AdvectionRHS(a)gu′′1u′′1(a)(b)gu′′1u′′1 Figure 1.17: Contributions of different terms in the R11 transport equations for multi-fluid A simulation averaged over 5 different realizations (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. Figure 1.18: Contributions of different terms in the R22 transport equations for multi-fluid A simulation averaged over 5 different realizations. (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. 31 k0x-202468-0.500.511.5Pressure-StrainPressure-VelocityDissipationProductionProduction (by mean pressure and stress gradient)Turbulent DiffustionViscous Diffusionk0x-202468-1-0.500.511.5AdvectionRHS(a)gu′′1u′′1(a)(b)gu′′1u′′1k0x-202468-0.6-0.5-0.4-0.3-0.2-0.10Pressure-StrainDissipationProductionProduction (by mean pressure and stress gradient)Turbulent DiffustionViscous Diffusionk0x-202468-0.6-0.5-0.4-0.3-0.2-0.10AdvectionRHS(a)gu′′2u′′2(a)(b)gu′′2u′′2 Figure 1.19: Normalized contributions of different terms in the R11 transport equations averaged over 5 different realizations for (a) multi-fluid A and (b) single-fluid simulations. Larsson et al., 2013). One difference is the magnitude of the dissipation term, which is relatively small compared to Larsson et al. (2013). To examine the difference, a preliminary simulation in a similar parameter range (Mt ≈ 0.15 and Ms = 1.5) was performed (not shown here) and the results were similar to Larsson et al. (2013). Therefore the difference can be attributed to different shock and turbulent Mach numbers. We note that the pressure-strain and pressure-velocity terms return to zero faster in the multi-fluid case, proving that the transient process for the multi-fluid case is shorter than that for the single-fluid case, which is in agreement with our previous observation. After the first large transient process of post-shock adjustment, we note that in both cases there is a secondary lower-level transient process that lasts very long in the post-shock region, with magnitude close to the dissipation term, especially for the multi-fluid case. This long-lasting small magnitude transient process makes the development of RANS models for these terms challenging. This observation is in agreement with Gréa et al. (2016), even though the jumps themselves can be well predicted by some of the current models (Schwarzkopf et al., 2016). Other than the time-averaged and space-averaged statistics, the instantaneous planar data in homogeneous directions could also provide valuable insights into the turbulence statistics and structure. Figure 1.21 shows the comparison of transverse kinetic energy spectra for single and 32 k0x012345678-2-101234567Pressure-StrainPressure-VelocityDissipationRHSProductionProduction (by mean pressure and stress gradient)Turbulent DiffustionViscous Diffusionk0x012345678-2-101234567Pressure-StrainPressure-VelocityDissipationRHSProductionProduction (by mean pressure and stress gradient)Turbulent DiffustionViscous Diffusion(a)gu′′1u′′1(a)gu′′1u′′1 Figure 1.20: Normalized contributions of different terms in the R22 transport equations averaged over 5 different realizations for (a) multi-fluid A and (b) single-fluid simulations. multi-fluid cases at different locations. The positions selected include: right before the shock (k0x ≈ −0.72) and at k0x ≈ 4.0. Again, because the shock surface is wrinkled, the location of the plane cannot be right before the shock everywhere, which will affect the accuracy of the calculations. In order to gather the spectrum information right before the shock, the instantaneous shock wave position and shape are first calculated. Then the spectra are taken at the plane closer overall to the shock wave. In figure 1.21, the normalized spectra are compared at positions before the shock and at k0x ≈ 4.0 for each case. Right before the shock wave, spectra for single- and multi-fluid cases are nearly identical. After passing through the shock wave, at k0x ≈ 4.0, the small scale energy is enhanced for both cases, agreeing well with early studies by Lee et al. (1993). However, the amplification of the small scales is stronger in the multi-fluid case at the same spatial location. 1.4.5 Vorticity Variance Vorticity is another important variable that is strongly affected by STI. Earlier studies on STI have 2ω(cid:48) 3ω(cid:48) 2 (or ω(cid:48) proven that the transverse vorticity variance ω(cid:48) 3) is amplified immediately after the shock wave, while the streamwise component ω(cid:48) 1ω(cid:48) 1 stays unchanged. This is then followed by 1 and a decrease in ω(cid:48) 1ω(cid:48) an increase in ω(cid:48) 2ω(cid:48) 2 as the turbulence returns to isotropy downstream of 33 k0x012345678-3-2-10123Pressure-StrainDissipationRHSProductionTurbulent DiffustionViscous Diffusionk0x02468-3-2-10123Pressure-StrainDissipationRHSProductionTurbulent DiffustionViscous Diffusion(a)gu′′2u′′2gu′′2u′′2(b) Figure 1.21: Two dimensional spectra of transverse turbulent kinetic energy. The comparison are made at right before shock wave and k0x ≈ 4.0. the shock. The amplification of ω(cid:48) 2ω(cid:48) 2 predicted by LIA has been reasonably captured by previous shock-resolving DNS studies (Ryu and Livescu, 2014). Current shock-capturing simulations also show good agreement on the amplification ratio of ω(cid:48) 2ω(cid:48) 2 with LIA (outside shock zone). When the effect of variable density (due to compositional change) is introduced, the modification of vorticity variance by shock can no longer be predicted by LIA. In figure 1.22, the vorticity variance for both single-fluid and multi-fluid cases are plotted. As demonstrated before, the amplification of transverse vorticity variance by the shock is greater when density variations are large. For the streamwise vorticity variance, both cases show values unaffected by the shock immediately after the shock, while the multi-fluid values increase faster downstream. The governing equation for vorticity variance can be derived from the momentum equation and is used here to identify the mechanisms responsible for the change in vorticity. After rearrangement, this equation takes the following form (Lee et al., 1993): 34 k100101102E22(k)+E33(k)10-710-610-510-410-310-210-1multi-fluidk0x≈−0.72multi-fluidk0x≈4.0single-fluidk0x≈−0.72single-fluidk0x≈4.0 Figure 1.22: Streamwise and transverse vorticity variance for multi-fluid and single-fluid case (ω uj ∂ ∂xj (cid:48)2 α ) = 2ω(cid:48) αω(cid:48) j j −ω (cid:48)2 α −ω (cid:48)2 α ∂uα ∂xj +2ω(cid:48) αω(cid:48) ∂u(cid:48) α ∂xj term II /ρ2 term I +2εi jk ω(cid:48) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) term V (cid:48)2 α u(cid:48) +2εi jk ω(cid:48) ∂uk ∂xk term III ∂u(cid:48) k ∂xk term IV ∂τkq ∂xq (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:18) 1 (cid:19) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) term VI ∂u(cid:48) +2ωαω(cid:48) α ∂xj αu(cid:48) −2ω(cid:48) − ∂ ∂xj ∂ωα ∂xj ∂p ∂xk ∂ ρ ∂xj ω α ∂ ∂xj α ρ α term VII term VIII term IX j j −2ωαω(cid:48) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) ∂u(cid:48) k ∂xk α term X (1.10) The terms on the LHS of Eq. (1.10) represents the convection or transport of vorticity variance by the mean velocity and the RHS terms are: (I) production by mean velocity; (II) production by turbulent stretching; (III) vorticity-mean dilatation; (IV) vorticity-turbulent dilatation; (V) baroclinic torque; (VI) viscous diffusion/dissipation; (VII) turbulent transport ; (VIII) mixed velocity-vorticity gradient; (IX) mixed vorticity-velocity gradient; and (X) vorticity-dilatation. 35 k0x-50510152005101520multi-fluidω′1ω′1multi-fluidω′2ω′2single-fluidω′1ω′1single-fluidω′2ω′2 Figure 1.23: Contributions of different terms in the transport equation for (a) ω(cid:48) 1 and (b) ω(cid:48) 1ω(cid:48) 2ω(cid:48) for multi-fluid simulation. 2 1 and (b) ω(cid:48) 1ω(cid:48) Figure 1.24: Contributions of different terms in the transport equation for (a) ω(cid:48) 2ω(cid:48) for single-fluid simulation. 2 The last three terms vanish due to the averaging in homogeneous directions, but they are shown in the equations and figures for completeness. The main contributing terms in Eq. (1.10) are shown in figures 1.23 and 1.24. Previous DNS study by Jamme et al. (2002) considered the vorticity variance budgets inside the shock wave. Here, the post-shock region is of interest. It can be seen that the dominating mechanisms in the post-shock development region of the flow are production by turbulent stretching (II) and viscous dissipation (VI) terms. These two terms have competing effects on vorticity variance. Turbulent stretching production adds a positive contribution while viscous dissipation decreases the vorticity 36 k0x0246810-200-150-100-50050100150200250300AdvectionTurbulent-strechingviscous diffusionRHSk0x0246810-400-300-200-1000100AdvectionTurbulent-strechingviscous diffusionRHS(a)ω′1ω′1(b)ω′2ω′2k0x0246810-30-20-10010203040506070AdvectionTurbulent-strechingViscous diffusionRHSk0x0246810-100-80-60-40-2002040AdvectionTurbulent-strechingViscous diffusionRHS(a)ω′1ω′1(b)ω′2ω′2 variance. The relative importance of these two terms describes different behaviors of ω(cid:48) 1ω(cid:48) 1 and 1ω(cid:48) ω(cid:48) 2. In ω(cid:48) 2ω(cid:48) 1 transport equation, term (II) is more significant than term (VI), resulting in a net increase in ω(cid:48) 1ω(cid:48) 1 . On the contrary, the magnitude of term (VI) is larger than term (II) for ω(cid:48) 2ω(cid:48) 2, causing the spanwise vorticity variance to decrease. It is worth mentioning that, other than the above dominating terms, terms like baroclinic torque (V) or mean velocity stretching (I) also make small contributions to the development of spanwise vorticity variance immediately after the shock, but they virtually have no effect on the streamwise vorticity. Comparison of vorticity variance budgets between multi- and single-fluid cases can be made to understand the effects of density variations. In both cases, production by turbulent stretching (II) and viscous dissipation (VI) terms are the main mechanisms for modeling vorticity in the post-shock region. Generally, the dominating terms in the multi-fluid case have much larger magnitudes than those in the single-fluid case. As for specific vorticity components, both multi-fluid and single fluid cases have increasing turbulent stretching (II) and decreasing viscous diffusion term (VI) close to the shock (figures 1.23(a) and 1.24(a)) for ω(cid:48) 1ω(cid:48) 1. The increase in the magnitude of viscous diffusion (VI) seems to be not as fast as that of turbulent stretching (II). Due to this, the multi- fluid A case shows greater contribution from the turbulent stretching (II) than from the viscous diffusion term (VI) in the vicinity of shock. The increase in the magnitude of the viscous term (VI) is closely connected to the increase of streamwise vorticity variance. For transverse vorticity variance budgets, the turbulent-stretching term (II) has a different behavior; in the multi-fluid case it increases slowly after the shock and in the single-fluid case it decreases almost immediately after shock. The ω(cid:48) 2ω(cid:48) 2 viscous term (VI) for the single-fluid case decreases in magnitude after the shock, which results from the decrease of the vorticity variance. In contrast, the same term for the multi-fluid case remains almost unchanged for a short distance before the decrease in magnitude, indicating a different transient process compared to the single-fluid case. Ribner (1954); Lee et al. (1993); Larsson and Lele (2009); Livescu and Ryu (2016) showed that STI leads to a two-dimensionalization of the flow immediately after the shock, i.e. the flow becomes locally axi-symmetric. The vortex stretching term is therefore the smallest immediately 37 Figure 1.25: Vortex-stretching term in the enstrophy equation, normalized by ω2ω2 and turbulence time scale TKE/ after the shock. This term then increases as the flow evolves away from the shock and returns to a 3-D structure. Following the expression of the vortex-stretching term in Livescu and Ryu (2016), the turbulent stretching term can be decomposed into contributions from eigenvectors and eigenvalues of strain rate. After normalizing the turbulent stretching term in the enstrophy equation using ω(cid:48) 2ω(cid:48) 2 and the turbulence time scale TKE/ (where  is the rate of dissipation), the effects of the eigenvectors and eigenvalues of the strain rate can be isolated. As shown in figure 1.25, the magnitude of turbulent stretching decreases to around zero across the shock wave. During the transient post-shock process, the multi-fluid case exhibits a faster return to 3-D isotropic turbulence structure, indicated by the faster increase in the normalized turbulent stretching term. Further downstream, the multi-fluid A case reaches its peak value much sooner than the single-fluid case. The behavior of the normalized turbulent stretching term indicates that the contributions to the return to 3-D turbulence not only come from increased enstrophy amplification, but also from the change of alignment between vorticity and strain rate tensor eigenvectors. This difference is hypothesized to be the result of the nonlinear effects introduced by large density variations. 38 k0x05101500.511.5multi-fluid A casesingle-fluid case 1.4.6 Mixing When turbulence passes through a shock wave, the enhancement of mixing by STI is usually expected (see figure 1.10). In this section, the mixing characteristics of multi-fluid turbulence interacting with the normal shock wave are studied in more details. When the variable density effect was introduced, the modification of turbulence by the shock was shown to be very different from that of single-fluid case. This leads to changes in basic features of the characteristics of mixing after the shock. Figure 1.26 compares the scalar Taylor micro scale of heavy fluid mole fraction or passive scalar (λφ = (φ(cid:48)φ(cid:48)/(∂φ/∂x1)2)0.5) between multi-fluid and single-fluid cases. Both cases are normalized using the values right before the shock. Immediately after the shock, the same decrease in λφ is observed. These results and those shown in figure 1.10(b) for the Batchelor scale indicate that the modification of scalar length scales by the shock wave is not dependent on the density variations in the pre-shock turbulence. Away from the shock, the multi-fluid flow exhibits a transient region with slowly decreasing λφ until k0x ≈ 3.0. After the transient region, λφ for the mole fraction starts to increase. On the contrary, λφ for the single-fluid case continuously increases after the shock. It was shown in figure 1.10 that the Batchelor scale λB for the multi-fluid case becomes even larger than that of the single-fluid case close to the end of the domain. This is not observed for λφ and the multi-fluid case λφ values remain small throughout the post-shock region. This suggests that the characteristics of the scalar field as affected by the shock and density variations are very different at small and intermediate length scales. The transport equation for Favre-averaged scalar (passive scalar or heavy fluid mass fraction) variance is considered here to study the characteristics of scalar mixing for multi-fluid turbulence interaction with the shock wave. This equation is rearranged in the following form (Livescu et al., 2000): 39 Figure 1.26: Plots of scalar Taylor micro scale for multi-fluid (red) and single-fluid (blue). ∂ ρ(cid:101)uj(cid:103)φ(cid:48)(cid:48)2 ∂xj ∂xj j φ(cid:48)(cid:48) ∂(cid:101)φ = −2ρu(cid:48)(cid:48) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) term I +2φ(cid:48)(cid:48) ∂Mj ∂xj term V −2M(cid:48) j (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) ∂φ(cid:48) ∂xj term II (ρu(cid:48)(cid:48) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) j φ(cid:48)(cid:48)2) − ∂ ∂xj term III (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) +2 ∂ ∂xj j φ(cid:48)) (M(cid:48) term IV (1.11) where Mj = (µ/Re0Sc)∂φ/∂xj. The LHS term in Eq. (4.3) represents the convection of the scalar variance by the mean velocity and the terms on the RHS are: (I) production; (II) scalar dissipation; (III) turbulent diffusion; (IV) viscous diffusion; and (V) compressibility. Figure 1.27 shows that for the multi-fluid case A the most important term is scalar dissipation (II). In the vicinity of the shock wave, there is a transient process where the production (I) and turbulent transport terms ((III)+(IV)) make relatively small contributions to the overall balance. These two terms reduce to zero around k0x ≈ 5.0 and after that the advection is only balanced by scalar dissipation. In contrast, the only dominating term for the single-fluid case is the scalar dissipation 40 k0x051015λφ0.20.40.60.81multi-fluidsingle-fluid Figure 1.27: Contributions of different terms in the scalar variance(cid:103)φ (cid:48)(cid:48)2 transport equations for multi-fluid simulation. (a) Plots of all the contributing terms in the transport equation and (b) balance of LHS and RHS of the transport equation. throughout the post-shock region as shown in figure 1.28(b). To make a better comparison between the two cases, the magnitude of scalar dissipation and advection terms are shown together in figure 1.29. Right after the shock, the same value for scalar dissipation is observed. As the flow moves away from the shock, the scalar dissipation for the multi-fluid case starts to decrease, while that of the single-fluid case keeps increasing. However, after k0x ≈ 3.0, the scalar dissipation in the multi-fluid case starts to increase with a faster rate than that in the single-fluid case. In the region before k0x ≈ 10.0, the scalar dissipation of the multi-fluid flow is higher in magnitude, which indicates enhanced mixing. Another important feature of scalar mixing in the STI configuration is that, unlike velocity field statistics, the direct effects of the shock on the scalar statistics are almost the same in the single- and multi-fluid cases. This can be explained by the fact that the Rankine-Hugoniot relation for the scalar is linear: φu=φd . This feature of scalar transport prevents direct changes in the scalar field by the shock. This trend exists so long as the pre-shock turbulent Mach number Mt and density ratio in the multi-fluid case is low. These criteria are satisfied in our reported simulations, so the modifications of scalar statistics across the shock remain the same regardless of the VD effects. Furthermore, if certain scalar statistics only uses the data in y-z plane, that statistics is continuous across the shock (e.g. scalar variance). For completeness, we need to mention that the evolution of 41 k0x051015-0.12-0.1-0.08-0.06-0.04-0.0200.020.04ProductionScalar DissipationCompressibility related termTransport termk0x051015-0.12-0.1-0.08-0.06-0.04-0.0200.020.04AdvectionRHS(a)gφ′′φ′′(b)gφ′′φ′′ Figure 1.28: Contributions of different terms in the(cid:103)φ (b) single-fluid simulations. (cid:48)(cid:48)2 transport equations for (a) multi-fluid and the scalar away from shock is indeed affected by the VD effects. The structure of the scalar field can be studied by examining the scalar spectra. Figure 1.30 shows the comparison of the scalar (passive scalar or heavy fluid mole fraction) variance spectra at three locations: right before the shock, at peak TKE, and at k0x ≈ 4.0. The spectra right before the shock wave match well with each other. However, downstream of the shock at k0x ≈ 4.0, the multi-fluid case shows more amplification of small-scale fluctuations (figure 1.30a). The more significant amplification of these fluctuations is accompanied by smaller scalar length scales λφ and λB and higher rate of scalar dissipation. Comparison of the scalar variance spectra at the peak TKE position shows almost the same profiles for the single- and multi-fluid cases (figure 1.21b). We already argued that the modification of scalar statistics (obtained by averaging over planes perpendicular to the mean flow direction) by the shock is not very sensitive to the pre-shock turbulence. The collapse of scalar variance spectra immediately after the shock at peak TKE position implies the scalar evolution away from the shock is mainly controlled by the post-shock turbulence. A mixing asymmetry in variable density turbulence was observed for Rayleigh-Taylor turbulence and its mechanisms are studied by Livescu and Ristorcelli (2008, 2009); Livescu et al. (2010). Due to the differential accelerations experienced by the fluids with different molecular weights, pure 42 k0x051015-0.12-0.1-0.08-0.06-0.04-0.0200.020.040.06ProductionScalar DissipationCompressibility related termTransport termRHSk0x051015-0.12-0.1-0.08-0.06-0.04-0.0200.020.040.06ProductionScalar DissipationCompressibility related termTransport termRHS(a)gφ′′φ′′(b)gφ′′φ′′ Figure 1.29: Comparison of scalar dissipation and advection terms obtained from multi-fluid and single-fluid simulations. Figure 1.30: 2D spectra of scalar. The comparison are made at various locations: (a) right before shock wave and k0x ≈ 4.0 and (b) right before shock wave and peak TKE position. 43 k0x051015-6-5-4-3-2-1012multi-fluid advectionmulti-fluid dissipationsingle-fluid advectionsingle-fluid dissipationgφ′′φ′′k100101102Sc(k)10-610-410-2multi-fluidk0x≈−0.72multi-fluidk0x≈4.0single-fluidk0x≈−0.72single-fluidk0x≈4.0k100101102Sc(k)10-610-410-2multi-fluidk0x≈−0.72multi-fluidpeakTKEsingle-fluidk0x≈−0.72single-fluidpeakTKE(a)(b) heavy fluid mixes slower than the pure light fluid in contrast to the Boussinesq approximation. Such non-Boussinesq effects in turbulence statistics are already observed for current multi-fluid STI results and a similar mixing asymmetry can be observed by comparing the PDFs of scalars (heavy fluid mole fraction, heavy fluid mass fraction and passive scalar) and density. Figures 1.31, 1.32, 1.33 show the PDFs of scalar and density as obtained from the spanwise 2D plane data at various streamwise positions for multi-fluid A (mole fraction), B (mass fraction), and single-fluid (passive scalar) cases. For the multi-fluid case A, figure 1.31 shows that the PDF of the inflow heavy fluid mole fraction is symmetrical. Before reaching the shock, the effect of multi-fluid density variations on the mixing asymmetry is negligible. After passing through the shock wave, however, the initially symmetric mole fraction PDF starts to become skewed and the skewness keeps increasing during the post-shock development. The skewed PDFs suggest that light fluid vanish faster than the heavy fluid. This is also shown in the density PDFs as the low density tails of the PDFs vanish faster. In the multi-fluid case B (figure 1.32), a similar behavior is observed for density PDF evolution. Initially, a symmetrical mass fraction PDF results in a skewed density PDF. The interaction with the shock itself has only a small effect on the asymmetry. Away from the shock, the skewness in the density PDF increases. This asymmetry in the multi-fluid cases is related to the larger kinetic energy in the light fluid regions (see figure 1.14a above). This preferential distribution of TKE suggests that the mixing is stronger in the low density regions. In contrast, figure 1.33 shows that the single-fluid scalar and density PDFs remain symmetrical throughout the domain. i/ρ = −u(cid:48)(cid:48) 1.4.7 Normalized Mass Flux and Density Specific Volume Covariance The normalized mass flux ai = ρu(cid:48) i and density specific volume covariance b = −v(cid:48)ρ(cid:48), where v = 1/ρ, are important quantities for the modeling of mixing in compressible turbulent flows with significant inhomogeneity and density variations (Schwarzkopf et al., 2016). Both of these two terms appear as a closure in the Favre-averaged Reynolds stress equations. Even though in the current study, the terms involving ai and b may not be important to Reynolds transport 44 Figure 1.31: PDFs of multi-fluid A case: (a) scalar (mole fraction) and (b) density fluctuation at various locations. The PDFs are calculated at specified y-z planes and then averaged over time. Figure 1.32: PDFs of multi-fluid B case: (a) scalar (mass fraction) and (b) density fluctuation at various locations. The PDFs are calculated at specified y-z planes and then averaged over time. on average, the contributions from these terms might be significant locally. More importantly, as the density ratio of two fluids increase, the interaction with the shock tends to generate strong inhomogeneity, particularly when the shock becomes unstable and breaks. In these conditions, the behavior and modeling of ai and b will be of great significance. Considering these factors, the behavior of normalized mass flux and density specific volume covariance are further investigated in this section. Figure 1.34 compares a1 and b for both single- and multi-fluid cases. The magnitude of these two terms in multi-fluid turbulence is comparably larger than those in the single-fluid case. In figure 45 xheavy00.20.40.60.8100.511.522.533.5inletrightbeforeshockpeakTKEk0x=5.0k0x=9.24ρ−ρ-0.6-0.4-0.200.20.40.600.511.522.533.5inletrightbeforeshockpeakTKEk0x=5.0k0x=9.24(a)(b)φ00.20.40.60.8100.511.522.5inletrightbeforeshockpeakTKEk0x=4.0k0x=7.1ρ−ρ-0.6-0.4-0.200.20.400.511.522.533.5inletrightbeforeshockpeakTKEk0x=4.0k0x=7.1(a)(b) Figure 1.33: PDFs of single-fluid case: (a) passive scalar and (b) density fluctuation at various locations. The PDFs are calculated at specified y-z planes and then averaged over time. 1.34(a), the pre-shock values of normalized mass flux are shown to be small. After passing through the shock wave, the streamwise normalized mass flux is greatly amplified, indicating a strong inhomogeneity in the streamwise direction, generated by the VD effects. This can be explained by the mismatch between the mean density jump in the simulated STI and the predicted value by R-H relation as observed in Larsson et al. (2013) for single-fluid flows. Immediately after the shock, the mean density jump is lower than the laminar R-H prediction. The mean density then slowly returns to the laminar R-H mean density as the turbulence decays downstream of the shock, introducing a small mean density inhomogeneity in the streamwise direction. This is strongly amplified in the multi-fluid simulations, as the wrinkling of the shock surface and variations in shock strength increase the deviation from the R-H values. Further downstream, a1 attains a very high magnitude for a short distance before decreasing, in agreement with previous observation that the mean density returns to R-H value and the flow becomes isotropic. In figure 1.34(b), the density specific volume covariance for both cases are compared. As b is a measure of the mixing state, its behavior is expected to be similar to that of scalar variance (figure 1.10a) outside the shock region, at least when the non-Boussinesq effects are not large. This is confirmed by figure 1.34(b), which shows that b slowly decreases before the shock and much faster after passing through the shock. For the single-fluid case, the density specific volume correlation remains statistically unimportant 46 φ00.20.40.60.8100.511.522.53inletrightbeforeshockpeakTKEk0x=5.0k0x=9.24ρ−ρ-0.1-0.0500.050.10510152025inletrightbeforeshockpeakTKEk0x=5.0k0x=9.24(a)(b) throughout the whole domain. The transport equations for the normalized mass flux components ai, which appear in unclosed form in the transport equations for Reynolds stresses, are (Livescu et al., 2009): (ρ(cid:101)ujai) = (ρai) + ∂ ∂t ∂ ∂xj term I bPi(cid:124)(cid:123)(cid:122)(cid:125) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) +ρv(cid:48) ∂p(cid:48) ∂xi term II (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) (aiaj) ∂ ∂xj +ρ ∂ ∂xj −ρaj (cid:102)(u(cid:48) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) i − ai) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) term III (ρ(cid:48)u(cid:48) j) + ρu(cid:48) iu(cid:48) −( ∂ ∂xj ∂u(cid:48) j ∂xj ) i term V term VI + (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) iu(cid:48) (ρ(cid:48)u(cid:48) j − Ri j) ∂ ρ ∂xj 1 ρ term IV (1.12) where Pi = ∂p/∂xi−∂τi j/∂xj. Figure 1.35 shows the terms in streamwise normalized mass flux equation. The contributing terms are normalized by its local value so that their relative magnitude can be shown. Right after the shock, terms (II) and (VI) are the leading order terms and have opposite sign. In the vicinity of the shock wave, positive net contributions from these two terms are observed, which correspond to the sharp increase of a1 after the shock wave. Further downstream, term (II) decreases and term (VI) increases. After k0x ≈ 2.0, term (II) becomes negative and keeps decreasing. Term (VI) approaches 0 after k0x ≈ 5.0, making term (II) the only dominating term. The overall behavior of a1 is dominated by term (II), while term (VI) contributes to the transient process after the shock interaction and affects the peak position of a1. The good balance of the advection term with all the RHS terms once again confirms the accuracy of the simulations. The density specific volume covariance b, is a measure of the mixing state. This is an unclosed term in the normalized mass flux equations. Its governing equation has the form (Livescu et al., 2009): b +(cid:101)uj ∂ ∂t ∂b ∂xj = 2aj (cid:124)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:125) ∂b ∂xj term I (cid:169)(cid:173)(cid:171)u(cid:48) (cid:170)(cid:174)(cid:172) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) j ρ(cid:48)v(cid:48) ∂ ∂xj +ρ ρ term III (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) +2ρv(cid:48) ∂u(cid:48) j ∂xj term IV (1.13) (cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125) −2aj(1 + b) ∂ ρ ∂xj 1 ρ term II 47 Figure 1.34: Plots of (a) normalized mass flux a1 and (b) density specific volume covariance − (cid:104)v ρ(cid:105) for multi-fluid and single-fluid simulations. Figure 1.35: Contributions of different terms in the normalized mass flux a1 transport equations for multi-fluid simulation. 48 k0x-505101520-0.02-0.0100.010.020.03multi-fluida1single-fluida1k0x-505101520-0.0100.010.020.030.040.05multi-fluidsingle-fluid(a)(b)k0x123456-20-10010203040AdvectionRHSIIIIIIIVVVI Figure 1.36: Contributions of different terms in the density specific volume covariance b transport equations for multi-fluid simulation. Figure 1.36 shows all terms in the transport equation of b (normalized by its local value). The only important term in this equation is term (IV). Terms (I), (II) and (III) have nonzero contributions close to the shock in the post-shock region, but their total contribution to b is negligible. The ‘transient’ region ends around k0 ≈ 3.0; after that term (IV) becomes the only nonzero term and remains nearly constant. It can be concluded that, the modeling of term (IV) is critical to the RANS simulation of variable density mixing in STI. Figure 1.36(b) shows the good balance between the advection and RHS terms. 1.5 Conclusions A hybrid numerical method, which combines a monotonicity-preserving scheme with a compact scheme, is used for high-order numerical simulations of isotropic multi-fluid turbulence interacting with a planar normal shock wave. A reference simulation of single-fluid turbulence is also conducted with similar conditions. The main objective is to study the effects of density variations on the shock- turbulence interaction and mixing. Convergence tests are conducted to establish the accuracy of numerical results using different meshes with a wide range of grid sizes. The computed statistics 49 k0x123456-10-8-6-4-202AdvectionRHSIIIIIIIV are found to be independent of the grid when the turbulence after the shock is well resolved and the scale separation between numerical shock thickness and turbulent scales is adequate. To further establish shock-capturing simulation as an effective tool in studying STI, LIA con- vergence tests are conducted to show that shock-capturing simulations exhibit similar converging trend to LIA predictions. As Reλ increases, both the streamwise and spanwise components of Reynolds stress tensor show converging trends to LIA predictions. Stronger Mt effects on the LIA convergence of the spanwise component are observed. The streamwise component, however, shows no converging trend toward LIA predictions at low Reλ. Further results show LIA limits for streamwise component can be approached using turbulence-resolving shock-capturing simulations with a minimum Reλ around 45 and a coarser grid compared to shock-resolving DNS. This advan- tage over DNS makes the shock-capturing simulations cheaper and at the same time able to achieve accurate predictions. The comparison between two multi-fluid cases and a single-fluid case is made to identify and explain the effects of major density variations on the flow structure and turbulence statistics in the STI configuration. Due to the strong nonlinear effects introduced by density variations, the modification of turbulence statistics by the normal shock is shown to be different from that of the single-fluid case and LIA prediction, with an increased amplification of turbulence variables (like TKE and vorticity variance) and a more significant reduction in turbulence length scales. Turbulent mixing enhancement by the shock is also more significant in the multi-fluid cases. The redistribution of turbulence statistics across the shock wave is another important feature of multi- fluid STI, as reflected in the changes of conditional expectations. For thermodynamic variables, the isentropic relations of fluctuations are satisfied for the single-fluid case before the shock. After passing through the shock wave, both acoustic and entropy modes are important for the correlations of thermodynamic properties. For the multi-fluid cases, the nonlinear effects and coupling between different modes complicates the analysis and further investigation is needed. It is also found that TKE has a preferential distribution on the pure fluid regions after the shock wave, while vorticity amplification is strongest in the mixed fluid areas. This difference is attributed to the different 50 roles density plays for these quantities. The effects from different density PDFs in the pre-shock turbulence are examined by comparing two multi-fluid cases. Both multi-fluid cases exhibit very similar behavior across the shock even though their density fields are different. Other parameters such as At and ks may also affect the interaction, but they are not the focus of current study. A detailed study of turbulent budgets for Reynolds stresses, vorticity variance, and scalar variance is conducted and the dominating mechanisms in the post-shock region are identified. The rapid increase in the streamwise Reynolds stress component after the shock wave is mainly caused by the pressure-strain and velocity-pressure correlation, whereas the decrease in the transverse Reynolds stress is mostly due to pressure-strain term, which are consistent with the observations made in the single-fluid simulation. Away from the shock wave, as transient processes vanish, the viscous dissipation becomes dominating and causes both Reynolds stresses to decrease. After being normalized using the corresponding local Reynolds stresses values, the contributing terms become similar in magnitude. The energy spectra of transverse kinetic energy are different at the same absolute position k0x ≈ 4.0, confirming multi-fluid and single-fluid cases have different post-shock length scales. For the vorticity variance, the dominating terms are turbulence stretching (positive contribution) and viscous diffusion (negative contribution). Turbulence stretching plays an important role in understanding the post-shock behavior: it attains small values right after the shock due to the two-dimensionalization of turbulence; during the return to a three-dimensional structure, the normalized turbulence stretching for the multi-fluid case increases faster and reaches the peak earlier compared to the single-fluid case. The mechanism that is responsible for enhanced mixing has been identified to be the increased scalar dissipation. The modification of scalar statistics, as measured by changes in variables like Batchelor scale, scalar Taylor micro scale, scalar dissipation across the shock wave is shown to be similar in single- and multi-fluid simulations, which is attributed to the linear Rankine-Hugoniot relation for the scalar (φu=φd). The mixing asymmetry is more pronounced after the interaction with the shock for the multi-fluid cases. This can be explained by the preferential distribution of TKE in the low density regions. Finally, the mechanisms that are responsible for the post-shock development of turbulence are identified and 51 used to assess the main characteristics of variable density compressible flows. 52 CHAPTER 2 DENSITY EFFECTS ON THE POST-SHOCK TURBULENCE STRUCTURE AND LAGRANGIAN FLOW DYNAMICS 2.1 Introduction In the first chapter, we have examined the turbulence statistics, turbulence budgets, conditional statistics, and energy spectrum in the multi-fluid STI and found that the nonlinear effects from the density variations significantly change the turbulence properties in both physical and spectral spaces. There have been few other detained study of multi-fluid or variable density STI in the literature. For example, Jin et al. (2015); Huete et al. (2017) considered a reactive shock wave in a premixed mixture and used LIA and shock-capturing simulations to study the detonation-turbulence interaction. However, there still exist many gaps in our knowledge of the variable density effects on the post-shock turbulence structure and flow topology. In this chapter, we will present a comprehensive study of the density effects on the post-shock turbulence structure from two perspectives: 1) the velocity gradient tensor and its dynamics and 2) Lagrangian statistics of non-inertial fluid particles. The properties of the velocity gradient tensor (VGT) determine a wide variety of turbulence characteristics, such as the flow topology, deformation of material volume, energy cascade, and intermittency. Understanding both the VGT field immediately after the shock-wave and its dynamics as the flow evolves away from the shock wave are also crucial to the development of subgrid-scale models that can accurately describe the shock interaction and return-to-isotropy effects. Chong et al. (1990) has proposed an approach to classify the local flow topology and structure using the invariants of VGT. The dynamical behavior of the VGT has been studied for incompressible flows using the Lagrangian evolution of the invariants along conditional mean trajectories (CMT) (Meneveau, 2011). The statistics regarding the invariants of VGT and their Lagrangian dynamics have been used to understand the structure of turbulence in many canonical flows, such as isotropic turbulence, turbulent boundary layer and 53 mixing layers (e.g. Chong et al., 1998; Ooi et al., 1999; Wang and Lu, 2012). Previous studies on single-fluid STI have examined some of the statistics of the PDF of VGT. Ryu and Livescu (2014); Livescu and Ryu (2016) took a step further to investigate the turbulence structure and vorticity dynamics based on the examination of VGT invariants. By taking advantage of the Shock-LIA procedure, they extracted the statistics of VGT and its invariants for a wide range of shock Mach numbers. On the other hand, the dynamics of VGT as the turbulence evolves away from the shock wave cannot be examined using the Shock-LIA procedure. Our preliminary numerical studies on variable density STI have revealed some important features of VGT (Tian et al., 2017c, 2018). The relation between velocity and scalar field has also been studied by Boukharfane et al. (2018) to better understand mixing in this flow configuration. However, these studies have not yet fully revealed the variable density effects on the post-shock turbulence/scalar structure. It is well-known that Lagrangian statistics used to understand the transport of a localized source of passive scalars. Essentially, if a group of marked fluid particles is released from a localized source of passive contaminants, it is well-known that the mean concentration field can be expressed in terms of "on-particle" statistics, and the concentration variance is similarly related to "two-particle" statistics of the relative motion between fluid particle pairs. These Lagrangian statistics have been extensively studied using DNS in the past few years due to the difficulty in experimental measurement of fluid particles. Yeung and Pope (1989) has used DNS database to extract single-particle Lagrangian statistics for isotropic turbulence. Squires and Eaton (1991) has also presesnted one-particle Lagrangian statistics for homogeneous shear flow. Particle-pair dispersion statistics have been reported by Yeung (1994); Ishihara and Kaneda (2002); Yeung and Borgas (2004); Biferale et al. (2005). These Lagrangian statistics have been used for the validation of Richardson’s theory on particle dispersion (Richardson, 1926) and the development of stochastic models (Durbin, 1980; Thomson, 1990; Komori et al., 1991; Pope, 2002), PDF methods (Pope, 1994) and so on. To extend the applicability of these models, researchers have been studying the particle dispersion in other anisotropic and inhomogeneous type of flows. Shen and Yeung (1997); Yeung (1997) studied the particle dispersion in homogeneous turbulent shear. Choi et al. (2004) 54 then reported the Lagrangian statistics in turbulent channel flow. The progress in theoretical and numerical progress in Lagrangian particle study has been reported in the review papers Sawford (2001); Yeung (2002); Salazar and Collins (2009); Toschi and Bodenschatz (2009). This study aims at using the database of the turbulence-resolving shock-capturing simulations of multi-/single- fluid STI generated in the first chapter to: 1) further study the density effects on the turbulence structure immediately after the shock wave, 2) perform the first Lagrangian simulations of this configuration to understand the dynamical behavior of VGT away from the shock, and 3) perform an Lagrangian analyses of the particle dispersion in the post-shock turbulence. The paper is organized as follows. Details of the simulations and the testing conducted to assess the accuracy of results are discussed in section 2.2. Results are presented in section 2.3 and concluding remarks are made in section 2.4. 2.2 Numerics and Accuracy In this section, we first briefly discuss the numerical approach used for shock-capturing turbulence-resolving simulations in our previous study (Tian et al., 2017a), from which we have extracted the VGT statistics addressed in this paper. The configuration of the variable-density STI extension is described next, followed by a discussion of the new Lagrangian simulations used to examine the VGT dynamics away from the shock. 2.2.1 Governing Equations and Numerical approach The conservative form of the dimensionless compressible Navier-Stokes equations for flows with two miscible species (i.e. continuity, momentum, energy, and species mass fraction transport equations) were solved numerically together with the perfect gas law using a high-order hybrid numerical method (Tian et al., 2017a). The inviscid fluxes for the transport equations were computed using the fifth-order Monotonicity Preserving (MP) scheme, as described in Suresh and Huynh (1997); Li and Jaberi (2012); Tian et al. (2017a). The molecular transport terms were calculated using the sixth-order compact scheme (Lele, 1992). The 3rd-order Runge-Kutta scheme 55 Figure 2.1: Instantaneous contours of vorticity and shock surface in isotropic turbulence interacting with a Mach 2 shock. (a) Vortex structure, identified by the Q criteria (Q = 2 (cid:104)Qw(cid:105)) and colored by the mole fraction of the heavy fluid. Fluid particles are initialized as a sheet that spans over the homogeneous directions at a given post-shock streamwise position and allowed to develop with the flow. (b) Visualized particle sheet, convected and distorted by the post-shock turbulence. The instantaneous shock surface is colored by the shock intensity across the shock for (c) single-fluid and (d) multi-fluid cases. was used for time advancement. 2.2.2 Numerical Setup of Eulerian Simulation The numerical setup of the Eulerian simulation has been explained in detail in the first chapter. 2.2.3 Numerical Detials of the Lagrangian Analyses For the current study, we have tracked more than 4.5 million particles that are initialized uniformly at various streamwise positions (cid:174)x0, and calculated various statistics concerning turbulence structure following their trajectories. The aim is to understand the evolution of flow structures following fluid particles. Figure 2.1 (a) marks with red lines a typical streamwise plane where particles are 56 initialized. The particles are then convected by the instantaneous turbulent velocity obtained by turbulence-resolving shock-capturing simulations (Tian et al., 2017a) and moved to another plane marked by the blue lines. At this stage, the initially flat particle sheet is distorted by the turbulence as shown in figure 2.1 (b). The fluid particles are non-inertial and follow the local flow velocity. The corresponding transport equations for particle positions x+ i are: i (t| (cid:174)x0, t0) dx+ dt i (t| (cid:174)x0, t0) = ui(x+ i , t) u+ = u+ i (t| (cid:174)x0, t0), (2.1a) i (t| (cid:174)x0, t0) can be obtained from the Eulerian velocity field ui(x+ (2.1b) i (t| (cid:174)x0, t0) represents the positions of the particles at time t that are initialized at (cid:174)x0 and time where x+ i , t) by t0. The particle velocity u+ interpolation. The interpolation is based on the cubic spline scheme, whose accuracy in predicting particle positions has been studied in Yeung and Pope (1988). In the STI configuration, there is a sharp change of the flow velocity at the shock, which deteriorates the interpolation accuracy. To achieve accurate interpolation of the particle velocity, the domain is partitioned into three different regions as shown in figure 2.1 (a): pre-shock, shock, and post-shock regions. The instantaneous shock surface is identified using the sensor: s = −θ/(|θ| + (cid:104)ωiωi(cid:105)0.5 yz ) > 0.5, where θ = ∂ui/∂xi is the dilatation, ωi = i jk ∂uk/∂xj is the vorticity, and (cid:104)(cid:105)yz represents the instantaneous average over the homogeneous directions. After the instantaneous shock region is identified, the pre- and post-shock turbulence fields can be separated for interpolation. Lagrangian dynamics of particles across the shock wave is not considered in this study. To extract single-particle dispersion statistics in the post-shock turbulence, we have sampled particles that are uniformly distributed at the y − z plane with k0x ≈ 0.64, and followed them to the end of the computational domain. This streamwise location is chosen such that the particles are away from the wrinkled shock wave so that the interpolation scheme is not affected by the shock jump. At the same time, we would like to ensure the particles are close enough to the shock to capture the effects of the post-shock transient process on particle dispersion. On top of that, the flow 57 has reached a statistically steady-state before the Lagrangian particle statistics are being computed. This allows us to initialize the particles in three homogeneous dimensions y, z and t with k0x being fixed. The total lifespan of the each of the particles, tend − tstart, can be approximately calculated as 2π/ud, where ud is post-shock mean velocity. For particle pair statistics, the separation vectors need to be clearly defined due to the anisotropy of the post-shock turbulence, in contrast to that in isotropic turbulence. Due to the homogeneity in the spanwise direction, the separation vectors can be divided into streamwise and spanwise directions. Here we denote the separation vector as (lx, ly, lz), where lx, ly, and lz are the relative distance between the particles. Again, the homogeneity in the spanwise direction allows us to rewrite the separation vector as (lx, lyzsinα, lyzcosα), where lyz is the relative distance in y − z plane and α is the angle between separation vector and z-axis. The only dominating parameters are lyz and lx. The angle α should be irrelevant in particle pair statistics and can be any random value. 2rad. The particle separation vectors can then be For convenience, we have only use α = 0 and π defined as (lx,0,0) and (0, lyz,0) or (0,0, lyz), where the latter two are equivalent. 2.2.4 Grid and Statistical Convergence The accuracy of the numerical results is addressed in this subsection through a series of convergence tests. To ensure that all the turbulence length scales are well resolved, a grid convergence test was conducted in Tian et al. (2017a). Here, we summarize these results for completeness, together with additional convergence results for small-scale quantities. Figure 2.2 shows the velocity dissipation rate ε and mass fraction dissipation rate εφ as a function of the normalized streamwise direction k0x for a series of meshes. The grey regions in the following figures indicate the unsteady shock region, inside which the results are affected by the shock wrinkling and unsteady shock movement. As the grid is refined in all three directions, both quantities display convergence, proving the accuracy of the turbulence database. Another aspect that needs to be considered is the scale separation between the numerical shock thickness δn and the Kolmogorov length scale η as suggested in our previous study (Tian et al., 2017a). With the finest mesh (512×512×2048), the scale separation ratio η/δn is 58 around 1.9, which is sufficient for resolving the interaction between the numerical shock wave and small-scale turbulent motions. Therefore, in the current study, we have obtained all the statistics from the turbulence field based on the finest grid to ensure accuracy. Secondly, LIA convergence tests were conducted following Ryu and Livescu (2014) to show that the shock-capturing simulations can capture the correct limits. Turbulent Mach number Mt and Taylor Reynolds number Reλ were varied for the canonical single-fluid simulations, covering a wide range of parameter space. The shock-capturing simulation results do converge to LIA predictions for individual Reynolds stress components as long as certain conditions are satisfied (Tian et al., 2017a). This was the first time that the asymptotic values for individual Reynolds stresses were approximated using shock-capturing simulations. Statistical convergence is another important factor that needs to be addressed before any further analysis. To reduce the statistical variability, all the results that are based on the Eulerian data are space-averaged over homogeneous directions and time-averaged for around two pass-over times. The averaging is performed after the flow has reached a statistically steady state to eliminate the effects of transient processes (Larsson et al., 2013). For the Lagrangian statistics, the number of fluid particles needs to be large enough for statistical convergence, especially for conditional averaged statistics. We define the conditional averaged quantity X(A, B), which depends on the random variables A and B using the following definition: (cid:104)X|(A = A0, B = B0)(cid:105) = (cid:104)X|(A0 − 1 (B0 − 1 1 2 ∆A), 2 ∆A) ≤ A < (A0 + 1 2 ∆B)(cid:105) 2 ∆B) ≤ B < (B0 + (2.2) where ∆A and ∆B are the bin sizes. The conditional statistics are obtained by ensemble averaging (denoted by (cid:104)(cid:105)) over all the fluid particles that fall into the bins. The number of samples needed to achieve statistical convergence will be examined for different bin sizes. Figure 2.3 and 2.4 show the Dt (cid:105)/(cid:104)Qw(cid:105)2 and convergence of two important conditional Lagrangian statistics (cid:104) DQ their standard deviation (see section 2.3.3 for definitions), depending on the number of particles Dt (cid:105)/(cid:104)Qw(cid:105)3/2, (cid:104) DR 59 Figure 2.2: Results of multi-fluid grid convergence tests at Reλ = 45 and Mt = 0.09. Streamwise development of (a) velocity dissipation rate ε and (b) mass fraction dissipation rate εφ is shown. The region of unsteady shock movement is marked in grey. Grid numbers for Grid 1 to 5 are 256×256×1024, 384×384×1024, 384×384×1536, 512×512×1536, 512×512×2048. in each bin. For the multi-fluid case, we note that the convergence of both conditional means and standard deviations can be achieved when using around 10,000 particles, larger than that needed for the canonical single-fluid case as shown in figure 2.4. This suggests that the variable density effects make the simulations more computationally demanding. The effects of the bin sizes are also examined by comparing three different set of bin numbers 30× 30, 40× 40 and 60× 60 in the (Q, R) phase plane at the same point (3.0,3.0). These bin numbers correspond to the following bin sizes: (1.3,1.3), (1.0,1.0) and (0.67,0.67). Our analysis indicate that the statistics converge to almost the same value when the sample size is large enough. In the present study, we uniformly sampled more than 4.5 million particles and made sure that there are at least 10,000 particles in each sample bin with the number of bins being 40 × 40 ((∆Q, ∆R) = (1.0,1.0)). 2.3 Results and Discussions The variable density effects on the post-shock turbulence structure and dynamics are examined in this section. The results obtained from the multi-fluid STI simulation are compared with those of a reference single-fluid case and standard isotropic turbulence. First, the post-shock turbulence state and its evolution away from the shock wave are examined to identify the variable density effects. 60 k0x01020ε050100150k0x01020εφ00.020.040.060.080.10.12Grid-1Grid-2Grid-3Grid-4Grid-5(a)(b) Figure 2.3: The statistical convergence for (a) (DQ/Dt)/(cid:104)Qw(cid:105)3/2 and (DR/Dt)/(cid:104)Qw(cid:105)2 and (b) their standard deviations conditioned at point (3.0,3.0) in the (Q, R) phase plane for multi-fluid case. The number of bins are 30 × 30 (solid), 40 × 40 (dashed) and 60 × 60 (dotted). Figure 2.4: The statistical convergence for (a) (DQ/Dt)/(cid:104)Qw(cid:105)3/2 and (DR/Dt)/(cid:104)Qw(cid:105)2 and (b) their standard deviations conditioned at point (3.0,3.0) in the (Q, R) phase plane for single-fluid case. The number of bins are 30 × 30 (solid), 40 × 40 (dashed) and 60 × 60 (dotted). 61 number of samples in each bin102103104-0.4-0.3-0.2-0.100.10.20.3number of samples in each bin1021031040.511.522.533.5(a)(b)σ(DRDt)/hQwi2σ(DQDt)/hQwi3/2DRDt/hQwi2DQDt/hQwi3/2number of samples in each bin102103104-1.5-1-0.500.5number of samples in each bin1021031040.811.21.41.61.82(a)σ(DQDt)/hQwi3/2σ(DRDt)/hQwi2DRDt/hQwi2DQDt/hQwi3/2(b) The results are based on temporal- and spatial-averaged statistics obtained from the Eulerian data. Then, the flow topology is analyzed to further understand the post-shock turbulence evolution. The dynamics that dominate the transient evolution of post-shock turbulence structure are examined using the Lagrangian equation of VGT and Lagrangian data collected for sample fluid particles. 2.3.1 Density Effects on the Post-shock Turbulence and Its Evolution Away From the Shock 2.3.1.1 Turbulence State Immediately After the Shock In our earlier study on multi-fluid STI, we have focused on turbulent length scales, Reynolds stresses, vorticity variance, scalar variance, and their budgets (Tian et al., 2017d,a) as the first quantities to examine and compare with the single-fluid STI predictions. It was found that, in the variable-density case, the small-scale turbulent motions are further enhanced after the interaction with the shock. In addition, the vortex stretching mechanism is more suppressed, but it experiences a faster increase during the development away from the shock. In this section, we take a step further to examine the post-shock turbulence structure in details. The PDFs of streamwise and spanwise (averaged over the two periodic directions right after the shock at k0x = 0.5) longitudinal velocity derivatives for both single-fluid and multi-fluid cases are shown in figure 2.5 alongside with the Gaussian distribution as a reference. The non-Gaussian nature of the velocity gradient PDFs and their connection to the energy cascade and intermittency are well documented in the turbulence literature. Here, we note that immediately after the shock wave, the PDF of the spanwise velocity gradient for both cases remains negatively skewed, as in isotropic turbulence. The streamwise component, however, becomes more symmetric and Gaussian-like due to the interaction with the shock wave. This indicates that the energy transfer to small scales is suppressed in the streamwise direction. We also note that the density has a relatively weak effect on the velocity derivatives PDFs since the single-fluid and multi-fluid cases have similar PDFs. The preferential amplification of the transverse components of the rotation and strain rate tensors is a well-known effect in STI and has been extensively studied for the canonical single-fluid flows 62 Figure 2.5: Comparison of the PDFs of the normalized post-shock velocity derivatives with a Gaussian distribution. (Mahesh et al., 1997; Ryu and Livescu, 2014; Livescu and Ryu, 2016). This amplification can lead to an increase in the correlation between the two quantities. When density variations are introduced, Tian et al. (2017a) showed that the transverse components of vorticity variance are further enhanced across the shock wave while the streamwise component retains almost the same value. To better understand the variable density effects on post-shock turbulence, the PDF of the strain-enstrophy angle, Ψ, is considered in figure 2.6. Ψ is calculated using Ψ =tan−1(Si jSi j/(Wi jWi j)), where Si j = 1/2(Ai j + Aji) and Wi j = 1/2(Ai j − Aji) are the strain and rotation tensors and Ai j = ∂ui/∂xj is the velocity gradient tensor. In isotropic turbulence, the PDF of Ψ peaks at 90◦ (Jaberi et al., 2000), indicating a strain dominated flow. In single-fluid post-shock turbulence, the PDF of Ψ exhibits a shift of the peak from 90◦ to around 45◦, as the shock Mach number increases. This has been observed by Livescu and Ryu (2016) and is interpreted as the increase in correlation of strain and rotation. However, in the multi-fluid case, the peak still occurs at relatively large angles and the increase in correlation is not as pronounced as that in the single-fluid case, at the same shock Mach number. Figure 2.6 implies that the rotation and strain are amplified differently by the shock when large density variations are present, which compromises the correlation between the above two quantities. 63 ∂u/∂x-8-6-4-20246810-410-310-210-1100Gaussian∂u/∂xpost-shock(m)∂v/∂ypost-shock(m)∂u/∂xpost-shock(s)∂v/∂ypost-shock(s) Figure 2.6: PDF of the strain-enstrophy angle Ψ in radian for post-shock turbulence. Figure 2.7: Conditional expectation of the magnitude of (a) rotation tensor and (b) strain rate tensor as a function of density before and after the shock. 64 Ψ(rad)00.511.5PDF00.20.40.60.81multi-fluidpost-shocksingle-fluidpost-shockρ−ρ-1-0.500.5111.522.533.544.55pre-shockWijWijpost-shockWijWijρ−ρ-1-0.500.51123456789pre-shockSijSijpost-shockSijSij(a)(b) To further understand the effects of the density variations on the strain and rotation tensors as the flow interacts with the shock wave, the pre- and post-shock values of the conditional expectations of the magnitudes of the two tensors are shown in figure 2.7 for the multi-fluid case. As observed, the correlation between vorticity (and strain rate) with density is completely different before and after the shock wave. Through the shock wave, the amplification of vorticity is stronger in the mixed fluid regions with near average density, but weaker in the pure fluid regions. This is not observed in the single-fluid simulation. One mechanism that might be responsible for this behavior is the baroclinic torque: (∇ρ × ∇p)/ρ2 in the vorticity transport equation. Consider a wrinkled shock surface extracted from the multi-fluid simulation as shown in figure 2.8 (a); a strong pressure gradient ∇p exists through the shock wave as marked in the figure; large density gradients ∇ρ also exist, especially in the mixed fluid regions. Since the pre-shock density field is isotropic, ∇ρ and ∇p can be locally misaligned, becoming a source of vorticity generation through the baroclinic torque. In the pure fluid regions or single-fluid simulation, however, the density gradients are much smaller, so that the cross product of ∇p and ∇ρ is also small. To confirm this, the PDF of the angle between the density gradient in the transverse direction and the vorticity vector is examined in figure 2.8 (b). Before the shock wave, it is noted that the PDF of the cosine of the angle between the density gradient and vorticity vectors is slightly tilted towards the 0-degree angle. After the shock wave, the multi-fluid case exhibits a stronger tendency of vorticity vector being perpendicular to the density gradient. In contrast, the alignment is completely different after the interaction with the shock wave for the single fluid case, as the two vectors are more likely to be parallel to each other. This provides evidence that density gradients and baroclinic torque play important roles in establishing the preferential deposition of vorticity across the shock wave. Figure 2.9 can help visualize the changes in the vortex structures across the shock wave. The vortex tubes are captured using the Q-criterion and are colored by their local density. Figure 2.9 (a) shows the vortex structures for pre-shock multi-fluid isotropic turbulence. For the visualized vortex tubes, there are no identifiable effects from the density variations; the vortex tubes are not preferentially distributed due to the density effects. However, the interaction with the shock has a 65 Figure 2.8: (a) Demonstration of the baroclinic torque generation across the shock wave and (b) PDF of the orientation between the vorticity vector and density gradient in y-z direction immediately after the shock wave. clear effect on the post-shock vortical structures (figure 2.9 b). Immediately behind the shock wave, the vortex tubes are aligned in the spanwise directions, which has been observed in previous STI studies. More importantly, the vortex tubes also get aligned with the density iso-surfaces, meaning that the vorticity becomes perpendicular to the density gradient. This is consistent with the earlier analysis of the baroclinic torque. As a consequence, the post-shock vorticity field will enhance the mixing between adjacent density regions which, in return, will weaken the transverse vortical field and accelerate the return to an isotropic state. This coupling is further explored in the next section. On the other hand, figure 2.7 (b) shows that the magnitude of the strain rate tensor tends to be stronger in the heavy fluid regions and weaker in the light fluid region. This trend is hypothesized to be related to the dependence of shock strength on the pre-shock density. Tian et al. (2017a) showed that shock strength is positively correlated with density. In the highest density regions, the shock compression is also stronger, while it is weaker in the smallest density regions, leading to the observed trend in the amplification of the strain rate tensor magnitude. This trend is different from that observed for the vorticity, which is explained above. As a result, the trend of the strain- enstrophy angle PDF peaking around 45◦, observed in the single fluid case at higher shock Mach numbers is weakened in the multi-fluid case. 66 Figure 2.9: Vortex structures captured using the Q-criterion, colored by density, for multi-fluid (a) pre-shock turbulence and (b) post-shock turbulence. 2.3.1.2 Evolution of Turbulence State Away from the Shock Wave The evolution of variable density turbulence away from the shock wave involves many coupled nonlinear processes, some of which have been discussed through the turbulence budgets analysis in Tian et al. (2017a). In this section, the evolution of turbulence structure will be examined in details. Figure 2.10 shows the development of some of the fundamental turbulence statistics. The evolution of these statistics helps in the understanding of the general characteristics of single- and multi-fluid STI. The enhancement of turbulent kinetic energy (TKE) and dissipation due to the interaction with the shock is a well-known effect in STI. With the introduction of strong density variations, the effect becomes larger, as shown in figure 2.10 (a). Away from the shock wave, TKE of the multi-fluid case exhibits a similar transient behavior compared to the single-fluid case, but it reaches its peak earlier (k0x ≈ 2.0 versus k0x ≈ 3.1). The mechanism that dominates this transient process has been shown to be related to the fast decay of the acoustic waves (Tian et al., 2017a). To highlight the decay of the acoustic waves, figure 2.10 (b) shows the pressure fluctuations variance as a function of the streamwise position. The amplification of pressure fluctuation across the shock wave is noted, agreeing with Sethuraman et al. (2018). The acoustic wave is stronger 67 in the multi-fluid case immediately after the shock wave. This is related to the shock intensity fluctuations induced by the strong density variations (Tian et al., 2017a). As a result, the decay of the acoustic field is also faster for the multi-fluid case, causing a faster increase of TKE. After the post-shock transient pressure adjustment, the multi-fluid case still exhibits larger absolute pressure fluctuations. However after normalizing with ρu(cid:48)u(cid:48), the pressure fluctuations become somewhat similar in magnitude in these two cases. In figure 2.10 (c), the vortex stretching term Σ is decomposed into streamwise and spanwise directions to explore the axisymmetric state and return-to-isotropy of post-shock turbulence. Pre- vious studies (Livescu and Ryu, 2016; Tian et al., 2017a) have shown that the normalized vortex stretching term reaches a low value after passing through the shock wave, indicating a tendency towards an axisymmetric state. Figure 2.10 (c) shows that, without normalization, the absolute values of the vortex stretching terms are magnified in both single- and multi-fluid cases, more so for the spanwise component. The two components then undergo a transient process, where they first increase and cross each other, before the flow returns to an isotropic state. k ω(cid:48) k − δi j In order to quantitatively study the evolution of turbulence anisotropy, we consider here the Reynolds stress anisotropy tensor defined as bi j = u(cid:48) ku(cid:48) j/u(cid:48) iu(cid:48) k − δi j 3 (Livescu et al., 2002). A similar iω(cid:48) anisotropy tensor, di j, can also be defined for the vorticity field, as di j = ω(cid:48) j/ω(cid:48) 3 . Due to the homogeneity in spanwise directions, the diagonal components of the anisotropy tensor are related by b22 = b33 = −0.5b11, so only b11 is discussed. The near zero value of b11 ≈ 0.0 is an indication that flow has reached an isotropic state, while b11 ≈ −1/3 means that the turbulent field has a tendency towards a 2D axisymmetric state. Figure 2.10 (d) shows that d11, a small-scale turbulent variable, reaches around -0.3 in the multi-fluid case, which is lower than that observed in the single-fluid case. This indicates that variable density intensifies the trend towards axisymmetric state for small-scale turbulence. On the other hand, the stronger turbulent stretching mechanism as observed in figure 2.10 (c), makes the return to isotropy much faster in the multi-fluid case as compared to that in the single-fluid case. For Reynolds stresses, which are large scale turbulence statistics, the multi-fluid flow reaches a quasi-isotropic state immediately after the shock wave 68 Figure 2.10: Development of (a) TKE and dissipation, (b) pressure variance, (c) vortex stretching, and (d) anisotropy (b11) of Reynolds stress and vorticity. (b11 ≈ 0.0), while single-fluid turbulence exhibits a tendency towards an axisymmetric state. Evidently, the variable density effects on the post-shock turbulence fields are different at small and large scales. Additionally, the quasi-isotropic state of the multi-fluid turbulence is not stable and is modified in the post-shock transition. Due to the energy transfer between the acoustic field and solenoidal turbulence field, R11 quickly increases, causing b11 to become larger than zero. The anisotropy reaches its maximum value around the peak TKE position (k0x ≈ 2.0) and then slowly decreases. For the single-fluid case, b11 keeps increasing till k0x ≈ 13.0, even though the acoustic effects almost vanish after peak TKE location (k0x ≈ π). In figure 2.11, the developments of skewness and the flatness of the longitudinal velocity gradients are examined before and after the interaction with the shock wave. They show how the non-Gaussian behavior of the velocity field and specifically VGT are affected by the combined shock and density effects. For isotropic turbulence, the skewness of the longitudinal velocity gradient 69 0510150246810ǫ(m)ǫ(s)TKE(m)TKE(s)05101501234p′p′(m)p′p′(s)k0x051015050100150200Σx(m)Σyz(m)Σx(s)Σyz(s)k0x051015-0.4-0.3-0.2-0.100.10.2Reynoldsstress(m)Reynoldsstress(s)vorticity(m)vorticity(s)0510151234(a)(b)(c)(d)prms/(ρu′u′) should be around -0.5, which is observed to be true in the pre-shock region for both single and multi-fluid cases for both streamwise as well as spanwise components (figure 2.11 a). Immediately after the interaction with shock, different components of the derivative skewness tensor are shown to be modified in different ways. The streamwise component for both single-fluid and multi-fluid cases approaches to values very close to 0.0, which is consistent with the tendency towards a two-dimensional axisymmetric state observed above. As the turbulence evolves away from the shock wave, the streamwise velocity derivative skewness decreases rapidly. Due to the strong density variations, the multi-fluid case exhibits a faster decrease before k0x = 5.0, after which it slowly increases towards the −0.54 value. The shock modification of the skewness of the transverse derivative is relatively small for the single-fluid case. For the multi-fluid case, the longitudinal transverse velocity derivative becomes less negatively skewed, with a value of around -0.25. This difference can be attributed to stronger shock intensity variations and shock wrinkling in the multi- fluid case. Away from the shock wave, for both cases, the skewness of ∂v/∂ y increases first until it reaches a peak and then slowly decreases. Comparably, the multi-fluid case exhibits a shorter but more intensified transition. At the end of the domain, however, the spanwise derivative skewness is still larger than −0.5, as the flow is still anisotropic. Figure 2.11 (b) shows the development of longitudinal velocity derivative flatness across and after the shock wave. Immediately after the shock, the flatness of the streamwise component decreases in value while that of the spanwise component increases. Similar to the skewness, the effect of density variations is relatively small on the flatness of streamwise component at the Atwood number considered in this study. On the other hand, the density variations in the multi-fluid case make the increase in transverse component less significant, with the pre- and post-shock values being almost the same. Away from the shock wave, the flatness of the longitudinal streamwise velocity derivative increases, returning to its pre-shock value, while the increase is much faster for the multi-fluid case. For the transverse longitudinal derivative component, the flatness slowly decreases after a small change. From the results above, it can be stated that the variable density effects are not strongly manifested immediately after the shock wave for some of the statistics, but they play an important 70 Figure 2.11: Development of (a) skewness, and (b) flatness, of the streamwise and transverse components of velocity derivatives. Figure 2.12: PDF of the density gradient at different streamwise locations for: (a) single-fluid case and (b) multi-fluid case. 71 k0x051015skewness-1-0.500.51∂u/∂x(m)∂v/∂y(m)∂u/∂x(s)∂v/∂y(s)k0x051015flatness34567(a)(b)∂ρ/∂xi-1-0.500.51012345∂ρ/∂xi-15-10-505101500.10.20.30.40.5∂ρ/∂x1post-shock∂ρ/∂x2post-shock∂ρ/∂x1peakTKE∂ρ/∂x2peakTKE(a)(b) role in the post-shock adjustment. To get an insight into this behavior, density gradient PDFs are examined at various streamwise positions in figure 2.12. Before the shock wave, the PDFs of the density gradients are symmetric in all three directions for both single- and multi-fluid cases (not shown). For the single-fluid case, after passing through the shock wave, the density gradients’ PDFs remain symmetrical, but the streamwise component PDF becomes wider due to the shock compression. As the turbulence develops away from the shock wave to the peak TKE position, the density gradients’ PDFs still remain symmetrical and become narrower, which is related to the fast decay of the acoustic waves. For the multi-fluid case, the density gradients’ PDFs are strongly amplified through the shock wave, but the changes away from the shock are relatively small, because the density variations are controlled by the compositional changes instead of the acoustic waves. More importantly, the streamwise component becomes negatively skewed. We have observed in our previous study (Tian et al., 2017a) that the post-shock streamwise velocity is correlated with the density. Thus, as the flow develops away from the shock wave, both the density and velocity gradients become negatively skewed. This will be further discussed and explained in the next subsection. The nonlinear effects of the shock wave, magnified by large density variations have been shown to cause the preferential deposition of turbulent quantities in regions with different densities (Tian et al., 2017a). This non-uniform turbulence distribution can further induce asymmetry in the development of post-shock turbulence, such as the skewed density distribution and asymmetric mixing. To understand the mechanism responsible for the skewness of the streamwise density gradient, we further examine here the orientation of the eigenvectors of strain rate tensor Si j. The PDFs of the cosines of the angles between the three eigenvectors with the streamwise direction, conditioned on regions with positive or negative density gradients, are plotted in figure 2.13. The eigenvalues of the strain rate tensor are γ1, γ2 and γ3 with the convention that γ1 < γ2 < γ3. The angles between these eigenvectors and streamwise direction are denoted by χ1, χ2 and χ3. For the multi-fluid case, in the positive density gradient regions, the extensive (γ3-) eigenvector is more likely to be aligned with the streamwise direction (figure 2.13 a). The intermediate (γ2-) 72 Figure 2.13: PDFs of the cosine angle between eigenvectors of the strain rate tensor and streamwise (x) axis for regions with (a) dρ/dx > 0 and (b) dρ/dx < 0 and for multi-fluid (solid lines) and single-fluid (dashed lines) cases. eigenvector is misaligned with the streamwise direction and the compressive (γ1-) eigenvectors have no preferential alignment. This implies that the density field is generally being stretched in the streamwise direction, making the magnitude of the density gradient smaller. On the other hand, the alignment of the γ1− and γ3− eigenvectors with the streamwise directions is reversed in the negative density gradient regions as shown in figure 2.13 (b). The density field is then compressed so that the magnitude of the density gradient is increased. This asymmetry in the alignment is induced by the VD nonlinear effects through the shock wave. This analysis serves as an explanation for the negatively skewed PDF of density gradient in the multi-fluid case. For the single-fluid case, the asymmetry in the eigenvector alignment is weak and vanishes quickly away from the shock wave. This implies that even though density variations do not affect some turbulence statistics in a straightforward manner, they modify the topology and structure of turbulence immediately after the shock and continue to manifest their effects in the post-shock turbulence evolution. 73 cosχi00.20.40.60.81PDF012345cosχ1cosχ2cosχ3cosχi00.20.40.60.81012345(a)(b) 2.3.2 Topological Analysis of the Post-shock Turbulence To further characterize the turbulence structure behind the shock wave, we have analyzed the invariant space of the VGT. The second and third invariants (denoted as Q∗ and R∗) of the anisotropic/deviatoric part of the VGT can reveal important features of the flow topology. The anisotropic part of the VGT is calculated using the formula A∗ = Ai j − θ/3I, where Ai j = ∂ui/∂xj i j is the VGT and θ = ∂ui/∂xi is the trace of the VGT of dilatation. Correspondingly, the second and third invariants can be calculated using the following formulas: Q∗ = −1 2 A∗ i j A∗ R∗ = −1 3 A∗ i j A∗ jk A∗ ki ji (2.3a) (2.3b) Similarly, the invariants of the symmetric and skew-symmetric parts of the anisotropic velocity i j, can also be calculated using the corresponding form of equation 2.3. s) and (Q∗ w) here. The following equations relate the above variables gradient tensor, S∗ They are denoted as (Q∗ for the anisotropic part of the velocity gradient tensor (Ooi et al., 1999): i j and W∗ s, R∗ w, R∗ Q∗ = Q∗ s + Q∗ w ∗ ∗ i S∗ s − ω i j ω j R∗ = R∗ (2.4a) (2.4b) w = ωi where ω∗ is the vorticity vector. The scalar variables Q∗ s and Q∗ w denote the local i dissipation rate (−Q∗ i jS∗ s = 1/2S∗ = 1/2W∗ i j) and enstrophy (Q∗ i jW∗ i j), respectively. Conse- quently, Q∗ represents the difference between enstrophy and dissipation (Chu and Lu, 2013). Similarly, R∗ ki represents the production of dissipation due to strain field and ω∗ i S∗ i j ω∗ j is the vortex stretching contribution to the enstrophy (Chu and Lu, 2013). Therefore, R∗ represents the difference between enstrophy production and dissipation production. Based on the local values of Q∗ and R∗, four types of local flow topologies can be identified: stable- focus/stretching (SFS), unstable-focus/contracting (UFC), stable-node/saddle/saddle (SN/S/S) and s = −1/3S∗ i jS∗ jk S∗ 74 unstable-node/saddle/saddle (UN/S/S). For isotropic turbulence, the joint PDF of (Q∗, R∗) has the well-known tear-drop shape. This has been further observed in other fully developed turbulent flows, such as boundary layers, mixing layers, and channel flows (Pirozzoli and Grasso, 2004; Wang et al., 2012). This type of distribution of Q∗ and R∗ is an indicator that the turbulence is more likely having a local topology of stable-focus/stretching or an unstable-node/saddle/saddle. In figure 2.14 (a), it is shown that the joint PDF of normalized second and third invariants, Q∗/(cid:104)Qw(cid:105) and R∗/(cid:104)Qw(cid:105)3/2, has the same tear-drop shape in the pre-shock flow. Using shock-LIA and DNS data, Ryu and Livescu (2014) showed that for single-fluid STI, the (Q∗, R∗) distribution is signifi- cantly modified by the shock wave, with a tendency towards symmetrization of the joint PDF. This indicates that the regions with stable-focus/compression and stable-node/saddle/saddle (first and third quadrant) are more likely to occur as turbulence develops a 2-D axisymmetric flow structure. To understand the variable density effects on this shock-induced symmetrization, the joint PDFs of (Q∗,R∗) for both single-fluid and multi-fluid post-shock turbulence are compared in figure 2.14 (b,c). Figure 2.14 (b) shows the joint distribution for the single-fluid post-shock turbulence. The dashed lines denote the locus of zero discriminant of A∗, where Q∗ and R∗ satisfy 27R∗2/4+Q∗3 = 0. Compared to the pre-shock joint PDF, there is a tendency towards symmetrization, with more points located in the first and third quadrants. This agrees very well with previous shock-LIA results (Ryu and Livescu, 2014). Similar to single-fluid STI, multi-fluid STI demonstrates a tendency towards symmetrization of the (Q∗, R∗) distribution. However, the multi-fluid distribution is slightly more symmetric and has a larger variance, with more points away from the axes. This implies that more extreme "events" exist in the post-shock multi-fluid turbulence. It also agrees with our previous study on multi-fluid STI that shows amplification of TKE to be stronger when there are significant density variations (Tian et al., 2017d,a). The density effects on the post-shock joint PDF of second and third invariants are further explored by comparing the conditional distribution, conditioned on regions with different densities, in figure 2.15 (a)-(c). Figure 2.15 (a) corresponds to regions with relatively high density (ρ > 75 Figure 2.14: Iso-contour lines of joint PDFs of normalized second and third invariants of the anisotropic part of the velocity gradient tensor, (Q∗,R∗), for (a) pre-shock, (b) single-fluid post-shock turbulence, and (c) multi-fluid post-shock turbulence. The lateral lines denote the locus of zero discriminant. (ρ + 90%ρ(cid:48) rms)), 2.15 (b) to regions with density around the post-shock mean value, and 2.15 (c) to low density regions (ρ < (ρ − 90%ρ(cid:48) rms)). For consistency check, the joint PDFs corresponding to these regions are also computed for the pre-shock flow (not shown) and found to be close to the single-fluid PDFs. After the shock wave, the joint PDFs demonstrate significant differences between regions with different densities. In regions with density closer to that of post-shock mean density, the distribution of invariants appears to be very similar to that shown in figure 2.14 (c). But for regions with higher density (figure 2.15 (a)), the joint PDF becomes more symmetric compared to the overall flow or single-fluid case. There is a much larger portion of data points having a local topology of stable-node/saddle/saddle, and fewer data points fall into the first and second quadrants, indicating larger strain-dominated regions. On the other hand, the post-shock regions with low- density values (figure 2.15 (c)) exhibit features similar to that of isotropic turbulence, with almost the same tear-drop shape, only with a larger variance or a wider distribution. The quantitative difference is hypothesized to be related to the higher shock strength variation in the multi-fluid case. It was observed in our previous studies (Tian et al., 2017b), that the local shock strength is 76 R∗-10010Q∗-10-50510R∗-4-2024-8-6-4-202468R∗-4-2024-8-6-4-202468(a)(b)(c) Figure 2.15: Iso-contour lines of post-shock (k0x ≈ 0.44) joint PDF of second and third invariants of the anisotropic part of the velocity gradient tensor, (Q∗,R∗), in regions with different densities. (a) regions with high density values, ρ > (ρ + 90%ρ(cid:48) post-shock mean value, and (c) regions with low density values, ρ < (ρ − 90%ρ(cid:48) rms), (b) regions with density around the rms). positively correlated with the pre-shock density. With a stronger shock, the two-dimensionalization effect on the post-shock turbulence should also appear stronger in the high-density regions (Livescu and Ryu, 2016). For low-density regions, the smaller two-dimensionalization effect reduces the symmetrization trend. Moreover, the relatively lower inertia in these regions leads to a faster response to the local strain field (Livescu et al., 2010), which could make a faster return to isotropic turbulence. The different characteristics of (Q∗, R∗) joint PDF in regions with different densities provide additional evidence for the previous argument made about the density role on the preferential amplification of the strain and rotation tensors. In figure 2.16, the planar distribution of the flow topologies are shown. Here, Qi refers to the quadrants on the joint PDF of (Q∗,R∗), which amounts to a representation of the local flow topology. Figure 2.16 (a) presents the 2D visualization of the flow topology in a typical x − z plane. The regions occupied by different quadrants are marked using different colors. Evidently, the vorticity-dominated regions (Q1 and Q2) cover a large portion of the flow and have more compact shapes. These regions are connected by UN/S/S areas (Q4), which are more elongated. 77 R∗-202Q∗-4-3-2-101234R∗-202-4-3-2-101234R∗-202-4-3-2-101234(a)(b)(c) Figure 2.16: Color illustration of the flow topology for the multi-fluid STI. The flow topology is represented by the quadrants (denoted as Qi) of the joint PDF of (Q∗,R∗). (a) 2D color contours in the x-z plane at y = 3.14 (half y-domain). The shock wave is located in the middle of the domain at k0x ≈ 0.0. The ratio of the fluid volume in different quadrants in the pre-shock region is Q1 : Q2 : Q3 : Q4=26.7%:38.7%:7.8%:26.8%. The 2D color contours in the homogeneous y-z plane at streamwise locations of: (b) k0x ≈ 0.2 (28.7%:34.4%:14.3%:22.6%), (c) k0x ≈ 2.0, peak TKE location (26.7%:36.9%:11.2%:25.2%) (d) k0x ≈ 4.0 (26.3%:37.9%:9.3%:26.2%). The SN/S/S (Q3) areas can be located either at the edge of Q4 regions or in-between different Q4 regions. These strain-dominated regions seem to be more fragmented than the compact vorticity- dominated regions. Figure 2.16 (b)-(d) show the 2D contours of Qi in the homogeneous (y − z) planes at different streamwise locations after the shock. These locations are also marked on figure 2.16 (a). Immediately after passing through the shock wave, the volume fractions of different quadrants are calculated to be Q1 : Q2 : Q3 : Q4 = 28.7% : 34.4% : 14.3% : 22.6%, indicating a trend towards symmetrization in the joint PDF. This can also be observed in the 2D contours in figure 2.16 (b). Moreover, the characteristic length scales associated with the regions occupied by 78 different quadrants are decreased across the shock wave. As the flow evolves away from the shock wave, the distribution slowly changes back to the pre-shock shape but still with smaller turbulence length scales. Most of fluid in different quadrants return to the pre-shock values at k0x = 4.0. The re-orientation of the flow structures into the streamwise direction is also noted in figure 2.16 (a), consistent with the return to isotropy of the flow. However, the rates at which different flow features return to an isotropic state are slightly different. The dynamics of flow and the return-to-isotropic turbulence process will be examined in detail in the next section using the Lagrangian statistics. In our previous study (Tian et al., 2017a), the (2-D) axisymmetric state of multi-fluid post-shock turbulence is explored by considering the normalized turbulent stretching term in the enstrophy budget equation. It was shown that this term is close to zero immediately after the shock and increases after the shock, more so in the multi-fluid case. In the current paper, we have shown that the multi-fluid joint PDF of (Q∗, R∗) becomes more symmetric. Below, we continue to explore how the vortex stretching rate is affected by the change in turbulence to an axisymmetric state and different flow topologies. w (Ooi et al., 1999). In figure 2.17, the joint PDF of (−Q∗ The rate at which the vorticity is stretched or contracted, i.e. the normalized vortex stretching rate, can be calculated based on the VGT invariants using the formula: Σ∗ = wiSi j wj = (R∗ s − R∗)/Q∗ s,Σ∗) is plotted to investigate the effects of the strain field on the vortex stretching rate. Positive and negative Σ∗ values correspond to the vortex being stretched or contracted. Figure 2.17 (a) shows the joint PDF of (−Q∗ s,Σ∗) for the isotropic turbulence. The results agree very well with those of Ooi et al. (1999), which indicates that the flow favors positive Σ∗ values or an overall vortex stretching, especially in the strong strain dominated regions. Here, we compare the results from isotropic turbulence to those from single- fluid and multi-fluid post-shock turbulence to understand the shock and variable density effects. We note that in figure 2.17 (b), the joint PDF becomes more symmetric around Σ∗ = 0.0 after passing through the shock wave. For the multi-fluid case, as shown in figure 2.17 (c), the joint PDF becomes almost fully symmetric, especially at lower −Q∗ s values. This symmetry has a strong effect on the overall vortex stretching rate for the multi-fluid post-shock turbulence because the positive 79 Figure 2.17: Iso-contour lines of joint PDF of (-Q∗ single-fluid and multi-fluid turbulence at post-shock position of k0x ≈ 0.44. s , Σ∗) for (a) isotropic box turbulence and (b,c) and negative Σ∗ values tend to cancel each other through averaging. Moreover, the variances of the stretching term are almost the same for single and multi-fluid cases, meaning that the lower stretching rate is mainly due to changes in the turbulence structure (especially more negative Σ∗ regions), and not simply the decrease in the magnitude of Σ∗. To understand the contribution from different topological states to the vortex stretching, the joint PDF of (−Q∗ s,Σ∗) is conditioned on different quadrants for the multi-fluid case. Figure 2.18 (a,b) shows the joint PDF for Q1 and Q2 regions, or areas with a local topology of SFS and UFC. It can be observed that in these rotation-dominated regions, the magnitude of the vortex stretching rate Σ∗ is relatively small. Moreover, Q1 is dominated by the negative Σ∗ and Q2 is dominated by positive Σ∗ areas, which can be inferred from their corresponding topologies. However, for Q3 and Q4 (figure 2.18 c,d), the magnitude of the vortex stretching rate is larger than that in the rotation- dominated regions (figure 2.18 a,b). In Q3, the joint PDF is relatively symmetric and seems to be slightly biased towards negative vortex stretching rates. Q4, on the other hand, is dominated by positive vortex stretching. Overall, the results explain the lower averaged vortex stretching rate in the multi-fluid case through enhancement of Q1 and Q3 regions. 80 Σ∗-202−Q∗s05101520Σ∗-20205101520Σ∗-20205101520(a)(b)(c) Figure 2.18: Iso-contour lines of joint PDF of (-Q∗ shock wave. (a) Q2, (b) Q1 ,(c) Q3 and (d) Q4. s , Σ∗) for different quadrants right after the 2.3.3 Lagrangian Dynamics of VGT The Lagrangian particle method has been widely used for studying many important aspects of turbulence, such as fluid particle dispersion, clustering, and mixing. Lagrangian analysis of the VGT is also proven beneficial in the study of turbulence (e.g. Chong et al., 1998; Ooi et al., 1999; Meneveau, 2011; Wang and Lu, 2012). Even though the turbulence structure at different post-shock locations can be extracted from the Eulerian data, the dynamics that dominate the development of VGT are still missing. In this section, we use non-inertial Lagrangian particles/tracers that move with the local flow velocity because their statistics can reflect important transient turbulent dynamics, which are difficult to study using Eulerian data (Yeung, 2002). More importantly, in the context of variable density flows, the Lagrangian statistics enable us to differentiate among particles 81 -3-2-10123−Q∗s051015-3-2-10123051015Σ∗-3-2-10123−Q∗s051015Σ∗-3-2-10123051015(a)(b)(c)(d) with different densities and investigate their dynamics separately. Lagrangian data are used here to perform a detailed analysis of the post-shock turbulence structure and its dynamics in the STI with and without significant density fluctuations. The first result considered is the timescale of particle motions related to different flow topologies. In figure 2.19, the percentage of fluid particles that remain in their starting quadrants are plotted over time so that we can identify the residence time of particles for different turbulence structures. In figure 2.19 (a), the percentages of fluid particles are plotted for decaying isotropic turbulence as a reference. It is noted that Q3 and Q4, which are the strain-dominated regions, have the smallest associated residence times among all the four quadrants, with Q3 time being the smaller of the two. For the rotation-dominated regions, the residence times are expectedly longer, especially for Q2. The residence times for single-fluid and multi-fluid simulations can be inferred from figure 2.19 (b,c). For both cases, Q3 always has the least residence time and Q2 has the largest one. Comparing all three cases, the particles in the multi-fluid and single-fluid post-shock turbulence are shown to evolve faster away from the original quadrant than the particles in isotropic turbulence. This is hypothesized to be due to the fact that post-shock turbulence still experiences a return-to-isotropy, while the isotropic turbulence has already reached an equilibrium among all the quadrants. Between the two post-shock turbulence fields, the multi-fluid case presents a relatively shorter residence time. This agrees with the fact that the multi-fluid turbulence exhibits a faster return-to-isotropy than the single-fluid turbulence (Tian et al., 2017a). The rotation-dominated regions are sometimes identified as coherent structures like stretched vortex tubes, as shown in figure 2.20 (a). When the particles travel through these structures, they tend to stay longer due to the coherency of the surrounding flow. This is confirmed by the circulating streamlines in such structures (black lines in figure 2.20 (a)). On the other hand, the strain-dominated regions are usually lacking such coherency, making the residence times smaller. The local streamlines pass through such structures instead of circulating around them. Figure 2.21 presents an example of the temporal development of the above-mentioned structures. The evolution of a vortex tube in the post-shock turbulence is tracked and visualized as it moves away from the 82 Figure 2.19: Percentage of fluid particles that stay in each quadrant following particles initialized uniformly in (a) isotropic turbulence and (b,c) single-fluid and multi-fluid turbulence at post-shock position of k0x = 0.44. shock wave in figure 2.21 (a). As expected, the vortical structure maintains its shape, except that it is being stretched and reoriented by the local flow field. Moreover, the vortex tube surface is almost parallel with the iso-surface of the density field, i.e. perpendicular to the density gradient. This is consistent with the discussion in section 2.3.1.1 regarding the bulk of vorticity generation across the shock wave. As the vortex tube evolves away from the shock wave, the reorientation of the density gradient by the vortex is also observed. In figure 2.21 (b), a strain-dominated structure is visualized using the iso-surface of negative Q∗. It can be clearly seen that such structures lack temporal coherency since they tend to be become fragmented as they evolve. In figure 2.22, the contributions to the normalized vortex stretching rate from particles that are initialized in each of the four quadrants are plotted following these particles. As expected, at t = 0.0, particles from Q2 and Q4 add positively to the vortex stretching rate, while those from Q1 have a negative vortex stretching rate contribution on average. This is in good agreement with the joint PDFs of (−Q∗ s,Σ∗), shown in figure 2.18. For Q3, the initial contribution is close to zero. As the fluid particles move with the turbulent flow, their contributions to the vortex stretching also change. It can be seen that the fast increase in vortex stretching can be mainly attributed to the 83 t00.5100.20.40.60.81Q1Q2Q3Q4t00.5100.20.40.60.81t00.5100.20.40.60.81(a)(b)(c) Figure 2.20: Visualization of turbulence structure using iso-surfaces of Q∗ colored by R∗ in the multi-fluid post-shock turbulence. (a) Q1 and Q2 (vorticity-dominated region), and (b) Q3 and Q4 (strain-dominated region). The black lines represent the instantaneous streamlines. Figure 2.21: Visualization of the temporal development (left to right) of the turbulence structure using iso-surfaces of Q∗ colored by density for the multi-fluid post-shock turbulence. These structures are captured immediately after the shock wave. (a) vorticity-dominated structure, and (b) strain-dominated structure. 84 Figure 2.22: Contributions to the vortex stretching rate from particles starting in each quadrant. The particles are initialized uniformly at the post-shock position k0x ≈ 0.44 and traced downstream till the vorticity returns to an isotropic state. (a) single-fluid and (b) multi-fluid cases. particles originating in Q1 and Q3. Particles starting in Q4 have an increasing vortex stretching contribution for a short period before their combined/average contributed value starts to decrease. The behavior is qualitatively similar for the single-fluid case, but the changes are smaller in this case. For both cases, the vortex stretching contribution from the initial Q2 particles decreases in time. To further understand this behavior, the Lagrangian equations of the VGT and its invariants are considered. The time evolution of Ai j for fluid particles can be obtained by taking the spatial derivatives of the Navier-Stokes equations. In dimensionless form, it can be written as (Chu and Lu, 2013): ∂ Ai j ∂t + uk ∂ Ai j ∂xk DAi j Dt + Aik Ak j = −Hi j + Ti j = −Aik Ak j − Hi j + Ti j (2.5a) (2.5b) with 85 ) = − 1 ρ2 ) ∂p ∂xi ∂σik ∂xk Hi j = Ti j = σi j = ∂ ∂xj ∂ ∂xj µ Re0 (1 ρ (1 ρ (cid:18) ∂ui ∂xj + ∂uj ∂xi − 2 3 ∂uk ∂xk δi j ∂ ρ ∂xj ∂p ∂xi 1 ρ ∂p2 ∂xi∂xj = Hb i j + Hp i j + (cid:19) (2.6a) (2.6b) (2.6c) where Re0 is the reference Reynolds number. From here, the dynamic equations for the three invariants of the VGT, P, Q, and R can be derived in the following form (Chu and Lu, 2013): DP Dt DQ Dt DR Dt + Hb ii − Tii + Ai j Hp = (P2 − 2Q) + Hp ii = (PQ − 3R) + (PHp ii = PR + (QHp ii +Ai j Ajk Hb ji) + (PHb + Ai j Ajk Hp ii + Ai j Hb ki) + (QHb + PAi j Hp ji ki) + (−QTii − PAi jTji − Ai j AjkTki) ji) + (−PTii − Ai jTji) ii + PAi j Hb ji where the three invariants of VGT are defined as: P = −tr(Ai j) Q = −1 2 R = −det(Ai j) (tr(Ai j)2 − tr(Ai j Ajk)) (2.7a) (2.7b) (2.7c) (2.8a) (2.8b) (2.8c) Here, tr(Ai j) and det(Ai j) denote the trace and determinant of a tensor. Note that instead of the deviatoric part of the VGT, the dynamic equations for the full VGT are considered. The reason is that due to the variable density effects and shock compression, the incompressibility condition is not satisfied especially when Mt and At become large. More importantly, the contributions from the isotropic part of the VGT to the flow dynamics are still unknown and need to be explored for STI. 86 Figure 2.23: PDFs of (a) (DQ/Dt)/(cid:104)Qw(cid:105)3/2 and (b) (DR/Dt)/(cid:104)Qw(cid:105)2 for fluid particles with different densities at streamwise location of k0x ≈ 0.5. The dynamical equations can be divided into contributions by four different parts: I) mutual- i j and IV) viscous term interaction among invariants, II) pressure Hessian, Hp Ti j. The statistics regarding these terms can be extracted from the Lagrangian data. i j, III) baroclinic, Hb Some general features of the Lagrangian dynamics of the VGT invariants are examined through the PDFs of their material derivatives. The variable density effects can be identified by comparing the PDFs corresponding to regions with different densities (figure 2.23). In the light fluid regions, the PDFs of DQ/Dt and DR/Dt have narrower tails, while the tails are larger in the heavy fluid regions. Another important observation is that the skewness of DQ/Dt is different in the light and heavy fluid regions. Heavy fluid particles have a positively-skewed PDF, similar to the overall flow. On the other hand, the DQ/Dt skewness resulting from light fluid particles is negative. This implies that heavy fluid particles are more likely to move towards rotation-dominated regions and vice versa. These differences can be attributed to different return-to-isotropy, experienced by fluid particles with different densities. The Lagrangian dynamics of the turbulence and the evolution of flow topology are further examined here by considering the conditional mean rate of change of Q and R in the invariants 87 (DQ/Dt)/hQwi3/2-20-100102010-410-310-210-1100(DR/Dt)/hQwi2-20-100102010-410-310-210-1100heavyfluidlightfluidmean(a)(b) Figure 2.24: Conditional mean rate of change vectors of (DQ/Dt/(cid:104)Qw(cid:105)3/2,DR/Dt/(cid:104)Qw(cid:105)2) in the (Q, R) plane for (a) isotropic turbulence, (b) single-fluid post-shock turbulence, and (c) multi-fluid post-shock turbulence at streamwise location of k0x ≈ 0.5. plane (Ooi et al., 1999). The rates of change are used to form a vector at each point in the invariants plane. The trajectories implied by these vectors can be followed to understand the return-to-isotropy. The procedure used to obtain the conditional mean vectors (CMVs) in this study is similar to that in Ooi et al. (1999). Based on the conditional averages introduced in equation 2.2, X(Q, R) represents a statistical quantity that is conditioned on Q and R. The statistical convergence concerning the bin sizes and the number of samples in each bin has been discussed in section 2.2.4. Figure 2.24 shows the normalized conditional mean vectors (DQ/Dt/(cid:104)Qw(cid:105)3/2,DR/Dt/(cid:104)Qw(cid:105)2) for different flows. The vectors obtained from isotropic turbulence data are shown in figure 2.24 (a) for reference. It can be seen that the CMVs exhibit a circulating behavior in the (Q, R) plot around the origin in the clockwise direction, indicating that the flow evolves from SFS to UFC, UN/S/S, SN/S/S then back to SFS on average. This circulating behavior represents the Lagrangian dynamics in fully developed turbulence that maintains the tear-drop shape of the (Q, R) distribution. This has been observed in many incompressible/compressible canonical turbulent flows (Ooi et al., 1999; Chu and Lu, 2013). The CMVs for single-fluid and multi-fluid post-shock turbulence are shown in figure 2.24 (b) and (c). Evidently, the joint PDF of (Q, R) becomes more symmetric due to shock compression. From the Lagrangian point of view, the circulating behavior as seen in figure 2.24 (a) isotropic turbulence disappears. The particles in Q2 tend to have an increasing Q and decreasing R, resulting in an overall trend of getting away from the original point, instead of circulating and 88 R/hQwi3/2-1-0.500.51Q/hQwi-2-1.5-1-0.500.511.52R/hQwi3/2-1-0.500.51-2-1012R/hQwi3/2-0.500.5-1-0.500.51(a)(b)(c) Figure 2.25: Conditional mean vectors in the (Q, R) invariants plane for (a) light fluid, (b) medium density fluid and (c) heavy fluid at streamwise location of k0x ≈ 0.5. then moving toward Q1. This trend in the second quadrant represents an increase of enstrophy. The particles in Q1 have similar dynamics as in isotropic turbulence and tend to move downward in the (Q, R) plane toward the zero discriminant curve. The particles in Q3 are more likely to move straight up towards Q2, while those in Q4 are likely to move away from the original point following the direction of the zero discriminant line. The overall behavior formed by these particles demonstrates the return-to-isotropy process, with an enlarging head in the second quadrant and elongating tail in the fourth quadrant, anticipating the formation of the classic tear-drop shape. The density effects can be further examined by conditioning the (DQ/Dt, DR/Dt) vector field on the local density. Figure 2.25 (a) shows the CMVs for the light fluid regions. The light fluid particles retain the circulating motion, except that the particles in Q3 and Q4 are likely to go straight left instead of following the zero discriminant line. In general, the flow dynamics in the light fluid regions are less affected by the shock wave. For the medium density fluid regions (figure 2.25 b), the circulating motion disappears. On the right side of the (Q, R) plane (R > 0), which is the strong dissipation-production region based on equation 2.4, the fluid particles tend to move downward, resulting in lower Q values. On the left side of the (Q, R) plane (R < 0), which is the enstrophy-production dominated region, the fluid particles tend to move to the left, indicating an increased enstrophy-production. The overall downward-moving behavior of the medium density 89 R/hQwi3/2-0.200.2Q/hQwi-1-0.500.51R/hQwi3/2-0.200.2-1-0.500.51R/hQwi3/2-0.200.2-1-0.500.51(a)(b)(c) Figure 2.26: Contributions to the transport equations of the VGT invariants by different terms for isotropic turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term, and (d) viscous term. fluid particles is indicative of decreasing vorticity. This is possibly due to the fact that vorticity is preferentially amplified in the medium density region across the shock wave. After passing the shock wave, the vorticity will decrease as the correlation between density and vorticity vanishes. Figure 2.25 (c) shows the CMVs for the heavy fluid regions. Interestingly, the heavy fluid particles exhibit counterclockwise motion. The heavy particles start from Q3 and move to Q4, Q1, and finally to Q2. This implies that they become vorticity dominated due to the fast depletion of strain. Evidently, density plays an important role in the development of the flow topology in the post-shock region, so special attention should be paid to the modeling of variable density STI. To better understand the underlying mechanisms that cause the behavior highlighted above, 90 -1-0.500.51Q/hQwi-2-1012-1-0.500.51-2-1012R/hQwi3/2-1-0.500.51Q/hQwi-2-1012R/hQwi3/2-1-0.500.51-2-1012(a)(b)(c)(d) the dynamic equations (2.7) governing the vector (DQ/Dt, DR/Dt) are examined. The right hand sides of these equations can be divided into four types of contributing terms: (a) mutual interaction among invariants (term I), (b) pressure Hessian (term II), (c) baroclinic (term III), and (d) viscous term (term IV). As a reference, these terms are shown in figure 2.26 for the isotropic turbulence. The variable Q tends to be amplified in the enstrophy-production dominated region due to the effects of vortex stretching mechanism and is decreased in the dissipation-production dominated region due to self-amplification of the strain rate tensor. On the other hand, the mutual effects on R are small because the first invariant P is usually small and the positive and negative values are likely to cancel each other. The contributions from the pressure Hessian (figure 2.26 b) tend to move the particles away from an asymptotic line, ending up amplifying the magnitude of R. For the current simulation, the asymptotic line starts from Q2 and ends in Q4 with a slope of around -2.5. This result agrees well with that observed in turbulent boundary layers (Chu and Lu, 2013). The baroclinic contributions are very small in the post-shock turbulence as shown in figure 2.26 (c). The viscous effects as shown in figure 2.26 (d) and as expected are reducing the magnitudes of Q and R and pushing the particles towards the origin. This has been observed in various types of turbulence (Ooi et al., 1999; Chu and Lu, 2013). The combined effects from the four above mechanisms determine the circulating behavior of the conditional mean of (DQ/Dt, DR/Dt) vectors. After interaction with the shock wave, the conditional mean vectors in the (Q, R) plane are different from those in the pre-shock isotropic turbulence. Figure 2.27 shows the results for the single-fluid case. By comparing it with figure 2.26, we note that even though the conditional means of (DQ/Dt, DR/Dt) vectors are different, the contributions from mutual interaction (figure 2.27 a), baroclinic term (figure 2.27 c) and viscous term (figure 2.27 d) are very similar. The only term that qualitatively distinguishes the post-shock turbulence from isotropic turbulence is the pressure Hessian term as shown in figure 2.27 (b). In the post-shock single-fluid turbulence, the asymptotic line that separates the vectors into two regions with different behaviors disappears. Instead, the pressure Hessian term tends to move the particles away from the origin in Q1 and Q2, thus increasing Q and |R| values of the particles. In Q3 and Q4, the vectors are parallel to the left 91 Figure 2.27: Contributions to the transport equations of the VGT invariants by different terms for single-fluid post-shock turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term, and (d) viscous term. zero discriminant line, making the particles move from Q3 to Q4, and then to Q1. For multi-fluid post-shock turbulence, the pressure Hessian term is also the only term that is qualitatively different than that in isotropic turbulence (figure 2.28). Despite the increased density and pressure gradient in the multi-fluid case, the baroclinic term is still considerably smaller than all the other terms. In Q1 and Q2, an asymptotic line similar to that in isotropic turbulence seems to exist, which repels the vectors away from it, causing an increase in |R| values. In Q3 and Q4, this term tends to move the fluid particles along the right zero discriminant line from Q4 to Q3. The further conditioned pressure Hessian term based on the local densities in figure 2.29 indicates that fluid particles with different densities exhibit very different behaviors in respond to the pressure Hessian dynamics. Specifically, the pressure Hessian generally moves the heavy particles toward 92 -1-0.500.51Q/hQwi-2-1012-1-0.500.51-2-1012R/hQwi3/2-1-0.500.51Q/hQwi-2-1012R/hQwi3/2-1-0.500.51-2-1012(a)(b)(c)(d) Figure 2.28: Contributions to the dynamics of the VGT invariants by different terms for multi-fluid post-shock turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term, and (d) viscous term. the regions with larger Q values. In Q3 and Q4, it also moves the heavy fluid particles towards the R > 0 plane. For the light fluid particles, the pressure Hessian term tends to make them move towards regions with larger |R| values in the first and second quadrant. In Q3 and Q4, the fluid particles move from Q4 to Q3. Last but not the least, the fluid particles with medium density seem to exhibit similar behavior to light fluid particles, except that in Q1, where the pressure Hessian contribution is moving the fluid particles towards the regions with large Q values. Examining figure 2.25 and figure 2.29 together, we observe that the different particle dynamics in the (Q, R) plane in regions with different densities are mainly due to differences in the pressure Hessian contribution. 93 -0.4-0.200.20.4Q/hQwi-1-0.500.51-0.4-0.200.20.4-1-0.500.51R/hQwi3/2-0.4-0.200.20.4Q/hQwi-1-0.500.51R/hQwi3/2-0.4-0.200.20.4-1-0.500.51(a)(b)(c)(d) Figure 2.29: Contributions from pressure Hessian to the dynamics of the VGT invariants in (a) light fluid region, (a) medium density fluid region and (c) heavy fluid region. 2.3.4 Lagrangian Particle Dispersion 2.3.4.1 Single Particle Statistics In this section, single particle statistics such as particle dispersion rate, autocorrelation of velocity will be presented. The Lagrangian statistics of fluid particles reveal important features of turbulent flow and are considered in this study for better understanding of post-shock flow dynamics. We will first consider single particle statistics, which are closely related to the mean concentration of a passive scalar. Stochastic models for particle dispersion has been proposed to predict various statistics. Due to the existence of mean flow, it would be reasonable to examine the particle displacement i (t) = relative to the mean flow. Thus, the relative particle dispersion can be calculated using y+ i (t) − ui(t − t0). Note that the only non-zero component of the mean velocity is the streamwise x+ component. The particle dispersion rate can then be calculated using Di j(t) =(cid:10)y+ diagonal components of the dispersion rate tensor Di j(t) are plotted in figure 2.30. Due to the homogeneity in the spanwise direction, D22 is equal to D33, so only D22 is shown in the figure. On top of that, the displacement variance is normalized using t2 to better visualize the temporal development. Figure 2.30 compares the single particle dispersion rate between single- and multi- α(t)(cid:11). The α(t)y+ 94 R/hQwi3/2-0.4-0.200.20.4Q/hQwi-1-0.500.51R/hQwi3/2-0.4-0.200.20.4-1-0.500.51R/hQwi3/2-0.4-0.200.20.4-1-0.500.51(a)(b)(c) fluid post-shock turbulence. For the multi-fluid case, both the streamwise and spanwise particle dispersion follows a t2 development. Later downstream, the streamwise particle dispersion rate increases and reaches a peak around t − t0 ≈ 0.5 and then decreases to a t scaling at larger time. On the other hand, the spanwise dispersion decreases from t2 to t scaling at larger time. Till the end of the domain, the streamwise particle dispersion is always larger than the that in the spanwise direction. For the single-fluid case, similar temporal scaling for each component of the particle dispersion is observed. However, the relative dispersion rate between the streamwise and spanwise direction are different. It is noted that in the single fluid simulations, the particle dispersion is initially faster in the spanwise direction. As the particles move downstream, D11 starts to grow and eventually becomes larger than D22. This shows that the dispersion anisotropy becomes different due to the variable density effects. The difference in the early time anisotropy can be related to the anisotropy in the Reynolds stress. The early time particle dispersion can be approximated using iu(cid:48) j(t − t0)2. Therefore, the early time anisotropy in the particle dispersion is strongly correlated with the Reynolds stress and this is also supported by the Eulerian results shown in figure 2.10 Di j =(cid:10)y+ α(t)(cid:11) ≈ u(cid:48) α(t)y+ Most of previous studies on Lagrangian particle dispersion are conducted for fully homogeneous Figure 2.30: Single particle dispersion rate normalized using t2 for both single and multi-fluid post-shock turbulence. 95 10-2100Dαα/t20.020.040.060.080.10.120.14D11/t2(m)D22/t2(m)D11/t2(s)D22/t2(s) turbulence, so the third-order moment of the particle displacement PDF is zero. For post-shock turbulence, the flow is only homogeneous in the spanwise directions, so the third-order moment or skewness of the PDF is expected to be zero in those directions. In the streamwise direction, the turbulence experiences complicated transient process as it evolves downstream, so it is not known what the skewness of the particle displacement PDF is. In figure 2.31, the skewness of particle displacement in the streamwise direction is shown for both cases. We note that the skewness of the dispersion µ3 is non-zero. Immediately after the shock wave, the multi-fluid case has a negative skewness of around -0.15. Then it increases to around -0.12 and remains at the value till the end of the domain. For the single-fluid case, however, the skewness is positive immediately after the shock wave. It then increases for a short period before decreases to a negative value. Due to the fact that many stochastic models assume a Gaussian process for particle dispersion, the non-Gaussian behavior of particle dispersion in streamwise direction needs to be addressed. It is well known that Lagrangian velocity autocorrelation plays an important role in the de- velopment of the particle dispersion (Yeung, 1994), which can be calculated using the formula: Figure 2.31: The temporal development of the skewness µ3 of the streamwise particle displacement PDF. 96 t00.511.52µ3-0.2-0.100.1µ3,x′1(m)µ3,x′1(s) (cid:68) (t0)(cid:69)1/2(cid:68) i (t0)u+ u+ u+2 j (cid:68) u+2 i j (t)(cid:69) (t)(cid:69)1/2 . For any homogeneous turbulence, the autocorrelation function ri j(t0, t) = will only depend on the time delay τ = t − t0. However, the flow is not homogeneous in the streamwise direction for post-shock turbulence, so the autocorrelation will depend on both t and t0. Here, we will focus on how the particle velocity autocorrelation from immediately after the shock wave. Figure 2.32 shows the diagonal components of the Lagrangian particle velocity autocorre- lation. Again, due to the homogeneity in the spanwise direction, only r11 and r22(t) are shown. Comparing single-fluid with multi-fluid case, we can clearly see that both velocity components in the single-fluid case have a stronger correlation than the multi-fluid case. For the multi-fluid case, the velocity autocorrelation functions for both velocity components are very similar to each other. For the single-fluid case, the spanwise velocity has a much stronger correlation than the streamwise component. Results from the first chapter have shown that the density variations in the multi-fluid mixture cause a differential distribution of turbulence statistics in regions with different densities. For the Lagrangian study, it would be reasonable to differentiate heavy and light fluid particles to understand the density effects. In figure 2.33, the Lagrangian velocity autocorrelation of heavy and light fluid Figure 2.32: Lagrangian velocity autocorrelation function for particles staring immediately after the shock wave. 97 t00.511.52rαα00.20.40.60.81r11(t0,t)(m)r22(t0,t)(m)r11(t0,t)(s)r22(t0,t)(s) particles are compared with the overall average. We note that for streamwise velocity, both heavy and light fluid particles have weaker correlations than the average. For the spanwise velocity, however, the autocorrelation functions for particles with different densities are very similar. 2.3.4.2 Particle Pair Statistics In this section, we will present the particle pair statistics such as the particle dispersion rate and particle pair velocity correlation to address the effects of density variation on particle pair dispersion. As described in section 3.2, the particle pair statistics depend strongly on the initial separation direction. For simplification, we will only address the initial separation direction as streamwise (lx = l0) or spanwise ( z = lyz = l0). (cid:113) + l2 l2 y +(2) j +(1) y i The first thing to examine is the diagonal components of the relative particle dispersion tensor i j(t) = Dr with different magnitudes and directions of the initial separation vectors. Figure 2.34 shows that the pair particle dispersion rate for particles that are initialized in the spanwise direction with initial separation distance of l0 = 0.25η, 4η and 64η. The results in y (cid:68) (cid:69) Figure 2.33: Lagrangian velocity autocorrelation function for particles staring immediately after the shock wave in the multi-fluid STI simulations. The heavy fluid particles (dashed line) and light fluid particles (dotted line) are shown for comparison. 98 00.511.52rαα(t0,t)00.20.40.60.81r11(t0,t)(m)r22(t0,t)(m) αα ≈(cid:68)(u+(1) α α − u+(2) )2(cid:69) (t − t0)2. On top of that, multi-fluid case (red lines) shows an figure 2.34 are normalized with t2 to better visualize the temporal scaling. At early times with t − t0 < tη, the particle dispersion for all the initial separation distances follows the t2 scaling. This is expected due to the fact that the early time dispersion is determined by the relative velocity variance, Dr overall larger dispersion rate than the single-fluid cases with the same initial separation distance. At intermediate times, interesting transient processes are observed for both single-fluid and multi-fluid cases. The anisotropy in the particle dispersion is mainly manifested in this transient process. For both single- and multi-fluid cases, the particle dispersion in three directions starts to have different behaviors. After this transient adjustment, the particle dispersion in three direction starts to converge to similar temporal scaling. In figure 2.34 (b), the initial particle separation distance is 4η. At intermediate times, different transient processes are also shown for the particle dispersion in three directions. Some similarity in the transient process between figure 2.34 (a) and (b) are observed. This is possibly due to the fact that the Reynolds number considered in this study is not Figure 2.34: Evolution of the pair particle dispersion rate normalized by t2 for both single-fluid (black) and multi-fluid simulations (red). The direction of the initial separation is in the spanwise direction with lyz = l0. The relative particle dispersion in three direction are separated into streamwise direction (dotted lines), longitudinal direction (solid line) and transverse direction (dashed lines). The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. 99 large (Reλ ≈ 40). For large initial separation distance (figure 2.34 (c)), we note that the temporal scaling of relative particle dispersion is very similar to the single-particle dispersion results (figure 2.30). This indicates that with large initial separation, the correlation between the two particles is too small, so the two-particle dispersion is very similar to that of single-particle dispersion. 22 ≈ Dr Figure 2.35 shows the same relative particle dispersion statistics for initial particle separation in the streamwise direction with lx = l0. Due to the homogeneity in the spanwise directions, we note that for all the cases in figure 2.35, Dr 33. At both early and large times, the temporal scaling agrees with previous observations. At the intermediate time, the anisotropy due to the different initialization direction starts to play a role in the relative particle dispersion. For both lx = l0 = 0.25η and 4η cases, we observe that the Dr 33) goes through a different transient process than Dr 11 during intermediate time. Through this transient process, the streamwise particle dispersion becomes larger than in the spanwise direction before they reach large time temporal scaling. For large initial separation distance, however, the transient process is less pronounced and the temporal scaling goes from t2 to t. For all the cases, it seems that the density variations only affect the magnitude of the particle dispersion rate and anisotropy, but no effects on the temporal scaling. 22 (Dr The Lagrangian particle relative velocity correlation is an important statistics that plays a central role in determining the relative particle dispersion (Salazar and Collins, 2009). Figure 2.36 shows relative velocity correlation for initial separation in the spanwise direction (lyz = l0) with initial separation distance l0 equals to 0.25η, 4.0η and 64η. For small initial separation, the velocities of the particle pair remain highly correlated till the end of the domain. Comparing each of the velocity components, the correlation is stronger for the streamwise velocity. For initial separation of l0 = 4η, the particle pair velocity correlation for all components are strong. As the particle evolve downstream, the particle pair velocities slowly de-correlate, but this process is even slower for the streamwise component. For large initial separation distance, the particle pair velocity correlation is very low. The streamwise velocity undergoes a short-lived transient process before slowing growing back. For the other velocity components, the correlation keeps decreasing after the very 100 short transient processes. Comparing single-fluid with the multi-fluid case, it seems that the density variations only affects the magnitude of velocity correlation. For initial separation in the spanwise direction (lx = l0), the particle pair velocity correlation is shown in figure 2.37. With small initial separation distance, the behaviors of the velocity correlation coefficient are very similar to those observed in figure 2.36, except the faster decorrelation. At intermediate initial separation, all velocity components experience transient processes of increasing and then decreasing correlation. On top of that, the length scales of the transient processes are larger for the streamwise velocity component. For large initial separation distance, all the correlation coefficients increase as the particles evolve away from the shock wave. The implies that at this scale, the separation distance is larger than the transient process length scale, so the transient process won’t affect the particle velocity correlation. 2.4 Conclusions Accurate shock-capturing turbulence-resolving simulations are conducted together with Eu- lerian and Lagrangian particle tracking post-processing methods to investigate the interaction of isotropic turbulence with a normal shock wave for both single-fluid and a mixture of different density fluids. The main objective is to develop a better understanding of the variable density and shock effects on the post-shock turbulence structure. Grid and LIA convergence tests are conducted to establish the numerical accuracy of the simulated data. The results show that the turbulence statistics are grid-converged and the LIA prediction can be approached, indicating good accuracy of the current computational method. Statistical convergence is also tested for both Eulerian and Lagrangian statistics. The density effects on the post-shock turbulence structure are first studied using Eulerian data. The non-Gaussianity of the velocity gradient tensor (VGT) is studied by examining the PDFs of velocity gradient components. The preferential amplification of rotation and strain rate tensors is found to be affected by the density variations, leading to a weaker correlation between the two tensors in the multi-fluid case. This is shown to be caused by different roles that density plays 101 on the modification of rotation and strain rate tensors across the shock wave. The skewness and flatness of VGT components before and after the shock wave are then examined to understand the evolution of VGT. It is shown that density effects are weak across the shock, but are stronger in the post-shock development. The density variations are also shown to cause the preferential alignment between eigenvectors of the strain rate tensor and density gradient vector, which then modifies the skewness of the velocity gradient and density gradient PDFs. The density effects on the turbulence topology are then examined by first comparing the joint PDF of the second and third invariants of the velocity gradient tensor. The pre-shock joint PDF has the classic tear-drop shape. However, after the shock wave, a tendency towards symmetrization of the joint PDF, as in single-fluid STI, is observed for the multi-fluid case, with more data points falling into the first and third quadrants. After conditioning the joint PDFs based on fluid density, large differences among heavy, medium, and light fluid regions are observed. In the heavy fluid regions, the joint PDF becomes almost completely symmetrical with an increasing portion of data in the third quadrant. In contrast, the majority of the light fluid data points have a similar distribution to that of isotropic turbulence. A connection between low vortex stretching and the joint Q∗,R∗ statistics is established for the post-shock turbulence, by considering the contribution to vortex stretching rate from each quadrant. Furthermore, Lagrangian fluid particles are used to track the development of the turbulence and VGT after the interaction with the shock. The Lagrangian dynamics of the VGT are also examined by using the conditional mean rate of change of the invariants of VGT. The results show that particles in Q3 have the least residence time, while those in Q2 have the longest residence times. Due to the return-to-isotropy experienced in the post-shock turbulence, the residence times are smaller than those in isotropic turbulence, especially in the multi-fluid case. It is also shown that particles starting in quadrants Q1 and Q3 play an important role in recovering of the vortex stretching rate. After interacting with the shock wave, the clockwise circulating behavior as observed in the isotropic turbulence disappears in both single- and multi-fluid cases. Our analysis highlights the mechanism through which post-shock turbulence recovers the classical tear-drop 102 shape, with an enlarging head in the second quadrant and elongating tail in the fourth quadrant. The contributions from different terms in the dynamic equations of VGT invariants compared with isotropic turbulence show that the pressure Hessian term is critical to the topological evolution of turbulence. The pressure Hessian term is also shown to be strongly dependent on the local density in the multi-fluid case, resulting in completely different dynamics in regions with different densities. Lagrangian statistics of non-inertial particles in post-shock turbulence show that the single-particle dispersion rate and anisotropy can be correlated to the development of post-shock Reynolds stress. The particle dispersion in the streamwise direction is found to be non-Gaussian, with the skewness of the dispersion PDF dependent on the density variations. Lagrangian integral time scales of fluid particles with different densities are obtained from velocity autocorrelation functions. Particle pair dispersion generally follows the temporal scaling for isotropic turbulence. 103 Figure 2.35: Evolution of the pair particle dispersion rate normalized by t2 for both single-fluid (black) and multi-fluid simulations (red). The relative particle dispersion in three direction are separated into streamwise direction (dotted lines), longitudinal direction (solid line) and transverse direction (dashed lines). The direction of the initial separation is in the spanwise direction with lx = l0. The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. Figure 2.36: Evolution of the pair particle velocity correlation for both single-fluid (black) and multi-fluid simulations (red). Velocities in the streamwise direction (solid lines), longitudinal direction (dashed line) and transverse direction (dotted line) are compared separately. The direction of the initial separation is in the spanwise direction with lyz = l0. The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. 104 Figure 2.37: Evolution of the pair particle velocity correlation for both single-fluid (black) and multi-fluid simulations (red). The direction of the initial separation is in the spanwise direction with lx = l0. Velocities in the streamwise direction (solid lines) and transverse directions (dotted and dashed lines) are compared separately. The corresponding initial separation distance is (a) l0 = 0.25η, (b) l0 = 4η and (c) l0 = 64η. 105 CHAPTER 3 EFFECTS OF DENSITY VARIATIONS ON A NORMAL SHOCK WAVE 3.1 Introduction In previous chapters, we studied the effects of shock wave on variable density flow. In this chapter, we will focus on the shock movement in a variable density environment. Extreme casses involving complex shock surface topology and shock breakup will not be considered. When a normal shock wave propagates through a non-uniform media, the strength and geometry of the shock wave, as well as the media itself, will be modified through their interactions. This is an important fundamental problem in the study of Inertial Confinement Fusion (ICF), supersonic and hypersonic combustion, astrophysics, and other applications (Ranjan et al., 2011; Hesselink and Sturtevant, 1988; Bisnovatyi-Kogan and Silich, 1995). In an early theoretical study, Chisnell (1955) investigated the one-dimensional propagation of shock waves through a non-uniform (variable density) media and derived an expression for the shock Mach number as a function of time. Later, Whitham (1958) developed a more general analytical solution for a medium with both varying density and pressure. These earlier formulations of shock propagation, even though useful, suffer from a number of difficulties, one of which originates from the fact that the waves generated from the post-shock region, i.e. re-reflected waves are neglected. The re-reflected waves effects are accounted for in the study conducted by Bird (1961), who performed numerical simulations with the method of characteristics and compared them with the Chisnell-Whitham (as referred to in Bird (1961)) results and asymptotic solutions. Soukhomlinov et al. (2002) further examined the effects of temperature gradients on shock propagation. Their results were in good agreement with Bird’s results. The above approaches can be applied to non-uniform flows with large density fluctuations. In the limit of weak fluctuations, this problem can be linearized to obtain the LIA prediction. To address the effects of pre-shock non-uniform density field, McKenzie and Westphal (1968) studied the interaction of isobaric entropy waves with an oblique shock. Mahesh et al. (1997) then treated 106 the 3D density fluctuations as a combination of vorticity and entropy waves and conducted a linear analysis. Griffond (2005); Huete et al. (2011) has considered the pre-shock non-uniform density field as a mixture of gases with different molecular weights and extended the linear analysis to multi- component flows. These LIA results are capable of predicting both unsteady shock propagation and post-shock fluctuations, but most of the discussions have been focused on the later. Numerical examination on this fundamental problem has been focused on the shock interaction with turbulence or Shock-Turbulence Interaction (STI). The canonical STI is an example of shock propagation in non-uniform media. Here, the non-uniformity in the pre-shock region is mainly due to turbulent motions. In this type of problem, the shock propagation can reach a quasi-steady- state for a range of parameters, so that the statistical behavior of turbulence through the shock or in the post-shock region is of interest. Such a flow has been studied using shock-capturing and shock-resolving simulations for limiting cases (e.g. Lee et al., 1993, 1997; Larsson and Lele, 2009; Larsson et al., 2013; Ryu and Livescu, 2014; Livescu and Ryu, 2016). On the other hand, there have been only a few numerical studies on shock propagation in non-uniform density field. Mahesh et al. (1997); Jamme et al. (2002) have studied the pre-shock vorticity-entropy fluctuations interaction with a normal shock wave, which is strongly related to shock/turbulent boundary layer interaction. Recently, our group has investigated the shock interaction with multi-fluid variable density turbulence, and has shown that shock front, shock intensity, and post-shock turbulence are strongly affected (figure 3.1 a,b) by the non-uniformities in density, even at low Atwood number At (Tian et al., 2017a,d). By further increasing the density ratio, the interaction exhibits a broken shock regime similar to the effect seen at high turbulent Mach number Mt values (Larsson et al., 2013) (figure 3.1 c). In previous studies, the interaction was studied in regimes where the shock motion was negligible (Ryu and Livescu, 2014) or the constant backpressure outflow condition was used to stabilize the shock motion even at high Mt values or density ratios (Larsson et al., 2013) so that the statistical analysis could be carried out. In this study, the transient shock propagation in variable density flow without boundary effects is of interest, as it resembles the conditions considered in the early theoretical work and practical applications of interest (e.g. ICF). We aim to test the available 107 theories on predicting the shock wave and post-shock flow under various conditions using accurate flow resolving shock resolving/capturing simulations. We will then extend the existing theories to account for the effects of the re-reflected waves for 1D problems, with the hope of laying the theoretical foundation for more complicating 3D problems. The rest of the paper is organized as follows. The numerical details will be presented in section 3.2, including the governing equations, numerical methods, and simulation setup. In this study, we will conduct both shock-capturing and shock-resolving simulations and assess their accuracy in predicting shock motions. The theoretical approaches will then be described in section 3.3. Both LIA and Chisnell-Whitham (CW) theory will be presented and discussed. In section 3.4, the numerical results will be compared with the theoretical results for various types of non-uniform density flows. The predictability of the CW theory will be assessed based on its agreement with numerical results. An improved model will then be proposed and the 3D effects on shock propagation will be addressed. The main findings and conclusions are summarized in the last section. Figure 3.1: Shock front in wrinkled shock regime (a,b) and broken shock regime (c). Shock intensity is represented by density jump ∆ρ. Shock intensity ratio is calculated as ∆ρ/∆ρ. 108 3.2 Numerical Methodology 3.2.1 Shock-capturing Simulation Numerical detials of the shock-capturing simulations have been described in detial in the first Chapter. 3.2.2 Shock-resolving Simulation Shock-resolving simulations are conducted by computing both the viscous and inviscid fluxes using the sixth-order compact scheme. To ensure numerical stability, the viscous terms are solved in the nonconservative form and the inviscid fluxes are solved in the skew-symmetric form (Johnsen et al., 2010). A newly developed six-order compact scheme with a fourth order boundary scheme is implemented for the calculation of first-order derivative (Brady and Livescu, 2018). On top of that, the nonconservative form of the viscous fluxes contains the second-order derivative of primitive variables, which will be calculated using the six-order compact scheme with a fourth Figure 3.2: (a) Demonstration of numerical setup of the 1D simulation for both shock-capturing and shock-resolving simulations. (b) Space-time diagram of the 1D shock propgation (thick solid line) in non-uniform density media. In the meantime, a reflected wave/forward-going wave will be generated (thin solid line), which will then interact with the post-shock non-uniform density media at the particle path (dashed line). As a result, a re-reflected wave/backward-going wave will be generated (dotted line), which will catch up with the shock wave. 109 order boundary scheme as described in Lele (1992). To fully resolve all flow scales, including shock width (δ), the mesh is further refined around the shock region. At least 12 grid points are used in the streamwise direction to resolve the shock profile, in contrast to the shock-capturing simulations, where the shock wave is numerical and is represented by only 2-3 grid points. 3.2.3 Simulation Setup The numerical methods discussed above were developed and tested for 3D simulations. In this study, we will only consider the 1D cases. As shown in figure 3.2 (a), the dimension of the computational domain is 4π in the streamwise direction. The normal shock wave is initialized in the middle of the domain at x = 2π. For the shock-resolving simulation, the shock is initialized using the laminar shock profile as derived in Thompson (1971). A moving reference of frame at the same speed as the laminar shock wave is used so that the shock stays relatively stationary throughout the simulations. A buffer layer in the streamwise direction is set at the end of the open end of the domain. Since in this study, we are only considering shock propagation with no boundary effects, the simulations will end before the outflow acoustic waves propagate to the shock wave. At the inflow, a non-uniform density field is superposed on a uniform mean flow that matches the laminar pre-shock condition for a shock wave with shock Mach number Ms. Here, we will use the overbar f to denote the mean and f (cid:48) to denote the fluctuations. The density non-uniformity ρ(cid:48) is introduced by correlating the mass fraction (or mole fraction) of the binary fluid to a scalar field. Different types of density non-uniformity will be considered and discussed in the results section. To quantify the non-uniform density field, we will use the peak wave number of the density spectrum, ks and Atwood number At = (MW2 − MW1)/(MW2 + MW1) to describe the length scales and magnitude of the non-uniform medium, where MW1 and MW2 are the molecular weights of binary mixture. 110 3.3 Theoretical Analysis In this section, the theoretical approach developed in previous studies and its extension to the non-uniform medium will be summarized. Here we will categorize the non-uniform density interaction with shock wave into two different regimes based on whether the linear interaction assumption holds true in the analysis. Here, we will present and discuss the theories for both interaction regimes, i.e. Linear Interaction Approximation for the linear regime and Chisnell- Whitham theory for the nonlinear regime. 3.3.1 Linear Interaction Approximation The interaction between non-uniform density medium and a shock wave in the linear regime has been studied by considering the pre-shock fluctuations as entropy fluctuations or entropy waves (Mahesh et al., 1997). The foundation of this analysis is the Kovasznay decomposition of compressible turbulence (Kovasznay, 1953): small-amplitude turbulence fluctuations can be decomposed into vorticity, acoustic and entropy waves. Both acoustic and entropy waves contribute to the density non-uniformity, but in this section, we will only focus on the shock interaction with entropy wave induced non-uniform density field due to its resemblance to the multi-component variable density effects. The interaction of these infinitesimally small waves with a shock wave is assumed to be linear, so that the linearized Rankin-Hogoniot jump condition can be used to derive the magnitude of post-shock waves. Here, we will re-examine the theoretical analysis of shock/entropy wave interaction by Mahesh et al. (1997) and extend it to include the multi-component density fluctuations. 111 The compressible multi-component Euler equations can be obtained from equation 1.1: + ∂ ρui ∂t ∂pui ∂xi + ∂p ∂t + ∂ ρ ∂t ∂ ρuiuj ∂ ρui = 0 ∂xi = − ∂p ∂xj ∂xi = −(γ − 1)p ∂ui ∂xi = 0 ∂ ρYiuj ∂xj + ∂ ρYi ∂t (3.1a) (3.1b) (3.1c) (3.1d) In the Euler equations, the energy equation is reorganized to form an equation for pressure. The mass fraction equations can be grouped to derive an equation for specific gas constant R by using the relation: R = RunivΣ : Yi MWi ∂ ρR ∂t + ∂ ρRuj ∂xj = 0 (3.2) The linearized form of equation 3.1 can be written for small amplitude fluctuations by ignoring , the linearized nonlinear terms. By introducing the entropy fluctuation (γ−1) s(cid:48) Euler equations for fluctuations can be expressed as: = p(cid:48) p − γ ρ(cid:48) µiY(cid:48) − Σ R T ρ i ∂ ρ(cid:48) ∂t ∂u(cid:48) i ∂t + ui + uj ∂s(cid:48) ∂t ∂R(cid:48) ∂t ∂ ρ(cid:48) ∂xi ∂u(cid:48) i ∂xj = −ρ = −1 ρ ∂s(cid:48) ∂xj ∂R(cid:48) ∂xj + uj + uj ∂u(cid:48) i ∂xi ∂p(cid:48) ∂xi = 0 = 0 (3.3a) (3.3b) (3.3c) (3.3d) Note that for multi-component media, the entropy fluctuations also include the chemical poten- . It can be easily seen that the linearized Euler equation for s(cid:48) is a linear ρ and tial of the species −Σ combination of the temperature contributions to the entropy fluctuations (γ − 1) s(cid:48) p − γ ρ(cid:48) = p(cid:48) µiY(cid:48) T i t R 112 µiY(cid:48) i T c R = −Σ the composition contributions (γ − 1) s(cid:48) . The composition part of the entropy wave has the same behavior as the specific gas constant R(cid:48), due to the fact they are both linear combinations of the mass fraction equations. For the convenience of further discussion on shock interaction, we will separate the two parts of the entropy fluctuations and only consider the temperature part of entropy wave: s(cid:48) t. Equations 3.3 can then be further reorganized to obtain the Kovasznay’s decomposition into different waves and its extension to multi-component flows. + uj + uj + uj ∂ω(cid:48) i ∂t ∂s(cid:48) t ∂t ∂R(cid:48) ∂t ∂2p(cid:48) ∂xi∂t ∂2p(cid:48) ∂t2 + 2ui = 0 = 0 ∂ω(cid:48) i ∂xj ∂s(cid:48) t ∂xj ∂R(cid:48) = 0 ∂xj ∂2p(cid:48) ∂xi∂xj = 0 + uiuj −a2 ∂2p(cid:48) (cid:114) ∂xi∂xi (3.4a) (3.4b) (3.4c) (3.4d) where ωi = i jk ∂uk/∂xj is vorticity and a = γp ρ is the sound speed. Four types of independent waves can be clearly observed from equations 3.4 with their corresponding wave speed. The vorticity wave (equation 3.4a) is convected by the mean field and is only related to the solenoidal part of the velocity field, which does not contribute to the density fluctuations. The density non- uniformity comes from the other three types of waves: entropy wave (equation 3.4b), ’composition wave’ (equation 3.4c) and acoustic wave (equation 3.4d). The acoustic waves consist of pressure fluctuations p(cid:48), density fluctuations ρ(cid:48) i,d and they propagate at local sound speed with respect to the flow. The entropy wave is isobaric and can be expressed as the combination of density and temperature fluctuations and they satisfy ρ(cid:48) et/ρ = −T(cid:48)/T. Another source of the density fluctuations is the composition fluctuations or ’composition waves’, which satisfy ρ(cid:48) c/ρ = −R(cid:48)/R. Equation 3.4 shows that in the linear and inviscid regime, the composition wave in multi-component flows has the same behavior as the entropy wave. Both a and the dilatational velocity fluctuations u(cid:48) 113 types of waves affect the density fluctuations and are convected by the mean flow. Combining all the waves, the linearized equation of state can be expressed as: p(cid:48) p p(cid:48) p = = ρ(cid:48) ρ ρ(cid:48) a ρ + + T(cid:48) T ρ(cid:48) et ρ + R(cid:48) R ρ(cid:48) c ρ + T(cid:48) T + + R(cid:48) R (3.5a) (3.5b) Following the procedure in Mahesh et al. (1997), we will express the pre-shock variable density field as a combination of different plane waves. In this study, we are interested in the shock interaction with non-acoustic non-uniform density variations, so we neglect pre-shock vorticty wave and acoustic wave and only considered entorpy wave and composition wave. The extensions of Mahesh’s analysis are marked by {}. u(cid:48) 1,u u1,u u(cid:48) 2,u u1,u ρ(cid:48) u ρu T(cid:48) u Tu p(cid:48) u pu R(cid:48) u Ru = 0 = 0 = Aeeik(mx+l y−u1,ut) Aceik(mx+l y−u1,ut)(cid:111) (cid:110) + = −Aeeik(mx+l y−u1,ut) (cid:110)−Aceik(mx+l y−u1,ut)(cid:111) = 0 = (3.6a) (3.6b) (3.6c) (3.6d) (3.6e) (3.6f) (3.6g) Here the subscript u represents the pre-shock variables. Ae and Ac are the magnitude of the entropy wave and composition wave. The variables m = cosφ and l = sinφ are related to the incident angle of the plane wave to the shock. The normal shock will be deformed by the waves and the displacement is denoted using xs = ξ(y, t). Correspondingly, the linearized Rankine-Hugoniot 114 (RH) jump condition can be extended to include the composition wave: u + RuT(cid:48) u) ( R(cid:48) u Ru M(cid:48) 1 = ) = 3 (cid:41) u(cid:48) 1,d − ξt u1,u = = = u(cid:48) 2,d u1,u ρ(cid:48) d ρd p(cid:48) d pd = R(cid:48) d Rd = − − R(cid:48) u Ru ) ) + ) ) 2 + 2 ) ξy  u(cid:48) 2,u u1,u + (γ + 1)M2 1 u(cid:48) (cid:40) 2(TuR(cid:48) M1 − γM1 1,u u1,u 2u1,u u(cid:48) (T(cid:48) M1 − M1 − M1 ) 1,u u 2 2 u1,u Tu (u(cid:48) (γ − 1)M2 1 − 2 1 − ξt (γ + 1)M2 u1,u 1 (T(cid:48) u Tu ( R(cid:48) u (γ + 1)M2 Ru 1 2(M2 1 − 1) (γ + 1)M2 1 (u(cid:48) 1 − ξt 4 (γ − 1)M2 u1,u 1 + 2 −(γ − 1)M2 (T(cid:48) 1 + 4 u (γ − 1)M2 Tu 1 + 2 ( R(cid:48) 1 + 4 u Ru 1 + 2 (u(cid:48) 1 − ξt u1,u (T(cid:48) u Tu ( R(cid:48) u Ru (γ − 1)M2 4γM2 1 1 − (γ − 1) 2γM2 1 1 − (γ − 1) 2γM2 1 1 − (γ − 1) −(γ − 1)M2  2γM2 2γM2 2γM2 ) ) )  ) (3.7a) (3.7b) (3.7c) (3.7d) (3.7e) (3.7f) Here, we assume that the specific heat ratio γ remains constant in the mixture. As we can see, the pre-shock entropy wave and composition wave will generate other types of waves, including 115 vorticity wave, acoustic wave and entropy wave. The post-shock variable fluctuations can therefore be written as: u(cid:48) 1,d u1,u u(cid:48) 2,d u1,u p(cid:48) d pd ρ(cid:48) d ρd T(cid:48) d Td R(cid:48) u Ru = Fei ˜k xeik(l y−mu1,ut) Geik(mr x+l y−mu1,ut) = Hei ˜k xeik(l y−mu1,ut) Ieik(mr x+l y−mu1,ut) ei ˜k xeik(l y−mu1,ut) K γ ei ˜k xeik(l y−mu1,ut) K γ +Qeik(mr x+l y−mu1,ut) = = Seik(mr x+l y−mu1,ut)(cid:111) (cid:110)−Seik(mr x+l y−mu1,ut)(cid:111) −Qeik(mr x+l y−mu1,ut) Kei ˜k xeik(l y−mu1,ut) γ − 1 γ (cid:110) + = = with r being the compression ratio of velocity: r = u1,u u1,d (3.8a) (3.8b) (3.8c) (3.8d) (3.8e) (3.8f) (3.9) The boundary conditions across the shock wave will yield the equations for shock velocity: = Leik(l y−mu1,ut) ξt u1,u ξy = − l m Leik(l y−mu1,ut) (3.10a) (3.10b) (3.10c) The magnitude of the post-shock waves can then be related to the pre-shock waves through the 116 linearized RH equations: F + G − L = −B1L − B2(Ae + Ac) + (Q + S) = −C1L − C2(Ae + Ac) K = −D1L − D2(Ae + Ac) K γ H + I = −E1 S = Ac l m L F = αK (3.11a) (3.11b) (3.11c) (3.11d) (3.11e) (3.11f) α, B1, B2, C1, C2, D1 and D2 are considered as constants that depend on γ and M1. Based on these extensions of Mahesh et al. (1997), we are able to theoretically calculate the magnitude of the post-shock waves. From equations 3.11, we notice some physical interpretations regarding the post-shock fluctuations. Firstly, due to the linearity of the R-H jump conditions for mass fraction, the magnitude of the post-shock composition wave S is the same as the pre-shock magnitude Ac. Composition waves are conserved through the shock and cannot be generated by other types of waves. Secondly, both entropy wave and composition wave can generate acoustic and vorticity field through interaction with the shock wave. On top of that, there exist similarities in the variable density effects from the two parts. The pre-shock variable density effects from both parts with a magnitude of Ae + Ac (equation 3.6 c,d) can be grouped together on the right hand side (RHS) of equation 3.11 (a-c). This implies that if the total pre-shock density fluctuations are fixed, i.e. Ac + Ae =constant, the post-shock pressure, density and velocity fluctuations will be the same, regardless of the individual magnitude of Ae and Ac. This has also been observed using transfer matrix by Griffond (2005). On the other hand, the effects both entropy wave and composition wave on post-shock entropy wave will be different depending on the values Ac and Ae. The magnitude of the composition wave is the same through the shock wave, while that of entropy wave is affected by both Ae and Ac. Note that the above analyses are based on the linearized Euler equations and the assumptions that specific heat ratio γ is constant. If scale separation between the shock wave and pre-shock disturbances is not satisfied, the viscous effects can play a role in the interaction, 117 which will make the post-shock field deviate from the LIA limit (Ryu and Livescu, 2014). Or if γ varies in the binary mixture, the analysis will also be further complicated (Griffond, 2005). 3.3.2 Chisnell-Whitham Theory In many of the engineering applications, At or the magnitude of the density non-uniformity is too large, which will invalidate the assumptions of LIA. For example, the density fluctuations ρ(cid:48)/ρ are O(1) in the hydrogen/air mixture in the scramjet engine. In the nonlinear regime, a different approach is needed to predict shock propagation in non-uniform density media. Early theoretical studies by Chisnell (1955); Whitham (1958) has considered the large density fluctuations as a sequence of weak contact discontinuities and derived a relation between the density change of dρ with the shock strength change dz, where z is the ratio of pressure across the shock wave. The underlying assumption is that the shock is quasi-equilibrium and the density change across the contact discontinuity is weak. Figure 3.2 (b) shows a wave diagram of the shock propagation in a non-uniform density media. As the shock wave encounters a density change dρ, its Mach number and propagation speed will be modified. In the meantime, a reflected wave/forward-going wave will be generated (marked as thin solid line), which will then interact with the post-shock non-uniform density media at the particle path (dashed line). As a result, a re-reflected wave/backward-going wave will be generated (dotted line), which will propagate in the same direction as the shock wave and eventually catch up the shock wave. The relation between upstream density jump dρu and dz can be written as (Chisnell, 1955): where λ =(cid:112)(γ − 1)/(γ + 1) is a constant. In equation 3.12, the re-reflected waves are neglected. ( 1 + λ2z z(1 + λ2))1/2 2 z − 1 1 ρu dρu dz λ2 + z = 2 z − 1 − 1 + (3.12) The differential equations can be integrated to obtain the explicit relation between shock strength 118 and density. ρ = A (z − 1)2 (λ2 + z) (cid:40)(cid:112) (cid:112) λ2 + 1/z + λ λ2 + 1/z − λ (cid:41)2λ/(cid:113)(1+λ2)(cid:40)√ √ λ2 + 1 −(cid:112) λ2 + 1 +(cid:112) λ2 + 1/z λ2 + 1/z (cid:41)2 (3.13) where constant A can be calculated from the initial condition. In this formulation, the pre-shock density non-uniformity is assumed to be either entropy fluctuations or composition fluctuations. This implies that in the nonlinear interaction regime, the density non-uniformity resulted from both temperature fluctuations and multi-component effects are also equivalent. On top of that, the shock strength z is a direct function of ρ, meaning that shock strength should be independent of the shape of density profile. This equation can then be numerically integrated over any given density profiles to obtain the temporal evolution of the post-shock medium, which can then be used as a boundary condition for the solution of the post-shock turbulent field. Due to the fact that the dependence of shock strength on density is implicit, z needs to be numerically solved using an iterative method. In this study, the numerical solution of z is considered to be converged when the error is smaller than 10−8. Here for simplicity, the dependence of z on ρ will be represented using implicit function z = g(ρ). Therefore, the shock position can be obtained using the following integration: Þ Þ (cid:112) Þ (cid:112) xs = = = dxs = usdt = γP/ρ(xs) γP/ρ(xs) Þ a(xs)M(xs)dt Þ ((z − 1)(γ + 1) + 1)dt ((h(xs) − 1)(γ + 1) 2γ (cid:115) (cid:115) + 1)dt 2γ (3.14) By integrating the above equations numerically, the shock position xs as a function of t can be obtained. The interaction between weak contact discontinuity with a shock wave will also generate Its magnitude can be expressed as acoustic waves that propagate away from the shock wave. (Chisnell, 1955): 119 z f = 1 + dz z (3.15) where dz can be related to the density jump through equation 3.12. When this elementary wave interacts with the post-shock non-uniform field, the re-reflect waves with be generated, which will propagate towards the shock wave and change the shock strength. Chisnell (1955) derived the magnitude of the re-reflected waves as: here ρ(z, z0) can be calculated using the following equation (Chisnell, 1955): zb = 1 − 1 1 4 α(z, z0) α(z, z0) = ( ρ(z,z0) dz z dz0 ρ(z,z0) ∂ ρ(z,z0) ∂z0 ρ(z,z) )1/4 Þ z (3.16) (3.17) (3.18) −1/4(z, z0) = ρ ρ −1/4(z0, z0) − 1 4γ −1/4(x, x) dx x ρ x=z0 where ρ(z, z) is the density immediately after the shock wave for a given shock intensity z. Based on these equations, the magnitude of the fast and slow acoustic waves can be calculated through numerical integration for a given density profile. In this study, we aim to test the applicability of these theoretical results to various types of problems with the hope of finding a universal and simple model for the shock propagation in non-uniform density field. 3.4 Results and Discussion In this section, the numerical results obtained from both shock-capturing and shock-resolving simulations will be examined and compared with theoretical results. 3.4.1 Shock Wave Propagation in Media with Linearly Varying Density In this section, the simple case of shock propagation in a linearly varying density field is considered. The aim is to verify the accuracy of current shock-capturing simulations in predicting shock 120 propagation in variable density media. Pseudo-1D simulations for linearly varying density profile in the weak shock limit of M0 = 1.1 are performed first. The ratio of the density at the beginning the field to that at the end of the field is 8 (for linearly decreasing density field) and 0.125 (for linearly increasing density field). Figure 3.3 (a) compares the time development of shock position obtained from numerical results and CW solution. The x-axis denotes time and the upper curve and lower curve represent cases with increasing and decreasing density field, respectively. As observed, in the weak shock limit, the results from current shock-capturing simulations are in excellent agreement with the CW solution for both decreasing and increasing density profiles. On the other hand, the match with the CW solution is less than satisfactory for a strong shock (M0 = 4). In this case, the numerical results for the case with increasing density gradually deviate from the CW solution in time (figure 3.3), resulting in a final difference of around 5%. The case with decreasing density shows a relatively better agreement at M0 = 4. To further understand the discrepancy, we compare the results from current simulations with those obtained by Bird (1961), who used the method of characteristics. The exact same conditions are used as in Bird (1961), for Mach number 4.0 and density ratios of 8 and 0.125. The comparison of the shock Mach number as a function of distance is shown in figure 3.4. Evidently, the current results agree better with those obtained by the method of characteristics for the increasing density case. On the other hand, for the linearly decreasing density case, all three methods give comparable predictions. It was argued in Bird (1961) that the influence of re-reflected waves, which is not considered in the CW solution, is different for different density profiles. This explains the different performance of the CW model in different situations. Our numerical results agree well with those obtained by the method of characteristics and also with CW theory when the re-reflected waves are not important. Note that the solutions are obtained by solving the Navier-Stokes instead of Euler. This implies that the diffusion effects are indeed week through the interaction for the parameter range considered in these cases. 121 3.4.2 Shock Wave Propagation in Media with Fluctuating Density The fundamental understanding of shock propagation in 1D fluctuating density is critical to the study of shock-turbulence/shock-boundary layer interactions. This problem provides a simplified configuration for the more complex 3D case with turbulence. The propagation of shock wave into non-uniform density with a sinusoidal shaped profile is called the "Shu-Osher" problem, which can be theoretical solved using the formulas in the linear interaction regime. This has been a classical problem for the validation of shock-capturing numerical schemes (Shu and Osher, 1989). In this section, we will use it as a simplified case for the study of shock propagation in randomly fluctuating density field, especially in the nonlinear regime. 3.4.2.1 The Similarity of Density Fluctuations from Entropy and Composition Fluctuations In section 3.3, the similarities and differences entropy wave and composition wave through shock wave have been discussed under certain conditions using theoretical derivation. In this section, we will use the numerical results to compare the density effects from both temperature fluctuations and composition fluctuations. Here, we will briefly describe the numerical setup for the "Shu-Osher" problem. The shock Figure 3.3: Shock wave propagation in the medium with linearly varying density for density ratios of 8 and 0.125. (a) M0 = 1.1 and (b) M0 = 4.0. 122 is initialized in a fluid with ρ = 1.0 at x ≈ π; the shock Mach number is 3.0; then a fluid with density profile given by ρ = 1.0 + Atsin(ksx) is superposed on the mean flow at the inlet, where At = 0.25 is the Atwood number. This density profile is obtained using two different methods: a) the composition of the flow is varied with fixed T and p and b) the temperature is varied with fixed Yi and p. The shock-resolving method is used for the simulations. Figure 3.5 shows the pressure, density and temperature in the streamwise direction. As shown in the figure, the post-shock pressure and density from both cases are exactly the same. This is in good agreement with the previous theoretical analyses. We also note that the post-shock density field has both short and long wavelength fluctuations. This is caused by the superposition of post-shock entropy/composition waves and acoustic waves and the short and long wavelength fluctuations are due to their different characteristic speed. The ratio of the wavelength is a function of the shock Mach number. Its dependence on Ms is plotted in figure 3.7 and confirmed with numerical results. On the other hand, the post-shock temperature fluctuations are different and the Figure 3.4: Evolution of shock Mach number in linearly varying density field for ρ1/ρ0 = 8 and 0.125. Initial shock Mach number M0 is 4.0. 123 x00.20.40.60.81Ms234567present workChisnell-WhithamBird (1961)ρ1=8M0=4ρ0=1ρ1=0.125 phase of the temperature fluctuations are reversed. This shows the difference between composition waves and entropy waves in the interaction with a normal shock wave. 3.4.2.2 Shock Propagation in Sinusoidal Density Media It was shown in the last section that the non-uniform density effects from both temperature or composition fluctuations through the shock wave are similar. Therefore in this section, we will only consider the variable composition cases. In figures 3.6 (a,b), the temporal development of shock position is shown for a Mach 2 shock and small density ratio of 2:1 (At = 0.33) at different wavenumbers of ks = 1 and 4. The CW solution is obtained by integrating equation 3.14. As shown in the figure, the CW solution predicts a shock motion that resembles the sine wave variation in the density profile, but with a constant mean velocity relative to the unperturbed shock that moves the shock wave upstream. The match between the numerical and theoretical solution is fairly good during the first (sine) wave cycle. However, after the first cycle, the shock wave acquires additional oscillatory motions Figure 3.5: Plots for flow variables in streamwise directions for a Mach 3 shock propagating in a sinusoidal shaped density field. (a) Pressure, (b) density and (c) temperature. 124 instead of moving upstream following the original sine wave only, as predicted by the theory. Thus, the current shock-capturing simulations are predicting long-wavelength oscillatory motions superposed on the original sinusoidal motion. Shock-resolving simulations are also conducted to verify this observation. As seen in figure 3.6 (b), shock-resolving results match very well with the shock-capturing results, proving that the deviation from CW solution and the long-wavelength oscillatory motions are both physical. The ratio of the wavenumber, kl, of the first harmonic of these motions to the sine wave wavenumber, ks, is denoted as rl = kl/ks ≈ 1/4. This ratio appears the to be the same for ks = 1 and 4. To explain these long-wavelength motions, we consider two velocities that are related to the behavior of shock: 1) the velocity at which the pre-shock medium approaches the shock wave vp, and 2) the velocity of the re-reflected waves approaching the shock wave, vr, relative to the shock speed. vp carries the density fluctuations responsible for moving the shock wave. vr can be found from the local sound speed behind the shock. The re-reflected waves can eventually overtake the initial shock and affect its strength and speed (Bird, 1961). The ratio of vr/vp is denoted as rv, and we find that rv ≈ rl. For a weaker shock, as shown in figure 5(c), rv and rl have small magnitudes: rv ≈ rl ≈ 1/12. Thus, the early time agreement with the CW model is improved and the long-wavelength variation in shock position is also clearer. It can be inferred that the re-reflected waves causing the long-wavelength motions of the shock (Bird, 1961) are responsible for deviation from the theoretical solution. For a fluctuating density profile, it seems that the re-reflected waves play a more important role in changing the shock strength and the shock motion, making the shock wave further deviate from the CW model. This hypothesis can be verified by theoretically deriving the wavelength ratio of the two effects. This is done by matching the frequency of the waves and then relate the wavelength to the wave speed: (cid:115) ( 1 s −(γ−1)/2 γM2 1+(γ−1)M2 s /2 − 1) rms = 2+(γ−1)M2 s (γ+1)M2 s (3.19) From equation 3.19, we can see that the wavelength ratio is only a function of shock Mach 125 number Ms and specific heat ratio γ. In this study, we have fixed γ to be 1.4, so the wavelength ratio will only be a function of Ms. To verify the applicability of this equation, we have conducted numerical simulations covering a wide of range of Ms. The numerical results (black circle) are compared with equation 3.19 (black solid line) in figure 3.7. As we can see, the theoretically derived wavelength ratio decreases as the Ms increases, reaching an asymptotic value ((cid:112)2γ + (cid:112)γ − 1)/(cid:112)γ − 1. The numerical results from shock-resolving simulations are marked in the figure using black circles. Results show that the agreement with the theoretical derivation is very good, which shows that the long-wavelength movement of the shock is caused by the re-reflected slow acoustic waves. It has been pointed out that in figure 3.6, the CW prediction of the first sinusoidal waves agrees very well with the numerical results. We have shown that the wavelength ratio decreases to an asymptotic limit as the shock Mach number Ms increases. For γ = 1.4, the asymptotic value is around 3.65, which still provides large enough separation between the two types of waves. This Figure 3.6: Shock propagation in a fluctuating density field. (a,b) M0 = 2.0, density ratio 2.0, and ks = 1 and ks = 4 and (c) M0 = 1.1, density ratio 2.0, and ks = 4. 126 means that the CW solution that neglects the re-reflected waves could provide a good estimation for the short wave-length shock oscillation. The magnitude of the short wavelength shock oscillation can be obtained by integrating equation 3.14 over the very first sin wave. The wavelength of the pre-shock density profile is fixed, so the results will be a function of shock Mach number Ms and Atwood number At. The shock position oscillation from both numerical simulations and numerical integration of CW solutions are compared in figure 3.8. The CW theoretical results show that as the At number increases, the shock oscillation becomes larger. For cases with different Ms, lower Ms case will have a relatively larger shock oscillation. The numerical results are also marked in the figures. It can be seen that the numerical results agree with the theoretical results. This confirms that the CW solutions provide as a reasonably good method for predicting the 1D shock oscillation in the nonlinear regime. Figure 3.7: Comparison of the wavelength ratios between numerical results and theoretical results. 127 3.4.2.3 CW Solution with Correction for Re-reflected Waves In the above analyses, the effects of the re-reflected waves are neglected in the CW solutions. Without taking the re-reflected waves into consideration, the theory still gives acceptable predictions for short wave-length shock oscillations. To further improve the applicability of the theory, the effects from the re-reflected waves need to be considered. Chisnell (1955) has considered the catching up of the re-reflected waves for the linearly increasing density profiles. However, the formulas are not closed, making it impossible to theoretically calculate the shock strength. In this section, we will extend some general formulas in Chisnell (1955), i.e. equations 3.15 and 3.16 to propose a theoretical model for the estimation of the re-reflected wave effects on shock propagation in fluctuating density profiles. The magnitude of the generated acoustic waves immediately after the shock wave can be estimated by numerically integrating equation 3.15. The magnitude of the wave can be expressed as: Figure 3.8: Comparison of the shock position fluctuations between numerical results and theoretical results. Results from shock-resolving simulations are maked using symbols for Ms = 1.1 (black circle), Ms = 1.3 (black cross), Ms = 1.5 (black star), Ms = 2.0 (red circle), Ms = 2.5 (red corss), Ms = 3.0 (red star). 128 Þ x 0 dz z(xs) z f (xs) − 1 = (3.20) Similarly, the magnitude of the re-reflected waves can be estimated by integrating 3.16 over dz0 Þ z 0 and then dz: zb(xs) − 1 = − 1 4 α(z, z0) 1 ρ(z, z0) ∂ ρ(z, z0) ∂z0 dz z dz0 (3.21) (λ2+z) 1+λ2z and ρ(z, z) is the post-shock density when the shock strength is z and can be calculated using the formula ρ(z, z) = ρ . From these equations, the profiles of the numerically integrated waves can be obtained. To estimate the applicability of these estimations, the results from shock-resolving simulations will be compared with the theoretical results. In order to decompose the post-shock waves from numerical simulations, a wave decomposition method as proposed by Payri et al. (1995) is used. The forward propagating (or downstream the shock wave) waves (p f )or backward propagating waves (pb) can be calculated using: p f = p0[1 2 (1 + ( p p0 ) γ−1 2γ (1 + γ − 1 2 u1 a )] 2γ γ−1 pb = p0[1 2 (1 + ( p p0 ) γ−1 2γ (1 − γ − 1 2 u1 a )] 2γ γ−1 (3.22) (3.23) Therefore, based on the definition of wave strength z, z f can be calculated to using p f ,max/p f ,min and zb can be calculated as pb,max/pb,min. The relative importance of the re-reflected waves can be estimated using the ratio of the magnitude of the two waves zb/z f . The ratio zb/z f from both theoretical derivation and numerical results are plotted in figure 3.9. The theoretical results predict that as the initial shock Mach number Ms increases, the ratio of zb/z f becomes smaller. And as the magnitude of density non-uniformity At increases, the relative importance of re-reflected waves zb also becomes weaker. Numerical results, as marked in the figure using symbols, matches very well with the theoretical results for both shock-capturing 129 simulations and shock-resolving simulations. This proves that the prediction of the magnitude of the post-shock re-reflected waves provided by Chisnell (1955) is accurate. On the other hand, the relative magnitude of the re-reflected wave decreases as the shock Mach number Ms increases and At increases. The correction to the shock propagation was proposed by (Chisnell, 1955) to compensate for the effects of the re-reflected waves. The strength of the re-reflected waves zb can be incorporated into the original equation. We have shown that the magnitude of the re-reflected waves can be predicted using the above formulation. Here, we will propose several simplifications of Chinsnell’s correction to the shock propagation. As shown in figure 3.10, the elementary re-reflected wave generated by a forward-going wave dz1 interacting with the non-uniform density field generated by dz0 can be expressed using equation 3.21. This elementary wave will propagate upstream and eventually changing the shock intensity at dz. Therefore, the modification of shock strength dz is not only affected by the pre-shock density change dρ, but also affected by the re-reflected waves along the negative characteristic lines, which Figure 3.9: Comparison of the relative magnitude zb/z f of the re-reflected waves.Results from shock-resolving simulations are maked using symbols for At = 0.1 (black circle), At = 0.2 (black cross), At = 0.333 (black star), At = 0.4 (red circle), At = 0.5 (red corss), At = 0.6 (red star), At = 0.9 (blue circle). 130 1 − dzÞ z0=z z0=0 1 − dzrere = is marked as the dash-dotted line in figure 3.10. From here, we will use dzre to represent the modification of the shock wave strength after considering the re-reflected waves. After getting the magnitude of the re-reflected wave dzrere that catches up with the shock wave, dzre can be estimated using: = ( 2 dρu ρu z−1( 1+λ2z z−1 − 1 + 2 λ2+z z(1+λ2))1/2dzrere z−1( 1+λ2z + 4 z(1+λ2))1/2)dzre (3.24) The only remaining question is the estimation of dzrere along the negative characteristic line. In Chisnell’s theory (Chisnell, 1955), the magnitude of the re-reflected wave which catches up the incident shock is: 1 z1 dz0 1 ρ(z1,z0) ∂ ρ(z1,z0) ∂z0 dz1 dz (3.25) 4 α(z, z1) 1 The only unknown in the integration is dz1 dz , which will make the Chisnell’s formulation not closed. This ratio dz1 dz represents the portion of the elementary re-reflected waves that will catch up and then modify the strength of the initial shock wave. This is determined by the relative speed of the shock wave and the re-reflected wave. Here we assume that the post-shock acoustic speed is constant so that the characteristic lines are straight, which implies dz1 = constant. The dz above assumption is reasonable when the shock wave is propagating in a fluctuating density field with a fixed mean value. By making the above assumption, one can easily integrate the re-reflected waves along the characteristic line without the need of solving for the full post-shock field. = ub−us us After obtaining the magnitude of the re-reflected waves, it can be substituted into equation 3.24 to obtain dzre, and therefore the modified shock strength zre(t). The modified shock position can then be calculated by replacing z(t) with zre(t) in the numerical integration of equation 3.14. The modified CW solution will be referred to as CW with re-reflected waves (CWRW). Figure 3.11 shows the results for shock-resolving simulation, CW, and CWRW model. It can be seen that by including the effects of re-reflected waves, the CWRW model can significantly improve 131 the prediction of shock propagation in a sinusoidally shaped density field. Evidently, both the mean deviation from the CW solution and the long-wavelength oscillation are predicted by the model. 3.4.2.4 Shock Propagation in Random Density Media The above analysis based on simplified density profiles can provide some basic assessment of the CW and CWRW models. However, due to the nonlinear nature of the problem, one should not expect that in a random realistic density media, the CWRW model solutions would still be valid. In this section, we will conduct an analysis on a randomly generated density profiles that satisfy a given spectrum and compare both CW and CWRW models with the numerical results. Here, we will use the same spectrum as used in Tian et al. (2017d) for the 1D density field. The characteristic length scale of the density field is modified by changing the peak wave number of the spectrum. The generated density field will then be used as the inflow condition for both numerical simulation and numerical integration of CW and CWRW models. Figure 3.12 shows Figure 3.10: Demonstration of the catchup of the re-reflected waves. 132 the shock propagation in a realization of the random density field with ks = 4.0 and At = 0.5. The numerical solutions are plotted using solid black lines and CW solutions are represented as dashed lines. As we can observe, the CW solutions predicted an overall similar profile as the numerical simulations. However, they still gradually deviate from the numerical solution, the same as we have observed in figure 3.6 (b). We also note that the long wavelength oscillation are not visible. On the other hand, the CWRW model shows improvement over the CW solution and agrees better with the shock-resolving simulation predictions. The results confirm that the CWRW model provides a significant improvement over the CW theory by including the re-reflected wave effects in 1D shock propagation. 3.4.3 Shock Wave Propagation in Three-dimensional Non-uniform Density Field In most real applications, the non-uniform density distribution/media is three-dimensional. The interaction between small amplitude 3D density fluctuations with a normal shock wave has been Figure 3.11: Shock propagation in sinusoidally varying density field. The temporal development of shock position is compared for the shock-resolving simulation, CW theory and the proposed CWRW theory. 133 theoretically studied by Griffond (2005); Huete et al. (2011) using linear theory. These studies have been focusing on the prediction of the post-shock field rather than the shock propagation and wrinkling. The amplitude of density fluctuations or At is usually too large in real-life applications, making the linear theory inappropriate for predicting shock motions. In the above discussions, we have considered the application of CW theory and CWRW model in predicting shock propagation in 1D non-uniform density media. The advantage of CWRW model is that it can be extended to the nonlinear regime. When the 3D spatial effects in the spanwise directions play a role, one may not simply treat each point on the 3D shock as individual 1D normal shock. In this section, we will present results and discussions to describe the 3D effects of shock propagation in non-uniform media. Figure 3.12: Shock propagation in random density media. Different realizations of the random density field are represented using different colors.(a) ks = 4.0 and (b) ks = 8.0. 134 3.4.3.1 3D Sinusoidal Density Field k2 y + k2 (cid:113) kyz , where kyz = The shock propagation in a realistic 3D random non-uniform media is understandably complicated. To simplify the problem, here we consider a 3D sinusoidal shaped density field that satisfies ρ(x, y, z) = 1.0 + Atsin(kx x + ky y + kzz). By using this 3D density field, one can prescribe At and kx, ky, kz to modify the intensity and length scale of the non-uniform density field. The anisotropy of the non-uniform density can also be changed by modifying the ratio of the streamwise to spanwise length scales rk = kx z . Using coordinate transformation, one may show that the 3D sinusoidal shaped density field is essentially a 2D plan wave with At, kx, and kyz being the important parameters. When the wavenumber ratio rk is equal to one, the density field becomes isotropic. When rk increases and approaches infinity, the density field becomes quasi-1D. On the other hand, when rk approaches zero, the spanwise shock surface wrinkling will dominate the shock propagation. Fundamental understandings of these two limits and how the shock behaves when the ratio rk is varied are critically important to the modeling of shock propagation in 3D non-uniform media. The following results and discussions on the 3D effects are divided into two parts. The limit when scale separation ratio rk approaches infinity will be studied first, followed by an examination of the limit when rk goes to zero. In figure 3.13, the instantaneous shock positions for several selected points on the shock surface are plotted as a function of time. We have verified that the selected points on the shock surface can represent the overall behavior of the shock wave, so the following discussions and conclusions may apply to the whole 3D shock surface. For all the cases in figure 3.13, At is very close to 0.5 and the length scale in the streamwise direction is fixed with kx = 4. The 3D effects are studied by varying ky and kz. Unless otherwise stated, ky is equal to kz. The results form CW theory and CWRW models are also plotted. Note that both CW theory and CWRW models are for 1D shock propagation, so in this 3D case, we have made the assumption that the selected point on the shock wave is only affected by the non-uniform density field right before that point. This essentially implies that the corresponding point on the shock is moving independently like a 1D normal shock. We note that when ky, kz is reduced to zero, the 3D simulation becomes quasi-1D. As expected, 135 the predicted values fall onto the CWRW model curve, which deviate from CW theory on average. As ky, kz are increased from 0 to 1 and 2, the results do not show significant deviations from the quasi-1D (ky = kz = 0) case and CWRW model. This indicates that with a wavelength ratio greater √ 2, the shock motion can be approximated by the 1D theory and the 3D effects are than or equal to generally negligible. Further increasing k y, kz to 3, however, significantly changes the behavior of the shock wave. In addition to the sinusoidal-shaped fluctuations, the shock wave moves forward with a mean velocity similar to that predicted by the CW theory. The sudden change in the shock 2) is hypothesized to be caused behavior when the wavenumber ratio is decreased from by the nonlinear effects induced by the 3D non-uniform density field. When the density field becomes more isotropic, the effects from the post-shock re-reflected waves seems to be weakened. Additionally, the 3D effects somehow further increase the mean shock propagation speed when ky, kz are further increased. We have conducted similar simulations for kx = 1 and 2 with same wavenumber ratio. The results (not shown) are found to be in good agreement with the results √ √ 2 to 4/(3 Figure 3.13: Shock propagation in 3D sinusoidal shaped density media. The temporal development of shock position on certain point on the shock surface is plotted. kx = 4 for all the cases and ky = kz unless the values are indicated by the legend. 136 t00.511.52Xs-1.2-1-0.8-0.6-0.4-0.20ky=0(1D)ky=1ky=2ky=3ky=3,kz=4ky=5,kz=0ky=4ky=5ky=6ky=10,kz=0ky=8ky=16Chisnell-WhithamCWRW above, implying that the wavenumber ratio is the only important parameter. To study the other limit (when rk approaches zero), the length scales in the spanwise directions are further decreased by increasing ky, kz. By doing so, we observe a smooth increase in the mean shock speed until ky = kz = 8 or a wavenumber ratio of 0.35. Results from a smaller wavenumber ratio have almost the same mean propagation speed. This suggests that when the wavenumber ratio is less than 0.35, the mean shock propagation speed reaches another limit and in this limit, the 3D effects in the streamwise direction do not play a role in the mean shock propagation. To better show the dependence of the shock propagation speed to the wavelength ratio, all the results for kx = 4 with varying ky and kz are plotted on figure 3.14. The shock speed is normalized using the mean sound speed in the variable density media. A hyperbolic tangent function is fit to the data points using least square fitting. The two limits on mean shock propagation for wavelength ratios rk = 0 and ∞ are not well captured by the fitted function but can be inferred from figure 3.14. For the limit where rk approaches infinity, the mean shock propagation speed reaches a small non- zero value, much smaller than that predicted by the CW theory. This limit is essentially a pseudo-1D case, which can be numerically approached by setting ky = kz = 0. For this limit, the improved Figure 3.14: The dependence of mean shock speed on wavenumber ratio for shock propagation in 3D non-unform density media with kx = 4. 137 rk01234Us/a00.020.040.060.080.10.120.140.16 CWRW model proposed earlier gives a better prediction of the mean shock propagation speed than the CW theory. For the other limit, results show that when rk approaches 0, the shock propagation speed becomes independent of the wavelength ratio and reaches a maximum speed around Mach 0.175. Dimensional analysis show that the mean shock speed will be a function of At and Ms, i.e. Us(At, Ms). For this limit, the smallest wavenumber ratio achieved in the numerical simulation is 0.044. Unlike the other limit, one can not simply set the streamwise wavenumber kx = 0, because this will result in a broken shock wave. It is also problematic to continuously increase the spanwise wavenumber ky, kz in simulations, because the wavelength of the fluctuation will keep decreasing and eventually reaching a scale where the diffusion effects become very important. In that case, the shock wave propagation speed will depend on Reynolds number Re. Previous simulations satisfy the scale separation criteria between the characteristic shock length scale and non-uniform density length scale, so the viscous effects can be neglected in the shock interaction. This essentially means that the above observations are for Reynolds number being infinitely large across the shock. When the small Re starts to affect the results, the predicted Us will not be the correct inviscid rk = 0 limit. In this study, the limit rk = 0 can only be approximated through simulations with rk = 0.044 with reasonably large length scale separation. The dependence of asymptotic value of mean shock propagation at the limit of rk ≈ 0 on Atwood number is also studied with fixed initial shock Mach number. The results are plotted in figure 3.15. We note that the mean propagation speed increases when the Atwood number is increased, which is mainly due to nonlinear effects. When At = 0, the mean shock propagation speed is zero, so the shock would remain stationary. As the Atwood number increases, the nonlinear effects introduced through the non-uniform density field will increase, causing the increase in mean shock propagation speed. In 1D cases, the re-reflected waves cancel out part of the nonlinear effects, making the shock propagation speed smaller than the 3D cases. By removing the mean propagation from the instantaneous shock position, one can obtain the magnitude of fluctuations in shock position at different spanwise locations, which can be used to describe the shock surface wrinkling. Figure 3.16 shows the magnitude of the fluctuations in 138 shock position (L2(∆xs)) as a function of the wavelength ratio. A least-square-fit to the function (c(1)rk + c(2)) × ec(3)r2 k − c(2) is also shown. Here we note that L2(∆xs) peaks around rk ≈ 0.57 and continuously decreases as the wavelength ratio increases. The asymptotic limit at rk = ∞, representing the 1D case, is not accurately captured by the fitted function but can be reasonably well predicted using the CWRW model. On the other hand, when the wavelength ratio decreases and approaches zero, the spanwise length scale of the non-uniform density field becomes smaller. Numerical results show that this will cause the shock oscillations to decrease. As the wavelength ratio approaches this limit, the fluctuations in the spanwise directions become more important to the shock wrinkling. Therefore, the decrease in length scale of the spanwise density fluctuations causes the decrease in the shock surface wrinkling. 3.4.3.2 3D Random Density Field For a 3D random density field, the temporal shock propagation becomes extremely difficult to predict, because there exists a wide range of spanwise and streamwise length scales. Despite the difficulty in the accurate prediction of shock position, previous results show that for a fixed At and Figure 3.15: The dependence of mean shock speed on wavenumber ratio for shock propagation in 3D non-uniform density media with kx = 4. 139 At0.20.30.40.50.60.7Xs00.10.20.30.40.50.60.7 Ms, the shock propagation speed should be bounded by the two limits obtained in the previous section. In this study, a shock-capturing simulation was conducted to study shock propagation in a randomly generated density field. The generation of the 3D density field has been described in the first chapter. Here, we have kept the peak wavenumber of the density spectrum to be ks = 4.0. The Atwood number is set to be 0.5. Figure 3.17 shows the mean shock propagation for many individual points on the shock surface. The two bounds for the Atwood number of 0.5 and shock Mach number of 2 are also shown. Interestingly, the results show that for each individual point on the shock, the mean propagation is bounded by the two limits, even though the density field is generated randomly. This confirms that even though predicting the instantaneous shock surface response to the random density fluctuations is very challenging, one can still provide bounds for the shock position at every point on the shock surface. 3.5 Conclusions In this study, accurate flow-resolving, shock-resolving/shock-capturing simulations are con- ducted together with theoretical analyses to study the shock propagation in non-uniform density Figure 3.16: The dependence of L2 norm of the shock fluctuations on wavenumber ratio for shock propagation in 3D non-unform density media with kx = 4. 140 rk012345L2(∆Xs)0.010.020.030.040.050.060.07 media. Theoretical analyses show that in both linear and nonlinear regime, there exist similarities between shock propagation in non-uniform density fields induced by entropy and composition fluc- tuations. In the nonlinear regime, the CW theory gives explicit solutions for the shock propagation in 1D non-uniform density media by neglecting the post-shock field effects, i.e. the re-reflected waves. Pseudo-1D numerical simulations for shock propagation in linearly varying density profiles show reasonably good agreement with the CW solutions as well as the results in the literature, which confirms the accuracy of our current numerical results. For oscillatory density profiles, the numerical results further deviate from the theoretical solutions and exhibit additional long- wavelength oscillations, which are shown to be related to the re-reflected waves. The magnitude of the re-reflected waves is numerically integrated by extending the original Chisnell’s formulation, which shows excellent agreement with the numerical results. A simplified model that includes the effects of re-reflected waves is then proposed to compensate for this deviation. The results obtained by the proposed model show significant improvement over the shock wave propagation in fluctuating density field predicted by the classical CW theory. The shock propagation in 3D non-uniform density media is also studied using a simplified 3D Figure 3.17: The mean propagation of shock wave in 3D random density media. Red lines represent points on the shock surface. The blue and black curves represent the bounds at the two asymptotic limit. 141 t00.511.52Xs-0.8-0.6-0.4-0.200.2 sinusoidal shaped density profile. The 3D effects are investigated by varying the wavenumber ratio of the streamwise and spanwise wavenumbers of the non-uniform density field. The asymptotic limit for wavelength ratio approaching zero and infinity are investigated and used to understand the 3D effects on the shock propagation. It is shown that when the wavenumber ratio is larger than 1.4, the results are very close to the proposed CWRW model, which implies one can treat each point on the shock wave as an independent 1D shock. Decreasing the wavenumber ratio will result in an increasing shock propagation speed. When the wavenumber ratio approaches zero, the shock propagation speed reaches another asymptotic limit. The shock propagation speed is shown to be bounded by the two limits for rk approaches 0 and infinity at any fixed At and Ms. The magnitude of the shock wrinkling is also found to be dependent on rk, reaching its maximum value at rk ≈ 0.57 for At = 0.5 and Ms = 2. Numerical study of shock propagation in 3D random density field indicates that the shock propagation speed on any point at shock surface is bounded by the two limits observed in our 3D simulations with sinusoidal density profile. 142 BIBLIOGRAPHY 143 BIBLIOGRAPHY Biferale, L., Boffetta, G., Celani, A., Devenish, B. J., Lanotte, A., and Toschi, F. (2005). Multipar- ticle dispersion in fully developed turbulence. Phys. Fluids, 17(11):111701. Bird, G. A. (1961). The motion of a shock-wave through a region of non-uniform density. J. Fluid Mech., 11(2):180–186. Bisnovatyi-Kogan, G. S. and Silich, S. A. (1995). Shock-wave propagation in the nonuniform interstellar medium. Rev. Mod. Phys., 67(3):661–712. Boukharfane, R., Bouali, Z., and Mura, A. (2018). Evolution of scalar and velocity dynamics in planar shock-turbulence interaction. Shock Waves, 28(6):1117–1141. Brady, P. T. and Livescu, D. (2018). High-order, stable, and conservative boundary schemes for central and compact finite differences. Comput. Fluids, 183:84–101. Brouillette, M. (2002). The Richtmyer-Meshkov instability. Annu. Rev. Fluid Mech., 34(1):445– 468. Budzinski, J. M., Zukoski, E. E., and Marble, F. E. (1992). Rayleigh scattering measurements of shock enhanced mixing. AIAA Paper 92-3546. Chisnell, R. F. (1955). The normal motion of a shock wave through a non-uniform one-dimensional medium. Proc. R. Soc. Lond. A, 232(1190):350–370. Choi, J. I., Yeo, K., and Lee, C. (2004). Lagrangian statistics in turbulent channel flow. Phys. Fluids, 16(3):779–793. Chong, M. S., Perry, A. E., and Cantwell, B. J. (1990). A general classification of three-dimensional flow fields. Phys. Fluids A: Fluid Dynam., 2(5):765–777. Chong, M. S., Soria, J., Perry, A. E., Chacin, J., Cantwell, B. J., and Na, Y. (1998). Turbulence structures of wall-bounded shear flows found using dns data. J. Fluid Mech., 357:225–247. Chu, Y. B. and Lu, X. Y. (2013). Topological evolution in compressible turbulent boundary layers. J. Fluid Mech., 733:414–438. Cook, A. W. and Riley, J. J. (1996). Direct numerical simulation of a turbulent reactive plume on a parallel computer. J. Comput. Phys., 129(2):263–283. Durbin, P. A. (1980). A stochastic model of two-particle dispersion and concentration fluctuations in homogeneous turbulence. J. Fluid Mech., 100(2):279–302. Gréa, B. J., Burlot, A., Griffond, J., and Llor, A. (2016). Challenging mix models on transients to self-similarity of unstably stratified homogeneous turbulence. J. Fluids Eng., 138(7):070904. 144 Griffond, J. (2005). Linear interaction analysis applied to a mixture of two perfect gases. Phys. Fluids, 17(8):086101. Hannappel, R. and Friedrich, R. (1995). Direct numerical simulation of a Mach 2 shock interacting with isotropic turbulence. Appl. Sci. Res., 54(3):205–221. Hesselink, L. and Sturtevant, B. (1988). Propagation of weak shocks through a random medium. J. Fluid Mech., 196:513–553. Huete, C., Jin, T., Martínez-Ruiz, D., and Luo, K. (2017). Interaction of a planar reacting shock wave with an isotropic turbulent vorticity field. Phys. Rev. E, 96(5):053104. Huete, C., Velikovich, A. L., and Wouchuk, J. G. (2011). Analytical linear theory for the interaction of a planar shock wave with a two-or three-dimensional random isotropic density field. Phys. Rev. E, 83(5):056320. Ishihara, T. and Kaneda, Y. (2002). Relative diffusion of a pair of fluid particles in the inertial subrange of turbulence. Phys. Fluids, 14(11):L69–L72. Jaberi, F. A., Livescu, D., and Madnia, C. K. (2000). Characteristics of chemically reacting compressible homogeneous turbulence. Phys. Fluids, 12(5):1189–1209. Jammalamadaka, A., Li, Z., and Jaberi, F. A. (2014). Numerical investigations of shock wave interactions with a supersonic turbulent boundary layer. Phys. Fluids, 26(5):056101. Jamme, S., Cazalbou, J.-B., Torres, F., and Chassaing, P. (2002). Direct numerical simulation of the interaction between a shock wave and various types of isotropic turbulence. Flow Turbul. Combust., 68(3):227–268. Jin, T., Luo, K., Dai, Q., and Fan, J. (2015). Simulations of cellular detonation interaction with turbulent flows. AIAA J., 54(2):419–433. Johnsen, E., Larsson, J., Bhagatwala, A. V., Cabot, W. H., Moin, P., Olson, B. J., Rawat, P. S., Shankar, S. K., Sjögreen, B., Yee, H. C., et al. (2010). Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves. J. Comput. Phys., 229(4):1213–1237. Kim, J. H., Yoon, Y., Jeung, I. S., Huh, H., and Choi, J. Y. (2003). Numerical study of mixing enhancement by shock waves in model scramjet engine. AIAA J., 41(6):1074–1080. Komori, S., Hunt, J., Kanzaki, T., and Murakami, Y. (1991). The effects of turbulent mixing on the correlation between two species and on concentration fluctuations in non-premixed reacting flows. J. Fluid Mech., 228:629–659. Kovasznay, L. S. G. (1953). Turbulence in supersonic flow. J. Aeronaut. Sci., 20(10):657–674. Larsson, J. (2010). Effect of shock-capturing errors on turbulence statistics. AIAA J., 48(7):1554– 1557. 145 Larsson, J., Bermejo-Moreno, I., and Lele, S. K. (2013). Reynolds- and Mach-number effects in canonical shock-turbulence interaction. J. Fluid Mech., 717:293–321. Larsson, J. and Lele, S. K. (2009). Direct numerical simulation of canonical shock/turbulence interaction. Phys. Fluids, 21(12):126101. Lee, S., Lele, S. K., and Moin, P. (1993). Direct numerical simulation of isotropic turbulence interacting with a weak shock wave. J. Fluid Mech., 251:533–562. Lee, S., Lele, S. K., and Moin, P. (1997). Interaction of isotropic turbulence with shock waves: effect of shock strength. J. Fluid Mech., 340:225–247. Lele, S. K. (1992). Compact finite difference schemes with spectral-like resolution. J. Com- put. Phys., 103(1):16–42. Li, Z. and Jaberi, F. A. (2012). A high-order finite difference method for numerical simulations of supersonic turbulent flows. Int. J. Numer. Meth. Fl., 68(6):740–766. Livescu, D., Jaberi, F. A., and Madnia, C. K. (2000). Passive-scalar wake behind a line source in grid turbulence. J. Fluid Mech., 416:117–149. Livescu, D., Jaberi, F. A., and Madnia, C. K. (2002). The effects of heat release on the energy exchange in reacting turbulent shear flow. J. Fluid Mech., 450:35–66. Livescu, D. and Ristorcelli, J. R. (2008). Variable-density mixing in buoyancy-driven turbulence. J. Fluid Mech., 605:145–180. Livescu, D. and Ristorcelli, J. R. (2009). Mixing asymmetry in variable density turbulence. In Eckhardt, B., editor, Adv. Turbul. XII, volume 132 of Springer Proceedings in Physics, pages 545–548. Springer. Livescu, D., Ristorcelli, J. R., Gore, R. A., Dean, S. H., Cabot, W. H., and Cook, A. W. (2009). High-Reynolds number Rayleigh–Taylor turbulence. J. Turbul., 10(13):1–32. Livescu, D., Ristorcelli, J. R., Petersen, M. R., and Gore, R. A. (2010). New phenomena in variable-density Rayleigh–Taylor turbulence. Phys. Scr., T142:014015. Livescu, D. and Ryu, J. (2016). Vorticity dynamics after the shock–turbulence interaction. Shock Waves, 26(3):241–251. Lombardini, M., Hill, D. J., Pullin, D. I., and Meiron, D. I. (2011). Atwood ratio dependence of Richtmyer–Meshkov flows under reshock conditions using large-eddy simulations. J. Fluid Mech., 670:439–480. Mahesh, K., Lee, S., Lele, S. K., and Moin, P. (1995). The interaction of an isotropic field of acoustic waves with a shock wave. J. Fluid Mech., 300:383–407. Mahesh, K., Lele, S. K., and Moin, P. (1997). The influence of entropy fluctuations on the interaction of turbulence with a shock wave. J. Fluid Mech., 334:353–379. 146 McKenzie, J. F. and Westphal, K. O. (1968). Interaction of linear waves with oblique shock waves. Phys. Fluids, 11(11):2350–2362. Meneveau, C. (2011). Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows. Annu. Rev. Fluid Mech., 43:219–245. Menon, S. (1989). Shock-wave-induced mixing enhancement in scramjet combustors. AIAA Paper 89-0104. Moore, F. K. (1953). Unsteady oblique interaction of a shock wave with a plane disturbance. NACA TN-2879. Morkovin, M. V. (1962). Effects of compressibility on turbulent flows. Mécanique de la Turbulence, pages 367–380. CNRS, Paris. In Favre, A., editor, Ooi, A., Martin, J., Soria, J., and Chong, M. S. (1999). A study of the evolution and characteristics of the invariants of the velocity-gradient tensor in isotropic turbulence. J. Fluid Mech., 381:141– 174. Payri, F., Desantes, J. M., and Torregrosa, A. J. (1995). Acoustic boundary condition for unsteady one-dimensional flow calculations. J. Sound Vib., 188(1):85–110. Pirozzoli, S. and Grasso, F. (2004). Direct numerical simulations of isotropic compressible turbu- lence: Influence of compressibility on dynamics and structures. Phys. Fluids, 16(12):4386–4407. Pope, S. B. (1994). Lagrangian pdf methods for turbulent flows. Annu. Rev. Fluid Mech., 26(1):23– 63. Pope, S. B. (2002). Stochastic lagrangian models of velocity in homogeneous turbulent shear flow. Phys. Fluids, 14(5):1696–1702. Quadros, R., Sinha, K., and Larsson, J. (2016). Turbulent energy flux generated by shock/homogeneous-turbulence interaction. J. Fluid Mech., 796:113–157. Ranjan, D., Oakley, J., and Bonazza, R. (2011). Shock-bubble interactions. Annu. Rev. Fluid Mech., 43:117–140. Ribner, H. S. (1954). Convection of a pattern of vorticity through a shock wave. NACA TR-1164. Richardson, L. F. (1926). Atmospheric diffusion shown on a distance-neighbour graph. Proc. R. Soc. (London) Ser. A, 110(756):709–737. Ristorcelli, J. R. and Blaisdell, G. A. (1997). Consistent initial conditions for the dns of compressible turbulence. Phys. Fluids, 9(1):4–6. Ryu, J. and Livescu, D. (2014). Turbulence structure behind the shock in canonical shock-vortical turbulence interaction. J. Fluid Mech., 756:R1. Salazar, J. P. and Collins, L. R. (2009). Two-particle dispersion in isotropic turbulent flows. Annu. Rev. Fluid Mech., 41:405–432. 147 Sawford, B. (2001). Turbulent relative dispersion. Annu. Rev. Fluid Mech., 33(1):289–317. Schwarzkopf, J. D., Livescu, D., Baltzer, J. R., Gore, R. A., and Ristorcelli, J. R. (2016). A two-length scale turbulence model for single-phase multi-fluid mixing. Flow Turbul. Combust., 96(1):1–43. Sethuraman, Y. P. M., Sinha, K., and Larsson, J. (2018). Thermodynamic fluctuations in canonical shock–turbulence interaction: effect of shock strength. Theor. Comp. Fluid Dyn., 32(5):629–654. Shen, P. and Yeung, P. K. (1997). Fluid particle dispersion in homogeneous turbulent shear flow. Phys. Fluids, 9(11):3472–3484. Shu, C. W. and Osher, S. (1989). Efficient implementation of essentially non-oscillatory shock- capturing schemes, ii. In Upwind and High-Resolution Schemes, pages 328–374. Springer. Soukhomlinov, V. S., Kolosov, V. Y., Sheverev, V. A., and Ötügen, M. V. (2002). Formation and propagation of a shock wave in a gas with temperature gradients. J. Fluid Mech., 473:245–264. Squires, K. D. and Eaton, J. K. (1991). Lagrangian and eulerian statistics obtained from direct numerical simulations of homogeneous turbulence. Phys. Fluids A: Fluid Dynamics, 3(1):130– 143. Suresh, A. and Huynh, H. T. (1997). Accurate monotonicity-preserving schemes with Runge-Kutta time stepping. J. Comput. Phys., 136(1):83–99. Thompson, P. A. (1971). Compressible-fluid dynamic. McGraw-Hill. Thomson, D. J. (1990). A stochastic model for the motion of particle pairs in isotropic high-reynolds- number turbulence, and its application to the problem of concentration variance. J. Fluid Mech., 210:113–153. Tian, Y., Jaberi, F., Li, Z., and Livescu, D. (2017a). Numerical study of variable density turbulence interaction with a normal shock wave. J. Fluid Mech., 829:551–588. Tian, Y., Jaberi, F., and Livescu, D. (2017b). Shock propagation in media with non-uniform density. In International Symposium on Shock Waves, pages 1167–1175. Springer. Tian, Y., Jaberi, F., and Livescu, D. (2018). Density effects on the flow structure in multi-fluid shock-turbulence interaction. In 2018 AIAA Aerospace Sciences Meeting, page 0374. Tian, Y., Jaberi, F., Livescu, D., and Li, Z. (2017c). Numerical study of shock–turbulence in- teractions in variable density flows. In Proceedings of TSFP-10 (2017) Chicago, volume 3 of International Symposium on Turbulence and Shear Flow Phenomena, pages 8C–3. Tian, Y., Jaberi, F. A., Livescu, D., and Li, Z. (2017d). Numerical simulation of multi-fluid shock- turbulence interaction. In AIP Conference Proceedings, volume 1793 of Shock Compression of Condensed Matter, page 150010. AIP Publishing. Toschi, F. and Bodenschatz, E. (2009). Lagrangian properties of particles in turbulence. Annu. Rev. Fluid Mech., 41:375–404. 148 Wang, J., Shi, Y., Wang, L. P., Xiao, Z., He, X. T., and Chen, S. (2012). Effect of compressibility on the small-scale structures in isotropic turbulence. J. Fluid Mech., 713:588–631. Wang, L. and Lu, X. Y. (2012). Flow topology in compressible turbulent boundary layer. J. Fluid Mech., 703:255–278. Whitham, G. B. (1958). On the propagation of shock waves through regions of non-uniform area or flow. J. Fluid Mech., 4(4):337–360. Yang, J., Kubota, T., and Zukoski, E. E. (1993). Applications of shock-induced mixing to supersonic combustion. AIAA J., 31(5):854–862. Yeung, P. K. (1994). Direct numerical simulation of two-particle relative diffusion in isotropic turbulence. Phys. Fluids, 6(10):3416–3428. Yeung, P. K. (1997). One-and two-particle lagrangian acceleration correlations in numerically simulated homogeneous turbulence. Phys. Fluids, 9(10):2981–2990. Yeung, P. K. (2002). Lagrangian investigations of turbulence. Annu. Rev. Fluid Mech., 34(1):115– 142. Yeung, P. K. and Borgas, M. S. (2004). Relative dispersion in isotropic turbulence. part 1. direct numerical simulations and reynolds-number dependence. J. Fluid Mech., 503:93–124. Yeung, P. K. and Pope, S. B. (1988). An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence. J. Comput. Phys., 79(2):373–416. Yeung, P. K. and Pope, S. B. (1989). Lagrangian statistics from direct numerical simulations of isotropic turbulence. J. Fluid Mech., 207:531–586. Zabusky, N. J. (1999). Vortex paradigm for accelerated inhomogeneous flows: Visiometrics for the Rayleigh-Taylor and Richtmyer-Meshkov environments. Annu. Rev. Fluid Mech., 31(1):495–536. 149