A bonded discrete element approach to simulate loading with hydraulic mining excavators

(1)

(2)

Preface

This master thesis is the last and final course within the Master Programme in Mechanical Engi-neering at Luleå University of Technology. The master thesis was going on between January 2021 and June 2021 and the work was carried out remotely in collaboration with Boliden.

Boliden has contributed with a lot of necessary material so that I was able to carry out the master thesis remotely, which I greatly appreciate. I have been accompanied by good and nice employees from Boliden and would like to thank them for their involvement in this work. I would like to give a big thanks to supervisor Andreas Svanberg. He has supported me throughout the master thesis and has always been available for questions. Also, thanks to Mattias Linnanki, the responsible contact person at Boliden, who helped me through the project by contributing to discussions that led the project forward. I would also like to give thanks to Simon Larsson for the same reason and for assisting me throughout the project as the examiner. I would also like to thank Jonas Edberg who has helped me with any questions related to the computer cluster. Finally, I would like to thank my friends and relatives who have supported me during my education at Luleå University of Technology.

(3)

Abstract

When operating hydraulic mining excavators the loading equipment is exposed to harsh conditions which lead to extensive wear of the equipment, especially the bucket and bucket teeth. Simulations are used to better understand the wear development and to evaluate new methods to operate excavators more efficiently. At the Aitik mine, operated by the high-tech metal company Boliden Mines, hydraulic excavators are used when loading the mined ore. One of the hydraulic excavators used at Aitik is the Komatsu PC7000.

In this master thesis, a simulation model for the hydraulic excavator Komatsu PC7000 was de-veloped with the simulation software LS-DYNA. This model consists of multi rigid body dynamics to describe the motion of the excavator and a granular material model to describe the rocks loaded into the bucket of the excavator. Simulations with two different types of granular material models have been utilized to study the wear development of the bucket. One of the models (bonded DE model) uses bonded discrete elements to describe the large rocks and single discrete elements are used to describe smaller rocks. This model is compared to the current FE-DE model which is being used today at Boliden. This model uses finite elements (FE) to model the larger rocks and discrete element spheres (DES) for smaller rocks. By using the bonded DE method a 71% reduction in simulation time could be achieved. This can be partly explained by the reduction of the number of elements included in the rock pile.

Archard’s wear law was used to numerically describe the wear development of the bucket. When simulating the wear a total of 30 bucket fillings were performed with the excavator. This was done with both the bonded DE method and the FE-DE method. In this wear study, the inside of the bucket was of interest. The resulting simulated wear map was compared to experimental measurements from which the plate thickness of the bucket had been measured two times to obtain the wear depth of some points inside the bucket. The experimental measurements and two 3D scanned point clouds were used to determine the wear depth inside the bucket. Results from the simulation showed that the wear is concentrated to the center of the bucket while less wear is concentrated to the sides of the bucket. With the bonded DE method the wear appeared to be more evenly distributed inside the bucket while the wear from the FE-DE method appeared in spots inside the bucket. The experimental results also showed that the wear was more extensive in the center of the bucket and also in the back of the bucket. Both simulation methods also showed that the wear was concentrated to the back of the bucket. From the simulations, it was also seen that the behavior of the material flow differed between the two methods. In the bonded DE method the material flow had more sliding behavior while the material flow in the FE-DE method had more rolling behavior. This could also be the reason why the bonded DE method captures the wear more evenly. The rolling behavior seen in the FE-DE method leads to more impact wear which is not captured by Archard’s wear law.

Overall, the bonded DE method leads to a big reduction in simulation time which is favorable when it comes to simulation. The larger rocks will have simpler shapes without sharp corners. However, the method allows for a more complex shape than just an ordinary sphere which is the simplest and most common shape to describe granular material. The bonded DE method also allows for easier configuration of contact definition since fewer contact interfaces must be added to the model. Furthermore, the post-processing of wear in LS-DYNA was facilitated since the wear does not have to be divided into two wear collectors for FE elements and DE elements.

(4)

Sammanfattning

När hydrauliska grävmaskiner används i gruvindustrin utsätts lastningsutrustningen för stora påfrestningar på grund av tuffa förhållanden vilket leder till slitage på utrustningen, särskilt på skopa och skoptänder. Simuleringar används för att bättre förstå utvecklingen av slitaget och för att utvärdera nya metoder som ska göra att grävmaskinerna kan manövreras mer effektivt. I Aitik-gruvan, som drivs av det högteknologiska metallföretaget Boliden Gruvor, används bland annat hydrauliska grävmaskiner när den brutna malmen ska lastas i gruvtruckar. En av de hydrauliska grävmaskiner som används i Aitik är Komatsu PC7000.

I detta examensarbete har en simuleringsmodell utvecklats för den hydrauliska grävmaskinen Komatsu PC7000. Simuleringsmodellen utvecklades med simuleringsprogrammet LS-DYNA. Mod-ellen använder sig av stelkroppsdynamik för att beskriva rörelsen av grävmaskinen och en gran-ulär materialmodell för att beskriva malmen som lastas i skopan på grävmaskinen. Simuleringar har genomförts med två olika granulära materialmodeller för att studera utvecklingen av slitaget på skopan som uppstår på grund av nötning från det granulära materialet. En av modellerna, klustrade DE modellen, använder sig av klustrade diskreta element för att beskriva de större ste-narna och enstaka diskreta element för att beskriva de mindre steste-narna. Denna modell jämförs med den nuvarande FE-DE modellen som används för simuleringar vid Boliden. Denna modell använder sig av finita element (FE) för att beskriva de större stenarna och enstaka diskreta element för att beskriva de mindre stenarna. Genom att använda modellen med klustrade DE kan simu-leringstiden minskas med ungefär 71% vilket bland annat kan förklaras av det minskade antalet element som ingår i stenhögen.

Archard’s nötningslag användes för att beskriva nötningen av skopan numeriskt. Totalt 30 lastningscykler simulerades med både klustrade DE metoden och med FE-DE metoden. Skopans insida var av intresse då nötningen simulerades. Den simulerade nötningskartan jämfördes mot experimentella mätningar där skopans plåttjocklek hade mätts upp två gånger för att bestämma nötningen av ett antal punkter i skopan. Mätningarna och två 3D skannade punktmoln av skopan användes för att bestämma nötningen inuti skopan. Resultatet från simuleringarna visade att nöt-ningen var koncentrerad till mitten av skopan och att nötnöt-ningen var mindre på sidorna i skopan. Med den klustrade DE modellen visade sig nötningen vara mer jämnt utbredd inuti skopan medan för FE-DE modellen visade sig nötningen uppstå punktvis inuti skopan. Det experimentella re-sultatet visade också att nötningen var större i mitten av skopan och även i den bakre delen av skopan. Båda simuleringsmodellerna visade att nötningen var koncentrerad till den bakre delen av skopan. Det visade sig också att flödet av materialet in i skopan betedde sig på två olika sätt för de två metoderna. I klustrade DE modellen gled materialet in i skopan medan materialet rullade in i skopan för FE-DE modellen. Det här kan också vara orsaken till att den klustrade DE modellen fångade nötningen mer jämnt. Det beror på att FE-DE modellen bidrar med nötning från stötar vilket inte fångas upp av Archard’s nötningslag som bara fångar nötning från glid.

Sammanfattningsvis visar resultaten att den klustrade DE modellen leder till en stor minskning av simuleringstiden vilket är en stor fördel när det kommer till simuleringar. De större stenarna får däremot en enklare form utan skarpa hörn men de får fortfarande en mer komplex form än formen av en sfär som är den enklaste och vanligaste formen för att beskriva granulära material. Med klustrade DE modellen är det också enklare att ställa in hur kontakt sker mellan partiklar och ytor då det är färre kontakter som måste definieras. Med klustrade DE modellen får man också enklare efterbehandling av nötning då all nötning fångas upp av en och samma behållare i simuleringsprogrammet LS-DYNA.

