Effect of Reynolds Number on Particle Interaction and Agglomeration in Turbulent Channel Flow

The work described in this paper employs large eddy simulation and a discrete element method to study turbulent particle-laden channel flows at low concentrations (particle volume fraction 10 -4 -10 -5 ), including particle dispersion, collision and agglomeration. Conventional understanding of such flows is that particle interactions are negligible, this work however demonstrates that such interactions are common at large Stokes numbers in turbulent flow. The particle-particle interaction model is based on the Hertz-Mindlin approach with Johnson-Kendall-Roberts cohesion to allow the simulation of cohesive forces in a dry air flow. The influence of different flow Reynolds numbers, and therefore the impact of fluid turbulence, on agglomeration behaviour is investigated. The agglomeration rate is found to be strongly influenced by the flow Reynolds number, with most of the particle-particle interactions taking place at locations close to the channel walls, aided by the higher turbulence levels and concentration of particles in these regions. modelled the agglomeration of rigid, dry and electrostatically neutral particles in turbulent gas flows. This work was based on a deterministic collision model in the framework of a coupled large eddy simulation-Lagrangian particle tracking approach for the description of the dispersed particle-laden flows considered. Based on their results, it was concluded that an enhanced agglomeration model, using the closely-packed sphere model for the arising agglomerates, realistically predicted the physical behaviour of the agglomeration process within particle-laden flows. Where the fluid phase has been modelled using RANS techniques, recent works Tong et al. who investigated the behaviour of agglomerate-agglomerate collisions based on a combined RANS and discrete element method (DEM) approach, aimed at developing a better understanding of the underlying mechanisms of loose drug agglomerate aerosolisation in dry powder inhalers. The collision between two mannitol agglomerates in a T-shape pipe was simulated, and the effects of key variables such as the agglomerate strength, air flow rate and collision angle were investigated. The results showed that the collision between agglomerates has a significant effect on the aerosolisation process.

Understanding the fundamental aspects of turbulent fluid-particle flows is of relevance to processes employed in a wide range of applications, such as oil and gas flow assurance in pipelines, powder dispersion from dry powder inhalers and particle re-suspension and processing in nuclear waste ponds and silos. Despite its importance, little is known about the influence of inter-particle collisions on the particle and fluid phase characteristics in the context of particle dispersion, collision, agglomeration and deposition in such turbulent, bounded flows laden with large particle numbers. For example, a major flow assurance problem encountered in oil and gas production is 'scale' which is formed by inorganic, sparingly soluble salts from aqueous brines (Cowan andWeintritt, 1976 andCarrell, 1987). The build-up of scale on pipe walls occurs under supersaturated conditions, for instance in the mixing of incompatible fluid streams, e.g. the formation water from the bottom hole and the injected seawater. The deposited scale adheres to the surfaces of the production well tubing and on parts of water handling equipment, where it accumulates over time and decreases flow rates in reservoirs, pumps, valves and topside facilities (Nasr-El-Din et al., 2001). The performance of equipment that involves heat transfer processes (e.g. boilers and heat exchangers) is further lowered due to a decrease in heat exchange rates (Karabelas, 2002). These phenomena are also encountered downstream (e.g. in distillation plants), where the build-up of mineral deposits damages equipment parts. In order to amend or replace these parts, operations have to be halted, and usually with associated high costs (Chernozubov et al, 1966).
The literature on fully coupled particle-laden flows can be divided into two groups based on the method used in the calculation of the fluid phase; flows where the fluid has been modelled using Reynolds-averaged Navier-Stokes (RANS) techniques, and where it has been simulated using large eddy simulation (LES) or direct numerical simulation (DNS). Where the fluid phase has been simulated, recent works include that of Alletto (2014) who first generalised a particle agglomeration model to take into account oblique collisions, thereby allowing relative tangential velocities at the contact point. In this work, a priori analysis and posteriori evaluation within a downward pipe flow at low Reynolds number were conducted giving reasonable results, although no validation procedure was included. Breuer and Almohammed (2015) modelled the agglomeration of rigid, dry and electrostatically neutral particles in turbulent gas flows. This work was based on a deterministic collision model in the framework of a coupled large eddy simulation-Lagrangian particle tracking approach for the description of the dispersed particle-laden flows considered. Based on their results, it was concluded that an enhanced agglomeration model, using the closely-packed sphere model for the arising agglomerates, realistically predicted the physical behaviour of the agglomeration process within particle-laden flows. Where the fluid phase has been modelled using RANS techniques, recent works include that of Tong et al. (2016) who investigated the behaviour of agglomerateagglomerate collisions based on a combined RANS and discrete element method (DEM) approach, aimed at developing a better understanding of the underlying mechanisms of loose drug agglomerate aerosolisation in dry powder inhalers. The collision between two mannitol agglomerates in a T-shape pipe was simulated, and the effects of key variables such as the agglomerate strength, air flow rate and collision angle were investigated. The results showed that the collision between agglomerates has a significant effect on the aerosolisation process.
For a given agglomerate, therefore, there was a threshold velocity above which aerosolisation was increased, and below which it was decreased. Analyses of the air flow field and the agglomerate properties indicated that aerosolisation performance was determined by two competing factors, i.e. interparticle cohesion and the total collision energy of the agglomerates. The work noted above highlights some of the research currently ongoing regarding the prediction of inter-particle collisions and their effect on the fluid and particle phase characteristics in particle-laden flows. It is important to note, however, that some predictive approaches currently employed are rudimentary as they either lack a soft sphere particle contact model for the particle phase, or the accuracy of LES or DNS in predicting the fluid phase. In this work, large eddy simulation is coupled with the discrete element method to provide a more detailed and broadly applicable description of flows containing solid particles. The work described builds on previous findings presented in Afkhami et al. (2015) which used LES and DEM to demonstrate that particle agglomeration in turbulent channel flows can occur at relatively low particle mass loadings.

