We investigate and model the initiation of motion of a single particle on a structured substrate within an oscillatory boundary layer flow, following a mechanistic approach. By deterministically relating forces and torques acting on the particle to the instantaneous ambient flow, the effects of flow unsteadiness are captured, revealing rich particle dynamics. Laboratory experiments in an oscillatory flow tunnel characterise the initiation and early stages of motion, with particle imaging velocimetry measurements yielding the flow conditions at the motion threshold. The experiments validate and complement results from particle-resolved direct numerical simulations, combining an immersed boundary method with a discrete element method that incorporates a static friction contact model. Within the parameter range just above the motion threshold, the mobile particle rolls without sliding over the substrate, indicating that motion initiation is governed by an unbalanced torque rather than a force. Both experimental and numerical results show excellent agreement with an analytical torque balance including hydrodynamic torque derived from the theoretical Stokes velocity profile, and contributions of lift, added mass, and externally imposed pressure gradient. In addition to static and rolling particle states, we identify a wiggling regime where the particle moves but does not leave its original pocket. Our deterministic approach enables prediction of the phase within the oscillation cycle at which the particle starts moving, without relying on empirical threshold estimates, and can be extended to a wide range of flow and substrate conditions, as long as turbulence is absent and interactions with other mobile particles are negligible.

We proposed a novel diffused interface approach for the interface normal and curvature calculation in the conservative level set method framework for both Newtonian and non-Newtonian multiphase flows. The standard benchmark of reversible vortex problem is used to test the interface capturing method for both standard and conservative level set method. In addition, the new approach is validated to accurately simulate oscillating droplet, bubble rising in a viscoelastic fluid, and droplet impact in a deep pool, in which a good agreement is obtained with analytical solutions and experimental measurements.

The dynamics of entrainment of finite-size heavy particles in a turbulent open channel flow over a smooth surface are analyzed. Three types of simulations, namely with freely moving, rotation-constrained, and spanwise-motion-constrained particles, were conducted using particle-resolved direct numerical simulations. With the aid of a relative velocity suitably defined in the vicinity of the finite-size particle, we decompose the hydrodynamic force into drag and lift contributions and evaluate the local wall-normal shear rate around the particles. By means of coherent structure eduction techniques, we investigate flow structures before and during lift-off events. Rotation-constrained simulations revealed the insignificance of particle rotation in the entrainment mechanism. Spanwise-motion-constrained simulations revealed the importance of particle location with respect to flow structures with apparent changes in entrainment frequency, duration of the entrainment process, wall-normal shear around the particles, and distance to the nearest vortical structures during lift-off. The contribution of lift to the wall-normal force is found to be responsible for the initiation of particle entrainment, which is induced by a high-shear event associated with fast-moving fluid. The presence of quasi-streamwise vortices is shown to be an important ingredient for the entrainment of particles into the bulk flow. The results show that, at marginal Shields number values, a high wall-normal shear rate and the proximity of an intense quasi-streamwise vortex are essential elements of the entrainment mechanism.

Particle-resolved direct numerical simulation (PR-DNS) is performed for turbulent open channel flow over a smooth horizontal wall with a vertical cylinder and a dilute set of mobile, heavy, spherical particles. At the chosen parameter point (which matches a previous study without a cylinder) the particles are mostly translating in the horizontal plane while remaining in contact with the wall. It is shown that the presence of the cylinder leads to the generation of intense vortical structures, enhanced turbulence intensity in the wake region, and to strong modifications of the local wall shear stress. These cylinder-induced perturbations have direct consequences for the average particle concentration: preferential accumulation/depletion in different parts of the wake region occurs, while the wall-normal transport of particles (against gravity) is significantly enhanced. A second simulation which adds roughness elements on the wall reveals an additional effect upon the wall-normal distribution of particles. It turns out that the configuration with wall-roughness and a wall-mounted cylinder features the largest fraction of entrained particles, even far from the wall.