(5)

1 Introduction

In today’s society, there is an ongoing transition to a world powered by renewable energy. This transition will require the necessary metals used in today’s innovation such as electric vehicles and wind turbines which is the key to a sustainable world.

At Aitik, an area situated outside the town of Gällivare in northern Sweden, ore is mined in the world’s most efficient open-pit copper mine. Aitik is one of currently 5 mines operated by the high-tech metal company Boliden Mines. In 2019 a total of 40 million tonnes of ore were mined at Aitik to provide the necessary metals for a sustainable world. An overview of the Aitik mine can be seen in Figure 1.

Figure 1: Aitik open pit mine

1.1 Background

This section will give a brief introduction of the underlying areas of the thesis such as open-pit mine environment, mining equipment and simulation methods for mining equipment.

1.1.1 Open-pit mine environment

Aitik is the world’s most efficient copper mine and this requires an optimal production process. Some parts of the production process are shown in Figure 2.

(8)

Figure 2: Overview of production at Aitik [3]

The process begins by drilling and blasting to break loose the ore. At Aitik blasting is done approximately one to two evenings per week. Ore is then loaded onto rock dumpers by using a combination of loading equipment such as rope shovel excavators, hydraulic excavators and wheel loaders. The ore is then transported with rock dumpers to the crusher. 100 000 tonnes of ore are processed per day at the crushers. The crushed ore is then transported to a storage using a conveyor. At the concentrator, the crushed ore is ground into sand and the metal is separated from the sand in the flotation process. With the dewatering and filtering process, the metal concentrate is separated from the sand. The metal concentrate is then transported to Rönnskär smelter using railroads. The sand is stored at a nearby sand magazine.

This master thesis will be based on the loading part of the production process where the blasted ore is loaded onto the rock dumpers. This part of the production is of high interest because if the loading stops then the whole production stops. Therefore it is valuable to reduce any downtime of the equipment used in the loading. For example, by creating a good schedule for the maintenance of the equipment. Some maintenance includes replacing worn bucket teeth of the excavators used in the loading process.

1.1.2 Mining equipment

In Aitik eight loaders are used when loading the mined ore into rock dumpers. Out of these, there are two rope shovel excavators, three hydraulic excavators and two wheel loaders.

One of the loaders used at the Aitik mine is the Komatsu PC7000. It is a hydraulic excavator with an operating weight ranging from 676-694 tonnes. Figure 3 shows the Komatsu PC7000 next to the dumper Caterpillar 795F at the Aitik mine.

(9)

Figure 3: Komatsu PC7000

A sketch of the Komatsu PC7000 with its moving components can be seen in Figure 4. The components include the upper body (2) that can rotate around the crawlers (1) with a revolute joint. The boom (3) is connected to the upper body with a revolute joint and the two piston and cylinder pairs (7-8) acts as a translational motor between the upper body and the boom. The stick (4) is connected to the boom with a revolute joint and the two piston and cylinder pairs (9-10) acts as a translational motor between the stick and the boom. The rear bucket (5) can rotate with a revolute joint around the stick and the rear bucket cylinder (12) respectively. Two piston and cylinder pairs (11-12) acts between the rear bucket and the boom as a translational motor. The two piston and cylinder pairs (13-14) holds the front bucket (6) in contact with the rear bucket and controls the opening and closing of the bucket. The front bucket rotates with a revolute joint around the rear bucket.

Figure 4: Components of the Komatsu PC7000 with front shovel attachment

The Komatsu PC7000 operates with a front shovel attachment with a bucket that has a capacity of loading 36 m3 _{rocks or a maximum weight of 68 tonnes (when taking into account the bulk}

density of 1.9 tonnes/m3_{). When in operation it takes an average of five bucket fillings to load}

(10)

Aitik has been followed and studied during a period of approximately 87 days. For this period the total load was 1289,4 ktonnes which is 14,9 ktonnes/day or approximately 300 fillings per day (taking into account the 48 tonnes of rock per filling). However, the excavator was tracked during its run-in period and therefore it was a lot of downtime of the excavator, for example, because of educational purposes. The master thesis will study the loading process of the Komatsu PC7000. 1.1.3 Simulation of mining equipment

The buckets used by the excavators are exposed to harsh conditions during the loading process. Not only due to the heavy amount of rocks carried by the buckets but also because of the many loading cycles performed per day. Such conditions lead to extensive wear of the equipment and therefore also an economical cost for Boliden. For example, through the downtime of the production because of maintenance of the equipment.

Simulations and experimental studies are performed at Boliden to understand the wear de-velopment better and to evaluate new methods to operate the excavators more efficiently during loading. The end goal is to have a schedule for when to replace parts that are exposed to extensive wear to prevent downtime in the production.

Simulation models at Boliden

Simulation models are developed by using commercial simulation software such as LS-DYNA. By using LS-DYNA a simulation model which simulates the complete loading process has been developed and is presented in Svanberg et al. [1]. This model simulates the loading process using the rope shovel excavator Bucyrus 495. The mined ore was simulated using finite elements and discrete elements (FE-DE) to model the granular material and the excavator was simulated with multi rigid body dynamics. This method is later referred to as the FE-DE method. FE was used to model the larger ore and DE was used to model the smaller ore. One bucket used by Bucyrus 495 has also been simulated using rigid FE elements for the bucket and discrete element spheres (DES) for the mined ore to evaluate the wear on the bucket teeth. This was done in Choudry [8] by using Archard’s wear law to numerically model the wear and the results from the simulation were then compared with experimental data of the wear which was determined from 3D scanning. Simulation models for granular material

Some applications of simulations of mining equipment require a granular material model to describe the mined ore. The granular material model can then be used together with the mining equipment such as the hydraulic excavator Komatsu PC7000 to simulate the loading process. From these simulations, it is possible to study the wear development of the bucket used in the loading process. These kinds of simulations could also be used in for example the crushing and the grinding process part of the production at Aitik which also would include a granular material model. Granular material models can also be used in the agricultural industry, for example, to describe granular flow in silos or grain hopper trailers. These models can also be used in the food industry or packaging industry to determine how grain typed content can be packaged efficiently.

One of the difficulties when it comes to simulations is to choose a model that represents reality well and a model that is also computationally efficient. There will always be pros and cons with the model chosen. If the model for the granular material was chosen as the exact representation of reality it would include a lot of small particles with complex shapes, however, it would lead to a computationally inefficient model. Therefore some trade-offs must be done when choosing the model.

(11)

One of the simulation models used for granular material today is the FE-DE method. By using this method it is possible to capture the complex and sharp shapes of larger rocks with FE. Smaller rocks are modelled with DES which just has a center point and a radius. By only using DES it is not possible to model sharp corners since the particle will always be spherical. The advantage however is the computational efficiency. By adding rolling friction or modelling shapes through continuous super quadratic functions this disadvantage can be reduced as discussed in Svanberg et al. [1]. However, it cannot alone model sharp corners.

By using a method referred to as clustering of particles or bonded DE the computational efficient spherical shape can be maintained while also achieving a more complex geometry. Several DES can be combined and form one single rock. However, to obtain sharp corners a very small radius is required and thus many particles to form a rock.

Clustering of particles/bonded DE method was used in Santos et al. [7] where the flow pattern was studied when simulating particles combined of one sphere, two spheres, three spheres or four spheres. In the case with two, three or four spheres clustering of particles was used. Bonded DE is also used to simulate granular material such as concrete. In Tannu [12] bonded DE was used to simulate the uniaxial compressive response of a concrete cylinder.

DEM can also be used together with polyhedral shapes by defining a number of flat polygonal faces and therefore obtaining a geometry with sharp corners but reducing the computational effi-ciency by not using spheres. This was evaluated in Smeets et al. [2] by using cubes to model the particles.

Simulation with LS-DYNA