Numerical simulation approach
In this work the fluid phase is calculated using large eddy simulation which is capable of accurately predicting complex fluid flow phenomena, with flow solutions provided by this method coupled to a discrete element method to predict particle motion and interaction. The bidirectional interaction of the particles with the flow and between the particles is considered, thus creating a four-way coupled methodology for simulating turbulent particle-laden flow. In the following the basic features of the methods used are described. For further information the reader is referred to a previous publication (Afkhami et al., 2015).

Large eddy simulation
In LES, only the largest and most energetic scales of motion are directly computed, whilst the small scales are modelled (Smagorinsky, 1963). Any function is decomposed using a localised filter function such that filtered values only retain the variability of the original function over length scales comparable with, or larger than, that of the filter width. The LES employed used a top-hat filter as this fits naturally into a finite-volume formulation. This is then used to decompose the continuity and momentum equations for an incompressible Newtonian fluid with constant properties into resolved and unresolved fields, bringing about terms which represent the effect of the sub-grid scale (SGS) motions on the resolved scale motions. The SGS stress model employed in this work was the dynamic model of Germano et al. (1992), applied using the approximate localisation procedure of Piomelli and Liu (1995). In this model the Smagorinsky constant is computed as a function of time and space based on the information provided by the resolved scales of motion. Computations were performed using the commercial CFD code ANSYS Fluent.
The code implements an implicit finite-volume incompressible flow solver using a co-located variable storage arrangement. Because of this arrangement, a procedure similar to that outlined by Rhie and Chow (1983) is used to prevent checker-boarding of the pressure field. In the solver, diffusion terms are discretised using a central differencing scheme. The governing equations are solved in a sequential (segregated) manner. The discretised algebraic equations are solved using a point-wise Gauss-Seidel iterative algorithm. An algebraic multi-grid method is employed to accelerate solution convergence. For temporal discretisation, the segregated solver uses a three-level, second-order scheme. Time advancement is performed via an implicit method for all transport terms, and the overall procedure is second-order accurate in both space and time. The code is parallel and uses a high performance and production quality implementation of the Message-Passing Interface standard for HP servers and workstations for Microsoft Windows operating systems (HP MPI). For further information on the mathematical model employed, and the numerical algorithm and its application, the reader is referred to Afkhami et al. (2015) and the ANSYS Fluent 15.0 Theory Guide (2013).

