DF Snow Contamination on Cars A DEM Study to Investigate the Influence of Ice Particle Adhesion on the Angle of Re- pose Master’s thesis in Applied Mechanics YESWANTH SAI TANNERU Department of Chemistry and Chemical Engineering CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2020 Master’s thesis 2020 Snow Contamination on Cars A DEM Study to Investigate the Influence of Ice Particle Adhesion on the Angle of Repose YESWANTH SAI TANNERU DF Department of Chemistry and Chemical Engineering Division of Chemical Engineering Chalmers University of Technology Gothenburg, Sweden 2020 Snow Contamination on Cars A DEM Study to Investigate the Influence of Ice Particle Adhesion on the Angle of Repose YESWANTH SAI TANNERU © YESWANTH SAI TANNERU, 2020. Supervisors: Tobias Eidevåg, David Kallin, VOLVO CAR CORPORATION Examiner: Anders Rasmuson, Chemistry and Chemical Engineering Master’s Thesis 2020 Department of Chemistry and Chemical Engineering Division of Chemical Engineering Chemical Engineering Design Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Pile of ice particles formed in an angle of repose simulation for poly-disperse spherical particles. Typeset in LATEX, template by David Frisk Printed by Chalmers Reproservice Gothenburg, Sweden 2020 iv Snow Contamination on Cars A DEM Study to Investigate the Influence of Ice Particle Adhesion on the Angle of Repose YESWANTH SAI TANNERU Department of Chemistry and Chemical Engineering Chalmers University of Technology Abstract Snow/ice particle adhesion on the sensors of an autonomous car is a very important phenomenon which determines the availability of safety systems controlled by these sensors. Therefore, it is necessary to identify the magnitude of snow impingement on the sensors. One of the key influences in snow impingement is the micro-mechanical behavior of snow such as the angle of repose, a bulk property. Angle of repose depends on the surface roughness and cohesion of the ice particles. Previous exper- imental studies have identified the temperature dependency of the angle of repose. It was observed that with the decrease in temperature of the snow the angle of re- pose also decreases i.e. colder snow will have a lower angle of repose. The current study focuses on modeling ice particles to identify different angle of repose values at different temperatures by varying work of adhesion between the ice particles, thereby correlating work of adhesion to temperature of snow. Work of adhesion is a cohesive (contact) parameter used in Johnson-Kendall-Roberts (JKR) model. The JKR model is used for modelling the adhesive part of the contact force calculation. Work of adhesion determines the magnitude of cohesiveness between ice particles. Along with the work of adhesion values, it is also important to correctly model other phenomena like velocity dependent visco-elastic damping, sliding and rolling of ice particles. Forces on all ice particles like body and contact forces are solved using DEM (dis- crete element modelling). The DEM simulations are performed in Altair EDEM, a commercial solver. Several investigations are made to identify the behavior of ice particles. They are: • Implementing drag force on the particles, which results in correct particle velocities of freely falling particles under the influence of still air. • Identifying and implementing the right velocity dependent visco-elastic damp- ing between the ice particles. • Implementing the suitable adhesion model to solve for contact forces. • Investigating several contact parameters like coefficient of restitution, work of adhesion, friction coefficients to get correct angle of repose values for ice particles at different temperatures. The results indicate the comparison of experimental trends done at Jokkmokk i.e., angle of repose versus temperature to the simulation trends done in Altair EDEM i.e., angle of repose versus work of adhesion. The influence of various contact param- eters like velocity dependent visco-elastic damping coefficient, friction coefficients, v work of adhesion on angle of repose are also discussed. The influence of fall height on the angle of repose is discussed. Keywords: Ice particles, angle of repose, work of adhesion, body force(s), contact force(s), cohesiveness, friction coefficients, visco-elastic damping, Altair EDEM, drag force. vi Acknowledgements I would like to thank Volvo Car Corporation for allowing me to work with an inter- esting master’s thesis. I would like to thank the Contamination and Core CFD team at Volvo Cars, especially Tobias Eidevåg and David Kallin for providing me valuable suggestions and discussions. Big thanks to Altair EDEM for the collaboration and providing us the software licenses. Thanks to Callum Bruce, Stephen Cole, David Curry from Altair EDEM for helping us by constantly providing the technical sup- port throughout the thesis. Many thanks to Professor Anders Rasmuson for taking the responsibility of being my examiner for this thesis. Finally, I thank my parents, friends and relatives for giving me moral and emotional support during the tough times of COVID-19. YESWANTH SAI TANNERU, Gothenburg, June - 2020 viii x Contents List of Figures xiii List of Tables xv 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Purpose & Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Assumptions and limitations . . . . . . . . . . . . . . . . . . . . . . . 4 2 Theory 5 2.1 Angle of Repose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Snow Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.1 Road/Untouched/Ground snow . . . . . . . . . . . . . . . . . 6 2.3 Discrete Element Method . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3.1 Hard sphere approach . . . . . . . . . . . . . . . . . . . . . . 6 2.3.2 Soft sphere approach . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Governing equations in DEM . . . . . . . . . . . . . . . . . . . . . . 7 2.4.1 Contact Forces . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.1.1 Other resistance models . . . . . . . . . . . . . . . . 12 2.4.2 Body (non-contact) Forces . . . . . . . . . . . . . . . . . . . . 12 2.4.3 Time-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 EDEM Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5.1 Contact radius (or radius of contact) in EDEM . . . . . . . . 14 3 Methodology 15 3.1 Experimental studies . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.2 Microscopic analysis . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.3 Other tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Particle size distribution . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.1 Creator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.2 Simulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3.3 Analyst . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4 Modeling coefficient of restitution for elastic damping . . . . . . . . . 22 3.5 Methods to reduce simulation run time . . . . . . . . . . . . . . . . . 24 3.5.1 Reducing Young’s modulus . . . . . . . . . . . . . . . . . . . . 24 xi Contents 3.5.2 Implementing dynamic domain (EDEM) . . . . . . . . . . . . 25 4 Results and discussion 27 4.1 Effect of rolling resistance on angle of repose . . . . . . . . . . . . . . 27 4.2 Effect of visco-elastic damping factor . . . . . . . . . . . . . . . . . . 28 4.3 Angle of repose for reduced Young’s modulus . . . . . . . . . . . . . . 29 4.4 Effect of sliding resistance on angle of repose . . . . . . . . . . . . . . 30 4.5 Effect of fall height on angle of repose . . . . . . . . . . . . . . . . . . 31 4.6 Variation of angle of repose with work of adhesion . . . . . . . . . . . 32 5 Conclusion 37 Bibliography 39 xii List of Figures 1.1 Snow impingement/accumulation on the rear-side of a Volvo car model. 1 1.2 Traditional way of doing an angle of repose experiment for non- adhesive particles, showing a cylinder filled with particles lifted-off to form a pile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1 Hard sphere approach (left) showing no deformation of the particles in contact, soft sphere approach (right) showing small overlap between the particles in contact. . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Work flow showing how forces on particles are updated/solved in Dis- crete Element Method (DEM). . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Representation of Hertz-Mindlin contact model between two spherical particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Hertz-Mindlin with JKR model showing contact and pull-off forces. The plot shows contact force (F) as a function of contact overlap (α). Picture credits: [30] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5 Work flow of simulation set-up in EDEM interface. . . . . . . . . . . 14 3.1 Experimental setup for angle of repose at Jokkmokk (Sweden): phys- ical setup (left), schematic view (center), sieves used as filters for particles (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Microscopic images of snow leaving the shaker. Snow which is ready to fall onto the platform. . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Sintering behaviour of snow after forming a pile. Snow pile untouched (left), no disruption in snow pile even if it is moved across the platform (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.4 Particle size distribution of ice particles classified as rounded in Lan- glois et al. [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.5 Hopper and platform geometries in EDEM. Particles falling from hop- per onto the platform (left), top view of platform and hopper (right). 19 3.6 Protractor tool used in EDEM to measure angle of repose of the snow pile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.7 Extrapolation trend-line for quasi-elastic coefficient of restitution given by equation 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 xiii List of Figures 3.8 Possible collision scenarios in the simulations. Particle-wall collision (left), particle-particle (same size) collision with one of the particles rebounding after collision with wall (middle), particle-particle (differ- ent size) collision with one of the particles rebounding after collision with wall (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.9 Impact tests to validate the equation 3.4. Test for normal impact validation (left), tangential impact validation (right). . . . . . . . . . 24 3.10 Dynamic domain implementation in the simulations. Implementation of dynamic domain with falling particles, no frozen particles (left), particles frozen outside dynamic domain (right). . . . . . . . . . . . 25 4.1 Effect of rolling resistance on angle of repose. Work of adhesion is kept constant i.e., 0.218 J/m2. . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Effect of velocity dependent visco-elastic damping on angle of repose. Angle of repose using modelled damping coefficients in section 3.4 (left), angle of repose with a higher damping coefficient of 0.9 (right). 29 4.3 Angle of repose for reduced Young’s modulus. Simulation with origi- nal Young’s modulus (left), reduced Young’s modulus [E∗ = E/100] (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4 Effect of sliding resistance on angle of repose for two work of adhesion values i.e., 0 J/m2 (blue) and 0.218 J/m2 (black). . . . . . . . . . . . 30 4.5 Angle of repose dependency on fall height. Experimental results at -17 ◦C for faceted snow (red), simulation results for spherical ice particles (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.6 Variation of angle of repose with change in work of adhesion val- ues.The angle of repose increases with an increase in work of adhesion. 32 4.7 Comparison of experimental and simulation trends for angle of repose versus temperature (red) and work of adhesion (blue, black) respec- tively. The figure also shows error bars while measuring angle of repose from both experiments and simulations. . . . . . . . . . . . . . 34 4.8 Three spherical particle used in simulations to evaluate the effect of non-sphericity. The model was adopted form [11], where sintered particles were represented as 3 sphere particles. . . . . . . . . . . . . 34 4.9 Comparison of trends of angle of repose versus work of adhesion for spherical particles (red) and equivalent non-spherical particles (black). 35 xiv List of Tables 3.1 Size distribution of spherical particles used in simulations. . . . . . . 19 3.2 Material properties of ice at -12 ◦C. . . . . . . . . . . . . . . . . . . . 19 3.3 Material properties of wall taken in the simulations. . . . . . . . . . . 20 3.4 Particle creation rate for each factory for an overall creation rate of 50, 000 particles per second. . . . . . . . . . . . . . . . . . . . . . . . 20 3.5 Critical velocities and quasi-elastic coefficient of restitution for parti- cle sizes in [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1 Initial simulation settings in EDEM. . . . . . . . . . . . . . . . . . . 27 4.2 Values of work of adhesion chosen for the simulations and the angle of repose values from them. . . . . . . . . . . . . . . . . . . . . . . . 33 xv List of Tables xvi 1 Introduction Contamination on road vehicles like cars is not an uncommon phenomenon. It may come in various forms like dust, snow, stone chipping, water etc, with each of them responsible for one or several effects on a car [17]. Contamination of snow especially road/ground snow occurs in cold regions when a car moves on the path covered by ice particles. An extreme case of snow contamination on rear-side of a car is shown in figure 1.1. In the case of autonomous cars, snow contamination can be very dangerous as the snow build-up may affect the working of autonomous systems like Camera, Radar, Lidar sensors. These sensors govern the availability of active safety systems in the car. Hence, it is necessary to know the snow accumulation on these sensors/active safety systems, to know their availability and develop necessary actions if they are not available. Identifying snow impingement on these active safety sensors like automatic braking, parking assistance (proximity sensors) helps in developing new technology to avoid the snow formation on these sensors. Figure 1.1: Snow impingement/accumulation on the rear-side of a Volvo car model. Development of new car models generally includes analysing for various forms of contamination that affect a car’s performance in several scenarios. At Contamina- tion and Core CFD department Volvo Car Corporation (VCC), continuous research on modelling and simulating the effects due to various contaminants on cars is be- 1 1. Introduction ing done. These simulations are then validated with experiments in field studies and physical tests in Wind tunnel(s) (Climate Wind Tunnel, Aerodynamic Wind Tunnel). The following master’s thesis is a part of a research project at Contamination and Core department at VCC in collaboration with the Chalmers University of Technol- ogy and Luleå University to model snow based on its micro-mechanical properties and use the results to model the snow contamination on cars. 1.1 Background Figure 1.2: Traditional way of doing an angle of repose experiment for non-adhesive particles, showing a cylinder filled with particles lifted-off to form a pile. Snow accretion on cars is a combined effect of both aerodynamic design and the degree of adhesion i.e., snow behaviour at various environments (different tempera- tures and humidity). It is important to know the adhesion properties of ice particles to correctly implement the coupling of snow physics with aerodynamics. Therefore, it is necessary to look at the micro-mechanical properties of snow at particle level 2 1. Introduction (ice particles). Some of the micro-mechanical properties of snow are the angle of repose, internal friction, packing by tapping, age hardening, capillary potential etc [1]. The steepest angle of descent relative to the horizontal plane to which a granular material can be piled is known as the angle of repose. It is one of the important granular properties that decide the degree of adhesion of snow. A traditional way of doing angle of repose can be seen in figure 1.2, where a cylindrical container filled with particles is lifted up to from the angle of repose. This method cannot be used for cohesive particles like ice particles in snow as they can stick together and to the test equipment during the preparation of test. Ice particles also sinter together by time which means that the experiments will depend on how long does the sample has been at rest before doing the test. Hence, it is necessary to do these tests in an alternate way for ice particles. One of the alternate ways of doing these experiments is pouring the particles from an altitude, which was adopted in the current work. Angle of repose experiments were carried out at Northern Sweden (Jokkmokk), to identify the physical behavior of different snow types at various temperatures. They gave an evidence that snow behavior (micro-mechanical properties) changes with change in size of particles in snow, type of particles in snow, temperature of snow, etc. The idea is to model a similar behaviour of snow using numerical simulations to approximately reproduce the experimental results. Discrete Element Modelling (DEM) [4] is chosen to model the contact physics between the particles. The DEM simulations were performed in Altair EDEM, a commercial DEM solver. 1.2 Purpose & Objectives This thesis focuses on performing simulations for Angle of Repose of snow by varying adhesion parameter(s). The idea is to investigate if the trend of angle of repose can be explained by changing adhesion parameter(s) and see if we capture this phenomena. The results can be applied to accurately predict the snow contamination on cars (especially autonomous sensors). Henceforth, new methods/technology can be developed to avoid the snow accumulation on these sensors. Main objectives of this work are: • Propose reasonable contact physics to model the ice particles in snow using DEM approach. • Apply accurate body force(s) on the falling particles. • Model the parameters/coefficients in various damping in the contact models. • Apply method(s) to reduce the simulation time thereby, making the simula- tions to run more fast and accurate. • Perform DEM simulations by varying different parameters to get various angle of repose values to correlate them with the behaviour of snow at different temperatures. 3 1. Introduction 1.3 Assumptions and limitations The current work focuses on the snow which is classified as untouched snow as described in section 2.2. Further investigations are needed to identify the properties of snow which actually causes contamination on cars. As described in [8], the snow which is lifted-off from the ground that is much smaller in size than the untouched snow is responsible for the snow impingement on the exterior of a car. Also, most simulations in current work are done using spherical particles by assuming the snow behaviour described in section 2.2. There are also particles that are irregular and not spherical in untouched snow, shown in figure 3.2. The idea is to compensate non-sphericity of particles by increasing the rolling resistance of spherical particles. A few simulations with non-spherical particles are done to see the effect of various parameters like shape, rolling resistance, sliding resistance of ice particles on the angle of repose. 4 2 Theory This chapter introduces the concepts of discrete element modeling, angle of repose, snow characteristics, and various terminology used in the thesis work. It also in- cludes several physical concepts that were used in the current work. 2.1 Angle of Repose Angle of repose is one of the micro-mechanical properties of granular materials. It can range from 0◦ to 90◦. Angle of repose experiments can be used to calibrate the magnitude of snow impingement on cars, which is the main motto of the current work. As mentioned in section 1.1 angle of repose of ice particles cannot be done traditionally i.e., the procedure used for non-adhesive particles. Kuroiwa et al. [1] has done investigations to evaluate micro-mechanical properties of snow. There they poured ice particles from a container to form a snow pile. It was found that the angle of repose increases with an increase in temperature i.e., colder snow has low angle of repose and vice-versa. There were also investigations made to check the influence of fall height on angle of repose and it was found that for lower fall heights it gave a higher angle of repose and vice-versa. In Kuroiwa et al. [1], the increased stickiness of ice particles was related to the increased liquid layer on the ice particles with an increase in temperature. Current work focuses on correlating temperature values to the work of adhesion to simulate the angle of repose experiments numerically. 2.2 Snow Characteristics Before going into the modelling of snow it is important to know some of the key phenomena happening with snow crystals and the classification of snow, to under- stand the assumptions made for this work. An ice crystal forms when a super-cooled water droplet comes into contact with aerosol particles like soot [25]. This process is known as the nucleation of ice and usually occurs in the troposphere. Water vapour around the ice crystal solidifies onto it and causes nearby particles to evaporate. This process continues and leads to the formation of snowflakes through a process called precipitation. The snowflakes start to grow and reshape and start to fall onto the earth surface. The snowflake shape depends on both temperature and humidity of the atmosphere, shown in Nakaya and Matsumoto [26]. The microstructure of the snowflake changes when it settles on the ground. This phenomenon is known as metamorphism [27, 11]. When there is a low-temperature gradient in the at- 5 2. Theory mosphere, dry snow experiences equilibrium metamorphism which leads in forming rounded grains. When the gradient is high kinetic growth metamorphism occurs which leads to forming hexagonal grains. Snow is a granular material made up of ice particles which come in various shapes and sizes [3]. As snow is a very complex material, several simplifications have to be done to model it. Snow properties like snow adhesion vary a lot with temper- ature, humidity etc. Change in different snowflake structures varies at different temperature and humidity [2]. 2.2.1 Road/Untouched/Ground snow Road snow is the type of snow which is settled on ground for longer periods, unlike the falling snow. The snow undergoes various phenomenon after getting settled on the earth’s surface. One of them is metamorphism which is discussed above. Another phenomenon is the sintering of ice particles [12, 13]. As explained in Colbeck [12], sintering in dry snow is a process where the grains grow slowly while they build inter- granular bonds. The snow that causes contamination on cars has a lower particle size distribution [8] than the road snow as the particles lifted-off from the ground are responsible for the contamination. In the current project, the road/ground snow is of more focus and has been analysed and simulated to find the angle of repose. This was chosen as the snow lifted-off from the ground has a very small particle size distribution which is computationally expensive to simulate. 2.3 Discrete Element Method Discrete Element Method (DEM) is a numerical method used to model or simulate bulk behaviour of granular materials. DEM theory was introduced by Alder and Wainwright in 1956 for molecular dynamics studies and later developed by Cun- dall and Strack [4] to describe the mechanical behaviour of assemblies of discs and spheres. This approach is used to analyze behaviour in a wide range of applications like materials engineering, mechanical engineering, pharmaceutical industry. DEM is useful in investigating phenomena at particle length scales. There are two main approaches in simulating particles in DEM. They are: • Hard sphere approach. • Soft sphere approach. 2.3.1 Hard sphere approach In the hard sphere approach, the particle interactions are impulsive i.e., the particles are in contact with each other for a very small time. As the particle interactions happen in a very short duration, the time-step in this approach is also very small and such simulations take very long simulation time. The particles only exchange momentum through collisions, forces between particles are not solved explicitly. 6 2. Theory 2.3.2 Soft sphere approach In the soft sphere approach, small overlaps are allowed to account for the contact deformations. The time-step in the soft sphere approach is usually larger than that in hard sphere approach. Along with momentum exchange other contact forces on the particles are also solved in this approach. Figure 2.1 exclusively shows the schematic difference between the hard and soft sphere approaches. Soft sphere approach is used to model the ice particles in this project, where the particles are in contact with each other for longer duration, and an overlap between the particles can be modeled. Figure 2.1: Hard sphere approach (left) showing no deformation of the particles in contact, soft sphere approach (right) showing small overlap between the particles in contact. 2.4 Governing equations in DEM The particles in DEM can exhibit two types of motion i.e. translation and rotational, as a result there are 6 degrees of freedom for each particle. The particle accelerations are calculated using Newton’s second law and are integrated numerically to evaluate velocities and positions of the particles in a DEM simulation. Equations 2.1, 2.2 are used to calculate the translation and rotational motions respectively. mdv dt = Fg + Fc + Fnc (2.1) Here, ’m’ is the mass of the particle, ’v’ is the linear velocity of the particle, ’Fg’ is the gravitational force acting on the particle, ’Fc’ and ’Fnc’ are the resultant contact and non-contact forces between the particles (or between particle and wall), ’t’ is time. 7 2. Theory Idω dt = M (2.2) Here, ’I’ is the moment of inertia of the particle, ’ω’ is the angular velocity of the particle, ’M’ is the resultant contact torque acting on the particle, ’t’ is time. The solver in DEM solves forces on each individual particle using Newton’s second law. Several time integration methods can be used to solve the equations. The work flow followed to solve particle forces in DEM can be seen in figure 2.2. Figure 2.2: Work flow showing how forces on particles are updated/solved in Discrete Element Method (DEM). The forces on the particles mentioned above can be split into two categories: • Contact forces. • Non-contact or body forces. 2.4.1 Contact Forces The Hertz-Mindlin model, i.e., Hertzian model (1882) with Mindlin’s no-slip model, is a default contact model to evaluate both normal and tangential contact forces between the particles effectively and accurately. There are damping forces in both normal and tangential directions where velocity dependent visco-elastic damping is related to coefficient of restitution. Each particle pair in contact is evaluated for normal and tangential forces based on the overlap as shown in figure 2.3. 8 2. Theory Figure 2.3: Representation of Hertz-Mindlin contact model between two spherical particles. The normal force is given as, Fn = 4 3E ∗ √ R∗δn 3 2 (2.3) Here, E∗ and R∗ are the equivalent Young’s modulus and equivalent particle radius respectively, defined in equations 2.13 & 2.12. Whereas, δn is the normal overlap between the particles. Damping force in normal direction is given as, Fn d = −2 √ 5 6β √ knm∗vn rel (2.4) Here, m∗ is the equivalent mass of the particles, vnrel is the normal component of relative velocity, kn is the normal stiffness given by kn = 2E∗ √ R∗δn (2.5) β = − ln e√ ln2 e+ π2 (2.6) Here, ’e’ is the coefficient of restitution. m∗ = m1m2 m1 +m2 (2.7) Here, m1, m2 are masses of spheres. Tangential force (Ft) depends on the tangential stiffness (kt) and tangential overlap (δt), given by Ft = −ktδt (2.8) kt = 8G∗ √ R∗δn (2.9) 9 2. Theory Here, G∗ is the equivalent shear modulus. And tangential damping is given by Ft d = −2 √ 5 6β √ ktm∗vt rel (2.10) Here, vtrel is the relative tangential velocity. The fact that the particles in the simulations are adhesive in nature, it is important to chose right adhesive model along with Hertz-Mindlin model. There are two most common models to model adhesive particles, which are Derjaguin-Muller-Toporov (DMT) and Johnson-Kendall-Roberts (JKR) [5]. Based on the Tabor number de- fined in equation 2.11, DMT and JKR models are applicable to different spectrum of particle physics. λT = ( R∗Γ2 E∗2δ3 )( 1 3) (2.11) Here, λT is the Tabor number, Γ is the work of adhesion, δ is the overlap between the particles, R∗ is the effective radius of the contact between the particles given by R∗ = R1R2 R1 +R2 (2.12) Here, R1, R2 are the radii of two particles in a contact. For particle-wall contact R∗ = R1 as R2 =∞. And E∗ is the effective Young’s modulus between the particles given by 1 E∗ = 1− ν1 2 E1 + 1− ν2 2 E2 (2.13) Here, E1, E2 are the Young’s Modulus of two particles, ν is the Poisson’s ratio. As given in Thornton [6] JKR theory is applicable to λT > 5 and DMT is applicable for λT < 0.1. In the current work all the values of Tabor numbers are over 5, therefore JKR model is chosen. A Hertz-Mindlin model with JKR is implemented in the current work. Figure 2.4 shows the variation of contact overlap (α) with the contact force (F) in a Hertz- Mindlin with JKR model. The force acting between the particles is zero from point A to point B. At point B, the particles physically come in contact with each other, where the normal contact force drops down to (8/9)Fc due to the presence of Van der Waals forces at point C. Here, Fc is the pull-off force given in equation 2.14. Upon loading the curve goes from point C to point D. After reaching point D the unloading occurs along the same path. The contact overlap becomes zero at point C, but the spheres stay in contact with each other. A minimum pull-off force shown at point E is required to break the contact and finally the particles separate at point F. 10 2. Theory Figure 2.4: Hertz-Mindlin with JKR model showing contact and pull-off forces. The plot shows contact force (F) as a function of contact overlap (α). Picture credits: [30] Fc = −3 2πΓR∗ (2.14) The contact force in normal direction for Hertz-Mindlin with JKR is given by Fn = 4E∗a3 3R∗ − ( 8πΓE∗a3 )0.5 (2.15) Here, Γ is work of adhesion (also called interfacial surface energy in EDEM), ’a’ is the contact radius. The contact overlap α shown in figure 2.4 is given as, α = a2 R∗ − ( 2πΓa E∗ )0.5 (2.16) As explained in [16] the work of adhesion Γ is given as, Γ = γ1 + γ2 − γ1,2 (2.17) Here, γ1 and γ2 are the surface energies of two spheres and γ1,2 is the interface surface energy. If the spheres have the same material then γ1,2 = 0 and Γ = 2γ. Along with the normal interactions, there are sliding and rolling resistance in tangen- tial and rotational directions respectively. Sliding resistance is applied on tangential force i.e., tangential force is limited by Coulomb friction given in equation 2.18. 11 2. Theory Fsliding = µsFn (2.18) Here, µs is the coefficient of static/sliding friction The rolling resistance is accounted for by applying torque to the contacting surfaces. It is given by, τrolling = −µrFnRiωi (2.19) Here, µr is the rolling friction coefficient, ’Ri’ is the distance of contact point from center of mass, ’ωi’ is the unit angular velocity of the object at contact point. Several other rolling resistance models that can be used for the similar application were described in [15] 2.4.1.1 Other resistance models There are other sliding and rolling resistance models that are used in [7] to correctly model Point ’O’ in figure 2.4. In the current rolling and sliding resistance equations given by 2.19 and 2.18 respectively, when the normal force Fn becomes zero the respective resistance also becomes zero. To compensate this, the equations 2.20 & 2.21 can be used instead. However, the current work was carried out using the aforementioned models i.e., equations 2.18 and 2.19. τrolling = −µr | Fn + 2Fc | Riωi (2.20) Fsliding = µs | Fn + 2Fc | (2.21) 2.4.2 Body (non-contact) Forces The particles in this work are falling freely from the sieves onto the platform, as shown in figure 3.1. Therefore it is enough to model drag resistance on the particles falling freely from the sieves. One of the assumptions in this work is using spherical particles, so the drag equations presented below are for spherical particles. The drag resistance is given by, FDrag = 0.5ρairCDAP | vair − vparticle | (vair − vparticle) (2.22) Here, ρair is the density of air taken as 1.351 kg/m3, vair and vparticle are velocities of air and particle respectively, CD is the drag coefficient given by, CD =  24 Re , if Re ≤ 0.5( 24 Re ) ( 1 + 0.15Re0.687 ) , if 0.5 < Re < 1000 0.44, if Re ≥ 1000 (2.23) Here, Re is the particle Reynolds’s number given by, 12 2. Theory Re = ρparticle | vair − vparticle | (2Rparticle) µair (2.24) Here, µair is the viscosity of air taken as 16.55 × 10−6 Pa s, AP in equation 2.22 is the projected area of the particle given by, AP = πR2 particle, for spheres (2.25) 2.4.3 Time-step Choice of time-step is key while implementing DEM simulations. Time-step deter- mines the stability of the simulations. Choosing improper time-steps will lead to numerical instabilities. As the movement of particle in a granular flow is affected not only by the contacts with its immediate neighbours but also depends on the disturbances by particles far away, which are called Rayleigh surface wave propaga- tion. These disturbance waves can be avoided by choosing small enough time-step. Time-step or Rayleigh time-step for the simulations used in this work is given by equation 2.26. Usually the use of a fraction of this time-step is recommended to ensure accurate force transmission and avoid numerical instabilities. It is common to use 20% of the Rayleigh time-step as mentioned in [29], which was chosen in the simulations. The particle radius mentioned in 2.26 usually is the smallest particle radius in the simulations, given that all particles have same properties. TRayleigh = πR( ρ G )0.5 0.1631ν + 0.8766 (2.26) Here, ’TRayleigh’ is the Rayleigh time-step, ’R’ is the particle radius, ’ρ’ is the density of the particle, ’G’ is the shear modulus of the particle, ’ν’ is the Poisson’s ratio of the particle. G = E 2(1 + ν) (2.27) Large/improper choice of time-step leads to excessive overlaps that may result in large contact forces and also cause Rayleigh waves in the system i.e. in a granular flow, the particle movement can also be affected by disturbance propagation from particles that are far away. 13 2. Theory 2.5 EDEM Overview Figure 2.5: Work flow of simulation set-up in EDEM interface. The contact or non-contact physics in EDEM can be modified or defined by the user using API’s (Application Programming Interface). 2.5.1 Contact radius (or radius of contact) in EDEM In EDEM contact radius which is greater than the physical radius of the particle must be defined. This enables the solver to account for the work of adhesion and allows the influence of negative overlap in force calculations. 14 3 Methodology This section explains the various methods used in the current work to set up accurate snow simulation in EDEM. The details about the assumptions in simulations are described and the reasoning for those assumptions is explained. The angle of repose experiments which are carried out at Jokkmokk are also described in this section. 3.1 Experimental studies The experimental work described in this section is not a part of this thesis project, the results from the experiments are taken to compare and correlate the experimental trends (angle of repose versus temperature) to the simulation trends (angle of repose versus work of adhesion). These experiments were done at Jokkmokk (Sweden). The equipment to do the angle of repose experiments is shown in figure 3.1. 3.1.1 Experimental setup Figure 3.1: Experimental setup for angle of repose at Jokkmokk (Sweden): physical setup (left), schematic view (center), sieves used as filters for particles (right). The experimental setup has a platform on which the particles fall and form a pile. The angle of repose is measured then for that pile. There is a shaker on the top 15 3. Methodology to which the particles are loaded. At the bottom of the shaker, there is a set of sieves, which act as filters to not allow lumps of snow to fall on the platform. When the experiment begins, the shaker starts to vibrate to allow the ice particles to fall on the platform. The vibrations are produced using an electric device, which are required to break the sintering bonds (section 2.2) if there are any. 3.1.2 Microscopic analysis Along with the angle of repose experiments, microscopic images of snow particles were also taken at Jokkmokk to evaluate the size and shape of the particles in snow. The microscopic pictures of ice particles can be seen in figure 3.2 Figure 3.2: Microscopic images of snow leaving the shaker. Snow which is ready to fall onto the platform. These images are necessary to observe and evaluate the characteristics of particles at different temperatures. 3.1.3 Other tests Several other tests were also done at Jokkmokk to see the characteristics and be- haviour of snow which include evaluation of bulk density of snow and test to observe sintering phenomena of particles after forming the pile. The bulk density of snow is measured by filling a container with snow and mass of snow filled into the container was measured. The bulk density is then calculated using both mass of the snow in the container and volume of the container. 16 3. Methodology Figure 3.3: Sintering behaviour of snow after forming a pile. Snow pile untouched (left), no disruption in snow pile even if it is moved across the platform (right). To observe sintering, the snow pile is pushed slowly to see if it gets destroyed due to the applied force. This can be seen in figure 3.3 where the pile stays put even after it is displaced and doesn’t disrupt. 3.2 Particle size distribution Figure 3.4: Particle size distribution of ice particles classified as rounded in Lan- glois et al. [14]. 17 3. Methodology The particle size distribution was adopted from Langlois et al. [14], where 3D image analysis of ice particles of various classifications were performed. The particle size distribution was taken from [14] as it is a more comprehensive study and the observed particle sizes agree with the what is observed in field studies (section 3.1.2). In Langlois et al. [14] a proper 3D image analysis was done for the similar types of ice particles to extract the particle size distributions with more number of particles in the samples. There are many types of particles classified in [14] like rounded, facetted, depth hoar, etc., along with eccentricity distribution for the particles. Particle size distribution of particles classified as rounded with zero eccentricity were chosen to perform the simulations in the current study. See figure 3.4. Total number of particles analysed under this classification i.e., rounded are 50,833 particles. The population density values shown in figure 3.4 are normalized with the total number of particles to get the number of particles at each diameter. Among all the particle sizes, the particles which have a higher population density were chosen to represent the entire size distribution. The rates at which particles fall from the hopper are given in table 3.4, which is calculated based on contribution of particles from this distribution. 3.3 Simulation setup The simulations are setup in EDEM a commercial solver for DEM simulations. This section provides an overview about different modules in EDEM like creator, simu- lator and analyst which are used to set-up, solve and post-process the simulations respectively. 3.3.1 Creator Creator module in EDEM allows the user to setup the simulation. Here, the user can do all the preprocessing work for the simulation i.e, setting-up geometries, particle distribution, properties, etc. To replicate the experiments that were explained in section 3.1.1 the geometries were created in EDEM accordingly, and are named platform and hopper. The ice particles falling from the hopper and are accumulated on the platform at the bottom as shown in figure 3.5. The platform is a cylinder with diameter 25mm and the hopper is a square plate 25mm × 25mm. The geometries are enclosed in a boundary; the particles which are outside this enclosure are not considered in the simulation and will be deleted. A periodic boundary can be set if necessary, which is not applicable for this case. The particles are defined using bulk material option in EDEM, which is used to create particles with different interactions, size, and shape. Five of such categories of spherical particles are created with different sizes as shown in table 3.1, properties of particles and wall are shown in tables 3.2, 3.3 respectively. The choice of particle size distribution is motivated in section 3.2. Spherical particles were chosen assuming the fact that the snow undergoes metamorphism as explained in section 2.2. The 18 3. Methodology Figure 3.5: Hopper and platform geometries in EDEM. Particles falling from hopper onto the platform (left), top view of platform and hopper (right). snow analysed in experiments is ground snow, which is settled for a longer duration. Radius of contact in EDEM is explained in section 2.5.1. Number Particle radius (µm) Radius of contact in EDEM (µm) 1 225 227 2 300 302 3 350 352 4 367 369 5 415 417 Table 3.1: Size distribution of spherical particles used in simulations. Name Property Reference Young’s Modulus, E [GPa] 9.7 [18] Poisson’s Ratio, ν [-] 0.31 [21] Density, ρ [kg/m3] 920 [18] Surface Energy, γ [J/m2] 0.109 [16, 18, 19, 20] Table 3.2: Material properties of ice at -12 ◦C. 19 3. Methodology The platform is defined as a geometry with material properties given in table 3.3. The hopper is defined as a virtual geometry which doesn’t have any material proper- ties assigned. Virtual geometries in EDEM can be assigned as factories. A factory (in EDEM) is a source from which the particles defined in table 3.1 are created. Various parameters like creation rate, position, orientation, velocity, angular veloc- ity, etc can be defined for the particles. In this case, the particles are assumed to be freely falling under the influence of a drag force, see section 2.4.2. Therefore, the particles are given no velocity or angular velocity, the positions and orientations of the particles were defined as random. An overall creation rate of 50, 000 particles per second is given to the 5 particle factories as shown in table 3.4. The fall height in most of the simulations was taken as 30mm. A 30mm fall height was chosen as it is sufficient to accommodate the angle of repose values ranging from 0◦ to 60◦ and is also allows the simulations to finish quicker i.e., a higher fall height would take more simulation time. Property ASA Polymer Reference Young’s Modulus, E [GPa] 2 [22] Poisson’s Ratio, ν [-] 0.35 [23] Surface Energy, γ [J/m2] 0.035 [24] Table 3.3: Material properties of wall taken in the simulations. Factory number Particle type (radius) (µm) Creation Rate (1/sec) 1 225 11, 836 2 300 11, 015 3 350 9, 952 4 367 9, 343 5 415 7, 854 Table 3.4: Particle creation rate for each factory for an overall creation rate of 50, 000 particles per second. Along with the properties, particle-particle and particle-wall interactions can also be defined under the particle and wall definitions. These properties are visco-elastic damping factor (also called coefficient of restitution), the coefficients of static and rolling frictions (section 2.4.1). The calculation of the coefficient of restitution for various interactions are formulated in section 3.4. Hertz-Mindlin with JKR model (section 2.4.1) is chosen to model the behaviour of cohesive/adhesive materials like snow. This model in EDEM requires interfacial surface energy, also called work of adhesion (Γ). To start with, the work of adhesion values were chosen as 0.218J/m2 for particle-particle interactions and 0.117J/m2 for particle-wall interactions. These values are taken from Eidevåg et al. [7] which are equivalent work of adhesion values to simulate colder snow (-12 ◦C). All particle- particle and particle-wall combinations were assigned with these values. To inves- tigate higher angle of repose values, the work of adhesion was increased to model snow at temperatures greater than -12 ◦C. 20 3. Methodology 3.3.2 Simulator Simulator allows the user to perform the simulation. A user can opt for different settings to run the simulations. Here, the user can set the time-step, total simulation run time, grid sizing, the solver configurations, etc. The time-step used in the simulations are calculated using Rayleigh formula men- tioned in section 2.4.3. It was suggested to use less than 20% of this time-step for the simulations as one of the best practices, to avoid the instabilities (Rayleigh waves) caused by improper defining of the time-step. The simulations were ran for an average time of 1.3 to 1.5 sec, which is when the pile formed on the platform stops growing further. However, some simulations with higher work of adhesion values required more simulation time as they required more particles in the pile. The grid size in EDEM is calculated based on the size of the smallest particle in the simulation. The best practice as mentioned in EDEM user guide [30] is to use 2 to 3 times the radius of the smallest particle (Rmin) in the domain. Using this makes the solver to detect the possible contacts faster. If too coarse grid size is chosen, it is possible that more particles can be in the same grid cell and the solver assumes that all particles in this cell have possible contact scenarios, which can make the simulation run slower. 3.3.3 Analyst Analyst is used to do all the post-processing work in EDEM. Determination of angle of repose, mass of the pile, plotting various parameters like velocities, angular velocities, number of particles etc., can be don in EDEM-analyst. For this work the protractor tool in the analyst is used to measure the angle of repose of the piles as shown in figure 3.6. Three points are marked on the GUI using the coordinate system of the geometries as reference. Figure 3.6: Protractor tool used in EDEM to measure angle of repose of the snow pile. 21 3. Methodology 3.4 Modeling coefficient of restitution for elastic damping The coefficient of restitution given in equation 2.6 is one of the input parameters in EDEM. This is an important input to accurately model the elastic damping during the contact. Higa et al. [10] proposed a correlation for the coefficient of restitution of ice particles based on the equivalent particle radius/diameter. Physical experiments with ice particles were performed to assess the critical (collision) velocities of the particles when they were made to collide with a wall. As explained in [10], each particle size had a different critical velocity after which the coefficient of restitution is calculated using equation 3.1. e = eqe ( vi vc )− log( vi vc ) , if vi ≥ vc (3.1) Here, ’e’ is the coefficient of restitution, eqe is the coefficient of restitution in quasi- elastic region i.e., e = eqe if vi < vc, vi is the impact velocity during collision, vc is the critical velocity based on particle size. Quasi-elastic coefficients of restitution for various particle sizes are given in table 3.5. However, the particle sizes analysed in [10] are much larger than the size distributions taken in the simulations. Therefore, the coefficient of restitution values were extrapolated for the particles used in the simulations. The extrapolation was done using equation 3.2 which was the best fit for the quasi-elastic coefficient of restitution trend in [10]. See figure 3.7 Particle radius (cm) Critical velocity (cm/s) eqe 3.6 22.7 0.95 1.5 40.6 0.89 0.8 54 0.86 0.4 70.2 0.79 0.14 124 0.71 Table 3.5: Critical velocities and quasi-elastic coefficient of restitution for particle sizes in [10]. e = 0.0744 ln(R∗) + 0.8611 (3.2) Here, R∗ is the equivalent particle radius (in cm). 22 3. Methodology Figure 3.7: Extrapolation trend-line for quasi-elastic coefficient of restitution given by equation 3.2. The experiments in Higa et al.[10] are done only for particle-wall interactions. There- fore, it is needed to calculate an equivalent particle radius using equation 2.12 to account for particle-particle interactions. A few collision scenarios (shown in fig- ure 3.8) were taken to see if any particle reaches a velocity higher than the critical velocities calculated by extrapolating the available values using equation 3.3. It was observed that the impact velocities of the particles did not exceed the critical velocities of the corresponding particles up to a particle radius of approximately 400µm. This implies that the relation specified in equation 3.1 may not be used in the current simulations i.e., the quasi-elastic coefficients of restitution are sufficient in this case. Vcritical = 46.232 (R∗)−0.506 (3.3) Here, R∗ is in cm. Figure 3.8: Possible collision scenarios in the simulations. Particle-wall collision (left), particle-particle (same size) collision with one of the particles rebounding after collision with wall (middle), particle-particle (different size) collision with one of the particles rebounding after collision with wall (right). 23 3. Methodology 3.5 Methods to reduce simulation run time 3.5.1 Reducing Young’s modulus One of the methods to reduce simulation run time is to increase the time-step of the simulation. As discussed in section 2.4.3 direct tweaking of time-step leads to numerical instabilities. Hærvig et al. [9] proposed a method of altering Young’s modulus (E) of the particles to increase the time-step as time-step in equation 2.26 depends on Young’s modulus. To compensate for the loss of Young’s modulus the work of adhesion (Γ) is reduced with the relation shown in equation 3.4. Γ∗ = Γ ( E∗ E ) 2 5 (3.4) Here, Γ∗ and Γ are the modified and original work of adhesion values respectively, E∗ and E are the modified and original Young’s modulus values respectively. A few single particle simulations with various sized particles were evaluated to val- idate the above-mentioned reduction. These simulations are done with no visco- elastic damping and with a work of adhesion value of 0.117 J/m2 so that the damp- ing in the simulations is only because of the work of adhesion. The coefficients of restitution were calculated using the rebound velocity after an impact with the wall as shown in 3.9. The rebound velocities are checked both for normal and tangential impacts. The tests with reduced Young’s modulus gave a better match for rebound velocities when compared with the tests with original Young’s modulus (with an insignificant error). Figure 3.9: Impact tests to validate the equation 3.4. Test for normal impact validation (left), tangential impact validation (right). When the Young’s modulus was reduced by 1000 times i.e., E∗ = E/1000 the time step has increased approximately by 100 times (for a particle size of 100µm). This implies that there is a significant potential to increase the time-step using this approach. 24 3. Methodology 3.5.2 Implementing dynamic domain (EDEM) Dynamic domain is an option provided in EDEM. The dynamic domain is a region in which all the particle velocities are considered and the particles with minimal displacement outside this region are considered as frozen i.e., the contact forces on these particles are frozen and not updated anymore. See figure 3.10. A geometry has to be defined as a dynamic domain whose only purpose is to act as a boundary. Figure 3.10: Dynamic domain implementation in the simulations. Implementation of dynamic domain with falling particles, no frozen particles (left), particles frozen outside dynamic domain (right). There are few solver options to implement this method namely, check interval, num- ber of checks and displacement of particle. Check interval is the simulation time interval (in seconds) which checks if the particle moves in the specified interval. If it doesn’t move it is flagged to be as frozen. The number of checks determines the number of times solver checks to see if a particle is frozen. Displacement of particle is given as a percentage of the radius of the smallest particle in the domain. This determines if the particle is stationary. If the displacement of the particle within the check interval is greater than this value then the particle is considered to be moving or else it is frozen and is removed from contact force and contact detection calculations. Avoiding contact force and contact detection calculation for at least few particles in the domain will decrease the simulation run time. Dynamic domain has to be used very carefully in the simulations. If a too low check interval with a higher displacement of particle is specified in the solver options, the particles will freeze rapidly producing incorrect results. In the current work a check interval of 0.01seconds is specified and the displacement of particle setting is defined 25 3. Methodology as 1% of minimum radius particle in the domain i.e., 2.25µm. Number of checks were limited to 3. 26 4 Results and discussion This section describes the simulation results and their comparison with experimental results. The influence of various contact parameters like velocity dependent visco- elastic damping coefficient, friction coefficients, work of adhesion on angle of repose are discussed. The influence of fall height on the angle of repose is also discussed. As discussed in section 3.3, the simulations were setup with different work of adhe- sion values starting with 0.218 J/m2. The coefficients of restitution were calculated using the method described in section 3.4 and were used as simulation inputs. The coefficient rolling friction is taken as 0.1. The sliding/static friction value was taken as prescribed in [28], which was 0.1. The initial simulation settings are shown in table 4.1. Parameter Value Young’s Modulus (particle) [GPa] 9.7 Young’s Modulus (wall) [GPa] 2 Work of Adhesion (particle-particle) [J/m2] 0.218 Work of Adhesion (particle-wall) [J/m2] 0.117 Rolling Friction Coefficient [-] 0.1 Sliding Friction Coefficient [-] 0.1 Table 4.1: Initial simulation settings in EDEM. 4.1 Effect of rolling resistance on angle of repose The effect of rolling resistance should be evaluated to correctly determine the angle of repose. A few simulations were performed by varying rolling friction coefficient to know the importance of rolling resistance. Figure 4.1 shows the effect of rolling resistance on angle of repose. The trend was made by taking a constant work of adhesion value of 0.218 J/m2. The trend in figure 4.1 is most likely because the rolling resistance is getting saturated with the increased in rolling friction coefficient and the particles that move tangentially start to slide instead of rolling. It can be seen that the angle of repose has increased with increase in rolling resistance initially and has settled at a constant value after rolling friction coefficient became 0.1. This motivated to continue the simulations with a rolling friction coefficient of 0.1. A simulation with higher rolling resistance (µr = 0.9) is also performed for each case with higher work of adhesion value(s) shown in table 4.2, but there was no change in angle of repose observed. 27 4. Results and discussion Figure 4.1: Effect of rolling resistance on angle of repose. Work of adhesion is kept constant i.e., 0.218 J/m2. 4.2 Effect of visco-elastic damping factor For all the simulations, the velocity dependent visco-elastic damping coefficient was modelled using the method described in section 3.4. The coefficients were assigned for each particle-particle and particle-wall interaction scenario. It is interesting to see the change in angle of repose when this value is changed. A test case was performed with an increased coefficient of restitution i.e., 0.9, for a simulation with three sphere particles (figure 4.8). It was compared with the simulation having a coefficient of restitution of 0.6 for particle-wall interaction and 0.548 for particle- particle interaction. The angle of repose results from both the simulations are shown in figure 4.2. 28 4. Results and discussion Figure 4.2: Effect of velocity dependent visco-elastic damping on angle of repose. Angle of repose using modelled damping coefficients in section 3.4 (left), angle of repose with a higher damping coefficient of 0.9 (right). It can be seen that when the damping was decreased (i.e., damping coefficient was increased) by a large extent, the effect on the angle of repose is very small. The angle of repose only decreased by 6.5◦ when the damping was decreased by 0.5 times. This implies that velocity dependent visco-elastic damping has a very little influence on the angle of repose. 4.3 Angle of repose for reduced Young’s modulus A method was proposed in section 3.5.1 to increase the time-step in the simulation, which in-turn reduces the simulation run time. To validate this method for angle of repose simulations, a simulation with original Young’s modulus value and a simula- tion with reduced Young’s modulus value are analysed. The angle of repose values of both the cases can be seen in figure 4.3. Figure 4.3: Angle of repose for reduced Young’s modulus. Simulation with original Young’s modulus (left), reduced Young’s modulus [E∗ = E/100] (right) . 29 4. Results and discussion It is evident that the simulations with reduced Young’s modulus were not suitable to simulate bulk simulations like angle of repose. Even though the normal and tangential interactions were validated for reduced Young’s modulus simulations, the sliding and rolling friction models were not verified. The bulk densities of these simulations match with the bulk density values measured for piles with original Young’s modulus. As explained in [9], the rolling resistance model should also be altered to compensate for the changed Young’s modulus and work of adhesion. The rolling resistance models were not modified in this work and the simulations are continued without changing the Young’s modulus and work of adhesion values as mentioned in section 3.5.1. 4.4 Effect of sliding resistance on angle of repose The effect of sliding resistance on angle of repose was evaluated using two work of adhesion values i.e., 0 and 0.218 J/m2 to see how angle of repose changes by changing sliding resistance. The influence of sliding resistance on angle of repose is shown in figure 4.4. Figure 4.4: Effect of sliding resistance on angle of repose for two work of adhesion values i.e., 0 J/m2 (blue) and 0.218 J/m2 (black). 30 4. Results and discussion We can see that the angle of repose increases with the increase in sliding resistance. However, these values of angle of repose shown in figure 4.4 with high sliding resis- tance seem to be too high and non-realistic for zero work of adhesion case. Therefore, the value of 0.1 was chosen as sliding friction coefficient as mentioned in [28]. 4.5 Effect of fall height on angle of repose To know the effect of fall height on the angle of repose, three simulations with fall heights of 30mm, 60mm, 100mm were performed and compared with the experi- mental results at various fall heights. Figure 4.5: Angle of repose dependency on fall height. Experimental results at -17 ◦C for faceted snow (red), simulation results for spherical ice particles (blue). The trend of the angle of repose with varying fall height value in simulations and experiments is shown in figure 4.5. It can be seen that both simulation and ex- perimental trends are decreasing almost linearly with increasing fall height. This happens because the particles in experiments/simulations with higher fall heights will attain higher velocities before colliding with platform/pile making them to build up a smaller pile and vice-versa happens at lower fall heights. The difference in val- ues of angle of repose between experiments and simulations at same fall height can 31 4. Results and discussion be motivated using several reasons. The snow type considered in simulations was rounded/spherical and in experiments it was considered as faceted. Another reason could be that the modelled sliding friction coefficient in simulations (µs = 0.1) is too low. As it was seen in figure 4.4 that the angle of repose increases with an increase in sliding resistance, the sliding resistance can be increased to exactly match the values on simulation and experimental trends. 4.6 Variation of angle of repose with work of ad- hesion The simulations were setup with an initial work of adhesion value of 0.218 J/m2 and the settings shown in table 4.1. The work of adhesion value is gradually increased to see the change in the angle of repose values. The work of adhesion values adopted in this work are presented in table 4.2. The coefficients of restitution were calculated using the method described in section 3.4 and were used as simulation inputs. The coefficient rolling friction is taken as 0.1, the choice of this value is motivated in section 4.1. The sliding/static friction value was taken as prescribed in [28], which was 0.1. The angle of repose piles with these settings are shown in figure 4.6. Figure 4.6: Variation of angle of repose with change in work of adhesion values.The angle of repose increases with an increase in work of adhesion. 32 4. Results and discussion The angle of repose was determined by taking the slope of particles formed at the surface of the pile not by measuring from the tip of the pile. Number Work of Adhesion (J/m2) Angle of Repose 1 0.218 37.2◦ 2 0.327 40◦ 3 0.436 49.21◦ 4 0.654 62.9◦ 5 0.872 90◦ Table 4.2: Values of work of adhesion chosen for the simulations and the angle of repose values from them. The simulation trend of angle of repose versus work of adhesion was compared with the experimental trend, which was between angle of repose and temperature of snow. Both trends are shown in figure 4.7. The comparison was plotted by equating the work of adhesion scale to the temperature scale at 0.218 J/m2 and -12 ◦C respectively. This was taken according to the information described in [7] and section 3.3. The lines in figure 4.7 show the experimental trend for a fall height of 100mm and the simulation trends for two fall heights i.e., 30mm and 100mm. From figure 4.7 the trend of both data i.e., angle of repose versus temperature and angle of repose versus work of adhesion at 100mm fall height seem to agree. Most of the simulations are done at 30mm fall height as explained in section 3.3. The diameter of platforms in simulation and experiments are also different taken as 25mm and 50mm respectively. The platform diameter was reduced in the simulations as the 50mm diameter platform requires more particles to form a pile when compared to the 25mm diameter platform. Simulating more particles requires more computational effort, in order to reduce the computational effort a lower diameter was chosen instead. Even though the simulation and experimental trends are similar, both the curves at similar fall height do not have the exact same values. The difference can be reasoned in several ways. As explained in section 4.5, the difference could be because of insufficiently modeled sliding resistance. The particles taken in experiments have different shapes and might have different snow types (faceted particles, depth hoar, etc.), but in simulations only spherical particles were considered. The bulk densities of the piles were measured by creating a box with dimensions 5mm × 5mm × 5mm inside the pile. The mass and volume of the box are then used to evaluate the bulk density, which varied in between 545 kg/m3 and 560 kg/m3. These values are close to the values from the experiments. To check the dependency of angle of repose on non-sphericity, an equivalent particle consisting is three spheres adopted from [11], shown in figure 4.8 was chosen to represent the entire particle distribution with non-sphericity and the simulations are done for the work of adhesion values mentioned in 4.2. Choosing a non-spherical particle will also indirectly change the rolling resistance value. 33 4. Results and discussion Figure 4.7: Comparison of experimental and simulation trends for angle of repose versus temperature (red) and work of adhesion (blue, black) respectively. The figure also shows error bars while measuring angle of repose from both experiments and simulations. Figure 4.8: Three spherical particle used in simulations to evaluate the effect of non-sphericity. The model was adopted form [11], where sintered particles were represented as 3 sphere particles. 34 4. Results and discussion It can be seen that the angle of repose values in simulations with the three sphere particles exactly matched with the trend of spherical particles with a minor difference at two points. This once again confirms that increasing rolling resistance by changing particle shape did not affect the angle of repose values. The trends of angle of repose versus work of adhesion for spherical and non-spherical particles can be seen in figure 4.9. Figure 4.9: Comparison of trends of angle of repose versus work of adhesion for spherical particles (red) and equivalent non-spherical particles (black). 35 4. Results and discussion 36 5 Conclusion The objective of this work was to correlate angle of repose of snow at different temperatures to the work of adhesion parameter described in JKR model. It was identified that the rolling and sliding resistance also influence the angle of repose values to a certain extent. It was observed that the simulation and experimental trends are matching. An exact match between data on the trend lines can be achieved by methods suggested in the previous sections i.e., increasing the sliding resistance of the particles. A similar simulation study can be made with the particles with different snow type like faceted snow or with a different non-spherical particle size distributions to see how these particles affect the angle of repose values. Alternate rolling and sliding resistance models which were prescribed in section 2.4.1 can be explored to see if adopting these models can change the angle of repose val- ues determined in the current study. An extrapolation method is used to model the velocity dependent visco-elastic damping, which was assumed to be accurate to model the coefficients of restitution for all the interactions. Also, the experiments in [10] are performed for particle-wall interactions, where the equivalent radius in all the interactions is the particle radius. For the particle-particle interactions the equivalent particle radius is taken as the harmonic mean of corresponding two par- ticles in contact. Further investigations are needed to validate the accuracy of the modeled visco-elastic damping in the current work. A method proposed to increase the time-step in the simulations to reduce the total simulation run time. This method worked well for single particle interactions in normal and tangential directions, but didn’t work well for bulk simulations i.e., the angle of repose simulations. However, the single particle simulations were not validated for proper working of rolling and sliding resistance models as there was limited project time available. The sliding and rolling resistance models can be checked for working correctly and necessary changes can be made if necessary to compensate for the increased time-step. The fall height trend seen in the simulations was same as in the experiments. The sliding resistance can be increased in the simulations to match the exact values of angle of repose in experiments. There can also be some other factors that influence the difference in angle of repose values for both the trends which can be transverse forces (wind currents) acting on the particles, change in velocities of the particles due to vibrations of the hopper. The purpose of this work was to predict the snow impingement/accumulation on the autonomous sensors of a car. The work of adhesion results from this work can be used to mimic the snow at different temperatures. These values of work of adhesion can 37 5. Conclusion therefore be used as an input in DEM-aerodynamic simulations with one or two way coupling to model the right degree of adhesion for the snow impingement. They can also be used in Lagrangian-Eulerian aerodynamic simulations, by using user-defined functions, which can have conditional statements of when the particle sticks. This can be done by correlating the velocities of particles achieved in angle of repose simulations with various fall heights and assigning the degree of adhesion/stickiness for each velocity value. 38 Bibliography [1] Kuroiwa D, Mizuno Y, Takeuchi M. Micromeritical properties of snow. Physics of Snow and Ice: proceedings. 1967;1(2):751-72. [2] Furukawa, Y., 2007. Snow and ice crystals. Physics today, 60(12), pp.70-71. [3] Hagenmuller, P., Chambon, G., Flin, F., Morin, S. and Naaim, M., 2014. Snow as a granular material: assessment of a new grain segmentation algorithm. Granular Matter, 16(4), pp.421-432. [4] Cundall, P.A. and Strack, O.D., 1979. A discrete numerical model for granular assemblies. geotechnique, 29(1), pp.47-65. [5] Marshall, J.S. and Li, S., 2014. Adhesive particle flow. Cambridge University Press. [6] Thornton, C., 2015. Granular dynamics, contact mechanics and particle system simulations. Springer. [7] Eidevåg, T., Abrahamsson, P., Eng, M. and Rasmuson, A., 2020. Modeling of dry snow adhesion during normal impact with surfaces. Powder Technology, 361, pp.1081-1092. [8] Abrahamsson, P., Eng, M. and Rasmuson, A., 2018. An infield study of road snow properties related to snow-car adhesion and snow smoke. Cold regions science and technology, 145, pp.32-39. [9] Hærvig, J., Kleinhans, U., Wieland, C., Spliethoff, H., Jensen, A.L., Sørensen, K. and Condra, T.J., 2017. On the adhesive JKR contact and rolling models for reduced particle stiffness discrete element simulations. Powder Technology, 319, pp.472-482. [10] Higa, M., Arakawa, M. and Maeno, N., 1998. Size dependence of restitution coefficients of ice in relation to collision strength. Icarus, 133(2), pp.310-320. [11] Colbeck, S.C., 1982. An overview of seasonal snow metamorphism. Reviews of Geophysics, 20(1), pp.45-61. [12] Colbeck, S.C., 1998. Sintering in a dry snow cover. Journal of Applied Physics, 84(8), pp.4585-4589. 39 Bibliography [13] Blackford, J.R., 2007. Sintering and microstructure of ice: a review. Journal of Physics D: Applied Physics, 40(21), p.R355. [14] Langlois, A., Royer, A., Montpetit, B., Roy, A. and Durocher, M., 2020. Pre- senting snow grain size and shape distributions in northern Canada using a new photographic device allowing 2D and 3D representation of snow grains. Frontiers in Earth Science, 7, p.347. [15] Ai, J., Chen, J.F., Rotter, J.M. and Ooi, J.Y., 2011. Assessment of rolling resistance models in discrete element simulations. Powder Technology, 206(3), pp.269-282. [16] Israelachvili, J.N., 2011. Intermolecular and surface forces. Academic press. [17] Gaylard, A.P., Pitman, J., Jilesen, J., Gagliardi, A., Duncan, B., Wanderer, J. and Konstantinov, A., 2014. Insights into rear surface contamination us- ing simulation of road spray and aerodynamics. SAE International Journal of Passenger Cars-Mechanical Systems, 7(2014-01-0610), pp.673-681. [18] Hobbs, P.V., 2010. Ice physics. Oxford university press. [19] Kloubek, J., 1974. Calculation of surface free energy components of ice accord- ing to its wettability by water, chlorobenzene, and carbon disulfide. Journal of Colloid and Interface Science, 46(2), pp.185-190. [20] Boinovich, L.B. and Emelyanenko, A.M., 2014, December. Experimental de- termination of the surface energy of polycrystalline ice. In Doklady Physical Chemistry (Vol. 459, No. 2, pp. 198-202). Pleiades Publishing, Ltd. [21] Gold, L.W., 1958. Some observations on the dependence of strain on stress for ice. Canadian Journal of Physics, 36(10), pp.1265-1275. [22] Parker, E.R., 1967. MATERIALS DATA BOOK. For Engineers and Scientists. [23] Professional Plastics, Mechanical Properties of Plastic Materials, [Online]. Available: https://www.professionalplastics.com/professionalplastics/MechanicalPropertiesofPlastics.pdf [24] Van Krevelen, D.W. and Te Nijenhuis, K., 2009. Properties of polymers: their correlation with chemical structure; their numerical estimation and prediction from additive group contributions. Elsevier. [25] Murray, B.J., O’sullivan, D., Atkinson, J.D. and Webb, M.E., 2012. Ice nucle- ation by particles immersed in supercooled cloud droplets. Chemical Society Reviews, 41(19), pp.6519-6554. 40 Bibliography [26] Nakaya, U. and Matsumoto, A., 1954. Simple experiment showing the existence of “liquid water” film on the ice surface. Journal of Colloid Science, 9(1), pp.41- 49. [27] Miller, D.A., Adams, E.E. and Brown, R.L., 2003. A microstructural approach to predict dry snow metamorphism in generalized thermal conditions. Cold regions science and technology, 37(3), pp.213-226. [28] Braghin, F., Cheli, F., Maldifassi, S., Melzi, S. and Sabbioni, E. eds., 2016. The engineering approach to winter sports. New York, NY: Springer. [29] Tamadondar, M.R., de Martín, L., Thalberg, K., Björn, I.N. and Rasmuson, A., 2018. The influence of particle interfacial energies and mixing energy on the mixture quality of the dry-coating process. Powder Technology, 338, pp.313- 324. [30] https://www.edemsimulation.com/edem-2020-0-documentation/ 41 Bibliography 42 List of Figures List of Tables Introduction Background Purpose & Objectives Assumptions and limitations Theory Angle of Repose Snow Characteristics Road/Untouched/Ground snow Discrete Element Method Hard sphere approach Soft sphere approach Governing equations in DEM Contact Forces Other resistance models Body (non-contact) Forces Time-step EDEM Overview Contact radius (or radius of contact) in EDEM Methodology Experimental studies Experimental setup Microscopic analysis Other tests Particle size distribution Simulation setup Creator Simulator Analyst Modeling coefficient of restitution for elastic damping Methods to reduce simulation run time Reducing Young's modulus Implementing dynamic domain (EDEM) Results and discussion Effect of rolling resistance on angle of repose Effect of visco-elastic damping factor Angle of repose for reduced Young's modulus Effect of sliding resistance on angle of repose Effect of fall height on angle of repose Variation of angle of repose with work of adhesion Conclusion Bibliography