During the master thesis, the simulation software LS-DYNA was used. This text will give a brief introduction to the LS-DYNA simulation tool. LS-DYNA is one of the leading tools used for simulation and is used in many industries such as automobile, aerospace and manufacturing [4]. The tool includes many capabilities to simulate for example nonlinear dynamics, rigid body dynamics and discrete element method among plenty of other applications. In LS-DYNA a simulation model is built up using keywords. By including the keywords in the model the user can specify things such as materials for parts, the motion of parts, contact between parts and which element type to represent the geometry of the simulated parts. These keywords are referenced with a starting * sign later in the work.

1.2 Aim and objective

The current simulation tool used by Boliden is today only considering the rope shovel excavator Bucyrus 495. However, at Aitik there is a total of 3 hydraulic excavators used during the loading process. Therefore, it is valuable to Boliden to develop a model simulating the loading process with a hydraulic excavator which is the topic of this thesis. Komatsu PC7000 is a hydraulic excavator used in Aitik and will be used as a reference for the project. This excavator was equipped with a front shovel attachment.

The main aim of this thesis is in the end to have a simulation model that simulates the loading process using hydraulic excavator Komatsu PC7000. The created simulation model can then be used together with numerical models to simulate wear to plan which bucket and ground engaging tools (GET) are the most suitable on-site in Aitik. The model could be used to optimize the digging procedure by studying the digging technique to see how the excavator should be operated

(12)

to minimize the wear development of, for example, the bucket teeth. The model could also be used to determine how to re-design the bucket or bucket teeth to reduce the wear development.

The first objective is to develop a multi rigid body simulation that simulates the motion of the hydraulic excavator Komatsu PC7000. The second objective of the master thesis is to develop a different method to simulate the granular material used when simulating the loading process with the Komatsu PC7000. This method will consider bonded DE to simulate the larger mined ore rather than using FE as in Svanberg et al. [1]. By utilizing a bonded DE method instead of FE a possible reduction in computation time could be achieved. A reduction in computational time is highly appreciated for all simulations, especially when used in industrial applications. The bonded DE method will be evaluated by comparison with the FE-DE method in terms of computational time, material flow and achieved wear pattern of the bucket. The third objective is to implement Archard’s wear law to numerically model the wear.

1.3 Problem description

The project will consist of the following:

1. From 3D scanning the geometry of the hydraulic excavator PC7000 will be retrieved together with the dimensions of important parts of the excavator. This information will be used to reverse engineer a CAD model for Komatsu PC7000 which will be used for a multi rigid body simulation in LS-DYNA.

2. The developed multi rigid body simulation will be used together with a FE-DE model for the granular material when simulating the loading process.

3. A new model for the granular material will be developed with bonded DE and will then be used with the multi rigid body simulation in order to simulate the loading process.

4. Both granular material models will be used together with Archard’s wear law in order to study how well the models capture the behavior of wear. This will be compared to the wear that has been experimentally measured on the bucket of the Komatsu PC7000.

5. A video recording of when the Komatsu PC7000 operates in Aitik will be used as valida-tion. This will be used to validate the motion of the excavator. During the validation, the mass loaded into the bucket will be measured and will be used to validate the mass in the simulation.

6. The two methods, FE-DE and bonded DE, will be evaluated and compared to each other but also to the experimental results.

1.4 Limitations

This section will describe the limitations of the thesis.

• The motion of the excavator in the simulation was limited to an estimated motion of how hydraulic excavators operate in reality. The validation was limited to a different motion than the motion used in the simulation.

• The granular material model used for simulation is limited to particles with a radius greater than 100 mm. If particles with a radius less than 100 mm would be used then the time step must be reduced as well as increasing the total number of elements in the rock pile which would lead to a too computationally expensive model.

(13)

• The bonded DE representing larger rocks in the bonded DE method was used with quite a few particles to represent the original FE rocks and was also limited to spheres with a radius greater than 100 mm to keep the same time step as in the FE-DE method.

• Archard’s wear law was used to numerically describe the wear. This law is easily implemented in LS-DYNA however it does not account for impact wear. Calibration is needed for the usage of wear laws. However, no calibration was made and instead the wear constant was used from an in-situ calibrated user-defined wear routine that takes both impact and sliding wear into account. When calibrating this wear routine similar element size was used for the particle representation of the rocks as well as the size of the finite element representation of the bucket. It was assumed that the wear constant from this wear routine is a good approximation for a wear routine that only uses Archard’s wear law. However, the wear constant does not matter since the wear results are normalized when making the comparison between simulation results and experimental results. The wear constant would have to be re-calibrated to be able to compare the magnitude of the wear.

• The validation of the simulated wear was performed against experimental data after loading approximately 1300 ktonnes of rocks. However, the simulation was limited to 30 loading cycles and a total of 1400 tonnes of rocks.

• It was not possible to validate the material flow into the bucket since the excavator was located in an inaccessible corner during the validation.

• The wear study has been limited to the inner bucket plate since there was no experimental data available for the wear of the bucket teeth.

(14)

2 Theory

In this section, the theory used for the thesis will be presented. Rigid body dynamics, the discrete element method (DEM), bonded discrete elements (bonded DE) and Archard’s wear law will be described in this section.

2.1 Rigid body dynamics

This section describes the theory of rigid body dynamics used in LS-DYNA. Rigid bodies are bodies that cannot be deformed and are implemented by the *MAT_RIGID keyword in LS-DYNA. The theory presented is based on the LS-DYNA Theory Manual [6]. Some parts of the derivation have been left out since this part should just give the reader a basic understanding of the principle of motion of rigid bodies.

2.1.1 Equation of motion for a rigid body

By using Newton’s second law of motion for a rigid body the following two equations can be formed. The first equation is:

Mρx = f¨ xρ (1)

where Mρis the physical mass, x is the location of the center of mass ( ˙x and ¨x is the velocity

and acceleration of the center of mass, respectively) and fx

ρ is the forces acting on the rigid body.

The second equation is:

Jρω + ω × J˙ ρω = fωρ (2)

where Jρ is the physical inertia tensor, ω is the angular velocity of the rigid body and fωρ is

the torques acting on the rigid body. The mass and inertia tensor is obtained by integration over the volume of the rigid body or directly from the *PART_INERTIA keyword where the user can specify the properties.

2.1.2 Equation of motion for slave nodes of rigid bodies

The rigid bodies are treated by processing the slave nodes of the rigid body. The slave nodes are defined by the finite element mesh of the rigid body or through nodes that have been added to the rigid body by the *CONSTRAINED_EXTRA_NODES keyword. The velocity of the coordinate of the i:th slave node ˙xi can be established as:

˙

xi= ˙x + ω × ri, (3)

and the rotional velocity of the i:th slave node ωi as:

ωi= ω (4)

where ri= xi− x. The two equations can be put into an equation system as:

" ˙ xi ωi # = " I −Ri 0 I # " ˙ x ω # (5) where I is the identity and Rir = ri× rfor an arbitrary vector r. These equations determines

(15)

The motion of each slave node with index i is also given by the equation of motion as:

mix¨i= fxi (6)

Jiω˙i+ ωi× Jiωi= fωi (7)

2.1.3 Generalized equation of motion for rigid bodies

The generalized equation of motion for rigid bodies can be established by combining the equation of motion for the rigid body (equation 1 and 2) with the equation of motion for the slave nodes (equation 6 and 7) while also satisfying equation 3 and 4. This is done by differentiating equation 3 and 4 with respect to time and inserting the result into equation 6 and 7. By doing this the following equation system can be obtained:

" mi −miRi 0 Ji # " ¨ x ˙ ω # = " fx_i − miω × ω × ri fω_i − ω × Jiω # (8) With the principle of virtual work and by using equation 5 to reduce the number of equations (equation of motion for rigid body and slave nodes) the following equation, which is known as the generalized motion of equation for the rigid body, can be obtained:

M ¨z = fx (9)

J ˙ω = fω (10)

where the mass M is given by

M = Mρ+

X

i∈S

mi (11)

and the generalized forces fx _{are given by:}

fx= fx_ρ+X

i∈S

fx_i (12)

where S is the set of slave nodes. The inertia tensor is given by: J = Jρ+ X i∈S Ji+ X i∈S miRTiRi− M RTz−xRz−x (13)

and the generalized torques as:

fω= −ω × J ω + fω_ρ +X

i∈S

(fω_i + ri× fxi) − rz−x× fx. (14)

This also introduces a variable substitution as: z = x + 1 M X i∈S miri (15) which gives ¨ x = ¨z + Rz−xω − ω × ω × r˙ z−x (16)

(16)

where rz−x= z − xand Rz−xr = rz−x× r

The rigid body accelerations can be determined from the generalized equations 9 and 10 as: ¨ x = f x M (17) ˙ ω = J−1fω. (18)

The translational and rotational increment ∆x and ∆θ can then be computed from the ac-celerations by using explicit time integration with timestep ∆t. The center of mass can then be updated from

xn+1= xn+ ∆x (19)

and the coordinates for each slave node can be updated as:

xn+1_i = xn_i + ∆x + (A − I)rn_i (20)

where A is the transformation matrix that can be computed with the Hughes-Winget algorithm. The velocity of the slave nodes is obtained from:

˙ xi =

xn+1_i − xn i

∆t . (21)

2.2 Discrete element method to model granular material

The Discrete Element Method (DEM) is an approach to model granular material by describing the physical properties of the particles in the material at an individual level and was originally formulated by Cundall and Strack [13]. Properties describing the particle are particle size, shape and stiffness properties. A common model to define these properties of the particle is the Discrete Element Spheres (DES) where the particle is given by a center point and a radius. The theory to describe DEM in Section 2.2.1 and 2.2.2 is based on Cundall and Strack [13].

2.2.1 Newton’s second law to model the motion of DEM particles

Newton’s second law can be used to describe the translational and rotational velocities of the particles. The translational velocity vi, where the subscript i denotes the index of the particle,

can be determined from:

mdvi dt =X j Fc_ij+X k Fnc_ik + Fg_i (22)

where m is the particle mass, Fc

ij is the contact force acting on particle i by particle j and F nc ik

is the non-contact forces acting on particle i from particle j. Fg

i is the gravitational force acting

on particle i and can be written as:

Fg_i = mig (23)

where g is the gravitational acceleration constant. When non-cohesive material is modeled the non-contact forces can be neglected. Non-contact forces are included when modeling capillary forces and other cohesive or adhesive forces between particles. The rotational velocity ωiof particle

(17)

Ii dωi dt = X j Mij (24)

where Ii is the inertia of the particle and Mij is the torque acting on particle i by particle j.

2.2.2 Contact modelling between DEM particles The contact forces denoted Fc

ij from Equation 22, can be determined by setting up a linear spring

and dashpot system between the particles. Fc

ij includes all normal and tangential contact forces.

Forces that are included are contributions from stiffness, damping and friction. A linear spring and dashpot system, between two particles, can be seen in Figure 5.

Figure 5: Linear spring and dashpot system between two particles where k denotes linear spring stiffness, c denotes viscous damping coefficient, n denotes normal direction, t denotes tangential direction, r denotes the particle radius and µ denotes friction coefficient

With the linear spring and dashpot definition the contact force attribution from the linear spring fc _{can be written as:}

fc= δk (25)

where δ is the overlap and k is the linear spring stiffness. The force contribution from the damping fd_{, which controls the energy dissipation, can be written as:}

fd= vc (26)

where v is the relative velocity of the two particles and c is the viscous damping coefficient. The force in the normal direction can then be determined as:

fn = δnkn+ vncn (27)

The force in the tangential direction is computed in the same way. The linear spring and dashpot system is one way of defining the contact in DEM. There are also other ways to define the contact, for example, through the Hertz-Mindlin contact model.

(18)

2.2.3 Contact modelling for particle-particle interaction in LS-DYNA

This section shows how friction, spring stiffness k and viscous dampening coefficient c are modeled and computed in LS-DYNA as stated in Karajan et al. [11].

Contact friction in LS-DYNA is modeled with Coulomb’s friction law as:

ff r= µfn (28)

where µ is the friction coefficient and fn_{is the sum of all contact forces in the normal direction.}

The normal forces can be computed from Equation 27 by only accounting to the normal direction illustrated in Figure 5. The rolling friction is modelled by introducing a moment Mras:

Mr= rµrollfn (29)

where µroll is the rolling resistance.

LS-DYNA computes the spring stiffness in the normal direction kn _{used in Equation 25 as:}

kn = κ1r1κ2r2 κ1r1+ κ2r2

Cnorm (30)

where ri is the radius of the particle, Cnorm is the normal penalty contact stiffness and κi is

given by the compression modulus and is computed as:

κ = E

3(1 − ν) (31)

where E is Young’s modulus and ν is Poisson’s ratio. The tangential spring stiffness kt _is

computed as:

kt= knCtan (32)

where Ctan is the tangential penalty contact stiffness. Also the viscous dampening coefficient

cfrom Equation 26, which is computed in the same way in normal cn and tangential ctdirection, can be computed as:

cn = ct= Dηcrit (33)

where D is the damping coefficient which is in the range of 0 ≥ D ≤ 1.0 and ηcrit is the critical

damping.

The parameters for the contact model for a particle to particle interaction are defined with the *CONTROL_DISCRETE_ELEMENT keyword in LS-DYNA. In this keyword, the user can spec-ify damping coefficients, friction coefficients and spring constants used in the DEM. The particle to surface interaction is defined with the *DEFINE_DE_TO_SURFACE_COUPLING keyword. Discrete element spheres (DES) are defined with the *ELEMENT_DISCRETE_SPHERE key-word. This keyword is generated automatically when DES are created with the DiscGen tool in LS-PrePost.

2.2.4 Time integration

Explicit time integration can be used to determine the time-dependent solution of Equation 22 and 24. The choice of time step ∆tc for this integration has been a discussed topic in order to achieve

numerical stability. In Malone and Xu [9] it is mentioned that the time step has been based on the natural frequency of the equivalent mass-spring system as:

(19)

∆tc= 2

r m

k (34)

From this, a common approach has been developed to use the time step as: ∆tc≤ C

r m

k (35)

where C is a constant, m is the particle mass and k is the stiffness. Malone and Xu [9] mentions that various constants C has been used such as 0.2, 0.4, 0.2π, 2 and 2π.

2.2.5 Bonded discrete elements

This section gives a brief review of the implementation of bonded discrete elements (bonded DE) in LS-DYNA as explained by Karajan et al. [10]. In LS-DYNA bonded particles are applied by creating bonds between nearby particles as illustrated by Figure 6.

Figure 6: Bonded particles. Bonds between particles are represented by the blue lines The bonds are subjected to tension, bending, shearing and twisting which is illustrated in Figure 7.

(a) Tension (b) Bending

(c) Shearing (d) Twisting

Figure 7: Bonds between two particles subjected to the possible force and moment transmission modes

In LS-DYNA bonded particles can be implemented with the *DEFINE_DE_BOND keyword [5]. A bond between two particles is created if the distance between the particles is less than the

(20)

specified distance MAXGAP. In this keyword, the user can specify the maximum normal stress (PBN_S) and maximum shear stress (PBS_S) of the bonds between particles. The bond breaks if the maximum stress is reached. However, the bond becomes unbreakable by setting the maximum stress to 0. In LS-DYNA, the normal bond force between two bonded particles is computed as:

∆fn=

P BN r1+ r2

A∆un (36)

where P BN is the bond modulus (in unit Pascal), A is the effective area of the bond, r1 and

r2 is the radius of the two particles and ∆un is the displacement in the normal direction. The