Discrete element method
A Lagrangian approach was used to model particle motion, with the particles tracked along their trajectories through the unsteady, non-uniform flow field. The particles can have two types of motion: translational and rotational. Their paths are computed based on Newton's second law for the translational and rotational accelerations. This is achieved by integrating the accelerations over a time-step, with particle velocities and positions updated. All particles were assumed to be soft spheres with particles much heavier than the fluid ( n / f >> 1). The most significant forces in such systems are the drag and buoyancy forces, although buoyancy was neglected due to the effect of gravity not being considered. The shear-induced Saffman lift force was taken into account as it assumes non-trivial magnitudes in the viscous sub-layer. A modified spherical, free-stream drag for calculation of the force on the particles was employed.
All fluid parameters were taken from the fluid finite-volume cell which contained the centre of the DEM particle. This treatment is therefore only valid for particles of the same size as, or smaller than, a fluid finite-volume cell (the maximum packing fraction used in this work was 0.95), or where the change in fluid parameters (velocity, density, viscosity, etc.) over the extent of a particle remain roughly constant. The particle-laden flow was assumed to be dilute (particle volume fraction 10 -4 -10 -5 ), and the method incorporated full coupling between the phases, i.e. interactions between the particles were considered, and the flow and particles were two-way coupled. Particle-wall collisions were assumed to be inelastic, with the coefficient of restitution set to 0.5.
Particle-particle interactions were modelled using the discrete element method incorporating the no-slip contact model of Herz-Mindlin with Johnson-Kendall-Roberts (JKR) cohesion to allow the simulation of the van der Waals forces which influence particle behaviour (Johnson et al., 1971). It has been shown that elastic adhesive models such as the JKR or van der Waals based models can recreate similar initial loose packings and, to a certain extent, capture the initial compression of the material. This is not the case for many of the other existing contact elastic models such as the DMT (Derjaguin, Muller and Toporov, 1994), EDEM linear cohesion model or capillary force models, which, may fail to capture the correct stress historydependent behavior (Morrissey and Thakur, 2014).
The approach only considered the attractive forces within the contact area, i.e. the attractive inter-particle forces are of infinite short range. JKR builds on the conventional Hertz model by incorporating an energy balance to extend it to cover two elastic-adhesive spheres. The contact area predicted by the JKR model is larger than that given by the Hertz model; this creates an outer annulus in the contact area which experiences tensile stress. This annulus surrounds an inner circular region over which a Hertzian compressive distribution acts (Thornton and Yin, 1991). When two spheres come into contact, the normal force between them immediately drops to a certain value (8/9 fc, where fc is the pull-off force (Thornton and Ning, 1998)) due to van der Waals attractive forces. The velocity of the spheres gradually reduces and some of the initial kinetic energy is radiated into the substrate as elastic waves. The loading stage is complete when the contact force reaches a maximum value and particle velocity drops to zero.
In the recovery stage, the stored elastic energy is released and converted into kinetic energy causing the spheres to move in opposite directions. All the work done during the loading stage has been recovered when the contact overlap becomes zero. At this stage, however, the spheres remain adhered to each other and further work (known as work of cohesion) is required to separate the surfaces. The contact breaks at a negative overlap, f, for a contact force 5/9 fc (Ning, 1995). The particle surface attractive force was altered by specifying the interface energy, , with the amount of interface energy influencing the cohesion of the material. The speed of disturbance waves was approximated by Rayleigh surface wave propagation based on the physical properties of the discrete medium. The time-step used must then be sufficiently less than the Rayleigh time-step in order to ensure realistic force transmission rates in the assembly and to prevent numerical instability (Ning and Ghadiri, 2006). In this work the DEM Solutions commercial software EDEM was used for the DEM simulations which is fully coupled with the ANSYS Fluent CFD software. For further details the reader is referred to Afkhami et al. (2015).