We analyse the long-time dynamics of trajectories within the stability boundary between laminar and turbulent square duct flow. If not constrained to a symmetric subspace, the edge trajectories exhibit a chaotic dynamics characterised by a sequence of alternating quiescent phases and intense bursting episodes. The dynamics reflects the different stages of the well-known near-wall streak-vortex interaction. Most of the time, the edge states feature a single streak with a number of flanking vortices attached to one of the four surrounding walls. The initially straight streak undergoes a linear instability and eventually breaks in an intense bursting event. At the same time, the downstream vortices give rise to a new low-speed streak at one of the neighbouring walls, thereby causing the turbulent activity to 'switch' from one wall to the other. If the edge dynamics is restricted to a single or twofold mirror-symmetric subspace, the bursting and wall-switching episodes become self-recurrent in time, representing the first periodic orbits found in square duct flow. In contrast to the chaotic edge states in the non-symmetric case, the imposed symmetries enforce analogue bursting cycles to simultaneously appear at two parallel opposing walls in a mirror-symmetric configuration. Both the localisation of turbulent activity to one or two walls and the wall-switching dynamics are shown to be common phenomena in marginally turbulent duct flows. We argue that such episodes represent transient visits of marginally turbulent trajectories to some of the edge states detected here.

Turbulent open channel flow is known to feature a multi-layer structure near the free surface. In the present work we employ direct numerical simulations considering Reynolds numbers up to $$\boldsymbol{R{e_\tau } = 900}$$ and domain sizes large enough ( $$\boldsymbol{{L_x} = 12\pi h}$$ , $$\boldsymbol{{L_z} = 4\pi h}$$ ) to faithfully capture the effect of very-large-scale motions in order to test the proposed scaling laws and ultimately answer the question: How far does the influence of the free surface extend? In the region near the free surface, where fluctuation intensities of velocity and vorticity become highly anisotropic, we observe the previously documented triple-layer structure, consisting of a wall-normal velocity damping layer that scales with the channel height $$h$$ , and two sublayers that scale with the near-surface viscous length scale $$\boldsymbol{{\ell _{\boldsymbol{V}}} = {\boldsymbol{Re}}_{\boldsymbol{b}}^{ - 1/2}h}$$ and with the Kolmogorov length scale $$\boldsymbol{{\ell _{\boldsymbol{K}}} = {\boldsymbol{Re}}_{\boldsymbol{b}}^{ - 3/4}h}$$ , respectively. The scaling laws previously proposed by Calmet and Magnaudet [J. Fluid. Mech. 474 , 355–378 (2003)] are found to hold with the following exceptions. The thin layer, where the intensity of surface-parallel components of the vorticity rapidly decreases to zero, is here found to scale with the Kolmogorov length scale $$\boldsymbol{{\ell _{\boldsymbol{K}}}}$$ rather than with the near-surface viscous scale $$\boldsymbol{{\ell _{\boldsymbol{V}}}}$$ . In addition, we argue that the Kolmogorov length scale is the relevant scale for the mean velocity gradient near the free surface. Both the mean velocity gradient and the fluctuation intensity of the surface-parallel component of vorticity decay to zero in the Kolmogorov sublayer $$\boldsymbol{{\delta _{\boldsymbol{K}}} \approx 20{\ell _{\boldsymbol{K}}}}$$ . On the other hand, the layer, where the wall-normal turbulence intensity decreases linearly to zero near the free surface, scales with $$\boldsymbol{{\ell _{\boldsymbol{V}}}}$$ rather than $$\boldsymbol{{\ell _{\boldsymbol{K}}}}$$ as suggested by Calmet and Magnaudet. The corresponding near-surface viscous sublayer measures $$\boldsymbol{{\delta _{\boldsymbol{V}}} \approx {\ell _{\boldsymbol{V}}}}$$ . Importantly, the streamwise turbulence intensity profile for $$\boldsymbol{{\boldsymbol{R}}{{\boldsymbol{e}}_\tau } \geq 400}$$ suggests that the influence of the free-slip boundary penetrates essentially all the way down to the solid wall through the appearance of enhanced very-large-scale motions ( $$\boldsymbol{{\delta _{{\boldsymbol{SIL}}}} \approx h}$$ ). In contrast, the layer where the surface-normal turbulence intensity is damped to zero is restricted to the free surface ( $$\boldsymbol{{\delta _{{\boldsymbol{NVD}}}} \approx 0.3h}$$ ). As a consequence, the partitioning of the surface-influenced region has to be expanded to a four-layer structure that spans the entire channel height $$\boldsymbol{h}$$ .

