A Comparative Study of Turbulence in Ramp-Up and Ramp-Down Unsteady Flows

DNS of a turbulent channel flow subjected to a step change in pressure gradient are performed to facilitate a direct comparison between ramp-up and ramp-down flows. Strong differences are found between behaviours of turbulence in the two flows. The wall shear stress in the ramp-up flow first overshoots, and then strongly undershoots the quasi-steady value in the initial stage of the excursion, before approaching the quasi-steady value. In a strongly decelerating flow, the wall shear stress tends to first undershoot but then overshoot the quasi-steady value. ??Slow?? response of turbulence as well as flow inertia is responsible for these behaviours. In the ramp-up flow, the response of turbulence is similar to that observed in uniformly accelerating flows from previous studies, exhibiting a three-stage development. However, the transition between the various stages is more gradual and the responding stage is much longer and slower in the flows considered here. It has been shown that the delay in the near wall region is longer than that in the buffer layer confirming that turbulence response first occurs at the location of peak turbulence production. In a strongly decelerating flow, the response of turbulence exhibits a two-stage development. In both ramp-up and ramp down flows, the energy distribution in the three components of turbulent kinetic energy deviates from that of the steady flow. In a ramp-up flow, more energy is in $u_1^\prime $ and less in $u_2^\prime $ and $u_3^\prime $ , whereas the trend is reversed in a ramp-down flow. This is a reflection of the redistribution of turbulence from $u_1^\prime $ to $u_2^\prime $ and $u_3^\prime $ .     

Towards the Development of an Evolution Equation for Flame Turbulence Interaction in Premixed Turbulent Combustion

Flame turbulence interaction is one of the leading order terms in the scalar dissipation \(\left (\widetilde {\varepsilon }_{c}\right )\) transport equation [35] and is thus an important phenomenon in premixed turbulent combustion. Swaminathan and Grout [36] and Chakraborty and Swaminathan [15, 16] have shown that the effect of strain rate on the transport of \(\widetilde {\varepsilon }_{c}\) is dominated by the interaction between the fluctuating scalar gradients and the fluctuating strain rate, denoted here by \(\overline {\rho }\widetilde {\Delta }_{c}= \overline {\rho {\alpha }\nabla c^{\prime \prime }S_{ij}^{\prime \prime }\nabla c^{\prime \prime }}\) ; this represents the flame turbulence interaction. In order to obtain an accurate representation of this phenomenon, a new evolution equation for \(\widetilde {\Delta }_{c}\) has been proposed. This equation gives a detailed insight into flame turbulence interaction and provides an alternative approach to model the important physics represented by \(\widetilde {\Delta }_{c}\) . The \(\widetilde {\Delta }_{c}\) evolution equation is derived in detail and an order of magnitude analysis is carried out to determine the leading order terms in the \(\widetilde {\Delta }_{c}\) evolution equation. The leading order terms are then studied using a Direct Numerical Simulation (DNS) of premixed turbulent flames in the corrugated flamelet regime. It is found that the behaviour of \(\widetilde {\Delta }_{c}\) is determined by the competition between the source terms (pressure gradient and the reaction rate), diffusion/dissipation processes, turbulent strain rate and the dilatation rate. Closures for the leading order terms in \(\widetilde {\Delta }_{c}\) evolution equation have been proposed and compared with the DNS data.     

Simultaneous measurement of velocity and pressure in a plane jet

