- Email: [email protected]
Three-dimensional simulation of multilayer particle deposition in an obstructed channel flow
Powder Technology 258 (2014) 134–143
Contents lists available at ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
Three-dimensional simulation of multilayer particle deposition in an obstructed channel flow Gregory Lecrivain a,⁎, Leopold Barry a, Uwe Hampel a,b a b
Institut für Fluiddynamik, Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany Technische Universität Dresden, Institut für Energietechnik, AREVA-Stiftungsprofessur für Bildgebende Messverfahren für die Energie- und Verfahrenstechnik, Germany
a r t i c l e
i n f o
Article history: Received 8 October 2013 Received in revised form 27 February 2014 Accepted 1 March 2014 Available online 10 March 2014 Keywords: Multilayer particle deposition Detached eddy simulation Self-organised criticality Dry granular mechanics Obstructed channel flow
a b s t r a c t A large variety of systems are subject to slow and lengthy processes of solid aerosol particle deposition in turbulent flows. As a result of a long exposure to deposition, a multilayer particle bed eventually forms over time. Notable examples are the formation of multilayer deposits in ventilation ducts, in nuclear reactors or on earth surfaces subject to atmospheric sedimentation. Simulations are of great importance to predict the multilayer deposition of solid aerosol particles. Theoretical models are quite limited since their complexity rapidly increases when the flow becomes turbulent and the surface geometry complex. The present study proposes a new three-dimensional approach to reproduce the growth of a multilayer deposit in a turbulent obstructed channel flow at Reynolds number Re = 10,000. Computational Fluid Dynamics and Computational Granular Dynamics are brought together to simulate 4 h of real deposition. A detached eddy simulation is employed to predict particle deposition while self-organised criticality is employed to reproduce the slow growth of the multilayer deposit. The three dimensional shape of the multilayer deposit matches remarkably well the experimental data. © 2014 Elsevier B.V. All rights reserved.
1. Introduction 1.1. Multilayer deposition of solid aerosol particles in turbulent flows Deposition, i.e. the process of attaching suspended aerosol particles from a fluid in turbulent motion to a wall surface, can either be monolayer or multilayer [1]. Monolayer deposition indicates that all deposited particles only adhere to a wall surface. Multilayer deposition indicates however a more complex particle arrangement: particles sit either on a wall surface or on top of other particles and form a multilayer particle bed, henceforth referred to as the multilayer deposit. A prolonged exposure to particulate matter inevitably leads to the formation of a multilayer deposit. The formation of multilayer deposit spans a wide range of industrial and non-industrial applications, such as particle deposition in ventilation ducts and pipelines [2] and particle deposition in graphite moderated High Temperature Reactors (HTR) [3,4]. There exists very scarce literature on multilayer deposition in turbulent flows. The complexity of the mechanisms involved renders the construction of computer models particularly difficult. Multilayer deposition involves the dynamic combination of two complex mechanisms: the transport of a particulate phase and the slow growth of a granular phase. Similarly to powders, the multilayer deposit exhibits a granular behaviour. One of the striking features of granular materials is the existence of a finite slope, given by the static angle of repose αs, below which ⁎ Corresponding author. Tel .: +49 351 260 3768; fax : +49 351 260 1 3768. E-mail address: [email protected] (G. Lecrivain).
http://dx.doi.org/10.1016/j.powtec.2014.03.019 0032-5910/© 2014 Elsevier B.V. All rights reserved.
the material is at rest [5]. A small stream of particulate matter falling down from a point source onto a flat surface will form a conical heap. The heap maintains under gravity a constant angle αs beyond which the excess material is lost through avalanches [6]. These two phenomena have largely been studied from two different standpoints. The review article of Guha [7] and the book of Pöschel and Schwager [8] testify to the increasing interest in the simulation of particle deposition and in the simulation of granular structures, respectively. An attempt is here made to combine computational fluid dynamics and computational granular dynamics to predict the complex growth of a multilayer deposit. The present work is a full three-dimensional extension of the work of Lecrivain et al. [9]. A fully unsteady CFD simulation and a new threedimensional granular pile model are brought together to simulate the multilayer deposition of solid aerosol particles in a turbulent obstructed channel flow over a fairly long period of time. 1.2. Granular pile models Similarly to sand piles the multilayer deposit obeys the laws of dry granular mechanics. Various methods each having cons and pros have been developed to simulate the heap formation (identical to sand pile formation), whose shape is characterised by the angle of repose αs. An extensive review of the current literature has shown that three approaches have been favoured in the recent years: the discrete element method, the differential approach and cellular automata. The discrete element method is a volume-based method. The heap is made of a large number of grains, which interact with one another
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
Nomenclature Greek αp αs αv Δh Δhth ΔP Δ δp ηd ν ∅ ρ τp ξ
Latin C0 dp e f g H h JS L Li = 1,2,3,4 m p Reb Rep S Ub ub ud uf up w yw
drag correction factor angle of repose airborne particle volume fraction local pile gradient stability threshold pressure drop cube root of cell volume inter-particle spacing deposited fraction kinematic viscosity of the fluid porosity of the multilayer deposit fluid density particle response time non-dimensional coefficient
135
introduced the concept of self-organised criticality. One of the first granular systems, for which cellular automaton was developed, involved the addition of one single grain at a time on a pile. These small perturbations to the system, i.e. the addition of grains, caused avalanche propagations. As illustrated in Fig. 1 a typical one-dimensional pile was regarded as a decreasing sequence of integers hi representing the height of the i-th column. The dynamics of the avalanches was simply described by two operators: adding a grain to the pile and relaxing the slope of the pile wherever the local gradient exceeded a critical stability threshold Δhth [17]. The stability of a column depended on the local gradient Δh = hi − hi + 1. If the local gradient exceeded the threshold Δh N Δhth, the grain in question toppled down in the next lower bin. The development of cellular automata has since been extended to simulate the three-dimensional heap formation on flat support [18]. Unlike the discrete element method and the differential approach, implementation of this surface-based method is relatively simpler. Self organised criticality will therefore be used in the present work to predict multilayer deposition in an obstructed turbulent channel flow. 2. Methods
airborne particle mass concentration particle diameter obstacle height friction coefficient gravity channel height local pile height particle mass flux channel length associated length of the i-th recirculation zone in the streamwise direction particle mass pitch Reynolds number based on bulk velocity Ub and channel height H particle Reynolds number particle-to-fluid density ratio bulk velocity of the fluid build-up velocity deposition velocity resolved component of instantaneous fluid velocity particle velocity kernel function distance to the nearest wall
Miscellaneous X Overbar denoting time or space averaging of arbitrary quantity X X+ Superscript denoting the normalisation of arbitrary quantity X with wall scaling (u⁎ for velocity variables, v/u⁎ for length variables, v/u⁎2 for time variables and u⁎3/v for acceleration variables)
through contact forces [10,11]. A fine tuning of the particle interaction model [12] and the relatively heavy computational cost [8] makes this method however unattractive for large systems. The differential approach is a surface-based method, which indicates that only the height h of the granular interface separating the multilayer deposit from the ambient gas is computed [13]. This analytical approach is however complex since it involves solving a system of partial differential equations for the dynamics of the surface height [14,15]. The emergence of cellular automaton to simulate avalanches on the surface of a granular heap goes back to 1988 when Bak et al. [16]
2.1. Model assumptions The present model offers a significant three dimensional extension of the work of Lecrivain et al. [9]. It therefore inherits all model assumptions previously made. 2.1.1. The flow is very dilute Experimental measurements in the wind tunnel revealed an airborne particle mass concentration of C0 = 2.75 g/m3, a value in line with real experiments previously performed in HTR [3]. With a particle-to-fluid density ratio equalled to S = 1800 and a mean particle diameter of dp = 5.3 μm, it led to a volume fraction of about αv = 1.2 × 10− 6 and a mean inter-particle spacing of δp/dp = 2500. For volume fraction αv b 10−6, Elghobashi [19] stated that the influence of particles on the gas phase and the particle interactions can be neglected.
2.1.2. The deposition process is very slow The multilayer deposit grows at some millimetres per hours while a particle crosses the channel at some metres per second. The build-up velocity ub, i.e. the velocity at which the multilayer deposit grows, is by several orders of magnitude smaller than the particle velocity up. During its life time in the channel, a particle sees a fairly constant shape of the multilayer deposit. The growth of the multilayer deposit is therefore a quasi-static process.
Fig. 1. Addition of one grain on a one-dimensional pile. The grain toppled in the next lower bin until it reaches a stable position.
136
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
2.1.3. Deposition is not impacted by the amorphousness of particles The experimental investigation of Kvasnak and Ahmadi [20] showed that the deposition rate of fibre-like particles in a horizontal duct is very close to that of spherical particles with identical response time. The deposition of amorphous graphite particles can thus be numerically reproduced using equivalent volume spheres without a significant loss of precision. 2.1.4. Particle detachment off the multilayer deposit does not take place Previous investigations performed by references [1] and [21] suggested that the diameter of re-entrained particles was up to an order of magnitude larger than that of the deposited particles, which suggests that particles resuspend off the multilayer deposit in the form of agglomerates rather than as sole entities. Agglomerates, because of their large size, are more easily re-entrained than single fine particles [22,23]. Lazaridis and Drossinos [24] demonstrated that the adhesion force between particles in the multilayer deposit is less than that between a spherical particle and a flat surface. Under very specific flow conditions, the two mechanisms of detachment and redeposition may even occur simultaneously [25,26]. To avoid the detachment of agglomerates, the friction velocity is low at the granular interface. Resuspension experiments performed by Reeks and Hall [22] showed that graphite particles in the size range dp = 6-30 μm initially at rest on flat plate started to reenter the flow beyond the critical friction velocity u∗c = 0.3 m/s. Barth et al. [27] detected particle detachment off a the multilayer deposit for friction velocities greater than u∗c = 0.15 m/s. In the present experiment the time-averaged friction velocity in the cavity did not exceed the value of u∗max = 0.1 m/s. 2.1.5. Depth compaction of the multilayer deposit does not occur Simulation of multilayer particle bed structures performed by Schmidt [28] and recent X-ray measurements of multilayer deposits performed by Barth et al. [27] under atmospheric conditions did not revealed any compaction. 2.2. Channel geometry The present simulation seeks to numerically reproduce the multilayer deposition experiment performed by reference [27]. The flow is wallbounded in a horizontal square duct having a rib-roughened bottom wall surface. The channel height, the rib height and the rib pitch are denoted by H, e and p, respectively. The rib pitch, defined as the distance from one leading edge of the rib to the next leading edge, is set to p = H. The height of the rib obstructs the channel height by e/H = 10 %. In the experimental tests 15 obstacles were placed on the bottom wall. To cut down on the computational cost the flow is however simulated around one single obstacle. The inlet and outlet surfaces coloured in grey in Fig. 2 are enforced a translational periodic boundary condition. 2.3. Numerical procedure A quasi-static approach is suggested to speed-up the three dimensional simulation. It would be extremely time-consuming to continuously simulate hours of deposition in the present obstructed channel flow. Similarly to the two-dimensional multilayer deposition model of Lecrivain et al. [9], the growth of the multilayer deposit is therefore discretised into a series of successive static events. Over the course of each event the geometry of the multilayer deposit remains constant. The simulation process throughout each event is threefold: 1. The transport and deposition of the particles by the gas phase is simulated for a few seconds, 2. The build-up velocity of the multilayer deposit is calculated from the monolayer deposition and linearly extrapolated to reproduce a larger time of exposure, and 3. The deposited mass of graphite is redistributed through avalanches to satisfy the granular properties of the multilayer deposit. The granular interface separating the multilayer deposit from the gas phase acts as a wall boundary in the simulation. The computational domain, in which the gas phase evolves during each event, is bounded
Fig. 2. Schematic of the obstructed channel. A periodic boundary condition is enforced on the two surfaces coloured in grey.
by the channel and the granular interface. Prior to each flow simulation, a Boolean operation is performed to subtract the volume beneath the granular interface from the original channel geometry. In the very first event deposition has not yet taken place. In that case the obstacle is simply subtracted from the channel. 2.4. Simulation of the gas phase Reynolds Averaged Navier-Stokes (RANS) models usually give acceptable results for most engineering applications. RANS simulations are fast, robust but often fail to accurately capture the unsteady flow structures, especially in regions of massive flow separation [29]. Large Eddy Simulation (LES) appears to be a more natural candidate for the unsteady simulation of complex flow separations. The potential accuracy of LES is generally well acknowledged in the literature. The great computational cost resulting from a fine mesh resolution near solid boundaries makes it however unattractive for practical applications. Detached Eddy Simulation (DES) overcomes these wall restrictions by considerably cutting down the cost of a computation. It is a hybrid method in which LES is employed in core and separated flow regions while near the wall RANS model has complete control over the solution. This hybrid concept has been successfully applied to simulate the turbulent flow in obstructed channels [30]. For the present investigation, the DES formulation available in the OpenFOAM library and similar to that proposed by Spalart et al. [31] is used. To accomplish the transition from nearwall RANS-based simulation to LES treatment of the interior flow in one formulation, the nearwall distance yw is replaced by the distance d, defined as e d ¼ minðyw ; C DES ΔÞ
ð1Þ
which act as a RANS model for Δ ≪ yw and as a sub-grid scale model for Δ ≫ yw. The symbol Δ is defined as the cube root of the cell volume and CDES is an empirical constant calibrated to a value of 0.65. More details on the exact implementation can be found in the work of Islam et al. [32]. 2.5. Simulation of the particulate phase The high particle-to-fluid density ratio S and a particle diameter greater than dp N 1 μm implies a negligible effect of the Saffman lift force [33] and of the Brownian force responsible for the molecular
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
137
diffusion [34], respectively. The acceleration of the particle is therefore only influenced by its drag and the gravity. The equation of motion for a point-like particle evolving in the obstructed channel flow reads
given by averaging the positions xi of its N connected sites. Mathematically the space-averaged position of a site is given by
dup u f −up ¼ þ g: dt τp =α p
x¼
ð2Þ
In above equation, up is the particle velocity and uf is the resolved component of instantaneous fluid velocity derived from the DES simulation. The contribution of the sub-grid scale velocity component is here left out since its effect on particle deposition is questionable [35]. The response time τp of a spherical particle is defined as τp ¼
Sd2p 18ν
ð3Þ
where ν is the fluid kinematic viscosity. The drag correction factor αp derived from the semi-empirical formula of Schiller and Naumann is a non-linear function of the particle Reynolds number Rep. Its formulation is given by αp ¼ 1 þ
0:687 0:15Rep :
ð4Þ
The term Rep is the particle Reynolds number based on the relative particle velocity to the fluid velocity and the particle diameter. The solid aerosol particles are spherical in the simulation. The in-house implementation of Eq. (2) is the same as that used in the work of Lecrivain and Hampel [36]. A particle deposits whenever the distance from the particle centre to the nearest wall is smaller than the particle radius. It is assumed that particle rebound does not occur here. In reality multiple rebounds upon collision with a wall surface occur without resulting in deposition, especially at low impact angle [37]. In a gravity-dominated deposition regime, as is the case here, the deposition velocity ud of a particle can be estimated with the formulation of Wood [38]. For a 20 μm graphite particle, the deposition velocity equals ud = g × τp = 0.03 m/s. This deposition velocity is far below the typical impact velocities used for the determination of the rebound coefficient of aerosol particles. Wall et al. [39] experimentally measured the rebound coefficient of 3– 7 μm spheres impacting a wall at velocities ranging from 1 to 20 m/s, for which an impact velocity of 1 m/s did not result in a rebound. Particle rebound, although possible, is unlikely. The low deposition velocity, the low friction velocity and the relatively strong adhesion forces tend to prevent this mechanism. 2.6. Simulation of the granular phase 2.6.1. Lattice construction An unstructured lattice is here used to represent the three dimensional shape of the granular interface. Each site of the lattice stores the height h of the multilayer deposit. At the initial time the shape of the lattice corresponds to that of the roughened channel surface, i.e. h = e for the obstacle else h = 0 for the channel floor. The procedure of Moukarzel and Herrmann [40] is employed to enforce a certain degree of disorder in the lattice, which will avoid an anisotropic distribution of the graphite during the simulation of avalanches. A chessboard-like structured square lattice is first defined. Within each square a random point with uniform distribution, henceforth referred to as a site, is generated. To create the connectivity of each site a constrained twodimensional Delaunay triangulation is performed. The Delaunay triangulation of these sites ensures that the circumcircle associated with each triangle contains no other site in its interior. Boundary sites lying on the bounding lines (walls and periodicity) and on the cliffs of the obstacle have constant coordinates either in the streamwise or spanwise direction. To avoid occurrence of high skewness elements in the lattice, a Laplacian smoothing operator is used. In this simple and yet effective smoothing technique, each internal site is moved to a new position x
N 1X x N i¼1 i
ð5Þ
where N is the number of adjacent sites connected. Collapse of edges smaller than an arbitrarily small value is then performed to finalise the lattice smoothing. The sites forming the obstacles are then shifted vertically by the distance e to form the obstacle. Following this extrusion, the sites located on the two cliffs of the obstacle show a discontinuous behaviour. Each node separating the obstacle from the channel surface now possesses two heights. The above mentioned steps necessary for the construction of smooth unstructured lattice are shown in Fig. 3. The number of sites was deliberately set to a low value for illustration purposes. 2.6.2. Boundary treatment of the lattice 2.6.2.1. Cliff of the obstacle. Each site located on one of the two cliffs of the obstacle must be dealt with cautiously because of a discontinuous behaviour. Each site is decomposed into a lower and upper component. The upper component belongs to the obstacle while the lower component belongs to the channel floor. Upon avalanche the upper component becomes passive and directly conveys the mass it receives to the lower site component. As a result the lower component moves upwards while the upper component remains immobile throughout the mass distribution. As soon as the height of the lower component exceeds that of the upper component by an arbitrarily small distance, the site loses its discontinuity and behaves exactly like any other regular site. 2.6.2.2. Periodicity. During the creation of the lattice, an additional constraint on the Delaunay triangulation is enforced to ensure one-to-one periodic connectivity. The spanwise position of the sites on the inlet line exactly corresponds to that of the sites on the outlet line. During the avalanche process the two paired sites with one-to-one connectivity behave as a single entity: they share the same neighbours and have the same height. 2.6.2.3. Channel walls. In reality, the two side walls of the duct prevent the graphite material from leaving the channel. No numerical treatment is actually needed to prevent the graphite from leaving the lattice, since the material can only be transferred from one site to another. 2.6.3. Avalanche A real function h(x, z) describing the local pile height is associated to each site, where the variables x and z respectively refer to the streamwise and spanwise directions of the lattice. The stability of a site in a given direction solely depends on the nearest neighbours, whose connectivities were created during the Delaunay triangulation. A site is said to be unstable or critical in the i-th direction if the difference in height Δh = h − hi between the site in question and the nearest site in the i-th direction exceeds a threshold Δhth i . The threshold is a function of the distance between the two sites and the angle of repose αs qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi th Δhi ¼ tan ðα s Þ ðx−xi Þ2 þ ðz−zi Þ2 :
ð6Þ
If the site is unstable in the i-th direction, a small amount of mass topples down the adjacent site. Mathematically it means ( th Δhi NΔhi ?
h hi
¼ ¼
h hi
th
þ ξ Δhi th ; − ξ Δhi
ð7Þ
where ξ is a non-dimensional coefficient ranging from 0 to 1. In the original model of Puhl [41] cubical sand grains were toppled (Fig. 1). It
138
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
a
b
c
d
Fig. 3. Lattice construction. (a) Random point generator and constraint lines, (b) constraint Delaunay triangulation, (c) Laplacian smoothing, (d) extrusion of obstacle.
would be however very time-consuming to simulate the toppling of each grain. The above coefficient is therefore introduced to reproduce a small mass transfer. In this case we had ξ = 0.1 which was found to yield a good compromise between accuracy and execution time. To simulate avalanche propagation, the algorithm successively updates the whole lattice until all sites become stable. The order, in which the sites and the nearest neighbours are tested, is fixed during each sweep. Random permutations of the directions did not show any improvements. 2.6.4. Test cases 2.6.4.1. Influence of lattice anisotropy. To emphasise the importance of the disorder in the lattice, the growth of heaps is reproduced on a flat plate. Two point sources, a first source feeding material on a structured lattice and a second feeding on an unstructured lattice, are simulated. The arrangement of the sites in the structured lattice is fairly simple.
Each site possesses nearest neighbourhoods in the four cardinal directions. Fig. 4 shows the two piles obtained. It can be seen that a pyramid is formed on the structured lattice while a conical pile is formed on the unstructured lattice. During the propagation of avalanches on the structured lattice the two directions of the structured lattice are preferred, which is why a pyramid is formed. One can clearly see that the disorder removes the two preferred directions. In reality, the shape the conical heap rarely presents sharp edges. Herrmann [42] showed that piles normally have a logarithmic tail. Therefore the Laplacian smoothing is used in the vertical direction with a very small number of iterations. It smooths the tail of the heap and removes the high frequency oscillations of the granular interface h. 2.6.4.2. Interaction of growing piles. The present granular pile model should satisfy the rules governing the growth of multiple piles [9,43]. For the sake of simplicity, let us assume two point sources feeding
Fig. 4. Effect of lattice structure on heap shape (αs = 45°).
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
139
Fig. 5. Heap growth on an unstructured lattice (αs = 45°). Top: two heaps growing at same velocity, middle: two heaps growing at different velocities, bottom: single heap growing on obstacle.
material with the corresponding fluxes J1 and J2. The two heaps initially grow independently as long as they are far apart and do not interact. However upon collision the final configuration will depend on the flux difference. If the two sources feed material at the same rate (J1 = J2) the configuration reaches an unstable equilibrium and the two heaps grow at the same velocity. If however one source feeds more material (J1 N J2), the faster-growing heap gradually encroaches the slowergrowing heap. One single heap centred at the position of the source J1 emerges from the initial configuration. The faster-growing heap eventually completely covers the slower growing heap. The evolution of this two-heap system is illustrated in Fig. 5. 2.6.4.3. Growth on obstacle. Experimental observations have shown that, when a pile reaches the border of a support, a runway forms and the grains falling from the point source move down that way to the edge and leave the support [18]. Let us assume a point source placed above the obstacle in the channel. A heap initially grows until it meets the edge of the obstacle. Upon contact with the edge, the heap stops growing and the material falls down the cliff. Through a “cascade effect” the
material is transferred to the channel bottom. A second heap starts to grow in the corner. When the summit of the offspring heap reaches the upper cliff, the two heaps join together and then behave as a single entity (Fig. 5). 2.7. Calculation of the build-up velocity The complexity in the present section lies in the transformation of a discrete deposition cloud derived from the particle-laden flow simulation into a continuous surface field of the build-up velocity ub. The build-up velocity, i.e. the velocity at which the multilayer deposit grows is defined as follows ub ¼
JS ; ρð1−ϕÞ
ð8Þ
where JS is the local particle mass flux on a wall surface, i.e. the mass of deposited particles per unit surface and per unit time. A simple method to work out the build-up velocity of each site of the lattice would be to
Fig. 6. Instantaneous and time-averaged flow field at the mid-section.
140
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
centre a circle of a fixed size and sum up the mass of each deposited particle within the circle. A much smoother distribution of the build-up velocity across the lattice can be achieved via a weighting function w (also referred to as kernel function) which assigns a weight to each deposited particle based on its distance from the circle centre [44]. The local particle mass flux of a site is calculated as follows J S ðx; zÞ ¼
N X 1 wðx; zÞmi ; 2 t s πr i¼1
ð9Þ
Table 2 Comparison of recirculation lengths with data from the literature.
DES (present study) Experiment [45] LES [46]
L1 e
L2 e
L3 e
L4 e
1.06 1–1.5 1.62
0.81 0.6–0.9 0.91
0.88 0.2–0.3 0.5
4.12 3.7–3.9 3.45
where mi is the mass of the i-th particle, r is the smoothing radius, ts is the duration of the particle-laden flow simulation and N is the total number of deposited particles on the wall surface. Several types of kernel functions are commonly used in statistics. We choose here to use the standard Gaussian density, with which the weights smoothly decrease towards zero as one moves away from the target point. For deposition on a surface the smoothing function reads
downstream of the obstacle (L3) and to the largest recirculation area downstream of the obstacle (L4), respectively. Table 2 compares the four recirculation lengths with the experimental data of Casarsa and Arts [45] and the LES of Lohász et al. [46]. All but the third associated length match relatively well the experimental data. It can be seen that the present DES offers a great improvement over a conventional RANS simulation, which was found unable to capture the second recirculation zone on top of the obstacle [9]. The time-average pressure drop over one pitch interval ΔP was also computed for validation purposes. The friction factor f is defined as
i 1 h 2 2 wðx; zÞ ¼ exp − 2 ðx−xi Þ þ ðz−zi Þ : 2r
f ¼
The weight function w is typically indexed by the smoothing radius r. A small value of r emphasizes high-frequency features whereas a greater value of r will result in a smoother distribution. Monolayer deposition on vertical walls does take place. However, the wall-induced adhesion forces are not strong enough to allow multilayer deposition on vertical walls. The graphite particles hitting the vertical walls fall off and therefore contribute to the build-up velocity of the sites in contact with the channel side walls and the cliffs (provided the site is discontinuous). For sites connecting the granular interface with a vertical wall, the flux is calculated slightly differently. The disc area associated with the kernel function in Eq. (9) is replaced with 2 ⋅ r ⋅ δ, where δ is the average spacing of two connected sites. To simulate the multilayer deposition, the monolayer deposition is extrapolated in time. The height h(x, z) of the granular interface is obtained from integration of the buildup velocity field as follows dh ¼ ub : dt
ð11Þ
3. Results and discussions 3.1. Turbulent flow field The computational domain, in which the gas phase evolves, has about 0.5 × 106 cells. The typical distance from the first cell centre to + the closest wall is about y+ w = 0.7 around the obstacle, yw = 1.5 on + the channel floor and yw = 3 on the side and top walls. While the resolution of the grid could be finer, the number of grid nodes was kept to a minimal value to ensure a relatively fast simulation of the particle laden flow. Fig. 6 shows the instantaneous and time-averaged flow field through the cross-sectional symmetry plane of the initial channel geometry. To obtain the time averaged flow field, the fluid statistics were averaged in time for 40 flow-through times. Four flow recirculation zones exhibiting a strong rotating nature can be observed near the obstacle. The lengths Li = 1 … 4 are associated to the recirculation areas, upstream of the obstacle (L1), on top of the obstacle (L2),
ð12Þ
where ρ the fluid density. The relative effect of obstacle on the pressure drop can be normalised with the friction factor f0 of smooth circular pipe. The Blasius expression correlates the friction factor with the bulk Reynolds number as follows f0 = 0.046 Re−0.2 . The present DES delivb ered a value of f/f0 = 2.91. It compares well with the value of f/f0 = 3.0 obtained by inserting the present values for obstacle dimension and spacing in the model of Tran [47]. The present ratio is also close to the experimental value of f/f0 = 3.2 taken from the work of Lai et al. [48] at similar Reynolds number. 3.2. Monolayer deposition 3.2.1. Deposition of monodisperse particles The present subsection seeks to compare the deposition rates with existing data available in the literature. While accurate literature data on particle deposition rates in complex geometries is scarce, a few attempts have been made to compute the deposition rates as a function of the particle size in non-obstructed channel flows. The theoretical deposited fraction ηd of fully mixed particles with uniform concentration in a channel bounded by four walls can be estimated using the formulation of Matida et al. [49] L ud ηd ¼ 1− exp −4 ; H Ub 1
0.8
ð13Þ
1
50 dp [μm]
10
Wood Eq. DES
0.6
0.4
0.2
Table 1 Input data for the simulation. Friction velocity Bulk Reynold number Porosity Angle of repose Particle size
H ΔP ; 2 p ρU 2b
ηd
ð10Þ
u* Reb ϕ αs dp
0.06 m/s 10,000 90% 35° 1..20 μm
0 0.001
0.01
0.1
1
τ+p Fig. 7. Deposited fraction against particle response time.
10
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
Discrete deposition cloud
141
Continuous deposition flux
0.8
6
0.6
4
0.4
2
0.2
z/e
8
J/Jmax
1
10
0
0 0
3
x/e
6
9
0
3
x/e
6
9
0.8
6
0.6
4
0.4
2
0.2
y/e
8
0
J/Jmax
1
10
0 5
6
7
8
9
0
1
2
3
4
x/e
5
6
7
8
9
0
1
2
3
4
x/e
Fig. 8. Deposition flux $J_S$ on bottom surface (top) and channel side (bottom).
where ud is the deposition velocity. Wood [38] derived from various experimental measurements a fairly simple expression for the estimation of the deposition velocity þ
þ2
þ þ
ud ¼ Kτ p þ g τ p :
ð14Þ
The first and second terms correspond, respectively, to the eddy diffusion impaction and to the gravitational sedimentation. These two deposition mechanisms where found to be predominant for the present study [50]. The constant K is an empirical value set to K = 4.5 ⋅ 10−4. The deposition is here simulated in an obstructed channel with a length-to-height ratio L/H = 10. To reproduce the length of the channel,
a counter is associated to each particle. Whenever a particle crosses the outlet periodic surface the counter is incremented by one. The simulation is run until all particles have either deposited or have shown a counter greater than 10. Particles with a counter greater then 10 are then removed from the cloud since they are assumed to have left the channel. For each class, about 15,000 particles of constant size are initially placed on a cross-sectional plane normal to the streamwise direction, placed at the position x = H/2. Particles are randomly injected across the planar surface stretching from the obstacle to the channel ceiling. The deposited fraction ηd = Nd/N0 is calculated at the end of the particle-laden flow simulation, where Nd is the number of deposited particles. The particle diameter is varied from 1 to 50 μm. It can be seen from Fig. 7 that the evolution of the deposited fraction compares well with the empirical expression of Wood [38]. The results of the present simulation are very satisfactory, since estimations of the deposition velocity in turbulent flows show large discrepancies often covering as much as one order of magnitude [51]. 3.2.2. Deposition of polydisperse particles
Fig. 9. Simulation of the granular interface after 4 h of deposition.
3.2.2.1. Deposition on channel floor. The size distribution of the particle cloud exactly corresponds to that used in the referenced experiment [27]. The particle equivalent aerodynamic diameter, whose mean value equals 5.3 μm, varies from 1 μm to 20 μm. Similarly to the deposition of monodisperse particles, the flow simulation is initially run during 10 flow-through times to remove artificial flow structures, originating from the initial flow perturbation. About 106 particles are then carried by the flow domain during 40 extra flow-through times. The cloud of deposited particles is used to compute the particle mass
142
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
y/e
0.4
t=1h
Sim t=2h Exp
t=3h
t=4h
0.2
0.0
y/e
0.6
0.4
0.2
0.0 0
3
6
0
3
6
x/e
9
x/e
Fig. 10. Streamwise build-up in the cavity ($z/e = 0.5$) — comparison between simulation and experiment.
flux on the channel bottom, on the obstacle surfaces and on the two channel sides. Fig. 8 illustrates the deposition cloud in the cavity. Not all deposited particles are displayed for illustration purposes. The transformation of the deposited cloud into a smooth surface particle mass flux using Eq. (9) reflects a very inhomogeneous flux in streamwise and spanwise directions across the channel floor. The effect of the recirculation zones on the deposition rate is here of particular importance. The channel bottom area coloured in dark blue and located downstream of the obstacle (x/e b 3) clearly indicates the lower deposition rates due to the largest flow separation. The channel bottom area coloured in dark red and located upstream of the obstacle (x/e N 8) is subject to a higher deposition flux. This is because of the recirculation area associated with the length L1, whose downward flow causes the particles to collide with the channel bottom. As one comes very close to a vertical wall the smoothing function w (Eq. (10)) normally shows a non-physical decreasing behaviour. It is a known limitation of the present smoothing technique because the smoothing circle of radius r goes beyond the lattice boundaries. Although ad-hoc techniques are available to correct the relatively small boundary values [52], their use is not justified here because the flux originating from the deposition on the vertical walls is much greater than the surface flux on the channel floor. For this reason, the surface flux very close to a vertical wall equals the closest flux value, whose value is not influenced by the boundary. 0.3
t=1h
3.2.2.2. Deposition on side wall. The multilayer deposition on vertical surfaces of the obstructions occurs neither in the simulation nor in the experiment. The mass depositing on vertical surfaces therefore solely contributes to the layer build-up on horizontal surfaces through a transfer of material (“cascade effect”). Therefore the deposition flux JS on the side walls solely varies in the streamwise direction. The large turbulence intensity next to the obstacle causes a heavy deposition flux on the vertical sides, while further away from the obstacle the deposition lowers. One may observe that the deposition flux is five times as large as that near the rib in far-off regions. 3.3. Multilayer deposition Simple measurements of the angle of repose for piles of the present graphite material delivered a value of αs = 35°. The experimental value worked out here agrees well with literature data. Typical values of the angle of repose for dry powders vary from 30° to 40° [53]. All input parameters used in the final simulation are summarised in Table 1. In this work, each event reproduces 30 min of real deposition which corresponds to a vertical geometric change of about 1% of the channel height. The discretisation step chosen here reaches a compromise between solution accuracy and computational cost. Each event involves the transport of about 106 particles in a time span of 4 s. Each particle-laden
t=2h
Sim Exp
y/e
0.2 0.1 0.0
t=3h
t=4h
0.4
y/e
0.3 0.2 0.1 0.0
0
1
2
3
z/e
4
0
1
2
3
4
z/e
Fig. 11. Spanwise build-up in the cavity ($x/e = 0.6$) — comparison between simulation and experiment.
5
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
flow simulation was performed on 2 × 64 processors of AMD 16-Core Opteron 2.3 GHz and required approximately 12 h of CPU time. Robustness of the granular model made possible the simulation of avalanches on a desktop computer. No more that a few hours were necessary for the redistribution of the deposited material across the lattice. All together about two weeks of simulation were necessary to reproduce 4 h of real multilayer deposition. Fig. 9 illustrates the shape of the lattice after 4 h of real deposition. The lattice is coloured by altitude. The height of the two black lines drawn on the granular interface in the streamwise direction (x) and in the spanwise direction (z) is compared with experimental data in Figs. 10 and 11. To observe the build-up progress the experiment was paused each hour and the height of the multilayer deposit in the cavity was measured using a laser scanner. The laser here is used to avoid direct contact that could potentially harm the multilayer deposit. It is found that the simulation matches remarkably well the experimental data. The shape of the multilayer deposit in the cavity and near the vertical walls is well captured. The slope of the granular interface near the vertical walls corresponds exactly to the value of the angle of repose. 4. Conclusions The present study has two novel aspects: 1. A fast and robust threedimensional granular pile model for multilayer growth on complex supports was suggested and 2. the dynamic coupling of the granular pile model was coupled with computational fluid dynamics. The combination of self-organised criticality and detached-eddy-simulation was successfully combined to simulate the multilayer deposition in a turbulent obstructed channel flow. The three dimensional shape of the multilayer deposit matched remarkably well the on-site experimental data. The outcome of this study has a direct impact on various applications such as atmospheric deposition and aerosol particle deposition in ventilation systems and in High Temperature Reactors. Acknowledgements This work was supported by the European Commission under the Grant 249337 of the “Thermal Hydraulics of Innovative Nuclear Systems” research project. References [1] H. Friess, G. Yadigaroglu, Modelling of the resuspension of particle clusters from multilayer aerosol deposits with variable porosity, J. Aerosol Sci. 33 (2002) 883–906. [2] Q. Li, J. Song, C. Li, Y. Wei, J. Chen, Numerical and experimental study of particle deposition on inner wall of 180° bend, Powder Technol. 237 (2013) 241–254. [3] R. Moormann, Fission product transport and source terms in HTRs: experience from AVR pebble bed reactor, Sci. Tech. Nucl. Install. 2008 (2008) 1–14. [4] M. Kissane, A review of radionuclide behaviour in the primary system of a veryhigh-temperature reactor, Nucl. Eng. Des. 239 (2009) 3076–3091. [5] H. Herrmann, Statistical models for granular materials, Phys. A 263 (1999) 51–62. [6] Y. Grasselli, H. Herrmann, On the angles of dry granular heaps, Phys. A 246 (1997) 301–312. [7] A. Guha, Transport and deposition of particles in turbulent and laminar flow, Annu. Rev. Fluid Mech. 40 (2008) 311–341. [8] T. Pöschel, T. Schwager, Computational Granular Dynamics: Models and Algorithms, Springer, 2010. [9] G. Lecrivain, S. Drapeau-Martin, T. Barth, U. Hampel, Numerical simulation of multilayer deposition in an obstructed channel flow, Adv. Powder Technol. 25 (2014) 310–320. [10] T. Shinbrot, M. Zeggio, F. Muzzio, Computational approaches to granular segregation in tumbling blenders, Powder Technol. 116 (2001) 224–231. [11] H. Zhu, Z. Zhou, R. Yang, A. Yu, Discrete particle simulation of particulate systems: theoretical developments, Chem. Eng. Sci. 62 (2007) 3378–3396. [12] J. Baxter, U. Tüzün, J. Burnell, D.M. Heyes, Granular dynamics simulations of twodimensional heap formation, Phys. Rev. E. 55 (1997) 3546–3554. [13] J. Bouchaud, M. Cates, J. Ravi Prakash, S. Edwards, A model for the dynamics of sandpile surfaces, J. Phys. I 4 (1994) 1383–1410. [14] L. Prigozhin, Variational model of sandpile growth, Eur. J. Appl. Math. 7 (1996) 225–235. [15] K. Hadeler, C. Kuttler, Granular matter in a silo, Granul. Matter 3 (2001) 193–197. [16] P. Bak, C. Tang, K. Wiesenfeld, Self organised criticality, Phys. Rev. A 38 (1988) 364–375.
143
[17] M. Aschwanden, Self-organized Criticality in Astrophysics: The Statistics of Nonlinear Processes in the Universe, Springer, 2011. [18] H. Puhl, On the modelling of real sandpiles, Phys. A: Stat. Mech. Appl. 182 (1992) 295–319. [19] S. Elghobashi, On predicting particle-laden turbulent flows, Appl. Sci. Res. 52 (1994) 309–329. [20] W. Kvasnak, G. Ahmadi, Fibrous particle deposition in a turbulent channel flow — an experimental study, Aerosol Sci. Technol. 23 (1995) 641–652. [21] S. Matsusaka, H. Masuda, Particle reentrainment from a fine powder layer in a turbulent air flow, Aerosol Sci. Technol. 24 (1996) 69–84. [22] M. Reeks, D. Hall, Kinetic models for particle resuspension in turbulent flows: theory and measurement, J. Aerosol Sci. 32 (2001) 1–31. [23] G. Ziskind, Particle resuspension from surfaces: revisited and re-evaluated, Rev. Chem. Eng. 22 (2006) 1–123. [24] M. Lazaridis, Y. Drossinos, Multilayer resuspension of small identical particles by turbulent flow, Aerosol Sci. Technol. 28 (1998) 548–560. [25] I. Adhiwidjaja, S. Matsusaka, H. Tanaka, H. Masuda, Simultaneous phenomenon of particle deposition and reentrainment: effects of surface roughness on deposition layer of striped pattern, Aerosol Sci. Technol. 33 (2000) 323–333. [26] S. Matsusaka, W. Theerachaisupakij, H. Yoshida, H. Masuda, Deposition layers formed by a turbulent aerosol flow of micron and sub-micron particles, Powder Technol. 118 (2001) 130–135. [27] T. Barth, M. Reiche, M. Banowski, M. Oppermann, U. Hampel, Experimental investigation of multilayer particle deposition and resuspension between periodic steps in turbulent flows, J. Aerosol Sci. 64 (2013) 111–124. [28] E. Schmidt, Simulation of three-dimensional dust structures via particle trajectory calculations for cake-forming filtration, Powder Technol. 86 (1996) 113–117. [29] H. Lübcke, S. Schmidt, T. Rung, F. Thiele, Comparison of LES and RANS in bluff-body flows, J. Wind. Eng. Ind. Aerodyn. 89 (2001) 1471–1485. [30] A. Viswanathan, D. Tafti, Detached eddy simulation of turbulent flow and heat transfer in a two-pass internal cooling duct, Int. J. Heat Fluid Flow 27 (2006) 1–20. [31] P. Spalart, W. Jou, M. Strelets, S. Allmaras, Comments of feasibility of LES for wings, and on a hybrid RANS/LES approach, International Conference on DNS/LES, Ruston, Louisiana, 1997. [32] M. Islam, F. Decker, E. de Villiers, A.e.a. Jackson, Application of Detached-Eddy Simulation for Automotive Aerodynamics Development, SAE Tech. Paper 2009. (SAE 2009-01-0333). [33] M. Soltani, G. Ahmadi, Direct numerical simulation of particle entrainment in turbulent channel flow, Phys. Fluids 7 (1995) 647–657. [34] M. Horn, H. Schmid, A comprehensive approach in modeling Lagrangian particle deposition in turbulent boundary layers, Powder Technol. 186 (2008) 189–198. [35] M. Salmanzadeh, M. Rahnama, G. Ahmadi, Effect of sub-grid scales on large eddy simulation of particle deposition in a turbulent channel flow, Aerosol Sci. Technol. 44 (2010) 796–806. [36] G. Lecrivain, U. Hampel, Influence of the Lagrangian integral time scale estimation in the near wall region on particle deposition, J. Fluid. Eng.-T. ASME 134 (2012) 1–6. [37] N. Konan, O. Kannengieser, O. Simonin, Stochastic modeling of the multiple rebound effects for particle–rough wall collisions, Int. J. Multiphase Flow 35 (2009) 933–945. [38] N. Wood, A simple method for the calculation of turbulent deposition to smooth and rough surfaces, J. Aerosol Sci. 12 (1981) 275–290. [39] S. Wall, W. John, H. Wang, S. Goren, Measurements of kinetic energy loss for particles impacting surfaces, Aerosol Sci. Technol. 12 (1990) 926–946. [40] C. Moukarzel, H. Herrmann, A vectorizable random lattice, J. Stat. Phys. 68 (1992) 911–923. [41] H. Puhl, Sandpiles on random lattices, Phys. A: Stat. Mech. Appl. 197 (1993) 14–22. [42] H. Herrmann, On the shape of a sandpile, Physics of Dry Granular Media, vol. 350, Springer, Netherlands, 1998, pp. 319–338. [43] C. Kuttler, On the competitive growth of two sand heaps, Math. Method Appl. Sci. 26 (2003) 1435–1449. [44] T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second edition Springer Series in Statistics, 2009. [45] L. Casarsa, T. Arts, Experimental investigation of the aerothermal performance of a high blockage rib-roughened cooling channel, J. Turbomach. T. ASME 127 (2005) 580–588. [46] M. Lohász, P. Rambaud, C. Benocci, Flow features in a fully developed ribbed duct flow as a result of MILES, Flow Turbul. Combust. 77 (2006) 59–76. [47] L. Tran, Turbulent friction and heat transfer calculations of a rectangular duct with repeated ribs on two opposite walls, ASME Winter Annual Meeting, New Orleans, USA, 1993. [48] A. Lai, M. Byrne, A. Goddard, Measured deposition of aerosol particles on a twodimensional ribbed surface in a turbulent duct flow, J. Aerosol Sci. 30 (1999) 1201–1214. [49] E. Matida, K. Nishino, K. Torii, Statistical simulation of particle deposition on the wall from turbulent dispersed pipe flow, Int. J. Heat Fluid Flow 21 (2000) 389–402. [50] T. Barth, G. Lecrivain, U. Hampel, Particle deposition study in a horizontal turbulent duct flow using optical microscopy and particle size spectrometry, J. Aerosol Sci. 60 (2013) 47–54. [51] P. Papavergos, A. Hedley, Particle deposition behaviour from turbulent flows, Chem. Eng. Res. Des. 62 (1984) 275–295. [52] M. Jones, Simple boundary correction for kernel density estimation, Stat. Comput. 3 (1993) 135–146. [53] G. Lumay, F. Boschini, K. Traina, S. Bontempi, J. Remy, R. Cloots, N. Vandewalle, Measuring the flowing properties of powders and grains, Powder Technol. 224 (2012) 19–27
Contents lists available at ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
Three-dimensional simulation of multilayer particle deposition in an obstructed channel flow Gregory Lecrivain a,⁎, Leopold Barry a, Uwe Hampel a,b a b
Institut für Fluiddynamik, Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany Technische Universität Dresden, Institut für Energietechnik, AREVA-Stiftungsprofessur für Bildgebende Messverfahren für die Energie- und Verfahrenstechnik, Germany
a r t i c l e
i n f o
Article history: Received 8 October 2013 Received in revised form 27 February 2014 Accepted 1 March 2014 Available online 10 March 2014 Keywords: Multilayer particle deposition Detached eddy simulation Self-organised criticality Dry granular mechanics Obstructed channel flow
a b s t r a c t A large variety of systems are subject to slow and lengthy processes of solid aerosol particle deposition in turbulent flows. As a result of a long exposure to deposition, a multilayer particle bed eventually forms over time. Notable examples are the formation of multilayer deposits in ventilation ducts, in nuclear reactors or on earth surfaces subject to atmospheric sedimentation. Simulations are of great importance to predict the multilayer deposition of solid aerosol particles. Theoretical models are quite limited since their complexity rapidly increases when the flow becomes turbulent and the surface geometry complex. The present study proposes a new three-dimensional approach to reproduce the growth of a multilayer deposit in a turbulent obstructed channel flow at Reynolds number Re = 10,000. Computational Fluid Dynamics and Computational Granular Dynamics are brought together to simulate 4 h of real deposition. A detached eddy simulation is employed to predict particle deposition while self-organised criticality is employed to reproduce the slow growth of the multilayer deposit. The three dimensional shape of the multilayer deposit matches remarkably well the experimental data. © 2014 Elsevier B.V. All rights reserved.
1. Introduction 1.1. Multilayer deposition of solid aerosol particles in turbulent flows Deposition, i.e. the process of attaching suspended aerosol particles from a fluid in turbulent motion to a wall surface, can either be monolayer or multilayer [1]. Monolayer deposition indicates that all deposited particles only adhere to a wall surface. Multilayer deposition indicates however a more complex particle arrangement: particles sit either on a wall surface or on top of other particles and form a multilayer particle bed, henceforth referred to as the multilayer deposit. A prolonged exposure to particulate matter inevitably leads to the formation of a multilayer deposit. The formation of multilayer deposit spans a wide range of industrial and non-industrial applications, such as particle deposition in ventilation ducts and pipelines [2] and particle deposition in graphite moderated High Temperature Reactors (HTR) [3,4]. There exists very scarce literature on multilayer deposition in turbulent flows. The complexity of the mechanisms involved renders the construction of computer models particularly difficult. Multilayer deposition involves the dynamic combination of two complex mechanisms: the transport of a particulate phase and the slow growth of a granular phase. Similarly to powders, the multilayer deposit exhibits a granular behaviour. One of the striking features of granular materials is the existence of a finite slope, given by the static angle of repose αs, below which ⁎ Corresponding author. Tel .: +49 351 260 3768; fax : +49 351 260 1 3768. E-mail address: [email protected] (G. Lecrivain).
http://dx.doi.org/10.1016/j.powtec.2014.03.019 0032-5910/© 2014 Elsevier B.V. All rights reserved.
the material is at rest [5]. A small stream of particulate matter falling down from a point source onto a flat surface will form a conical heap. The heap maintains under gravity a constant angle αs beyond which the excess material is lost through avalanches [6]. These two phenomena have largely been studied from two different standpoints. The review article of Guha [7] and the book of Pöschel and Schwager [8] testify to the increasing interest in the simulation of particle deposition and in the simulation of granular structures, respectively. An attempt is here made to combine computational fluid dynamics and computational granular dynamics to predict the complex growth of a multilayer deposit. The present work is a full three-dimensional extension of the work of Lecrivain et al. [9]. A fully unsteady CFD simulation and a new threedimensional granular pile model are brought together to simulate the multilayer deposition of solid aerosol particles in a turbulent obstructed channel flow over a fairly long period of time. 1.2. Granular pile models Similarly to sand piles the multilayer deposit obeys the laws of dry granular mechanics. Various methods each having cons and pros have been developed to simulate the heap formation (identical to sand pile formation), whose shape is characterised by the angle of repose αs. An extensive review of the current literature has shown that three approaches have been favoured in the recent years: the discrete element method, the differential approach and cellular automata. The discrete element method is a volume-based method. The heap is made of a large number of grains, which interact with one another
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
Nomenclature Greek αp αs αv Δh Δhth ΔP Δ δp ηd ν ∅ ρ τp ξ
Latin C0 dp e f g H h JS L Li = 1,2,3,4 m p Reb Rep S Ub ub ud uf up w yw
drag correction factor angle of repose airborne particle volume fraction local pile gradient stability threshold pressure drop cube root of cell volume inter-particle spacing deposited fraction kinematic viscosity of the fluid porosity of the multilayer deposit fluid density particle response time non-dimensional coefficient
135
introduced the concept of self-organised criticality. One of the first granular systems, for which cellular automaton was developed, involved the addition of one single grain at a time on a pile. These small perturbations to the system, i.e. the addition of grains, caused avalanche propagations. As illustrated in Fig. 1 a typical one-dimensional pile was regarded as a decreasing sequence of integers hi representing the height of the i-th column. The dynamics of the avalanches was simply described by two operators: adding a grain to the pile and relaxing the slope of the pile wherever the local gradient exceeded a critical stability threshold Δhth [17]. The stability of a column depended on the local gradient Δh = hi − hi + 1. If the local gradient exceeded the threshold Δh N Δhth, the grain in question toppled down in the next lower bin. The development of cellular automata has since been extended to simulate the three-dimensional heap formation on flat support [18]. Unlike the discrete element method and the differential approach, implementation of this surface-based method is relatively simpler. Self organised criticality will therefore be used in the present work to predict multilayer deposition in an obstructed turbulent channel flow. 2. Methods
airborne particle mass concentration particle diameter obstacle height friction coefficient gravity channel height local pile height particle mass flux channel length associated length of the i-th recirculation zone in the streamwise direction particle mass pitch Reynolds number based on bulk velocity Ub and channel height H particle Reynolds number particle-to-fluid density ratio bulk velocity of the fluid build-up velocity deposition velocity resolved component of instantaneous fluid velocity particle velocity kernel function distance to the nearest wall
Miscellaneous X Overbar denoting time or space averaging of arbitrary quantity X X+ Superscript denoting the normalisation of arbitrary quantity X with wall scaling (u⁎ for velocity variables, v/u⁎ for length variables, v/u⁎2 for time variables and u⁎3/v for acceleration variables)
through contact forces [10,11]. A fine tuning of the particle interaction model [12] and the relatively heavy computational cost [8] makes this method however unattractive for large systems. The differential approach is a surface-based method, which indicates that only the height h of the granular interface separating the multilayer deposit from the ambient gas is computed [13]. This analytical approach is however complex since it involves solving a system of partial differential equations for the dynamics of the surface height [14,15]. The emergence of cellular automaton to simulate avalanches on the surface of a granular heap goes back to 1988 when Bak et al. [16]
2.1. Model assumptions The present model offers a significant three dimensional extension of the work of Lecrivain et al. [9]. It therefore inherits all model assumptions previously made. 2.1.1. The flow is very dilute Experimental measurements in the wind tunnel revealed an airborne particle mass concentration of C0 = 2.75 g/m3, a value in line with real experiments previously performed in HTR [3]. With a particle-to-fluid density ratio equalled to S = 1800 and a mean particle diameter of dp = 5.3 μm, it led to a volume fraction of about αv = 1.2 × 10− 6 and a mean inter-particle spacing of δp/dp = 2500. For volume fraction αv b 10−6, Elghobashi [19] stated that the influence of particles on the gas phase and the particle interactions can be neglected.
2.1.2. The deposition process is very slow The multilayer deposit grows at some millimetres per hours while a particle crosses the channel at some metres per second. The build-up velocity ub, i.e. the velocity at which the multilayer deposit grows, is by several orders of magnitude smaller than the particle velocity up. During its life time in the channel, a particle sees a fairly constant shape of the multilayer deposit. The growth of the multilayer deposit is therefore a quasi-static process.
Fig. 1. Addition of one grain on a one-dimensional pile. The grain toppled in the next lower bin until it reaches a stable position.
136
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
2.1.3. Deposition is not impacted by the amorphousness of particles The experimental investigation of Kvasnak and Ahmadi [20] showed that the deposition rate of fibre-like particles in a horizontal duct is very close to that of spherical particles with identical response time. The deposition of amorphous graphite particles can thus be numerically reproduced using equivalent volume spheres without a significant loss of precision. 2.1.4. Particle detachment off the multilayer deposit does not take place Previous investigations performed by references [1] and [21] suggested that the diameter of re-entrained particles was up to an order of magnitude larger than that of the deposited particles, which suggests that particles resuspend off the multilayer deposit in the form of agglomerates rather than as sole entities. Agglomerates, because of their large size, are more easily re-entrained than single fine particles [22,23]. Lazaridis and Drossinos [24] demonstrated that the adhesion force between particles in the multilayer deposit is less than that between a spherical particle and a flat surface. Under very specific flow conditions, the two mechanisms of detachment and redeposition may even occur simultaneously [25,26]. To avoid the detachment of agglomerates, the friction velocity is low at the granular interface. Resuspension experiments performed by Reeks and Hall [22] showed that graphite particles in the size range dp = 6-30 μm initially at rest on flat plate started to reenter the flow beyond the critical friction velocity u∗c = 0.3 m/s. Barth et al. [27] detected particle detachment off a the multilayer deposit for friction velocities greater than u∗c = 0.15 m/s. In the present experiment the time-averaged friction velocity in the cavity did not exceed the value of u∗max = 0.1 m/s. 2.1.5. Depth compaction of the multilayer deposit does not occur Simulation of multilayer particle bed structures performed by Schmidt [28] and recent X-ray measurements of multilayer deposits performed by Barth et al. [27] under atmospheric conditions did not revealed any compaction. 2.2. Channel geometry The present simulation seeks to numerically reproduce the multilayer deposition experiment performed by reference [27]. The flow is wallbounded in a horizontal square duct having a rib-roughened bottom wall surface. The channel height, the rib height and the rib pitch are denoted by H, e and p, respectively. The rib pitch, defined as the distance from one leading edge of the rib to the next leading edge, is set to p = H. The height of the rib obstructs the channel height by e/H = 10 %. In the experimental tests 15 obstacles were placed on the bottom wall. To cut down on the computational cost the flow is however simulated around one single obstacle. The inlet and outlet surfaces coloured in grey in Fig. 2 are enforced a translational periodic boundary condition. 2.3. Numerical procedure A quasi-static approach is suggested to speed-up the three dimensional simulation. It would be extremely time-consuming to continuously simulate hours of deposition in the present obstructed channel flow. Similarly to the two-dimensional multilayer deposition model of Lecrivain et al. [9], the growth of the multilayer deposit is therefore discretised into a series of successive static events. Over the course of each event the geometry of the multilayer deposit remains constant. The simulation process throughout each event is threefold: 1. The transport and deposition of the particles by the gas phase is simulated for a few seconds, 2. The build-up velocity of the multilayer deposit is calculated from the monolayer deposition and linearly extrapolated to reproduce a larger time of exposure, and 3. The deposited mass of graphite is redistributed through avalanches to satisfy the granular properties of the multilayer deposit. The granular interface separating the multilayer deposit from the gas phase acts as a wall boundary in the simulation. The computational domain, in which the gas phase evolves during each event, is bounded
Fig. 2. Schematic of the obstructed channel. A periodic boundary condition is enforced on the two surfaces coloured in grey.
by the channel and the granular interface. Prior to each flow simulation, a Boolean operation is performed to subtract the volume beneath the granular interface from the original channel geometry. In the very first event deposition has not yet taken place. In that case the obstacle is simply subtracted from the channel. 2.4. Simulation of the gas phase Reynolds Averaged Navier-Stokes (RANS) models usually give acceptable results for most engineering applications. RANS simulations are fast, robust but often fail to accurately capture the unsteady flow structures, especially in regions of massive flow separation [29]. Large Eddy Simulation (LES) appears to be a more natural candidate for the unsteady simulation of complex flow separations. The potential accuracy of LES is generally well acknowledged in the literature. The great computational cost resulting from a fine mesh resolution near solid boundaries makes it however unattractive for practical applications. Detached Eddy Simulation (DES) overcomes these wall restrictions by considerably cutting down the cost of a computation. It is a hybrid method in which LES is employed in core and separated flow regions while near the wall RANS model has complete control over the solution. This hybrid concept has been successfully applied to simulate the turbulent flow in obstructed channels [30]. For the present investigation, the DES formulation available in the OpenFOAM library and similar to that proposed by Spalart et al. [31] is used. To accomplish the transition from nearwall RANS-based simulation to LES treatment of the interior flow in one formulation, the nearwall distance yw is replaced by the distance d, defined as e d ¼ minðyw ; C DES ΔÞ
ð1Þ
which act as a RANS model for Δ ≪ yw and as a sub-grid scale model for Δ ≫ yw. The symbol Δ is defined as the cube root of the cell volume and CDES is an empirical constant calibrated to a value of 0.65. More details on the exact implementation can be found in the work of Islam et al. [32]. 2.5. Simulation of the particulate phase The high particle-to-fluid density ratio S and a particle diameter greater than dp N 1 μm implies a negligible effect of the Saffman lift force [33] and of the Brownian force responsible for the molecular
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
137
diffusion [34], respectively. The acceleration of the particle is therefore only influenced by its drag and the gravity. The equation of motion for a point-like particle evolving in the obstructed channel flow reads
given by averaging the positions xi of its N connected sites. Mathematically the space-averaged position of a site is given by
dup u f −up ¼ þ g: dt τp =α p
x¼
ð2Þ
In above equation, up is the particle velocity and uf is the resolved component of instantaneous fluid velocity derived from the DES simulation. The contribution of the sub-grid scale velocity component is here left out since its effect on particle deposition is questionable [35]. The response time τp of a spherical particle is defined as τp ¼
Sd2p 18ν
ð3Þ
where ν is the fluid kinematic viscosity. The drag correction factor αp derived from the semi-empirical formula of Schiller and Naumann is a non-linear function of the particle Reynolds number Rep. Its formulation is given by αp ¼ 1 þ
0:687 0:15Rep :
ð4Þ
The term Rep is the particle Reynolds number based on the relative particle velocity to the fluid velocity and the particle diameter. The solid aerosol particles are spherical in the simulation. The in-house implementation of Eq. (2) is the same as that used in the work of Lecrivain and Hampel [36]. A particle deposits whenever the distance from the particle centre to the nearest wall is smaller than the particle radius. It is assumed that particle rebound does not occur here. In reality multiple rebounds upon collision with a wall surface occur without resulting in deposition, especially at low impact angle [37]. In a gravity-dominated deposition regime, as is the case here, the deposition velocity ud of a particle can be estimated with the formulation of Wood [38]. For a 20 μm graphite particle, the deposition velocity equals ud = g × τp = 0.03 m/s. This deposition velocity is far below the typical impact velocities used for the determination of the rebound coefficient of aerosol particles. Wall et al. [39] experimentally measured the rebound coefficient of 3– 7 μm spheres impacting a wall at velocities ranging from 1 to 20 m/s, for which an impact velocity of 1 m/s did not result in a rebound. Particle rebound, although possible, is unlikely. The low deposition velocity, the low friction velocity and the relatively strong adhesion forces tend to prevent this mechanism. 2.6. Simulation of the granular phase 2.6.1. Lattice construction An unstructured lattice is here used to represent the three dimensional shape of the granular interface. Each site of the lattice stores the height h of the multilayer deposit. At the initial time the shape of the lattice corresponds to that of the roughened channel surface, i.e. h = e for the obstacle else h = 0 for the channel floor. The procedure of Moukarzel and Herrmann [40] is employed to enforce a certain degree of disorder in the lattice, which will avoid an anisotropic distribution of the graphite during the simulation of avalanches. A chessboard-like structured square lattice is first defined. Within each square a random point with uniform distribution, henceforth referred to as a site, is generated. To create the connectivity of each site a constrained twodimensional Delaunay triangulation is performed. The Delaunay triangulation of these sites ensures that the circumcircle associated with each triangle contains no other site in its interior. Boundary sites lying on the bounding lines (walls and periodicity) and on the cliffs of the obstacle have constant coordinates either in the streamwise or spanwise direction. To avoid occurrence of high skewness elements in the lattice, a Laplacian smoothing operator is used. In this simple and yet effective smoothing technique, each internal site is moved to a new position x
N 1X x N i¼1 i
ð5Þ
where N is the number of adjacent sites connected. Collapse of edges smaller than an arbitrarily small value is then performed to finalise the lattice smoothing. The sites forming the obstacles are then shifted vertically by the distance e to form the obstacle. Following this extrusion, the sites located on the two cliffs of the obstacle show a discontinuous behaviour. Each node separating the obstacle from the channel surface now possesses two heights. The above mentioned steps necessary for the construction of smooth unstructured lattice are shown in Fig. 3. The number of sites was deliberately set to a low value for illustration purposes. 2.6.2. Boundary treatment of the lattice 2.6.2.1. Cliff of the obstacle. Each site located on one of the two cliffs of the obstacle must be dealt with cautiously because of a discontinuous behaviour. Each site is decomposed into a lower and upper component. The upper component belongs to the obstacle while the lower component belongs to the channel floor. Upon avalanche the upper component becomes passive and directly conveys the mass it receives to the lower site component. As a result the lower component moves upwards while the upper component remains immobile throughout the mass distribution. As soon as the height of the lower component exceeds that of the upper component by an arbitrarily small distance, the site loses its discontinuity and behaves exactly like any other regular site. 2.6.2.2. Periodicity. During the creation of the lattice, an additional constraint on the Delaunay triangulation is enforced to ensure one-to-one periodic connectivity. The spanwise position of the sites on the inlet line exactly corresponds to that of the sites on the outlet line. During the avalanche process the two paired sites with one-to-one connectivity behave as a single entity: they share the same neighbours and have the same height. 2.6.2.3. Channel walls. In reality, the two side walls of the duct prevent the graphite material from leaving the channel. No numerical treatment is actually needed to prevent the graphite from leaving the lattice, since the material can only be transferred from one site to another. 2.6.3. Avalanche A real function h(x, z) describing the local pile height is associated to each site, where the variables x and z respectively refer to the streamwise and spanwise directions of the lattice. The stability of a site in a given direction solely depends on the nearest neighbours, whose connectivities were created during the Delaunay triangulation. A site is said to be unstable or critical in the i-th direction if the difference in height Δh = h − hi between the site in question and the nearest site in the i-th direction exceeds a threshold Δhth i . The threshold is a function of the distance between the two sites and the angle of repose αs qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi th Δhi ¼ tan ðα s Þ ðx−xi Þ2 þ ðz−zi Þ2 :
ð6Þ
If the site is unstable in the i-th direction, a small amount of mass topples down the adjacent site. Mathematically it means ( th Δhi NΔhi ?
h hi
¼ ¼
h hi
th
þ ξ Δhi th ; − ξ Δhi
ð7Þ
where ξ is a non-dimensional coefficient ranging from 0 to 1. In the original model of Puhl [41] cubical sand grains were toppled (Fig. 1). It
138
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
a
b
c
d
Fig. 3. Lattice construction. (a) Random point generator and constraint lines, (b) constraint Delaunay triangulation, (c) Laplacian smoothing, (d) extrusion of obstacle.
would be however very time-consuming to simulate the toppling of each grain. The above coefficient is therefore introduced to reproduce a small mass transfer. In this case we had ξ = 0.1 which was found to yield a good compromise between accuracy and execution time. To simulate avalanche propagation, the algorithm successively updates the whole lattice until all sites become stable. The order, in which the sites and the nearest neighbours are tested, is fixed during each sweep. Random permutations of the directions did not show any improvements. 2.6.4. Test cases 2.6.4.1. Influence of lattice anisotropy. To emphasise the importance of the disorder in the lattice, the growth of heaps is reproduced on a flat plate. Two point sources, a first source feeding material on a structured lattice and a second feeding on an unstructured lattice, are simulated. The arrangement of the sites in the structured lattice is fairly simple.
Each site possesses nearest neighbourhoods in the four cardinal directions. Fig. 4 shows the two piles obtained. It can be seen that a pyramid is formed on the structured lattice while a conical pile is formed on the unstructured lattice. During the propagation of avalanches on the structured lattice the two directions of the structured lattice are preferred, which is why a pyramid is formed. One can clearly see that the disorder removes the two preferred directions. In reality, the shape the conical heap rarely presents sharp edges. Herrmann [42] showed that piles normally have a logarithmic tail. Therefore the Laplacian smoothing is used in the vertical direction with a very small number of iterations. It smooths the tail of the heap and removes the high frequency oscillations of the granular interface h. 2.6.4.2. Interaction of growing piles. The present granular pile model should satisfy the rules governing the growth of multiple piles [9,43]. For the sake of simplicity, let us assume two point sources feeding
Fig. 4. Effect of lattice structure on heap shape (αs = 45°).
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
139
Fig. 5. Heap growth on an unstructured lattice (αs = 45°). Top: two heaps growing at same velocity, middle: two heaps growing at different velocities, bottom: single heap growing on obstacle.
material with the corresponding fluxes J1 and J2. The two heaps initially grow independently as long as they are far apart and do not interact. However upon collision the final configuration will depend on the flux difference. If the two sources feed material at the same rate (J1 = J2) the configuration reaches an unstable equilibrium and the two heaps grow at the same velocity. If however one source feeds more material (J1 N J2), the faster-growing heap gradually encroaches the slowergrowing heap. One single heap centred at the position of the source J1 emerges from the initial configuration. The faster-growing heap eventually completely covers the slower growing heap. The evolution of this two-heap system is illustrated in Fig. 5. 2.6.4.3. Growth on obstacle. Experimental observations have shown that, when a pile reaches the border of a support, a runway forms and the grains falling from the point source move down that way to the edge and leave the support [18]. Let us assume a point source placed above the obstacle in the channel. A heap initially grows until it meets the edge of the obstacle. Upon contact with the edge, the heap stops growing and the material falls down the cliff. Through a “cascade effect” the
material is transferred to the channel bottom. A second heap starts to grow in the corner. When the summit of the offspring heap reaches the upper cliff, the two heaps join together and then behave as a single entity (Fig. 5). 2.7. Calculation of the build-up velocity The complexity in the present section lies in the transformation of a discrete deposition cloud derived from the particle-laden flow simulation into a continuous surface field of the build-up velocity ub. The build-up velocity, i.e. the velocity at which the multilayer deposit grows is defined as follows ub ¼
JS ; ρð1−ϕÞ
ð8Þ
where JS is the local particle mass flux on a wall surface, i.e. the mass of deposited particles per unit surface and per unit time. A simple method to work out the build-up velocity of each site of the lattice would be to
Fig. 6. Instantaneous and time-averaged flow field at the mid-section.
140
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
centre a circle of a fixed size and sum up the mass of each deposited particle within the circle. A much smoother distribution of the build-up velocity across the lattice can be achieved via a weighting function w (also referred to as kernel function) which assigns a weight to each deposited particle based on its distance from the circle centre [44]. The local particle mass flux of a site is calculated as follows J S ðx; zÞ ¼
N X 1 wðx; zÞmi ; 2 t s πr i¼1
ð9Þ
Table 2 Comparison of recirculation lengths with data from the literature.
DES (present study) Experiment [45] LES [46]
L1 e
L2 e
L3 e
L4 e
1.06 1–1.5 1.62
0.81 0.6–0.9 0.91
0.88 0.2–0.3 0.5
4.12 3.7–3.9 3.45
where mi is the mass of the i-th particle, r is the smoothing radius, ts is the duration of the particle-laden flow simulation and N is the total number of deposited particles on the wall surface. Several types of kernel functions are commonly used in statistics. We choose here to use the standard Gaussian density, with which the weights smoothly decrease towards zero as one moves away from the target point. For deposition on a surface the smoothing function reads
downstream of the obstacle (L3) and to the largest recirculation area downstream of the obstacle (L4), respectively. Table 2 compares the four recirculation lengths with the experimental data of Casarsa and Arts [45] and the LES of Lohász et al. [46]. All but the third associated length match relatively well the experimental data. It can be seen that the present DES offers a great improvement over a conventional RANS simulation, which was found unable to capture the second recirculation zone on top of the obstacle [9]. The time-average pressure drop over one pitch interval ΔP was also computed for validation purposes. The friction factor f is defined as
i 1 h 2 2 wðx; zÞ ¼ exp − 2 ðx−xi Þ þ ðz−zi Þ : 2r
f ¼
The weight function w is typically indexed by the smoothing radius r. A small value of r emphasizes high-frequency features whereas a greater value of r will result in a smoother distribution. Monolayer deposition on vertical walls does take place. However, the wall-induced adhesion forces are not strong enough to allow multilayer deposition on vertical walls. The graphite particles hitting the vertical walls fall off and therefore contribute to the build-up velocity of the sites in contact with the channel side walls and the cliffs (provided the site is discontinuous). For sites connecting the granular interface with a vertical wall, the flux is calculated slightly differently. The disc area associated with the kernel function in Eq. (9) is replaced with 2 ⋅ r ⋅ δ, where δ is the average spacing of two connected sites. To simulate the multilayer deposition, the monolayer deposition is extrapolated in time. The height h(x, z) of the granular interface is obtained from integration of the buildup velocity field as follows dh ¼ ub : dt
ð11Þ
3. Results and discussions 3.1. Turbulent flow field The computational domain, in which the gas phase evolves, has about 0.5 × 106 cells. The typical distance from the first cell centre to + the closest wall is about y+ w = 0.7 around the obstacle, yw = 1.5 on + the channel floor and yw = 3 on the side and top walls. While the resolution of the grid could be finer, the number of grid nodes was kept to a minimal value to ensure a relatively fast simulation of the particle laden flow. Fig. 6 shows the instantaneous and time-averaged flow field through the cross-sectional symmetry plane of the initial channel geometry. To obtain the time averaged flow field, the fluid statistics were averaged in time for 40 flow-through times. Four flow recirculation zones exhibiting a strong rotating nature can be observed near the obstacle. The lengths Li = 1 … 4 are associated to the recirculation areas, upstream of the obstacle (L1), on top of the obstacle (L2),
ð12Þ
where ρ the fluid density. The relative effect of obstacle on the pressure drop can be normalised with the friction factor f0 of smooth circular pipe. The Blasius expression correlates the friction factor with the bulk Reynolds number as follows f0 = 0.046 Re−0.2 . The present DES delivb ered a value of f/f0 = 2.91. It compares well with the value of f/f0 = 3.0 obtained by inserting the present values for obstacle dimension and spacing in the model of Tran [47]. The present ratio is also close to the experimental value of f/f0 = 3.2 taken from the work of Lai et al. [48] at similar Reynolds number. 3.2. Monolayer deposition 3.2.1. Deposition of monodisperse particles The present subsection seeks to compare the deposition rates with existing data available in the literature. While accurate literature data on particle deposition rates in complex geometries is scarce, a few attempts have been made to compute the deposition rates as a function of the particle size in non-obstructed channel flows. The theoretical deposited fraction ηd of fully mixed particles with uniform concentration in a channel bounded by four walls can be estimated using the formulation of Matida et al. [49] L ud ηd ¼ 1− exp −4 ; H Ub 1
0.8
ð13Þ
1
50 dp [μm]
10
Wood Eq. DES
0.6
0.4
0.2
Table 1 Input data for the simulation. Friction velocity Bulk Reynold number Porosity Angle of repose Particle size
H ΔP ; 2 p ρU 2b
ηd
ð10Þ
u* Reb ϕ αs dp
0.06 m/s 10,000 90% 35° 1..20 μm
0 0.001
0.01
0.1
1
τ+p Fig. 7. Deposited fraction against particle response time.
10
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
Discrete deposition cloud
141
Continuous deposition flux
0.8
6
0.6
4
0.4
2
0.2
z/e
8
J/Jmax
1
10
0
0 0
3
x/e
6
9
0
3
x/e
6
9
0.8
6
0.6
4
0.4
2
0.2
y/e
8
0
J/Jmax
1
10
0 5
6
7
8
9
0
1
2
3
4
x/e
5
6
7
8
9
0
1
2
3
4
x/e
Fig. 8. Deposition flux $J_S$ on bottom surface (top) and channel side (bottom).
where ud is the deposition velocity. Wood [38] derived from various experimental measurements a fairly simple expression for the estimation of the deposition velocity þ
þ2
þ þ
ud ¼ Kτ p þ g τ p :
ð14Þ
The first and second terms correspond, respectively, to the eddy diffusion impaction and to the gravitational sedimentation. These two deposition mechanisms where found to be predominant for the present study [50]. The constant K is an empirical value set to K = 4.5 ⋅ 10−4. The deposition is here simulated in an obstructed channel with a length-to-height ratio L/H = 10. To reproduce the length of the channel,
a counter is associated to each particle. Whenever a particle crosses the outlet periodic surface the counter is incremented by one. The simulation is run until all particles have either deposited or have shown a counter greater than 10. Particles with a counter greater then 10 are then removed from the cloud since they are assumed to have left the channel. For each class, about 15,000 particles of constant size are initially placed on a cross-sectional plane normal to the streamwise direction, placed at the position x = H/2. Particles are randomly injected across the planar surface stretching from the obstacle to the channel ceiling. The deposited fraction ηd = Nd/N0 is calculated at the end of the particle-laden flow simulation, where Nd is the number of deposited particles. The particle diameter is varied from 1 to 50 μm. It can be seen from Fig. 7 that the evolution of the deposited fraction compares well with the empirical expression of Wood [38]. The results of the present simulation are very satisfactory, since estimations of the deposition velocity in turbulent flows show large discrepancies often covering as much as one order of magnitude [51]. 3.2.2. Deposition of polydisperse particles
Fig. 9. Simulation of the granular interface after 4 h of deposition.
3.2.2.1. Deposition on channel floor. The size distribution of the particle cloud exactly corresponds to that used in the referenced experiment [27]. The particle equivalent aerodynamic diameter, whose mean value equals 5.3 μm, varies from 1 μm to 20 μm. Similarly to the deposition of monodisperse particles, the flow simulation is initially run during 10 flow-through times to remove artificial flow structures, originating from the initial flow perturbation. About 106 particles are then carried by the flow domain during 40 extra flow-through times. The cloud of deposited particles is used to compute the particle mass
142
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
y/e
0.4
t=1h
Sim t=2h Exp
t=3h
t=4h
0.2
0.0
y/e
0.6
0.4
0.2
0.0 0
3
6
0
3
6
x/e
9
x/e
Fig. 10. Streamwise build-up in the cavity ($z/e = 0.5$) — comparison between simulation and experiment.
flux on the channel bottom, on the obstacle surfaces and on the two channel sides. Fig. 8 illustrates the deposition cloud in the cavity. Not all deposited particles are displayed for illustration purposes. The transformation of the deposited cloud into a smooth surface particle mass flux using Eq. (9) reflects a very inhomogeneous flux in streamwise and spanwise directions across the channel floor. The effect of the recirculation zones on the deposition rate is here of particular importance. The channel bottom area coloured in dark blue and located downstream of the obstacle (x/e b 3) clearly indicates the lower deposition rates due to the largest flow separation. The channel bottom area coloured in dark red and located upstream of the obstacle (x/e N 8) is subject to a higher deposition flux. This is because of the recirculation area associated with the length L1, whose downward flow causes the particles to collide with the channel bottom. As one comes very close to a vertical wall the smoothing function w (Eq. (10)) normally shows a non-physical decreasing behaviour. It is a known limitation of the present smoothing technique because the smoothing circle of radius r goes beyond the lattice boundaries. Although ad-hoc techniques are available to correct the relatively small boundary values [52], their use is not justified here because the flux originating from the deposition on the vertical walls is much greater than the surface flux on the channel floor. For this reason, the surface flux very close to a vertical wall equals the closest flux value, whose value is not influenced by the boundary. 0.3
t=1h
3.2.2.2. Deposition on side wall. The multilayer deposition on vertical surfaces of the obstructions occurs neither in the simulation nor in the experiment. The mass depositing on vertical surfaces therefore solely contributes to the layer build-up on horizontal surfaces through a transfer of material (“cascade effect”). Therefore the deposition flux JS on the side walls solely varies in the streamwise direction. The large turbulence intensity next to the obstacle causes a heavy deposition flux on the vertical sides, while further away from the obstacle the deposition lowers. One may observe that the deposition flux is five times as large as that near the rib in far-off regions. 3.3. Multilayer deposition Simple measurements of the angle of repose for piles of the present graphite material delivered a value of αs = 35°. The experimental value worked out here agrees well with literature data. Typical values of the angle of repose for dry powders vary from 30° to 40° [53]. All input parameters used in the final simulation are summarised in Table 1. In this work, each event reproduces 30 min of real deposition which corresponds to a vertical geometric change of about 1% of the channel height. The discretisation step chosen here reaches a compromise between solution accuracy and computational cost. Each event involves the transport of about 106 particles in a time span of 4 s. Each particle-laden
t=2h
Sim Exp
y/e
0.2 0.1 0.0
t=3h
t=4h
0.4
y/e
0.3 0.2 0.1 0.0
0
1
2
3
z/e
4
0
1
2
3
4
z/e
Fig. 11. Spanwise build-up in the cavity ($x/e = 0.6$) — comparison between simulation and experiment.
5
G. Lecrivain et al. / Powder Technology 258 (2014) 134–143
flow simulation was performed on 2 × 64 processors of AMD 16-Core Opteron 2.3 GHz and required approximately 12 h of CPU time. Robustness of the granular model made possible the simulation of avalanches on a desktop computer. No more that a few hours were necessary for the redistribution of the deposited material across the lattice. All together about two weeks of simulation were necessary to reproduce 4 h of real multilayer deposition. Fig. 9 illustrates the shape of the lattice after 4 h of real deposition. The lattice is coloured by altitude. The height of the two black lines drawn on the granular interface in the streamwise direction (x) and in the spanwise direction (z) is compared with experimental data in Figs. 10 and 11. To observe the build-up progress the experiment was paused each hour and the height of the multilayer deposit in the cavity was measured using a laser scanner. The laser here is used to avoid direct contact that could potentially harm the multilayer deposit. It is found that the simulation matches remarkably well the experimental data. The shape of the multilayer deposit in the cavity and near the vertical walls is well captured. The slope of the granular interface near the vertical walls corresponds exactly to the value of the angle of repose. 4. Conclusions The present study has two novel aspects: 1. A fast and robust threedimensional granular pile model for multilayer growth on complex supports was suggested and 2. the dynamic coupling of the granular pile model was coupled with computational fluid dynamics. The combination of self-organised criticality and detached-eddy-simulation was successfully combined to simulate the multilayer deposition in a turbulent obstructed channel flow. The three dimensional shape of the multilayer deposit matched remarkably well the on-site experimental data. The outcome of this study has a direct impact on various applications such as atmospheric deposition and aerosol particle deposition in ventilation systems and in High Temperature Reactors. Acknowledgements This work was supported by the European Commission under the Grant 249337 of the “Thermal Hydraulics of Innovative Nuclear Systems” research project. References [1] H. Friess, G. Yadigaroglu, Modelling of the resuspension of particle clusters from multilayer aerosol deposits with variable porosity, J. Aerosol Sci. 33 (2002) 883–906. [2] Q. Li, J. Song, C. Li, Y. Wei, J. Chen, Numerical and experimental study of particle deposition on inner wall of 180° bend, Powder Technol. 237 (2013) 241–254. [3] R. Moormann, Fission product transport and source terms in HTRs: experience from AVR pebble bed reactor, Sci. Tech. Nucl. Install. 2008 (2008) 1–14. [4] M. Kissane, A review of radionuclide behaviour in the primary system of a veryhigh-temperature reactor, Nucl. Eng. Des. 239 (2009) 3076–3091. [5] H. Herrmann, Statistical models for granular materials, Phys. A 263 (1999) 51–62. [6] Y. Grasselli, H. Herrmann, On the angles of dry granular heaps, Phys. A 246 (1997) 301–312. [7] A. Guha, Transport and deposition of particles in turbulent and laminar flow, Annu. Rev. Fluid Mech. 40 (2008) 311–341. [8] T. Pöschel, T. Schwager, Computational Granular Dynamics: Models and Algorithms, Springer, 2010. [9] G. Lecrivain, S. Drapeau-Martin, T. Barth, U. Hampel, Numerical simulation of multilayer deposition in an obstructed channel flow, Adv. Powder Technol. 25 (2014) 310–320. [10] T. Shinbrot, M. Zeggio, F. Muzzio, Computational approaches to granular segregation in tumbling blenders, Powder Technol. 116 (2001) 224–231. [11] H. Zhu, Z. Zhou, R. Yang, A. Yu, Discrete particle simulation of particulate systems: theoretical developments, Chem. Eng. Sci. 62 (2007) 3378–3396. [12] J. Baxter, U. Tüzün, J. Burnell, D.M. Heyes, Granular dynamics simulations of twodimensional heap formation, Phys. Rev. E. 55 (1997) 3546–3554. [13] J. Bouchaud, M. Cates, J. Ravi Prakash, S. Edwards, A model for the dynamics of sandpile surfaces, J. Phys. I 4 (1994) 1383–1410. [14] L. Prigozhin, Variational model of sandpile growth, Eur. J. Appl. Math. 7 (1996) 225–235. [15] K. Hadeler, C. Kuttler, Granular matter in a silo, Granul. Matter 3 (2001) 193–197. [16] P. Bak, C. Tang, K. Wiesenfeld, Self organised criticality, Phys. Rev. A 38 (1988) 364–375.
143
[17] M. Aschwanden, Self-organized Criticality in Astrophysics: The Statistics of Nonlinear Processes in the Universe, Springer, 2011. [18] H. Puhl, On the modelling of real sandpiles, Phys. A: Stat. Mech. Appl. 182 (1992) 295–319. [19] S. Elghobashi, On predicting particle-laden turbulent flows, Appl. Sci. Res. 52 (1994) 309–329. [20] W. Kvasnak, G. Ahmadi, Fibrous particle deposition in a turbulent channel flow — an experimental study, Aerosol Sci. Technol. 23 (1995) 641–652. [21] S. Matsusaka, H. Masuda, Particle reentrainment from a fine powder layer in a turbulent air flow, Aerosol Sci. Technol. 24 (1996) 69–84. [22] M. Reeks, D. Hall, Kinetic models for particle resuspension in turbulent flows: theory and measurement, J. Aerosol Sci. 32 (2001) 1–31. [23] G. Ziskind, Particle resuspension from surfaces: revisited and re-evaluated, Rev. Chem. Eng. 22 (2006) 1–123. [24] M. Lazaridis, Y. Drossinos, Multilayer resuspension of small identical particles by turbulent flow, Aerosol Sci. Technol. 28 (1998) 548–560. [25] I. Adhiwidjaja, S. Matsusaka, H. Tanaka, H. Masuda, Simultaneous phenomenon of particle deposition and reentrainment: effects of surface roughness on deposition layer of striped pattern, Aerosol Sci. Technol. 33 (2000) 323–333. [26] S. Matsusaka, W. Theerachaisupakij, H. Yoshida, H. Masuda, Deposition layers formed by a turbulent aerosol flow of micron and sub-micron particles, Powder Technol. 118 (2001) 130–135. [27] T. Barth, M. Reiche, M. Banowski, M. Oppermann, U. Hampel, Experimental investigation of multilayer particle deposition and resuspension between periodic steps in turbulent flows, J. Aerosol Sci. 64 (2013) 111–124. [28] E. Schmidt, Simulation of three-dimensional dust structures via particle trajectory calculations for cake-forming filtration, Powder Technol. 86 (1996) 113–117. [29] H. Lübcke, S. Schmidt, T. Rung, F. Thiele, Comparison of LES and RANS in bluff-body flows, J. Wind. Eng. Ind. Aerodyn. 89 (2001) 1471–1485. [30] A. Viswanathan, D. Tafti, Detached eddy simulation of turbulent flow and heat transfer in a two-pass internal cooling duct, Int. J. Heat Fluid Flow 27 (2006) 1–20. [31] P. Spalart, W. Jou, M. Strelets, S. Allmaras, Comments of feasibility of LES for wings, and on a hybrid RANS/LES approach, International Conference on DNS/LES, Ruston, Louisiana, 1997. [32] M. Islam, F. Decker, E. de Villiers, A.e.a. Jackson, Application of Detached-Eddy Simulation for Automotive Aerodynamics Development, SAE Tech. Paper 2009. (SAE 2009-01-0333). [33] M. Soltani, G. Ahmadi, Direct numerical simulation of particle entrainment in turbulent channel flow, Phys. Fluids 7 (1995) 647–657. [34] M. Horn, H. Schmid, A comprehensive approach in modeling Lagrangian particle deposition in turbulent boundary layers, Powder Technol. 186 (2008) 189–198. [35] M. Salmanzadeh, M. Rahnama, G. Ahmadi, Effect of sub-grid scales on large eddy simulation of particle deposition in a turbulent channel flow, Aerosol Sci. Technol. 44 (2010) 796–806. [36] G. Lecrivain, U. Hampel, Influence of the Lagrangian integral time scale estimation in the near wall region on particle deposition, J. Fluid. Eng.-T. ASME 134 (2012) 1–6. [37] N. Konan, O. Kannengieser, O. Simonin, Stochastic modeling of the multiple rebound effects for particle–rough wall collisions, Int. J. Multiphase Flow 35 (2009) 933–945. [38] N. Wood, A simple method for the calculation of turbulent deposition to smooth and rough surfaces, J. Aerosol Sci. 12 (1981) 275–290. [39] S. Wall, W. John, H. Wang, S. Goren, Measurements of kinetic energy loss for particles impacting surfaces, Aerosol Sci. Technol. 12 (1990) 926–946. [40] C. Moukarzel, H. Herrmann, A vectorizable random lattice, J. Stat. Phys. 68 (1992) 911–923. [41] H. Puhl, Sandpiles on random lattices, Phys. A: Stat. Mech. Appl. 197 (1993) 14–22. [42] H. Herrmann, On the shape of a sandpile, Physics of Dry Granular Media, vol. 350, Springer, Netherlands, 1998, pp. 319–338. [43] C. Kuttler, On the competitive growth of two sand heaps, Math. Method Appl. Sci. 26 (2003) 1435–1449. [44] T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second edition Springer Series in Statistics, 2009. [45] L. Casarsa, T. Arts, Experimental investigation of the aerothermal performance of a high blockage rib-roughened cooling channel, J. Turbomach. T. ASME 127 (2005) 580–588. [46] M. Lohász, P. Rambaud, C. Benocci, Flow features in a fully developed ribbed duct flow as a result of MILES, Flow Turbul. Combust. 77 (2006) 59–76. [47] L. Tran, Turbulent friction and heat transfer calculations of a rectangular duct with repeated ribs on two opposite walls, ASME Winter Annual Meeting, New Orleans, USA, 1993. [48] A. Lai, M. Byrne, A. Goddard, Measured deposition of aerosol particles on a twodimensional ribbed surface in a turbulent duct flow, J. Aerosol Sci. 30 (1999) 1201–1214. [49] E. Matida, K. Nishino, K. Torii, Statistical simulation of particle deposition on the wall from turbulent dispersed pipe flow, Int. J. Heat Fluid Flow 21 (2000) 389–402. [50] T. Barth, G. Lecrivain, U. Hampel, Particle deposition study in a horizontal turbulent duct flow using optical microscopy and particle size spectrometry, J. Aerosol Sci. 60 (2013) 47–54. [51] P. Papavergos, A. Hedley, Particle deposition behaviour from turbulent flows, Chem. Eng. Res. Des. 62 (1984) 275–295. [52] M. Jones, Simple boundary correction for kernel density estimation, Stat. Comput. 3 (1993) 135–146. [53] G. Lumay, F. Boschini, K. Traina, S. Bontempi, J. Remy, R. Cloots, N. Vandewalle, Measuring the flowing properties of powders and grains, Powder Technol. 224 (2012) 19–27