We analyse the dynamics within the stability boundary between laminar and turbulent square duct flow with the aid of an edge-tracking algorithm. As for the circular pipe, the edge state turns out to be a chaotic attractor within the edge if the flow is not constrained to a symmetric subspace. The chaotic edge state dynamics is characterised by a sequence of alternating quiescent phases and regularly occurring bursting episodes. These latter reflect the different stages of the well-known streak-vortex interaction in near-wall turbulence: The edge states feature most of the time a single streak with a number of flanking quasi-streamwise vortices attached to one of the four surrounding walls. The initially straight streak undergoes the classical linear instability and eventually breaks in an intense bursting event due to the action of the quasi-streamwise vortices. At the same time, the vortices give rise to a new generation of low-speed streaks at one of the neighbouring walls, thereby causing the turbulent activity to `switch' from one wall to the other. When restricting the edge dynamics to a single or twofold mirror-symmetric subspace, on the other hand, the outlined bursting and wall-switching episodes become self-recurrent in time. These edge states thus represent the first periodic orbits found in the square duct. In contrast to the chaotic edge states in the non-symmetric case, the imposed symmetries enforce analogue bursting cycles to simultaneously appear at two parallel opposing walls in a mirror-symmetric configuration. Both localisation of the turbulent activity to one or two walls and wall switching are shown to be a common phenomenon in low Reynolds number duct turbulence. We therefore argue that the marginally turbulent trajectories transiently visit the identified edge states during these episodes, so that the edge states become actively involved in the turbulent dynamics.

The present work extends the direct-forcing immersed boundary method introduced by Garcia-Villalba et al. (2023), broadening its application from spherical to arbitrarily-shaped particles, while maintaining its capacity to address both neutrally-buoyant and light objects (down to a density ratio of 0.5). The proposed method offers a significant advantage over existing methods regarding its simplicity, in particular for the case of neutrally- buoyant particles. Three test cases from the literature are selected for validation: a neutrally-buoyant prolate spheroid in a shear flow; a settling oblate spheroid; and, finally, a rising oblate spheroid.

The long-time behaviour of trajectories embedded in the edge manifold between the laminar and turbulent basins of attraction is investigated in order to identify simple invariant solutions to the Navier-Stokes equations in straight ducts with square cross-section. With the aid of the iterative bisection procedure of Toh and Itano ( J. Fluid Mech. , vol. 481, 2003, pp. 67-76), we detect three travelling waves that represent attracting ‘edge states’ within the separatrix and saddles w.r.t. the full state space. All three travelling waves are presumably connected to solutions known in the community. In other situations, the edge states to which trajectories are attracted seem to be invariant periodic cycles that feature a rich dynamics: Along each cycle, the flow oscillates periodically between two flow states, in each of which streamwise vorticity is predominantly residing near one set of parallel walls, while the vortical activity associated to the other pair of opposing walls is markedly reduced. This ‘switching’ dynamics is accompanied by intermittent bursting-like high-dissipation excursions. While we cannot rigorously prove that the identified periodic cycle is a periodic orbit rather than a finite-precision approximation to a heteroclinic cycle, sound arguments for the former supposition are given. Interestingly, the different phases of the periodic cycle qualitatively resemble a similar ‘switching’ behaviour observed in (chaotic) marginal turbulence, which suggests that the found periodic cycles could help understanding the characteristic dynamics of duct turbulence at low Reynolds numbers.

We conducted direct numerical simulations of turbulent open channel flow (OCF) and closed channel flow (CCF) of friction Reynolds numbers up to Re_τ≈ 900 in large computational domains up to L_x× L_z=12π h × 4π h to analyse the Reynolds number scaling of turbulence intensities. Unlike CCF, our data suggests that the streamwise turbulence intensity in OCF scales with the bulk velocity for Re_τ≳ 400 . The additional streamwise kinetic energy in OCF with respect to CCF is provided by larger and more intense very-large-scale motions in the former type of flow. Therefore, compared to CCF, larger computational domains of L_x× L_z=12π h× 4π h are required to faithfully capture very-large-scale motions in OCF—and observe the reported scaling. OCF and CCF turbulence statistics data sets are available at https://doi.org/10.4121/88678f02-2a34-4452-8534-6361fc34d06b .

In the present chapter we focus on the fundamentals of non-grid-conforming numerical approaches to simulating particulate flows, implementation issues and grid convergence vs. available reference data. The main idea is to avoid adapting the mesh (and - as much as possible - the discrete operators) to the time-dependent fluid domain with the aim to maximize computational efficiency. We restrict our attention to spherical particle shapes (while deviations from sphericity are treated in a subsequent chapter). We show that similar ideas can be successfully implemented in a variety of underlying fluid flow solvers, leading to powerful tools for the direct numerical simulation of large particulate systems.

We discuss computational methods to model collisions and lubrication at the particle scale in particle-laden flows. We review the well-established approaches to compute the Lagrangian tracking of spherical and non-spherical rigid particles with collisional and lubrication-dominated interactions. This should equip the readers with a variety of efficient methods to choose from depending on the specific features of the problem they are interested in. While non-spherical particle collisions and lubrication between a pair of spheres are (reasonably) well understood, much less is known about lubrication between two non-spherical particles. We also provide the readers with assorted ways to speed up the implementation of the collision (also referred to as granular dynamics) solver in a computer program. We put the numerical models to the test in simple and more challenging flow configurations. Finally, we summarize the challenges and open questions lying ahead to improve the modeling of short-range collisional and hydrodynamic particle–particle interactions in particle-laden flows.

We present a simple modification of the direct-forcing immersed boundary method (IBM) proposed by Uhlmann [J. Comput. Phys, 2005] in order to enable it to be applied to particulate flows with solid-to-fluid density ratios around unity. The main difference with respect to the original formulation lies in the particle velocity update which is performed directly after the preliminary velocity field has been computed in the absence of any IBM volume forcing term. In addition, we apply the forcing term to the entire space occupied by the immersed solid object (instead of to the vicinity of its interface only). The present approach requires the evaluation of integrals of the velocity field over the volume occupied by the solid particle, which are evaluated efficiently as sums over the respective quantities available at particle-attached force points. The resulting method can be used seamlessly for density ratios down to ρ_p/ρ_f>0.5. The new formulation has been validated using three configurations: (i) lateral migration of a neutrally buoyant circular particle in two-dimensional Couette flow; (ii) the release from rest of a neutrally buoyant sphere in a free stream; (iii) the release of a particle in a free stream after an initial phase in which it is translationally fixed with an imposed angular velocity. In all three test cases the present IBM formulation yields a very good agreement with the available reference data. Thus, the proposed approach is a cost-efficient and accurate modification of the original method which allows for the simulation of fluid systems involving density-matched solid particles.

We review some of the results obtained to date with the aid of the PR-DNS approach to turbulent particulate flows. It is shown that the method has matured to a point which allows to apply it successfully to a wide variety of fluid/particle configurations, albeit still at a relatively large computational cost. Due to the availability of high-fidelity space-and-time-resolved data, a number of challenging open questions have already been addressed in unprecedented detail.

We have performed particle-resolved direct numerical simulations of many heavy non-spherical particles settling under gravity in the dilute regime. The particles are oblate spheroids of aspect ratio 1.5 and density ratio 1.5. Two Galileo numbers are considered, namely 111 and 152, for which a single oblate spheroid follows a steady vertical and a steady oblique path, respectively. In both cases, a strongly inhomogeneous spatial distribution of the disperse phase in the form of columnar clusters is observed, with a significantly enhanced average settling velocity as a consequence. Thus, in contrast to previous results for spheres, the qualitative difference in the single particle regime does not result in a qualitatively different behavior of the many-particle cases. In addition, we have carried out an analysis of pairwise interactions of particles in the well-known drafting-kissing-tumbling setup, for oblate spheroids of aspect ratio 1.5 and for spheres. We have varied systematically the relative initial position between the particle pair and we have considered free-to-rotate particles and rotationally-locked ones. We have found that the region of attraction for both particle shapes, with and without rotation, is very similar. However, significant differences occur during the drafting and tumbling phases. In particular, free-to-rotate spheres present longer drafting phases and separate quickly after the collision. Spheroids remain close to each other for longer times after the collision, and free-to-rotate ones experience two or more collision events. Therefore, we have observed a shape-induced increase in the interaction which might explain the increased tendency to cluster of the many-particle cases.

Ice crystal number concentrations were often found to be orders of magnitude higher than the number concentration of ice nucleating particles; a finding that indicated the presence of secondary ice production (SIP). Although 6 mechanisms of SIP have been both discovered and theorized, it is still not fully understood and the recent studies have been inconclusive in identifying the dominant process in real conditions. This lack of constraint of ice multiplication adds to the uncertainty of cloud simulations in climate models. Studying SIP is challenging due to the various interfering factors involved.In this study, we attempt to further our knowledge in understanding the SIP mechanisms using two different but complementary approaches. The first consists of using remote sensing tools such as Himawari-8 and MODIS retrievals in addition to the SOCRATES in-situ data to identify the presence of SIP and categorize the possible mechanism involved.The second approach utilizes numerical simulations to further understand these mechanisms that are potentially responsible for SIP, but through the study of the dynamics of the different particles (ice crystals, supercooled droplets, graupel...) involved in these processes. In this approach we focus on the characteristics of the particles, such as their diameter and concentration, as well as the presence of turbulence, that are crucial in describing their movement and the feasibility of the mechanisms under study.

We present direct numerical simulation results of turbulent open channel flow at bulk Reynolds numbers up to 12000, coupled with (passive) scalar transport at Schmidt numbers up to 200. Care is taken to capture the very large scale motions which appear already for relatively modest Reynolds numbers. The transfer velocity at the flat, free surface is found to scale with the Schmidt number to the power "-1/2", in accordance with previous studies and theoretical predictions for uncontaminated surfaces. The scaling of the transfer velocity with Reynolds number is found to vary, depending on the Reynolds number definition used. To compare the present results with those obtained in other systems, we define a turbulent Reynolds number at the edge of the surface-influenced layer. This allows us to probe the two-regime model of Theofanous [Turbulent mass transfer at free, gas-liquid interfaces, with applications to open-channel, bubble and jet flows. Int. J. Heat Mass Transfer 19, 613--624, 1976], which is found to correctly predict that small-scale vortices significantly affect the mass transfer for turbulent Reynolds numbers larger than 500. It is further established that the root-mean-square of the surface divergence is, on average, proportional to the mean transfer velocity. However, the spatial correlation between instantaneous surface divergence and transfer velocity tends to decrease with increasing Schmidt number and increase with increasing Reynolds number. The latter is shown to be caused by an enhancement of the correlation in high-speed regions, which in turn is linked to the spatial distribution of surface-parallel vortices.

We have numerically investigated the turbulent flow and sediment grain motion in an open-channel flow configuration over a subaqueous sediment bed featuring two-dimensional transverse ripples at moderate Reynolds number and super-critical Shields number values. The simulation data, which were generated by means of particle-resolved direct numerical simulation, are the same as in our previous work (Kidanemariam & Uhlmann, J. Fluid Mech., vol. 818, 2017, pp. 716–743). By carefully choosing the computational box sizes, we were able to accommodate single ripple units which form over an initially flat sediment bed at a wavelength equal to the domain length. The ripples then evolve into their asymmetric shape relatively quickly and eventually migrate downstream steadily while maintaining their shape and size. In the present study, using a ripple-conditioned phase-averaging procedure, we are able to obtain novel insights into the evolution of the turbulent flow and particle motion over the bedforms, in particular the spatial structure of the basal shear stress and its relation to the particle flow rate. Our analysis confirms that the boundary shear-stress maximum is located upstream of the ripple crest, while the particle flow rate is essentially in phase with the ripple topology, with an average phase difference between the two in the range of 18–19 particle diameters for the considered parameter values. We were further able to confirm the link between the sediment flux relaxation behaviour and the observed shear-stress/geometry lag, by direct evaluation of the saturation length scale.

The question whether or not the presence of warm hydrometeors in clouds may play a significant role in the nucleation of new ice particles has been debated for several decades. While the early works of Fukuta and Lee (1986) and Baker (1991) indicated that it might be irrelevant, the more recent study of Prabhakaran et al. (2019) [arXiv:1906.06129] suggested otherwise. In this work, we are aiming to quantify the ice-nucleating potential using high-fidelity flow simulation techniques around a single hydrometeor and use favorable considerations to upscale the effects to a collective of ice particles in clouds. While we find that ice nucleation may be enhanced in the vicinity of a warm hydrometeor by several orders of magnitude and that the affected volume of air is much larger than previously estimated, it is very unlikely that this effect alone causes the rapid enhancement of ice nucleation observed in some types of clouds, mainly due to the low initial volumetric ice concentration. Nonetheless, it is suggested to implement this effect into existing cloud models in order to investigate second-order effects such as ice nucleus preactivation or enhancement after the onset of glaciation.

We have performed spectral/spectral-element simulations of a single oblate spheroid with small geometrical aspect ratio settling in an unbounded ambient fluid, for a range of Galileo numbers covering the various regimes of motion (steady vertical, steady oblique, vertical periodic and chaotic). The high-fidelity data provided includes particle quantities (statistics in the chaotic case), as well as flow profiles and pressure maps. The reference data can be used as an additional benchmark for other numerical approaches, where a careful grid convergence study for a specific target parameter point is often useful. We further describe an extension of a specific immersed boundary method (Uhlmann, J. Comput. Phys, 209(2):448–476, 2005) to enable the tracking of non-spherical particles. Finally, the reference cases are computed with this immersed boundary method at various spatial and temporal resolutions, and grid convergence is discussed over the various regimes of spheroidal particle motion. The cross-validation results can serve as a guideline for the design of simulations with the aid of similar non-conforming methods, involving spheroidal particles with Galileo numbers of O(100).