Effect of gas-particle interaction on roller spreading process in additive manufacturing

Article history: Received 10 April 2020 Received in revised form 27 May 2020 Accepted 31 May 2020 Available online 06 June 2020 Powder spreading in Additive Manufacturing (AM) has been analysed extensively by the Discrete Element Method, but without considering the presence of ambient gas. For fine particles, as commonly used in AM, the gas drag could affect the quality of spread layer. Here, we consider the dynamics of powder spreading by a roller for a gas-atomised metal powder and analyse the combined effects of gas-particle interaction and interparticle adhesion on the particle flow in the heap and spread layer uniformity. In the presence of gas, the convection and circulation of particles within the heap are slowed down, and the heap repose angle becomes steeper. The amount of particles spread on the base is reduced, as compared to the case in which gas drag is not considered, but surprisingly particles with larger interparticle adhesion form a more uniform spread layer with larger total particle volume when gas drag is taken into account. © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
Additive Manufacturing (AM) based on powder spreading has received great attention in recent years in wide ranging applications [1][2][3][4][5][6]. In a number of methods for additive manufacturing, fine powder is spread as a thin layer for building up an object layer by layer. Spreading fine and cohesive powders poses great challenges for producing a uniform layer as the interparticle attractive forces and ambient gas drag have adverse effects on spreading. They could cause transient jamming, defects in the spread layer, and even formation of empty patches [7], thereby adversely affecting final product quality [8,9]. The spread layer should be very thin for selective melting by a radiative energy beam or other bonding processes, therefore requiring a narrow gap between the spreader and layering base, typically a few multiples of particle diameter, and the spreading speed should be high to promote a rapid manufacturing. However, lack of sufficient understanding of the dynamics of this kind of near-boundary flow under a high shear strain rate hinders further advancement in this technology and introduction of new materials for high quality structures.
Recently, the Discrete Element Method (DEM) has been used extensively to describe the mechanical behaviour of the powder spreading system in AM [7,[10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Parteli and Pöschel [10] and Haeri et al. [11] studied the effect of spreading conditions on the bed quality in the roller spreading system, and showed that a high translational velocity of the spreader could lead to low bed quality. Herbold et al. [18], Mindt et al. [19], Zhang et al. [20], Lee et al. [21], Fouda and Bayly [22], Desai and Higgs [24] also identified the spreading conditions affecting the layer quality, but the particles in their DEM simulations were spheres without taking account of interparticle adhesion, which are not representative of real powders used in AM. Also the use of mono size spheres in the work of Zhang et al. [20] and Fouda and Bayly [22] brings about a crystalline order to the structure, which is unrealistic, as slip occurs more readily on selective planes similar to cleavage planes of crystals. Haeri [13] identified the optimum blade tip shape to produce a spread layer with volume fraction and surface roughness comparable to a roller under the actual operation conditions. Chen et al. [12,23], Meier et al. [15,16] and Han et al. [17] studied the effect of interparticle adhesion on the bed quality in the blade spreading system, where spherical particles with tuned mechanical properties were used for DEM simulation in the blade spreading process, and found that large particle adhesion could lead to poor quality of spread layer. Nan et al. [7] characterised all relevant physical properties and interaction parameters of 316 L steel particles used in DEM simulations, and showed for the first time that transient jamming occurred in narrow spreader gaps and was responsible for the formation of empty patches over the work surface in the blade spreading process. This was later validated by experiment work on the same powder system [26]. Nan et al. [14] examined the shear band of the near-boundary particle flow in the spreading process, and showed that the velocity of the particles in front of the blade could well be described by Gauss error function. Nan et al. [25] recently identified the particle convection/circulation of particles within the heap in the roller spreading system, and found that the segregation in the spread layer was related to the jamming events of the near-boundary flow between the rough roller and base. However, in additive manufacturing processes, fine powders in the micrometre size range are commonly used, which implies that they are not only affected by interparticle adhesion but also by fluid drag [27,28]. Nevertheless, the effect of gas-particle interaction on the spreading process has not so far been analysed in detail. This is especially the case for roller spreading system, where large roller rotational speeds exert a notable drag on the particles and ambient gas, thereby affecting the spreading process. We therefore make a first attempt in this work to analyse the combined effects of gas particle interactions and interparticle adhesion, which give rise to bulk cohesion.
The effect of gas-particle interaction on the particle flow in the roller spreading process is analysed by numerical simulations by the combined approach of Discrete Element Method (DEM) and Computational Fluid Dynamics (CFD) using the full four-way coupling method [29][30][31][32]. A heap of gas-atomised stainless 316 L steel particles in the presence of argon gas is simulated and subjected to the translational and anticlockwise rotation of a roller, as it moves from left to right. The particles are spread onto a rough base through a narrow gap between the roller and base. The dynamics of the particles in the heap and the quality of spread layer are analysed. This provides a step forward in our understanding of the role of ambient gas interaction and interparticle adhesion on particle flow in the spreading process, which is a key factor affecting the production quality.

Methods
The spreading process of particles is simulated by the coupled DEM-CFD approach, in which the particle and fluid dynamics interact with each other, i.e. the so called 'four-way coupling' [29][30][31][32]. The DEM-CFD numerical simulation platform used here is based on EDEM™ software provided by DEM Solutions, Edinburgh, UK and Fluent software provided by ANSYS.

Discrete element method
The particles are modelled as discrete entities and their translational and rotational motions are tracked individually by solving Newton's laws of motion [33][34][35]: where m i , I i , v i and ω i are the mass, moment of inertia, translational velocity and angular velocity, respectively; f pf,i is the fluid-particle interaction force on the particle; F c,i is the contact force, originating from its interaction with neighbouring particles or walls; M c,i is the contact torque, arising from the tangential and normal contact forces; R i is the rotation matrix from the global to the local coordinate system in which the calculation of the rotation expressed by Eq. (2) is accomplished. As introduced by Favier et al. [36], the non-spherical particles used in the spreading process are described by the overlapping multi-sphere model. Thus, the interactions between any two non-spherical particles can be simplified as that of spherical particles. In this work, the elastic contact force is described by Hertz-Mindlin contact model [34], and the adhesive interaction is accounted for by JKR model [37]: where Г is the interfacial surface energy; E* is the equivalent Young's modulus; R* is the equivalent radius; a is the contact radius. Based on Ferrari's solution [38], the contact radius a could be analytically calculated from the normal overlap α through this equation: In the unloading process, the normal contact force F n is not zero when the normal overlap α is negative, as further work is required to separate the cohesive contact. More features and further information of the contact model are given by Thornton [34] and Pasha et al. [39], which are not shown here for brevity.

Continuum model of fluid phase
The motion of the gas phase is described by the fluid continuity and Navier-Stokes equations based on the local mean variables over a computational cell. The governing equations for the incompressible fluid accounting the effects of fluid-particle interaction are given as [40]: where ρ f is the density of fluid; ε f is the volume fraction of fluid in the computational cell with the volume of ΔV; u and p are respectively the velocity and pressure of fluid; τ is the shear stress tensor of fluid; F pf is the volumetric fluid-particle interaction force acting on the fluid phase, given as: where k c is the number of particles in the fluid cell; f d,i is the drag force on the particle.

Fluid-particle interaction
The dominant forces acting on a particle due to the fluid flow are the pressure gradient force and drag force, given as [40]: where V p,i is the volume of an individual particle. For simplicity, the subscript i is omitted in following description. The drag force f d is calculated from the local particle volume fraction and relative velocity between the fluid and particle, given as: where β is the fluid-particle exchange coefficient, given as [41]: where d V is the diameter of the sphere which has the same volume as the considered particle, and C d is the drag coefficient of particles, given as [41]: where Re p is the particle Reynolds number.

Simulation conditions
The particle properties used in the simulation are based on gasatomised 316 L stainless steel particles, provided by Sandvik Osprey Ltd., Neath, UK, and have previously been characterised by Nan et al. [7]. They are summarised in Table 1 for ease of reference. The particles have a size distribution in the range 15-55 μm, obtained by image analysis of their projected area, obtained by scanning electron microscopy. Their characteristic measures of the particle size distribution, i.e. number-based equivalent-circle diameters D 10 , D 50 and D 90 , are 20 μm, 32 μm and 45 μm, respectively. They are classified into four main size classes based on their equivalent-circle diameter of the projected area. The number frequencies for the size classes of 15-25 μm, 25-35 μm, 35-45 μm and 45-55 μm are 29.6%, 40.8%, 23.9% and 5.7%, respectively. For each size class, several types of particles are randomly selected and their shapes are reconstructed by the clumped sphere method, with a total of 24 types of particle shapes used in the simulation, as shown in Fig. 1. To make the computational time practical, Young's modulus in the simulation is 100 times smaller than the one measured in the experiment, i.e. E sim = 2.1 GPa, and the surface energy is correspondingly scaled down in the simulation, i.e. Г sim = 1.4 mJ/m 2 , according to the criterion proposed by Behjani et al. [42], Haervig et al. [43] and Washino et al. [44]: where Г and E are the experimental values shown in Table 1.
The simulation system comprises a spreading roller with diameter of 4 mm, and a rectangular base with length of 400D 90 , as shown in Fig. 2. The front and rear boundaries (i.e. in the Y direction) of the simulation domain are treated as periodic boundaries for particle flow. Both the roller and base have the same width as the simulation domain in the Y direction, i.e. 10D 90 . To mitigate the bulk sliding of the particles, the base and roller are made up of clumped cylinders with axes along the Y direction. The initial particle bed is prepared by using the poured packing method, where approximately 21,000 particles are generated. The roller spreader is then placed at a specified position, forming a vertical gap of δ (i.e. 2.5D 90 = 112.5 μm) between the roller and base. The previous work showed that a uniform spread layer could be obtained by this gap height [25]. As the spreading process begins, the roller moves along the X direction (i.e. from left to right) with a constant translational speed of U (i.e. 0.08 m/s) and rotates anti-clockwise with a constant angular speed of ω (i.e. 60 rad/s), giving a tip speed of 0.2 m/s, by which the particles are spread onto the rough base. As the roller starts to move, its traction with the particle heap drags the particles (Re p < 5) upward and with them the surrounding gas, by which a laminar gas flow around the roller is induced, as shown in Fig. 2(c). The effect of roller movement on the gas flow is considered by using the method of dynamic mesh. For DEM-CFD simulations, the gas of argon with density of 9.7 kg/m 3 and viscosity of 2.125 × 10 −5 Pa·s is used. The time steps in DEM and CFD are 10 −8 s and 10 −6 s, respectively, and the mean grid size in CFD is 2.5D 90 . In the powder spreading process in some AM machines, there is an inert gas flow in the lateral direction [45,46]. However, its effect on the particles being laid on the base is small to avoid entraining the particles away from the process zone. Furthermore, it does not intervene with the gas flow that is induced in the heap in front of the roller by the roller rotation. Thus, the background lateral gas flow is omitted in this work and the DEM-CFD simulation starts from a quiescent flow.
To explore the effect of extent of gas-particle interaction on the heap and spread layer, the spreading process is simulated with increased fluid drag on the particles, corresponding to smaller particle sizes, d s (1/1.5, 1/2, and 1/3 times of the original size class, i.e. D = 15-55 μm shown in Fig. 1), whilst keeping the actual particle size constant. So, for the standard PSD under consideration, i.e. D = 15-55 μm, the gasparticle interaction force is based on Eqs. (8)(9)(10). To account for larger fluid drag on smaller particles, i.e. for d s = 1/1.5D, 1/2D and 1/3D, the gas-particle interaction force is scaled up by the scale factor SF according to following criterion [47,48]:  (45-55 μm). The number frequency of each particle shape is 5% for the first 19 boxes (i.e. 15-45 μm) and 1% for the last 5 boxes (i.e. 45-55 μm), resulting in number frequencies for four size classes of 30%, 40%, 25% and 5%, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 Physical and mechanical properties of single particle characterised by Nan et al. [6].
where F f DEM-CFD = SF×F f,s is the scaled-up gas-particle interaction force on the particles in simulation; F f,s is the gas-particle interaction force calculated based on Eqs. (8)-(10) and d s ; m s and m are the mass of the particles with d s and D, respectively. This procedure speeds up the simulation, as the time-step is proportional to the smallest size of sphere elements used to describe the shape of particles shown in Fig.1.
It should be noted that for the particles with d s , their surface energy and other physical properties (e.g. friction coefficient, Young's Modulus) and spreading conditions (i.e. spreading speed and gap) in the simulation are kept the same as the one with D (i.e. 15-55 μm), which are shown in Table 1 and Fig. 2. Thus, the ratio of adhesive force to particle weight of the smaller particles (d s ) is the same as the one for particles of size D. As the gas drag on the particles is increased to correspond to smaller particle sizes, the ratio of fluid drag to particle weight increases, thus affecting the dynamics of spreading profoundly. Therefore, we first examine the effect of gas-particle interaction, where the surface energy is kept the same in all the simulations, and then investigate the coupled effect of gas-particle and adhesive interactions to account for their influence on small particles, which are more affected by both gas drag and adhesion. This is done for a case with a larger value of surface energy.
Archimedes number (Ar), terminal velocity (u t ) and Stokes number (Stk) of a single particle with the diameter of d s,90 are given in Table 2. u t is calculated based on the drag coefficient in Eq. (10). The ratio of roller tip speed (i.e. u tip = 0.2 m/s) to the terminal velocity (u t ) is used as the non-dimensional number to quantify the effects of roller-particle-gas interactions on particle dynamics. In the case here, Stk describes the drag of the particles on the gas, as they are moved by the roller. Stk may be used interchangeably with u t here, because the particles are small so that the dependence of both u t and Stk on particle size scales with d s 2 . The traction of the rough roller on the heap drags the particles, and with them the gas. For large particles and hence large Stk values, the particles cannot drag the gas up with them as effectively as the small ones do, corresponding to small Stk values. Nevertheless, the effect of the particles dragging the gas with them can be equally expressed by either of them.

Particle flow within the heap
The snapshots of the heap with particles coloured, based on their velocity magnitude, are shown in Fig. 3, and the corresponding probability distribution of particle velocity is shown in Fig. 4. For the case in which the gas drag is not considered, hereinafter referred to as 'vacuum' (i.e. Fig. 3(a)), the heap is mainly classified into four velocity blocks, i.e. A, B, C and D. They are all in cascading style, where the boundaries between them are in parabolic form. This feature is attributed to the rotational motion of the roller, which has been detailed discussed in Nan   et al. [25]. Due to the action of roller rotation, some particles are cast far away from the heap and the topmost part of the heap does not touch the roller surface. Correspondingly, as shown in Fig. 4, the probability distribution of the particles with velocity of 0.02-0.12 m/s is almost flat with the total probability of 0.7. As the particle shear band caused by roller shearing is narrow (i.e. around 3-5 particle diameter), and due to the base being rough, the probability of particle velocity close to maximum tip speed of roller (i.e. 0.2 m/s) is very small. This is consistent with the observations shown in Fig. 3. As the gas-particle interaction in the spreading process is considered for the cases in the presence of gas (i.e. u tip /u t = 0.7 and 2.2 shown in Fig. 3), the heap shows several distinct differences from the case of vacuum: (i) the heap is mainly classified into blocks B and C, i.e. 0.04-0.08 m/s and 0.08-0.12 m/s; (ii) the probability distributions of particle velocity become narrower with larger peaks, as shown in Fig. 4; (iii) the slope of the heap at the right hand side approaches closer to the roller and has a larger tilt angle (i.e. tending to vertical) for finer particles (i.e. at larger u tip /u t ); (iv) interestingly, there are almost no particles cast away from the heap, consistent with the previous statement and the topmost point of the heap is attached to the roller surface. This implies that the particles close to the slope of heap at the right hand side are pushed back on to the heap. The intriguing observation that the repose angle becomes steeper when the gas drag on the particles is increased and also steeper than the case of vacuum is clearly indicative of the influence of gas drag during the spreading process, which cannot be ignored for fine particles (i.e. at large u tip /u t ) and requires further probing. The drag force on the particles within the heap is in the opposite direction of the roller spreading, due to larger particle velocity than the induced gas flow (i.e. u gas −u p < 0).
To illustrate the influence of gas-particle and particle-particle/roller interactions, the region with significant shear straining within the heap is tracked, i.e. the arc close to the roller, which is in orange colour in Fig. 5(b). It is divided into 9 bins with an arc angle of 5°. Fig. 5 (a) shows the distribution of gas-particle interaction force F f and particle contact force F c , which are normalised by the weight of particle and then averaged for all particles over time in the considered bin. For all the cases, the contact force F c is larger than the maximum adhesion force F c,JKR = 1.5πГD 90 , while the gas-particle interaction force F f is less than 1.5πГD 90 , indicating that the motion of particles around roller is very much dominated by the contact interaction of particle-particle/ roller. The roller lifts the particles in a narrow band and drags the gas with them. Compared to the case of vacuum and with the increase of u tip /u t , the contact force F c in the presence of gas is larger. This is also the case with the increase of contact number within the band and with the roller, which is not shown here for simplicity. However, the extent of increase of F c is much less than that of gas-particle interaction force F f , resulting in the decrease of the ratio of F c to F f , as shown in Fig. 5(b). This suggests that finer particles (i.e. corresponding to larger u tip /u t ) are more affected by the gas drag, as intuitively expected, and also as observed in Figs. 3-4. To further illustrate the effect of gas-particle interaction on the particle flow in the heap, the spatial positions of the particles in four cuboid cells with size 2D 90 in both X and Z directions and full depth in Y direction (i.e. 10D 90 ), are tracked from t = 0.04 s (i.e. the position shown in Fig. 6) to t = 0.15 s (i.e. the end of spreading process) for the two cases of spreading in vacuum (a1 to l1) and in the presence of argon gas with u tip /u t = 2.2 (a2 to l2). The tracked particles in the four cells are in red, green, black and orange colours, respectively, while other particles are transparent. For both cases of vacuum and considering gas-particle interaction in the presence of ambient gas, particles in cell 1 (i.e. red colour) are spread onto the rough base quickly, but particles in other  cells show extensive convection and circulation within the heap, as they lifted up by the roller and cascade down on the heap slope. However, comparing the two cases of vacuum and ambient gas, distinct differences prevail, as the convection and circulation of particles in the latter case are slowed down. For example, the particles in cell 2 (i.e. green colour) could be completely spread onto the base in the case of vacuum, but in the case of u tip /u t = 2.2, they are still cascading down on the slope of heap at the right hand at the end of spreading process. To better illustrate this feature, the averaged trajectories of the particles in cell 3 (i.e. black colour) are recorded, as shown in Fig. 7, where the abscissa is the relative position of the particle centre x with respect to the roller centre x c . For the case of vacuum, the particles could finish a whole cycle of convection/circulation. As the gas-particle interaction is considered  and with the increase of u tip /u t , the trajectories become shorter, resulting in less entrainment of particles by roller movement. For example, the particles in the case of u tip /u t = 4.6 only reach the top of the heap at the end of spreading process, implying that the gas drag slows the particles being sheared upwards, as intuitively expected.

Particle spread layer
When the spreading process is finished, a thin layer of particles is formed on the base, as shown in Fig. 8. The total volume of particles V layer within the spread layer is normalised as: where V p is the volume of an individual particle; L is the length of the spread layer (i.e. in the X direction); W is the width of the spread layer (i.e. in the Y direction). For the case of vacuum, as shown in Fig. 8(a), the base is almost fully covered by particles, and the total particle volume of particle V layer is 0.61, as shown in Fig. 9. As the gas-particle interaction is considered and with the increase of u tip /u t , the layer becomes patchier, as shown in Fig. 8, resulting in a decrease of V layer shown in Fig. 9. For example, V layer in the case of u tip /u t = 0.7 is a little lower than that of vacuum, and it decreases to almost half of that of vacuum in the case of u tip / u t = 4.6, implying finer particles (i.e. at larger u tip /u t ) are severely affected by gas drag. In the case of vacuum, the convection and circulation of particles in the heap are strong due to the shearing effect of the roller on the heap, and this makes the particles more easily entrained by the roller movement, as previously analysed by Nan et al. [25]. In contrast, in the spreading process in the presence of gas, as considered in this work, the particles are less entrained by the roller movement (i.e. slow particle convection/circulation) as the fluid drag is increased, corresponding to smaller particle sizes. As discussed in Section 3.1, the presence of the gas appears to produce a heap having a steep repose angle, and with no particles getting scattered around the heap slope, somewhat reminiscent of flow of partially-wet sand under negative pore pressure. This implies a lower spreadability, which explains the visual Fig. 9. Variations of the total particle volume of spread layer in vacuum and the presence of gas. Fig. 10. Ratio of shear to normal stresses with inertial number in the gap region (x = 0-6D 90 ) for the cases in vacuum and the presence of gas.  differences of Fig. 8(a-d), i.e. with increasing the fluid drag on the particles (expressed by increasing u tip /u t ), the spread layer becomes less uniform and patchier. Particle flow through the gap is responsible for forming the spread layer on the base and is examined in terms of particle stress and velocity in the gap region. The former is related to the rheological behaviour of the particles in the gap, whilst the latter describes particles leaving the gap. Here, we consider the gap region x = 0-6D 90 , where x = 0 is the x position of roller centre, and x = 6D 90 is in front of roller vertical centre line in the spreading direction. The stress tensor of the particle flow within the gap region is given as: where V is the volume of the gap region; m p is the mass of particle p; δv i and δv j are the fluctuation velocities of particle p; f ij is the contact force at contact c and r ij is the corresponding branch vector between mass centre of particle i and that of particle j. Based on the stress tensor, three principal stresses could be calculated: major one σ 1 , intermediate one σ 2 and minor one σ 3 . The normal stress σ and shear stress τ are then given as: Similar to the rheological analysis of dense particle flow, the inertial number I [49] is used to describe the shear stress regime: The ratio of shear to normal stresses τ/σ of particle flow in the gap region is shown in Fig. 10, indicating resistance to flow due to bulk friction and cohesion. As the gas-particle interaction is considered and with the increased fluid drag on the particles (expressed by the increase of u tip /u t ), inertial number decreases, which is due to the increase of normal stress (the particle size is kept constant), and τ/σ creeps down as the fluid drag starts to affect the particle motion in the gap region. τ/σ decreases linearly with the decrease of inertial number, indicating that the particle flow in the gap region is in immediate regime.
The particle velocity in the spreading direction in the gap region is shown in Fig. 11. As the gas-particle interaction is taken into account and with the increase of u tip /u t , particle velocity v x increases. This is attributed to the positive gas drag, i.e. in the same direction as roller spreading, as in the gap region (u gas -v x ) > 0 where v x is retarded by the frictional traction of the rough base. Thus, the particles are more dragged by the induced gas flow to move with the roller in this region, resulting in smaller leaving velocity (i.e. U-v x ) if taking the coordinate system on roller centre. The particles in the gap region follow the moving roller more easily, and resulting in the decrease of total particle volume of spread layer.

Combined effect of particle adhesion and gas-particle interaction
To explore the combined effect of particle adhesive interaction and gas-particle interaction, the spreading processes of particles with surface energy of Г sim = 11.2 mJ/m 2 in vacuum and also in the presence of gas for the case of u tip /u t = 2.2 are simulated with the same spreading Fig. 13. Normalised contact force F c and gas-particle interaction force F f on the particles around the roller in the cases of vacuum (solid symbol) and u tip /u t = 2.2 (other symbols).  The total particle volume of layer, V layer , for the two values of surface energy is shown in Fig. 12. For the case of vacuum, as Г sim is increased from 1.4 mJ/m 2 to 11.2 mJ/m 2 , V layer decreases from 0.61 to 0.57. In contrast, in the case of the gas-particle interaction the normalised particle volume that is spread actually increases as the surface energy is increased, although the presence of gas reduces the spread volume. Considering the case of u tip /u t = 2.2, V layer is smaller than that of vacuum for both surface energy values. However, the extent of the decrease in the case of Г sim = 11.2 mJ/m 2 is less than that of the case of Г sim = 1.4 mJ/m 2 , indicating that the interparticle adhesion reduces the negative impact of gas-particle interaction on the total particle volume of spread layer. Meanwhile, in the presence of gas, V layer of Г sim = 11.2 mJ/m 2 is larger than that of Г sim = 1.4 mJ/m 2 , suggesting that surprisingly more adhesive particles could have a better spreadability (in terms of total particle volume of spread layer) if the gas-particle interaction is considered. This is consistent with the results of the work of Ghadiri et al. [50], which contrasted the mechanisms of particle flowability and spreadability. The latter is related to the particle flow in narrow gaps with a thickness of a few particle diameter, whilst in the case of the former, the information obtained by traditional flowability testing methods is not easily extendable to describe spreadability. To depict the underlying mechanisms, the combined effects of interparticle adhesion and gas-particle interaction on the particle flow within the heap and gap are explored, as shown in Figs. 13-15. For the cases in vacuum, as Г sim is increased from 1.4 mJ/m 2 to 11.2 mJ/m 2 , the contact force increases (see Fig. 13), resulting in larger impact of the roller movement on the heap and more entrainment of particles by the roller movement. Meanwhile, the ratio of shear stress to normal stress increases, indicating larger resistance to flow through the gap region (see Fig. 14). Therefore, with the increase of particle surface energy, less particles could enter into the gap region and they have more difficulty to flow through the gap region, resulting in the decrease of the total particle volume of spread layer.
For addressing the effect of gas-particle interaction, the case of u tip /u t = 2.2 is considered. As Г sim is increased from 1.4 mJ/m 2 to 11.2 mJ/m 2 , the contact force increases substantially, while the gas-particle interaction force changes very little, resulting in a large ratio of F c to F f , as shown in Fig. 13. Thus, the effect of the gas-particle interaction on the heap is lessened by increasing interparticle adhesion. Meanwhile, the ratio of shear stress to normal stress is larger for Г sim = 11.2 mJ/m 2 as compared to the case of Г sim = 1.4 mJ/m 2 , as shown in Fig. 14. Thus, interparticle adhesion has a contributory effect to increased stress ratio, but to a lesser effect as compared to the case of vacuum. Interestingly, the combined effect of gas-particle interaction and adhesion is more notable on the particle velocity in the gap region, as shown in Fig. 15. As Г sim is increased from 1.4 mJ/m 2 to 11.2 mJ/m 2 , particle velocity in the spreading direction in the presence of gas is lower for the latter case (i.e. 11.2 mJ/m 2 ), indicating a better spreading. However, gas drag has a deleterious effect on spreading as the velocity in the spreading direction is larger than that of the case of vacuum. This is supported by the spread layer particles in Fig. 8.

Discussions
The presence of gas drag on the particles in roller spreading has a remarkable influence on the physical and mechanical features of the heap and spreading dynamics. The steeper dynamic repose angle of the heap due to gas drag is indicative of rheological changes in the heap, deserving further analysis of the distribution of the apparent viscosity and bulk friction angle within the heap, considering that the shear straining is not uniform. The convection and circulation of particles within the heap are  clearly slowed down in the presence of gas, and the amount of particle spread on the base is reduced. But surprisingly, increasing the bulk cohesion (i.e. increasing the surface energy of the particles) counteracts the adverse effect of gas drag and provides a more uniform spread layer with larger total particle volume. It appears that the combined effect of gas drag and interparticle adhesion imparts a bulk behaviour that is reminiscent of flow of partially-wet sand, which is under negative pore pressure, based on the appearance of the heap surface viewed by video playback of the simulation. Further analysis of the heap dynamics to explore the stress state therein as a function of the roller rotational and translations speeds will be informative. Some interesting macroscopic trends are depicted in Figs. 16 and 17. Clearly the resistance against the translational motion of the roller increases with increased gas drag, as shown in Fig. 16. Also the downward force and shear traction on the base, normalised by the heap weight, increase notably as the gas drag is increased, shown in Fig. 17(a & b). Interestingly, the downward force increases by 35% over the weight of the heap for the case in which the gas drag is the greatest, implying the heap is pressed down, despite the lifting action of the counter-rotating roller. This in turn increases the shear traction on the base, as shown in Fig. 17(b). The above trends are therefore indicative of the influence of gas drag during the spreading process, which cannot be ignored for fine particles (i.e. at large u tip /u t ).

Conclusions
The particle spreading process with a roller spreader has been analysed by DEM-CFD simulations, taking account of particle-fluid drag and using realistic physical and mechanical properties of particles as measured for single particles in the previous work. The effect of gas-particle interaction on the spreading has been analysed and quantified in terms of the particle flow in the heap and the quality of spread layer. The combined effect of interparticle adhesion and gas-particle interaction is also analysed. The main results are summarised as follows: 1) Taking account of the gas drag has notable influence on the particle behaviour in the heap. The convection and circulation of particles within the heap are slowed down, and the heap repose angle becomes steeper. The gas drag has a deleterious effect on the quality of the spread layer. 2) Particle entrainment in the heap by the drag exerted by the counter rotation of the roller, the stress ratio in the gap, and the velocity of particles leaving the gap are influenced by the coupled effect of particle adhesion and gas-particle interaction. 3) With the increase of gas drag (represented by u tip /u t ), particles in the gap more easily follow the moving roller without spreading on the base, resulting in the decrease of total particle volume of spread layer.

4) Increasing the bulk cohesion (by increasing interparticle adhesion)
counteracts the adverse effect of gas drag and provides a more uniform spread layer with larger total particle volume.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.