Simulation-Based Study of Particle Dynamics in Homogeneous Isotropic Turbulence and Open Channel Flows A DISSERTATION SUBMITTED TO THE FACULTY OF THE UNIVERSITY OF MINNESOTA BY Yuanqing Liu IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY Lian Shen August, 2024 © Yuanqing Liu 2024 Acknowledgements This thesis could not have been completed without the invaluable support of numerous individuals. First and foremost, I am deeply grateful to my advisor, Prof. Lian Shen, whose guidance was crucial throughout my Ph.D. studies. Prof. Shen introduced me to the field of computational fluid dynamics, allowing me the freedom to choose my research topics and encouraging me continuously. His broad knowledge and connections greatly facilitatedmy networkingwith experts inmyfield. Fromhim, I learned essential academic skills such as paperwriting, generating research ideas, and crafting proposals. Prof. Shen’s exemplary presentation skills and enthusiasm for his work have significantly influenced and motivated me. I extend my gratitude to Prof. Filippo Coletti, with whom I collaborated during the early years of my Ph.D. on the research topic of tracking coherent clusters. Our discussions, both at the university and during APS meetings after he left U.S., have always been inspiring and encouraging. Special thanks to Prof. Chris Paola, who introduced me to the study of sediment trans- port over vibrating beds. His rigorous approach to research and constructive discussions have been extremely helpful. He also encouraged me to value my work and maintain a positive outlook on my achievements. i I am thankful for Prof. Jiarong Hong, Prof. Sungyon Lee, and Prof. Chris Paola for serving on my final committee and for their time spent reviewing my thesis. I am fortunate to have collaborated with talented colleagues in our research group. Dr. Anqing Xuan, in particular, provided insightful discussions that greatly enhanced my research. I am also grateful to Dr. Zixuan Yang, Dr. Han Liu, Dr. Sida He, Dr. Qiang Gao, and Dr. Yadong Zeng for their assistance throughout my Ph.D. studies. Thanks also go to the other members of the Fluid Mechanics Lab for their friendship and support both in research and in life. Finally, I cannot express enough gratitude to my family, whose unwavering support has been my cornerstone throughout my Ph.D. journey. ii Dedication This thesis is dedicated to my family and friends who encouraged me to pursue the doctoral degree. iii Abstract Particle-laden turbulent flows are fundamental in both natural environments and industrial applications. The dynamics of the particles are governed by the particle properties themselves and the fluid flow field. In this work, we studied particle dynamics in three different scenarios: (i) particle clusters in homogeneous isotropic turbulence; (ii) sediment transport of particleswith different shapes; (iii) sediment transport over vibrating beds. In the study of particle clusters in homogeneous isotropic turbulence, we proposed a novel method to study the temporal behaviors of particle clusters. The results show that large particle clusters tend to have longer lifetimes. During the span of a particle cluster’s lifetime, the size of the cluster, enstrophy, and turbulent intensity experienced by the particles remain constant. In the study of sediment transport of particles with different shapes, we performed particle-resolved simulations of spherical particles, oblate particles, and the mixture of these two particles. Our results show that the mean velocity profiles are highly related to the solid volume fraction. Spherical particles prefer rolling motion in the sediment bed, while oblate particles tend to slide, and this trend persists even in the mixed particle case. Preferential orientations for oblate particles have also been observed. In the study of sediment transport over vibrating beds, varying oscillation conditions in horizontal and/or vertical directions are considered. Time-averaged results show that the mean velocity profile and Reynolds normal stresses are decreased in cases with oscillations. Phase-averaged results suggest that fluid velocities and Reynolds normal stresses are significantly influenced by the oscillation phases. Comparisons with the iv pulsating turbulent channel flows are made, and similarities and differences are discussed. Phase-averaged sediment transport rates are reported, and different sediment transport patterns are observed, underscoring the critical role of phase differences in sediment transport processes. Overall, this thesis presents significant advancements in the field of particle-laden flow dynamics, providing a deeper understanding of the complex interactions at play, which are crucial for improving the modeling of environmental processes and engineering systems. v Contents Acknowledgements i Dedication iii Abstract iv List of Tables ix List of Figures x 1 Introduction 1 1.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Numerical method 7 2.1 Immersed boundary method for simulation of particles . . . . . . . . . 7 2.2 Collision model for particles . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Validation cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4 Particle cluster identification . . . . . . . . . . . . . . . . . . . . . . . 24 vi 3 Life and death of inertial particle clusters in turbulence 28 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Numerical Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2.1 Identification and tracking of clusters . . . . . . . . . . . . . . 32 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3.1 Lifetime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3.2 Birth and death . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3.3 Relative dispersion of clustered particles . . . . . . . . . . . . . 38 3.3.4 The role of the turbulent activity . . . . . . . . . . . . . . . . . 40 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4 Interface resolved simulation of particles with different shapes in sediment transport 45 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2.2 Simulation set-up . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2.3 Sediment height and pressure gradient . . . . . . . . . . . . . . 50 4.2.4 Simulation initialization . . . . . . . . . . . . . . . . . . . . . 52 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3.1 Fluid flow properties . . . . . . . . . . . . . . . . . . . . . . . 56 4.3.2 Particle properties . . . . . . . . . . . . . . . . . . . . . . . . 59 4.3.3 Shear stress decomposition . . . . . . . . . . . . . . . . . . . . 66 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 vii 5 Particle-resolved simulation of sediment transport over vibrating beds 70 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2.1 Reference frame . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2.2 Governing equation and computational frameworks . . . . . . . 74 5.2.3 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2.4 Flow property decomposition, sediment height, and Reynolds number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3.1 Basic statistics of the flow and particles . . . . . . . . . . . . . 81 5.3.2 Comparison with pulsating turbulent flows . . . . . . . . . . . 94 5.3.3 Sediment transport results . . . . . . . . . . . . . . . . . . . . 101 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6 Conclusing Remarks 114 6.1 Contribution of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.2 Future studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.2.1 GPU-accelerated simulations . . . . . . . . . . . . . . . . . . . 115 6.2.2 Expanding particle diversity in simulations . . . . . . . . . . . 116 6.2.3 Sediment dynamics under varied oscillation conditions . . . . . 117 Bibliography 118 viii List of Tables 2.1 Parameters for single spherical particle sedimentation in quiescent fluid. 20 2.2 Parameters for particle-wall collision case. . . . . . . . . . . . . . . . . 22 3.1 Non-dimensional parameters characterizing the different simulations, and the characteristic lifetime (expressed in Kolmogorov time-scales) for each case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1 Simulation cases and parameters, #( is the number of moving spherical particles and #$ is the number of moving oblate particles. . . . . . . . 50 5.1 Simulation cases and parameters. . . . . . . . . . . . . . . . . . . . . . 78 ix List of Figures 2.1 Example of Lagrangian marker points on a particle surface. . . . . . . . 9 2.2 Particle collision notations. . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Sketch of collision between spheroid particles. . . . . . . . . . . . . . . 18 2.4 Particle settling velocity versus time for case 1. . . . . . . . . . . . . . 21 2.5 Comparison of particle trajectory between simulation and experiment for normal particle-wall collision: (a) (C = 27, (b) (C = 152. . . . . . . . . 22 2.6 Definition on spheroid particles, (a) shows the definition of polar and equatorial radius, (b) defines the angles \ and q: \ is the angle between the symmetry axis and the z-axis, and q is the angle between the x-axis and the projection of the symmetry axis on the x-y plane. . . . . . . . . 23 2.7 Angular velocity of spheroid particle with aspect ratio 0/1 = 2.5 in simple shear flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.8 Example of two-dimensional Voronoï diagram. . . . . . . . . . . . . . 25 2.9 An example of PDF of Voronoï cell volume normalized by the average cell volume compared with Γ distribution. . . . . . . . . . . . . . . . . 25 2.10 Cluster surface area versus cubic root of cluster volume. . . . . . . . . . 26 x 3.1 (a) Schematic illustration of the method of recognizing connections be- tween clusters. The numerical labels identify the same particles through two successive time steps. (b) Forward connection distribution of case (C = 4, ({ = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 PDF of the cluster lifetimes for the different cases. . . . . . . . . . . . . 36 3.3 (a) PDF of cluster size conditioned on the range of cluster lifetime for case St = 4, Sv =37. (b) Cluster size versus lifetime for the different cases. 37 3.4 (a) Relative probability of the combinations of different types of cluster birth and death for the case St=1, Sv =1.2. For each combination, the average cluster lifetime (b) and size (c) are normalized by their respective maximum values in the population. . . . . . . . . . . . . . . . . . . . 38 3.5 Evolution of the radius of gyration '6 of the particle clouds normalized by its average during C;8 5 4: (a) before birth; (b) during lifetime; and (c) after death, with lines indicating exponential best-fits. In panel (b), the temporal abscissa is normalized by the average cluster lifetime for each case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.6 Evolution of the enstrophy experienced by the particle clouds normalized by its average during C;8 5 4: (a) before birth; (b) during lifetime; and (c) after death. In panel (b), the temporal abscissa is normalized by the average cluster lifetime for each case. . . . . . . . . . . . . . . . . . . . 41 4.1 Visualization for (a) fixed layer of spherical particles and (b) initial par- ticle arrangement for case mix. . . . . . . . . . . . . . . . . . . . . . . 53 4.2 Mean streamwise velocity profiles for three cases. The markers indicate the location of Sediment bed height �B43 . . . . . . . . . . . . . . . . . 55 xi 4.3 Mean solid fraction profiles for three cases. . . . . . . . . . . . . . . . 55 4.4 Reynolds normal stresses normalized with the bulk velocity: (a) stream- wise components, (b) wall-normal components, and (c) spanwise com- ponents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.5 Solid volume fraction contribution from different particles in case mix: (a) spherical and oblate particle solid fraction, (b) ratio of spherical to oblate particle solid fraction. Red dashed line indicates the sediment bed height. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.6 Mean particle velocities comparison: (a) comparison between simulated cases, (b) comparison within case mix. The figure only shows the region where the particles are present. . . . . . . . . . . . . . . . . . . . . . . 60 4.7 Spanwise angular velocity comparison: (a) comparison between simu- lated cases, (b) comparison within case mix. The figure only shows the region where the particles are present. . . . . . . . . . . . . . . . . . . 60 4.8 Spanwise angular velocity comparison for spherical and oblate particles. 61 4.9 Orientation of oblate particles: (a) angle between the major axis and wall normal direction, (b) orientation of the major axis. . . . . . . . . . . . . 61 4.10 Snapshots of the sediment bed for case (a) mix, (b) oblate, and (c) sphere. 65 4.11 Spatial-temporal evolution of sediment bed interface in case (a) mix, (b) oblate, and (c) sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.12 Shear stress and its different components: (a) total shear stress, (b) viscous shear stress, (c) Reynolds shear stress, and (d) particle-induced shear stress. 66 5.1 Sketch of laboratory experiment setup. . . . . . . . . . . . . . . . . . . 71 5.2 Profiles of (a) solid fraction, (b) mean flow velocity obtained for each case. 81 xii 5.3 Snapshot of the waveform of sediment bed obtained in case H+V. . . . . 82 5.4 Profiles of (a) time-averaged Reynolds normal stress D′D′ (—), {′{′ (- -),|′|′ (– · –), and (b) time-averaged Reynolds shear stress −D′{′. . . 82 5.5 Profiles of particle-induced stress g? for the first four cases. . . . . . . . 83 5.6 Profiles of (a) mean particle velocity and (b) mean particle angular ve- locity in the spanwise direction for case N (The red dashed line indicates H = �B43), also here we only show the regions where the particle statistics are available. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.7 Profiles of the phase-averaged (a) fluid velocity and (b) solid volume fraction for case V. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.8 Profiles of phase-averaged r.m.s. values of velocity fluctuations for case V D ′ D ′ (—), {′{′ (- -),|′|′ (– · –). . . . . . . . . . . . . . . . . . . . . . 87 5.9 Profiles of phase-averaged particle velocity for case V. The dashed lines indicate the time-averaged particle velocity. . . . . . . . . . . . . . . . 88 5.10 Profiles of phase-averaged (a) fluid velocity, (b) particle velocity, (c) relative velocity of the fluid with respect to the particles for case H. The dashed lines indicate time-averaged values. . . . . . . . . . . . . . . . 89 5.11 Profiles of phase-averaged r.m.s. values of velocity fluctuations for case H 91 5.12 Profiles of phase-averaged (a) fluid velocity, (b) particle velocity, (c) relative velocity of the fluid with respect to the particles for case H-V. The dashed lines indicate time-averaged values. . . . . . . . . . . . . . 91 5.13 Profiles of phase-averaged (a) fluid velocity, (b) particle velocity, (c) relative velocity of the fluid with respect to the particles for case H+V. The dashed lines indicate time-averaged values. . . . . . . . . . . . . . 92 xiii 5.14 Profiles of phase-averaged solid volume fraction for case H+V. . . . . . 92 5.15 Comparison of velocity amplitude and phase lag of the oscillating velocity with the Stokes solution for case H (a, b) and case H-V (c, d). . . . . . . 97 5.16 Comparison of velocity amplitude and phase lag of the oscillating velocity with the Stokes solution for case H+V (a, b) . . . . . . . . . . . . . . . 98 5.17 Phase-averaged turbulence intensity modulation in case H. The red lines indicate the sediment bed height. . . . . . . . . . . . . . . . . . . . . . 100 5.18 Phase-averaged turbulence intensity modulation in case H-V. The red lines indicate the sediment bed height. . . . . . . . . . . . . . . . . . . 100 5.19 Time evolution of sediment transport rate in case N. . . . . . . . . . . . 102 5.20 Time evolution of sediment transport rate in, (a) case V, (b) case H, (c) case H-V, (d) case H+V. . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.21 Spatial–temporal evolution of sediment bed interface in (a) case N, (b) case V, (c) case H, (d) case H-V. . . . . . . . . . . . . . . . . . . . . . 105 5.22 Spatial–temporal evolution of sediment bed interface in case H+V. . . . 106 5.23 Time evolution of a the r.m.s. of sediment bed height fluctuations and b the mean wavelength _ℎ in case H+V. . . . . . . . . . . . . . . . . . . 106 5.24 Time evolution of sediment transport rate in, (a) (c) (e) (g) horizontal directions, (b) (d) (f ) (h) vertical directions. . . . . . . . . . . . . . . . 109 xiv Chapter 1 Introduction 1.1 Motivations Particle-laden turbulent flows are ubiquitous in natural and engineering settings, com- monly seen in river flow, fluidized beds, and dust storms. Fascinating behaviors of the particles have been reported, including preferential concentration of particles (Eaton & Fessler, 1994) and preferential alignment of non-spherical particles (Voth & Soldati, 2017). Rapid development of experimental and numerical tools has been made in recent years to study different aspects of particle dynamics in turbulent flows (Brandt & Coletti, 2022). Our study on particle dynamics is motivated by the various unresolved scien- tific questions across different flow conditions, specifically, the clustering of particles in homogeneous isotropic turbulence and sediment transport in open channel flow. Particle clustering, the phenomenon where particles aggregate to form clusters, has been a subject of research for decades (Eaton & Fessler, 1994; Monchaux et al., 2010). Both numerical and simulation works have shown that the cluster is most pronounced 1 when the Stokes number of the particles, (C, is around 1 (Aliseda et al., 2002; Eaton & Fessler, 1994; Maxey, 1987). While individual particles may have negligible impact on the turbulent flow, significant accumulations can alter the flow dynamics considerably (Elghobashi, 1994), and in turn, the dynamics of the particles can be affected (Ardekani et al., 2017). Numerous efforts have been made to characterize and identify these clusters (Aliseda et al., 2002; Monchaux et al., 2010; Baker et al., 2017). Despite these advances, the temporal evolution of clusters remains poorly understood. By studying the temporal evolution of the clusters, we can better understand the formation and destruction of the clusters and the role of turbulence during the process, which would be helpful for future modeling of the particle-laden flow. Sediment transport is commonly seen in natural settings, with themost frequent exam- ples being the transport of sand on land, in rivers, and under seas. Sediment transport is a complex phenomenon, with interactions between the particles and the fluid and massive collision events between the particles themselves. Factors such as particle electrification also impact sediment transport, as shown in studies on wind-blown sand (Zheng et al., 2003). Particle motion in sediment transport is generally categorized into rolling, sliding, and saltating movements (Dey & Ali, 2017). The lifted particles in saltating movements will eventually deposit on the sediment bed, and under certain conditions, this behavior can amplify small perturbations on the flat sediment bed, leading to the formation of ripples and dunes (Kidanemariam & Uhlmann, 2014a). Various numerical tools have been used to study sediment transport, including Eulerian descriptions (Zedler & Street, 2006; Yang et al., 2017), Lagrangian descriptions (Schmeeckle, 2014; Sun &Xiao, 2016; Zheng et al., 2021), and particle-resolved simulations (Kidanemariam&Uhlmann, 2017; Jain et al., 2019). Particle-resolved simulations offer the most detailed view by directly 2 resolving particle geometries and interactions without relying on empirical formulas, and have gained popularity with the rapid increase in computational power. Particle-resolved simulations have shed light on phenomena such as saltation (Ji et al., 2014; Zhu et al., 2022), the generation of bedform patterns (Kidanemariam & Uhlmann, 2014a, 2017), and sediment transport under waves (Mazzuoli et al., 2019; Vittori et al., 2020). However, there are still many unknowns in sediment transport. Two topics have caught our attention. The first one is the sediment transport of non-spherical particles. Non-spherical particles have exhibited complex behaviors in turbulent flows (Voermans et al., 2017). In the context of sediment transport, Allen (2012) observed that spherical particles tend to roll, while disc-like particles tend to slide, and Krumbein (1942) found the largest traveling velocities for spherical particles among investigated shapes, with disc-like particles having the smallest velocity. Unfortunately, current studies using particle-resolved simulation on sediment transport primarily consider spherical particles. Limited literature can be found for particle-resolved simulations of non-spherical particles in sediment transport, and there is an urgent need to fill the gap in our understanding of sediment transport involving non-spherical particles, enriching the dataset for particle- resolved simulations and offering new insights into their behaviors. The second topic is related to surface-sediment remobilization during seismic events. Earthquakes are a major mechanism for the initiation of deep-sea turbidity currents, and seismic history is believed to be relevant to underwater avalanches (Adams, 1990; Goldfinger et al., 2007; Ashi et al., 2014). Little is known about the remobilization of surface particles during seismic events. As mentioned by (Atwater et al., 2014), even earthquakes of identical magnitude in the same locale can produce distinct turbidity signatures. We believe that performing the first particle-resolved simulations on this 3 topic can provide valuable insights into the interplay between fluid and particles under seismic events and can inform future studies and applications in civil engineering and environmental management. 1.2 Thesis overview This study aims to provide in-depth insights into particle dynamics within homogeneous isotropic turbulence and open channel flow, exploring the interplay between particles and fluid dynamics under various conditions. In chapter 2, we delineate the numerical method employed in our research. The simulations on open channel flow are based on particle-resolved direct numerical simula- tions (PR-DNS), utilizing the immersed boundary method proposed by Uhlmann (2005); Breugem (2012). Particle-particle and particle-wall collisions are modeled using the soft sphere model introduced by Costa et al. (2015); Ardekani et al. (2016). We present detailed descriptions of the numerical models, supplemented by canonical validation cases that confirm the correct implementation of our numerical tools. Additionally, we describe the particle cluster identification method, which is pivotal for our analysis of particle clusters. In chapter 3, we delve into the study of particle clusters within homogeneous isotropic flow. While previous research has primarily focused on the transient characteristics of particle clusters, our study introduces a novel methodology for tracking clusters over time, enabling an analysis of their temporal behaviors. We analyze a broad range of particle data, covering various Stokes numbers and Froude numbers, observing consistent characteristics of cluster dynamics. Our findings suggest that the exponential distribution 4 provides a reasonable approximation of particle lifetime probabilities, and we found that larger clusters tend to persist longer within the turbulence. The study shows that large, long-lived clusters are usually generated by the recombination of other substantial clusters and are likely to recombine into new, large clusters. The size of the cluster during its life span is relatively constant. The study also examines the role of turbulence in the clustering process, noting that enstrophy and turbulent intensity tend to remain constant during the lifespan of the particles. In chapter 4, we focus on particle-resolved simulations of sediment transport in open channel flow, specifically examining spherical and oblate particles with an aspect ratio of 0.5. We maintain a constant flow rate across all simulation cases. Our analysis indicates that mean velocity profiles correlate strongly with the solid volume fraction, with a larger sediment bed height reported for cases with only spherical particles. Notably, cases involving oblate particles exhibit higher normal and spanwise Reynolds normal stresses. The dynamics within the sediment bed show that spherical particles typically exhibit rolling motions, while oblate particles tend to slide—a trend that persists even in the mixed particle case. Preferential orientations for oblate particles are observed, with similar mean orientation angles maintained across various cases. Additionally, the sediment transport rates calculated in this chapter show good agreement with empirical formulas. The analysis of shear stress decomposition reveals that particle-induced shear stress is highly correlated with solid volume fraction, as well as particle shape. In Chapter 5, we investigate sediment transport over a vibrating bed, in an effort to study the remobilization of particles during seismic events. We consider varying os- cillation conditions in horizontal and/or vertical directions, with time-averaged analyses 5 revealing smaller mean velocities and reduced Reynolds normal stresses for oscillation- involved cases. Phase-averaged results suggest that fluid velocities and Reynolds normal stresses are significantly influenced by the oscillation phase, except in cases with only vertical oscillations. Our simulations with horizontal oscillations fall into the current- dominated, very-high-frequency regime for pulsating channel flow, showing good agree- ment with the Stokes solution in these cases, while phase lags continue to increase within the sediment bed. For turbulent intensity modulations, we observe collective changes of turbulent intensity in all three directions. The phase-averaged sediment transport rates show that the highest sediment rates typically coincide with near-zero vertical sediment transport rates. Comparative analyses between cases with the same vertical and horizon- tal oscillations, but with a phase difference of c, exhibit varying transport rate patterns, highlighting the critical role of phase differences in sediment transport processes. 6 Chapter 2 Numerical method 2.1 Immersed boundary method for simulation of parti- cles The Immersed Boundary (IB) method, first proposed by Peskin (1982), has been de- veloped over the years. For simulation of particles, the numerical algorithm have been modified by many researchers (Mohd-Yusof, 1996; Uhlmann, 2008; Kim & Choi, 2006; Lucci et al., 2010; Breugem, 2012; Tschisgale et al., 2018). For our simulation study, we primarily follow the IB method proposed by Uhlmann (2005) and Breugem (2012). This IB method consists of two grids, a fixed uniform grid for the fluid phase and a uniform Lagrangian grid attached to and moving with the particles. The Navier–Stokes equations for fluid phase are: ∇ · u = 0, (2.1) 7 mu mC + u · ∇u = − 1 d 5 ∇? + ` d 5 ∇2u + f, (2.2) here, u is the fluid velocity, ? is the pressure, d 5 is the fluid density, ` is the dynamic viscosity, and f is a volume force term representing the force from particles to the fluid field. For each particle, the governing equations are: d?+? 3u? 3C = d 5 ∮ m+ (3 · n p)3� + (d? − d 5 )+?g + Lc, (2.3) O? 38? 3C = ∮ m+ r × (3 · n p)3� + Zc, (2.4) here, u? is the velocity of the particle centroid, and 8? is the angular velocity of the particle. d? and +? are the density and volume of the particle, respectively. g = −?O + ` 5 (∇u + ∇u) ) is the stress tensor, n is the outward-pointing normal vector of the particle surface, g is the gravitational acceleration, Lc and Zc are the contact force and torque, respectively. O? is the moment of inertia of the particle, and r is the position vector from the centroid to the point on the particle surface. The calculation of collision force and torque will be introduced in § 2.2. Each particle is represented by numerous uniformly distributed points on the particle surface, known as Lagrangian marker points. As shown in figure 2.1, there points are uniformly distributed on a spheroid. These marker points form a shell around the particle with a thickness equivalent to the grid spacing, ΔG, as described by Uhlmann (2005). Each point represents a volume, Δ+ , ideally equivalent to the volume of a single grid cell in the uniform Eulerian mesh, ΔG3. This equivalency governs the number of Lagrangian points on the particle surface. 8 Figure 2.1: Example of Lagrangian marker points on a particle surface. For example, for a spherical particle with radius r, imagine a shell on the particle surface with radii A − ΔG/2 and A + ΔG/2. The shell volume can be calculated as: +Bℎ4;; = cΔG 3 (12A2 + ΔG2), (2.5) and to make the Lagrangian point cell volume equal to the single grid cell volume, the number of Lagrangian points on the particle surface can be calculated as: # ≈ c 3 (12 A2 ΔG2 + 1), (2.6) To uniformly distribute N particles on a sphere with radius r, Saff & Kuijlaars (1997) proposed the following method: ℎ: = −1 + 2(: − 1) # − 1 , : = 1, 2, ..., #, (2.7) 9 \: = arccos(ℎ: ), : = 1, 2, ..., #, (2.8) q: = (q:−1 + 3.6 √ # 1√ 1 − ℎ2 : ) (<>32c), : = 2, 3, ..., # − 1, q1 = q# = 0. (2.9) and the Lagrangian point locations can be calculated as: G: = A sin(\: ) cos(q: ), (2.10) H: = A sin(\: ) sin(q: ), (2.11) I: = A cos(\: ). (2.12) For non-spherical particles, there is no straightforward method to directly calcu- late particle distributions. To achieve uniformly distributed Lagrangian points on non- spherical particle surface, we introduce the Riesz s-energy, which is defined as: �B = ∑ 8≠ 9 1 |x8 − x 9 |B , (2.13) Hardin & Saff (2005) proved that for a 3D geometry, when B ≥ 2, the minimum s-energy indicates the uniform distribution of points on the geometry surface. Thus, for Lagrangian points on a non-spherical particle surface, we first calculate the number of Lagrangian points on the surface. Then, we randomly generate points on the particle surface. These points move under the influence of the Riesz s-energy until the total displacement of all Lagrangian points is smaller than a threshold. Finally, these points represent the uniformly distributed points on the particle surface. Having obtained the Lagrangian grid points on the particle surface, we proceed to calculate the Immersed Boundary force between the fluid and the particle. Prior to this, it 10 is essential to revisit the solution process for the Navier–Stokes equations on the Eulerian grid. The fractional-step method (Kim & Moin, 1985) is commonly employed for this purpose, which involves the following steps: 1. Predict the intermediate velocity field: u∗ = u= + AℎB=ΔC, (2.14) 2. Solve the Poisson equation for pressure: ∇2?=+1 = d 5 ΔC ∇ · u∗, (2.15) 3. Correct the velocity field: u=+1 = u∗ − ΔC d 5 ∇?=+1. (2.16) in the first step, ΔC is the time step, AℎB= is the right-hand side of the momentum equation AℎB = −u · ∇u + 1 d 5 ∇? − ` d 5 ∇2u, and = is the time step index. For the IB method, the fractional-step method can still be utilized, with the primary difference being the inclusion of the IB force, resulting in a second prediction step for the velocity field: u∗∗ = u∗ + fΔC, (2.17) where f is the force from the particles to the fluid field. Then step 2 and 3 of the fractional-step method are performed, with the change of u∗ to u∗∗. 11 For any point X at a particle surface, the velocity U at this point can be expressed as: U = up + 8 p × (X − X?), (2.18) where X? is the location of the particle centroid. The boundary condition at particle surface is a non-slip and non-penetrable condition, u = U(X). The difference between the fluid velocity and the particle surface velocity is the source for the IB force. The IB force can be calculated by: 1. Interpolating the first prediction velocity (u∗) to the Lagrangian marker points: U∗; = ∑ 8 9 : u∗8 9 :X3 (x8 9 : − X;)ΔG3, (2.19) 2. Calculating the force on the Lagrangian marker points: F; = U? (X;) − U∗ ; ΔC , (2.20) 3. Interpolating the force back to the Eulerian grid: f8 9 : = ∑ ; F;X3 (x8 9 : − X;)Δ+, (2.21) here, the uppercase letter represents quantities on the Lagrangian marker points, while lowercase letters represent Eulerian grid quantities. The subscript (8 9 :) denotes the position on the Eulerian grid with index (8, 9 , :). U∗ ; is the interpolated first prediction velocity on the Lagrangian marker points, and X; is the position of the Lagrangian point with index l on the Eulerian grid. X3 is the regularized Dirac delta function, and Δ+ is the 12 volume of the Lagrangian grid cells. Following Uhlmann (2005) and Breugem (2012), we use the regularized Dirac delta function of Roma et al. (1999), which has the form: X= (A) =  1 6 ( 5 − 3|A | − √ −3(1 − |A |)2 + 1 ) , if 0.5 ≤ |A | ≤ 1.5 1 3 ( 1 + √ −3A2 + 1 ) , if |A | ≤ 0.5 0, if |A | > 1.5 (2.22) where A = |x8 9:−X; | ΔG . For the particle governing equation (2.3), there is a surface integral term, ∮ m+ (3 · n p)3�, this term can be written as: ∮ m+ (3 · n p)3� = −d 5 ∫ +? f3+ + d 5 3 3C ( ∫ +? u3+). (2.23) here the first term is related to the sum over all IB forces: − d? ∫ +? f3+ = −d 5 ∑ 8 9 : fΔG3, (2.24) and the second term of equation 2.23 accounts for the change of fluid inertia within the particle, and we can directly evaluate the volume integrals by means of a second-order accurate midpoint rule as proposed by Kempe et al. (2009), in which the momentum integral is computed as: ∫ +? u3+ = ∑ 8 9 : u8 9 :U8 9 :ΔG3. (2.25) here, U8 9 : is the volume fraction of the particle in the Eulerian grid cell (8, 9 , :), and this can be calculated by looking at a level-set function q, which is the signed distance 13 function from the particle surface. The value of q is negative when it is inside the particle and positive otherwise. The solid volume fraction can be evaluated by the 8 corner nodes of the grid cell in the following form: U8 9 : = ∑8 ==1 −q=� (−q=)∑8 ==1 |q= | . (2.26) where � is the Heaviside step function. This approach ensures accurate and efficient computation of the interactions between the fluid and the particle within the Eulerian grid framework. Using the same analysis, the flow-induced torque on the particle can be calculated by: ∮ m+ r × (3 · n p)3� = −d 5 ∫ +? r × f3+ + d 5 3 3C ( ∫ {? r × u3+). (2.27) 2.2 Collision model for particles The simulation of particle interactions in fluid dynamics relies critically on the choice of an appropriate collision model. Three primary models are prevalent in particle-resolved simulations: Repulsive Potential Model: This model is widely applied in studies of dilute flows due to its simplicity (Glowinski et al., 2001; Uhlmann, 2008). However, it does not account for energy dissipation or tangential forces, which limits its usage in collision- intensive systems. Hard Sphere Model: Based on binary, quasi-instantaneous collisions, this model ensures that particles do not overlap and their velocities are updated immediately (Veera- mani et al., 2009; Derksen, 2011; Jain et al., 2019). It faces challenges in scenarios where 14 particles remain in contact for extended periods, such as in sediment transport. Soft SphereModel: This approach extends the collision duration numerically, allow- ing particles to overlap. The collision force is calculated based on the overlap distance and the relative velocity between particles (Kempe & Fröhlich, 2012; Costa et al., 2015; Biegert et al., 2017), making it well-suited for dense particle systems. For our study, we employed the soft sphere model proposed by Costa et al. (2015), where the collision can be seen as a linear spring–dashpot system. For the collision between two spherical particles, labeled 8 and 9 , the force is determined by two primary parameters: the relative velocity at the contact point and the extent of overlap. The relative velocity at the contact point, u8 9 , is defined by the equation: u8 9 = (u8 + '888 × n8 9 ) − (u 9 + ' 98 9 × n 98), (2.28) where u8 and u 9 represent the velocities of the centers of particles 8 and 9 , respectively; 88 and 8 9 are the angular velocities; '8 and ' 9 are the radii; and n8 9 and n 98 are the unit vectors from the center of one particle to the other. These vectors are defined as: n8 9 = X 9 − X8 |X 9 − X8 | , n 98 = −n8 9 , (2.29) with X8 and X 9 denoting the particle center locations. The relative velocity can be decomposed into normal and tangential components, u8 9 ,= and u8 9 ,C , respectively. The overlap between the particles, denoted as %8 9 ,=, is calculated as: %8 9 ,= = ('8 + ' 9 − ‖x 9 − x8‖)n8 9 , (2.30) 15 The normal force F8 9 ,= resulting from the collision is then expressed as: F8 9 ,= = −:=%8 9 ,= − [=u8 9 ,=, (2.31) where := and [= are the normal stiffness and damping coefficients, respectively. These two coefficients ensures the physical characteristics of the collision response. The notations for the particle collision are shown in figure 2.2 In soft sphere model, the collision duration extends over several computational time steps, ΔC. Typically, the collision time, )2, is set as )2 = #ΔC, with N being 10 in most simulations. To accurately simulate the collision dynamics, we define two key coefficients: the normal coefficient of restitution, 4=,3 , and the tangential coefficient of restitution, 4C,3 . These are expressed as: 4=,3 = D>DC,= D8=,= , 4C,3 = D>DC,C D8=,C , (2.32) where the subscripts 8= and >DC denote the velocities before and after the collision, respectively. The stiffness coefficient := and damping coefficient [= are then determined such that after collision, the particles’ velocities conform to the defined coefficients of restitution. These coefficients are calculated as: := = <4 (c2 + log2(4=,3)) (#ΔC)2 , [= = − 2<4 log(4=,3) (#ΔC) , (2.33) here, <4 represents the effective mass given by <4 = ( 1 <8 + 1 < 9 )−1 , where <8 and < 9 are the masses of particles 8 and 9 . 16 Figure 2.2: Particle collision notations. The tangential force during particle collisions can be modeled using a linear spring– dashpot system. However, this force is constrained by the maximum normal force. Therefore, the tangential force, F8 9 ,C is formulated as: F8 9 ,C = min ( ‖ − :C%8 9 ,C − [Cu8 9 ,C ‖, ‖ − `2F8 9 ,=‖ ) t8 9 , (2.34) where %8 9 ,C denotes the tangential displacement, and t8 9 is the unit vector in the tangential direction. In this model, F8 9 ,C is limited by the Coulomb friction coefficient, `2. The tangential stiffness and damping coefficients, :C and [C , are calculated as: :C = <4,C (c2 + log2(4=,C)) (#ΔC)2 , [C = − 2<4,C log(4C,3) (#ΔC) . (2.35) where <4,C is defined as <4,C = 2 7<4. 17 Figure 2.3: Sketch of collision between spheroid particles. For particle–wall collisions, the wall can be seen as a particle with infinite mass and infinite radii, which makes a particle–wall collision equivalent to a particle–particle collision. The parameters for the particle–wall collision can be calculated in the same way. For modeling collisions involving non-spherical particles, we adopt the method pro- posed by Ardekani et al. (2016). They suggest that the collision dynamics of spheroid particles can be effectively approximated by treating them as spherical particles. This approximation uses spheres with equivalent mass and radii corresponding to the local curvature at the contact points between particles. As shown in figure 2.3, the collision between the two spheroids are approximated by the collision between the two dashed spheres. An important thing to notice is that the normal force can also apply torque on the particle, since the normal force does not necessarily pass through the center of the particle. 18 In scenarios where particles are in close proximity, particularly when the gap between two particles or between a particle and a wall is less than one grid width, the IB method fails to predict forces accurately due to insufficient spatial resolution. To mitigate this limitation, a lubrication force model is employed, enhancing force estimation under narrow gap conditions. In Costa et al. (2015), the lubrication force can be written as: ΔF;D1 = −6c`'D8 9 ,= [_(n) − _(nΔG)], (2.36) where the Stokes amplification factor _ is given by: _?? (n) = 1 2n − 9 20 ln(n) − 3 56 n ln n, (2.37) _?| (n) = 1 n − 1 5 ln(n) − 1 21 n ln n, (2.38) for particle–particle collision and particle–wall collision. Here, n = X8 9 ,=/', where ' is the largest particle radius. The lubrication force is only applied when the gap between particles is less than a threshold, nΔG . 2.3 Validation cases Here we provide several numerical cases to validate our implementation of the Immersed Boundary method and the collision model. 19 Case d?/d 5 6 a 5 '4? Experiment '4? 1 2.56 9.81 1.04 × 10−3 377 362.70 2 7.71 9.81 2.68 × 10−3 283 280.42 Table 2.1: Parameters for single spherical particle sedimentation in quiescent fluid. Case 1: Sedimentation of a spherical particle in a quiescent flow To validate our IB method for a moving particle, we perform the simulation of a spherical particle sedimenting in a quiescent fluid with different density ratios. The original experiments are performed byMordant&Pinton (2000), and it is often used as a validation case for IB method (Uhlmann, 2005). The computational domain is !G × !H × !I = (8 × 100 × 8)�, with � = 20ΔG (3ΔG is the grid spacing) in current simulation, and a grid of #G × #H × #I = 160 × 2000 × 160 is used. The particle is initially at G = !G/2, H = 0.95!H, I = !I/2 and released under gravity from rest at t=0. Two different cases are considered, and the parameters are reported in table 2.1. The particle quickly reaches a terminal velocity, and the particle Reynolds number can be calculated as: '4? = 2' |D? | a 5 , where D? is the particle velocity. The calculated '4? and '4? from the experiments are also reported in table 2.1, where a good agreement is observed. Furthermore, figure 2.4 shows the comparison of particle settling velocity with experiment results, where a good agreement is observed. Case 2: Normal particle-wall collision To validate the collision model, we simulate the bouncing motion of a single spherical particle in a viscous fluid. The experiment is performed by Gondret et al. (2002). Here, two cases with different Stokes number, (C, are considered. The same simulation has been 20 Figure 2.4: Particle settling velocity versus time for case 1. used to validate collision models in several works (Costa et al., 2015; Biegert et al., 2017). The experiment parameters of the two cases are listed in table 2.2. In this simulation, the particle is resolved by � = 20ΔG, and the coefficient of restitution is set to 0.97, the same as the particle properties in the experiment. Periodic boundary conditions are applied in the streamwise and spanwise directions, and a no-slip boundary condition is applied at the top and bottom surface. Initially, the particle is place at G = !G/2, H = !H − 1.5�, I = !I/2 with zero velocity. The particle falls under gravity and reaches its terminal velocity, and then it will collide with the wall. X= denotes the distance between the particle and wall, and it is normalized by the particle diameter�. The particle trajectories are shown in figure 2.5, with the red squares showing the experiment results. As we can see, our simulation results are in good agreement with the experiment. 21 Figure 2.5: Comparison of particle trajectory between simulation and experiment for normal particle-wall collision: (a) (C = 27, (b) (C = 152. St '4? � ? (<) d?/d 5 a 5 (<2/B) 6( !2 are considered as coherent clusters. The coherent clusters are our focus in the study in chapter 3. 27 Chapter 3 Life and death of inertial particle clusters in turbulence 3.1 Introduction The tendency of particles to aggregate while interacting with fluid flows, termed clus- tering, is evident in countless situations, including atmospheric clouds, fluidized bed reactors, swimming microorganisms, just to name a few. Here we focus on the classic case of small heavy particles in turbulence, which are known to cluster when their re- sponse time is comparable to the time scale of the flow. This phenomenon has attracted the attention of the fluid mechanics community for decades (Eaton & Fessler, 1994; Mon- chaux et al., 2012). This is justified by the potentially important role played by clusters This chapter was published as the paper "Life and death of inertial particle clusters in turbulence" (Yuanqing Liu, Lian Shen, Remi Zamansky, Filippo Colltti). My major contribution is formulating ideas, analyzing data and writing the manuscript. 28 in triggering and enhancing the interaction between particles, as well as their collective backreaction on the carrier fluid. As such, significant effort has gone into quantitatively characterizing clustering, with tools ranging from box-counting (Aliseda et al., 2002), to radial distribution functions (Sundaram & Collins, 1997) and advanced topological descriptors (Calzavarini et al., 2008). A method to characterize clustering that has gained increasing popularity is the Voronoï tessellation. Since its introduction to the field of particle-laden turbulence by Monchaux et al. (2010), it has been widely used for both heavy particles and light bubbles (Tagawa et al., 2012), homogeneous and wall-bounded turbulence (Wang et al., 2019), point-like and finite-size particles (Fiabane et al., 2012), monodisperded and polydis- persed particles (Lian et al., 2019), in one-way and two-way coupled regimes (Monchaux & Dejoan, 2017), mostly in hydrodynamically forced but also in buoyancy-driven tur- bulence (Zamansky et al., 2016) and even in freely sedimenting systems (Uhlmann & Doychev, 2014). An advantage of this technique is that it provides the local concentration associated to each particle (through the inverse of the Voronoï cell size), which in turn allows the identification of individual clusters as sets of neighboring particles satisfying a threshold concentration. This feature was used by Baker et al. (2017) to show that clusters of heavy particles in homogeneous turbulence are self-similar objects with a broad range of sizes, and that they have significantly different fall speeds compared to non-clustered particles. The approach has been used to characterize individual clusters, including their shape, orientation, velocity and acceleration, in homogeneous turbulence (Petersen et al., 2019; Momenifar & Bragg, 2020) and in wall-bounded flows (Fong et al., 2019). While the presence of clusters in particle-laden turbulence is not disputed, their origin and dynamics are still debated. The classic explanation refers to a centrifugingmechanism 29 that pushes particles out of vortex cores and in high-strain regions (Eaton & Fessler, 1994). This has been challenged by alternative and not mutually compatible explanations, notably: the path-history effect (Bragg & Collins, 2014), i.e., particles retaining memory of the flow fluctuations they experienced; and the sweep-stick mechanism (Goto & Vassilicos, 2008), i.e., particles sticking to zero-acceleration points swept and clustered by large-scale motions. Naturally, information on the temporal evolution of clusters appears crucial to reach a predictive understanding of their dynamics. There is, however, a dire scarcity of such information in the literature. Rare exceptions are represented by Tagawa et al. (2012) and Lian et al. (2019), who used the temporal autocorrelation of the Voronoï cell sizes to identify the time scales during which particles remain clustered. Indeed, the Voronoï tessellation method is naturally apt to obtain temporal information on clustering dynamics. This was one of the original motivations of Monchaux et al. (2010), who already showed the feasibility of tracking the cell area along particle trajectories. To date, no study has addressed the temporal evolution of individual clusters, nor investigated systematically their lifecycle. As such, several fundamental questions remain unanswered: What is the lifetime of a cluster? How does it depend on the particle properties, and how is it linked to the cluster size? How do clusters merge with and split from others? Which turbulent events lead to the formation and disruption of clusters? These have deep ramifications for how such objects are modeled and how they may affect the carrier fluid flow. In this study we attempt to answer these questions by following clusters in space and time. We use numerical simulations of homogeneous turbulence laden with small inertial particles, with and without the effect of gravity, and introduce a novel methodology to track clusters across successive realizations of the flow. This is done by adding a temporal dimension to the cluster identification method introduced in 30 Baker et al. (2017), to establish whether a cluster survives over successive time steps. 3.2 Numerical Cases We carry out direct numerical simulations of homogeneous isotropic turbulence in a periodic box of length 2c at a resolution of 5123 grid points. The continuity and momentum equations are solved using a pseudo-spectral method (see Zamansky et al., 2016). Steady state is achieved by imposing a large-scale forcing in Fourier space that results in constant mean dissipation (Carbone et al., 2019), obtaining a Taylor-microscale Reynolds number '4_ = 127.4. The time integration follows the second-order Adams– Bashford method. The flow contains 1 million particles, assumed to be spherical, with a density d% much larger than the fluid density, a diameter �% smaller than the Kolmogorov scale, and a response time g% = d%�2 % /(18`) (where ` is the fluid dynamic viscosity). We consider dilute conditions in which the particles do not affect the flow or each other, and Stokes drag and gravity (when present) are the only forces acting on them. We use Lagrangian tracking to obtain the evolution of the particle positions, where the fluid velocity is evaluated by cubic spline interpolation. The time advancement for the particle transport also uses the second-order Adams–Bashford method. In table 3.1 we list the considered cases, generated by a matrix of two Stokes numbers (C = g%/g[ and three Froude numbers �A = 0[/6, where g[ and 0[ are the Kolmogorov time scale and acceleration, respectively, and 6 is the gravitational acceleration. The cases with �A = ∞ correspond to zero-gravity conditions. ({ = (C/�A compares the still-fluid terminal velocity of the particles to the Kolmogorov velocity D[. 31 3.2.1 Identification and tracking of clusters At each time step, we identify the clusters via criteria defined in Baker et al. (2017). Briefly, from the particle positions we perform Voronoï tessellation of the domain and identify sets of adjacent cells smaller than a threshold. Following Monchaux et al. (2010), the latter is given by the intersection of the probability density function (PDF) of the Voronoï cell volumes with the distribution that one would obtain if the particles were randomly distributed (closely approximated by a Γ curve). Of the resulting clusters, we consider those in a size-range that displays a power-law decay (Baker et al., 2017). This translates to a minimum threshold on the cluster size (calculated as the cubic root of the cluster volume) between 4.8[ and 5.8[, the precise value of which does not affect the conclusions of the study. The average inter-particle distance is around 6[, for which qualitative biases in the characterization of the clusters are not expected (Momenifar & Bragg, 2020). The tracking strategy is inspired by the work of Lozano-Durán& Jiménez (2014), who identified and tracked coherent flow structures in turbulent channel flows and analyzed their splitting and merging. However, the present method is significantly different as it leverages the Lagrangian nature of the simulations. Considering two clusters identified in two consecutive time steps, we take both to be successive realizations of the same cluster if the number of particles they share is above a given threshold. Here the time step refers to the period of g[ at which the data is stored and analyzed. The shared particles across clusters in successive time steps are termed connections. We consider forward-in-time and backward-in-time connections, and apply thresholds on the fraction of connected particles over the total number of particles in each cluster. This is illustrated in figure 3.1(a). Cluster A (identified in time step 1) shares all its particles with cluster C (identified 32 in time step 2), while C shares 2/3 of its particles with A. Therefore, the fractions of forward and backward connections between A and C are 1 and 0.67, respectively. On the other hand, B shares 2/3 of its particles with C, and C shares 1/3 of its particles with B. Thus, the forward and backward connections between B and C are 0.67 and 0.33, respectively. In figure 3.1(b) we display the PDF of the fraction of forward connections for the case (C = 4, ({ = 0, representative also of the other cases. The PDFs of the fraction of forward and backward connections are virtually indistinguishable, which is expected given the nature of the definition and the large number of considered clusters. While the occurrence of clusters sharing 90% or more of their particles is rare, the majority share between 30% and 60% of their particles with clusters in previous/following time steps. Considering the above, we adopt the following definition: two clusters in consecutive time steps are identified as the same cluster when their fractions of backward and forward connections are both above 0.5. This criterion eliminates ambiguities, as no cluster can be recognized in two different clusters in past/future instances. Small variations of the threshold do not qualitatively alter the trends reported below. In the example of figure 3.1, A and C are recognized as the same cluster, which therefore is “alive” in both time steps. The cluster lifetime is defined as the time elapsed between birth (the first instance a cluster is identified) and death (the last time it is recognized). We verify that a time step g[/2 produces the same results, as expected since we consider particle response times equal or larger than g[. 33 (C �A ({ );8 5 4/g[ 1 +∞ 0 3.10 1 0.82 1.2 3.38 1 0.11 9.2 5.11 4 +∞ 0 4.78 4 0.82 4.9 5.01 4 0.11 37 8.71 Table 3.1: Non-dimensional parameters characterizing the different simulations, and the characteristic lifetime (expressed in Kolmogorov time-scales) for each case. 0 0.2 0.4 0.6 0.8 1 fraction of forward connections 0 0.05 0.1 0.15 0.2 P D F (b) 3 5 1 2 8 4 6 7 9 10 11 12 13 14 Time 1 Cluster A Cluster B Cluster C Time 2 10 11 12 13 3 5 1 2 8 4 6 7 (a) Figure 3.1: (a) Schematic illustration of the method of recognizing connections between clusters. The numerical labels identify the same particles through two successive time steps. (b) Forward connection distribution of case (C = 4, ({ = 0. 34 3.3 Results 3.3.1 Lifetime We begin the analysis by considering the distribution of the cluster lifetime C;8 5 4 for the various cases, displayed in figure 3.2. Most clusters live for only a few Kolmogorov time-scales, but the stretched tails of the distributions indicate significant variability. An exponential distribution provides a reasonable approximation of the probability %(C;8 5 4), fromwhich a characteristic lifetime);8 5 4 can be defined, i.e., %(C;8 5 4) ∼ 4G?(−C;8 5 4/);8 5 4). A least-square fit returns the values listed in table 3.1. The values are consistent with the time-scale ∼ 4g[ obtained by Tagawa et al. (2012) from the Lagrangian autocorrelation of the Voronoï volume sizes for inertial particles in homogeneous turbulence. They did not report significant variation of the decorrelation times with increasing St, and did not consider the effect of gravity. Our results indicate that the cluster lifetimes (which depend on the relative motion and position of a large number of particles) do increase with St. Moreover, gravitational settling also contributes to extending the cluster lifetime. One can expect the lifetime to be related to other cluster properties, in particular the size !� . The latter is taken as the cubic root of the cluster volume averaged over the lifetime. Figure 3.3(a) shows the PDF of !� conditioned on C;8 5 4 for the case (C = 4, ({ = 37 (the other cases displaying analogous trends). Clearly, the cluster size and lifetime are positively correlated. The trend is confirmed by figure 3.3(b) that displays the mean !� conditioned on C;8 5 4 for all cases. In general, size and lifetime grow with increasing importance of both inertia (St) and gravity (Sv). The trend with Sv is consistent with the argument of Ireland et al. (2016) that, for (C & 1, gravity reduces the relative impact of turbulence dispersion of inertial particles: gravitational drift causes them 35 Figure 3.2: PDF of the cluster lifetimes for the different cases. to cross fluid trajectories, leading to shorter particle-eddy interactions and a reduction of dispersion compared to tracers (Squires & Eaton, 1991, Wang & Stock, 1993, and more recently Parishani et al., 2015 and Mathai et al., 2018). Therefore, at least in the considered range of (C, gravity allows for larger and longer-lived clusters. This will be also corroborated by the relative dispersion results presented in §3.3.3. An exception is represented by the case (C = 4, ({ = 37, which shows relatively small and long-lived clusters. In this case, the massive gravitational drift does not allow sufficient time for the turbulence to engender large clusters; small groups of nearby particles are registered as clusters and travel together across the domain, with minor changes in their mutual position and therefore attaining long lifetimes. Considering the previous observation on the effect of inertia and gravity on the lifetime, the correlation between !� and C;8 5 4 is consistent with the notion that particles with higher St and Sv generally form larger clusters (Petersen et al., 2019). We remark 36 5 10 15 20 L C / 0 0.1 0.2 0.3 0.4 0.5 P D F lifetime = 1-2 lifetime = 3-4 lifetime = 5-6 lifetime = 7-8 (a) (b) Figure 3.3: (a) PDF of cluster size conditioned on the range of cluster lifetime for case St = 4, Sv =37. (b) Cluster size versus lifetime for the different cases. that !� was shown to be in approximately linear relation with the number of particles belonging to a cluster #%� (Baker et al., 2017), and therefore the relation between C;8 5 4 and !� is similar to that between C;8 5 4 and #%� . 3.3.2 Birth and death We now examine the type of events leading to the formation and destruction of clusters. A simple categorization is proposed, based on the particle origin and destination before and after death. In particular, births can originate from three types of events: coagulation, separation, and recombination. Coagulation refers to a cluster born from a majority of particles that do not belong to any clusters in the previous time step. For clusters which are not born by coagulation, we term separation the case in which a “single parent” cluster contributes more than half the particles of the newborn cluster. The remaining cases originate from the recombination of parts of multiple parent clusters. The causes of deaths are categorized analogously: disintegration, if the majority of the particles do not 37 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 disintegration absorption splitting (a) (b) (c) coagulation separation recombination Figure 3.4: (a) Relative probability of the combinations of different types of cluster birth and death for the case St=1, Sv =1.2. For each combination, the average cluster lifetime (b) and size (c) are normalized by their respective maximum values in the population. belong to any new cluster in the next time step; absorption, if one new cluster receives the majority of the particles; and splitting otherwise. Figure 3.4(a) graphically shows the probability of the different combinations of birth and death types, while figure 3.4(b) and 3.4(c) compare lifetime and size in those scenar- ios, respectively. The case (C = 1, ({ = 1 is displayed. The trends for the other St and Sv cases are qualitatively the same, being somewhat more accentuated with increasing inertia and gravity effect. The most common scenario is birth by coagulation and death by disintegration of small and short-lived clusters. Birth by recombination and death by splitting are relatively rare, but their occurrence is not inconsequential as they are typically associated to large and long-lived clusters. The close similarity between figure 3.4(b) and 3.4(c) corroborates the connection between cluster size and lifetime. 3.3.3 Relative dispersion of clustered particles Beyond the specific birth mechanisms, the formation of a cluster inherently implies that, in a mean sense, particles approach each other from some further distance. Likewise, 38 0 0.5 1 t/t life 1 2 3 4 5 6 7 (c)(a) (b) Figure 3.5: Evolution of the radius of gyration '6 of the particle clouds normalized by its average during C;8 5 4: (a) before birth; (b) during lifetime; and (c) after death, with lines indicating exponential best-fits. In panel (b), the temporal abscissa is normalized by the average cluster lifetime for each case. a cluster destruction implies that particles move away from one another. To investigate the rate at which these processes happen, we consider the set of particles belonging to each cluster and track them in time before birth and after death. Because we follow them beyond the time during which they belong to a cluster, we generally refer to this set of particles as a cloud. At each time instant, we quantify the characteristic size of the cloud using the radius of gyration '6, i.e. the mean particle distance from the cloud center of mass. In figure 3.5 we plot '6 as a function of t, during approximately 10g[ before birth (figure 3.5a) and after death (figure 3.5c). The results are ensemble-averaged over all clusters with C;8 5 4 ≥ 4g[ (a threshold which emphasizes but does not qualitatively alter the trends) and normalized by the average value during cluster’s lifetime. The latter, denoted as '6,0{6, is of the same order as the cluster sizes !� displayed in figure 3.3(b). These span a significant range of scales, and are generally outside of the dissipative range of the turbulent motions. At such scales, the flow is not differentiable and an approach based 39 on Lyapunov exponents in not applicable (Bec et al., 2007). Also, figure 3.5(b) shows the normalized '6 during the cluster lifetime, plotted against the normalized time C/C;8 5 4. Several points appear noteworthy. First, the clouds contract and expand exponentially in time before and after the cluster lifetime, respectively. (The exponential trend is apparent in semi-logarithmic versions of the plots, not shown.) Second, the contraction rate before birth is about twice as high as the expansion rate after death. This is consistent with the notion that relative particle dispersion is faster backward-in-time than forward-in-time (Sawford et al., 2005), and that such asymmetry is more pronounced for inertial particles (Bragg et al., 2016). Third, both contraction rates and expansion rates generally decrease with increasing St and Sv. The effect of inertia on particle relative dispersion has been explored only recently (Bec et al., 2010; Chang et al., 2015; Bragg et al., 2016): its impact on the particle dispersion rate compared to fluid tracers was found to be non-monotonic in time and strongly dependent on the initial separation. Because here we consider clusters with complex shapes and a broad range of sizes, a direct comparison with results obtained for particle pair dispersion is not straightforward. We notice, however, that the slower separation in presence of gravity is qualitatively consistent with the results of Chang et al. (2015). Fourth, and perhaps most striking, the radius of gyration shows only marginal variations during the cluster lifetime, indicating that the relative distance between the particles remain approximately unchanged during the lifetime. This is true even for the higher St and Sv cases, which have relatively long-lived clusters. 3.3.4 The role of the turbulent activity We finally turn our attention to the role played by turbulence in the clustering process. We use the local enstrophy l2 to characterize small-scale turbulence activity (Carter 40 0 0.5 1 t/t life 1 2 3 4 (a) (b) (c) Figure 3.6: Evolution of the enstrophy experienced by the particle clouds normalized by its average during C;8 5 4: (a) before birth; (b) during lifetime; and (c) after death. In panel (b), the temporal abscissa is normalized by the average cluster lifetime for each case. & Coletti, 2018). We interpolate l2 at the location of the particles belonging to the clouds considered in the previous subsection, and (similar to figure 3.5) we plot the ensemble-averaged values (normalized by their average during C;8 5 4) before the birth (figure 3.6a), during the lifetime (figure 3.6b), and after the death of the clusters (figure 3.6c). The fluid enstrophy sampled by the cloud decays approximately linearly during the time before birth, grows at a somewhat faster rate after death, and remains relatively constant during the lifetime. The rate of change of the sampled enstrophy before birth and after death is much steeper for particles with St and Sv of order unity. Consider the case (C = 1, ({ = 1.2. During lifetime, the clustered particles see an average enstrophy level around 24% of the unconditional average in the computational domain. About 10g[ before birth (after death), they sample fluid with enstrophy around 75% (98%) of the unconditional average. This provides a strong indication of the mechanistic influence of small-scale turbulence activity on the formation, survival, and destruction of inertial particle clusters: these are formed when a cloud of particles senses decreasing turbulence 41 activity, persist while their environment is relatively quiet, and are disrupted by the next rise of local turbulence fluctuations. This behavior appears to be most pronounced for the cases in which clustering is most intense, but is still present for particles with more inertia. We also note that using the local strain-rate instead of the enstrophy leads to quantitatively analogous results (not shown). Therefore, the trends are not the consequence of the oversampling of low-enstrophy regions, but reveal a more fundamental tendency of the clusters to survive in times when the small-scale turbulence activity is weak. 3.4 Conclusions We have carried out the first study of the temporal coherence of inertial particle clusters in turbulence. We have focused on a range of St (1 to 4) and Sv (0 to 37) leading to clusters of considerable size. The picture that emerges reveals several novel and unexpected aspects of particle-laden turbulence. The cluster lives have typical durations of a few Kolmogorov time-scales, consistent with the idea that clustering is primarily driven by small-scale turbulence. Cluster lifetimes are strongly related to their size, i.e. large clusters tend to be long-lived, and vice versa. Accordingly, clusters formed by particles with more inertia (which are on average larger) last longer in time. Gravitational settling also increases lifetime, because (for the present range of St) it reduces the influence of turbulent dispersion. Despite the strong link between size and lifetime, we remark that the former follows a power-law probability distribution (Monchaux et al., 2010; Monchaux & Dejoan, 2017; Baker et al., 2017; Petersen et al., 2019; Momenifar & Bragg, 2020, among others), strongly suggestive of a self-similar process; while the latter follows an approximately exponential distribution, indicating a dominant time scale. Therefore, the 42 relation between both quantities is non-trivial and its mathematical modeling requires further investigation. By tracking individual particles before and after the formation of a cluster, we gain insight into its birth, survival, and death. The most common scenario is the formation of small clusters from the coagulation of previously non-clustered particles, quickly followed by disintegration into parts mostly consisting of non-clustered particles. Large and long-lived clusters, on the other hand, are born most frequently by recombination of other similarly large clusters, and upon their death their components also recombine into new, large, and long-lived clusters. This implies that, albeit rarely, the same particles may remain in a “clustered state” for a time much longer than even the longest cluster lifetimes. This may have important consequences for the probability of particle-particle interaction. Particle clouds contract exponentially in time before giving birth to a cluster, and after its death they expand exponentially at a lower rate, consistently with the known asymmetry between forward-in-time and backward-in-time pair dispersion. This may also be related to the effective compressibility of the inertial particle field (Maxey, 1987). However, as recently pointed out by Oujia et al. (2020), the relation between the local concentration and the sign of the particle velocity divergence is not trivial. In future studies, this point may be investigated coupling our approach to the method of Oujia et al. (2020), who tracked the rate of change of the Voronoï cell volumes. Both contraction rates and expansion rates generally decrease with increasing St and Sv. For all considered cases, the cluster size is remarkably constant during its lifetime, and the particles experience anomalously low enstrophy and strain-rate of the fluid turbu- lence. Specifically, the cluster formation happens during a phase in which the small-scale turbulence activity decays, and its disruption coincides with the end of such quiet state. 43 Thus, in contrast with the common view that clustering is associated to intense turbulence, the clusters need a relatively quiescent environment in order to survive for an extended amount of time. As weak turbulence fluctuations are associated to reduced inter-scale energy transfers rates (Carter & Coletti, 2018), this points to a link between cluster for- mation and the turbulence cascade which deserves more attention. For similar reasons, further research is warranted to investigate higher Reynolds numbers for which intermit- tency is more intense, as this may have direct impact on the cluster lifetime. Moreover, the large-scale spatial organization of small-scale turbulence activity (also dependent on the Reynolds numbers) is likely to impact the life cycle of the clusters, that are shown to form and evolve over inertial–range scales. In the future, a more articulated view of the merging and splitting of clusters can be gained by the applications of graph theory, which has been successfully used to investigate the evolution of clusters of vortices in wall turbulence (Lozano-Durán & Jiménez, 2014). Also, light bubbles are known to produce even more concentrated and long-lived clusters (Tagawa et al., 2012; Mathai et al., 2020), and therefore the present methodology is expected to provide similar insight in bubble-laden turbulence. In general, the approach can be applied to any system of clustering elements that can be tracked in a Lagrangian framework. 44 Chapter 4 Interface resolved simulation of particles with different shapes in sediment transport 4.1 Introduction Sediment transport represents a complex interaction of hydrodynamic forces, gravity, and collision forces between particles and boundaries, particularly under strong shear flow. It is a commonly observed natural phenomenon extensively studied due to its ecological and geological significance (Geyer et al., 2001; Prosser et al., 2001). Over the years, significant efforts have been made to develop comprehensive models of sediment transport (Merritt et al., 2003; Durán et al., 2011; Kok et al., 2012; Valance et al., 2015; Rana et al., 2021), highlighting the ongoing research interest in this field. 45 Numerical simulations have become a pivotal tool in sediment transport studies, em- ploying various representations of particles, from Eulerian descriptions of sedimentation (Yang et al., 2017) to Lagrangian point models (Schmeeckle, 2014; Sun & Xiao, 2016; Zheng et al., 2021). These methods fail to capture direct interactions between particles and fluid, utilizing models where particle sizes are smaller than the fluid scales. Re- cently, particle-resolved direct numerical simulation (PR-DNS) has gained popularity for its detailed study of particle dynamics, utilizing the Immersed Boundary (IB) method. The immersed boundary method was initially developed by Peskin (1982) for flexible boundaries in flow fields and has been extensively adapted to accommodate interactions between flow and rigid particles (Mohd-Yusof, 1996; Uhlmann, 2005; Kim&Choi, 2006; Lucci et al., 2010; Breugem, 2012; Tschisgale et al., 2018), allowing for a high-resolution capture of the particle geometry and its interactions with the fluid. Collision models are also employed to accurately describe particle–particle and particle–wall interactions (Kempe & Fröhlich, 2012; Costa et al., 2015; Biegert et al., 2017; Jain et al., 2019). PR-DNS has facilitated in-depth studies on various aspects of sediment transport. For instance, Zhu et al. (2022) used PR-DNS to study saltation events and developed models describing key parameters of this process. Kidanemariam & Uhlmann (2014a, 2017) observed the formation of sediment patterns, noting that a computational domain length of at least 100 times the particle diameter is necessary for these patterns to emerge. Scherer et al. (2022) investigated the role of turbulent large-scale streaks in generating spherical sediment ridges, confirming that these ridges can arise from instabilities between the mobile sediment bed and turbulent flow. For non-spherical particles, studies have reported complex behaviors, such as non- uniform settling rates and orientation-dependent drag, which affect their trajectories in 46 unpredictable waysVoermans et al. (2017). Baker&Coletti (2022) found that in turbulent boundary layers, disks and fibers oversample high-speed zones near the wall and respond differently to near-wall turbulence compared to spheres, and the preferential orientation of the particles is also observed in the near wall region, aligning with the numerical observation (Shin & Koch, 2005; Mortensen et al., 2008; Marchioli et al., 2010, 2016; Challabotla et al., 2015a,b). In terms of particle-resolved simulations, Ardekani et al. (2017) performed simulations of oblate particles in turbulent channel flow and observed that near the wall, oblate particles experience significantly lower angular velocities and act as barriers damping turbulent fluctuations. Fornari et al. (2018) studied the settling of many oblate particles and noted increased settling speeds when oblate particles form clusters. For sediment transport, the particle shape impact on the fluid flow is a complex process. The wall shear stress is important for sediment transport, which is determined by the channel roughness (Blois et al., 2014; Patil & Fringer, 2022), and the channel roughness is influenced by the shape of the particle, since different particle shapes can result in different packing structures with varying permeability (Zhou et al., 2011). The drag and life coefficients are different for different particle shapes, even for the same shape, the drag and lift coefficients can be different for different orientations (Fröhlich et al., 2020), which can affect the trajectories of particles (Schmeeckle et al., 2001). Allen (2012) observed that spherical particles prefer rolling motion, while disc-shaped particles prefer sliding motion. As for sediment transport rate, Cassel et al. (2021) reported different bedload transport rates for particles with different shapes, again highlighting the importance of particle shape in sediment transport. 47 Given the complex interactions between non-spherical particles and fluid flow, PR- DNS is an ideal tool for studying the sediment transport of non-spherical particles. However, research on this topic is extremely limited, only Jain et al. (2021) performed a PR-DNS study of the sediment transport of four different particles, which shows that fluid flow, particle properties, and bedforms are highly influenced by the particle shapes. To contribute to the understanding of sediment transport of non-spherical particles, in this study, we perform PR-DNS of different shapes of particles, particularly spherical and oblate particles. The mixture of both particles is also considered, and the behaviors of the particles are compared. This study should provide another valuable dataset for the study of the sediment transport of non-spherical particles. 4.2 Simulation details 4.2.1 Numerical method In this study, particle-resolved direct numerical simulations are conducted using the immersed boundarymethod, as proposed byUhlmann (2005); Breugem (2012). The fluid flow is governed by the Navier-Stokes equations, which are solved using a second-order central difference scheme on a staggered Cartesian grid. Time integration is performed with the second-order Runge-Kutta (RK2) method. Each particle is individually resolved by Lagrangian marker points on its surface, and their dynamics are governed by equations that account for hydrodynamic forces, gravity, and collision forces. Particle–particle interactions are modeled using the soft sphere model proposed by Costa et al. (2015), while the approach of Ardekani et al. (2016) is applied to non-spherical particles. Further numerical details are available in chapter 2. 48 4.2.2 Simulation set-up In this study, the fluid field is modeled as open channel flow. The computational domain, defined by dimensions !G × !H × !I = (100 × 22.5 × 27)�eq, is consistent across all cases examined. Here, �eq is defined as the equivalent diameter, which is the diameter of a sphere that has the same volume as the irregularly shaped particle under study. The dimensions !G, !H, and !I correspond to the channel lengths in the streamwise (x), wall- normal (y), and spanwise (z) directions, respectively. For the current simulation, �eq is set to 163G, where 3G denotes the grid spacing. Consequently, the domain comprises 1600 × 360 × 432 uniform grid points. Boundary conditions are defined as follows: periodic in the G- and I-directions, free-slip on the top wall, and non-slip on the bottom wall. Similar to previous studies on particle-resolved sediment transport (Vowinckel et al., 2014; Kidanemariam & Uhlmann, 2014b, 2017; Jain et al., 2021), this study maintains a constant bulk Reynolds number '41 = *1� 5 /a 5 across all cases. Here,*1 represents the bulk velocity of the flow, and � 5 is the height of the fluid part, essentially maintaining a consistent flow rate & 5 since*1 = & 5 /� 5 . The sediment height is denoted as �B43 , and the total channel height !H is the sum of � 5 and �B43 . As the sediment bed evolves over time, both �B43 and consequently � 5 change, making � 5 a function of time. Therefore, these values are calculated a posteriori. Typically, channel flow is driven by a constant pressure gradient. However, to sustain the targeted bulk Reynolds number under dynamic conditions, a time-varying pressure gradient is applied. The methodology for calculating this pressure gradient will be detailed in § 4.2.3. In our simulation, we consider two types of particles: spherical and oblate particles with an aspect ratio of 0.5. Both particle types have the same equivalent diameter, �4@, 49 Case #( #$ '41 '4g �B43/� ΔG+ (ℎ &?/&?,A4 5 sphere 12000 0 3000 271.67 5.81 1.023 0.251 0.417 oblate 0 12000 3000 273.27 5.41 0.999 0.240 0.397 mix 6000 0 3000 271.57 5.45 0.995 0.238 0.444 Table 4.1: Simulation cases and parameters, #( is the number of moving spherical particles and #$ is the number of moving oblate particles. which is used throughout. To replicate the randomness of a natural sediment bed, a fixed layer of spherical particles with a diameter of �4@ is established at the channel’s bottom, as described in Jain et al. (2017). This fixed layer remains consistent across all simulation cases. Three different scenarios are examined, each involving 12,000 movable particles: Case sphere: all movable particles are spherical; Case oblate: all movable particles are oblate; Case mix: there is an equal number of spherical and oblate particles. Further details about these cases are provided in table 4.1. The initial arrangement of the particles in case mix is shown in figure 4.1. In our current study, all the resolved particles have the same relative density, d ′ = d?/d 5 − 1 = 1.5, which corresponds to the sand density in water. The galilei number �0 = (d ′�3 4@6)1/2/a 5 = 32.66 in all cases. The coefficient of restitution for collisions in the normal direction has a value of 4=,3 = 0.97 (Ge & Monroe, 2019), while that in the tangential direction is 0.3. The coefficient for friction `B = 0.2 is applied to all the cases. 4.2.3 Sediment height and pressure gradient To determine the fluid–sediment interface, we employ a solid phase indicator function q(x, C), as proposed byKidanemariam&Uhlmann (2014b). This function assigns a value 50 of unity to point G located within any particle, and zero otherwise. The instantaneous sediment height at every streamwise location, denoted as ℎB (G, C), corresponds to the height where the spanwise-averaged q(x, C) reaches a threshold of 0.1. The average height across the streamwise distance is then denoted as 〈ℎB〉G . The time-averaged sediment bed height is then calculated as: �B43 = 1 )0{4 ∫ C0+)0{4 C0 〈ℎB〉G (C)3C, (4.1) here, )0{4 signifies the duration over which the calculation is conducted. Fluid statistics are calculated by averaging the fluid part in the horizontal directions, and for a fluid quantity 5 , we have: 〈 5 〉(H) = 1 )0{4+ 5 ∫ C0+)0{4 C0 ∫ +0 5 (x, C) (1 − q(x, C))3+3C, (4.2) where+ 5 is the fluid part of the volume+0, andwe are sampling locations with q(x, C) = 0. + 5 can be calculated as: + 5 = 1 )0{4 ∫ C0+)0{4 C0 ∫ +0 (1 − q(x, C))3+3C. (4.3) In order to maintain a constant flow rate, the pressure gradient is adjusted every time step. At each time step, sediment height at every streamwise location (〈ℎB〉G) is calculated first, then the bulk velocity is defined as: *1 (C) = 1 + 5 ∫ !H 〈ℎB〉G ∫ +0 D(x, C) (1 − q(x, C))3+3H, (4.4) 51 where D(x, t) is the instantaneous streamwise velocity of the fluid field. The targeted bulk velocity, *target, is defined as *target = & 5 /(!H − 〈ℎB〉G). According to Vowinckel et al. (2014), to achieve the desired flow rate in the simulation, a PI-controller can be utilized to control the instantaneous pressure gradient, which can be calculated by: 3? 3G (C) = % (*target −*1 (C)) + � ∫ C C0 (*target −*1 (C))3C, (4.5) here % and � are coefficients for the proportional and integral terms respectively. In practice, we have % = 8 and � = 0.3, which can ensure the maximum deviation from the desired flow rate is less than 1%. The mean total shear stress is calculated as: gC>C = d 5 〈 3? 3G 〉 ∫ !H H (1 − 〈q〉G)3H, (4.6) where 〈 3? 3G 〉 is defined the time-averaged pressure gradient. The wall shear stress is then defined as g| = gC>C (H = �B43). Consequently, the friction velocity Dg can be calculated based on the value of g| in each case. Then the friction Reynolds number for each case can be calculated as '4g = Dg (!H − �B43)/a 5 . The shield number (ℎ = D2 g/(d ′ 6�4@) and grid resolution in wall unit ΔG+ = DgΔG/a 5 are also calculated for each case and reported in table 4.1. 4.2.4 Simulation initialization In this study, we applied the same initialization method as used by Vowinckel et al. (2014) and Jain et al. (2021), which has been shown to reduce the time required to reach 52 (a) (b) Bottom layer of fixed spherical particles Sphere Oblate Initial arrangement of mix case Figure 4.1: Visualization for (a) fixed layer of spherical particles and (b) initial particle arrangement for case mix. 53 a statistically steady state. This method begins with the simulation of open-channel flow over a fixed particle bed layer without any particle movement. Once turbulence is fully developed, a snapshot of this state serves as the initial flow field for subsequent simulations. Particles are then loosely packed into the simulation domain, ensuring no interactions between them; their initial placements remain consistent across all three cases studied. Figure 4.1(b) illustrates this setup. For case mix, a random number generator decides whether a location is filled with a sphere or an oblate particle. In case oblate and mix, oblate particles are randomly oriented, as shown in figure 4.1(b). Then the simulation starts, and the fluid flow is driven by the pressure gradient calculated as described in § 4.2.3, with particles settling and advancing simultaneously. The initial phase involves the consolidation of the loosely packed bed, which, as observed, takes approximately 400�4@/*1 time units for all cases. Contrary to expectations, the bed height slightly increases after this phase. This observation aligns with Jain et al. (2021), who attribute the phenomenon to the development of bedforms and the resultant sediment agitation. The simulation then runs for a total of )C>C0; = 8000�4@/*1 time units, and the statistics are calculated for the last )0{4 = 3500�4@/*1 time units. The development time for the sediment bed for each case is )34{ = 4500�4@/*1 time units, and as shown in Jain et al. (2021), this is enough time for the simulation to reach a statistically steady state. 54 Figure 4.2: Mean streamwise velocity profiles for three cases. The markers indicate the location of Sediment bed height �B43 . Figure 4.3: Mean solid fraction profiles for three cases. 55 4.3 Results 4.3.1 Fluid flow properties Our analysis begins by examining the mean streamwise velocity profiles for the three simulated cases, as shown in figure 4.2. This figure illustrates the comparison of mean velocities, with markers indicating the sediment height in each case. The velocity profiles for cases oblate and casemix develop at a lower region compared to case sphere. However, at higher regions, the mean velocity profiles for all three cases nearly converge. The primary differences in velocity profiles occur around the sediment bed interface, where the existence of particles significantly influences the mean velocity. Notably, case oblate and case mix exhibit similar trends in mean streamwise flow across the entire domain. Figure 4.3 shows the mean solid fraction profiles for the three cases. Within the sediment region, the solid fraction for spherical particles is approximately 0.6, aligning with experimental findings reported by Boyer et al. (2011),Aussillous et al. (2013), and Voermans et al. (2017). The mixed case exhibits a slightly higher solid fraction than the spherical case, while the oblate case has the highest solid fraction. This increase is attributed to the more compact packing of oblate particles due to their shape. The three mean solid fraction profiles are similar to each other near the sediment bed interface, suggesting that for the current simulation, the shape of the particles does not significantly affect the sediment structures near the sediment bed height �B43 . Although in the work of Jain et al. (2021), substantial differences in the solid fraction profiles between case oblate and others can be observed. This is because in their work, the oblate particles form high porosity structures at the bottom, causing a higher sediment bed compared with other cases, however, such high porosity packing is not observed in the current simulation. The 56 difference can be attributed to the difference in oblate particle aspect ratio, in their study, the aspect ratio is 0.29, while in the current study, the aspect ratio is 0.5. Different sediment heights, resulting from variations in solid fraction profiles, help explain the observed differences in mean streamwise velocity profiles. Both case oblate and case mix, which have lower sediment height (�B43) and consequently lower solid fractions near this region compared to case sphere, show faster developments in their streamwise velocity profiles. The similarity in the solid fraction profiles between case oblate and case mix suggests that a corresponding similarity in mean streamwise velocity profiles can be expected. Figure 4.4 presents the Reynolds normal stresses for the three cases under investi- gation. The peak of the streamwise stress components occurs just above the sediment height (�B43), while the peak values for the wall-normal and span-wise components are observed to be 2 − 3�4@ higher than �B43 . Across all cases, the maximum values in the streamwise direction are similar. However, the oblate case exhibits higher maximum values in both the wall-normal and spanwise directions compared to the other cases. The peak value location for the streamwise component follows the �B43 , where case oblate reaches the peak value first, followed by case mix, and finally case sphere. The same trend is observed for the wall-normal and spanwise components. The similarity of Reynolds normal stresses between case mix and case sphere can be partly explained by the differing solid fractions of the particles in case mix. Figure 4.5(a) illustrates the solid volume fractions for each type of particle in case mix. At the bottom, the fixed layer of spherical particles contributes solely to the solid fraction. Above this region, oblate particles add to the solid fraction as they settle on the fixed bed. Just above the fixed bed, oblate particles contribute more significantly to the solid fraction due to 57 (a) (b) (c) Figure 4.4: Reynolds normal stresses normalized with the bulk velocity: (a) streamwise components, (b) wall-normal components, and (c) spanwise components. 58 (b)(a) Figure 4.5: Solid volume fraction contribution from different particles in case mix: (a) spherical and oblate particle solid fraction, (b) ratio of spherical to oblate particle solid fraction. Red dashed line indicates the sediment bed height. their shapes, which make it easy to find locations between spherical particles. From this region upwards, the solid fraction of oblate particles begins to decrease rapidly, while that of spherical particles declines more gradually. Figure 4.5(b) displays the ratio of the solid fraction of spherical particles to that of oblate particles, highlighting that above the sediment bed, a greater presence of spherical particles results in the similar behaviors of Reynolds normal stresses observed in case mix and case sphere. 4.3.2 Particle properties Figure 4.6(a) compares the mean particle velocities in each case. At the same height, case mix exhibits the highest mean particle velocity, closely followed by case oblate, with case sphere displaying the lowest velocity. This trend correlates with the comparison of the mean streamwise velocity in figure 4.2, where case sphere demonstrates the lowest mean streamwise velocity near the sediment bed interface, aligning with the observed mean 59 (b) Figure 4.6: Mean particle velocities comparison: (a) comparison between simulated cases, (b) comparison within case mix. The figure only shows the region where the particles are present. (b)(a) Figure 4.7: Spanwise angular velocity comparison: (a) comparison between simulated cases, (b) comparison within case mix. The figure only shows the region where the particles are present. 60 Figure 4.8: Spanwise angular velocity comparison for spherical and oblate particles. (b)(a) symmetry axis Figure 4.9: Orientation of oblate particles: (a) angle between the major axis and wall normal direction, (b) orientation of the major axis. 61 particle velocities. Interestingly, the mean particle velocity for case mix is slightly higher than that for case oblate, despite case oblate showing a higher mean streamwise velocity. This counterintuitive result is explained by the shape of the particles, where spherical particles tend to follow the fluid flow better, which elevates their velocity compared to oblate particles. Given the similarities in mean streamwise velocities between case oblate and case mix as observed in figure 4.2, the slightly elevated mean particle velocity in case mix is to be expected. Further insights are provided in figure 4.6(b), which details the mean particle velocities within case mix itself. At the same height, oblate particles exhibit lower velocities than spherical particles, reinforcing the assertions made regarding the comparative analysis between case mix and case oblate. Figure 4.7(a) illustrates the mean spanwise angular velocities for each case, revealing that for both spherical and oblate particles, angular velocity increases with height starting from the region where particles begin to move, reaching a peak around H = �B43 + (1 ∼ 2)�4@. Beyond this height, the angular velocity begins to decline. For case oblate, the decay rate is much smaller compared with case sphere. Case sphere exhibits the highest spanwise angular velocity, followed by case mix and case oblate. This observation aligns with findings from Allen (2012) and Jain et al. (2021), which suggest that spherical particles tend to engage in rolling motions within the sediment, whereas oblate particles predominantly slide, exhibiting minimal rotational motion. Figure 4.7(b) further explores case mix, showing that spherical particles maintain higher angular velocities than oblate particles, consistent with the overall trend observed in figure 4.7(a). This consistency underscores the distinct movement behaviors between particle shapes. Continued investigation into angular velocities of identical particle types across dif- ferent scenarios is presented in figure 4.8. A comparative analysis of oblate particles 62 between case mix and case oblate shows similar behaviors, with angular velocities in case mix being marginally higher even beyond the sediment bed region. The trend of slowly decreasing angular velocities above H = �B43 + (1 ∼ 2)�4@ is consistent with the behaviors observed in case oblate. For spherical particles, case mix follows the trend of case sphere up to its sediment bed height. As case sphere exhibits a higher �B43 , rolling motions extend into higher regions. Conversely, at the same height in case mix, the angular velocity of spherical particles begins to decrease. The comparison indicates that the intrinsic behaviors of the particles remain consistent even in complex scenarios. Non-spherical particles tend to exhibit preferential orientation in the flow. In DNS simulations of turbulent channel flow incorporating inertial point particles, preferential alignment near the wall has been observed (Marchioli et al., 2010; Challabotla et al., 2015a; Zhao et al., 2015). This phenomenon has also been documented in experimental studies(Baker & Coletti, 2022). Additionally, (DiBenedetto et al., 2018) observed pref- erential orientations for anisotropic particles under wave conditions. In nature, clasts of preferred orientation are often observed (Fairbridge & Bourgeois, 2006). In the current study, to investigate the preferential orientation of the particles, we calculate the mean angle between the major axis of the oblate particles and wall-normal direction, as defined in figure 4.9(a). The result of the mean angle value in figure 4.9(b) shows that for oblate particles inside the sediment bed, the mean orientation of the particles is around 125◦, for both oblate particles in case oblate and case mix. Above the sediment bed, the preferred orientation is less pronounced and becomes increasingly random, and the trend for the angle change is the same for oblate particles in both cases. This orientation of the oblate particles resembles the orientation of disk-like stones near the river bed regions. According to Kidanemariam & Uhlmann (2017), the instantaneous volumetric flow 63 rate of the particles, &?, can be calculated as: &? (C) = c�3 6!G!I #?∑ ;=1 D;? (C), (4.7) this formula sums all streamwise particle velocities multiplied by the particle volume and is divided by the product of domain length and width. To interpret its value,&? is usually normalized by a reference value, &?,A4 5 = (d ′ 6�3)1/2, where d ′ = d?/d 5 − 1 represents the relative density. The normalized mean volumetric flow rate are reported in table 4.1. Comparing to the empirical formula proposed by Wong & Parker (2006), the normalized mean volumetric flow rate shows reasonably good agreement with the empirical formula, which is around 0.4. The agreement with of the sediment transport with empirical laws is also observed in Jain et al. (2021), except the oblate case, where higher sediment transport rate is observed, which is caused by the different packing patterns of the oblate particles in that case. Our computational domain extends 100�4@, which, according to Kidanemariam & Uhlmann (2017), is the minimal length required for dune-like structures to emerge. In our simulations, wave-like bedforms were generated in all cases, as depicted in the top view of the sediment bed shown in figure 4.10. Each case features a large streamwise bedform structure. Comparing case oblate (figure 4.10(b)) with case sphere (figure 4.10(c)), the spanwise distribution of oblate particles is less uniform. The bedforms exhibit temporal movement. To illustrate the dynamic evolution of the bedforms, we plotted the spatio- temporal evolution of the sediment bed interface, ℎB (G, C), in figure 4.11. Heights are normalized by equivalent particle diameter (�4@), and each figure captures the evolution over a 4000�4@/*1 period. The movement of the bedforms is clearly observable in 64 (a) (b) (c) Figure 4.10: Snapshots of the sediment bed for case (a) mix, (b) oblate, and (c) sphere. (b)(a) (c) Figure 4.11: Spatial-temporal evolution of sediment bed interface in case (a) mix, (b) oblate, and (c) sphere. 65 (a) Figure 4.12: Shear stress and its different components: (a) total shear stress, (b) viscous shear stress, (c) Reynolds shear stress, and (d) particle-induced shear stress. each case. Given the limited length of our computational domain, caution is advised when interpreting the wavelength and speed of these bedforms. However, compared to the spherical and mixed cases, the structure of the bedforms in the oblate case shows more pronounced variations over time. Additionally, the differences in bedform structure highlighted in figure 4.10 underscore the significant impact of particle shape on bedform dynamics. 4.3.3 Shear stress decomposition To understand how the shear stress varies in different simulation cases, shear stress decomposition can be performed. The total shear stress in balance by the sum of viscous 66 shear stress, based on the formulation proposed in Nikora et al. (2007), we have: gC>C = d 5 〈 m? mG 〉 ∫ !H H 〈q(H)〉3H = d 5 a 5 3〈D〉 3H − d 5 〈D ′ { ′〉 + g?, (4.8) where g? denotes the particle-induced stress. The contribution of each component are illustrated in figure 4.12. The viscous shear stress (figure 4.12(b)) displays similar values across all cases away from the sediment interface. However, notable differences are observed below this interface, highlighting the significant impact of particle shapes on the flow dynamics. The Reynolds shear stress, detailed in figure 4.12(c), also follows similar trends across all scenarios, with a slightly higher peak value observed for case oblate. Figure 4.12(d) presents the particle-induced shear stress, which strongly correlates with the solid fraction in each case. Typically, particle-induced shear stress is modelled solely as a function of the solid fraction. However, in our simulations at the sediment bed interface (�B43), where a solid fraction of 0.1 is achieved across all cases, we observed differences in the particle-induced shear stress, indicating that particle shape influences the particle-induced shear stress. These findings suggest that particle shape should be considered in future modelling efforts to accurately predict shear stress distributions in sediment transport. 4.4 Conclusions Particle-resolved simulations for sediment transport involving spherical and oblate parti- cles were conducted across three distinct cases – case sphere, case oblate and case mix – each featuring 12,000 mobile particles. These simulations provided extensive statis- tical data, revealing both similarities and differences in particle and fluid dynamics as 67 influenced by particle shape. Our simulations demonstrated similar mean streamwise velocity profiles across all cases, with notable differences occurring near the sediment bed interface. Particle shape influenced sediment packing, altering the solid volume fractions and consequently af- fecting the mean velocity near the sediment bed interface. This variance in velocity was primarily observed in this region, whereas further from the wall, the impact of particles diminished, resulting in consistent velocity trends in the upper flow regions due to the implementation of a constant flow rate. For Reynolds normal stress, significant increases were observed in the wall-normal and spanwise directions when comparing case oblate with case sphere. Case mix exhibited behavior akin to case sphere, likely due to a higher concentration of spherical particles just above the sediment bed height. At the same height, case sphere displayed the lowest velocities due to higher solid volume fractions. Interestingly, in case mix, spherical particles consistently demon- strated higher velocities than oblate particles at comparable heights. Angular velocity profiles were similar across all cases for oblate particles throughout the domain, whereas for spherical particles, profiles converged only within the sediment bed. Additionally, preferential orientation was observed for oblate particles. Wave-like bedforms were generated in all simulated cases. By examining the spatio- temporal evolution of the sediment bed interface, distinct differences in transport pro- cesses were identified. Shear stress decomposition revealed that particle-induced shear stress is strongly correlated with the solid fraction of the particles. Further investigations are warranted to elucidate the effects of particle shape on sediment transport more comprehensively. In the work of Jain et al. (2021), large porosity is observed in oblate particle simulations with an aspect ratio of 0.29, which led 68 to significant differences in the mean velocity profile. While our simulations reported similar porosity and mean velocity profiles across all cases, this suggests a potential critical aspect ratio that could dictate the behavior of oblate particles. Moreover, the much-increased sediment transport rate for this case is reported in Jain et al. (2021), while in our simulation, the sediment transport rate shows good agreement with the empirical formulation in Wong & Parker (2006). This encourages a deeper exploration of the impact of aspect ratio on sediment transport, potentially contributing to empirical formulations for estimating sediment transport rates. Overall, our simulations provide a valuable dataset for understanding sediment trans- port dynamics involving different particle shapes, serving as a crucial reference for future numerical or experimental studies on this topic. 69 Chapter 5 Particle-resolved simulation of sediment transport over vibrating beds 5.1 Introduction Earthquakes are a major mechanism for the initiation of deep-sea turbidity currents. Un- derwater avalanches of sediment are believed to be relevant to the records of seismic history, offering glimpses of past earthquakes that have shaped our world (Adams, 1990; Goldfinger et al., 2007, 2012; Ashi et al., 2014). Meanwhile, other mechanisms are also responsible for initiating deep-sea turbidity currents, such as storm surges, floods and volcanic eruptions (Goldfinger et al., 2012; Ikehara et al., 2020). Among these trig- gers, surface-sediment remobilization emerges as a pivotal player in earthquake events (McHugh et al., 2016; Moernaut et al., 2017). Our understanding of surface-sediment remobilization remains tantalisingly incomplete. As Atwater et al. (2014) pointed out, even earthquakes of identical magnitude in the same locale can yield distinct turbidity 70 Figure 5.1: Sketch of laboratory experiment setup. signatures, underscoring the importance of more studies on surface-sediment remobiliza- tions. To understand the mechanism of surface-particle remobilization in seismic events, this study performs numerical simulations of sediment in turbulent channel flows with different oscillation conditions. Meanwhile, a series of laboratory experiments are per- formed by our collaborators. The experiments are performed in a channel, as shown in figure 5.1, affixed to an oscillator capable of imparting vertical oscillations at specified frequencies and amplitudes. Precise control over the flow rate within the channel is achieved, allowing for the manipulation of hydraulic conditions. The channel is outfitted with a variety of sediment types commonly found on the seabed, including coarse sand, fine sand and mud. This setup enables the emulation of vertical oscillations akin to those observed during seismic events. The effects of different factors on sediment transport 71 dynamics can be studied by varying the flow rate and the parameters of vertical oscillation. In our simulations, we investigate both vertical and horizontal oscillations, as well as their combinations. Specifically, the horizontal oscillation scenario involves applying an oscillating acceleration, which effectively renders the setup analogous to pulsating turbulent channel flows. Pulsating turbulent flows over smooth walls have been studied extensively both experimentally (e.g. Tu & Ramaprian, 1983; Mao & Hanratty, 1986; Tardu & Binder, 1993; Tardu et al., 1994; Lodahl et al., 1998) and numerically (e.g. Scotti & Piomelli, 2001; Manna et al., 2012, 2015; Papadopoulos & Vouros, 2016; Weng et al., 2016; Cheng et al., 2020), for channel and pipe flows, in which various combinations of governing parameters have been studied. Recently, pulsating flows over rough and bumpy walls have also been studied. Jelly et al. (2020) compared pulsatile and non-pulsatile rough-wall turbulent pipe flows and found that the major difference between these two cases in the current-dominated very-high-frequency regime is confined to the near-wall region. Ghodke & Apte (2016) and Patil & Fringer (2022) studied pulsating flows in turbulent channels with rough walls. They found enhanced drag for rough walls. They also examined the phase-dependent turbulent kinetic energy (TKE) budget, finding that there is a strong attenuation in the net production for the pulsating rough wall cases. The present study employs particle-resolved direct numerical simulation (PR-DNS) to investigate sediment transport behaviours in oscillatory channels. PR-DNS, rooted in the immersed boundary (IB) method initially proposed by Peskin (1982), for flexible boundary simulation in a flow field, has been further developed by many researchers (e.g. Mohd-Yusof, 1996; Uhlmann, 2005; Kim & Choi, 2006; Lucci et al., 2010; Breugem, 2012; Tschisgale et al., 2018) to study flow-particle interactions. The interactions among particles are described by collision models, which typically fall into three categories: 72 repulsive potential models (Glowinski et al., 2001; Uhlmann, 2008), hard sphere models (Veeramani et al., 2009; Derksen, 2011; Jain et al., 2019) and soft sphere models (Kempe & Fröhlich, 2012; Costa et al., 2015; Biegert et al., 2017). PR-DNS has been increasingly utilised in sediment transport research in recent years. Kidanemariam & Uhlmann (2014b) examined the pattern formation in flows over erodi- ble beds, observing distinct dune formation under laminar and turbulent flow conditions. Subsequently, Kidanemariam & Uhlmann (2017) extended their study on pattern forma- tion and discovered that the computation domain length needs to exceed 100 times the particle diameter to observe at least one unstable mode. Recent studies by Mazzuoli et al. (2019) and Vittori et al. (2020) delved into sediment transport in oscillatory flows. They noted ripple formation under such conditions. They also proposed refined empirical formulas for a more accurate estimation of the sediment transport rate. In the present work, we use PR-DNS to investigate sediment transport behaviours in channel flows under various oscillation conditions. We apply periodic oscillations in vertical and horizontal directions to simulate different seismic conditions. Our objective is to examine flow patterns under different oscillation conditions. Particle behaviours in different oscillatory states are also reported. Our goal is that, starting from the canonical oscillation conditions, our results can elucidate the intricate interplay between fluid dynamics and sediment transport, paving the way for understanding sediment transport dynamics under natural seismic conditions. 73 5.2 Simulation details 5.2.1 Reference frame Our numerical simulations aim to study the sediment transport behaviour under the various conditions of the oscillating bed. As shown in figure 5.1, we fix the reference frame, denoted by (, with the oscillating system. Without loss of generality, for a typical vertical oscillation, the displacement H in the Earth frame (′ is described as H(C) = � sin(lC), (5.1) and the acceleration is given by 0(C) = −�l2 sin(lC) = 0H, (5.2) The frame ( is a non-inertial frame. For governing equations of velocities observed in this non-inertial frame, an inertial force C = d 5 m? mG ∫ � H (1 − q̄(H))3H, (5.15) The wall shear stress is then defined as g| = gC>C (H = �B43). Consequently, the frictional velocity Dg can be calculated based on the value of g| in each case. Then the friction Reynolds number for each case can be calculated as '4g = Dg (� − �B43)/a 5 . The sediment height and friction Reynolds number for each case are listed in table 5.1. 80 Figure 5.2: Profiles of (a) solid fraction, (b) mean flow velocity obtained for each case. 5.3 Results 5.3.1 Basic statistics of the flow and particles We begin our analysis of the simulation results by examining the basic fluid and particle statistics for the different cases. The results are presented in two forms, time-averaged and phase-averaged. Time-averaged results We first compare the time-averaged results among the five cases. Figure 5.2 shows the profiles of time-averaged solid fraction, q̄(H), and mean flow velocities, D̄(H). As illustrated in figure 5.2(a), the first four cases exhibit similar solid fractions within the sediment bed, approximately 0.6, consistent with previous experimental results (Boyer et al., 2011; Aussillous et al., 2013; Voermans et al., 2017). At the top of the sediment bed, case V has slightly higher solid fractions than the other three cases. Conversely, case H+V demonstrates a markedly different solid fraction distribution. Further analysis 81 x z y Figure 5.3: Snapshot of the waveform of sediment bed obtained in case H+V. Figure 5.4: Profiles of (a) time-averaged Reynolds normal stress D′D′ (—), {′{′ (- -),|′|′ (– · –), and (b) time-averaged Reynolds shear stress −D′{′. 82 Figure 5.5: Profiles of particle-induced stress g? for the first four cases. in §5.3.3 reveals that this unique distribution in case H+V is due to the generation of a large amplitude waveform. Figure 5.2(b) shows the mean velocity profile for each case. Case N exhibits the highest mean velocity. Despite the presence of vertical oscillations in case H-V and their absence in case H, these two cases have similar mean velocities, which are slightly lower than those in case N. This finding aligns with Patil & Fringer (2022), who noted that a pulsating flow over a rough wall results in reduced mean velocity compared to non-oscillating flow. Interestingly, in case V, the vertical oscillations reduce the mean flow to a level lower than that in case H and case H-V. This reduction is partially due to the slightly higher solid fraction near the sediment surface and the increased sediment transport rate (shown in table 5.1 and analysed in §5.3.3) in case V, which slows down the mean flow. Similar to the solid fraction distribution, the mean velocity in case H+V is distinctly different from the other cases, displaying an almost linear profile within the sediment region. The effect of the different solid fractions goes beyond the sediment region, we can observe a significantly slower growth of the mean velocity profile in case 83 H+V. Even in the outer region where no particles presented, it shows a slower growth rate compared to other cases. The above time-averaged results indicate that caseH+V is distinct from the other cases. In our simulation, we observed that the waveform (shown in figure 5.3) reaches the fixed bottom layer of random sediment particles, and the height of the sediment bed can reach above a height of H = 0.35. As shown in Ghodke &Apte (2016), for oscillation flows over rough walls, the effects of the rough wall can extend to more than twice the highest point of the rough wall. If the same criterion is applied here, the effects of the sediment bed penetrationwould be higher than our simulation domain, so the fluid structures induced by this large sediment structures necessitate amuch larger computational domain. Therefore, the fluid flow results for case H+V should be interpreted with caution. Additionally, figure 5.4 shows the Reynolds stress components for the first four cases. In each case, the maximum value of streamwise velocity fluctuations occurs slightly above the sediment bed height, whereas the maximum fluctuations in the other two velocity components occur at higher locations, consistent with the findings of Jain et al. (2021). Compared to case N, the streamwise fluctuations are higher within the sediment bed for the cases with oscillations and lower away from the sediment bed. For fluctuations in the other two velocity components, the values in case N are smaller than those in the other cases. Figure 5.4(b) compares the Reynolds shear stress of the first four cases. Case V exhibits the highest values within the sediment bed. At the location slightly higher than the sediment bed, case N and case H have the largest values. Following the examination of the profiles of the fluid flow’s mean velocity and Reynolds shear stress, it is logical to assess the particle contribution to the mean shear stress next. The driving force of the flow should be balanced by the total shear stress 84 in the system, which includes viscous shear stress, Reynolds shear stress and particle- induced stress (Kidanemariam & Uhlmann, 2017). Although the horizontal acceleration of the non-inertial frame alters the flow momentum, its time-averaged contribution is zero. Thus, the total shear stress can still be calculated using the following equation proposed by Nikora et al. (2007): d 5 m? mG ∫ !H H q̄(H)3H = d 5 a 5 3D̄ 3H − d 5 D′{′ + g?, (5.16) where g? denotes the particle-induced stress. Figure 5.5 shows g? for the first four cases. Form q̄ shown in figure 5.2(a), we can see that just below the sediment bed height, case V has a smaller solid volume fraction, followed by case H. This trend is reflected in the particle-induced stress shown in figure 5.5. However, despite that case N and case H-V have similar q̄ profiles, the particle-induced stress for case H-V is larger than that for case N, indicating increased stress due to oscillations. Away from the sediment bed height, the particle-induced stress is higher for the cases with higher solid fraction. Additionally, figure 5.6 shows basic sediment properties for case N. In this case, only a few particles are present in the top layer, resulting in a non-smooth profile above the height H/� = 5. However, despite the limited number of moving particles, figure 5.6(a) shows that below the sediment bed height �B43 , even though the fluid velocity is not zero, the particle velocity is zero due to the packing of the sediment bed. Above �B43 , the mean sediment particle velocity increases rapidly with fluid velocities. Figure 5.6(b) shows the mean angular velocity of the particles. Around the height H/� = 4 − 5, the increase in the angular velocities indicates the rolling motion of the particles in this region. 85 (a) (b) Figure 5.6: Profiles of (a) mean particle velocity and (b) mean particle angular velocity in the spanwise direction for case N (The red dashed line indicates H = �B43), also here we only show the regions where the particle statistics are available. Phase-averaged results Figure 5.7(a) shows the phase-averaged fluid velocity for case V. We observe that the fluid velocity does not change significantly across different phases. Figure 5.7(b) presents the phase-averaged solid volume fraction, indicating a clear elevation at the top layer; however, the solid volume fraction value there is small and does not significantly affect the fluid motions at higher locations. Figure 5.8 displays the phase-averaged profiles of root-mean-square(r.m.s.) of velocity fluctuations for case V. There is a minimal difference among phases, especially away from the sediment bed, further confirming that the fluidmotion at higher locations is not affected by the slight increase in solid volume fraction. For the streamwise velocity fluctuations, although small changes in the peak value are observed, the location of the maximum value remains fixed. Near the sediment bed, phase \ = 0 exhibits the largest streamwise velocity fluctuations and the smallest fluctuations in the other two velocity components, while phase \ = c/2 shows the opposite pattern. 86 Figure 5.7: Profiles of the phase-averaged (a) fluid velocity and (b) solid volume fraction for case V. Figure 5.8: Profiles of phase-averaged r.m.s. values of velocity fluctuations for case V D ′ D ′ (—), {′{′ (- -),|′|′ (– · –). 87 Figure 5.9: Profiles of phase-averaged particle velocity for case V. The dashed lines indicate the time-averaged particle velocity. Although the fluid field does not change much across different phases, sediment movements strongly correlate with the phase. Figure 5.9 compares the phase-averaged particle velocities with the time-averaged particle velocity. From phase \ = 0 to \ = 3c/4, the mean velocity of the sediments increases at all height, while from phase \ = c to \ = 7c/4, the mean velocity decreases. Additionally, the velocity profiles indicate that particles can reach different heights at different phases. At \ = c, particles reach the highest location, causing the elevated solid fraction observed in figure 5.7(b). Since the mean fluid velocity does not vary significantly across the different phases, the relative velocity between the fluid and particles, 〈DA〉 = 〈D〉 − 〈D?〉, depends solely on the particle velocity. Thus, from phase \ = 0 to \ = 3c/4, 〈DA〉 increases and from phase \ = c to \ = 7c/4, the relative velocity decreases. Figure 5.10(a) shows the phase-averaged fluid velocity for case H. Due to the large oscillation of the horizontal acceleration, the phase-averaged fluid velocity oscillates 88 (a) (b) (c) Figure 5.10: Profiles of phase-averaged (a) fluid velocity, (b) particle velocity, (c) relative velocity of the fluid with respect to the particles for case H. The dashed lines indicate time-averaged values. 89 around its time-averaged profile starting from the sediment bed height. From phase \ = 0 to c, the mean velocity accelerates, while from \ = c to 2c, it decelerates. The solid volume fraction for case H also does not significantly across different phases (results not plotted here). Figure 5.10(b) displays the phase-averaged particle velocity at different phases. Sedi- ment particle velocity increases during the fluid accelerating phases and decreases during the fluid decreasing phases. Figure 5.10(c) illustrates the relative velocity between the particles and fluid, with the dashed line representing the time-averaged relative velocity. It should be noted that D̄A ≠ 〈DA〉, because the height of the particles can change at different phases, resulting in sampling only at high fluid velocity phases. From phase \ = 0 to \ = c/2, despite the increasing particle velocity, the relative velocity 〈DA〉 still rises (figure 5.10(c)) as the fluid velocity increases more than the particle velocity increases. The relative velocity decreases during the other phases. This behaviour contrasts with case V, where no horizontal acceleration present and particles can only gain momentum from the fluid flow. In case H, horizontal acceleration acts as a driving force for both the fluid flow and sediment particles. The phase-average values for velocity fluctuations of case H are shown in figure 5.11. Unlike the results of case V shown in figure 5.8, substantial changes in the velocity fluctuations are observed, especially for the horizontal fluctuations. Both the maximum values and their locations vary across different phases. A detailed analysis of how the velocity fluctuations change over different phases is provided in § 5.3.2. For case H-V, figure 5.12(a) demonstrates the phase-averaged fluid velocity, which is very similar to the fluid profiles for case H shown in figure 5.11(a). This similarity is expected since the horizontal accelerations are the same, and there is little change 90 Figure 5.11: Profiles of phase-averaged r.m.s. values of velocity fluctuations for case H (a) (b) (c) Figure 5.12: Profiles of phase-averaged (a) fluid velocity, (b) particle velocity, (c) relative velocity of the fluid with respect to the particles for case H-V. The dashed lines indicate time-averaged values. 91 (a) (b) (c) Figure 5.13: Profiles of phase-averaged (a) fluid velocity, (b) particle velocity, (c) relative velocity of the fluid with respect to the particles for case H+V. The dashed lines indicate time-averaged values. Figure 5.14: Profiles of phase-averaged solid volume fraction for case H+V. 92 observed in the sediment bed during the simulation. Figure 5.12(b) presents the phase- averaged particle velocity. In the region near the sediment bed, there is a good agreement with the time-averaged velocity for most cases. In contrast, in cases V and H only a few phases have good agreement (figures 5.9(b) and 5.10(b)). This difference indicates that the particle movement behaviour for case H-V is distinct from the previous cases V and H. Additionally, although it is not obvious in figure 5.12(b), slightly negative mean particle velocities are seen for some phases. Figure 5.12(c) shows that the fluid relative to particle velocity increases during the fluid accelerating phases and decreases during the fluid decelerating phases, similar to case H. The velocity fluctuations in case H-V also align with case H and are analysed in § 5.3.2. Figure 5.13 shows the results for case H+V. With horizontal accelerations, velocity oscillations away from the sediment bed appear as expected. Additionally, the combined horizontal and vertical accelerations contribute to sediment motions at the same phases, resulting in much larger sediment transport for this case. Figure 5.13(a) presents the mean velocity profiles. It shows that the height of the linear velocity zone varies at different phases, corresponding to the solid volume fraction changes over different phases, as illustrated in figure 5.14. Figure 5.14 shows a much larger elevation of the solid volume fraction in case H+V than in case V (figure 5.7(b)). We also observe a clear reverse flow close to the bottom of the domain. Figure 5.13(b) demonstrates the phase-averaged sediment particle velocities. Unlike other cases with noticeable particle velocities at all phases, particle velocities are close to zero at phase \ = 0 and \ = 7c/4, showing a much-confined particle movement at those phases. Additionally, at phase \ = 3c/2, negative particle velocities at lower heights are 93 observed. Figure 5.13 (c) shows the fluid relative to particle velocities. For the time- averaged velocities, the time-averaged particle velocities are larger than the time-averaged fluid velocity at the bottom and at higher locations. Near the bottom, during the reverse flow phases, the particles have zero velocity, resulting in a positive relative velocity at phase \ = 0 and \ = 7c/4. At higher locations, close to H/� = 7, particles can only reach this height during high fluid velocity phases. Therefore, sampling the particle velocities in that region yields higher values than the time-averaged fluid velocities (figure 5.13 (c)). However, the phase-averaged relative velocity shows that particle velocities in these regions close to the sediment bed height are always smaller than the fluid velocities (figure 5.13(c)). 5.3.2 Comparison with pulsating turbulent flows In the present study, cases involving horizontal accelerations exhibit similarities to pul- sating turbulent flows. In this section, we compare our results with the results of pulsating turbulent flows reported in the literature. Amplitude and phase shifts of the velocity oscillations Pulsating pipe and channel flows can be categorized as either current-dominated or wave- dominated, depending on the ratio between the oscillation amplitude and the bulk velocity. In our simulations, this ratio is less than unity, indicating current-dominated flows. Scotti & Piomelli (2001) found that current-dominated flows are primarily influenced by the frequency rather than the amplitude of the oscillation. Current-dominated flows have been extensively studied and are typically categorized into different regimes based on the frequency parameter. In our study, we use the 94 inner scale frequency l+ = la 5 /D2 g and characteristic length scale ;+B = ;BDg/a 5 (where ;B = √ 2a 5 /l denotes the Stokes layer thickness) as proposed by Tardu et al. (1994). Based on the value of l+, pulsating flows are generally classified into five regimes (Tu & Ramaprian, 1983; Manna et al., 2012; Papadopoulos & Vouros, 2016). 1. Quasi-steady regime (l+ → 0): The flow operates like a steady flow, with no velocity overshoot or phase lag in the ensemble-averaged quantities compared to steady values with a constant pressure gradient. 2. Low-frequency regime (l+ < 0.005): While the time-averaged flow is not significantly different from the quasi-steadyflow, the oscillation affects the phase-averaged velocity across the entire shear layer. 3. Intermediate-frequency regime (0.005 < l+ < 0.02): There exists interactions between the imposed oscillation and turbulent structures, and the interactions are more significant with increasing frequency. 4. High-frequency regime (0.02 < l+ < 0.04): The effects of oscillations are strong but confined to the viscous sublayer. Inside this region, the amplitude of the oscillations and the phase lag compared with centerline velocity follow the Stokes laminar values. Outside this region, the bulk flow behaves like a solid body movement. 5. Very-high-frequency regime (l+ > 0.04): The interaction between the oscillation and turbulence is very strong, though confined in a small region. This regime is also characterised by the fact that the oscillation frequency is similar to the bursting frequency of the steady turbulent flow. In our simulations, all cases with horizontal acceleration have l+ > 0.04, falling in the very-high-frequency category. For this regime, past experiments (Brereton & Hwang, 1994) and simulations (Scotti & Piomelli, 2001; Manna et al., 2012, 2015; Papadopoulos 95 & Vouros, 2016) have demonstrated that the oscillatory flow near the wall adheres to the Stokes solution. Jelly et al. (2020) examined the oscillatory flow in rough pipes. They found that roughness height can significantly affect the comparison between the oscillatory flow and Stokes solution in the near-wall region, whereas the comparison is not affected away from thewall. It would be intriguing to ascertain if the same comparison can be observed in the present problem. The Stokes solution can be written as: DB (A) = D<'4 [ 8 ( �0(83/2 √ 2A/;B) �0(83/2 √ 2�/;B) − 1 ) 48lC ] , (5.17) where �0 is the Bessel function of the first kind of order zero, 8 is the imaginary unit and '4[] is the real part of the argument. The phase-averaged velocity can be considered as the combination of a mean part and an oscillating part. Tardu et al. (1994) expanded the phase velocity 〈D〉 using Fourier series: 〈D〉(H) = D̄(H) + �D̃ (H) cos(lC + \D̃ (H)) + =∑ 8=2 �D̃8 (H) cos(8lC + \D̃8 (H)), (5.18) where �D̃ and \D̃ are the amplitude and phase of the fundamental mode, which describe the oscillating components of the quantities that follow the driving frequency. After the decomposition using 5.18, the fundamental mode at each wall normal distance can be obtained, and then �D̃ (H) and \D̃ (H) are compared with the centerline velocity component. In our results, �D̃ (�) and \D̃ (H) are obtained at the free-slip boundary. Comparison of our simulations with the Stokes solution begins at the sediment height location, and the comparison is performed using the height coordinates denoted as (H − �B43)/;+B . Figure 5.15(a) and (c) depict the amplitude ratio comparison for cases H and 96 (b) Figure 5.15: Comparison of velocity amplitude and phase lag of the oscillating velocity with the Stokes solution for case H (a, b) and case H-V (c, d). 97 Figure 5.16: Comparison of velocity amplitude and phase lag of the oscillating velocity with the Stokes solution for case H+V (a, b) H-V, revealing a reasonably good agreement with the Stokes solution. At a distance far from the sediment bed, the amplitude ratio equals 1, consistent with the Stokes solution. In our simulations, sediment particles ascend in each period, impacting the fluid velocities periodically. Figure 5.15(a) and (c) demonstrate that in the region approximately (H−�B43)/;+B = 2−6, the velocity amplitude ratio exceeds unity, indicating a longer region compared to the Stokes solution. Furthermore, unlike the rapid decline of the amplitude for the oscillating velocity component near the wall as described by the Stokes solution, our simulation exhibits a much slower decrease into the sediment bed. Figures 5.15(b) and (d) illustrate the phase lag relative to the free-slip boundary, revealing only minor differences even down to the sediment bed height location. Notably, the comparison for phase lag extends into the sediment bed, where it is observed that the phase lag continues to increase, until the region where the amplitude ratio vanishes. Figure 5.16 presents the comparison with the Stokes solution for case H+V, which stands out as the only case demonstrating large waveform development (figure 5.3) under 98 horizontal accelerations. This unique characteristic underscores its significance. While the amplitude comparison in case H+V resembles that in cases H and H-V above the sediment bed height, case H+V displays an even slower decrease within the sediment bed. Additionally, compared to cases H and H-V, the phase lag in case H+V notably decelerates around and below the sediment height. Similar trends were observed by Jelly et al. (2020) in rough-wall pulsatile turbulent pipe flows. The above results suggest that the formation of bed forms profoundly influences the flow near the surface sediment, which in turn affects the evolution of particle patterns in this region. Turbulence intensity modulation In this section, we examine the modulation of the turbulent velocity fluctuations by the oscillatory bed. The oscillatory modulation is the phase-dependent component, which can be calculated using formula 5.10. Figures 5.17 and 5.18 show the phase-averaged turbulence intensity modulation for case H and case H-V, respectively. The red dashed line in each figure indicates the sediment bed height. For case H (figure 5.17), the fluctuations in all directions increase around the sediment bed height during the fluid accelerating phase (\ = 0 to \ = c) while decreasing during the fluid decelerating phase (\ = c to \ = 2c). A similar trend is observed for case H-V in figure 5.18(a) and (c) for �√ D ′ D ′ and �√ | ′ | ′. However, for �√ { ′ { ′ the modulation remains at a high level during the first half of the fluid decelerating phase (\ = c to \ = 3c/2) due to the effect of vertical oscillations. As described by Manna et al. (2015), the streamwise intensity modulations behave as a wave propagating away from the wall region. In particular, the wave is generated at the beginning of the decelerating phase. Compared with the results shown in Manna 99 (a) (c)(b) Figure 5.17: Phase-averaged turbulence intensity modulation in case H. The red lines indicate the sediment bed height. (a) (b) (c) Figure 5.18: Phase-averaged turbulence intensity modulation in case H-V. The red lines indicate the sediment bed height. 100 et al. (2015), our turbulence intensity modulation exhibits a similar wave propagating pattern, but it is generated at phase \ = c/2 instead of \ = c. However, as pointed out by citetpapadopoulos2016pulsating, different cases might have different response time of the turbulent fluctuations to the oscillations, causing phase differences among the cases. For pulsating turbulent pipe flows, the streamwise velocity fluctuation modulation starts at the wall, and the maximum value of the modulation is achieved at locations around H/;+B = 2 (Papadopoulos & Vouros, 2016). In our comparison of mean veloc- ity oscillations with the Stokes solution, we observed good agreement with the Stokes solution above the sediment bed region. However, here for the turbulence velocity fluctu- ation, we observe a different pattern compared with the pulsating pipe flow results. The streamwise fluctuation modulations start at the height of H = �B43 −�, which is one sed- iment diameter below the defined sediment height. Additionally, the highest modulation happens at the sediment height location rather than at locations higher than that. 5.3.3 Sediment transport results In this section, we report results regarding the sediment transport process. Sediment transport rate According to Kidanemariam & Uhlmann (2017), the instantaneous volumetric flow rate of the particles, &?, can be calculated as: &? (C) = c�3 6!G!I #?∑ ;=1 D;? (C), (5.19) 101 Figure 5.19: Time evolution of sediment transport rate in case N. this formula sums all streamwise particle velocities multiplied by the particle volume and is divided by the product of domain length and width. To interpret its value,&? is usually normalised by a reference value, &?,A4 5 = (d ′ 6�3)1/2, where d ′ = d?/d 5 − 1 represents the relative density. Figure 5.19 shows the sediment transport rate for case N over 70 oscillation periods. The average value for the sediment transport rate is &?/&?,A4 5 = 0.006, which is small. The Shields number, (ℎ = (D2 g)/(d ′ 6�), serves as a common metric for assessing sediment particle mobility. For case N, we have (ℎ = 0.035. According to Wong & Parker (2006), the critical Shields number for the onset of sediment motion is around 0.03, indicating that a low sediment transport rate is expected in case N. 102 (a) (b) (c) (d) Case V Case H Case H-V Case H+V Figure 5.20: Time evolution of sediment transport rate in, (a) case V, (b) case H, (c) case H-V, (d) case H+V. 103 For the cases with oscillations, figure 5.20 illustrates the rhythmic nature of sediment transport rates. For each case, higher sediment transport rates are observed at the beginning, attributed to the initial arrangement of the stationary sediment bed. However, after approximately 20 periods, the sediment transport rates stabilise. For cases V and H, the sediment transport rates remain positive throughout, while cases H-V and H+V exhibit negative sediment transport rates at some times. This phenomenon is consistent with the results reported in § 5.3.1, where negative particle velocities are observed. The time-averaged sediment transport rate is computed for each case, with the results summarised in table 5.1. Compared to case N with no oscillations, the mean transport rate has increased across all cases with varying oscillations. Notably, case H+V displays a remarkable tenfold increase in the mean sediment transport rate, attributed to the simultaneous contributions of oscillations to sediment transport from both directions. Direct observation of sediment transport As defined above, the instantaneous sediment height is denoted as ℎB (x, C). Similar to the figures from Kidanemariam & Uhlmann (2014b), we plot the spatial–temporal evolution of the sediment bed interface in figures 5.21 and 5.22. The height is normalised by the particle diameter. Each figure depicts 30 periods of oscillations. In figure 5.21(a), corresponding to case N, vertical lines of the same colour are noticeable, representing stationary top particles throughout the simulation. Additionally, several yellow stripes intersect these vertical lines, depicting the movement of individual top layer particles. Figure 5.21(c), which is for case H, also shows these vertical lines. Figure 5.21(b)−(d) are characterised by vertical blocks, indicative of the pronounced rhythmic motion of sediment particles. Specifically, figure 5.21(b) (case V) showcases a 104 (a) (c) (b) (d) Case N Case V Case H Case H-V Figure 5.21: Spatial–temporal evolution of sediment bed interface in (a) case N, (b) case V, (c) case H, (d) case H-V. 105 Figure 5.22: Spatial–temporal evolution of sediment bed interface in case H+V. (a) (b) Figure 5.23: Time evolution of a the r.m.s. of sediment bed height fluctuations and b the mean wavelength _ℎ in case H+V. 106 dune-like pattern characterised by limited height variations between ridges and troughs. Figure 5.21(d), on the other hand, shows amore averaged sediment height in the simulation domain. Figure 5.22 illustrates a significant wave-like pattern formed during the oscillation process for case H+V, attributed to the high sediment transport rate. As mentioned above, due to the limited particle layer in the present simulation, the trough of the wave-like structures reaches the fixed irregular layer. It would be of interest to examine how the wave-like structures develop in case H+V. According to Kidanemariam&Uhlmann (2017), the sediment height fluctuations can be defined as: ℎ ′ 1 (G, C) = ℎB (G, C) − ℎB G (C), (5.20) the wall normal dimension of the patterns is characterised by the r.m.s. of the sediment bed, fℎ, defined as: f2 ℎ (C) = ℎ ′ 1 (G, C)ℎ′ 1 (G, C) G , (5.21) To determine the average wavelength of the pattern, a two-point correlation coefficient of the bed-height fluctuations is calculated as follows: 'ℎ (XG, C) = ℎ ′ 1 (G, C)ℎ′ 1 (G + XG, C) G f2 ℎ (C) , (5.22) the average bedform wavelength _ℎ is defined as: _ℎ = 2Gmin, (5.23) where Gmin is the global minimum of 'ℎ (XG, C). 107 Figure 5.23 shows fℎ (C) and _ℎ (C). Similar to the sediment transport rate, fℎ (C) oscillates with the bed oscillating period and shows increasing fluctuations when the bedform is moving. However, the mean wavelength is not as sensitive to the oscillation duration. At approximately C/)| = 10−40, the sediment bed maintains its fluctuation and wavelength at a relatively stable stage. However, due to the limitation of our simulation domain length, the unstable wavelength is restricted, ultimately leading the sediment bed to adopt the entire domain length !G as the wavelength. Although the waveform changes over time, figure 5.20(d) indicates that the sediment transport rate is not affected by these changes, suggesting that the sediment transport rate is more closely related to oscillation properties such as direction, amplitude and frequency. Phase-averaged sediment transport rate Figure 5.20 above shows the rhythmic nature of the sediment transport rate in various oscillations, suggesting that it would be informative to examine how sediment transport rate changes within one period. In addition to the typical calculated sediment transport rate in the streamwise direction, the rhythmic movement allows us to examine sediment movement in the vertical direction, denoted as &H,?, which is also normalised by &?,A4 5 . Although the mean of &H,? is zero, the changes within one period can provide insights into the sediment transport pattern. Figure 5.24 demonstrates the phase-averaged sediment transport rate in the streamwise and vertical directions. In each figure, the x-axis indicates the normalised time within one period, and the pink area is the standard deviation over 30 periods. For&G,?, the time mean value is plotted using the dashed line in each figure. The small standard deviation values for both directions in each case indicate that the difference between periods is 108 (a) (c) (d) (b) (e) (f) (g) (h) Figure 5.24: Time evolution of sediment transport rate in, (a) (c) (e) (g) horizontal directions, (b) (d) (f ) (h) vertical directions. 109 minimal. Notably, case H+V has the largest sediment transport rate but the smallest relative standard deviation, which again suggests that the large sediment transport rate is closely related to the bed oscillation properties. For case H with only horizontal oscillations, figure 5.24(c) shows that the peak sediment transport rate occurs in themiddle of the period with \ = c, where themaximum fluid velocity is reached (figure 5.10(a)), with nearly zero vertical transport rate as shown in figure 5.24(d). For case V and case H+V, figure 5.24(a) and (g) indicate similar sediment transport rate variation patterns during each period. Case H+V exhibits a negative overshot in the sediment transport rate (figure 5.24(g)), which occurs when the sediment particles fall back to the bottom. For vertical oscillations, the maximum sediment transport rate happens slightly earlier than the phase \ = c. Similar to case H, figures 5.24(a) and (g) show that at the time of maximum streamwise sediment transport. The vertical transport rate is approximately zero. In case H-V, where vertical and horizontal accelerations contribute to sediment trans- port at different phases, a distinct transport pattern is observed in figures 5.24(e) and (f). In this case, the value of the maximum positive and negative sediment transport rate is about the same, and the maximum negative sediment transport happens closely after the maximum positive sediment transport happens. Unlike the other cases, the maxi- mum positive streamwise sediment transport happens at the time when vertical sediment transport is negative maximum, and the maximum negative streamwise sediment trans- port happens at the time when vertical sediment transport is positive maximum. Figure 5.24(e) shows that the minimum sediment transport is observed around \ = c, even though the maximum particle velocities are observed at the same phase (figure 5.12 (b)). This discrepancy arises because, at this time, the fluid velocity is at its maximum, and 110 sampling of particles at higher locations captures the highest velocities, but the number of particles at those heights is small. Additionally, at time \ = 3/2c, vertical oscillations facilitate particle movement, but the horizontal acceleration causes a reverse flow and negative particle velocities just below the sediment bed height, resulting in a negative sediment transport rate. The above comparison between cases H-V and H+V highlights how the phase difference between horizontal and vertical oscillations can significantly impact the sediment transport rate, indicating a need for further study of different phase difference groups. 5.4 Conclusion This study presents PR-DNS simulations of sediment transport over a vibrating bed, with earthquake vibrations modelled as sinusoidal oscillations in horizontal and/or vertical directions. By analysing time-averaged and phase-averaged results for the fluid flow and sediment transport under different oscillation conditions, we have identified both similarities and differences among the various cases. The time-averaged results show that caseswith oscillations have smallermean velocity profiles and smaller velocity fluctuations than those of the case without oscillations. A comparison among different cases indicates a higher particle-induced stress for case V. The phase-averaged results reveal that the profiles of the mean velocity and velocity fluctuations do not change significantly with different phases for case V compared to the cases with horizontal oscillations. However, for all the cases with oscillations considered, the particle velocities exhibit high correlationswith the phase. In the cases with horizontal oscillations, the fluid flow velocity increases much faster than the particle velocities, 111 resulting in increasing relative velocities during the fluid acceleration phase. Comparing the horizontal oscillation cases with the pulsating turbulent flows reported in the literature, we observe distinct behaviours for different quantities. Our simulation cases fall into the current-dominated very-high-frequency regime. For the mean velocity above the sediment bed, our results show good agreement with the Stokes solution, similar to the pulsating turbulent flow in this regime. However, inside the sediment bed, we observe a slower decrease of the oscillation velocity amplitude and a continuation of phase lag. Additionally, we report the effects of the bedform. For the turbulent intensity modulation, we observe collective changes of turbulent intensity in all three directions during different phases. The streamwise turbulent intensity modulation originates at a distance of one particle diameter below the sediment bed height, with the highest modulation occurring at the sediment bed height. Sediment transport over vibrating beds has the same oscillation period, and an in- creased mean sediment transport rate is observed for all oscillation cases. Case H+V exhibits a large waveform structure; however, the sediment transport rate is insensitive to this structure and is more related to the bed oscillation properties. The phase-averaged sediment transport rate shows similar patterns for cases V, H and H+V, with the highest sediment rate coinciding with the near zero vertical sediment transport rate. In contrast, case H-V shows a distinct phase-averaged sediment transport rate pattern, indicating the importance of phase differences between different oscillation components in the sediment transport process. Further investigations are needed to provide amore comprehensive picture of sediment transport during seismic events. For instance, we have shown the different transport patterns with oscillations in vertical and horizontal directions when the phase difference 112 is either zero or c. However, more phase angles should be studied to complete the comparison. Additionally, our simulation considers coarse grains, but the seabed can also be formed by fine grains and mud; the sediment transport for those conditions should not be overlooked. Furthermore, the oscillation frequencies in our simulation are the same in both directions. It would be interesting to examine different frequencies in different directions and their effects on sediment transport patterns. Overall, our simulations are an initial step towards a deeper understanding of sed- iment transport under seismic-like oscillations, revealing the critical role of oscillation in different directions. These insights contribute to a better understanding of sediment dynamics in natural environments affected by seismic activities and can inform future studies and applications in civil engineering and environmental management. 113 Chapter 6 Conclusing Remarks 6.1 Contribution of the thesis In this study, we have explored particle dynamics across a broad spectrum of topics, offering insights into the intricate interactions between particles and fluid flow. We introduced a novel method for tracking coherent particle clusters within turbulent flows. Utilizing this method, we conducted the first study on the temporal coherence of inertial particle clusters in turbulence, discovering a strong correlation between the size of a cluster and its lifespan. Our findings indicate that larger clusters predominantly emerge from the recombination of other substantial clusters and, upon dispersal, their components generally integrate into new, large, and long-lived clusters. We also demonstrated that a cluster’s size remains relatively constant throughout its existence, with the clusters residing in regions of the fluid characterized by low enstrophy and strain rates. Our investigations into sediment transport involved different particle shapes in open channel flow, where distinct behaviors of various particles were noted, and bedform 114 generation was observed in all cases. Despite similar mean velocity profiles across the simulations, we found notable differences near the sediment bed interface due to the particle shape impact on sediment packing. We noted increased Reynolds normal stresses, particularly in the wall-normal and spanwise directions for cases with oblate particles. We confirmed that spherical particles predominantly exhibit rolling motions within the sediment bed, whereas oblate particles tend to slide. These dynamics were consistent even in mixed-particle cases, underscoring that particle shape significantly influences both bedform and sediment transport. To study particle remobilization during seismic events, we performed studies of particle transport over a vibrating bed under varied oscillation conditions. We have observed increased sediment transport rates in all simulation cases, and we reported the time-average and phase-average results for fluid flow and particle movement. Comparing our results with those from pulsating channel flow, we found good agreement with the Stokes solution above the sediment bed, though phase lag increased within the sediment bed. Turbulent modulation was also examined. Our research highlights the rhythmic nature of sediment transport under oscillations and analyzes varying sediment transport rate behaviors across different cases. 6.2 Future studies 6.2.1 GPU-accelerated simulations The computational cost of particle-resolved simulations is substantial, particularly in scenarios involving a large number of particles. Mazzuoli et al. (2019) reported that their simulation required 10 million CPU hours and spanned approximately one and a half 115 years. Similarly, Vittori et al. (2020) used 60 million CPU hours for a simulation that ran for two years. In our case, typical simulations require between 300,000 to 800,000 CPU hours, presenting significant challenges in exploring a broader range of cases and scenarios due to these high computational demands. Recently, the use of GPUs has emerged as a promising solution to accelerate CFD simulations. Weymouth & Font (2023) reported speed increases of 10 to 100 times with GPU acceleration. Our research group has observed similar improvements, indicating that GPU technology holds substantial promise for enhancing the efficiency of particle- resolved simulations. 6.2.2 Expanding particle diversity in simulations In most particle-resolved simulations, a single particle type is typically modeled within each case. However, in the simulation detailed in § 4, the case mix includes two different particle types, expanding the complexity of the modeled system. Despite this enhance- ment, real-world sediment transport processes involve a significantly broader variety of particle types. Particle sizes in natural settings can vary widely, ranging from a few mi- crometers to several millimeters. Exploring the impacts of this size diversity on sediment transport could provide deeper insights into the dynamics of these systems. Furthermore, studying the combined effects of particle shape and size on sediment transport could yield crucial understanding relevant to environmental engineering and geophysical research. 116 6.2.3 Sediment dynamics under varied oscillation conditions In our simulation detailed in § 5, we observed a waveform generated when vertical and horizontal oscillations were applied in-phase, resulting in a significant increase in sedi- ment transport rate. Conversely, when the phase angle between the two directions was c, the sediment bed interface remained relatively flat, accompanied by a relatively smaller increase in sediment transport rate. Currently, our simulations utilize sinusoidal oscil- lations at uniform frequencies. Future studies should explore more realistic oscillation conditions. Specifically, examining a broader range of phase angles could provide a com- prehensive understanding of how these oscillations interact. Additionally, advancements in computational capabilities may allow for simulations over larger domains, enabling the complete capture of large bedforms and ensuring that the simulations are not constrained by boundary effects. 117 Bibliography Abhyankar, S., Brown, J., Constantinescu, E. M., Ghosh, D., Smith, B. F. & Zhang, H. 2018 Petsc/ts: A modern scalable ode/dae solver library. arXiv preprint arXiv:1806.01437 . Adams, J. 1990 Paleoseismicity of the cascadia subduction zone: Evidence from tur- bidites off the oregon-washington margin. Tectonics 9 (4), 569–583. Aliseda, A., Cartellier, A., Hainaux, F. & Lasheras, J. C. 2002 Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence. J. Fluid Mech. 468, 77–105. Allen, J. R. 2012 Principles of physical sedimentology. Springer Science & Business Media. Ardekani, M. N., Costa, P., Breugem, W. P. & Brandt, L. 2016 Numerical study of the sedimentation of spheroidal particles. Int. J. Multiph. Flow 87, 16–34. Ardekani, M. N., Costa, P., Breugem, W.-P., Picano, F & Brandt, L. 2017 Drag reduction in turbulent channel flow laden with finite-size oblate spheroids. J. Fluid Mech. 816, 43–70. 118 Ashi, J., Sawada, R., Omura, A. & Ikehara, K. 2014 Accumulation of an earthquake- induced extremely turbid layer in a terminal basin of the nankai accretionary prism. Earth Planets Space 66 (1), 1–9. Atwater, B. F., Carson, B., Griggs, G. B., Johnson, H.P. & Salmi, M.S. 2014 Rethinking turbidite paleoseismology along the cascadia subduction zone. Geology 42 (9), 827–830. Aussillous, P., Chauchat, J., Pailha, M., Médale, M. & Guazzelli, E. 2013 Inves- tigation of the mobile granular layer in bedload transport by laminar shearing flows. J. Fluid Mech. 736, 594–615. Baker, L., Frankel, A., Mani, A. & Coletti, F. 2017 Coherent clusters of inertial particles in homogeneous turbulence. J. Fluid Mech. 833, 364–398. Baker, L. J. & Coletti, F. 2022 Experimental investigation of inertial fibres and disks in a turbulent boundary layer. J. Fluid Mech. 943, A27. Bec, J., Biferale, L., Cencini, M., Lanotte, A., Musacchio, S. & Toschi, F. 2007 Heavy particle concentration in turbulence at dissipative and inertial scales. Phys. Rev. Lett. 98 (8), 084502. Bec, J., Biferale, L., Lanotte, A. S., Scagliarini, A. & Toschi, F. 2010 Turbulent pair dispersion of inertial particles. J. Fluid Mech. 645, 497–528. Biegert, E., Vowinckel, B. & Meiburg, E. 2017 A collision model for grain-resolving simulations of flows over dense, mobile, polydisperse granular sediment beds. J. Com- put. Phys. 340, 105–127. 119 Blois, G., Best, J. L., Sambrook S., Gregory H. & Hardy, R. J. 2014 Effect of bed permeability and hyporheic flow on turbulent flow over bed forms. Geophys. Res. Lett. 41 (18), 6435–6442. Boyer, F., Guazzelli, É. & Pouliquen, O. 2011 Unifying suspension and granular rheology. Phys. Rev. Lett. 107 (18), 188301. Bragg, A. D. & Collins, L. R. 2014 New insights from comparing statistical theories for inertial particles in turbulence: I. Spatial distribution of particles. New J. Phys. 16 (5), 055013. Bragg, A. D., Ireland, P. J. & Collins, L. R. 2016 Forward and backward in time dispersion of fluid and inertial particles in isotropic turbulence. Phys. Fluids 28 (1), 013305. Brandt, L. & Coletti, F. 2022 Particle-laden turbulence: progress and perspectives. Annu. Rev. Fluid Mech. 54 (1), 159–189. Brereton, G. J. & Hwang, J. L. 1994 The spacing of streaks in unsteady turbulent wall-bounded flow. Phys. Fluids 6 (7), 2446–2454. Breugem, W. P. 2012 A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 231 (13), 4469–4498. Calzavarini, E., Kerscher, M., Lohse, D. & Toschi, F. 2008 Dimensionality and morphology of particle and bubble clusters in turbulent flow. J. Fluid Mech. 607, 13–24. 120 Carbone, M., Bragg, A. D. & Iovieno, M. 2019 Multiscale fluid–particle thermal interaction in isotropic turbulence. J. Fluid Mech. 881, 679–721. Carter, D. W. & Coletti, F. 2018 Small-scale structure and energy transfer in homoge- neous turbulence. J. Fluid Mech. 854, 505–543. Cassel, M., Lavé, J., Recking, A., Malavoi, J.-R. & Piégay, H. 2021 Bedload transport in rivers, size matters but so does shape. Sci. Rep. 11 (1), 508. Challabotla, N. R., Zhao, L. & Andersson, H. I. 2015a Orientation and rotation of inertial disk particles in wall turbulence. J. Fluid Mech. 766, R2. Challabotla, Niranjan Reddy, Zhao, Lihao & Andersson, Helge I 2015b Shape effects on dynamics of inertia-free spheroids in wall turbulence. Phys. Fluids 27 (6). Chang, C., Qiao, F., Bo, J., Peng, D. & Li, Q. 2023 The influence of seismic frequency spectrum on the instability of loess slope. Sci. Rep. 13 (1), 10949. Chang, K., Malec, B. J. & Shaw, R. A. 2015 Turbulent pair dispersion in the presence of gravity. New J. Phys. 17 (3), 033010. Cheng, Z., Jelly, T. O., Illingworth, S. J., Marusic, I. & Ooi, A. S. H. 2020 Forcing frequency effects on turbulence dynamics in pulsatile pipe flow. Int. J. Heat Fluid Fl. 82, 108538. Costa, P., Boersma, B. J., Westerweel, J. & Breugem, W. P. 2015 Collision model for fully resolved simulations of flows laden with finite-size particles. Phys. Rev. E. 92 (5), 053012. 121 Derksen, J. J. 2011 Simulations of granular bed erosion due to laminar shear flow near the critical shields number. Phys. Fluids 23 (11). Dey, S. & Ali, S. Z. 2017 Mechanics of sediment transport: Particle scale of entrainment to continuum scale of bedload flux. J. Eng. Mech. 143 (11), 04017127. DiBenedetto, M. H., Ouellette, N. T. & Koseff, J. R. 2018 Transport of anisotropic particles under waves. J. Fluid Mech. 837, 320–340. Durán, O., Claudin, P. & Andreotti, B. 2011 On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws. Aeolian Res. 3 (3), 243–270. Eaton, J. K. & Fessler, J. R. 1994 Preferential concentration of particles by turbulence. Int. J. Multiph. Flow 20, 169–209. Elghobashi, S. 1994 On predicting particle-laden turbulent flows. Appl. Sci. Res. 52, 309–329. Fairbridge, R. W. & Bourgeois, J. 2006 Encyclopedia of Sedimentology. Springer Berlin, Heidelberg. Ferenc, J. S. & Néda, Z. 2007 On the size distribution of poisson voronoi cells. Physica A 385 (2), 518–526. Fiabane, L., Zimmermann, R., Volk, R., Pinton, J.-F. & Bourgoin, M. 2012 Clustering of finite-size particles in turbulence. Phys. Rev. E 86 (3), 035301. Fong, K. O., Amili, O. & Coletti, F. 2019 Velocity and spatial distribution of inertial particles in a turbulent channel flow. J. Fluid Mech. 872, 367–406. 122 Fornari, W., Ardekani, M. N. & Brandt, L. 2018 Clustering and increased settling speed of oblate particles at finite reynolds number. J. Fluid Mech. 848, 696–721. Fröhlich, K., Meinke, M. & Schröder, W. 2020 Correlations for inclined prolates based on highly resolved simulations. J. Fluid Mech. 901, A5. Gao, Q., Deane, G. & Shen, L. 2021 Bubble production by air filament and cavity breakup in plunging breaking wave crests. J. Fluid Mech. 929, A44. Ge, J. & Monroe, C. A. 2019 The effect of coefficient of restitution in modeling of sand granular flow for core making: Part i free-fall experiment and theory. Int. J. of Metalcast. 13, 753–767. Geyer, W. R., Woodruff, J. D. & Traykovski, P. 2001 Sediment transport and trapping in the hudson river estuary. Estuaries 24, 670–679. Ghodke, C. D. & Apte, S. V. 2016 Dns study of particle–bed–turbulence interactions in an oscillatory wall-bounded flow. J. Fluid Mech. 792, 232–251. Glowinski, R., Pan, T. W., Hesla, T. I., Joseph, D. D. & Periaux, J. 2001 A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow. J. Comput. Phys. 169 (2), 363–426. Goldfinger, C., Morey, A. E., Nelson, C. H., Gutierrez-Pastor, J., Johnson, J. E., Karabanov, E., Chaytor, J., Eriksson, A. & Party, S. S. 2007 Rupture lengths and temporal history of significant earthquakes on the offshore and north coast segments of the northern san andreas fault based on turbidite stratigraphy. Earth Planet Sci Lett. 254 (1-2), 9–27. 123 Goldfinger, C., Nelson, C. H., Morey, A. E., Johnson, J. E., Patton, J. R., Kara- banov, E. B., Gutierrez-Pastor, J., Eriksson, A. T., Gracia, E., Dunhill, G., Enkin, R. J., Dallimore, A. & Vallier, T. 2012 Turbidite event history–methods and implications for holocene paleoseismicity of the cascadia subduction zone. Tech. Rep.. US Geological Survey. Gondret, P., Lance, M. & Petit, L. 2002 Bouncing motion of spherical particles in fluids. Phys. Fluids 14 (2), 643–652. Goto, S. & Vassilicos, J. C. 2008 Sweep-stick mechanism of heavy particle clustering in fluid turbulence. Phys. Rev. Lett. 100 (5), 054503. Hardin, D. P. & Saff, E. B. 2005Minimal riesz energy point configurations for rectifiable d-dimensional manifolds. Adv. Math. 193 (1), 174–204. He, S., Liu, H. & Shen, L. 2022 Simulation-based study of turbulent aquatic canopy flows with flexible stems. J. Fluid Mech. 947, A33. Ikehara, K., Usami, K. & Kanamatsu, T 2020Repeated occurrence of surface-sediment remobilization along the landward slope of the japan trench by great earthquakes.Earth Planets Space 72 (1), 1–9. Ireland, P. J., Bragg, A. D. & Collins, L. R. 2016 The effect of reynolds number on inertial particle dynamics in isotropic turbulence. Part 1. Simulations without gravita- tional effects. J. Fluid Mech. 796, 617–658. Jain, R., Tschisgale, S. & Fröhlich, J. 2019 A collision model for dns with ellipsoidal particles in viscous fluid. Int. J. Multiph. Flow 120, 103087. 124 Jain, R., Tschisgale, S. & Fröhlich, J. 2021 Impact of shape: Dns of sediment transport with non-spherical particles. J. Fluid Mech. 916, A38. Jain, R., Vowinckel, B. & Fröhlich, J. 2017 Spanwise particle clusters in dns of sediment transport over a regular and an irregular bed. Flow, Turbul. Combust. 99, 973–990. Jeffery, G. B. 1922 The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. Lond. A 102 (715), 161–179. Jelly, T. O., Chin, R. C., Illingworth, S. J., Monty, J. P., Marusic, I. & Ooi, A. S. H. 2020 A direct comparison of pulsatile and non-pulsatile rough-wall turbulent pipe flow. J. Fluid Mech. 895, R3. Ji, C., Munjiza, A., Avital, E., Xu, D. & Williams, J. 2014 Saltation of particles in turbulent channel flow. Phys. Rev. E. 89 (5), 052202. Kempe, T. & Fröhlich, J. 2012 Collision modelling for the interface-resolved simulation of spherical particles in viscous fluids. J. Fluid Mech. 709, 445–489. Kempe, T., Schwarz, S. & Fröhlich, J. 2009 Modelling of spheroidal particles in vis- cous flows. In Proceedings of the Academy Colloquium Immersed Boundary Methods: Current Status and Future Research Directions (KNAW, Amsterdam, The Netherlands, 15–17 June 2009), , vol. 845. Kidanemariam, A. G. & Uhlmann, M. 2014a Direct numerical simulation of pattern formation in subaqueous sediment. J. Fluid Mech. 750, R2. 125 Kidanemariam, A. G. & Uhlmann, M. 2014b Interface-resolved direct numerical simu- lation of the erosion of a sediment bed sheared by laminar channel flow. Int. J. Multiph. Flow 67, 174–188. Kidanemariam, A. G. & Uhlmann, M. 2017 Formation of sediment patterns in channel flow: minimal unstable systems and their temporal evolution. J. Fluid Mech. 818, 716–743. Kim, D. & Choi, H. 2006 Immersed boundary method for flow around an arbitrarily moving body. J. Comput. Phys. 212 (2), 662–680. Kim, J. & Moin, P. 1985 Application of a fractional-step method to incompressible navier–stokes equations. J. Comput. Phys. 59 (2), 308–323. Kok, J. F., Parteli, E. J. R., Michaels, T. I. & Karam, D. B. 2012 The physics of wind-blown sand and dust. Rep. Prog. Phys. 75 (10), 106901. Krumbein, W. C. 1942 Settling-velocity and flume-behavior of non-spherical particles. Eos, Trans. Am. Geophys. Union 23 (2), 621–633. Lian, H., Chang, X. Y. & Hardalupas, Y. 2019 Time resolved measurements of droplet preferential concentration in homogeneous isotropic turbulence without mean flow. Phys. Fluids 31 (2), 025103. Liu, H., He, S., Shen, L. & Hong, J. 2021 Simulation-based study of covid-19 outbreak associated with air-conditioning in a restaurant. Phys. Fluids 33 (2). Lodahl, C. R., Sumer, B. M. & Fredsøe, J. 1998 Turbulent combined oscillatory flow and current in a pipe. J. Fluid Mech. 373, 313–348. 126 Lozano-Durán, A. & Jiménez, J. 2014 Time-resolved evolution of coherent structures in turbulent channels: characterization of eddies and cascades. J. Fluid Mech. 759, 432–471. Lucci, F., Ferrante, A. & Elghobashi, S. 2010 Modulation of isotropic turbulence by particles of taylor length-scale size. J. Fluid Mech. 650, 5–55. Manna, M., Vacca, A. & Verzicco, R. 2012 Pulsating pipe flow with large-amplitude oscillations in the very high frequency regime. part 1. time-averaged analysis. J. Fluid Mech. 700, 246–282. Manna, M., Vacca, A. & Verzicco, R. 2015 Pulsating pipe flow with large-amplitude oscillations in the very high frequency regime. part 2. phase-averaged analysis. J. Fluid Mech. 766, 272–296. Mao, Z.-X. & Hanratty, T. J. 1986 Studies of the wall shear stress in a turbulent pulsating pipe flow. J. Fluid Mech. 170, 545–564. Marchioli, C., Fantoni, M. & Soldati, A. 2010 Orientation, distribution, and deposi- tion of elongated, inertial fibers in turbulent channel flow. Phys. Fluids 22 (3). Marchioli, C., Zhao, L. & Andersson, H. I. 2016 On the relative rotational motion between rigid fibers and fluid in turbulent channel flow. Phys. Fluids 28 (1). Mathai, V., Huisman, S. G., Sun, C., Lohse, D. & Bourgoin, M. 2018 Dispersion of air bubbles in isotropic turbulence. Phys. Rev. Lett. 121 (5), 054501. Mathai, V., Lohse, D. & Sun, C. 2020 Bubble and buoyant particle laden turbulent flows. Annu. Rev. Condens. Matter Phys. 11. 127 Maxey, M. R. 1987 The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. J. Fluid Mech. 174, 441–465. Mazzuoli, M., Kidanemariam, A. G. & Uhlmann, M. 2019 Direct numerical simula- tions of ripples in an oscillatory flow. J. Fluid Mech. 863, 572–600. McHugh, C. M., Kanamatsu, T., Seeber, L., Bopp, R., Cormier, M. H. & Usami, K. 2016 Remobilization of surficial slope sediment triggered by the ad 2011 mw 9 tohoku-oki earthquake and tsunami along the japan trench. Geology 44 (5), 391–394. Merritt, W. S., Letcher, R. A. & Jakeman, A. J. 2003 A review of erosion and sediment transport models. Environ. Model. Softw. 18 (8-9), 761–799. Moernaut, J., Van Daele, M., Strasser, M., Clare, M. A., Heirman, K., Viel, M., Cardenas, J., Kilian, R., de Guevara, B. L., Pino, M., R., Urrutia & De Batist, M. 2017 Lacustrine turbidites produced by surficial slope sediment remobilization: a mechanism for continuous and sensitive turbidite paleoseismic records. Mar. Geol. 384, 159–176. Mohd-Yusof, J. 1996 Interaction of massive particles with turbulence. PhD thesis. Cornell University. Momenifar, M. & Bragg, A. D. 2020 Local analysis of the clustering, velocities, and accelerations of particles settling in turbulence. Phys. Rev. Fluids 5 (3), 034306. Monchaux, R., Bourgoin, M. & Cartellier, A. 2010 Preferential concentration of heavy particles: A Voronoï analysis. Phys. Fluids 22 (10), 103304. 128 Monchaux, R., Bourgoin, M. & Cartellier, A. 2012 Analyzing preferential con- centration and clustering of inertial particles in turbulence. Int. J. Multiph. Flow 40, 1–18. Monchaux, R. & Dejoan, A. 2017 Settling velocity and preferential concentration of heavy particles under two-way coupling effects in homogeneous turbulence. Phys. Rev. Fluids 2 (10), 104302. Mordant, N. & Pinton, J.-F. 2000 Velocity measurement of a settling sphere. Eur. Phys. J. B 18, 343–352. Mortensen, P. H., Andersson, H. I., Gillissen, J. J. J. & Boersma, B. J. 2008Dynamics of prolate ellipsoidal particles in a turbulent channel flow. Phys. Fluids 20 (9). Nikora, V.and McEwan, I., McLean, S., Coleman, S., Pokrajac, D. & Walters, R. 2007 Double-averaging concept for rough-bed open-channel and overland flows: Theoretical background. J. Hydraul. Eng. 133 (8), 873–883. Oujia, T., Matsuda, K. & Schneider, K. 2020 Divergence and convergence of inertial particles in high reynolds number turbulence, arXiv: 2005.00525. Papadopoulos, P. K. & Vouros, A. P. 2016 Pulsating turbulent pipe flow in the current dominated regime at high and very-high frequencies. Int. J. Heat Fluid Fl. 58, 54–67. Parishani, H., Ayala, O., Rosa, B., Wang, L.-P. & Grabowski, W. W. 2015 Effects of gravity on the acceleration and pair statistics of inertial particles in homogeneous isotropic turbulence. Phys. Fluids 27 (3), 033304. 129 Patil, A. & Fringer, O. 2022 Drag enhancement by the addition of weak waves to a wave–current boundary layer over bumpy walls. J. Fluid Mech. 947, A3. Peskin, C. S. 1982 The fluid dynamics of heart valves: experimental, theoretical, and computational methods. Annu. Rev. Fluid Mech. 14, 235–259. Petersen, A. J., Baker, L. & Coletti, F. 2019 Experimental study of inertial particles clustering and settling in homogeneous turbulence. J. Fluid Mech. 864, 925–970. Prosser, I. P., Rutherfurd, I. D., Olley, J. M., Young, W. J., Wallbrink, P. J & Moran, C. J/ 2001 Large-scale patterns of erosion and sediment transport in river networks, with examples from australia. Mar. Freshwater Res. 52 (1), 81–99. Rana, S., Anderson, W. & Day, M. 2021 An entrainment paradox: How hysteretic salta- tion and secondary transport augment atmospheric uptake of aeolian source materials. J. Geophys. Res. 126 (10), e2020JD033493. Roma, A. M., Peskin, C. S. & Berger, M. J. 1999 An adaptive version of the immersed boundary method. J. Comput. Phys. 153 (2), 509–534. Saff, E. B. & Kuijlaars, A. B. J. 1997 Distributing many points on a sphere. Math. Intell. 19, 5–11. Sawford, B. L., Yeung, P.-K. & Borgas, M. S. 2005 Comparison of backwards and forwards relative dispersion in turbulence. Phys. Fluids 17 (9), 095109. Scherer, M., Uhlmann, M., Kidanemariam, A. G. & Krayer, M. 2022 On the role of turbulent large-scale streaks in generating sediment ridges. J. Fluid Mech. 930, A11. 130 Schmeeckle, M. W. 2014 Numerical simulation of turbulence and sediment transport of medium sand. J. Geophys. Res. 119 (6), 1240–1262. Schmeeckle, M. W., Nelson, J. M., Pitlick, J. & Bennett, J. P. 2001 Interparticle collision of natural sediment grains in water. Water Resour. Res. 37 (9), 2377–2391. Scotti, A. & Piomelli, U. 2001 Numerical simulation of pulsating turbulent channel flow. Phys. Fluids 13 (5), 1367–1384. Shin, M. & Koch, D. L. 2005 Rotational and translational dispersion of fibres in isotropic turbulent flows. J. Fluid Mech. 540, 143–173. Squires, K. D. & Eaton, J. K. 1991 Preferential concentration of particles by turbulence. Phys. Fluid A 3 (5), 1169–1178. Sun, R. & Xiao, H. 2016 Sedifoam: A general-purpose, open-source cfd–dem solver for particle-laden flow with emphasis on sediment transport. Comput. Geosci. 89, 207–219. Sundaram, S. & Collins, L. R. 1997 Collision statistics in an isotropic particle-laden turbulent suspension. Part 1. Direct numerical simulations. J. FluidMech. 335, 75–109. Tagawa, Y., Mercado, J., Prakash, V. N., Calzavarini, E., Sun, C. & Lohse, D. 2012 Three-dimensional Lagrangian Voronoï analysis for clustering of particles and bubbles in turbulence. J. Fluid Mech. 693, 201–215. Tardu, S. F. & Binder, G. 1993 Wall shear stress modulation in unsteady turbulent channel flow with high imposed frequencies. Phys. Fluids 5 (8), 2028–2037. 131 Tardu, S. F., Binder, G. & Blackwelder, R. F. 1994 Turbulent channel flow with large-amplitude velocity oscillations. J. Fluid Mech. 267, 109–151. Tosi, P., Sbarra, P. & De Rubeis, V. 2012 Earthquake sound perception. Geophys. Res. lett. 39 (24). Tschisgale, S., Kempe, T. & Fröhlich, J. 2018 A general implicit direct forcing im- mersed boundary method for rigid particles. Comput. Fluids 170, 285–298. Tu, S. W. & Ramaprian, B. R. 1983 Fully developed periodic turbulent pipe flow. part 1. main experimental results and comparison with predictions. J. Fluid Mech. 137, 31–58. Uhlmann, M. 2005 An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209 (2), 448–476. Uhlmann, M. 2008 Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime. Phys. Fluids 20 (5). Uhlmann, M. & Doychev, T. 2014 Sedimentation of a dilute suspension of rigid spheres at intermediate Galileo numbers: the effect of clustering upon the particle motion. J. Fluid Mech. 752, 310–348. Valance, A., Rasmussen, K. R., El Moctar, A. O. & Dupont, P. 2015 The physics of aeolian sand transport. C. R. Phys. 16 (1), 105–117. Veeramani, C., Minev, P. D. & Nandakumar, K. 2009 Collision modeling between two non-brownian particles in multiphase flow. Int. J. Therm. Sci. 48 (2), 226–233. 132 Vittori, G., Blondeaux, P., Mazzuoli, M., Simeonov, J. & Calantoni, J. 2020 Sediment transport under oscillatory flows. Int. J. Multiph. Flow 133, 103454. Voermans, J. J., Ghisalberti, M. & Ivey, G. N. 2017 The variation of flow and turbulence across the sediment–water interface. J. Fluid Mech. 824, 413–437. Voth, G. A. & Soldati, A. 2017 Anisotropic particles in turbulence. Annu. Rev. Fluid Mech. 49 (1), 249–276. Vowinckel, B., Kempe, T. & Fröhlich, J. 2014 Fluid–particle interaction in turbulent open channel flow with fully-resolved mobile beds. Adv. Water Resour. 72, 32–44. Wang, G., Fong, K. O., Coletti, F., Capecelatro, J. & Richter, D. H. 2019 Inertial particle velocity and distribution in vertical turbulent channel flow: a numerical and experimental comparison. Int. J. Multiph. Flow 120, 103105. Wang, L.-P. & Stock, D. E. 1993 Dispersion of heavy particles by turbulent motion. J. Atmos. Sci. 50 (13), 1897–1913. Weng, C., Boij, S. & Hanifi, A. 2016 Numerical and theoretical investigation of pulsatile turbulent channel flows. J. Fluid Mech. 792, 98–133. Weymouth, G. D. & Font, B. 2023 Waterlily. jl: A differentiable fluid simulator in julia with fast heterogeneous execution. arXiv preprint arXiv:2304.08159 . Wong, M. & Parker, G. 2006 Reanalysis and correction of bed-load relation of meyer- peter and müller using their own database. J. Hydraul. Eng. 132 (11), 1159–1168. Yang, Z., Deng, B.-Q. & Shen, L. 2018 Direct numerical simulation of wind turbulence over breaking waves. J. Fluid Mech. 850, 120–155. 133 Yang, Z., Lu, X.-H., Guo, X., Liu, Y. & Shen, L. 2017 Numerical simulation of sediment suspension and transport under plunging breaking waves. Comput. Fluids 158, 57–71. Zamansky, R., Coletti, F., Massot, M. & Mani, A. 2016 Turbulent thermal convection driven by heated inertial particles. J. Fluid Mech. 809, 390–437. Zedler, Emily A & Street, Robert L 2006 Sediment transport over ripples in oscillatory flow. J. Hydraul. Eng. 132 (2), 180–193. Zhao, L, Challabotla, N. R., Andersson, H. I. & Variano, E. A. 2015 Rotation of nonspherical particles in turbulent channel flow. Phys. Fluids 115 (24), 244501. Zheng, X., Feng, S. & Wang, P. 2021 Modulation of turbulence by saltating particles on erodible bed surface. J. Fluid Mech. 918, A16. Zheng, X. J., Huang, N. & Zhou, Y.-H. 2003 Laboratory measurement of electrification of wind-blown sands and simulation of its effect on sand saltation movement. J. Geophys. Res. 108 (D10). Zhou, Z.-Y., Zou, R.-P., Pinson, D. & Yu, A.-B. 2011Dynamic simulation of the packing of ellipsoidal particles. Ind. Eng. Chem. Res. 50 (16), 9787–9798. Zhu, Z., Hu, R., Lei, Y., Shen, L. & Zheng, X. 2022 Particle resolved simulation of sediment transport by a hybrid parallel approach. Int. J. Multiph. Flow 152, 104072. 134 Acknowledgements Dedication Abstract List of Tables List of Figures Introduction Motivations Thesis overview Numerical method Immersed boundary method for simulation of particles Collision model for particles Validation cases Particle cluster identification Life and death of inertial particle clusters in turbulence Introduction Numerical Cases Identification and tracking of clusters Results Lifetime Birth and death Relative dispersion of clustered particles The role of the turbulent activity Conclusions Interface resolved simulation of particles with different shapes in sediment transport Introduction Simulation details Numerical method Simulation set-up Sediment height and pressure gradient Simulation initialization Results Fluid flow properties Particle properties Shear stress decomposition Conclusions Particle-resolved simulation of sediment transport over vibrating beds Introduction Simulation details Reference frame Governing equation and computational frameworks Simulation setup Flow property decomposition, sediment height, and Reynolds number Results Basic statistics of the flow and particles Comparison with pulsating turbulent flows Sediment transport results Conclusion Conclusing Remarks Contribution of the thesis Future studies GPU-accelerated simulations Expanding particle diversity in simulations Sediment dynamics under varied oscillation conditions Bibliography