effective area of the bond A can be computed as:

A = πr2_{ef f} (37)

where the effective radius ref f is the minimum radius of the two particles times a scale factor

SF Aand can be written as:

ref f = min(r1, r2)SF A. (38)

The shear force is computed as:

∆fs= P BS

P BN r1+ r2

A∆us (39)

where P BS is the stiffness ratio (shear stiffness/normal stiffness) and ∆us is the displacement

from shearing.

2.3 Archard’s wear law

Archard’s wear law is a common approach to numerically model the wear. This is done with two different keywords in LS-DYNA. The first keyword *DEFINE_DE_TO_SURFACE_COUPLING models the wear of surfaces in contact with DE rocks. Archard’s wear depth h for DE is computed in LS-DYNA as:

˙h =W EARCfnvt

A (40)

where fn is the normal contact force exerted by the DE, vtis the tangential sliding velocity of

the DE relative to the surface element and A is the area of the surface element.

The second approach is used for surfaces in contact with FE rocks and is applied with the *CON-TACT_ADD_WEAR keyword where the user specifies the contact ID (CID) for the contact surfaces that should include wear. Archard’s wear depth w is then described by:

˙ w = kp ˙d

H (41)

where k is a dimensionless scale factor, p is the contact interface pressure, d is the relative sliding velocity of the contact points and H is the hardness of the surface [5].

(21)

3 Method

This section will describe the method used for the master thesis and will go through the steps of creating the simulation model.

3.1 3D scanning and reverse engineering

A 3D scanning of the Komatsu PC7000 was completed and a point cloud was obtained which can be seen in Figure 8. The main goal of the 3D scan was to capture all dimensions of the moving objects involved in the digging process such as bucket, stick, cylinders and boom. 3D scanning was made with the 3D laser scanner Leica RTC 360.

Figure 8: Point cloud of Komatsu PC7000

A file with stl surfaces for the excavator was obtained from the point cloud and can be seen in Figure 9a. From the stl surfaces a reverse-engineered CAD model was created which can be seen in Figure 9b. The CAD software Siemens NX was used for reverse-engineering.

(22)

(a) STL surfaces from point cloud (b) Reverse engineered CAD model

Figure 9: The retrieved stl-surfaces from the point cloud and reverse-engineered CAD model of Komatsu PC7000

Some areas of the excavator were not captured using 3D scanning. Assumptions were then made to reverse engineer the geometry. A sketch from the service manual was used as a reference when reverse-engineering the bucket dump mechanism. The reversed engineered piston and cylinder pair for the bucket dump can be seen in Figure 10.

(23)

3.2 Multi rigid body simulation

This section describes the rigid body simulation. The section describes the rigid bodies included in the simulation, how the joints are created between rigid bodies and how the motion of the excavator is defined.

3.2.1 Rigid bodies

All parts of the assembly in Figure 9b were meshed individually in NX Pre/Post and were then exported to LS-DYNA format. These parts are referred to as the rigid bodies and a list of the parts included in the simulation can be seen in Table 1. Every rigid body has its own part ID (PID) in the simulation.

Table 1: Rigid bodies

PID Part name

1 Crawler 2 Upper_body 3 Boom 4 Piston_Boom_L 5 Piston_Boom_R 6 Cylinder_Boom_L 7 Cylinder_Boom_R 8 Stick 9 Cylinder_Stick_L 10 Cylinder_Stick_R 11 Piston_Stick_L 12 Piston_Stick_R 13 Rear_Bucket 14 Cylinder_Rear_Bucket_L 15 Cylinder_Rear_Bucket_R 16 Piston_Rear_Bucket_L 17 Piston_Rear_Bucket_R 18 Cylinder_Bucket_Dump_L 19 Cylinder_Bucket_Dump_R 20 Front_Bucket 21 Piston_Bucket_Dump_L 22 Piston_Bucket_Dump_R

To control the mass of the rigid bodies the element type was set to shell for the parts where the mass would influence the results. For the remaining rigid bodies where the mass is not of interest, the element type was set to solid. For example, crawler and upper body. The cylinders and pistons had a reasonable weight using the solid element type and therefore was kept as a solid element type. The weight of the rigid bodies with a shell mesh can be seen in Table 2. To obtain the correct weight of the rigid body the thickness of the shell was adjusted until the weight matched the weight stated in the service information document for Komatsu PC7000.

(24)

Table 2: Weight of rigid bodies with shell mesh PID Part name Weight (kg)

3 Boom 40 647

8 Stick 22 791

13 Rear_Bucket 23 568 20 Front_Bucket 23 163

The meshed rigid body assembly can be seen in Figure 11. Coarse mesh size was given to rigid bodies that will not be in contact with the granular material. Thus a finer mesh was given to the front bucket and the rear bucket. The mesh size is 150 mm on the side and bottom surfaces of the bucket. A mesh size of 75 mm was set to the inside surface of the bucket as well as to the bucket teeth.

Figure 11: Meshed rigid bodies of Komatsu PC7000 in LS-DYNA. The geometry of the upper body is simplified and for example, does not include the cab of the excavator

3.2.2 Defining joints and motors in LS-DYNA

In LS-DYNA a revolute joint is defined with the *CONSTRAINED_JOINT_REVOLUTE key-word by inputting four nodes along the line of revolution. An example is illustrated in Fig-ure 12. Two node pairs are created, 1,2 and 3,4 on the red dashed line of revolution. Node 1, 3 is assigned to rigid body A and node 2, 4 is assigned to rigid body B using the *CON-STRAINED_EXTRA_NODES keyword. The node pairs must be coincident. To create a rota-tional motor a third node pair must be created and assigned to the rigid bodies. This is illustrated as the node 5 and 6 in the figure on the green dashed line. A motor can then be defined with the *CONSTRAINED_JOINT_ROTATIONAL_MOTOR keyword.

(25)

Figure 12: Two rigid bodies A and B which can revolute around the red dashed line A cylindrical joint can be defined using the *CONSTRAINED_JOINT_CYLINDRICAL key-word in LS-DYNA. This is done similarly as a revolute joint however four nodes are created along the centerline instead. An example is illustrated in Figure 13. The four nodes consist of two coincident node pairs 1,2 and 3,4 that must be created as well as assigned to the rigid bodies. Node 1, 3 is assigned to rigid body A and node 2, 4 is assigned to rigid body B. A third node pair is created and assigned to the rigid bodies to create a translational motor with the *CON-STRAINED_JOINT_TRANSLATIONAL_MOTOR keyword.

(26)

Figure 13: The cylindrical joint between rigid body A and rigid body B allows rotation around the red dashed center line as well as translation along the centerline

3.2.3 Joints and motors in the rigid body simulation

Joints and motors were defined to simulate the motion of the hydraulic excavator. These consisted of revolute and cylindrical joints together with rotational and translational motors. The positioning of the joints can be seen in Figure 14.

Figure 14: Joints of Komatsu PC7000. The red arrows indicate a revolute joint and the blue arrows indicate a cylindrical joint

(27)

joint number. The table includes the part name of the two bodies (A and B) for which the joint acts between.

Table 3: Joints and motors in multi rigid body simulation

Nr Joint Motor Part A Part B

1 Revolute Rotational Crawler Upper_body

2 Revolute Upper_body Boom

3 Revolute Upper_body Piston_Boom

4 Revolute Boom Cylinder_Boom

5 Cylindrical Translational Piston_Boom Cylinder_Boom

6 Revolute Stick Boom

7 Revolute Cylinder_Stick Boom

8 Revolute Piston_Stick Stick

9 Cylindrical Translational Piston_Stick Cylinder_Stick

10 Revolute Cylinder_Rear_Bucket Rear_Bucket

11 Revolute Stick Rear_Bucket

12 Revolute Piston_Rear_Bucket Boom

13 Cylindrical Translational Piston_Rear_Bucket Cylinder_Rear_Bucket