A new technique was developed for the simultaneous measurement of velocity and pressure in turbulent flows. To accomplish this objective, a new probe (hereafter called the combined probe) that consists of an X-type hot-wire probe and a newly devised pressure probe was developed. The pressure probe was miniaturized by the MEMS fabrication process and by using a 0.1-in. microphone as a pressure sensor for improving the spatial resolution. This pressure probe was placed between two hot-wire sensors of which the X-type hot-wire probe was composed. The pressure probe was given a hemispherical tip, like that of a pitot tube, because an earlier pressure probe with a conical tip suffered from a reduction in spatial resolution. The spatial arrangement of the pressure probe and the hot-wire probe for the combined probe was carefully determined, because there was a risk that the measurement accuracy of one probe will be influenced by disturbances caused by the other probe when the two probes were placed very close to each other. Therefore, the combined probe was arranged to engender no noticeable interference between the velocity data and the pressure data measured by their respective probes. As one application of this combined probe, simultaneous measurements of pressure and two components of instantaneous velocity were performed in a plane jet. The turbulent energy budget and the cross-correlation coefficient of velocity and pressure in the intermittent region of the plane jet were estimated. The results show that the mean streamwise velocity, velocity fluctuation, and pressure fluctuation profiles were consistent with those measured individually using the X-type hot-wire probe or pressure probe. Moreover, it was shown that the integral value of the diffusion term (which should theoretically be equal to zero) in the turbulent energy transport equation was closer to zero than previous reports (Bradbury in J Fluid Mech 23(Part 1):31–64, 1965). In addition, the time variation of the cross-correlation coefficient in the intermittent region supports the vortex structure model predicted in previous studies (Browne et?al. in J Fluid Mech 149:355–373, 1984; Tanaka et?al. JSME Int J Ser B 49(4):899–905, 2006; Sakai et?al. J Fluid Sci Technol 2(3):611–622, 2007).     

Hypervelocity shock standoff on spheres in air

To provide data for the validation of computational fluid dynamics models, measurements of the shock standoff distance on spheres in hypervelocity flows have been made. Test flows of air at 8.7 and 9.7 km/s were generated in the X2 expansion tunnel fitted with a Mach 10 nozzle. High-speed video images were analysed with a least-squares shape-fitting algorithm. Assuming a spherical shock shape near the nose enabled increased resolution measurements beyond the native pixel size. Normalised shock standoff distances, $\Delta $ / $D$ , in the range 0.03–0.04 were measured, with sphere diameters, $D$ , of 40, 60 and 80 mm.     

An Updated Portrait of Transition to Turbulence in Laminar Pipe Flows with Periodic Time Dependence (A Correlation Study)