Results and discussion
The treatment of discrete particle motion amongst the large scale turbulence predicted by the LES leads to questions concerning the performance of the numerical methods used in the fluid flow solver, and the accuracy of the interpolation scheme adopted to calculate the instantaneous fluid velocity at the particle location. All particles were individually tracked and their velocity calculated at the particle centre using interpolated fluid velocity values. In this approach, parameters such as the finite-volume grid resolution and the time-step used in solving the governing flow balance equations are important.
Moreover, the level of fluid turbulence in wall-bounded turbulent flows has an influence on particle-particle interactions and subsequent particle dispersion and deposition phenomena.
Modelling these physical processes is extremely important, particularly when using techniques that lack the detail of DNS. The difficulty is associated with the complicated interactions between non-homogenous turbulent structures in the wall-normal direction and the inertia of the particles. Small scale turbulent structures are not solved for in any way by the LES approach adopted, eliminating their effect on tracked particles which may, nevertheless, be influenced by such motions. It is therefore necessary to assess the level of mesh resolution used and the ability of the SGS model employed to predict accurately the selective response of different inertial particles. That said, for the simulations considered below, the particles were large with a relaxation time greater than that of the smallest fluid time scales, therefore the influence of the unresolved fluctuating velocities in the LES on particle motion may be deemed unimportant (Pozorski and Apte, 2009).
The aim of the work described is to predict the behaviour of particles in a turbulent two-phase channel flow, with the potential to lead to physical insights in to how particles disperse and agglomerate in such flows. The ability of the LES approach within the FLUENT platform has been assessed in Afkhami et al. (2015), and was found to be capable of providing accurate predictions of single-phase low Reynolds number, turbulent channel flows. It is important to note, however, that such computational fluid dynamic codes, even though frequently exploited in predicting high Reynolds number flows in complex geometries, require further validation of their ability to reliably predict multiphase flows. At present, however, this is difficult due to the lack of appropriate experimental data, and in the absence of DNS-based predictions over a wide enough range of Reynolds numbers and particle sizes, and the fact that they generally do not accommodate flow-particle interaction at the level of four-way coupling with agglomeration.

Flow configuration and initial conditions
The flow into which particles were introduced was a turbulent channel flow of air; Figure 1, gives a schematic diagram of the channel geometry and co-ordinate system. The flow is described using a three-dimensional Cartesian co-ordinate system (x, y and z) representing the streamwise, spanwise and wall-normal directions, respectively. The velocities in the x, y and z-directions are Ux, Uy and Uz, respectively. The boundary conditions for the momentum equations were set to no-slip at the channel walls, and the instantaneous flow field was considered to be periodic along the streamwise and spanwise directions, with a constant mass flux through the channel in the streamwise direction maintained by a dynamically adjusted pressure gradient used to drive the flow. The rectangular channel considered was of dimensions Figure 1). The length of the channel in the streamwise direction was sufficiently long to capture the streamwise-elongated, near-wall turbulent structures that exist in wall-bounded shear flows; such structures are usually shorter than ~ 1000 wall units (Robinson, 1991). Some variables reported in this work are in dimensionless form, represented by the superscript (+), and expressed in wall units, with the latter obtained by combining variables with u , and f .