14 Revolute Front_Bucket Rear_Bucket

15 Revolute Cylinder_Bucket Rear_Bucket

16 Revolute Piston_Bucket Front_Bucket

17 Cylindrical Translational Piston_Bucket Cylinder_Bucket 3.2.4 Defining motion of joints

The motion of the bucket is controlled through the displacement of four translational joints and one rotational joint. The displacement is given by a curve in LS-DYNA with the *DEFINE_CURVE keyword. These curves can be seen in Figure 15. To determine how the excavator operates several videos have been studied and have been used as inspiration to create these curves [14]. The operation starts by obtaining the correct initial position of the cylinders. This part includes withdrawing the rear bucket cylinder between 0 to 5 s. After 5 s the digging operation begins by extending the stick cylinder between 5 to 15 s. The rear bucket cylinder starts extending between 10 and 15 s. When the stick and rear bucket cylinder stops the boom cylinder extends between 15 to 20 s to remove the bucket off the rock pile and to obtain a position of the bucket above the dumper. Between 20 s and 25 s, the excavator rotates until the bucket is positioned above the back of the dumper. The bucket dump cylinder withdraws between 30 to 35 s to dump the rocks onto the dumper. The bucket dump cylinder extends between 35 to 40 s to close the bucket. The loading process is completed and the excavator rotates back to its initial position between 40 and 45 s. The rest of the cylinders also moves to their initial position

(28)

Figure 15: Displacement curves for the motion of the four translational joints and one rotational joint

3.3 Granular simulation

In this section, the two different methods to simulate granular material will be described. The granular material is simulated as rigid with the following material parameters: density 3700 kg/m3_,

Young’s modulus 0.45 GPa and Poisson’s ratio 0.30. In Table 4 the parameters that control the behavior of the DES are listed. The parameters have remained the same as in Svanberg et al. [1] in which they performed calibration of these parameters. The parameters can be specified in the *CONTROL_DISCRETE_ELEMENT keyword.

Table 4: Parameters for discrete elements

Parameter Description Value

NDAMP Normal damping coefficient 0.9 TDAMP Tangential damping coefficient 0.75

FRICS Friction coefficient 0.9

FRICR Rolling friction coefficient 0.05

NORMK Normal spring constant 0.05

SHEARK Shear spring constant 0.27

3.3.1 FE-DE method to model granular material

In Svanberg et al. [1] a variety of rocks were 3D scanned to obtain the general geometry of rocks in a rock pile at Aitik. The FE representation of the larger rocks can be seen in Figure 16a. The number of each rock was created such that it follows the particle size distribution at Aitik. The empty spaces between the rocks were filled with smaller rocks represented by DES with a radius of 100 mm. This is seen in Figure 16b.

(29)

(a) FE representation of larger rocks (b) DES representation of smaller rocks

Figure 16: FE and DES representation for common shapes of rock fragmentation found at Aitik The calibrated rock pile for the FE-DE method can be seen in Figure 17. The rock pile was calibrated to an angle of repose of 35 degrees which is the general angle of repose for rock piles in Aitik. The same calibration method as for the bonded DE method, described in Figure 23, was used to initialize the rock pile for the FE-DE method.

Figure 17: Rock pile with FE-DE representation. Some rocks have been removed under the pile as well as on the top to reduce computational time. The bottom layer of DES has been fixed to prevent the pile from collapsing. This pile includes 76 416 DE elements and 198 237 FE elements 3.3.2 Contact definition for FE-DE method

For the FE-DE method, there are five different contact definitions to be specified. Between rocks to bucket, rocks to rocks, rocks to walls, rocks to ground and rocks to dumper. In this section, FS and FD are referred to as the static coefficient of friction and dynamic coefficient of friction, respectively. SFS and SFM are the slave penalty stiffness of the slave and master surface, respectively.

The contact definition between FE dumper to DE rocks, FE walls to DE rocks and FE rocks to DE rocks was established with the *CONTACT_AUTOMATIC_NODES_TO_SURFACE key-word. The parameters for these contact definitions are listed in Table 5.

(30)

Table 5: Contact parameters with nodes to surface keyword

Parts FS/FD SFS/SFM

FE dumper to DE rocks 0.85 0.5

FE walls to DE rocks 0 0.5

FE rocks to DE rocks 0.85 0.1

The *CONTACT_AUTOMATIC_SINGLE_SURFACE keyword was used to establish contact between FE rocks themselves. It was used with a pinball segment-based contact (SOFT = 2) and friction coefficients were set to 0.85. Contact definitions between FE rocks to FE dumper, FE rocks to FE ground, FE rocks to FE walls and FE rocks to FE bucket were established with the *CONTACT_AUTOMATIC_SURFACE_TO_SURFACE keyword. The parameters for these contact definitions are listed in Table 6.

Table 6: Contact parameters with surface to surface keyword

Parts FS/FD SFS/SFM SOFT

FE rocks to FE dumper 0.9 1.0 2

FE rocks to FE ground 0.9 1.0 2

FE rocks to FE walls 0.9 1.0 2

FE rocks to FE bucket 0.53 0.53 0

For the FE rocks to FE bucket contact a penalty formulation was used (SOFT = 0) to allow contact wear between the parts. Wear was added to the model by using the *CON-TACT_ADD_WEAR keyword. This interface is shown in Figure 18.

Figure 18: Keyword for wear between FEM rocks and bucket. WTYPE = 0 sets the wear law to Archard’s, P1 is the dimensionless wear factor k, P2 is the hardness of the slave surface and P3 is the hardness of the master surface

The wear coefficients were obtained from a calibrated model which uses a user-defined wear routine that takes both impact wear and sliding wear into account. This model was calibrated to an in-situ experiment with a Hardox 450 wear plate. This gave the wear coefficient W EARC =

k

H = 3.96·10

−10_{. The dimensionless scale factor k is then applied as k = W EARC ·H = 1.98·10}−7_.

The hardness H was set to 500. If k is inserted into Equation 41 then the hardness H is canceled out and thus the equations are only dependent on the calibrated WEARC coefficient and it does not matter what H is set to.

The *DEFINE_DE_TO_SURFACE_COUPLING keyword was used to establish contact be-tween DE rocks to FE bucket and DE rocks to ground. The parameters for this keyword are listed in Table 7. To add wear to the simulation the wear coefficient WEARC was set to 3.96e-10 in the keyword for the bucket contact.

(31)

Table 7: Contact parameters with surface coupling keyword

Parts FricS FricD

DE rocks to FE ground 0.9 0.3 DE rocks to FE bucket 0.53 0.05 3.3.3 Bonded DE method to model granular material

The FE rocks from Figure 16a were used to create bonded DE that matched the shape of the rocks. The radius was kept to a minimum of 100 mm to be able to keep the same time step for the simulation as in the FE-DE method. DES was fitted to the rocks with the DiscGen tool in LS-PrePost and the bonded DES can be seen in Figure 19.

Figure 19: Bonded DES fitted to FE rocks. The FE rocks are replaced by one very large rock with eight spheres in a quad pattern, one large rock with six spheres in a triangular pattern, one large rock with three spheres in a straight pattern, three smaller rocks with six spheres each in a triangular pattern and 13 smaller rocks with three spheres each in a straight pattern

The bonded DE can be seen in Figure 20a and the empty spaces between the bonded DES were filled with 134 normal DES with a radius of 100 mm. This is seen in Figure 20b.

(a) Bonded DES representation of larger rocks (b) DES representation of smaller rocks together with the bonded DES.

Figure 20: Bonded DES and DES representation for common shapes of rock fragmentation found at Aitik

(32)

