Angefertigt am
Department of Particulate Flow Modelling
Erstbeurteiler
Assoz. Univ.-Prof. Dr. Stefan Pirker - Johannes Kepler Uni-versitaet
Zweitbeurteiler
Ass.Prof. Dipl.-Ing. Dr.techn. Michael Harasek - Technis-che Universitaet Wien
March 2016
JOHANNES KEPLER UNIVERSIT ¨AT LINZ
Altenbergerstraße 69
Die Filling of Cohesive
Powders: Material
Char-acterization,
Numerical
Simulation and
Experi-mental Validation
Dissertation
zur Erlangung des akademischen Grades
Doktor der Technischen Wissenschaften
im Doktoratsstudium der
Abstract iii Zusammenfassung v Publications ix List of Figures xi List of Tables xv 1 Introduction 1 1.1 General background . . . 1
1.2 Powder flow in die filling process . . . 2
1.2.1 Physical phenomena of powder flow . . . 2
1.3 Cohesion in granular media . . . 5
1.4 DEM in die filling process . . . 7
1.5 Coarse grain modelling (CGM) . . . 8
1.6 Granular flow through silo . . . 10
1.7 Wall stress in silos . . . 11
1.8 Granular flow under variable g . . . 13
1.9 Aim of this thesis . . . 15
2 Discrete Element Method (DEM) 17 2.1 Concept of DEM simulations . . . 17
2.2 Contact models . . . 19
2.2.1 Hooke model . . . 19
2.2.2 Hertz model . . . 20
2.2.3 JKR model . . . 24
2.2.4 Capillary model . . . 25
2.3 Coarse graining of contact laws . . . 26
2.3.1 Coarse graining of cohesion models . . . 29
2.3.2 Box Shearing . . . 32
2.3.3 Box Filling . . . 39
2.4 Concluding remarks . . . 45
3 Experimental Investigation 47 3.1 Conventional calibration experiments for non-cohesive powders . . . 47
3.1.1 Angle of repose . . . 48
3.1.2 Static angle . . . 49
3.2 Conventional calibration experiments for very cohesive powders . . . 51
3.3.1 Experimental apparatus „Sandy“ . . . 55
3.3.2 Silo Discharge Centrifuge . . . 59
3.4 Concluding remarks . . . 69
4 Numerical Simulation I 71 4.1 Numerical analysis of calibration experiments . . . 71
4.1.1 Static angle . . . 71
4.1.2 Angle of repose . . . 73
4.1.3 Shear cell . . . 74
4.2 Simulation results of Sandy . . . 79
4.3 Numerical analysis of the centrifuge . . . 86
4.3.1 Model parameters sensitivity . . . 87
4.3.2 Numerical Simulation - cohesionless material . . . 89
4.3.3 Similarities analysis . . . 90
4.3.4 Radial and axial velocities . . . 93
4.3.5 Effect of bulk density . . . 97
4.3.6 Effect of polydispersity . . . 100
4.3.7 Effect of cohesion . . . 102
4.4 Concluding remarks . . . 104
5 Numerical Simulation II - Industrial Application 107 5.1 Process description . . . 107
5.2 Numerical modelling . . . 108
5.3 Results . . . 111
5.4 Concluding remarks . . . 129
6 General Conclusions and Suggestions 131 6.1 General conclusion . . . 131
6.2 Suggestions for further works . . . 132
Bibliography 133
This work deals with numerical simulation and experimental validation of cohesive and cohesionless granular material. In the first chapter a review of the die filling process, which is the motivation of this work, is described. Also a review on the main aspects of cohesion in granular media and material flow under variable gravitational force conditions is given. The second chapter describes the Discrete Element Method (DEM), methodology chosen in this work to simulate granular materials. Also the Coarse Graining Model (CGM), a modelling technique to describe the behaviour of fine particles by simulating coarser particles, is described from literature and our contributions for contact and cohesion models are added to the method. Chapter 3 describes a series of experiments used for DEM calibration, which are necessary to obtain physically correct results, such as angle of repose, shear cell and static angle. Furthermore, an experiment that mimics die filling process codenamed “Sandy”, and a centrifuge used to analyze cohesion and cohesionless material flow under increased gravitational force conditions, are described along with their respective obtained data. Chapter 4 describes numerical simulations performed to calibrate our DEM model and validation through comparison of simulations to experiments performed with “Sandy”. Material flow through a hopper is simulated and data is compared to experimental results. Finally, Chapter 5 describes a real size industrial application using DEM model calibrated for the very cohesive Molybdenum powder. Multiple shaking modes are applied to the die and their effect on the final density distribution is analyzed.
Diese Arbeit handelt von numerischer Simulation und experimenteller Validierung des Verhaltens von kohäsiven und nicht kohäsiven granularen Materialien- Im ersten Kapitel wird der Prozess der Formfüllung beschrieben, der die Grundlage für die vorliegende Arbeit darstellt. Weiters wird ein Überblick über die wichtigsten Aspekte des Verhaltens kohäsiver Granulate unter variabler Gravitiation gegeben. Das zweite Kapitel beschreibt die Discrete Element Method (DEM), die in dieser Arbeit gewählt wurde um granulare Materialien zu simulieren. Das Coarse Graining Model (CGM), eine Modellierungstechnik in der das Verhalten kleiner Partikel durch die Simulation größerer Partikel angenähert wird, wird ebenfalls beschrieben, zusammen mit unserer Erweiterung des CGM um ein Kohäsionsmodell. Kapitel 3 handelt von einer Reihe von Experimenten für die Kalibrierung der DEM-Simulationen: Schüttwinkel, Scherzelle und Reibungswinkel. Weiters werden ein Experiment namens "Sandy", das den Formfüllungsprozess nachstellen soll, und eine Zentrifuge in der das Ausfließen von kohäsivem und nicht kohäsivem Material aus einem Silo unter erhöhter Schwerkraft untersucht wird, vorgestellt und Ergebnisse der beiden Experimente werden präsentiert. In Kapitel 4 wird die Kalibrierung der Simulationen mittels der zuvor erwähnten Experimente beschrieben. Weiters werden die Simulationen durch Vergleich mit dem Experiment "Sandy" validiert. Zusätzlich werden die Messungen an der Zentrifuge mit Simulationen verglichen. Kapitle 5 beschreibt Simulationen im industriellen Maßstab mit dem zuvor kalibrierten und validierten DEM-Modell: Eine mit hoch kohäsivem Molybdänpulver gefüllte Form wird auf verschiedene Arten geschüttelt, und der Einfluss auf die Dichteverteilung in der Form wird untersucht.
It took me three and a half years to complete my PhD project, and now I am finally compiling this thesis. I would like to take the completion of this work as a chance to express my gratitude to all the people who have accompanied me during that time.
First of all, I would like to thank Assoz. Univ.-Prof. Dr. Stefan Pirker, head of the Particulate Flow Modelling department, who was a great supervisor throughout my PhD time, not only because of his expertise and enthusiasm about numerical modelling, particulate flow and fluid flow, but also because of his humanity and a good balance in guiding his staff. I would like to thank Ass.Prof. Dipl.-Ing. Dr.techn. Michael Harasek from Vienna University of Technology, Head of Research Division in the Institute of Chemical Engineering for reviewing this thesis.
Also, I would like to thank all my current and former colleagues, especially DI Dr. Christoph Kloss, DI Dr. Christoph Goniva, DI Luca Benvenuti, Mahdi Saeedipour MSc., DI Dr. Stefan Puttinger, DI Philippe Seil, DI Gerhard Holzinger, DI Bernhard König and DI Dr. Thomas Lichtenegger for their support and providing always a very friendly atmosphere.
Of course, I would like to thank those who made this work possible from a financial point of view: Plansee SE, Austria, with special thanks going to Dr. Arno Plankensteiner, for his persistent, under-standing and professional guidance, as well as to Dr. Christian Grohs, who was the main contact person, for sharing his expertise and his feedback always very helpful.
A special thanks to my wife Julia Nasato for her patience, support and for having accepted to go through this new step of our lives with me. A special thanks also to my parents Nilton and Rosana, and to my brother Marcelo for their support.
Journal
Schiochet Nasato D., Goniva C., Pirker S., Kloss C. (2015) Coarse Graining for Large-Scale DEM
Sim-ulations of Particle Flow – An Investigation on Contact and Cohesion Models Procedia Engineering,
Volume 102, Page(s) 1484–1490, 2015.
Proceedings
Schiochet Nasato D., Pirker S., Lichtenegger T., Kloss C. (2015) Numerical simulation of porosity of
co-hesive powders after die filling process in dependence of different pouring and shaking recipes Proceedings
of EURO PM 2015 Congress, 2015.
Schiochet Nasato D., Goniva C., Pirker S., Kloss C. (2014) Coarse Graining for Large-Scale DEM
Sim-ulations of Particle Flow – An Investigation on Contact and Cohesion Models Proceedings of The 7th
World Congress on Particle Technology (WCPT7), 2014.
Pirker S., Schiochet Nasato D., Benvenuti L., Kloss C., Schneiderbauer S. (2014) Simulation
partikelbe-ladener Strömungen: Diskrete, kontinuierliche und hybride Modellierungsansätze Danninger H., Kestler
H., Kolaska H. (Eds.): Neue Horizonte in der Pulvermetallurgie - Werkzeuge, Produkte und Verfahren, Series Pulvermetallurgie in Wissenschaft und Praxis, Volume 30, 2014.
Schiochet Nasato D., Goniva C., König B., Pirker S., Kloss C. (2013) Die Filling Process Simulation
Using Discrete Element Method (DEM) Proceedings of PARTICLES 2013, 2013
Grohs C., Plankensteiner A., Schiochet Nasato D., Kloss C. (2013) Numerical Simulation of Refractory
Metals and Cemented Carbides in the Regime of Powder Filling and Powder Transfer Proceeding on
Plansee Seminar 2013, 2013.
Presentations
Schiochet Nasato D. (2015) Numerical simulation of porosity of cohesive powders after die filling process
in dependence of different pouring and shaking recipes EURO PM 2015 Congress, October 06, Reims,
France
Schiochet Nasato D. (2014) Coarse Graining for Large-Scale DEM Simulations of Particle Flow – An
Investigation on Contact and Cohesion Models 7th World Congress on Particle Technology (WCPT7),
April 20, Beijing, China
Schiochet Nasato D. (2013) Die Filling Process Simulation Using Discrete Element Method (DEM) PAR-TICLES 2013, September 18, Stuttgart, Germany.
1.1 Types of flow in die filling process . . . 2
1.2 Effect of liquid in sand grains . . . 5
1.3 Effect of different amounts of liquid in granular media . . . 6
1.4 Slice element (bulk solid) in the vertical section of the silo . . . 11
1.5 Cylinder equipped with ultrasonic sensor at the top to detect air-grains interface . . . 13
2.1 Illustration of liquid bridges between two grains of different sizes . . . 25
2.2 Coarse graining model . . . 26
2.3 Particle collision test proposed to evaluate normal overlap . . . 30
2.4 Normal overlap of particles using high stiffness . . . 31
2.5 Normal overlap of particles using low stiffness . . . 31
2.6 Characterization of different flow regimes by means of pressure measurements in shear flow 32 2.7 DEM model of the shear box simulation with Lees-Edwards boundary conditions . . . 33
2.8 Dimensionless pressure in the form P Dp/k is depicted for different parcels diameters . . . 34
2.9 Dimensionless shear stress in the form τDp/k is depicted for different parcels diameters . 35 2.10 Pressure for different parcels diameters using Hertz model . . . 37
2.11 Shear stresses for different parcels diameters using Hertz model . . . 38
2.12 Simple box test proposed to evaluate CGM . . . 40
2.13 Coordination and porosity for standard Hertz . . . 40
2.14 Coordination and porosity for Hooke . . . 41
2.15 Coordination and porosity for SJKR2 . . . 42
2.16 Coordination and porosity for JKR . . . 43
2.17 Coordination and porosity for Capillary . . . 44
3.1 Angle of repose test using CFY sample . . . 48
3.2 Static angle experiment - experimental and DEM simulation . . . 49
3.3 Static angle experiment using CFY sample . . . 50
3.4 Angle of repose test using Molybdenum sample . . . 52
3.5 Principle of shear deformation in a translational shear tester . . . 53
3.6 Principle of rotational shear tester . . . 54
3.7 Ring shear tester data from Molybdenum . . . 54
3.8 Details of the testbench proposed to mimic die filling process . . . 55
3.9 Illustration of the “Sandy” testbench steps . . . 56
3.10 Testbench “Sandy” designed to mimic die filling process . . . 56
3.11 Picture of CFY sample tested in “Sandy” . . . 57
3.12 Picture of the Mo sample tested in Sandy . . . 58
3.13 3D CAD model of the centrifuge designed for silo mass flow tests . . . 59
3.14 Centrifuge used for experimental measurements . . . 60
3.15 Details of centrifuge opening mechanism . . . 60
3.19 Measured mass flow dry and wet 900 µm particles . . . 64
3.20 Measured mass flow for Molybdenum sample with 4.6 µm diameter . . . 65
3.21 Large Mo blocks formed in centrifuge operating at 20 g . . . 66
3.22 Large Mo blocks formed in centrifuge operating at 65 g . . . 66
3.23 Large Mo blocks formed in centrifuge operating at 70 g . . . 67
3.24 Measured mass flow for Molybdenum 4.6 µm compared to Beverloo . . . 68
3.25 Mass flow Molybdenum 4.6 µm compared to Carleton . . . 68
4.1 Static angle simulation steps . . . 72
4.2 Results for calibration process of CFY using Sandy the angle of repose test. . . 73
4.3 Shear cell sketch and DEM model . . . 74
4.4 Mo shear cell simulation - friction only . . . 75
4.5 Mo shear cell simulation - coefficient of internal friction . . . 75
4.6 Mo shear cell simulation - servo wall . . . 76
4.7 Mo shear cell simulation - no friction step . . . 77
4.8 Mo shear cell simulation - friction and cohesion . . . 78
4.9 Numerical simulation and experimental measurements of CFY sample using “Sandy” . . . 79
4.10 Numerical simulation of Mo with over predicted cohesion . . . 80
4.11 Numerical simulation results of Sandy testbench for calibrated Mo material. Images of the different steps adopted in the simulation are depicted. . . 81
4.12 Numerical simulation of Mo with calibrated cohesion . . . 81
4.13 Results for calibration process of Molybdenum using Sandy. . . 82
4.14 Experimental data with adjusted Coulomb cohesion . . . 83
4.15 Mo shear cell simulation - calibrated data for shear force . . . 84
4.16 Mo shear cell simulation - calibrated data for coefficient of internal friction . . . 85
4.17 DEM model of the centrifuge . . . 86
4.18 Calibration of dry 900 µm particles - rolling friction . . . 87
4.19 Calibration of dry 900 µm particles - static friction . . . 88
4.20 Mass flow simulations of dry 900 µm . . . 89
4.21 Similarities in flow at 1 g and at increased g-force . . . 90
4.22 Similarities in flow at 1 g and at increased g-force - zoom in dilated region . . . 91
4.23 Dilated region similarities for granular flow under different g-forces . . . 91
4.24 Dilated region similarities for granular flow under different g-forces . . . 92
4.25 Bins for radial velocity . . . 93
4.26 Details of velocity in X, Y and Z direction . . . 93
4.27 X velocity in radial direction . . . 94
4.28 Y velocity in radial direction . . . 94
4.29 Z velocity in radial direction . . . 95
4.30 Bins for axial velocity . . . 95
4.31 Z velocity in axial direction . . . 96
4.32 Solid fraction and g-force for dry 900 µm particles . . . 97
4.33 Mass flow with solid fraction adjustment for dry 900 µm particles . . . 98
4.34 Mass flow - numerical and experimental result for dry 900 µm particles . . . 99
4.35 Mass flow data for mono and polydisperse for dry 900 µm particles . . . 100
4.36 Mass flow data for polydisperse for dry 900 µm particles and experimental measurements 101 4.37 Mass flow under increase gravity for different cohesion forces . . . 102
4.38 Mass flow under increase gravity for highly cohesive forces . . . 103
5.1 Dimensions of the cuboid die . . . 108
5.5 Porosity distribution for the no shaking case . . . 112
5.6 Porosity distribution for the no shaking case - calibrated and uncalibrated models . . . . 112
5.7 Free surface details no shaking . . . 113
5.8 Numerical simulation results for sinusoidal shaking of 5 Hz . . . 113
5.9 Numerical simulation results for sinusoidal shaking of 10 Hz . . . 114
5.10 Free surface profiles for the sinusoidal case . . . 114
5.11 Porosity distribution for the sinusoidal shaking case . . . 115
5.12 Modulated movement of a square wave function . . . 116
5.13 Numerical simulation results for modulated shaking of 5 Hz . . . 117
5.14 Numerical simulation results for modulated shaking of 10 Hz . . . 117
5.15 Free surface profiles for the modulated shaking case . . . 117
5.16 Porosity distribution for the modulated shaking case . . . 118
5.17 Particle size effect for the no shaking case . . . 119
5.18 Wall pressure for the no shaking case . . . 121
5.19 Wall pressure sinusoidal shaking case . . . 121
5.20 Wall pressure modulated shaking case . . . 122
5.21 Wall pressure comparison for different particle sizes . . . 122
5.22 Subcases from no shaking case . . . 123
5.23 Subcases 1/2 from sinusoidal 5 Hz shaking case . . . 124
5.24 Subcases 2/2 from sinusoidal 5 Hz shaking case . . . 124
5.25 Subcases 1/2 from sinusoidal 10 Hz shaking case . . . 125
5.26 Subcases 2/2 from sinusoidal 10 Hz shaking case . . . 125
5.27 Subcases 1/2 from modulated 5 Hz shaking case . . . 127
5.28 Subcases 2/2 from modulated 5 Hz shaking case . . . 127
5.29 Subcases 1/2 from modulated 10 Hz shaking case . . . 128
3.1 CFY powder - measured data through angle of repose and static angle calibration
experi-ments. . . 49
3.2 Mo powder - measured data through calibration experiments. . . 53
3.3 Results for tests using “Sandy” experiment and CFY sample. . . 57
3.4 Results for tests using “Sandy” experiment and Mo sample. . . 57
3.5 Data used in Beverloo equation for large granular material . . . 61
3.6 Data used in Beverloo equation for comparison with experimental data obtained for Mo . 67 3.7 Data used in Carleton equation for comparison with experimental data obtained for Mo . 67 4.1 Calibrated parameters for CFY using angle of repose and static angle tests . . . 73
4.2 Parameters used in cases 1 to 4 depicted in Figure 4.4 . . . 76
4.3 Parameters used in cases 1 to 6 depicted in Figure 4.7 . . . 76
4.4 Parameters used in cases 1 to 6 depicted in Figure 4.8 . . . 77
4.5 Experimental and simulation data for CFY using Sandy test bench . . . 79
4.6 Final calibrated data obtained for Mo after calibration using “Sandy” experimental test rig 84 4.7 Data used for DEM numerical simulations of the centrifuge . . . 88
Introduction
1.1
General background
Powder metallurgy (PM) is a near net shape process, which can produce components with high precision at low cost since there is little waste of material. Also, more and more complex multileveled engineering components are being manufactured using PM technology. This increasingly challenges quality control of the finished part. Consequently, proper handling techniques throughout the entire PM process are required to manufacture consistent and uniform products [83].
Powder compaction and sintering represent the primary industrial manufacturing steps of powder metallurgical processing routes for refractory metals and cemented carbides. Increasing demands in their respective industrial applications encourage net-shape manufacturing aiming at no machining being nec-essary for, e.g., chromium based interconnector plates in solid oxide fuel cell applications, and cemented carbide tools for cutting applications [23].
Inhomogeneous density distributions after die filling are a ubiquitous problem in powder technological part production [10]. In general, a consistent and uniform die filling process is desirable. Heterogeneity during die filling can propagate through the subsequent processes and finally lead to serious product defects, such as cracking, low strength, distortion and shrinkage. Furthermore, a fast filling process is preferable in order to improve the productivity [26].
Material density distribution in die region is influenced by a number of process and materials pa-rameters. Powder characteristics (particle size, size distribution, density, shape and surface properties), apparatus features (shoe and die design) and operating conditions (shoe kinematics, the absence or pres-ence of air, suction, vibration, agitation, aeration, humidity and temperature) are known to influpres-ence the material density distribution. [26]. Understanding the die filling process can provide essential guidelines to optimize materials, equipment and processes in terms of uniformity and productivity [26].
For a better understanding of die filling process, Discrete Element Method (DEM) has been used as a tool for simulating the discrete behaviour of granular materials involved in such process. DEM modelling can provide some micro-dynamic information, e.g. individual particle trajectories and transient inter-action forces, which are usually impossible to obtain in the physical experiments. Such micro-dynamic information is essential to understand the underlying physics of granular matter. Compared to experi-mental investigations, it is much easier to undertake the parametric study in the DEM modelling because the system parameters (particle properties, die geometry, shoe speed etc.) can be prescribed freely [26].
1.2
Powder flow in die filling process
Powder compaction and sintering are important techniques for the mass production of geometrically complex parts. Powder is poured from a reservoir into the feeding shoe, which then passes the cavity one or more times thereby delivering powder into it. The next step is an uniaxial compaction of the powder creating a relatively brittle green body. Finally, the green body is ejected from the cavity and sintered in a furnace where thermal activation below the melting point produces a fully dense structure. Necks form and grow between adjacent grains thereby eliminating the porosity of the manufactured part [10].
1.2.1
Physical phenomena of powder flow
Powder flow from a feeding shoe into a die can be classified in 3 major types: nose, bulk and intermittent flow [26, 72]. When a fixed mass of powder is placed in a shoe, the initial acceleration of the shoe and friction between the powder mass and the base plate forces the powder towards the back of the shoe, forming a nose shaped profile. As the tip of this nose translates across the die opening, material can flow over the surface of the nose to the tip and avalanche into the die. This is referred as nose flow (Figure 1.1 left). At high speeds, or if the die opening is small, the tip of the nose rapidly moves across the opening and powder is delivered into the die by detaching from the bottom free surface of the powder mass. This is referred as bulk flow (Figure 1.1 center).
For some materials, flow occurs as a result of a series of discrete instabilities, which releases either small clusters of particles or large chunks of agglomerated powder into the die. This is a random, of-ten infrequent, process and is referred as intermitof-tent flow (Figure 1.1 right). The intermitof-tent flow is illustrated as the random detachment of individual particles. This flow is especially depicted in strongly cohesive powders which tend to agglomerate.
Figure 1.1: Representation of types of flows. Nose flow (left), bulk flow (center), intermittent flow (right). Schneider et al. [72]
1.2.1.1 Influence of particle properties
The influence of the different material and process parameters in the material density distribution is reported in the literature. Guo [26] reported a lower critical shoe velocity for small irregular shaped particles when compared to regular shaped particles. Critical shoe velocity is defined as the maximum shoe velocity at which the die can be completely filled. This lower critical shoe velocity for irregular shaped particles can be attributed to strong interaction and interlocking of such particles.
Wu and Cocks [83] reported that the addition of liquid reduced flow rates due to increased higher shear resistance between the particles as they flow into the die. However they depicted denser packing for systems with addition of lubricants since inter-particle friction is reduced in the presence of lubricants.
Wu and Cocks [83] also found through numerical simulations that lower inter-particle friction leads to a higher flow rate – particle can slide easily through each other and have a better flowability. Particles with shape that can tessellate (rectangle, hexagon) have a lower critical velocity and filling rate – due to interlocking - than those that can not tessellate (pentagon, circle). Also polydisperse system were reported to have a higher critical velocity and filling rate than monodisperse system [26].
1.2.1.2 Influence of air
Studies found in literature (Guo et al. [25], Schneider et al. [72], Wu et al. [84]) compare results of physical die filling experiments in the presence of air and in vacuum. In general a smoothly and rapidly powder flow is depicted in vacuum except for some cohesive pharmaceutical powders [72]. Schneider et al. [72] reported that air pressure builds up in the die as it is filled when die filling process is conducted in the presence of air. This opposes the flow of powder and reduces the rate at which powder is delivered to the die. Air is also reported to lead in some cases to a transition from bulk to intermittent flow.
Wu et al. [84] suggested the air present during die filling has four major effects on the process: i) the expansion of the powder as it falls from the shoe into the die creates an adverse pressure gradient opposing the motion of the particles; ii) pressure built up by air trapped in the die also opposes the flow of powder into the die; iii) air flowing between the particles provides a lubricant effect, allowing the particles to effectively slide past each other more easily; iv) the drag force exerted by the air as it escapes from the die also slows down the process with some particles even being carried out with the upward flow or air in some cases.
Air has also effects on the segregation of powders. Guo et al. [25] reported that vertical segregation occurred with a high concentration of fine particles at the bottom of the die when filling in vacuum was performed. This is explained by the fine particles being sieved by the coarser particles. The presence of air reduced the extend of this segregation by suppressing the percolation of fines through the voids.
1.2.1.3 Influence of shoe
Shoe kinematics have an important role in die filling. When a low shoe velocity is used nose flow dom-inates the filling process ([72, 84]). This is known to promote fast air evacuation and rapid powder delivery. However powder flows from the feeding shoe only through a small portion of the die opening. For the case when high shoe velocity is used, the tip of the nose moves across the die quickly and the bulk flow dominates the process. The net filling rate in the case decreases as shoe velocity increases [26]. Air-replacement shoe are reported by Hashimoto et al. [28], Sawayama and Seki [71]. The idea is to facilitate the rapid air evacuation of the die and a pipe is positioned in the bottom of the die to perform this task. Such system significantly improved the filling rate ([28, 71]). Hashimoto et al. [28] found im-provements in the filling rate when an agitating device was used in the shoe. Depending on the rotating speed, agitating speed and die geometry used, density variation was reduce up to 30% for agitated feed shoe when compared with conventional feed shoe.
Zahrah et al. [85] conducted experiment using a fluidized shoe. Fluidized shoe causes a lubricant effect in particles and can considerably improve flowability, especially in powders with poor flow properties. As a result higher filling rates were obtained. Maximum compact density and height variation values were significantly reduced, especially in fast filling. However fluidization was reported to cause segregation of powder blends [85].
Reports from Hjortsberg and Bergquist [32] suggest that the powder density increases as the shoe speed increases for a ring shaped die. During die filling, the kinetic energy of moving particles can be partially transferred into the die cavity, in which it causes rearrangement and compaction of the particles. Rice and Tengzelius [64] also reported that the packing density in the die was found to increase with the
number of shoe passages until a limit was reached. This is due to the interactions of the moving particles in the shoe and the stationary particles in the die when the shoe moves over the die, resulting in a denser packing rearrangement of particles in the cavity.
1.3
Cohesion in granular media
Dry and wet particulate media might be greatly influenced by cohesive forces. In case of fine dry granular media van der Waals forces often dominate interactions and adhesion between fine particles. This, in turn, decisively influence the bulk behaviour of powders when they are in close proximity and throughout contact [45, 74]. Parteli et al. [57] reported that packing behavior of powders can be strongly influenced by inter-particle attractive forces of different types, such as adhesion and non-bonded van der Waals forces. For very fine particles porosity can reach values of 0.8 due to cohesive forces.
The van der Waals forces at solid interfaces that occurs at the atomic level were integrated by Hamaker [27] to predict the attraction between two macroscopic non-deformable bodies. For particles which are not touching, the attractive force between two particles whose surfaces are a distance z apart is calculated using Hamaker’s expression [74]:
FHam=
AR∗
12z2 (1.1)
where A is the Hamaker constant which depends on the chemical properties of the solid, and R∗ is the
effective radius defined as [74]:
1 R∗ = 1 R1 + 1 R2 (1.2)
In Equation (1.1) there are no terms which take particles’ deformation into calculations. Johnson et al. [39] proposed a model that considers the Hertzian elastic behaviour of particles - the so called JKR model. The JKR model calculates the surface force acting over the contact area, causing deformation at the contact point. This model will be discussed in details in the next chapter.
Another very important source of cohesive force is the presence of moisture in granular media. In these cases the liquid bridges formed between particles may have an important role depending on the particle size. Even humidity in the air may result in a tiny liquid bridge at a contact point, which introduces cohesion. The effect of cohesive capillary forces in granular media can be easily depict by adding small amount of water (less than 4% w/w) to sand grains as for example in Figure 1.2.
Figure 1.2: Effect of liquid in sand grains. (a) Dry sandpile with a well-defined surface angle. (b) Wet sandpile with a tunnel Mitarai and Nori [54].
The presence of water causes different effects in granular materials according to the amount of liquid in the system. Cohesion occurs in wet granular material unless the system becomes over wet, i.e. the
granular medium is completely immersed in a liquid. Four states of liquid content can be identified in granular media: pendular, funicular, capillary, and slurry state. They can be distinguished in the following:
• Pendular state: Particles are held together by liquid bridges at their contact points.
• Funicular state: Some pores are fully saturated by liquid, but there still remain voids filled with air.
• Capillary state: All voids between particles are filled with liquid, but the surface liquid is drawn back into the pores under capillary action.
• Slurry state: Particles are fully immersed in liquid and the surface of liquid is convex, i.e. no capillary action at the surface.
Figure 1.3 details the effect of the amount of liquid in the granular media and a description of the physical effects.
Figure 1.3: Effect of different amounts of liquid in granular media. In the schematic diagrams in the third column, the filled circles represent the grains and the grey regions represent the interstitial liquid. Mitarai and Nori [54]
1.4
DEM in die filling process
Discrete Element Method (DEM) is a relatively new technique for simulating the die filling process. A general goal when studying die filling is to address the filling density in the die region, as this directly affect the quality of the final product. DEM can model particle size distribution, complex particle geom-etry, complex cavity geomgeom-etry, cavity vibration and other key parameters known to affect the material density distribution. In literature we found the work of Wu and Cocks [83], who used DEM for cavity filling simulations to investigate the influence of escaping air and cavity geometry. Guo et al. [24] studied the effect of air entrapment during the filling process using CFD and DEM coupled simulations. He also analyzed the effect of polydispersity and adhesion in the flow rate of the die filling. Guo [26] investigated the effect of suction applied in cavity region and how size distribution induces segregation during filling. DEM methodology simulate the behaviour of each individual particle in the system and is, therefore, limited to the computational resources available, as the increase in the number of particles also increase the computational costs. A typical die filling process can result, in many cases, in a prohibitive number of particles to be simulated using real particle size. To surpass this limitation a coarse graining technique can be applied for DEM (Bierwisch [11]). The idea behind this technique is to enlarge particles and, as consequence reduce the total number of particles to be simulated. However care has to be taken to guarantee that same physical phenomena is depicted in the scaled and in the original system by the numerical models used. Bierwisch [11] also studied the effect of non spherical particles and cohesion in the material density distribution at the cavity region using DEM.
1.5
Coarse grain modelling (CGM)
Coarse grain model is a term found when referring to multiscale modeling. The term multiscale mod-eling is widely used to describe a hierarchy of simulation approaches to treating systems which bridges across different scales. For a given scale of interest, one picks a simulation method capable of simulating systems at the length and time scales of interest. When one moves to larger scales, a coarse grained model is required to enable simulations. In terms of modelling it is necessary to develop ways of thinking about coarse-graining, and how to link coarse-grained models to the underlying small-scale physics that drive larger-scale dynamics. Information through scales can be propagated in a bottom-up approach where fundamental physical principles at the more detailed scale are used to parametrize a model at a coarse-grained scale. Or in a top-down approach where the behavior at larger scales is used to inform the interactions at more detailed scales [61, 68].
The conceptual idea of coarse grain models is to obtain information of one size or time scale by means of simulations of a different scale. Different applications and approaches for coarse grained modelling are found in literature. Coarse-grained molecular dynamics (CG-MD) are models where groups of atoms are clustered into beads that interact through an effective force-field. This allows simulating mesoscale physical processes while retaining the molecular detail of the system [66]. Many CG-MD models differ by the number of atoms they incorporate. A model reported by [14] maps 4 to 6 atoms in 1 coarse grained bead, while the model reported by [34] maps 200 atoms to a single coarse grained stick with a determined orientation. There are also lattice coarse grained models [44] and others off-lattice models [7]. Ruiz et al. [66] proposed a coarse grained model for molecular dynamics simulation of graphene.
In the Discrete Element Method (DEM) our major interest is in the reduction in the number of parti-cles necessary to simulate a determined process. Since DEM tracks each individual particle in the system and calculate their collisions, the computational power required is proportional to the amount of particles to be simulated.
CGM specifically for DEM is found, e.g., in the work of Sakai and Koshizuka [69]. They studied the pneumatic conveying of millimeter size powders. Sakai and Koshizuka [69] consider that kinetic energy of the coarse grained particle agrees with that of the original particles. They also consider that when the bi-nary collision of the coarse grain particle occurs, the bibi-nary collisions due to all the original particles occur simultaneously. Through a balance of forces they propose the usage of stiffness and damping parameters from the original system to calculate the collisional forces of the CGM system. Forces are calculated for the original system and further applied to the CGM system by means of a scale factor (l3), where l in this
case is the scaling factor (CGM divided by original diameter). Forces applied to the CGM systema are then the forces calculated for the original system multiplied by the scale factor (l3). The disadvantage
in this case is the usage of same time step for the original and CGM particles limiting gains in terms of simulation speed up. Further validations of this CGM model using a fluidized bed are also reported in [70]. Similar to [69] the work of Hilton and Cleary [30] considers the stiffness and damping parameters of the fine particles when calculating the collisions of coarse particles. Hilton and Cleary [30] extend their analysis for tangential collisional forces in a similar way, with the elastic and dissipative constants calculated from the fine particles and extended to larger particles by the addition of a CGM factor.
Radl et al. [62] propose scaling directly the contact law to turn contact forces independent of par-ticle size. Scaling is based on conserving energy from the original to the scaled system. In their work translational and rotational velocity are assumed to remain invariant as well as friction and restitution parameters. [62] start their analysis from the equation of the normal overlap and dimensionless numbers are used to make the normal overlap equation scale independent. Elastic constant should scale with the particle radius and dissipative constant should scale with the square of the particle radius to satisfy previously conditions. Radl et al. [62] found these scaling conditions work fine for quasi-static dense
regime but for inertial dilute regime stresses increase with the increase in particle size. [62] proposes then a modification in the contact law based on the flow regime that is being simulated. One main advantage of this approach is the increase in maximum time step as enlarged particle radius are used to calculate forces during collisions.
1.6
Granular flow through silo
A hopper is the conical, or converging, section of a powder storage vessel. The bin is the parallel sided section, usually cylindrical or rectangular, and the word silo is used to cover the entire vessel. In mass flow type silo discharge, the flow pattern is often described as: first in, first out and in core or funnel flow, the pattern is last in, first out. A silo may acceptably operate under funnel or core flow conditions so long as piping, or rat holing doesn’t occur. These two terms refer to a silo in which the storage capacity consists substantially of stationary powder with just a hole within the silo, taking newly deposited material from the top straight to the bottom and discharge [33].
In contrast with the discharge of a Newtonian fluid for which the flow velocity depends on the height of fluid column, the discharge rate of granular silos is constant. Numerous correlations have been de-veloped over the years to predict the discharge rate of bulk particulate materials from hoppers [6]. The independence of the discharge rate on the filling height is best demonstrated by the well-known Beverloo scaling given by [6, 9]:
Q= Cρb
√
g(L − kd)(5/2) (1.3)
where Q is the mass flow rate, ρb is the bulk density, d is the mean particle diameter and C and k are
constants. The Beverloo equation Equation (1.3) relates the mass flow rate (Q) to the outlet size (L) in a very robust manner for flat bottom silos. Berverloo scaling is commonly accepted as the evidence of a screening effect responsible for a constant and low value of the pressure (i.e., lower than the hy-drostatic prediction) in the outlet area. The traditional physical explanation for this screening effect resorts to Janssen’s analysis: the friction forces mobilized at walls reduce the apparent weight of the ma-terial in the silo and prevent the bottom area to sense the pressure, now partly sustained by the walls [76]. From inspection of Equation (1.3) the dependence of the mass flow with the gravitational force in the form√gis depicted. Frictional effects in the Beverloo equation are incorporated to the constants C and
k. Also as reported by [6] the particle-wall friction has none or very little effect in the mass flow. Beverloo correlation does not consider the drag force in particles in the region near the outflow. This way Berveloo is applicable only for particles with dp ≥500µm. For smaller particles other correlations
are necessary. Particle beds need to dilate (increase distance between particles) before the powder can flow. This means air must penetrate into the bed through the bottom surface of the hopper as the powder moves through the constriction formed by the conical walls. For fine particles the pore diameters in the powder bed are small and there is a significant amount of air drag that resists the powder motion.
Carleton [13] proposed a correlation that considers the effect of air in the vicinity of the outflow region. This is given by Equation (1.5):
4V2 0sinθ B + 15 ρ1/3µ2/3V04/3 ρpd 5/3 p = g (1.4) QCar = ρ0AV0 (1.5)
where V0 is the average velocity of solids discharging (m/s), A and B are geometrical parameters. For
cylindrical hoppers A is the outflow area (πD2
out/4) and B is the outflow diameter Dout. ρ and µ are the
1.7
Wall stress in silos
The classical method to calculate stresses in a silo’s vertical section (active state of stress) was derived by Janssen [36]. Janssen proposed a model that establishes the relation between the average stress at the bottom and the powder bed height in a cylindrical column, based on three major hypotheses: i) at the measurement scale, the granular bed is considered as a continuous medium; ii) the vertical stress is supposed to be horizontally uniform; iii) the grains at the wall are at the sliding threshold [47].
Janssen started his derivation by considering a slice element, i.e., an elemental section of the vertical section, of infinitesimal height, dz as depicted in Figure 1.4.
Figure 1.4: Slice element (bulk solid) in the vertical section of the silo (Schulze [73]).
With the assumptions of constant vertical stress, σv, acting across the cross-section, and constant
bulk density, ρb, equilibrium of forces in z-direction is given by Equation (1.6) [73].
Aσv+ gρbAdz= A(σv+ dσv) + τwU dz (1.6)
Introducing the wall friction angle (assumed to be fully mobilized):
tanϕx= τw/σh (1.7)
and the lateral stress ratio:
K= σh/σv (1.8)
one obtains a differential equation for the vertical stress, σv [73]:
dσv
dz + σvK U
Atanϕx= gρb (1.9)
Assuming constant values of parameters ρb, ϕx, and K , this is a first order differential equation which
can be solved analytically. With the condition that a vertical surcharge stress, σv0, is acting on the top
surface (z = 0), integration of the differential equation yields:
σv= gρbA KtanϕxU +σv0− gρbA KtanϕxU e−KtanϕxUzA (1.10)
For the case where the stress on the top surface is zero ( σv0 = 0 at z = 0), Equation (1.10) yields the
well-known “Janssen equation” given by Equation (1.11)
σv= gρbA KtanϕxU h1 − e−KtanϕxUz A i (1.11)
σhand τw follow from Equation (1.10) by substituting K and tanϕx, respectively, using Equation (1.7)
and Equation (1.8) [73]: σh= gρbA tanϕxU + Kσv0− gρbA tanϕxU e−KtanϕxUzA (1.12) τw= gρbA U + Ktanϕxσv0− gρbA U e−KtanϕxUzA (1.13)
For large values of z the exponential functions in Equation (1.10) to Equation (1.13) approach zero. Thus, the expression in front of the brackets is the final (asymptotic) stress, which is attained for z → ∞. The final value of the vertical stress is:
σv∞= σv(z → ∞) =
gρbA
KtanϕxU (1.14)
σv∞ is independent of silo height and surcharge stress, σv0. It depends on the bulk solid properties and
the ratio A/U, which is A/U = D/4 in a cylinder of diameter D.
Janssen’s work was important because it showed that stress is not transmitted in a similar way to hydraulic head, and wall friction has a very significant influence on the internal stresses of the granu-lar media. However, the assumption of a constant coefficient linking the vertical and horizontal stresses (Equation (1.8)) has no theoretical justification [33]. Also, arches can be formed among particles, suggest-ing that the surface of interest is not planar and that the stress in a plane is not uniform. Nevertheless, it provides a useful semi-theoretical analysis of stress inside a hopper [33].
1.8
Granular flow under variable g
Gravity is the determinant factor in most granular shear flows from geologic situations like landslides and sand piles on a beach to industrial applications (processing of grains, ores, and pharmaceuticals). However, the effect of changing the gravitational acceleration is still largely unexplored [12]. In the work of Brucks et al. [12] they explore the effect of a tumbler connected to a centrifuge in environments with an acceleration of up to 25 g. Angle of repose at the center of the tumbler is depicted and a decrease in its value with the increasing levels of g-force could be seen. They found that scaling relations for granular surface flows developed on Earth at 1 g could be extended to other gravitational accelerations when scaled according to Froude number [12].
In the work of Dorbolo et al. [18] they designed a 2D and a 3D experiment to measure mass flow through a hopper under increased gravity. 3D system consists of a cylindrical hopper filled with non cohesive material and mass flow is obtained through measuring the free surface level of the material inside the container. As the outflow is constant, the free surface position h(t) can be fitted by a linear law h(t) = a(t0− t) where a and t0are fitting parameters. Mass flow is then calculated from Q = aπR2c.
An schematic of the mass flow compartment and free surface level is depicted in Figure 1.5.
Figure 1.5: Cylinder equipped with ultrasonic sensor at the top. These sensors are able to detect the air-grains interface. The receptacle to collect material is a mobile part and is depicted as the cylinder in red. Dorbolo et al. [18].
Dorbolo et al. [18] reported experimental results for 2D and for 3D conditions and for low and high apparent gravity environment, the square root dependence of the flow on the gravity can be considered as a fact (as predicted by Beverloo equation).
Another interesting result reported by Dorbolo et al. [18] was that the apparent gravity does not significantly increase the minimum aperture size below which the flow is clogged. To increase the flow rate out of a silo, the silo can be centrifuged without any risk of clogging.
More recently Mathews and Wu [51] studied the effect of gravity on silo discharge and internal flow patterns using a geotechnical centrifuge. They reported similar findings to [18] with the discharge rate proportional to the square root of the gravity. Mathews and Wu [51] studied 3 different materials - two sand samples with average particle size of 0.4 mm and 0.85 mm and angle of internal friction of 35°. the other sample is a mixture of glass beads with particles of 1.45 mm and 3.15 mm and an angle of internal friction of 22°. To be comparable mass flow is adimensionalised through Equation (1.15):
Qequivalent= Qobserved/(QBev/
√
g) (1.15)
All samples had a good agreement with Beverloo equation with respect to the square root dependence of the flow on the gravity. Differences are found in the absolute value of mass flow and the Beverloo equation. However these differences are expected as the constants in Beverloo equation were not fitted to each material tested [51].
1.9
Aim of this thesis
This thesis focuses on the experimental and numerical simulation of the die filling process using cohesive metallic powder. Very fine particles (< 5µm) of Molybdenum are object of our study. At this size scale cohesive forces have a dominant effect in particle dynamics as well as in static systems such as packed systems [57]. To perform physically acceptable numerical simulations of these fine particles we extend the coarse graining modelling proposed by Radl et al. [62] for non linear contact model (Hertz-Mindling) as well as for cohesion models - JKR and capillary models.
Finding correct DEM parameters for frictional and cohesive data is not trivial specially for highly cohesive material and is also part of this thesis. To our best knowledge modelling die filling process with highly cohesive powders is still an open topic. In addition to calibration experiments for DEM we built an experimental apparatus that mimics the die filling process. This is then used as a validation experiment for cohesive and non cohesive powders.
We also proposed the investigation of material behaviour under increased gravity through granular silo flow. A centrifuge was built to investigate cohesive and non cohesive material flow under acceleration forces of up to 70 g. Molybdenum powder is also experimentally investigated. Cohesive and non cohesive material are then simulated through Discrete Element Method (DEM) and the sensitivity of mass flow with respect to frictional parameters as well as bulk density is investigated.
Finally we simulated an industrial die filling process using Molybdenum as filling material. It was investigated how the different shaking modes affect the porosity distribution in a cuboid shape die. Furthermore the different particle sizes and how they affect porosity distribution and wall pressure were investigated.
Discrete Element Method (DEM)
2.1
Concept of DEM simulations
Discrete element method (DEM) was first developed by Cundall and Strack [16] in 1979. The method aims at tracking individual particles or parcels representing a number of particles in the flow domain. Particles interactions are modeled using a soft-sphere model where rigid spheres are allowed to overlap each other at the contact point.
Equations of motion govern the behavior of classical particulate systems by determining the particles’ trajectories. Changes in positions and velocities of the particles are calculated from the integration of Newton’s equation of motion. The governing equations for the translational and rotational motion of particle i with mass mi and moment of inertia Ii can be written as [86]:
mi dvi dt = X j Fcij+X k Fncik + Ffi+ Fgi (2.1) Ii d dtωi= X j Mij (2.2)
where vi and ωi are the translational and angular velocities of particle i, respectively, Fcij and Mij are
the contact force and torque acting on particle i by particle j or walls, Fnc
ik is the non contact force acting
on particle i by particle k or other sources, Ff
i is the particle–fluid interaction force on particle i, and F g i
is the gravitational force.
The behavior of single particles and global properties are of interest in most granular systems, however it is not possible to find analytical expressions for the phase space trajectories of a system. Therefore in DEM simulations the problem is discretized in time as this will be shown following. For simplicity only the force equation that describes the translational movement is considered and the torque equation, which describes the rotational movement, can be treated analogously:
˙ri= vi (2.3)
miv˙i= F (ri, vi, t) (2.4)
DEM comprises various methods to numerically integrate these equations for the case of solid parti-cles, i.e. for short-ranged, dissipative forces. To this end time derivatives are approximated with finite
differences, e.g. forward differences: ˙r ≈ r(t + ∆t) − r (t) ∆t (2.5) ˙ v ≈ v(t + ∆t) − v (t) ∆t (2.6)
Equations (2.5) and (2.6) are equivalent to use Taylor’s expansion and truncating after the second term: r(t + ∆t) = r (t) + ∆t˙r (t) + 1 2∆t 2¨r(t) + O(∆t3) (2.7) = r(t) + ∆tv(t) +12∆t2 a(t) + O(∆t3) (2.8)
Different integration schemes can be adopted to solve this system of linear equations such as Verlet scheme, the velocity Verlet scheme [82], the Newmark beta scheme [20], among others. The one adopted in LIGGGHTS®[41] and in our simulations is the velocity Verlet scheme, which is derived from the Verlet
algorithm. Written in it’s most compact form the velocity Verlet algorithm is given by:
r(t + ∆t) = r(t) + ∆tv(t) +1
2∆t
2a(t) + O(∆t3) (2.9)
v(t + ∆t) = v(t) +∆t
2 a(t) + a(t + ∆t) + O(∆t2) (2.10) A common way to implement Equation (2.9) and Equation (2.10) consists of introducing half-steps for the velocities:
v t+∆t 2 = v(t) +∆t2 a(t) (2.11) r(t + ∆t) = r(t) + ∆tv t+∆t 2 (2.12) v(t + ∆t) = v t+∆t 2 +∆t 2 a(t + ∆t) (2.13) The velocity Verlet algorithm reads [40]:
1. Calculate v t + ∆t 2
from Equation (2.11) 2. Calculate r(t + ∆t) from Equation (2.12)
3. Evaluate interaction forces and calculate a(t + ∆t) 4. Calculate v(t + ∆t) from Equation (2.13)
2.2
Contact models
When two particles collide energy is dissipated and a resulting repulsive force gives origin to a post colli-sional translational and rotational velocity. There are different models to compute energy dissipation and repulsive forces. Contact models are based in a normal and a tangential force with a repulsive (spring) and a dissipative (damping) term. Normal forces contribute to the translational velocity and tangential forces contribute to the rotational velocity as depicted in Equations (2.1) and (2.2).
When two spherical particles i and j with radii ri and rj, respectively, are in contact, they interact
and are allowed to slightly overlap each other. From this overlap one can calculate the contact force and dissipated energy originated from this collision [46]:
δ= (ri+ rj) − (ai− aj) · n (2.14)
δ >0 means that particles are in contact, with the unit vector n = nij = (ai− aj)/|ai− aj| pointing
from j to i. The positions of the particles i and j are ai and aj. The force on particle i, from particle j, at
contact c, can be decomposed into a normal and a tangential part as Fc
i = Fnn+ Ftt, where n · t = 0.
The tangential force leads to a torque like rolling and torsion do [46].
Our work is based in purely elastic models, i.e. after separation of the contacting particles, they recover their initial shape - there is no plastic deformation [31]. Other models consider plastic deformation after some threshold is surpassed. These models usually tune stiffness or cohesion parameters in a way that original position of the particles are not restored after a cycle of loading/unloading is applied. Some examples of these models are found in literature [37, 46, 49, 79]. In general the form of the normal and tangential spring/damping system is described as:
m∗¨δ = knδn+ cn ˙δn + ktδt+ ct˙δt (2.15)
2.2.1
Hooke model
Models differentiate in the form of how the elastic and dissipative constants kn, kt, cnand ctare calculated.
In LIGGGHTS® and the model considered in further studies as the linear model, a linear spring dashpot
model based on Hooke’s Law is adopted. This model is called Hooke model[1] and it considers a linear increment in repulsive force with particle overlap. In Hooke model the elastic constants are function only of material properties and characteristic velocity and are given by [1]:
kn =16 15 √ R∗Y∗ 15m∗V2 char 16√R∗Y∗ (2.16) cn= v u u t 4m∗k n 1 + π ln(e) 2 ≥0 (2.17) where Vchar is the characteristic impact velocity magnitude in the system. For Hooke model tangential
spring and damping coefficients are the same as normal coefficients (kt = kn and ct = cn). Effective
Young’s modulus Y∗, effective radius R∗ and effective mass m∗ are given as:
1 Y∗ = 1 − ν2 1 Y1 + 1 − ν22 Y2 (2.18) 1 R∗ = 1 R1 + 1 R2 (2.19) 1 m∗ = 1 m + 1 m (2.20)
ν is the Poisson’s ratio. One can see that Hooke’s model results in a linear dependency from the normal
overlap when repulsive force is calculated.
2.2.2
Hertz model
More complex models take into account the geometry of the contact and are mostly based in Hertz studies [29]. Hertz first developed a theory of frictional elasticity of a sphere in contact with a wall for the normal contact. Further Mindlin and Deresiewicz [53] complemented the theory adding a force component for the tangential direction. This model called Hertz-Mindlin is a non-linear contact model and the repulsive force has a dependency of δ3/2.
The Hertz-Mindlin is the model adopted in most of our simulations. The choice is due to its non-linear nature and because most cohesion models are also derived for this model. In the work of [5] there are interesting comparisons of linear and Hertz-Mindlin (non-linear) model. In general better results are found when using a non-linear model, but in some particular cases the opposite effect was reported [5]. For the Hertz-Mindlin model the spring and damping parameters are calculated as [1, 5, 80]:
kn= 4 3Y∗ p R∗δ n (2.21) cn = −2 r 5 6β p Snm∗≥0 (2.22) kt= 8G∗ p R∗δ n (2.23) ct= −2 r 5 6β p Stm∗≥0 (2.24)
with the Sn, St and β parameters calculated from:
Sn= 2Y∗ p R∗δ n (2.25) St= 8G∗ p R∗δ n (2.26) β= ln(e) pln2(e) + π2 (2.27)
here Stis a function of δn. There is a dependency of the tangential damping with respect to the normal
overlap. G∗is the equivalent shear modulus given by the following equation:
1 G∗ = 1 − ν2 1 G1 + 1 − ν22 G2 (2.28)
The shear modulus is related to the Young’s modulus in the form Y = 2G(1 + ν). Tangential force is limited by the coefficient of friction µsthrough the Coulomb criteria in the form [60]:
2.2.2.1 Rayleigh time step
DEM uses an explicit methodology in time and the stability of the solution is dependent of the chosen time step. One has also to consider the force propagation in granular media for the correct choice of the time step. When particles collide, this collision is given between particle and neighboring particles and/or walls. However, the propagation of disturbance waves generated far from the particles’ neighboring also affect the particle movement [81].
This problem is solved by selecting a suitably small value for the time step such that, during a single time step, a disturbance can only propagate from a particle to other particles in contact with it. These disturbances propagate in a form of Rayleigh waves along the surface of the solids [78]. The simulation time step adopted for a physically correct behaviour and a stable solutions is a part of Rayleigh time [59], which is taken by energy wave to transverse the smallest element in the particular system. The critical time step is given by the following equation [81]:
∆tc= πRmin vR = πRmin λ r ρ G (2.30)
where Rminis the minimum particle radius, ρ is the particle density, G is the particle shear modulus,
vR is the Rayleigh wave speed and λ can be obtained from [81]:
(2 − λ2)4= 16(1 − λ2) 1 − λ2 1 − 2ν 2(1 − ν) (2.31) which can be approximated by:
λ= 0.8766 + 0.1631ν (2.32)
where ν is the Poisson’s ratio of the particle. The Rayleigh critical time equation is commonly found in literature in the following form:
∆tr=
πRPpGρ
2.2.2.2 Rolling friction
When introducing rotational degrees of freedom to the particles, changes of their angular velocity occur due to the torque that results from the tangential forces. When unconfined, purely spherical models cannot capture formation of stable angles of repose and the spheres behave fluid-like in that sense [10]. However the static behaviour can be restored to some extent by incorporating a torque against the relative rotation between the two contact entities. We need to introduce a coefficient of rolling resistance µrthat
is defined as a dimensionless parameter:
µr= tan(β) (2.34)
where β is the maximum angle of a slope on which the rolling resistance torque counterbalances the torque produced by gravity acting on the body. This coefficient of rolling friction is incorporated to suppress our impossibility to model real particle shapes and surface effects, e.g., particle roughness.
There are different models to describe the rolling resistance torque. The so called Constant Directional Torque (CDT) are a category of models that apply a constant torque on a particle to represent the rolling friction. The direction of the torque is always against the relative rotation between the two contact [3]. The model is given by Equations (2.35) and (2.36) [1]:
Mrf = −µrknδnR∗
ωrel
|ωrel|
(2.35)
ωrel= ωi− ωj (2.36)
where ωi and ωj are the angular velocities of the particles i and j, respectively, ωrel denotes the relative
angular velocity and R∗ is the effective contact radius.
Another category of rolling friction models are the Elastic-Plastic Spring-Dashpot (EPSD) models. The torque consists of two components: a mechanical spring torque and a viscous damping torque. The mechanical spring torque is dependent on the relative rotation between the two contacting entities. The total rolling resistance torque Mr consists of a spring torque Mkr and a viscous damping torque Mdr in
this model [3]:
Mr= Mrk+ M d
r (2.37)
This model is presented by Ai et al. [3] and is based on the work of Jiang et al. [38]. It is appropriate for both one way rolling and cyclic rolling cases. For EPSD model the torque due to the spring Mk
r is calculated as [1, 3]: ∆Mk r = −kr∆θr (2.38) Mr,t+∆tk = Mr,tk + ∆Mrk (2.39) Mr,t+∆tk ≤ Mrm (2.40) Mrm= µrR∗Fn (2.41)
Here krdenotes the rolling stiffness and ∆θris the incremental relative rotation between the particles.
The spring torque is limited by the full mobilization torque Mm
r that is determined by the normal force
Fn and the coefficient of rolling friction µr. The viscous damping torque Mrd is [1, 3]:
Mr,t+∆td = −Cr˙θr if M k r,t+∆t < M m r −f Cr˙θrif M k r,t+∆t =M m r (2.42)
The damping coefficient Crmay be expressed as [1, 3]: Cr= ηrCrcrit (2.43) Crcrit= 2pIrkr (2.44) Ir= 1 Ii+ miri2 + 1 Ij+ mjr2j !−1 (2.45) Here Ii, Ij, mi and mi are the moment of inertia and masses of the particles i and j, respectively.
An alternative elastic-plastic spring-dashpot model called EPSD2 is implemented into LIGGGHTS®
and is the standard model adopted in our simulations. This model is similar to the EPSD model, but in contrast to the original model the rolling stiffness kris defined as [1, 35]:
kr= ktR∗2 (2.46)
where kt is the tangential stiffness obtained from the contact model (will differ from Hooke or Hertz
model). Furthermore, the viscous damping torque Md
r is disabled at all. For this model there is no
need to define a viscous rolling damping ratio as in the case of the EPSD model. This reduce the number of parameters necessary to model the rolling friction coefficient. A more complete description on elastic-plastic spring-dashpot models and the constant directional models as well as comparison cases are depicted in Ai et al. [3].
2.2.3
JKR model
When fine particles are involved and/or moisture exists, contact and non-contact interparticle forces may affect the packing and flow behaviour of particles significantly. Often the non-contact forces involve a combination of three fundamental forces, i.e. the van der Waals force, capillary force and electrostatic force, which can act concurrently or successively to different extents [86]. For contact the adhesive ther-modynamic forces can be relevant and they will be now described in terms of JKR Model [39].
The classical Hertz contact theory is for the elastic deformation of bodies in contact, but neglects the adhesive forces. The theory of Johnson et al. [39], referred to as JKR model, is used to describe elastic isotropic adhesive spheres in contact [58]. This theory assumes that the attractive forces are confined within the area of contact and are zero outside. JKR model extends the Hertz model to two elastic-adhesive spheres by using an energy balance approach.
The contact area predicted by the JKR model is larger than that by Hertz. When two spheres come into contact, the normal force between them will immediately drop to a certain value 8
9F
jkr
c0 , where F jkr c0
is the pull-off force) due to van der Waals attractive forces. The particle velocity will reduce to zero at a point where the contact force reaches a maximum value and the loading stage is complete. In the recovery stage, the stored elastic energy is released and is converted into kinetic energy, which causes the spheres to move in the opposite direction.
All the work done during the loading stage has been recovered when the contact overlap becomes zero. However, at this point, the spheres remain adhered to each other and further work (known as work of adhesion) is required to separate the surfaces[50]. The contact breaks at a negative overlap with the contact force being 5
9F
jkr
c0 . Pull-off force, F jkr
c0 , is the maximum tensile force the contact experiences
and is given by [39].
Fc0jkr=3
2πR∗ξ (2.47) where ξ is the surface energy. The JKR interaction force is given by[50]:
Fnjkr= − 4Y∗ 3R∗a 3−4pπξY∗a3/2 (2.48) Here, the radius of the contact surface, a, is defined and related to the overlap by the expression[50]:
δij = a2 R∗ − 2πaξ E∗ 1/2 (2.49)