Single-phase flow
For a smooth wall, suitable scaling parameters for the inner layer include the viscosity, v, and the friction velocity,  u . Inner layer scaling then requires that the relationship given below holds for the mean streamwise fluid velocity, x U : where Ux + is the non-dimensional mean streamwise fluid velocity, f is a universal function (independent of Reynolds number) and  z is the dimensionless distance from the wall. The mean velocity profile in the inner and outer layers of a turbulent channel flow at high Reynolds number may be represented using the expressions given by von Karman (1930): The region near the wall can be divided into two sections, the viscous sub-layer between 5 0    z and the buffer layer between In this section some of the most relevant statistics for the fluid phase are presented and discussed to benchmark the performance of the LES approach. It is important to mention here that all fluid velocity statistics presented in this paper refer to a fully developed flow. This was deemed to have been achieved when the first-and second-order moments (specifically, the mean streamwise velocity, the root mean square (rms) of the fluctuating velocity components in all three co-ordinate directions and the shear stress) remained constant with time. To achieve smooth profiles the fluid statistics were both spatial-and time-averaged over hundreds of thousands of time steps. It is always beneficial to compare numerical predictions against experimental data, however, in its absence the present LES results are compared with DNS predictions. The results generated by the LES for the fluid phase were therefore compared with DNS predictions for shear Reynolds flows of Re = 150, 300 and 590 obtained by Marchioli et al. (2008), Marchioli and Soldati (2007), and by Moser et al. (1999), respectively, as well as from other authors, as noted below. The LES results for Re = 300 have been published in Afkhami et al. (2015) where they showed satisfactory agreement with DNS, and for that reason will not be repeated in this paper. providing reasonable agreement overall. The LES clearly predicts the viscous sub-layer region to a high degree of accuracy and quantitatively tends towards (2) as this region is approached.
It is seen that the LES slightly over predicts the DNS in this region, although the log scale used emphasises any discrepancies close to the wall and therefore highlights those differences. The logarithmic law given by (3) is shown for the region z + > 30, based on the values suggested by Kim et al. (1987), with the LES results seen to over predict this analytical profile and the DNS results in this region, although the various approaches come in line at the centre of the channel.
In this region of the channel, the flow characteristics are dominated by large energetic scales of motion and, given that these scales are directly computed by the LES, the predicted profiles should match those of the DNS, with the differences observed likely due to the lack of resolution in the LES. Overall, however, the streamwise mean velocity generated by the LES is in acceptable agreement with the DNS. For the normal stresses in Figure 3, the LES results are in good agreement with the DNS for the U x,rms + component, with the magnitude and position of the peak and centre-line values of this profile predicted reasonably well. The U y,rms + and U z,rms + profiles follow the trend of the DNS, although qualitative and quantitative differences are observed in some regions. For U y,rms + , an over prediction by the LES increases through the buffer layer where it reaches a maximum before decreasing towards the log region. For U z,rms + , this difference increases through the log region where it reaches a maximum, before decreasing towards the outer layer.
Agreement between the LES and DNS in the channel centre and close to the walls is good for all the profiles given in Figure 3. the constraint that particles must be the same size as, or smaller than, any fluid finite-volume cell.