Two different node sets were created for the DE elements. One node-set containing all the nodes of the bonded DE and one node set for all normal DE together with the bonded DE. Bonded DE was created with the *DEFINE_DE_BOND keyword in LS-DYNA with the inputs as shown in Figure 21. PBN and PBS specifies the parallel-bond modulus [Pa] and parallel-bond stiffness ratio (shear stiffness/normal stiffness), respectively. PBN S and PBS S specifies the maximum normal stress and maximum shear stress, respectively. ALPHA defines numerical damping and MAXGAP defines the maximum gap between two bonded spheres. If the distance between two spheres is less than MAXGAP, a bond between the spheres is created.

Figure 21: Inputs for bonded DE in LS-DYNA

PBN S and PBS S were set to 0 which indicates infinite maximum stress and thus the bonds will be unbreakable. PBN and PBS were modified until a reasonable behavior between two bonded spheres was achieved i.e. such that the particles do not oscillate relative to each other. PBN was set to 1000 and PBS to 1. Numerical damping was introduced to the system by setting ALPHA to 0.9. Numerical damping was introduced to achieve a stable behavior of the particles. MAXGAP was set to -20 (mm). A value of MAXGAP less than 0 indicates that an absolute value will be used. The node-set containing all DES to be bonded (displayed in Figure 20a) was chosen as SID. The rock pile for the bonded DE method can be seen in Figure 22.

(33)

Figure 22: Rock pile with bonded DE representation. Bonded DE represented by blue color and single DES by orange color. Rocks have been removed on the top and in the bottom of the pile to reduce computational time. The bottom layer of DES is fixed to prevent the pile from collapsing. This pile includes a total of 142 292 DE elements out of which 51 534 are bonded DE elements 3.3.4 Contact definition for bonded DE method

Contacts were defined by using the *DEFINE_DE_TO_SURFACE_COUPLING keyword for the bonded DE method. The parameters for the contacts are listed in Table 8. FricS and FricD are the friction coefficient and rolling friction coefficient, respectively. SFP is the scale factor on contact stiffness.

Table 8: Contact definitions for bonded DE method

Parts FricS FricD SFP

DE rocks to FE dumper 0.9 0.3 1.0 DE rocks to FE ground 0.9 0.3 1.0 DE rocks to FE walls 0.9 0.3 1.0 DE rocks to FE bucket 0.53 0.05 0.01

The SFP was set to 0.01 for the DE rocks in contact with the FE bucket. This was done to achieve a stable behavior for the contact. When the SFP was set to 1.0 the bucket could "shoot" away the particle because the bucket teeth came in between two bonded particles. By lowering the SFP to 0.01 this problem was avoided.

3.3.5 Calibration and initialization of bonded DE rock pile

The rock piles have been initialized with the same approach in both the bonded DE method as well as the FE-DE method. This approach is described by Figure 23. At t = 0s the rocks are

(34)

spread out evenly in different segments. At this stage, the rocks in the bonded DE method have the distribution as shown in Figure 20b. For the FE-DE method, the rocks have the distribution as shown in Figure 16b. The rocks fall because of gravity landing on the rigid wall which has an angle relative to the horizontal ground. This rigid wall begins to move up vertically at t = 2s. It moves 11 m under a period of 18s. At t = 20s the rigid wall becomes stationary. The rocks continue to fall to their equilibrium position until t = 45s when the simulation ends. The aim was to obtain an angle of repose for the pile at 35 degrees which is a good estimate for the piles after rock blasting in Aitik.

(a) t=0s (b) t=2s

(c) t=20s (d) t=45s

Figure 23: Calibration method of rock piles

(35)

Figure 23d an initialization of the rock pile must be made, for the bonded DE method. This is because it is not possible to start the simulation at t=45s since all bonded particles would bond to other bonded particles and form a big unmovable rock. This is because LS-DYNA has not implemented any possibility to save information of bonded discrete elements.

Therefore all new simulations are made from t = 0s in Figure 23a and initialized to t=45s by using the keyword *BOUNDARY_PRESCRIBED_FINAL_GEOMETRY. At time t=0s there is enough space between bonded particles to prevent them from bonding to other bonded particles. The inputs for this keyword are shown in Figure 24.

Figure 24: Inputs to the *BOUNDARY_PRESCRIBED_FINAL_GEOMETRY keyword The (x, y, z) coordinates of each node in the rock pile (from Figure 23d) are inserted into the keyword. LCID and DEATH for each node are kept as 0 so it will automatically choose LCIDF and DEATHD as input. LCIDF should be defined as a curve between 0 to 1. At 0 the nodes will have their initial coordinates and at 1 the nodes will reach their specified position. The nodes will then move in a straight-line trajectory to their specified coordinates from their initial coordinates. At time DEATHD all nodes will be able to move freely and are not controlled by the keyword. *CONTROL_DISCRETE_ELEMENT is turned off and is given a birth time after the initialization has been made. Gravity was also turned off during the initialization. This was done to prevent contact forces and bond forces between the rocks. The initialization process is seen in Figure 25.

(36)

(a) Initial position. LCIDF curve=0, t=0s

(b) Nodes moving in a straight line trajec-tory. 0<LCIDF-curve<1, t=0.3s

(c) Nodes reached speci-fied coordinates. LCIDF curve=1, t=0.5s

Figure 25: Initialization of rock pile with *BOUNDARY_PRESCRIBED_FINAL_GEOMETRY with DEATHD=1s. Blue spheres are bonded DE and orange spheres are single spheres

3.3.6 Calibration of bulk density

The method to determine the bulk density of the granular material is described in Figure 26. Granular material is filled into a box with a known volume. The total mass of the rocks occupied in the box is then computed and divided by the known volume of the box to obtain the bulk density. In this case, a 5x5x5 m box was used. The density of the rocks is varied until a bulk density of 1900 kg/m3 _{has been obtained. For the bonded DE method, the rock density was determined as}

3100 kg/m3_{. The density is lower compared to the density used for the FE-DE method which is}

(37)

(a) t=0s. Granular material spaced out evenly in the same way as in the calibration of the rock pile

(b) t=10s. Both the funnel and box have been filled with granular material due to gravity act-ing on the rocks

(c) t=12s. The funnel is being moved away from the box, de-activating any rocks outside of the box

(d) t=20s. Box has been filled and all rocks outside of the box has been removed from the simulation

(38)

4 Results

The results of the master thesis are presented in this section.

4.1 Multi rigid body simulation with bonded DE granular material

model

The complete simulation process of the loading with the bonded DE granular material model can be seen in Figure 27. The figure shows snapshots of a complete cycle beginning at t = 0s and ends at t = 40s.

(a) t=0s (b) t=5s

(c) t=10s (d) t=15s

(e) t=19s (f) t=25s

(g) t=30s (h) t=40s

Figure 27: Simulation of Komatsu PC7000 with multi rigid body dynamics and bonded DE granular material model. Bonded particles in blue and single DES in orange. Cycle time 40 seconds

(39)

4.2 Multi rigid body simulation with FE-DE granular material model

The complete simulation process of the loading with the FE-DE granular material model can be seen in Figure 28.

(a) t=0s (b) t=5s

(c) t=10s (d) t=15s

(e) t=19s (f) t=25s

(g) t=30s (h) t=40s

Figure 28: Simulation of Komatsu PC7000 with multi rigid body dynamics and FE-DE granular material model

(40)

4.3 Validation of simulation results

In Figure 29 the Komatsu PC7000 can be seen operating in the Aitik mine. The time between the first snapshot to the last snapshot was 36 seconds which is the cycle time.

(a) t=0s (b) t=6s

(c) t=12s (d) t=17s

(e) t=20s (f) t=23s

(g) t=26s (h) t=36s

(41)

4.4 Comparison between granular material models

Figure 30 shows the bucket filling process from a side view for the bonded DE method and FE-DE method.

(a) t=0s (b) t=5s (c) t=10s (d) t=15s

(e) t=0s _{(f) t=5s} _{(g) t=10s} _{(h) t=15s}

Figure 30: Bucket filling using bonded DE method (a)-(d) and FE-DE method (e)-(h) illustrated with transparent bucket

