Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 3, pp. 643-657, Warsaw 2007 LES OF TURBULENT CHANNEL FLOW AND HEAVY PARTICLE DISPERSION Jacek Pozorski Tomasz Wacławczyk Mirosław Łuniewski Institute of Fluid Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland e-mail: jp@imp.gda.pl; twacl@imp.gda.pl; mirmir@imp.gda.pl TheEulerian-Lagrangianapproachhasbeen applied to two-phase turbu- lent flows with dispersed heavy particles. First, numerical computations of the continuous phase (fluid) have been performed using the Direct Numerical Simulations (DNS) and Large-Eddy Simulations (LES) for the case of a fully-developed channel flow. Parallelisation efficiency of flow and particle solvers has been estimated. Residual turbulent kinetic energy has been found from filtered DNS (also called a priori LES); a model for this quantityhas beenassessed.Then, heavyparticles havebe- en tracked in the LES velocity field. Statistics of particle motion (mean velocity, intensity of velocity fluctuations) and the profile of cross-stream number density aswell as preliminary results for particlewall deposition have been obtained. Key words: turbulent flow, large eddy simulation, particle dispersion 1. Introduction Two-phase polydispersedflows (with solid particles, droplets or bubbles) com- monly occur in various industrial applications, including power engineering (coal/spray combustion, wet steam flow), chemical and process engineering (catalytic reactions in two-phase flow, particle separators, etc.). These flows are often physically complex because of flow turbulence, particle polydispersi- ty, phase changes and/or chemical reactions (cf. Moin and Apte, 2006). Con- sequently, an improved modelling of polydispersed two-phase flows is largely perceived as an open issue. Although major progress has been achieved sin- ce the advent of Computational Fluid Dynamics (CFD), there is still a need formodel developments, both fundamental and for practical applications. Ge- nerally speaking, three main issues in CFD can be identified as: numerics 644 J. Pozorski et al. (including appropriate discretisation of governing equations, effective solution algorithms, parallelisation of numerical computations), geometrical complexi- ty (including grid generation, i.e. suitable discretisation of the flow domain, often differing for external and internal flows), and physical complexity of the flow phenomena. Out of the three, the numerical and geometrical issues have largely been overcome to date; however, soundmodelling of physical complexity remains to beapacing item in theprogress ofCFD.Quite often, the turbulencemodelling is perceived as the central issue, since a physically-sound description of the underlying flow field is a sine qua non condition for a correct and reliable computation of two-phase dispersed flows, including heat andmass transfer. As far as the turbulence modelling is concerned, various levels of descrip- tion have been developed, starting from the statistical resolution of the flow (RANS, PDF), through computation of larger flow structures (LES), up to the full resolution of three-dimensional, unsteady flow field (DNS). For seve- ral reasons, nowadays, the LES stands out as a promising approach to deal with turbulent polydisperse flows. One of the reasons is that particle motion is often dominated by energetic, large-scale flow structures that are (or should be) adequately resolved in this approach and the effect of Sub-Grid Scales (SGS) remains limited (cf. Armenio et al., 1999). Yet, the LES still needs fur- ther developments, specially forwall-bounded turbulence, towork out a viable compromise between its computational cost and accuracy. Improved models need to be put forward for fluid-particle interactions, including the subgrid- scale particle dispersion (cf. Pozorski et al., 2004; Kuerten andVreman, 2005; Shotorban andMashayek, 2006), the effects of preferential particle concentra- tion (Squires andEaton, 1991), possibly also the effect of particles on the fluid SGS models, etc. Another reason for the increasing popularity of LES is the rapid growth of available computational power that allows for simulations of more andmore complex flow cases. Although the LES is governed by deterministic equations, it still conserves a flavour of the statistical approach, since the resulting flow realisations can be compared to experimental orDNS results via the ensemble- or phase-averaged statistics, rather than via single realisations. Yet, contrary to the statistical turbulencemodelling, LES allows one to exactly resolve a spectrum of length scales in the flow (down to themesh size) and their interactions. This picture getsmore complex for two-phasedispersedflows, since thepresenceof particles creates new scales of time and length due to inertia, gravity effects, collisions, interaction with walls, evaporation, chemistry, etc. The two above facts, i.e. the practical importance of dispersed turbulent flows on the one hand and modelling difficulties involved in their sound de- scription on the other hand, provide themainmotivation for the presentwork. We report here the LES computation with Lagrangian discrete particle trac- LES of turbulent channel flow... 645 king of turbulent, particle-laden channel flow. The paper constitutes the first step towards amore comprehensive study of such flows. 2. LES of turbulent flow 2.1. Governing equations and their closure The assumed starting point are the Navier-Stokes (N-S) equations for the incompressible flow; after spatial smoothing (with a given filter function G of a length-scale ∆), they become the filtered equations for the resolved scales, or large eddies (cf. Pope, 2000) ∂Ui ∂t +Uj ∂Ui ∂xj =− 1 ρ ∂p ∂xi +ν∇2Ui− ∂τij ∂xj (2.1) According to the idea of LES, flow variables are decomposed into resolved (large-scale) and residual (subgrid-scale) parts: U = U +u′, p = p+ p′, etc., where the symbol (·) stands for spatial filtering (smoothing), being a convolution: U = G∗U. In the filtered field, the flow scales smaller than a certain cut-off length (related to ∆) are eliminated. In Eq. (2.1), τij is the residual (subgrid-scale) stress tensor, τij =UiUj−UiUj; the divergence of τij represents the effect of small-scale velocity on the resolved flow. As seen from Eq. (2.1), the large-eddy fluid dynamics is closed once a su- itablemodel for the residual stress tensor is provided.We note that additional models, e.g. for the residual kinetic energy, kr = u ′ iu ′ i/2, may also be needed for some purposes. For the present computation, the dynamic SGS model of Germano and Lilly (cf. Pope, 2000) has been applied. The deviatoric part of the residual stress tensor, τdij = τij − (2/3)krδij, is modelled as τdij =−2νrSij (2.2) where Sij is the resolved strain rate (the symmetric part of the velocity gra- dient tensor). The residual, or sub-grid scale, viscosity νr is determined as νr = CG∆ 2 |S| where |S| = √ 2SijSij. The dynamic procedure is applied for theGermanomodel coefficient CG, assuming a certain scale-similarity around the cut-off length CG = 1 2 〈LijMij〉av 〈MklMkl〉av (2.3) where Lij and Mij are defined as Lij = ÛiUj − ÛiÛj Mij =∆ 2 |̂S|Sij − ∆̂ 2 |Ŝ|Ŝij (2.4) 646 J. Pozorski et al. In Eq. (2.3), 〈·〉av stands for a subsequent averaging (e.g., over flow homoge- neity directions, whenever possible) to yield a smoother CG field. Moreover, (·), (̂·) in Eq. (2.4) stand for spatial filtering with two different filter sizes: ∆ and ∆̂, respectively. We note that for the computation of CG, the opera- tional formula for the effective filter length ∆̂ is needed. It is readily found for a spectral filter: ∆̂ = ∆̂, and for a Gaussian filter: ∆̂ = √ ∆̂2+∆ 2 (cf. Pope, 2000). Yet, a lingering question remains about the working meaning of the effective length for the top-hat filter (often used for implicit filtering, also here). The point is that the double filtering in this case is not equivalent to a single filtering with ∆̂ and the same, top-hat filter shape. To end this point, let us recall two often-used concepts. If the DNS of the N-S eqs. is run with filtering of the solution at each time step (U = G∗U, but U is not used to compute the flow dynamics), such a procedure is called a priori LES. It is often used to validate the LES closures (also for dispersed two-phase flows). As opposed to this, a true LES computation, where Eqs. (2.1) are solved, is called a posteriori LES.Both a priori and a posteriori LES computations are reported in the following. 2.2. Estimation of residual kinetic energy The issue of the residual or subgrid scale (SGS) velocity field, its influence on the statistics of particle motion and their preferential concentration has received only limited attention in the literature. Amodel for the SGS partic- le dispersion seems to be crucial for coarser LES where a considerable part of the turbulent kinetic energy (say, 20% or more) is not explicitly resolved. Such a model has to be validated first, possibly through a priori LES tests. Wall-bounded flows are of particular importance because of the LES resolu- tion difficulties, the agglomeration of particles in the near-wall zone and their separation on the walls. Although the residual turbulent kinetic energy is not needed for LES flow computation (save for SGS stress models based on the transport equation for kr), estimation of this quantity may be necessary to model some small-scale phenomena. In particular, the residual kinetic energy is needed in SGSparticle dispersionmodels. In the present work, we estimate the residual kinetic energy kr through a dynamic procedure. This is certainly appealing since we have gone for the dynamic procedure to close the SGS stress term, too. The starting point of this approach is the assumption that (cf. Shotorban andMashayek, 2006) kr =CI∆ 2 |S|2 (2.5) LES of turbulent channel flow... 647 where the proportionality parameter CI is not really a constant but varies in space and time like the Germano coefficient CG. Analogously to CG, the dynamic procedure with double filtering is applied to solve for CI CI = 1 2 〈ÛkUk− ÛkÛk〉av 〈∆ 2 |̂S|2− ∆̂ 2 |Ŝ|2〉av (2.6) 2.3. Flow and particle solvers We consider turbulent particulate flows in the Eulerian-Lagrangian approach where the three-dimensional, unsteady simulations of the carrier phase: DNS, a priori LES (filteredDNS) and true LES have been performed. The LES has been coupled with the particle tracking. As the flow solver, we have used an open-source academic code FASTEST3D, currently being developed by the group of Prof. Michael Schäfer (TU Darm- stadt,Germany). It is afinitevolume, second-order accuracy, block-structured, full-multigrid solver for incompressible flows, with the SIMPLE pressure- velocity coupling scheme, suitable for parallel computations of general CFD applications. As the particle solver, we have further developed our in-house PTSOLVmo- dule for the Lagrangian particle tracking. It provides the necessary interpo- lation of flow variables at the particle locations, solves a system of ordinary differential equations for particle positions and velocities, and computes sta- tistics. In the present work, the dispersed phase has been assumed dilute; consequently, the one-waymomentumcoupling (fluid-to-particles) is adequate and no inter-particle collisions are considered. 2.4. Parallelisation issues The LES computations using FASTEST3D for fluid flow and our PTSOLVmodule for particle motion have been parallelised using theMPI library. The parallel code has been prepared and tested locally; only afterwards the computational task has been submitted to amainframe server for a production run. The flow domain has been divided into a number of equal-size parallel blocks to be run on separate processors. For the particle solver, each proces- sor has been assigned an approximately equal number of particles, indepen- dently of their spatial locations. This choice seems suitable for various flow cases, including those with highly non-uniform particle distributions, such as spatially-developing particle-laden jets, or wall-bounded flows where particles can accumulate in the near-wall regions due to turbophoresis.Themain draw- back of the chosen approach is the need tomake the whole fluid velocity field available for each processor at each flow time step to perform the interpolation 648 J. Pozorski et al. of flow velocity to current particle locations. This may become problematic, or even infeasible, for cases where a huge number of computational cells are used for fluid. An alternative choice would be to assign to each of the pro- cessors the particles currently located in the flow subdomain served by that processor. For flow cases where there is nomajor risk of load unbalancing due to a non-uniform particle distribution in space and for bigger computational problems, that choice would bemore efficient. Yet, for this alternative option there would be a need to exchange some particles between processors at each time step. Fig. 1. Total elapsed time (arbitrary units) vs. the number of processors; particle solver (•), flow solver on a fine grid (◦), on a coarser grid ( ); dashed line: ideal scalability (slope−1 on log-log plot) In order to find an optimum number of processors to be used in simula- tions, we have performed scalability tests from parallel computations. First, we ran the flow solver alone (no particles) on 44×44×32 and 128×128×128 meshes for several time steps using a number of processors varying from 1 up to 32. Then, we ran the particle solver only (in a frozen fluid velocity field) using 3.2M Lagrangian fluid elements. The results of scalability tests are qu- ite instructive; they are shown in Fig.1. First observation is that the switch from a single-processor configuration to a two-processor parallel computation is only marginally effective; this is explained by internal computer architec- ture (fast cache memory available for single-processor runs, no need for data exchange). Then, the parallel runs become quite efficient when passing to 4 and 8 processors. However, a shift from 8 to 16 processors for fluid compu- tation is rather useless, probably because of the increase of ”area-to-volume” ratio of each flow sub-domain. As expected, we note a better scalability for a finer fluidmesh (more grid points) than for a coarser one. On the other hand, the particle solver alone shows a good scalability for up to 32 processors; only beyond, the inter-processor data exchange (necessary to compute the fluid ve- locity seen by particles) makes a further increase of the number of computing LES of turbulent channel flow... 649 nodes problematic. In conclusion, for typical mesh sizes (643 and 1283) and a typical number of particles (around 106) of our problem, we have decided to stay with 8-processor runs. 3. Computation results for single-phase turbulence The computations of the channel flowcase have beenperformed for Reτ =150 (the Reynolds number based on the friction velocity uτ and the channel half- width h) on the 1283 grid for the DNS, and on two coarser grids for the LES. The size of the flow domain in the streamwise (x), wall-normal (y) and spanwise (z) directionswas 4πh×2h×2πh. The flowwas assumed periodic in the streamwise and spanwise directions, and driven by the constant pressure gradient. Parallel computations have been run on 8 processors. To validate the LES and to gather necessary data for a priori tests, we ran the DNS of the channel flow first. The mesh size (in wall units) was ∆x+ = 14.7 and ∆z+ = 7.35. In the wall-normal direction, ∆y+ varied from 0.05 at thewall up to 6.75 at the channel centerline. The simulation time step (in viscous units) was ∆t+ = 3.6. The necessary time until convergence was about 150 hours (total wall clock elapsed time) or 30 FTT (flow-through times, Lx/UCL). Then, flow statistics were gathered for another 50 hours (or 10 FTT). All obtained velocity statistics compare very well with results from another DNS (Marchioli and Soldati, 2006, priv. comm.), performed using a different approach (spectral resolution in the x and z directions) and solver. Next, we ran a priori LES, i.e. we filtered theDNS field at each time step, to have a look at results of the dynamic procedure. In particular, we obtained an estimate of the residual turbulent energy kr, Eq. (2.5). This quantity will be a crucial ingredient of a SGS particle dispersion model (currently tested by the authors), but is also interesting as such. The cross-stream profile of CI is shown in Fig.2a. Its behaviour is qualitatively similar to that of CG (not shown): it decreases to zero at the walls; yet, a single maximum is noticed in the core of the flow. Using thevalues of CI computed fromthea prioriLESof channel flow, the residual turbulent kinetic energy is estimated from Eq. (2.5). It qualitatively agrees with the profile of kr resulting from the filtering of DNS results; both are plotted in Fig.2b. The results have been obtained from a single snapshot (iteration) of theDNS velocity field with averaging over x and z. It is readily noticed that the residual turbulent kinetic energy kr predicted by Eq. (2.5) is over-estimated since it attains up to 10% of the total turbulent energy k, and the DNS filtering for this case (done on 1283 to the 643 grid) results in up to 3% of k removed. 650 J. Pozorski et al. Fig. 2.A priori LES of channel flow; (a) profile of the CI parameter; two different effective filter width coefficients: Gaussian (1.25 – solid lines) and spectral (1.0 – dashed lines), (b) turbulent kinetic energy: total (k, 1283 DNS – solid line), resolved (kf , 64 3 filtered DNS – dashed line), residual (kr ≡ ksg, two CI coefficients – lower lines) Following the DNS computation, we performed the LES of channel flow. There, two coarser mesheswere chosen, non-uniform in thewall-normal direc- tion: amesh of 643 with the cell size varying from ∆y+ =0.22 to ∆y+ =13.5 and a finer mesh with 86 cells in the wall-normal direction. The wall-clock ti- me until convergence was about 48 hours, or 20 FTT; the flow statistics were gathered for another 10 FTT (or 24 hours of the elapsed time). Themean velocity and the turbulent shear stress are shown inFig.3 in the near-wall scaling. The results for LES on 643 mesh are not accurate enough whereas those on afinermesh compare quitewellwith theDNS referencedata. As reported in other LES studies, to obtain good velocity statistics, the flow resolution in the wall-normal direction has to be close to that of the DNS. Fig. 3. (a) The mean velocity profile, (b) the turbulent shear stress; LES – lines; DNS – symbols The resolved and residual turbulent kinetic energy in the LES of the chan- nel flow on the 643 mesh are shown in Fig.4. As compared with Fig.2b, both LES of turbulent channel flow... 651 statistics are consistently over-estimated, probably due to insufficient LES re- solution in the near-wall region. Fig. 4. Profiles of turbulent kinetic energy in LES computation on the 643 mesh (single iteration, averaging over x,z); solid line – resolved energy k+, lower dashed lines – residual energy (estimated k+r , two different CI coefficients) 4. Heavy particles in turbulent flow 4.1. Particle equation of motion In case of heavy particles (ρp ≫ ρf), the governing equations are sim- plified and (in the absence of gravity) contain only the drag force, based on the particle velocity Up and the fluid velocity along the particle trajectory U ∗ f =Uf(xp, t) dxp dt =Up (4.1) dUp dt = 3 4 ρf ρp CD dp |U∗f −Up|(U ∗ f −Up) The drag coefficient CD is taken from a well-established correlation CD = 24 Rep (1+0.15Re0.687p ) (4.2) where Rep = dp|U ∗ f −Up|/ν is the particle Reynolds number (based on the particle diameter dp, the relative particle velocity, and the kinematic viscosity of the carrier fluid ν). Let us also define the particle relaxation time τp = ρp ρf d2p 18ν (4.3) 652 J. Pozorski et al. In the limit of small Rep, the drag term in Eq. (4.1)2 becomes (U ∗ f−Up)/τp. In a general case (e.g. for the particle density comparable to that of fluid), the equation of motion should also contain other terms (like the added-mass force, the lift force or the history term) that may be significant. In the statistical description of turbulence only the mean fluid velocity is available, and a reconstruction of the fluctuating field is needed. Propo- sals for the statistical approach were put forward by Pozorski and Minier (1999). They consist in formulating the joint PDF of fluid and particle velo- city, f(U∗f,Up;x, t) and then stochastic modelling of the fluid velocity ”se- en” by particles: dU∗f = Adt+BdW . On the other hand, the simplest tre- atment of the particulate phase in LES is to use in Eq. (4.1)2 the filtered (large-scale) field instead of the instantaneous velocity field ”seen” by par- ticles: U∗f → U ∗ f = Uf(xp, t) with no account for the residual scales. This approach has been used here. In practical implementation, one needs to inter- polate the fluid velocity (known on a discrete mesh) at particle locations. In the present work, a tri-linear interpolation has been used. This was deemed sufficientbecause of the spatial integration scheme in thefluid solver.However, a sensitivity analysis with a higher-order interpolation scheme is foreseen. The simplest numerical integration scheme for the particle equations of motion is the first-order explicit scheme, for small Rep written as U n+1 p =U n p +(U ∗n f −U n p) ∆t τp over a time step ∆t = tn+1 − tn; yet, the scheme becomes unstable for ∆t > O(τp). The same remains true also for higher-order Runge-Kutta me- thods based on this predictor. Hence, an acceptable time step (specially for small particles) can become smaller than the time step used in the LES solver for fluid. In this case, an unconditionally stable, first-order scheme may be preferred. Assuming that U∗f is frozen over ∆t, the exact analytical solution yields an ”exponential” scheme U n+1 p =U ∗n f +(U n p −U ∗n f )exp (−∆t τp ) Higher-order schemes of this kind were proposed by Barton (1996) and Pe- irano et al. (2006). We have checked some statistics of particle motion and found them independent of the time step used; however, a second-order in- tegration scheme would allow us to increase the time step, in particular for small particles. 4.2. Particle statistics in channel flow A convenient measure of particle inertia in the flow is the Stokes number, here defined as St = τ+p = τp/(νf/u 2 τ). In a statistically-steady LES of fluid LES of turbulent channel flow... 653 velocity, several particle groups (”swarms”) of varying inertia, characterised by the Stokes numbers of St = 2, 10, 50, 250, have been tracked for time t+ =1400 (non-dimensionalised with viscous units) with perfect-reboundwall boundary conditions. In computations, particle trajectories are approximated by line segments, and implementation of this boundary condition for the first- order integration scheme of particle equations is straightforward. Themean velocity profiles of particles and fluid across the channel and the intensities of streamwisevelocity fluctuationsare shown inFig.5.Note that the mean particle velocity tends to zero at the walls for all particle classes except the largest ones. At the channel center, the particlemean velocity tends to lag behind that of the fluid. It also transpires fromFig.5a that themean velocity of smallest particles computed here (St= 2) is quite close to that of the fluid; however, the differences become more evident for the second-order moments, as seen in Fig.5b. This may be due to a net wallward motion and inertia, specially for larger particles. Fig. 5. Profiles of mean velocity and intensity of streamwise velocity fluctuations: fluid and particles Themean concentration profiles, i.e. theparticle numberdensity across the channel at t+ =1400 are shown in Fig.6. The phenomenon of turbophoresis (the tendency of heavy particles to move off the highest turbulence intensity regions) is readily observed; the non-uniformity of particle number density across the channel increases with particle inertia. We have also checked that the concentration profile for smaller particles (St = 2) has already achieved its steady state at the simulation time consi- dered (t+ = 1400); the concentration profile for larger particles (St = 10) is also close to stationary and shows stronger gradients, whereas for the largest particles (St = 250) the non-uniformity of number density profile across the 654 J. Pozorski et al. Fig. 6. Particle number density across the channel. Initial distribution (statistically uniform) – dotted line; distribution at t+ =1400 for particles of St=2 – solid line, St= 10 – dot-dashed line, and St= 250 – dashed line channel further develops and ends upwith a considerable level of particle con- centrations close to the walls, ultimately undermining the assumption on the dispersed phase being dilute. These trends are noticed in the following Fig.7. There, the time record of thenon-uniformityof particle concentration σ across the channel is shown for various particle classes. The value of σ is measured as the r.m.s. departure of the particle number density from its initial value c0 (that of the spatially-uniform distribution) and normalised with c0. At the simulation time considered, the number-density profile of smallest particles is already statistically-steady, whereas for larger particles the non-uniformity still develops. Fig. 7. Time records of the non-uniformity of number density profiles for particle Stokes number: 2 – solid line, 10 – dotted line, 50 – dashed line, and 250 – dot-dashed line 4.3. Particle separation on walls In this Section, we consider particle motion in a turbulent flow with the ab- sorbing boundary condition at the wall. Themodelling of particle separation, LES of turbulent channel flow... 655 or deposition, on solid walls presents a major difficulty for the statistical, one-point approaches, because of the crucial role of the fluid fluctuating velo- city in the process, at least for moderate-inertia particles. In their case, the mechanisms of deposition involve complex interactions with near-wall vorti- cal structures. Consequently, it should in principle be easier to achieve better results for particle deposition in a well-resolved LES. The mass flux Jw of particles separating on the channel walls, divided by the particle bulk density ρdp (computed out of the current number of particles ntot in the flow domain), is quantified as the particle separation velocity (cf. Uijttewaal andOliemans, 1996). The non-dimensional deposition velocity is written as V+ dep = Jw ρdpuτ = ndep ntot h 2uτ∆t where ndep is the number of particles that deposit on the walls in the time interval ∆t. Fig. 8. Separation velocity for heavy particles: LES results (2) vs. experimental data of Liu and Agarwal (1974) at Re=10000 (•) and Re=50000 (N) To estimate the particle deposition, the LES computation of the channel flow at Reτ =150 (Re≈ 3000) with particle tracking (no SGS dispersionmo- del) has beenperformedwith the perfectly absorbingwall boundarycondition. Results for particle mass flux separating on the wall are presented in Fig.8 as the non-dimensional deposition velocity V+ dep for a selection of particle Stokes numbers and compared to the experimental data for turbulent pipe flow (Liu and Agarwal, 1974). Although the experiment has been performed at higher Reynolds numbers, the experimental data itself (• and N) suggests that the particle separation is rather weakly dependent on Re. The separation velocity is visibly underestimated in our LES, specially for particles of small inertia, probably due to the neglected effect of subgrid-scale flow scales and/or the lift force. An addition of the sub-grid scale particle dispersionmodel considerably improves the results (work in progress). 656 J. Pozorski et al. 5. Conclusion In the paper, we have described our DNS and LES of the turbulent channel flow followed by the LES studies of the particle-laden case where the flow so- lver has been coupledwithdiscrete particle simulations of the dispersedphase. The outcome of parallel tests has been presented. The estimation of residual turbulent kinetic energy from the dynamic procedure has been shown toge- ther with a priori LES results. Reasonable flow statistics have been obtained from the LES on a finer grid. Particle velocity statistics have been compu- ted and differences with respect to fluid ones have been shown, including the phenomenon of turbophoresis in the near-wall region. The next step will be the computation of particle wall deposition (sepa- ration from flow) with the account of subgrid-scale fluid motion on particle dynamics modelled with the use of residual kinetic energy. Moreover, it se- ems advantageous to implement a second-order numerical scheme for particle motion. At a longer term, the promising idea about the dispersed two-phase turbulent flows (including reactive flows or combustion) seems to be a hybrid approachwhere large-scale turbulent structureswill becomputed fromLES(or POD), and sub-grid processes (particle dispersion, species transport/mixing, chemical reactions) will bemodelled using the FDF (filtered density function) formulation. Acknowledgments The work has been partly supported by the European Community under the auspices of the High Performance Computing (HPC-Europa) programme. J.P. and T.W. wish to express their warmest thanks to Prof. Alfredo Soldati and Dr. Cristian Marchioli (University of Udine, Italy) for their hospitality (June 2006). Also, sincere thanks go to Dr. Giovanni Erbacci and Andrea Tarsi from CINECA (Bologna) for their interest and kind assistance in our using the computing resources (CLXparallel server). Some simulation runswere also performed at the regionalComputing Centre (TASK) at Gdańsk. The licence agreement of Professor Michael Schäfer (Dept. of NumericalMethods inMechanical Engineering,TUDarmstadt,Germany) to use and further develop the research code FASTEST3D is gratefully acknowledged. References 1. Armenio V., Piomelli U., Fiorotto V., 1999, Effect of the subgrid scales on particle motion,Phys. Fluids, 11, 3030-3042 2. Barton I.E., 1996, Exponential-Lagrangian tracking schemes applied to Sto- kes law, J. Fluids Engng., 118, 85-89 3. Kuerten J.G.M., Vreman A.W., 2005, Can turbophoresis be predicted by large-eddy simulation?Phys. Fluids, 17, 011701 LES of turbulent channel flow... 657 4. Liu B.Y.H., Agarwal J.K., 1974, Experimental observation of aerosol depo- sition in turbulent flow,Aerosol Sci., 5, 145-155 5. Moin P., Apte S.V., 2006, LES ofmultiphase reacting flows in complex com- bustors,AIAA J., 44, 698-708 6. Peirano E., Chibbaro S., Pozorski J., Minier J.P., 2006, Mean- field/PDF numerical approach for polydispersed turbulent two-phase flows, Prog. Energy Combust. Sci., 32, 315-371 7. Pope S.B., 2000,Turbulent Flows, Cambridge University Press 8. Pozorski J., Apte S., RamanV., 2004,Filtered particle tracking for disper- sed two-phase turbulent flows,Proceedings of the Summer Program, Center for Turbulence Research, Stanford University, 329-340 9. Pozorski J., Minier J.P., 1999, PDF modeling of dispersed two-phase tur- bulent flows,Phys. Rev. E, 59, 855-863 10. Shotorban B., Mashayek F., 2006, A stochastic model of particle motion in large-eddy simulation, J. Turbulence, 7, art. 18 11. Squires K.D., Eaton J.K., 1991, Preferential concentration of particles by turbulence,Phys. Fluids A, 3, 1169-1178 12. Uijttewaall W.S.J., Oliemans R.V.A., 1996, Particle dispersion and de- position in DNS and LES,Phys. Fluids, 8, 2590-2604 Zastosowanie metody dużych wirów do przepływu turbulentnego w kanale i dyspersji cząstek Streszczenie W pracy zastosowano podejście Eulera-Lagrange’a do wyznaczenia ruchu fazy dyspersyjnej (cząstki, krople) w turbulentnym przepływie dwufazowym. Obliczenia ruchu fazy ciągłej (płynu) przeprowadzono za pomocą rozwiązania pełnych równań przepływu (DNS) orazprzyużyciumetodydużychwirów(LES) dla przypadku rozwi- niętegoprzepływuwkanalepłaskim.Oszacowanoefektywnośćzrównolegleniaobliczeń numerycznych.Określono poziom energii kinetycznej skal podsiatkowych; porównano ją z wynikami uzyskanymi z modelu dynamicznego.Wyznaczono trajektorie cząstek fazy dyspersyjnej w polu prędkości LES. Określono statystyki ich ruchu: prędkość średnią, intensywność fluktuacji prędkości, profil koncentracji cząstek w poprzek ka- nału; uzyskano wstępne wyniki dla separacji cząstek na ściankach. Manuscript received January 22, 2007; accepted for print May 14, 2007