Effects of Reynolds number on particle agglomeration
The initial particle positions were distributed randomly throughout the channel, with their initial velocity set to zero and with the particles coming in-line with local flow velocities with time. It is worth noting that the particle runs did not achieve steady state since the timedependent agglomeration of particles was the focus of the study. Particles were assumed to interact with turbulent eddies over a certain period of time, that being the lesser of the eddy lifetime and the particle transition time. Particles that moved out of the rectangular channel in the streamwise and spanwise directions were re-introduced back into the computational domain using periodic boundary conditions. All of the particles considered in this paper had identical intensive physical properties, as given in Table 1. The particle size, particle number and densities used were 150 m, 20,000 and 1000 kg m -3 , respectively, unless stated otherwise. The particle relaxation time is given by p = 24 p dp 2 /18 fCDRep, where p is the particle density, dp is the particle diameter, and CD is the drag coefficient (Shirolkar et al., 1996). The particle Stokes number, St, is defined as the non-dimensional particle relaxation time p + , i.e. St This section investigates the effects of turbulence on particle agglomeration, with three different flow Reynolds numbers considered, Re = 150, 300 and 590. The particle surface energies selected were 0.05 and 0.5 J m -2 due to their practical relevance (Afkhami et al., 2015), with the corresponding particle relaxation times, Stokes number and other relevant parameters given in Table 2. The integration time-step used for the particle simulations was 5.2 × 10 -7 s. This behaviour suggests that the higher flow turbulence in the Re = 590 case is initially responsible for creating a larger number of particle-particle interactions compared to the Re = 300 and 150 flows. The subsequent decline in the rate of particle contacts for the Re = 590 case is then indicative of an increase in the rate of particle contact breakage. This behaviour can be attributed to the initial conditions employed; as the particles accelerate and their velocity comes in line with that of the fluid, the greater flow turbulence causes the particles to encounter more fluid resistance (due to the drag forces acting in the opposite direction to the relative motion of the particle moving with respect to the surrounding fluid), with this increased resistance responsible for the increased rate of particle contact breakage in the higher Reynolds velocity to such an extent that the flow turbulence now causes particle-particle interactions. A linear increase in particle contact numbers then continues to about t = 0.05 s, after which an increasing divergence is seen between the various Reynolds number flows. This behaviour indicates a mechanism within the flow that advantages the particles exposed to higher Reynolds numbers in the formation of agglomerates. This occurs as a result of regions of high particle concentration and low particle mean velocity near the channel walls; in such regions the number of contacts formed is proportionally higher for particles of higher Reynolds number as the particles migrate to these regions at a faster rate. Moreover, the increased shear in the high Reynolds number flows increases the intensity of these turbulent regions, and therefore the particle velocity fluctuations and hence their number of interactions. Further analysis is desirable to establish a quantitative relationship between the particle fluctuating velocity and its impact on the formation of successful contacts. This, along with the dispersing behaviour of the particles and the regions in which contacts are formed, are discussed further below. For the Re = 150, 300 and 590 flows, respectively, there are 528, 635 and 1524 particle contacts in the flow at the end of the simulation for the higher surface energy case. These figures further demonstrate that increases in the flow Reynolds number dramatically enhance turbulence, and as a result particle agglomeration. It is thus again clear that the effects of turbulence are significant in creating successful particle-particle contacts, and that the flow Reynolds number is a key factor in determining particle agglomeration. respectively. In all three cases, particle agglomeration near to the walls can be attributed to the high particle concentration and the high turbulence levels in this region. Further scrutiny of the results shows that for the highest Reynolds number flow, particle agglomeration is roughly double that of the other flows at the channel centre. Close to the walls, the number of agglomerates is approximately equal, although the highest Re case tends to have the largest number of contacts in this region, with the number of contacts at all locations significantly greater than in the low Reynolds number case. This particle behaviour reflects the higher turbulence levels in the Re =590 flow, which drives the particles to regions of lower fluid mean velocity. Throughout the Re =590 flow, particle agglomeration is enhanced through high fluctuating velocities which cause a large number of particle-particle interactions, with peak levels approximately 30 wall units away from the solid boundaries. This effect is therefore most evident in the results for columns 2 and 3, which contain the highest agglomerate number.
These results indicate that flow Reynolds number is important to the occurrence of particleparticle contacts for high surface energies, as might be anticipated. In contrast, and at all Reynolds numbers, a lower surface energy leads to few agglomerates because of an increase in contact breakage due to the effects of flow turbulence. Figure 7 shows the instantaneous location of individual particles and contacts in the wallnormal direction, and their number at each location, at time t = 0.2 s. These plots provide more focus on locations very close to the wall, with results shown for equally spaced regions across half the channel height that are equivalent in size to the particle radius, with particle statistics combined within each of the zones of fluid considered. As already noted, generally the number of particles and contacts increases towards the walls. In the results of Figure 7, however, a sharp decrease is seen in the number of particles and contacts at the wall itself. Further scrutiny of the results shows that the point at which the particle number and contacts peak for the Re = 150 case is further away from the wall compared to the other two flows. This is related to the fact that the streamwise mean fluid velocity gradient is less steep towards the walls in the former case, resulting in a thicker boundary layer and lower levels of wall-normal turbulence intensity.
It is important to highlight that not all particle-particle interactions lead to the formation of agglomerates. The contact forces between colliding particles are based on the concept of contact mechanics, which takes plastic deformation of the particles into consideration. In the work presented only elastic deformation occurs, since the maximum stress does not reach the yield strength of the colliding particles (Bitter, 1963). The numerical model in DEM predicts the critical sticking, rebound and removal velocities, which are important parameters in determining the formation of agglomerates. During the collision of particles the normal contact force, Fn, and the tangential contact force, Ft, are considered. Friction is modelled using Coulomb's law of friction, with both the static friction coefficient and the kinetic (or sliding) friction coefficient accounted for. In general, for low shear forces there is little relative motion between the particles, therefore, the particles stick, whilst for high shear forces there is a more significant relative motion causing the particles to slip. The rolling friction was also taken into consideration in the present work. The normal contact force acts along the line joining the centres of the colliding particles, whilst the tangential contact force acts perpendicularly to that line. The contact force is defined by the collision phase and the relative velocity of the colliding particles. Collisions between particles can be divided into two consecutive phases, the approach and the restitution phase. The approach phase ends when the two bodies have a relative normal velocity equal to zero as a result of impact. According to Thornton and Ning (1998) the work required to break the contact between two particles is given by: If energy losses due to elastic wave propagation are neglected, the only work dissipated during a collision is the work done in separating the surfaces, We. Therefore: If the rebound velocity Vr = 0, then the impact velocity becomes equal to the sticking velocity, i.e. Vi = Vs, and accounting for the coefficient of restitution, e, and the critical velocity below which sticking occurs from (4) and (5), the sticking criterion becomes: If Vi > Vs then bounce occurs and (5) may be written as: The sticking velocity has been calculated using Error! Reference source not found.), and values for three levels of surface energy are given in Table 3. Note that here, the largest surface energy of 5 J m -2 corresponds to extremely cohesive particles that are of less practical relevance than the cohesive (0.05 J m -2 ) and very cohesive particles (0.5 J m -2 ) considered thus far. i.e. when two particles collide it is registered as one collision, regardless of how long the particles stay in contact for, with results collected over the duration of the collision. It should be noted that collisions can occur in between data write-outs and never register as contacts.
The normal impact velocity of two particles in a collision is Vr1 -Vr2; this value was calculated relative to the contact points and not the particle centres. In agglomerating systems the magnitude of the impact velocities may decrease until a steady state is reached in the total contact number. Data sampling should therefore ideally continue until this point is reached to give results that are fully representative of any quantitative differences. However, this requires very long computer run times, and is not considered relevant to the analysis presented below.
For Re =150 (Figure 8(a)), it is seen that the impact velocities range from less than 0.1 m s -1 to a maximum of 1.5 m s -1 , although with few collisions in the upper range. The velocity regions in which sticking occurs have been highlighted in this figure for different surface energies.
Based on these cut-off points the number of successful collisions (forming contact) are 45 and one order of magnitude, from 0.05 to 0.5 J m -2 , therefore results in twice the number of contacts.
For Re =300 (Figure 8(b)), the impact velocities range from around 0.002 m s -1 to a maximum of 3.1 m s -1 . The number of successful collisions is now 49 and 188 for the 0.05 and 0.5 J m -2 cases, respectively. Increasing the particle surface energy by one order of magnitude in this flow therefore results in almost an eight-fold increase in contact number. Lastly, for Re =590 (Figure 8(c)), the results show that the relative velocities range from less than 0.1 m s -1 to a maximum of 7.1 m s -1 . Based on the cut-off points the number of successful collisions are 15 and 227 for the 0.05 and 0.5 J m -2 particles, respectively. Increasing the particle surface energy by one order of magnitude in this case therefore results in a 55-fold increase in contact number.
Furthermore, it is seen that for low surface energy particles, the number of collisions within the range required for sticking is very low compared to that for the higher surface energy particles, and this explains the low number of agglomerates formed in the former case. From this analysis it is clear that the number of contacts formed is a function of the number of collisions, in relation to the sticking velocity between the particles, which varies with Reynolds number. Figure 9 shows the mean streamwise particle velocity, Vx, for the 0.05 J m -2 surface energy particles at all three shear Reynolds numbers at t = 0.2 s, obtained by spatially averaging over 100 non-uniform bins in the channel cross-section, and it is clear that the particle velocity increases with the flow Reynolds number, as would be anticipated. Profiles obtained for the 0.50 J m -2 particles are not shown as these gave very similar results and would not add to the discussion. Figure 10 shows the relationship between the particle position in the wall-normal direction for all three Reynolds numbers, plotted against the rms of the particle fluctuating velocity in the wall-normal (Vz,rms), spanwise (Vy,rms) and streamwise (Vx,rms) directions, again for 0.05 J m -2 particles at t = 0.2 s and averaged as indicated above. The locations of the points are plotted relative to the lower wall at 0, and the centre of the channel at 1. In the region closest to the walls, the particle streamwise, spanwise and wall-normal velocity fluctuations show peak values of approximately 0.44, 0.97 and 1.86 m s -1 , 0.08, 0.11 and 0.22 m s -1 , and 0.06, 0.08 and 0.18 m s -1 , respectively, with increasing Reynolds number. These results clearly illustrate the dramatic increase in fluctuating velocities with Reynolds number in regions where preferential agglomeration occurs. In the next region moving away from the wall (z/h = 0.01 -0.11), the sum of these fluctuating velocities is approximately maintained, despite the decrease in the streamwise fluctuations at the highest Reynolds number. The range in particle velocity fluctuations in the highest Re case demonstrates the significant influence flow turbulence can have on particle agglomeration in both these regions. Relating the results of Figure 9 and Figure   10, the difference in particle agglomeration between the various Reynolds number flows in the regions noted can be attributed to a combination of both the enhanced particle mean velocity and the particle velocity fluctuations. Finally, at the channel centre, the particle velocity fluctuation peak values are significantly reduced at all Reynolds numbers, thereby explaining the lower levels of particle collision and hence agglomeration in this region.
Lastly, Figure 11 shows the time evolution of the maximum value of the particle number density, np max , near the wall. The rationale for monitoring this quantity lies in the fact that the concentration close to the wall takes the longest time to reach a steady state within these flows.
The binning procedure used for obtaining np max was as follows. First, the channel cross-section was divided into 533 uniform regions equal to the particle radius. Next, the number of particles within each region was divided by the volume of that region at each time step to give the local concentration, np = np(s). Last, the maximum value of np amongst all regions was selected to give np max , with this value normalised by the total number of particles within the channel volume np total . The particle number density distribution will then be greater than 1 in crosssectional regions of the channel where the particles tend to accumulate, and ≈ 1 in regions where the particles are uniformly dispersed. The results clearly show that, starting from an initial distribution corresponding to a flat profile centered around np max = 1, the values for the Re = 150, 300 and 590 flows initially increase roughly linearly with time. In the case of the Re = 590 flow, the profile subsequently asymptotes to an almost constant value at around 0.19 s. In contrast, the other two profiles (Re = 150 and 300) continue to increase roughly linearly with time, although at later times it may be anticipated that they would also achieve steady values. Interestingly, little difference is seen between results obtained for the different particle surface energies at all Reynolds numbers. Overall, this behaviour suggests that the turbulence in the highest Reynolds number flows accelerates the particles at a faster rate in all directions, including towards the walls. Figure 12 gives values of the mean particle velocity in the streamwise direction for the 0.05 J m -2 surface energy case in the highest Reynolds number flow at a number of simulation times. In line with previous results, these figures provide a visual representation of how the particle velocity first increases with time up to approximately 0.05 s at it comes in line with the fluid velocity, and then decreases with time close to the walls, with the associated accumulation of particles near the channel solid surfaces.

Conclusions
Particles with identical physical parameters have been simulated in three channel flows with different levels of flow turbulence, achieved by increasing the Reynolds number of the flow, using a fully coupled LES-DEM approach. The particle diameter and surface energies selected were 150 µm, and 0.05 and 0.5 Jm -2 . Results suggest that the rate of agglomeration is strongly influenced, and increases, with the intensity of the flow turbulence, with most of the particleparticle interactions taking place at locations close to the channel walls and in regions of high turbulence where their agglomeration is aided both by the high levels of turbulence and the high concentration of particles. Additionally, the three different flows show that the rate of agglomeration is strongly influenced for high surface energy particles by, and increases with, the intensity of the flow turbulence. In contrast, for lower surface energy particles, the rate of agglomeration diminishes with an increase in flow turbulence intensity. It can be concluded that a combination of the effects of both the surface energy and fluctuating velocities is most significant in determining successful particle agglomeration in channel flows.

Acknowledgement
The authors would like to thank the Engineering and Physical Sciences Research Council of the UK for providing financial support through grant EP/P505593/1.     open red symbols = DNS (Moser et al., 1999)).