(42)

(a) t=0s _{(b) t=0s}

(c) t=5s _{(d) t=5s.}

(e) t=10s (f) t=10s

(g) t=15s (h) t=15s

Figure 31: Bucket filling with bonded DE method to the left and FE-DE method to the right The two granular material models were compared through a video of the loading with a section view similar to Figure 30. With this video, it was seen that granular material tends to "roll" into the bucket for the FE-DE method. With the bonded DE method it was seen that the granular material slides into the bucket and does not roll as much as seen in the FE-DE method.

4.5 Simulation of wear

For this part, the rotational motion of the excavator was removed which means that the excavator never rotates to the dumper to unload the rocks. Instead, the rocks are unloaded back into the rock pile. This prevents empty spaces in the rock pile and thus it is easier to fill the bucket repetitively. Five bucket fillings were carried out during the first simulation. The amount of rocks per filling for the respective method is listed in Table 9.

(43)

Table 9: Mass of rocks for each bucket filling in unit tonnes

Filling Nr Method

Bonded DE FE-DE Experimental

1 45,3 46,8 55 2 48,5 48,1 59 3 49,8 46,8 61 4 48,4 50,6 57 5 50,1 49,2 52 Average 48,4 48,3 56.8 Total 242,1 241,5 284

The bonded DE simulation time for five repeated fillings was 29 hours and 27 minutes with parallel execution on 32 CPUs. For the FE-DE method, the simulation time was 99 hours and 4 minutes with the same execution. In Figure 32 the five filled buckets for both methods can be seen.

(a) (b) (c) (d) (e)

(f) (g) (h) (i) (j)

Figure 32: (a)-(e) shows bucket filling 1-5 with the bonded DE method and (f)-(j) shows bucket filling 1-5 with the FE-DE method

The wear pattern after five repeated fillings with the bonded DE method and FE-DE method can be seen in Figure 33.

(44)

(a) Bonded DE method (b) FE-DE method

Figure 33: Absolute wear pattern achieved on bucket after five bucket fillings in unit mm The focus of the wear study was on the inside of the bucket and therefore the wear plate rows inside the bucket were extracted from the previous figure. The wear on these rows can be seen in Figure 34a and 34b after five bucket fillings for the bonded DE and FE-DE method, respectively.

(45)

(a) Bonded DE method

(b) FE-DE method

Figure 34: Absolute wear pattern achieved on inner bucket after five bucket fillings in unit mm A comparison for the wear after 30 loading cycles was also completed. This was done by adding the wear from six different simulations with each simulation consisting of five bucket fillings. A total of 1376 tonnes with an average of 45,9 tonnes per bucket was filled with the FE-DE method. A total of 1405 tonnes with an average of 46,9 tonnes per bucket was filled with the bonded DE method. Simulations were made when the bucket was positioned higher up in the rock pile or lower down in the rock pile. For each of the two locations, a simulation was made in the left, middle and right part of the rock pile. The mass for each bucket filling per location can be seen in Appendix A and B for the FE-DE and bonded DE method, respectively. From the 6 simulations, the average simulation time for the bonded DE method and FE-DE method was 31h, 12min and 102h respectively.

The resulting wear map from the bonded DE and FE-DE method can be seen in Figure 35a and 35b, respectively.

(46)

(b) FE-DE method

Figure 35: Absolute wear pattern achieved on inner bucket after 30 bucket fillings in unit mm The wear result from Figure 35 was normalized and any outliers were removed with the rmoutlier function in MATLAB. The rmoutlier function removes any outlier from the data if it is more than three standard deviations from the mean. This was done because in the wear results there were some nodes with unreasonable wear compared to the neighboring nodes. The unreasonable wear in some nodes seems to register due to irregularities in the model, for example, the upper left red spot in Figure 35a was caused by a single DES penetrating the bucket in one of the 30 bucket fillings. The reason for this might be that the contact between particle and surface is not stiff enough. If the data were normalized without removing the outliers the normalized data would be biased to the outliers. From this criteria, a total of 21 out of 1758 nodes were neglected from the bonded DE method and 32 out of 1758 nodes from the FE-DE method. Normalization was done to better compare the results independent of the range. Wear results with normalized data and outliers removed can be seen in Figure 36.

(47)

(b) FE-DE method

Figure 36: Normalized wear pattern after 30 bucket fillings with range 0 to 1

By setting the lower bound on the color scale to 0.3 it is more clear to see where the wear is most extensive. This is shown in Figure 37

(48)

(b) FE-DE method

Figure 37: Normalized wear pattern after 30 bucket fillings with range 0.3 to 1

From the FE-DE method, the average nodal wear contribution from the DE rocks and FE rocks was 1,26E-06 mm and 5,63E-06 mm, respectively per five bucket fillings. This information can also be found in Appendix A.

From the bonded DE method, the average nodal wear from the DE rocks was 17,8E-06 mm or about 2.5 times more than the total average nodal wear in the FE-DE method. This information can also be found in Appendix B.

4.6 Validation of wear

The Komatsu PC7000 was followed for 86,67 days when operating in Aitik. At the end of the period, the bucket was examined with an ultrasonic thickness gauge (KARL DEUTSCH ECHOMETER 1075 Basic) to obtain the thickness of some points inside the bucket. For this period the total

(49)

tonnes of rocks loaded was 1289,4 ktonnes. The bucket had an initial plate thickness of 25mm. When subtracting the measured thickness from the initial plate thickness the wear depth as shown in Figure 38 is obtained.

Figure 38: Wear depth at some points in the bucket in unit mm

A surface fitting z(x,y), where z is the wear depth, was applied to the normalized measuring points. Normalization was done because the wear ranges from 0-11.5 mm in the experimental data and 0-1.0 µm in the simulation. With normalization, the map becomes independent of the range and the range after normalization is 0-1. From this, the wear map of the experimental results in Figure 39 was obtained. The x-axis is given by the width position inside the bucket and the y-axis is given by the depth position inside the bucket. Coordinate (x,y) = (0,0) is the inner left corner of the bucket when looking from the front of the bucket.

(50)

Figure 39: Normalized wear map obtained from experimental measurements. The black points indicate the position of the measurement points

The measuring points from Figure 38 were used to align two 3D scanned point clouds with each other. These 3D scans included a scan of the original bucket plate with an initial 25 mm plate and the second scan was the bucket plate after loading 1300 ktonnes. After the point clouds are aligned with each other, the wear map in Figure 40 can be established. The results from the middle of the bucket are uncertain where the color shows blue/green because there was water located inside the bucket during the first scan. From the first scan, no points were obtained of the bucket where the water was located. Instead, a CAD stand-in was inserted, where the water was located, together with the remaining point cloud. For the second scan, there was no water inside the bucket. This is also the reason why there is a sudden shift in the color scale in the figure. The color goes from yellow/red to green/blue with no transition phase. Therefore, the middle area with the CAD stand-in should not be compared to the back and the front of the bucket and should only be compared to the area where the CAD stand-in is included.

Figure 40: Experimental wear depth inside the bucket obtained from a comparison of two 3D scanned point clouds in unit mm

A bonded discrete element approach to simulate loading with hydraulic mining excavators

Preface

Abstract

Sammanfattning

Contents

1

Introduction

1.1

Background

1.2

Aim and objective

1.3

Problem description

1.4

Limitations

2

Theory

2.1

Rigid body dynamics

2.2

Discrete element method to model granular material

2.3

Archard’s wear law

3

Method

3.1

3D scanning and reverse engineering

3.2

Multi rigid body simulation

3.3

Granular simulation

4

Results

4.1

Multi rigid body simulation with bonded DE granular material

model

4.2

Multi rigid body simulation with FE-DE granular material model

4.3

Validation of simulation results

4.4

Comparison between granular material models

4.5

Simulation of wear

4.6

Validation of wear