Transition to turbulence in axially symmetrical laminar pipe flows with periodic time dependence classified as pure oscillating and pulsatile (pulsating) ones is the concern of the paper. The current state of art on the transitional characteristics of pulsatile and oscillating pipe flows is introduced with a particular attention to the utilized terminology and methodology. Transition from laminar to turbulent regime is usually described by the presence of the disturbed flow with small amplitude perturbations followed by the growth of turbulent bursts. The visual treatment of velocity waveforms is therefore a preferred inspection method. The observation of turbulent bursts first in the decelerating phase and covering the whole cycle of oscillation are used to define the critical states of the start and end of transition, respectively. A correlation study referring to the available experimental data of the literature particularly at the start of transition are presented in terms of the governing periodic flow parameters. In this respect critical oscillating and time averaged Reynolds numbers at the start of transition; Re _os,crit and Re _ta,crit are expressed as a major function of Womersley number, $\sqrt {\omega ^\prime } $ defined as dimensionless frequency of oscillation, f. The correlation study indicates that in oscillating flows, an increase in Re _os,crit with increasing magnitudes of $\sqrt {\omega ^\prime } $ is observed in the covered range of $1<\sqrt {\omega ^\prime } <72$ . The proposed equation (Eq. 7), ${\rm{Re}}_{os,crit} ={\rm{Re}}_{os,crit} \left( {\sqrt {\omega ^\prime } } \right)$ , can be utilized to estimate the critical magnitude of $\sqrt {\omega ^\prime }$ at the start of transition with an accuracy of ±12?% in the range of $\sqrt {\omega ^\prime } <41$ . However in pulsatile flows, the influence of $\sqrt {\omega ^\prime }$ on Re _ta,crit seems to be different in the ranges of $\sqrt {\omega ^\prime } <8$ and $\sqrt {\omega ^\prime } >8$ . Furthermore there is rather insufficient experimental data in pulsatile flows considering interactive influences of $\sqrt {\omega ^\prime } $ and velocity amplitude ratio, A ₁. For the purpose, the measurements conducted at the start of transition of a laminar sinusoidal pulsatile pipe flow test case covering the range of 0.21<?A ₁?<0.95 with $\sqrt {\omega ^\prime } <8$ are evaluated. In conformity with the literature, the start of transition corresponds to the observation of first turbulent bursts in the decelerating phase of oscillation. The measured data indicate that increase in $\sqrt {\omega ^\prime } $ is associated with an increase in Re _ta,crit up to $\sqrt {\omega ^\prime } =3.85$ while a decrease in Re _ta,crit is observed with an increase in $\sqrt {\omega ^\prime } $ for $\sqrt {{\omega }'} >3.85$ . Eventually updated portrait is pointing out the need for further measurements on i) the end of transition both in oscillating and pulsatile flows with the ranges of $\sqrt {\omega ^\prime } <8$ and $\sqrt {\omega ^\prime } >8$ , and ii) the interactive influences of $\sqrt {\omega ^\prime } $ and A ₁ on Re _ta,crit in pulsatile flows with the range of $\sqrt {\omega ^\prime } >8$ .     

The Scalar Gradient Alignment Statistics of Flame Kernels and its Modelling Implications for Turbulent Premixed Combustion

The effects of mean flame curvature on reaction progress variable gradient, $\nabla c$ , alignment with local turbulent strain rate are studied based on three-dimensional Direct Numerical Simulation (DNS) data of turbulent premixed flame kernels with different initial radii under decaying turbulence. A statistically planar flame is also considered in order to compare the results obtained from the kernels with a flame of zero mean curvature. It is found that the dilatation rate effects diminish with decreasing kernel radius due to defocusing of heat in the positively curved regions. This gives rise to a decrease in the extent of reaction progress variable gradient alignment with most extensive principal strain rate with decreasing kernel radius. The modelling implications of the statistics of the alignment of $\nabla c$ with local strain rate have been studied in terms of scalar dissipation rate transport. A new modelling methodology for the contribution of the scalar-turbulence interaction term in the transport equation for the mean scalar dissipation is suggested addressing the reduced effects of dilatation rate for flame kernels and the diminished value of turbulent straining at the small length scales at which turbulence interacts with small flame kernels. The performance of the new models is found to be satisfactory while comparing to DNS results. The existing models for the dilatation contribution and the combined chemical reaction and molecular dissipation contributions to the transport of mean scalar dissipation, which were originally proposed for statistically planar flames, are found to satisfactorily predict the corresponding quantities for turbulent flame kernels.     

The influence of gas phase velocity fluctuations on primary atomization and droplet deformation

The effects of grid-generated velocity fluctuations on the primary atomization and subsequent droplet deformation of a range of laminar liquid jets are examined using microscopic high-speed backlit imaging of the break-up zone and laser Doppler anemometry of the gas phase separately. This is done for fixed gas mean flow conditions in a miniature wind tunnel experiment utilizing a selection of fuels, turbulence-generating grids and two syringe sizes. The constant mean flow allows for an isolated study of velocity fluctuation effects on primary atomization in a close approximation to homogeneous decaying turbulence. The qualitative morphology of the primary break-up region is examined over a range of turbulence intensities, and spectral analysis is performed in order to ascertain the break-up frequency which, for a case of no grid, compares well with the existing literature. The addition of velocity fluctuations tends to randomize the break-up process. Slightly downstream of the break-up region, image processing is conducted in order to extract a number of metrics, which do not depend on droplet sphericity, and these include droplet aspect ratio and orientation, the latter quantity being somewhat unconventional in spray characterization. A turbulent Weber number $We^{\prime}$ which takes into account gas phase fluctuations is utilized to characterize the resulting droplet shapes, in addition to a mean Weber number <We _d>. Above a $We^{\prime}>0.05$ a clear positive relationship exists between the mean aspect ratio of droplets and the turbulent Weber number where $We^{\prime}$ is varied by altering all relevant variables including the velocity root mean square, the initial droplet diameter, the surface tension and the density.     

A spectral chart method for estimating the mean turbulent kinetic energy dissipation rate

We present an empirical but simple and practical spectral chart method for determining the mean turbulent kinetic energy dissipation rate $ \left\langle \varepsilon \right\rangle $ in a variety of turbulent flows. The method relies on the validity of the first similarity hypothesis of Kolmogorov (C R (Doklady) Acad Sci R R SS, NS 30:301–305, 1941) (or K41) which implies that spectra of velocity fluctuations scale on the kinematic viscosity ν and $ \left\langle \varepsilon \right\rangle $ at large Reynolds numbers. However, the evidence, based on the DNS spectra, points to this scaling being also valid at small Reynolds numbers, provided effects due to inhomogeneities in the flow are negligible. The methods avoid the difficulty associated with estimating time or spatial derivatives of the velocity fluctuations. It also avoids using the second hypothesis of K41, which implies the existence of a ?5/3 inertial subrange only when the Taylor microscale Reynods number R _λ is sufficiently large. The method is in fact applied to the lower wavenumber end of the dissipative range thus avoiding most of the problems due to inadequate spatial resolution of the velocity sensors and noise associated with the higher wavenumber end of this range.The use of spectral data (30?≤?R _λ?≤?400) in both passive and active grid turbulence, a turbulent mixing layer and the turbulent wake of a circular cylinder indicates that the method is robust and should lead to reliable estimates of $ \left\langle \varepsilon \right\rangle $ in flows or flow regions where the first similarity hypothesis should hold; this would exclude, for example, the region near a wall.     

A two-parameter model of turbulence,and its application to free jets

A model of turbulence is investigated in which the Reynolds stress appearing in the momentum equation is calculated from the expression \(\overline {u'\upsilon '} = - \sqrt k L(\partial U/\partial y)\) ; the kinetic energy, k, and the length scale,L, of turbulence are determined from differential transport equations for these quantities. These equations are solved for various free-jet situations, and the empirical constants involved are adjusted so as to give best agreement between predictions and experimental results. The plane mixing layer, the plane jet and the radial jet can be predicted with a single set of constants; for the round jet a different set has to be used. This suggests a dependence of the otherwise universal constants on the ratio ofL to the radiusr. Comparisons are presented of predicted and measured rates of spread, profiles forU, k and \(\overline {u'\upsilon '} \) , and energy balances. For most cases the agreement is within the experimental accuracy.     

Wall-Driven Versus Pressure Gradient-Driven Flows in Porous Media

The heat transfer characteristics of two boundary layer flows past an isothermal plane surface adjacent to a saturated Darcy–Brinkman porous medium is compared to each other in this paper. The flows are driven either by a stretching of the adjacent plane boundary, or by an external pressure gradient. It is found that below a threshold value $\tilde{P}r_{*} $ of the modified Prandtl number $\tilde{P}r$ , the Nussselt number in case of the pressure gradient-driven flow is larger than in case of the wall- driven flow, while for $\tilde{P}r>\tilde{P}r_{*} $ the flow driven by the moving wall provides a more efficient heat transfer mechanism. The dependence of $\tilde{P}r_{*} $ on the Darcy number is also discussed in detail.     

Plane jet excited by disturbances with spanwise phase variations

Evolution of the near-field structures of a plane jet excited by temporal periodic disturbances with spanwise phase variations was investigated with stereoscopic particle image velocimetry. The three-dimensional vorticity distributions were reconstructed by using Taylor’s frozen field hypothesis. When ?, the temporal phase difference of disturbances in the spanwise direction was π; chain-link-fence type structures were formed. The $\Uplambda$ vortices in the chain-link-fence structures were then distorted into an $\Upomega$ shape, and the head of the vortex was detached and reconnects to form a vortex ring, or reconnects to the adjacent V-shaped vortices to form an A-shaped vortex. After the reconnection stage, the flow field was occupied by uniformly distributed fine scale eddies. Here, the overall turbulent kinetic energy and shear stress were suppressed, and the jet width was narrower than that of the unexcited case and other forced cases. In the case of ? = π/2, spanwise rollers and rib structures were formed near the nozzle exit after the first vortex pairing. However, further vortex pairing did not occur downstream, and the rate at which the jet widened was reduced.     

Accounting for Polydispersion in the Eulerian Large Eddy Simulation of the Two-Phase Flow in an Aeronautical-type Burner

A major issue for the simulation of two-phase flows in engines concerns the modeling of the liquid disperse phase, either in the Lagrangian or the Eulerian approach. In the perspective of massively parallel computing, the Eulerian approach seems to be more suitable, as it uses the same algorithms as the gaseous phase solver. However taking into account the whole physics of a turbulent spray, especially in terms of polydispersity, requires an additional modeling effort. The Mesoscopic Eulerian Formalism (MEF) [13] accounts for the effect of turbulence on the disperse phase, and was extended to the Large Eddy Simulation framework [41], but is limited to monodisperse flows. In [38], the influence of polydispersity on resolved and unresolved turbulent motions of the disperse phase was highlighted, and a first model was proposed, based on size-conditioned statistics. Starting from this idea, a coupling between the MEF and the Multifluid Approach (MA) [30] is proposed. The MA decomposes the Eulerian phase into several fluid classes called sections, and corresponding to size intervals. Each section uses then size-conditioned closures. The original idea of this work is to use the MEF closures independently in each section, taking into account the mean droplet size of this section. This new approach, called Multifluid Mesoscopic Eulerian Formalism (MMEF), is then able to capture polydispersion with associated size-conditioned turbulent dynamics. First, the importance of polydispersity and the ability of MMEF to capture it are highlighted with a 0D evaporation case and a 2D vortex case, showing its impact on dynamics in both size and physical spaces. Then, the MMEF is applied to the MERCATO configuration of ONERA [18]. Results are compared to monodisperse Eulerian, Lagrangian and experimental results.     

Dynamics of Isotropic Homogeneous Turbulence with Linear Forcing Using a Lattice Boltzmann Method

The Lattice Boltzmann Method (LBM) has proved to be a promising approach to solve the Navier–Stokes equations, especially for incompressible and isothermal cases. For turbulent flows, the quality of the predictions has been previously studied considering standard spectral forced (ten Cate et al., Comput Fluids 35:1239–1251, 2006) statistically homogeneous isotropic turbulence. In the present contribution, a recently proposed linear forcing scheme working in physical space (Lundgren 2003; Rosales and Meneveau, Phys Fluids 17(9):095106–1,8, 2005) has been integrated in a three-dimensional fifteen-velocity LBM formulation. Results have been analyzed, with special attention to the dynamics of the flow through the invariants of the velocity tensor. This topic had not been studied yet for the linear forcing, regardless of the nature (spectral or LBM) of the numerical method. Results fully agree with standard pseudo-spectral direct numerical simulations, results proving the validity of the LBM with linear forcing in real space to study this kind of turbulent flows.     

Multi-component Reynolds-averaged Navier–Stokes simulations of Richtmyer–Meshkov instability and mixing induced by reshock at different times

Turbulent mixing generated by shock-driven acceleration of a perturbed interface is simulated using a new multi-component Reynolds-averaged Navier–Stokes (RANS) model closed with a two-equation $K$ – $\epsilon $ model. The model is implemented in a hydrodynamics code using a third-order weighted essentially non-oscillatory finite-difference method for the advection terms and a second-order central difference method for the gradients in the source and diffusion terms. In the present reshocked Richtmyer–Meshkov instability and mixing study, an incident shock with Mach number $M\!a_{\mathrm{s}}=1.20$ is generated in air and progresses into a sulfur hexafluoride test section. The time evolution of the predicted mixing layer widths corresponding to six shock tube test section lengths are compared with experimental measurements and three-dimensional multi-mode numerical simulations. The mixing layer widths are also compared with the analytical self-similar power-law solution of the simplified model equations prior to reshock. A set of model coefficients and initial conditions specific to these six experiments is established, for which the widths before and after reshock agree very well with experimental and numerical simulation data. A second set of general coefficients that accommodates a broader range of incident shock Mach numbers, Atwood numbers, and test section lengths is also established by incorporating additional experimental data for $M\!a_{\mathrm{s}}=1.24$ , $1.50$ , and $1.98$ with $At=0.67$ and $M\!a_{\mathrm{s}}=1.45$ with $At=-0.67$ and previous RANS modeling. Terms in the budgets of the turbulent kinetic energy and dissipation rate equations are examined to evaluate the relative importance of turbulence production, dissipation and diffusion mechanisms during mixing. Convergence results for the mixing layer widths, mean fields, and turbulent fields under grid refinement are presented for each of the $M\!a_{\mathrm{s}}=1.20$ cases.     

Development of a feedback model for the high-speed impinging planar jet

This article experimentally investigates the self-excited impinging planar jet flow, specifically the development and propagation of large-scale coherent flow structures convecting between the nozzle lip and the downstream impingement surface. The investigation uses phase-locked particle image velocimetry measurements and a new structure-tracking scheme to measure convection velocity and characterize the impingement mechanism near the plate, in order to develop a new feedback model that can be used to predict the oscillation frequency as a function of flow velocity ( $U_o$ ), impingement distance ( $x_o$ ) and nozzle thickness ( $h$ ). The resulting model prediction shows a good agreement with experimental tone frequency data.     

Statistical behavior of post-shock overpressure past grid turbulence

When a shock wave ejected from the exit of a 5.4-mm inner diameter, stainless steel tube propagated through grid turbulence across a distance of 215 mm, which is 5–15 times larger than its integral length scale \(L_{u}\) , and was normally incident onto a flat surface; the peak value of post-shock overpressure, \(\Delta P_{\mathrm{peak}}\) , at a shock Mach number of 1.0009 on the flat surface experienced a standard deviation of up to about 9 % of its ensemble average. This value was more than 40 times larger than the dynamic pressure fluctuation corresponding to the maximum value of the root-mean-square velocity fluctuation, \(u^{\prime }= 1.2~\hbox {m}/\hbox {s}\) . By varying \(u^{\prime }\) and \(L_{u}\) , the statistical behavior of \(\Delta P_{\mathrm{peak}}\) was obtained after at least 500 runs were performed for each condition. The standard deviation of \(\Delta P_{\mathrm{peak}}\) due to the turbulence was almost proportional to \(u^{{\prime }}\) . Although the overpressure modulations at two points 200 mm apart were independent of each other, we observed a weak positive correlation between the peak overpressure difference and the relative arrival time difference.     

Computational study of a supersonic free jet rotating perpendicular to the streamwise direction

We study the flow structure of supersonic jets rotating perpendicular to the streamwise direction using RANS simulations, and we assess the performance of different turbulence model rotation corrections. The Coriolis and centrifugal terms were added to the equations of motion to perform calculations in this non-inertial (rotating) frame of reference. An explicit, cell-centred, finite-volume numerical method, coupled to a k?ε turbulence model, was used for the computations. The turbulence model rotation corrections of Howard et al. (1980), Park and Chung (1999), and Cazalbou et al. (2005) were attempted. In the absence of experimental data for jets rotating perpendicular to the streamwise direction, the rotation corrections were examined against the available measurements of a swirling jet; the comparison of the numerical and experimental data indicates that the Cazalbou et al. (2005 Cazalbou, J.B. 2005. Two-equation modeling of turbulent rotating flows. Physics of Fluids, 17(5): 1–14. [Crossref], [Web of Science ®] , [Google Scholar]) and Park and Chung (1999 Park, J.Y. and Chung, M.K. 1999. A model for the decay of rotating homogeneous turbulence. Physics of Fluids, 11(6): 1544–1549. [Crossref], [Web of Science ®] , [Google Scholar]) corrections improve the performace of the turbulence model. Simulations were then run of a supersonic jet rotating perpendicular to the stream direction at 0, 50, 100 and 150 rad/s, using no turbulence model rotation correction, and using the three rotation corrections. The results indicate that the Cazalbou et al. (2005 Cazalbou, J.B. 2005. Two-equation modeling of turbulent rotating flows. Physics of Fluids, 17(5): 1–14. [Crossref], [Web of Science ®] , [Google Scholar]) correction is more physical than the other two, as it yields results that are qualitatively consistent with the known effects of rotation: that turbulence is enhanced and suppressed on the concave and convex sides of a rotating jet centreline, respectively, and that the effect of rotation saturates as the rotation rate increases. The findings are in qualitative agreement with the available literature.     

High-order Large Eddy Simulations of Confined Rotor-Stator Flows

In many engineering and industrial applications, the investigation of rotating turbulent flow is of great interest. In rotor-stator cavities, the centrifugal and Coriolis forces have a strong influence on the turbulence by producing a secondary flow in the meridian plane composed of two thin boundary layers along the disks separated by a non-viscous geostrophic core. Most numerical simulations have been performed using RANS and URANS modelling, and very few investigations have been performed using LES. This paper reports on quantitative comparisons of two high-order LES methods to predict a turbulent rotor-stator flow at the rotational Reynolds number Re(=?Ωb ²/ν)?=4 × 10⁵. The classical dynamic Smagorinsky model for the subgrid-scale stress (Germano et al., Phys Fluids A 3(7):1760–1765, 1991) is compared to a spectral vanishing viscosity technique (Séverac & Serre, J Comp Phys 226(2):1234–1255, 2007). Numerical results include both instantaneous data and post-processed statistics. The results show that both LES methods are able to accurately describe the unsteady flow structures and to satisfactorily predict mean velocities as well as Reynolds stress tensor components. A slight advantage is given to the spectral SVV approach in terms of accuracy and CPU cost. The strong improvements obtained in the present results with respect to RANS results confirm that LES is the appropriate level of modelling for flows in which fully turbulent and transition regimes are involved.     

A Fast Algorithm for Invasion Percolation

We present a computationally fast Invasion Percolation (IP) algorithm. IP is a numerical approach for generating realistic fluid distributions for quasi-static (i.e., slow) immiscible fluid invasion in porous media. The algorithm proposed here uses a binary-tree data structure to identify the site (pore) connected to the invasion cluster that is the next to be invaded. Gravity is included. Trapping is not explicitly treated in the numerical examples but can be added, for example, using a Hoshen–Kopelman algorithm. Computation time to percolation for a 3D system having $N$ total sites and $M$ invaded sites at percolation goes as $O(M \log M)$ for the proposed binary-tree algorithm and as $O(M N)$ for a standard implementation of IP that searches through all of the uninvaded sites at each step. The relation between $M$ and $N$ is $M = N^{D/E}$ , where $D$ is the fractal dimension of an infinite cluster and $E$ is Euclidean space dimension. In numerical practice, on finite-sized cubic lattices with invasion structures influenced by the injection boundary and boundary conditions lateral to the flow direction, we observe the scaling $M = N^{0.852}$ in 3D (valid through the second decimal place) instead of $M= N^{0.843}$ based on the infinite cluster fractal dimension $D=2.53$ .