CET Volume 86 DOI: 10.3303/CET2186136 Paper Received: 26 October 2020; Revised: 10 February 2021; Accepted: 8 April 2021 Please cite this article as: Marchelli F., Di Felice R., 2021, On the Influence of Contact Models on Friction Forces in Discrete Element Method Simulations, Chemical Engineering Transactions, 86, 811-816 DOI:10.3303/CET2186136 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 On the Influence of Contact Models on Friction Forces in Discrete Element Method Simulations Filippo Marchelli*, Renzo di Felice Department of Civil, Chemical and Environmental Engineering, University of Genova, Via all’Opera Pia 15a, 16145 Genova, Italy filippo.marchelli@edu.unige.it The discrete element method (DEM) is an important tool to simulate granular systems with high accuracy. Depending on the application, it is often unclear which model is more appropriate for the calculation of collision forces: the Hertzian approach is generally considered more accurate, but it makes the simulation significantly slower than the Hookean one. In this work, these two approaches are compared in two different situations: the stress distribution of static particles in a cylindrical column (DEM) and the onset of water fluidisation for a completely segregated mixture (CFD-DEM). In both cases, particle contact forces are of great relevance to determine the output. It is found that in the first case the Hookean approach does not produce the expected asymptotical stress trend and does not even respond satisfyingly to friction mobilisation. Conversely, in the fluidisation simulations the results are virtually identical, pointing out that the more complex Hertzian approach may be unnecessary in that case. 1. Introduction Granular materials are of the utmost importance in both natural and industrial processes. In recent decades, several methodologies have been proposed to simulate them at various degrees of accuracy and detail. The popularity of these methodologies has been aided by the great enhancement in computational capacity of modern computers. One of the most common is the Discrete Element Method (DEM), which was proposed by Cundall and Strack (1979). The method relies on solving the Newtonian equations of motion for each involved particle, which is achieved through numerical methods. The higher the number of the simulated particles is, the longer the simulation will need to reach the results. Typically, no more than a few hundred thousand particles are simulated, although there have been cases with higher numbers. Despite its limitations, the method is nowadays popular (Guo and Curtis, 2015), also for its coupling with computational fluid dynamics (CFD) to simulate multiphase gas- or liquid-solid flows (Zhu et al., 2008). The methodology has been successfully applied for several applications, both in its pure granular form (Elfeky et al., 2020; Tian et al., 2019; Zhao et al., 2018) or for multiphase flow units (Golshan et al., 2017; Pietsch et al., 2018), including for geological applications (Wang et al., 2019). In spite of the previous positive points, the DEM still has some open issues and limitations (Berger and Hrenya, 2014; Liu and Hrenya, 2014). Among these, the role and correct choice of some sub-models are still unclear, including that of the particle-particle contact equations (Di Renzo and Di Maio, 2004; Lupo et al., 2019). For elastic collision forces, either a linear Hookean approach or a Hertzian non-linear approach is generally employed. While the former makes the simulation less computationally intensive, the latter is generally deemed as more accurate (Paulick et al., 2015). Overall, there are no established guidelines on when it is recommendable to employ one model over the other. In this framework, this work aims at comparing the effect of the two contact models in two different situations: the stress distribution of static particles in a cylindrical column and the water fluidisation of segregated particles of different densities. Both comparisons are based on previously published experimental studies (Di Felice and Scapinello, 2010a, 2010b). Such a direct comparison between the two is quite lacking in the literature (Coetzee, 2017). As will be explained below, these situations present peculiar and possibly unexpected features that are caused by friction forces, which are calculated as a function of the normal contact force. 811 2. Simulation methodology The simulations were performed with the program MFIX (Multiphase Flow with Interphase eXchanges) version 20.2.1 (mfix.netl.doe.gov). It is developed by the National Energy Technology Laboratory of the U.S. Department of Energy and allows simulating granular materials through various methodologies. Most notably, it is open-source, freely available to anyone and transparent to the user. It is also provided with a graphical user interface (GUI) that makes it more user-friendly. 2.1 Governing equations The simulations were based on the DEM, which determines the position and velocity of each computational particle through its Newtonian equations of motion. Every collision between two particles or a particle and the wall is handled by letting them slightly overlap, through what is known as the soft-sphere approach. The magnitude of the overlap determines the intensity of the collision force. For the sake of brevity, the equations of the model are not reported here, but can be found on the MFIX guide (Garg et al., 2012) or in prior publications (e.g. Yurata et al., 2020). The focus of this work is on the modelling of the elastic contact force that, for two colliding particles i and j and in the normal direction, is given as: , = − (1) in which kn is the spring stiffness, δn is the normal overlap and nij is the normal versor. The subscript n indicates the normal direction, and an analogous equation is used for the calculation of the force in the tangential direction. When employing the linear Hookean approach, kn is given as a constant and is provided directly by the user. Conversely, in the Hertzian approach, kn is calculated at every time step and for every collision pair as a function of the overlap and particle properties: , = 43 ∗1 − + (1 − ) , (2) In this, E is the Young modulus, σ is the Poisson’s ratio, rij is the effective collision radius, and G is the shear modulus. This approach clearly entails larger computational requirements due to its additional equations and non-linearity and hence is avoided if unnecessary. Regardless of the approach, the value of the particle stiffness determines the collision time and, therefore, the simulation time step and duration. Employing a realistic value of the stiffness leads to very large computational time. Quite often, researchers employ much lower values, as this has been shown not to alter the simulation results significantly above a certain threshold (Marchelli et al., 2019). For the second batch of simulations, the coupled CFD-DEM methodology was employed to consider the effect of the flow of water on the particles. While the DEM equations are the same, the local mass and momentum balance equations for the fluid phase are added and solved simultaneously. The coupling between the fluid and solid phase is achieved through the drag force. Since these simulations focused on a simple case, the Gidaspow drag model (1994) was employed. It is perhaps the most commonly employed and it often provides acceptable results, although there often are better options (Marchelli et al., 2020). 2.2 Simulation conditions The simulations focused on two different physical situations, but both involved a vertical cylindrical container and are based on previously published experimental data (Di Felice and Scapinello, 2010a, 2010b). Table 1 contains the employed values for the parameters, while Fig. 1 presents snapshots from the two simulations. Table 1: Simulation parameters (the contact parameters are for both interparticle and particle-wall collisions) Parameter Value Parameter Value Spring constant (k) 1000 N/m Average particle filling flux (DEM simulations) 3.18 kg/(cm2 s) Young modulus (E) 107 Pa Poisson’s ratio (σ) 0.29 DEM time step 1/50 of the collision time Friction coefficient 0.3 CFD time step 10-4 s (adaptive) Restitution coefficients 0.9 Water density 1000 kg/m3 Particle diameter 1.7 mm Water viscosity 0.001 Pa s Particle density 2500 kg/m3 (glass) 10800 kg/m3 (lead) Initial bed height (CFD- DEM simulations) 0.12 m 812 Figure 1: (a) Particles falling and stabilising in the DEM simulations; (b) onset of fluidisation in the CFD-DEM simulations (blue is glass while red is lead) The pure DEM simulations focused on the stress distribution of 1.7 mm static glass particles in a 28-mm-diameter cylindrical column. Particles were gradually let into the column through a funnel, until they settled under the effect of gravity. After the system had stabilised, the stress distribution was analysed through the values of the contact forces. To force the particles into reaching the active state, the lateral wall was then let move upwards for 1 s with a velocity of 1 mm/s. After this, the stress distribution was analysed again. The CFD-DEM simulations regarded the water fluidisation of 1.7 mm glass and lead particles in a 30-mm-diameter column. The geometry was discretised in a cut-cell cartesian mesh with a size of 4 times the particle diameter, in accordance with common heuristics (Boyce et al., 2015). At the beginning, particles were settled at the bottom of the column and completely segregated, with the layer of lead particles placed above the layer of glass particles. Water was let flow from the bottom of the column at gradually increasing velocity. The pressure drop at the bottom of the column was measured after it had stabilised, for each inlet velocity value. 3. Results and discussion 3.1 DEM simulations The first group of simulations focused on static particles in a cylindrical column. As known from the literature, and especially from the pioneering work by Janssen (1895), the pressure exerted on the bottom of such a column does not increase linearly with the number of particles, but tends to an asymptotical value. This effect is caused by the particle-wall friction, which sustains part of the particles’ weight. If the forces are properly calculated, the same effect should be visible from DEM simulations. In recent years, this topic has experience moderate popularity in the DEM field (Ramírez-Gómez, 2020) and we focused on it in greater detail in a recent publication (Marchelli and Di Felice, 2020). Fig. 2a shows the trend of the obtained solids pressure at the bottom of the column for different bed heights, from the experiments and simulations, obtained with the two approaches. The results are very evident: while with the Hertzian approach the curve correctly shows the asymptotical behaviour, this is not achieved with the Hookean approach. With the latter, the force exerted on the bottom of the column almost coincides with the weight of the particles, with no significant friction action exerted by the walls. The values are not perfectly matched even with the Hertzian approach, but this is probably due to the impossibility of perfectly reproducing all the phenomena involved in the settling phase, including the filling methodology. In the experiments it is expected that walls slightly expand in response to the force exerted by particles, in addition to the bottom possibly moving downwards if a balance is employed. This leads particles to move towards the so-called active state (Nedderman, 1992). As these effects could not be captured in the simulations, we let the lateral wall slightly move upwards (this is called ‘friction mobilisation’ (Horabik et al., 2018)) to achieve the same effect. Fig. 2b displays the results of this procedure, normalised dividing them by the pressure in a frictionless column. The Hertzian approach leads to the expected result: the vertical 813 movement of the wall increases the particle-wall friction force, and thus the pressure exerted on the bottom decreases. This happens as well with the Hookean approach, but the extent of the pressure decrease is almost negligible compared to the very abrupt change observed with the Hertzian approach. The problem is also not caused by the relatively low value of the spring constant: even with a 100-fold increase of its value, the results are nearly identical. With a further analogous increase, the program stopped functioning properly, with particles keeping bouncing instead of stabilising at the bottom of the column. On the one hand, this suggests that the particle stiffness has a low influence on the results, but on the other hand it clearly proves that the Hookean method is unable to capture the studied physical behaviour. Since higher values of k lead to lower time steps, this probably does not play a role in determining the observed effect. The incorrect friction calculation is probably due to the formation of the settled particles’ configuration: since in the Hertzian approach the value of k is modified iteratively, for every collisional pair and at every time step, it may lead to a more accurate prediction of the particles’ positions and, therefore, of the force network. (a) (b) Figure 2: (a) Effect of the contact model on the vertical solids stress at various bed heights; (b) effect of the contact model and particle stiffness on the bottom solids pressure before and after mobilisation (bed height 7.2 cm). The experimental data are from Di Felice and Scapinello (2010a). 3.2 CFD-DEM simulations The second group of simulations regarded the water fluidisation of a heterogeneous mixture of particles, starting from a completely segregated configuration with the heavier ones at the top. This phenomenon has in common with the previous one the fact that is heavily affected by the friction force, and has indeed been described with the same approach (Di Felice and Scapinello, 2010b). Contrarily to the fluidisation of standard homogeneous mixtures, in this case the bed of solid particles remains in a fixed state even after the pressure drop at the inlet becomes larger than the weight of the particles. This is caused by the particle-wall friction force, which impedes the movement of the heavy particles by balancing the upward force exerted by the lower layer of light particles. Due to the similarity with the previous case, it interesting to assess whether the role of the contact model is the same. Fig. 3 depicts the obtained trend of the water pressure drop as a function of its inlet velocity. Figure 3: Trend of the pressure drop of water as a function of its inlet velocity. The values of the Ergun equation (Ergun, 1952) are calculated with the parameters of Table 1 and a porosity of 0.4. 814 From the Figure, it is evident that both approaches produce nearly identical values of the water pressure drop, and both are in accordance with the value provided by the Ergun equation (Ergun, 1952) for a fixed bed of such composition. As expected in this peculiar case, particles remain in a fixed state even after the pressure drop has exceeded their weight, and only start getting fluidised when the inlet water velocity exceeds a value of 0.0625 m/s. This value, defined as the ‘critical velocity’, is the same through both approaches, and no visible difference is observed. The data are also in agreement with the experimental value: in the experiments, the calculated critical velocity for this mixture was 4.7 % less than that of pure lead, while these simulations led to a value that is 3.8 % lower than the experimental one. Given the similarity of the results, in this case the use of the Hertzian approach appears nonsensical, provided that, even with this moderate number of particles, it made the simulations about 26 times slower. The different performance compared to the first case is probably due to the preponderance of the drag force in the current system: since it is the main force, particles respond to it and the program correctly predicts the force network even if the calculation of the contacts is less accurate. 4. Conclusions This work focused on comparing the linear Hookean and non-linear Hertzian contact modelling approaches in DEM and CFD-DEM simulations through the open-source program MFIX. Two different systems were simulated: static particles in a cylindrical column and the onset of water fluidisation for a completely segregated binary mixture. In both systems, particles are packed and interparticle contact forces play a key role in determining the final output. The simulations show that, in the first case, the Hookean approach does not permit to yield the asymptotic behaviour of the vertical solids stress that is known from practice and often described through the model of Janssen. Regardless of the bed height, the pressure exerted on the bottom of the column coincides with the weight of the particles, and this is not even markedly modified by an upward movement of the lateral wall. Conversely, through the Hertzian approach, the expected asymptotic behaviour is correctly displayed. The bottom pressure becomes almost constant after the ratio between bed height and the column’s diameter exceeds one. Moreover, upon a wall increase, the bottom pressure drastically decreases, becoming less than one third of the original value in response to the increase of friction forces. In the second batch of simulations, employing the CFD-DEM methodology, a similar disparity is not observed: the water pressure drop during fluidisation is nearly identical with both methods, following the values of the Ergun equation. With the layer of heavier particles placed at the bottom, the simulation reproduces what is observed experimentally: the pressure drop visibly exceeds the weight of the particles before a critical value of inlet water velocity is reached and fluidisation starts. This is caused by the particle-wall friction force, which hinders the movement of heavy particles. Considering the similarity of the results, in this case it appears pointless to employ the Hertzian approach, given that it slows down the simulation by 2600 % with no apparent advantage. The difference between the first and the second case may be due to the prominence of the drag force compared to contact forces when water flows. Despite the importance of contact forces in hindering fluidisation, the program correctly predicts the configuration of particles and the force network even with the Hookean approach. In the first case, instead, the Hookean approach leads to a wrong determination of the collisional forces, and particles do not manage to adapt their configuration in response to the falling of other particles. In conclusion, these simulations point in favour of the theory that the Hertzian approach is only vital when particles are closely packed, and gravity is the main force acting on them. Acknowledgments Filippo Marchelli’s current research grant is co-funded by the Liguria Region under the Programma Operativo Por FSE Regione Liguria 2014–2020, grant number RLOF18ASSRIC/42/1. References Berger, K.J., Hrenya, C.M., 2014. Challenges of DEM: II. Wide particle size distributions. Powder Technol. 264, 627–633. https://doi.org/10.1016/j.powtec.2014.04.096 Boyce, C.M., Holland, D.J., Scott, S.A., Dennis, J.S., 2015. Limitations on Fluid Grid Sizing for Using Volume- Averaged Fluid Equations in Discrete Element Models of Fluidized Beds. Ind. Eng. Chem. Res. 54, 10684–10697. https://doi.org/10.1021/acs.iecr.5b03186 815 Coetzee, C.J., 2017. Review: Calibration of the discrete element method. Powder Technol. 310, 104–142. https://doi.org/10.1016/J.POWTEC.2017.01.015 Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Géotechnique 29, 47–65. https://doi.org/10.1680/geot.1979.29.1.47 Di Felice, R., Scapinello, C., 2010a. On the interaction between a fixed bed of solid material and the confining column wall: The Janssen approach. Granul. Matter 12, 49–55. https://doi.org/10.1007/s10035-009-0157-z Di Felice, R., Scapinello, C., 2010b. The transition from the fixed to the fluidised state of binary-solid liquid beds. Chem. Eng. Sci. 65, 5187–5192. https://doi.org/10.1016/j.ces.2010.06.022 Di Renzo, A., Di Maio, F.P., 2004. Comparison of contact-force models for the simulation of collisions in DEM- based granular flow codes. Chem. Eng. Sci. 59, 525–541. https://doi.org/10.1016/j.ces.2003.09.037 Elfeky, K.E., Mohammed, A.G., Wang, Q., 2020. Numerical study on the performance of thermocline tank for hybrid solar tower power plants using DEM-CFD coupled method. Chem. Eng. Trans. 81, 505–510. https://doi.org/10.3303/CET2081085 Ergun, S., 1952. Fluid Flow Through Packed Columns. Chem. Eng. Prog. 48, 89–94. Garg, R., Galvin, J., Li, T., Pannala, S., 2012. Documentation of Open-Source MFIX–DEM Software for Gas- Solids Flows. Gidaspow, D., 1994. Multiphase flow and fluidization : continuum and kinetic theory descriptions. Academic Press. Golshan, S., Esgandari, B., Zarghami, R., 2017. CFD-DEM and TFM Simulations of Spouted Bed. Chem. Eng. Trans. 57, 1249–1254. https://doi.org/10.3303/CET1757209 Guo, Y., Curtis, J.S., 2015. Discrete Element Method Simulations for Complex Granular Flows. Annu. Rev. Fluid Mech. 47, 21–46. https://doi.org/10.1146/annurev-fluid-010814-014644 Horabik, J., Parafiniuk, P., Molenda, M., 2018. Stress profile in bulk of seeds in a shallow model silo as influenced by mobilisation of particle-particle and particle-wall friction: Experiments and DEM simulations. Powder Technol. 327, 320–334. https://doi.org/10.1016/j.powtec.2018.01.003 Janssen, H.A., 1895. No Title. Zeitschr. D. Vereines Dtsch. Ingenieure 39, 1045. Liu, P., Hrenya, C.M., 2014. Challenges of DEM: I. Competing bottlenecks in parallelization of gas-solid flows. Powder Technol. 264, 620–626. https://doi.org/10.1016/j.powtec.2014.04.095 Lupo, M., Sofia, D., Barletta, D., Poletto, M., 2019. Calibration of DEM simulation of cohesive particles. Chem. Eng. Trans. 74, 379–384. https://doi.org/10.3303/CET1974064 Marchelli, F., Di Felice, R., 2021. A Discrete Element Method Study of Solids Stress in Cylindrical Columns Using MFiX. Processes 9, 60. https://doi.org/10.3390/pr9010060 Marchelli, F., Hou, Q., Bosio, B., Arato, E., Yu, A., 2020. Comparison of different drag models in CFD-DEM simulations of spouted beds. Powder Technol. 360, 1253–1270. https://doi.org/10.1016/j.powtec.2019.10.058 Marchelli, F., Moliner, C., Bosio, B., Arato, E., 2019. A CFD-DEM sensitivity analysis: The case of a pseudo- 2D spouted bed. Powder Technol. 353, 409–425. https://doi.org/10.1016/j.powtec.2019.05.035 Paulick, M., Morgeneyer, M., Kwade, A., 2015. Review on the influence of elastic particle properties on DEM simulation results. Powder Technol. 283, 66–76. https://doi.org/10.1016/J.POWTEC.2015.03.040 Pietsch, S., Kieckhefen, P., Heinrich, S., Müller, M., Schönherr, M., Kleine Jäger, F., 2018. CFD-DEM modelling of circulation frequencies and residence times in a prismatic spouted bed. Chem. Eng. Res. Des. 132, 1105–1116. https://doi.org/10.1016/J.CHERD.2018.01.013 Ramírez-Gómez, Á., 2020. The discrete element method in silo/bin research. Recent advances and future trends. Part. Sci. Technol. 38, 210–227. https://doi.org/10.1080/02726351.2018.1536093 Tian, X., Yang, J., Guo, Z., Wang, Q., Sunden, B., 2019. Numerical study of flow and heat transfer in gravity- driven particle flow around a circular or elliptical tube. Chem. Eng. Trans. 76, 235–240. https://doi.org/10.3303/CET1976040 Wang, Y., Zhou, L., Yang, Q., 2019. Hydro-mechanical analysis of calcareous sand with a new shape- dependent fluid-particle drag model integrated into CFD-DEM coupling program. Powder Technol. 344, 108–120. https://doi.org/10.1016/J.POWTEC.2018.12.008 Yurata, T., Piumsomboon, P., Chalermsinsuwan, B., 2020. Effect of contact force modeling parameters on the system hydrodynamics of spouted bed using CFD-DEM simulation and 2k factorial experimental design. Chem. Eng. Res. Des. 153, 401–418. https://doi.org/10.1016/j.cherd.2019.10.025 Zhao, H., An, X., Wu, Y., Qian, Q., 2018. DEM modeling on stress profile and behavior in granular matter. Powder Technol. 323, 149–154. https://doi.org/10.1016/j.powtec.2017.10.006 Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2008. Discrete particle simulation of particulate systems: A review of major applications and findings. Chem. Eng. Sci. 63, 5728–5770. https://doi.org/10.1016/j.ces.2008.08.006 816