Full-scale simulation and validation of bucket filling for a mining rope shovel by using a combined rigid FE-DEM granular material model

Rope shovels and other heavy mining equipment used for loading fragmented rocks to extract minerals from the earth are used in almost every open pit mine. The optimization of the loading process is of enormous value due to the extremely large amount of material turn over. In this work, a full-scale numerical model of the loading process is developed. Granular material of copper ore is modeled in a combination of rigid finite elements for larger particles with complex shapes, and the discrete element method for smaller particles. A multi rigid body dynamic model, discretized with finite elements are used to model the rope shovel. Calibration of the numerical model for the granular material is performed via a new and unique experimental full-scale approach of analyzing waste rock pile angles with a height of approximately 15 m. In situ experimental data acquisition is performed during the loading process for validation of the model. After model validation, the influence of several loading variables such as bucket rake angle, velocity, and position from the pile are investigated and evaluated. When comparing the numerical model results with experimental mass measurement an excellent agreement was observed. Also, drone camera video recordings of the mass flow behavior and the numerical mass flow behavior are in agreement. Small adjustments of dig variables show a significant effect on the average dig force as well as the bucket fill factor.


Introduction
The loading process of ore and waste rock is a very critical process in the production flow in any open-pit mine. Usually, rope shovels or other mining excavators such as hydraulic mining shovels are used for excavating the material onto haul trucks. The material is then transported to a crusher or waste rock dump site depending on if the material is ore with profitable mineral content or waste rock. A typical mining rope shovel has a bucket volume capacity of around 30-60 m 3 and fills a dump truck with a payload of over 300 tons with 3-4 buckets. Normally, a rope shovel fills around 1000-3000 buckets corresponding to 30-70 kton of ore and waste rock daily. Naturally, the bucket of a rope shovel is the part of the machine that is subjected to the immediate B Andreas Svanberg Andreas.Svanberg2@Boliden.com; andreas.svanberg@ltu.se 1 Boliden Minerals AB, Aitik Mine, 98292 Sakajärvi, Sweden 2 Division of Mechanics of Solid Materials, Luleå University of Technology, 97187 Luleå, Sweden 3 Boliden Minerals AB, 93632 Boliden, Sweden interaction with the granular material. Several factors are of significant importance to ensure that the rope shovel operates efficiently and that the machine downtime is kept to an absolute minimum. It has been found that several factors have a great influence on how efficient a bucket is filled. Therefore it is of high interest to obtain a better understanding of how the material flows into the bucket in order to make operational decisions. Improved knowledge and understanding are important in order to optimize and evaluate the process to perform at maximum efficiency, e.g., to make decisions on what equipment is best, how to operate the rope shovel for best bucket fill factor, and identify critical damage and wear areas. A slight increase in bucket fill factor can save thousands of bucket loads annually, and thus lead to increased productivity and decreased maintenance costs.
Numerical modeling and simulation of granular material have in the last decades been extensively studied with different methods and for different applications, with the aim to capture the very complex nature of granular material and its interaction with structures. Numerical models and simulations can provide insight into processes that are difficult or impossible to study experimentally. Hence, in order to further improve the knowledge and the understanding of granular processes, simulations are a powerful tool. It should be stressed that validation is of utmost importance in order to use a numerical tool for industrial decision making.
Perhaps the most common approach to model granular materials is by the use of the discrete element method (DEM). The DEM was originally developed by Cundall and Strack [13,14] for the study of granular assemblies. Since then, the DEM has been used directly or in combination with other numerical methods such as the particle finite element method (PFEM) or smoothed particle hydrodynamics (SPH) for many applications such as safety berms for haul trucks in mines [44], material dumping and loading process [16,17,20], blast furnaces and granular column collapse [27,34], and grinding in tumbling and stirred media mills [4,[22][23][24]28,29,48] Traditionally, DEM describes each particle in a granular assembly with a spherical geometry. The fact that each granular particle is described by a center point and a radius simplifies the numerical computations significantly. A perfect sphere has a very different behavior when in movement and in contact with another particle compared to, e.g., a rock with large aspect ratio and more angularity, and is the reason for research on particle shape affects [7,26]. The most common solution to this in DEM-simulations is to add artificial friction, usually called rolling friction, which can be applied to the spheres in order to change their behavior at contact. The influence and accuracy of utilizing rolling friction have been studied and it has been shown that the rolling friction itself can not capture the non-spherical behavior alone [46]. The DEM has the possibility of modeling other geometries than spheres. One way of modeling the shapes is through the use of continuous super quadratic functions [8]. Other geometries than spheres lead to greater computational demand. The continuous functions have a limited selection of geometries since the geometry needs to be described by a function. Also, sharp corners are almost impossible to obtain through continuous functions. Another common method to obtain more complex particle shapes is by using clusters of particles. This method uses several discrete element spheres to obtain one particle, by such the simple sphere geometry can still be used, and thus be computationally efficient, but still generating a particle shape that is not a sphere. The method with particle clusters has been used for different applications, [42]. Once again, due to the spherical shape, e.g., sharp corners and small features require very small radii. The explicit time step in DEM is directly dependent on the particle radius and it is most often not feasible with too complex geometries due to this time step dependency. Also, with smaller particles, the particle count will drastically increase, and thus increase simulation time. Some researchers have used DEM particles with polyhedral shapes [36,43]. In the polyhedral DEM the particle shape can arbitrarily be chosen due to that the parti-cle geometry is defined by a number of flat polygonal faces. For this reason, arbitrary shapes can be modeled, and the resolution of the particle shape can be defined by adjusting the number of faces or the mesh representing the geometry. The polyhedral DEM has quite recently (compared to the spherical DEM) received more attention. Some commercial software utilizes the polyhedral DEM and perhaps the most common in the mining industry is the Rocky DEM. However, the majority of available DEM softwares does not have any option for polyhedral particle shapes, though implementation in open-source codes may be performed [35]. Latham et. al. [30] modeled particles with combined FEM and DEM (FDEM) for a range of applications and also gave a review on currently available methods for particle shape representation. Zhong et. al. [51], did a thorough review of recent developments and methods to model a non-spherical particulate system. Gustafsson et al. [19] showed that the multi-particle finite element method can be used to simulate granular material. This method is especially powerful for complex-shaped particles and when the particle itself is studied in detail.
In both experimental and numerical research aimed for industry, it has been very common to perform the research on smaller scales than the actual process, e.g., in pilot or lab scale. There is a clear benefit of doing laboratory-scale research, due to the controlled environment, and due to laboratory-scale research, some of the most important scientific discoveries have been revealed. The disadvantage of the pilot or laboratory-scale studies is that the industrial process such as the loading process is at a very large scale, and it is not adequate to have several 1-m large rocks inside any traditional laboratory. Since the translation from laboratory to full scale is not always straight forward, and sometimes not even possible, large scale research is of importance in order to study mechanisms at their real scale. In the last decades, the numerical capacity has increased and algorithms have been developed. Full and large scale simulations with DEM of up to 128 million particles have been performed [8,45].
Numerical and experimental studies in mining-related granular flows have been performed previously. Rasimarzabadi et al. [40] did a thoroughly experimentally study, but also a numerical study [39] where different factors such as rake angle, velocities, and size through the filling of a scaled rope shovel bucket with crushed limestone in a controlled environment. The bucket was transparent, to make it possible to study the flow into the bucket. Curley et al. [15] used a 2D DEM model to study the impact of particle size characteristics of the experimental work from Rasimarzabadi. Another research group studied how the digging forces and power consumption on a scaled rope shovel could be predicted experimentally but also through utilizing DEM with clustered particles [3,47]. Cleary, [5,6,8] created a numerical model to simulate the dragline bucket filling both in 2D and 3D to investigate the flow into buckets and also how it is affected by e.g. boulders (large rocks) in the bulk material or muck-pile. Coetzee, [10][11][12] have been studying the bucket filling with both discrete (DEM) and continuum methods such as material-point method (MPM). Nezami [36] used DEM with polyhedral elements to model the bucket filling of a wheel loader in order to characterize force feedback during soil and bucket interaction.
One objective of this paper is to develop experimental methods that better support calibration and validation of granular models for industrial use. A second objective is to utilize a numerical model for an industrial full-scale study in order to better understand the loading process and in the extent, with the aim to load ore more efficient. In this work, a Bucyrus rope shovel used in the Boliden's Aitik copper mine located in northern Sweden was studied during the loading operation of copper ore and waste rock. Movement data, geometry, and relative positions were recorded experimentally. The granular ore was modeled with a unique combination of DEM and rigid FE for simulating the smaller, and larger particles respectively with realistic geometries. The geometry of particles was 3D laser scanned and a semidynamic full-scale calibration of the model was performed. A multi rigid body dynamic FE-model of the rope shovel, with 3D laser-scanned bucket and the relevant positions of bucket and muck pile were used to load the granular material. Furthermore, the loading process was validated at full scale with video recording and mass measurement. Finally, the model was utilized to investigate and evaluate critical loading parameters.

Materials
In Aitik, the active mining area was several kilometers in all directions, thus the geology varied depending on where in the mine the ore and waste rock were extracted from. The copper, gold and silver content in the ore in Aitik is in the range of 2000-3000, 0.2 and 2 g/ton respectively, where the remaining is defined as waste rock, which thus makes the difference between waste rock and ore on a bulk level hard to distinguish between due to the very small metal content. In the present study, it is assumed that the copper, gold, and silver ore as well as the waste rock behaves similarly on a material flow level, and thus can be treated as one material. The material is later referred to as the granular material or granular ore.
The rock types in Aitik were mainly pegmatite, muscovite schist, banded hornblende, and amphibole-biotite. The mineralogical composition in the rocks varies but the largest contents (approximately 92 % by volume) are quartz, plagioclase, K-feldspar, biotite, and muscovite [2,50]. Original PSD curve adopted from [37] Particle size distribution (PSD) of the material is shown in Fig. 1. PSD in the Aitik mine was previously studied [37,38], and will be applied in the present study as well, although the position in the mine was not identical. The material had a bulk density of approximately 1800-2000 kg/m 3 , and particle densities varying from 2600 to 2970 kg/m 3 between the rock types [50]. The rock shapes were mostly angular and blocky. The Aitik mine utilized a mining method, where benches of 15 m were used. After blasting, the muck-pile height was approximately 18 m, due to the swell factor.
Arbitrary rocks from a muck-pile were 3D-scanned with a Leica RTC 360, seen in Fig. 2. To obtain a granular model containing a variation in both size and shape, seven different rocks were randomly picked from a loading site.

Experimental procedures
A Bucyrus 495 rope shovel was experimentally observed during the loading process. In Fig. 3a, a picture of the studied rope shovel is given with its main components. Four fundamental components move on a rope shovel, all moving with the intention to maneuver the bucket. Swinging allows for rotation between the upper and lower body. Propelling is caused by the rotation of the crawlers. The hoist ropes can be extended and retracted in order to lift and sink the bucket respectively. The crowd arm has the possibility to be translated along the saddle block, and the saddle block can pivot around an axis corresponding to the position of number 8 in Fig. 3. The bucket is tied to the crowd arm, which makes the movement of the bucket identical to the movement of the crowd. Figure 4a presents the bucket assembly with the main components. The width of the bucket is defined as the distance between the two sidewalls, and the length is defined as the  (7), and pivoting point (8) between crowd arm and boom distance from the front lip to the door. Height is defined as the distance between the bottom and top of the bucket. Bucket rake angle is defined as the angle between one line, drawn along the centerline of the dipper handle, and another line, drawn from the very bottom of the bucket heel, toward the tip of the teeth. The rake angle is presented in Fig. 4b. The bucket fill factor is determined as the ratio of achieved volume and the volumetric rating (struck volume) of the bucket.
Initially, the geometry of the rope shovel and the bucket was 3D scanned in high resolution with more than 10 million scan points. The 3D scanner was of model Leica RTC 360 and the bucket was an ESCO master dipper of 42 m 3 (55 cubic yard) volume capacity. The geometry of the muck-pile and the position of the rope shovel relative to each other were also detected through 3D scanning right before the loading started. During the loading process, data was recorded from an HMI (Human machine interface) installed on the rope shovel. Displacement of the crowd, hoist ropes, and rotational velocity of the swing were recorded. The sampling frequency was 0.5-1 Hz for all the outputs. The shovel operator was instructed to dig in a controlled manner and to try to repeat the same dig operation for every bucket filling. This was done to keep the operational variables to a minimum. During the loading, a drone (Mavic 2 Enterprise DUAL) with a mounted camera of 12 megapixels resolution and recording frequency of 30 fps was recording the bucket filling from a bird's eye view. The mass in the bucket was measured when it was dumped on the haul trucks of type, Caterpillar 793. Mass on the truck bodies was measured with strut pressure sensors on the haul trucks. According to the manufacturer, the weight results can vary depending on how well the scales are calibrated; however, the mass measurements should always be within 5% accuracy. The data acquisition was active for a total of 9 bucket fills.
Measurements for full-scale calibration of the granular material were performed. The static angle at which material rests after unloading at waste rock dumpsites was measured, and dynamically the dumping sequence was video recorded with the Mavic 2 drone. Photographs of the profiles of the waste rock piles were taken from a distance in order to capture the full pile. Figure 5 shows an example of one of the photographs. The photographs were further imported to the open-source software ImageJ to measure the angle. A total of 6 waste rock piles were photographed and analyzed. Photographs were taken approximately perpendicularly to the plane of the pile angle.

Numerical model and simulation setup
In the following sections, a description of the numerical methods and the simulation model of the rope shovel loading process setup is presented.

Granular material model
The granular material model was a combined finite element (FE) and discrete element (DE) method. For the larger particles, a type of finite element rigid body dynamics was used. Each particle is discretized with an arbitrary number of finite elements and hence the geometry is defined. Further, all of the elements within a rock are specified to be rigid itself and also rigid relative to each other. Solution of the rigid bodies is performed via rigid body dynamics equation of motion, based on the mass and inertia tensors from the element. Since all the elements are assumed to be rigid bodies, the material properties such as Young's modulus and Poisson's ratio are of no interest for the material itself, however the parameters are used in the contact algorithm (described in Sect. 3.1.1) and thus need to be correctly specified. This is the way rigid bodies are modeled in LS-DYNA, and the full algorithm is described in detail by Benson and Hallquist [1]. The benefit of utilizing this method compared to DEM with spheres was the possibility of modeling arbitrary shapes of the particles. However, since the geometrical description with several elements for a shape requires computational effort due to the larger amount of contact possibilities the traditional spherical DEM was utilized for the smallest rocks in the granular material model.

Discrete element method
The DEM was originally developed by Cundall and Strack, [13,14], and it is a discrete method where individual granular particles and their interactions are modeled. The interactions between the particles then govern the behavior of the bulk material when subjected to external forces. DEM has its basic foundation in the use of Newton's second law of motion to determine particle translations and rotations. The DEM treats each granular particle as a discrete and usually spherical element. Each element acts with its own individual size, shape, and stiffness properties, and the interaction with other discrete elements is governed by some contact algorithm.
In a DEM model, each particle i has mass m i , radius r i and velocity v i . Through Newton's second law of motion, the translation and rotational motion can be calculated. The translation of a particle is calculated as where F i c is the sum of forces acting on the particle due to normal and tangential contact with other particles and g is the gravitational acceleration. The rotational motion is calculated by where I i is the moment of inertia for the particle, ω i is the angular velocity, and M i the sum of torque acting on the particle.
Equations 1 and 2 constitute the governing equations for a DEM model and there are two unknown motion variables to solve for. The translation and rotation are calculated through numerical integration of the governing equations. In most cases, an explicit time integration is performed, and by such the time step size should be kept small enough to ensure numerical stability. Furthermore, the contact between particles happens on a certain time scale, and thus the time step needs to be kept small enough to detect and resolve the contacts with adequate accuracy. The general suggestion for selecting the DEM time step size is usually on the form Δt = c Δt √ m/k, where c Δt is a constant and k is the stiffness of the particle. The choice of the constant c Δt has been discussed in the literature and a summary can be found in [21].

Contact model for DEM
A common method for treating contact between particles or bounding surfaces, is to model the contact as a linear spring and dash-pot system. This model makes it possible to treat the particle-particle or particle-wall contact in the same way as for a damped harmonic oscillator. The contact force in the normal direction is calculated as [25,33] where K n is the spring normal stiffness, δ n is the contact overlap in the normal direction, v n is the relative velocity in normal direction of the two bodies in contact, and c n is a viscous normal damping coefficient reducing a specific part of the relative kinetic energy at contact. K n for a particleparticle contact is calculated as where k i , r i and C norm is the particle bulk modulus, particle radius and stiffness scaling coefficient in normal direction respectively. The bulk modulus is given as with Young's modulus E i and Poisson's ratio ν i . The tangential force is calculated in the same way as for the normal contact force in Eq. 3, however replacing the subscript n with t. The tangential spring stiffness is defined as K t = K n C tan where the constant C tan is a scaling coefficient. The viscous damping coefficient for both tangential and normal direction is given as, where D is the damping coefficient (0 ≤ D ≤ 1.0) and η crit the critical damping. Forces from friction between particles is modeled with Coulomb's friction law given as where μ fr and F n is the coefficient of friction and sum of forces in the normal direction [25]. The frictional forces govern the resistance between the particles in the tangential direction and thus is contributing to the tangential forces at contact. Rolling friction or rolling resistance is defined through introducing a momentum M r for the particle, that resists the rolling motion where μ roll is the rolling friction coefficient.

Contact model for FE rigid bodies
For the rigid finite element rocks, two types of penalty-based contacts for shell and solid finite elements implemented in LS-DYNA were used [32]. A segment-based contact, which means that the contact detects segment-segment penetration were used for the rigid FE rocks self contact and for rigid FE rocks and all other rigid FE bodies such as the bucket. The segment-segment contact stiffness for both normal and tangential direction is given by where C F E 1 is a penalty stiffness scale factor for the contact, m 1 and m 2 are the segment masses and Δt 0 is the time step. The contact is activated in LS-DYNA by using SOFT = 2. The major benefit of the segment-based contact is when surfaces are not smooth. The penalty contact model where the contact detects nodesegment penetrations were used for contact between DEspehers and rigid FE rocks. The contact stiffness for normal and tangential direction is calculated as with material bulk modulus k, penalty scale factor C F E 2 , segment area A and volume V . This contact model does not contain any rolling friction. It is activated in LS-DYNA by using SOFT = 0. The frictional force for the FE rigid parts are calculated by Coulomb's friction according to Eq. 7. Damping for the

Generating particles
In Table. 1, the particle size distribution used for the granular material model is presented. The PSD in Fig. 1. was used as a base for the data in Table 1. The PSD, given in terms of mass percentage was translated into number of particles percentages in order to implement it into the model. Particles smaller than 200 mm are neglected in the present study to maintain feasible computational times. The shape of the rocks was defined in the CAD software Siemens NX, by using the 3D-scanned point data of the rocks shown in Fig. 2. as a reference. Furthermore, the CAD model of the rocks was discretized with 4-noded solid tetrahedron finite elements and each rock was assumed rigid. In Fig. 6. the different geometries and sizes of the FE-rocks are seen. DEM-spheres were filled in the voids between the FE-rocks with the LS-PrePost packing routine to obtain the correct particle count according to Table 1.

Model calibration
The granular material model was calibrated through simulating a dynamic generation of a waste rock pile that should end up with a static representation of a waste rock pile where the angle can be measured and compared to experimental values. The dynamic scales occurring at the waste rock piles were assumed adequate for the loading process. In Fig. 7, the geometrical setup for the calibration simulation model is shown. A total of approximately 145,000 DEM particles and 300,000 finite elements distributed on 20,600 rigid FE particles were used. Consistent static pile angle θ and general dynamic behavior of the simulation model was obtained by adjusting the sliding friction, rolling friction, damping, and stiffness parameters between rigid FE-FE, DEM-FE, and DEM-DEM contacts in the model. Initially, the friction and damping coefficients were set to values chosen with judgment by the authors and further tuned in through iterating the simulations until simulated angle and general material flow behavior was consistent with the average angle of the experiments.
An explicit time integration scheme was used, and the stiffness controlling the contact of the ore was scaled to a lower value (with scale factors C norm , C tan , C F E 1 and C F E 2 ) in order to maintain a stable contact and at the same time obtain an adequate time step to keep reasonable computational times. Previous research shows that Young's modulus, which is related to the contact stiffness (see Eqs. 3, 4 and 5) may be scaled within a reasonable range without affecting the material bulk behavior [31,49], however in some cases were compression, large penetration and shearing behavior occurs, the authors suggest that users should keep in mind that bulk behavior may be affected.
In the calibration process, it was observed that the damping coefficients between all parts had to be large and that the contact stiffness had to be down-scaled in order to make a soft and smooth contact. A too stiff contact was leading to unstable contact behaviour. However, as the damping and stiffness coefficients were in a stable range their influence did not have a significant affect on the bulk behavior. In the model, the contact stiffness was scaled between 0.01 and 0.1 times the contact stiffness without observing any effect on the bulk behavior; however, when the stiffness was higher, the contact became unstable. When a stable contact was established, the calibration of the friction coefficients was performed. The first step was to calibrate the DEM-spheres and rigid FE rocks separately to keep the calibration parameters to a minimum. The rolling friction and the sliding friction coefficient for the DEM particles were adjusted until it matched the experiment. Since the real rocks had quite large aspect ratios and sharp corners the rolling friction was applied. The rolling resistance between the ground and the DE particles was sensitive, since the particles rolled away when more particles were added with a too small rolling resistance. The same procedure was then performed with the rigid FE particles. For the rigid FE particles, the calibration was quite simple since there are no rolling friction to consider due to the shapes. When both the FE and DE-particles were behaving relatively close to the experiment, they were combined in the model. In the last step where the full calibration model was used, the main adjustments were performed on the FE-DE friction coefficients in order to receive a correct behavior, even though some smaller adjustments were done on frictions between DE and FE.
Bulk density could not directly be given as an input to the numerical model, instead the particle density was given. Hence, a simulation was performed where a box with a defined volume of 28 m 3 was filled to the top with granular ore. Furthermore, the weight of the material inside of the box was measured in the model with the initial particle density to obtain the simulated bulk density. In order to receive the correct bulk density in the simulation, the initial particle density was scaled by the ratio between experimental and ini-  . 7 Illustration of the simulation domain for calibration. The grey area illustrates the granular material. The angled support sheet is moved upwards with a velocity v z = 0.61 m/s for a time of 20 s. The material is flowing as the sheet is moved upwards, and finally rests at a simulated pile angle of θ when the sheet stops moving. The initial height and depth of the domain pile was approximately 35 and 12 m respectively tial simulated bulk density. The scaled particle density was the same for both the FE and DEM particles.
Friction between ore material and bucket was measured from an experimental pin on disc tests between a Hardox plate (bucket material) and two rocks from the Aitik mine. The average of the tests was used as the friction coefficient between the bucket and the granular material in the model, see Table 2.

Rope shovel model
A CAD model was created with the basis of a 3D scanning of the Bucyrus rope shovel (Fig. 3). In the present study, the focus was on the bucket, and hence, a very coarse model of the rope shovel and a more detailed CAD model of the bucket was performed. Further, the CAD model was used to create an FE-mesh. The numerical representation of the Esco bucket is seen in Fig. 8, the bucket has a total width, depth, and height of approximately 4, 4, and 3 m respectively. Bucket fill volume capacity was given as 42 m 3 by the manufacturer. The rope shovel is numerically represented as a rigid multi-body system where each body consists of several finite elements. Since rigid bodies were used, the mesh size is of little interest for all parts except the bucket. Dynamic behavior is applied through penalty-based joints. In Sect. 2.2 the possible modes of movements were described, the numerical representation of the movements are summarized in Table 3. Prescribed motions were given as translation-time curves for the hoist ropes and crowd arm. Curves were extracted from the data acquisition and implemented in the model. Some of the waiting time occurring in experiments are removed from the curves to minimize simulation time, Fig. 9. Rotational velocity-time curves for the doors and upper body rotations were estimated by looking at the video recordings during the data acquisition. Totally, there were 14 joints in the model. A rigid constraint between crowd arm and bucket existed. The crawlers were constrained in all 6 possible degrees of freedom. Since the model does not include any bearings, lubrication, etc. which in reality contributes to damping in the system, artificial damping and stiffness were introduced in the joints as well as mass damping for some of the rigid parts. Joint stiffness parameters were scaled up to obtain a behavior similar to infinitely strong joints for the motors in order to follow the prescribed motion as closely as possible. The prismatic joint stiffness around its rotational axis was adjusted with reasonable values in order to obtain correct model behavior.
The two hoist ropes, connecting to the bucket at the lower ends and runs through the cylindrical joint at the boom top were modeled separately. Rigid beams are used for approximately 70 % of the rope length on the upper part where it passes through the top of the boom. The lower part of the ropes is modeled with short elastic beam elements. The elements are free to rotate in the end nodes and only takes force in tension, thus, making it behave like a realistic rope. In Fig. 10, an illustration of the rope model is presented. Rope modeling was applied in order to obtain similar behavior as was observed in the real loading process.

System model
A system simulation model was set up, consisting of the granular ore model (Sect. 3.1), and the multi-body dynamic rope shovel model (Sect. 3.2). In Table 2 the parameters used in the model are presented. Experimentally known values are particle and bulk densities of the ore material described in Sect. 2.1, and measured friction between granular material and bucket material. The Young's modulus and Poisson's ratio is set to 45 GPa and 0.3 respectively for all particles. The Young's modulus value is chosen without specific testing on the current rock types and is based on similar rock types. However, since the contact stiffness, which is dependent on the Young's modulus is down-scaled with approximately order of magnitude 1, the actual value of the modulus will not play a significant role. The ground is modeled as a FE-plane, hence it does not contain any surface roughness or unevenness which were seen in reality. For this reason, the friction coefficient was set to 0.9 between the ground and particles.
If not else specified in Table 2, the damping coefficient is 1 and the stiffness scale factor is 0.05. The remaining material properties were calibrated. In Fig. 11a the numerical setup of the system is shown. A 3D scan of the rope shovel and the muck pile was used in order to achieve the correct positions for the crowd arm, hoist ropes relative to the rope shovel and the whole rope shovel relative to the muck-pile. This positioning ensured that the model was geometrically set up in the same way as during the experiment. Since the muck pile was generated when there was no rope shovel present, the rocks occupying the space where the bucket was were manually moved slightly or removed in order to avoid initial penetrations. A rigid haul truck in approximately the same size as a CAT 793 was also included in the system model. In Fig. 11b, a model is presented where some of the particles were removed in the areas where the resultant velocity of particles was 0. Particles on a boundary of approximately 0.5 m were constrained in all directions. This simplification was done to reduce simulation time further and is later referred to as a simplified model. The simplified model was used to study different dig variables such as rake angle, position from the muck-pile, and dig velocity. In Fig. 11c, a profile of the 3D scanned muck-pile is positioned on top of the muck pile generated by cutting the final calibration model at a height of 10 m. The muck-pile was created by using the same model as the one used for calibration; however, to reduce the computational effort, the height of the muck-pile was reduced.
During the simulations, the contact between particles and bucket were activated during digging and dumping. The contact was deactivated when the rope shovel was moving back to the dig position. This was done in order to avoid unnecessary contact between bucket and muck pile to reduce computational time.

Results and discussion
In the following sections, the results from the numerical model are presented. First, the validity of the numerical model is checked by comparing experimental/real behavior and measurements with simulated results. Subsequently, the characterization of the bucket filling is performed, and finally, a number of simulations where three rope shovel parameters that are critical to the loading efficiency are studied in order to investigate how the loading process is affected by these. For all simulations, unless else stated, the time step was set to 4.5e-5 s, and simulations process time was approximately 140 seconds in order to capture 3 bucket fillings. Simulations were performed using the multi-physics software LS-DYNA (version MPP R.11.0) with 16 2.6 GHz CPUs. The simulation time for the full simulation model to CPU time was approximately 1:4190, meaning that 1 simulation time unit requires 4190 CPU time units. For the simplified model, the same ratio was 1:2590.

Calibration
Calibration of the granular ore model was performed according to the description in Sect. 3.1.5. Figure 12. shows 4  Table 4. The static pile angle agreed well with the average of the experimental pile angles. Measurements of the pile angles were however not a straight forward process, and the result varied slightly depending on how the angle was measured, hence the angle θ from experiments and simulations may vary ± 1 • . Fig. 12 Simulation results from the final calibration test described in 3.1.5. The figure shows snapshots from the calibration simulation model at the different times, t = 0 s, at start, t = 9 s, when the sheet has moved approximately halfway upwards, t = 18 s when the sheet is almost at a stop, (sheet stops at t = 20 s), and finally at t = 40 s the pile is at rest. Repose angle θ is measured to 34.8 • and the length, l and height, h are 20 and 16 m, respectively

Validity of numerical model
Validation of the numerical model was checked by comparing the results from the simulation with the experimental results.
The results are based on the loading of one truck body, which corresponds to three bucket fillings. In Table 5 a comparison between simulated and experimental results of ore mass in the bucket for each pass is presented. Results from both the full and simplified simulation models are seen. Simulated and experimental mass for the two first passes is differing with less than 3 % which are within the error tolerance that the manufacturer of the measuring equipment gives. The third pass gives a somewhat larger difference of around 8%. A possible reason for the difference might be because of several factors. The mass from simulations is measured when the ore is in the bucket, whereas the experimental mass is recorded on the truck body. From the drone recording, it was observed that the amount of spillage from the truck body for the last pass was significantly larger than for the first and second pass. Another factor that might give this deviation could be due to the muck-pile shape. Since the initial shape of the muck-pile is only known for the first pass, the change of the muck-pile shape is only governed by the granular ore model. Since the granular ore model was cut at a height of approximately 10 m, the "refilling" of material down into the muck-pile might vary. It was clear that the experimental muck-pile was higher, but the area where the actual digging occurs (along the bucket trajectory) was not as high, hence the muck-pile in the model was cut. Most of the material in the experimental muck-pile that is on a higher level than Results from the full and the simplified simulation model is presented. Values inside parenthesis corresponds to the percent error between simulated and experimentally measured mass 10 m is very saturated/compressed due to that the material has been blasted and further compressed after the material has been free-falling, also freezing temperatures could have affected this. During the experiment, the material collapsed from the muck-pile in batches more than a continuous flow.
These observations were the cause of the simplified muckpile, but may however affect the refilling of material into the muck-pile. Another explanation might be due to the PSD in the model. Particles smaller than 200 mm are not included in this model, and if the number of fines are smaller in the last pass, the model cannot recognize this, and i.e. give a larger difference in mass. A larger amount of fines will decrease the total air void ratio and cause a larger mass in the bucket. From the experimental drone recording, a comparison of the bucket filling sequence for all the passes was done. Figure 13a, b shows a sequence of snapshots from experiments and simulations for the first pass. The drone recording was captured at an angle, giving that the whole bucket was not visible at all times. However, when studying the video recording from simulation, the granular ore model behaved very similar to the behavior observed in experiments. Smaller particles were as mentioned previously ignored in the simulation and flow from fines was not detected in the model. However, the bulk behavior did not illustrate any major difference between experiments and simulations. Resultant velocity field of the particles seen in Fig. 13c shows that the particles near the boundaries are 0 or close to zero which shows that the assumption of the width of the muck-pile was reasonable.
Considering the results from mass measurements in Table  5 and the granular material flow comparison in Fig. 13 the numerical results shows good agreement with the experimentally observed ore material behavior during the loading process.

Bucket filling analysis
The particle flow mechanism into the bucket was experimentally studied by Rasimarzabadi [40] by using a scaled version of a transparent bucket. From the experimental observation, a five-stage particle flow pattern was proposed. The observations were also captured in the numerical model proposed in this paper. For definitions of the different terms that will be used in this section, see Fig. 4. In Fig. 14, 6 snapshots from simulations are presented to describe the filling process, in accordance with the 5-stage flow pattern described in [40]. At t = 0 s in Fig. 14, the bucket was positioned in front of the muck-pile and was right before it starts to move. At t = 2 s; stage 1, the front lip of the bucket engaged the material and was cutting the muck-pile. Particles above the front lip were seen entering the very front of the bucket sliding along the bucket and rolling into the bucket. Rasimarzabadi refereed to the particles in a rotational and sliding motion as the "flow zone". As the bucket engaged the muck-pile further, it approached stage 2, at t = 5 s. The majority of this stage was to develop a thicker "flow zone" that went along the whole width of the bucket. Particles were rolling and sliding over one another and along the length of the bucket. Stage 3, occurred approximately at t = 5 s. In stage 3, the flow zone reached the rear of the bucket and was increasing in thickness. After some thickness growth, the particles closest to the bottom would interlock with each other and along the bottom. Hence, particles lost their sliding/rotational motion and i.e. the flow zone, and entered the so-called bulk zone, where particles had 0 or a very small relative velocity. The "bulk zone" was located underneath the flow zone. The "flow zone" was still active as a layer of material under rotational and sliding movement on the top. As the bucket was moving along its trajectory, the bulk zone was increasing in thickness due to the interlocking with the flow zone layer. Due to gravity, the particles were moved backward, to the rear of the bucket. In stage 4, occurring at approximately t = 9 s, the bucket was further penetrating the muck-pile. In the front of the bucket, a large build-up of material was observed. Plenty of deadload was pushed upward due to the build-up of material in the front, even though there was a fair bit of empty volume towards the rear of the bucket. Stage 4 was directly dependent on the depth of cut, i.e. how far the bucket penetrates into the muck-pile. For some depths of cut, the bucket is filled to the top, which creates a bulk zone of all the material inside the bucket. In the case of the results seen from these simulations, this did not occur since the bucket was not filled to the top at this stage. This showed that a relatively small depth of cut was used. This leads to a fairly large flow zone trough out the whole bucket filling. At t = 13-20 s, stage 5 is taking place, which is the final step of the bucket filling process. The bucket was moving out from the muck-pile, causing the Column (c) shows a fringe plot illustrating the resultant velocity of the particles, with units m/s. For the fringe plot sequence, the whole muck-pile is visible in order to see how the bulk material is affected by the rigid walls. Each row represents a specific timestamp, starting from 0, which corresponds to the time right before the bucket starts to move. Time units are in seconds piled up material in front of the bucket to fall away from the bucket, and a pile like shape, or heaped material was seen at the highest point of the material inside of the bucket. As the bucket got ready for dumping the material into the truck, it changed direction and was subjected to gravity. The material was then moved backward to the rear of the bucket, which is seen at t = 20 s.

Rope shovel parameter influence on bucket filling
A number of simulations were performed in order to detect trends and understand how different rope shovel parameters affect the loading process. Three rope shovel loading parameters that are of interest for loading efficiency were investigated, the rope shovel rake angle, position of the rope shovel, and the dig velocity. The evaluation was performed by measuring the mass of material inside the bucket and the digging force, while different parameters were varied. Dig force was calculated as the resultant of the force from the hoist and the crowd arm the crowd and hoist forces were recorded from the simulation models. For all the simulations in this section, the first pass of the simplified simulation model was utilized as a baseline case in order to maintain a relatively fast computational cycle.

Rake angle
Bucket rake angle, defined in Sect. 2.2 was varied. 5 different angles were chosen. The baseline case had a rake angle of 61.8 • . From the rake angle of the baseline case, it was varied  Fig. 15 the rake angle, material mass in the bucket, and the mean digging force are presented. In Fig.  16 the material inside the bucket is shown for the different rake angles. It is clear from the figures that a slight change in the rake angle significantly affects both the digging force required, and also how well the bucket was filled. Rake angle dependency seems to follow a trend of increased mass gives increase dig force, which seems logical. However, at the rake angle of 62.8 • , the average dig force was smaller than for the rake angle of 61.8 • , even though the mass in the bucket was larger. The following finding shows that there is an optimum rake angle (with respect to average dig force and mass in the bucket) found somewhere around 62.8 • , meaning that the bucket is penetrating the muck-pile in a feasible way, where the bucket is filled sufficiently. From the rake angle variation, it was found, that a change of the rake angle 2 • gives a change of mass in the bucket of almost 4 tons. Being able to fill the bucket with 4 extra tons every pass leads to savings of up to millions of tons on the truck bodies annually. However, it is important to keep in mind that other factors such as wear on the outside and bottom of the bucket may increase with a larger rake angle, and for a too small rake angle, the wear on the teeth and front lip may be excessive. Since wear usually occurs due to resistive forces, the average dig force should give an indication of the resistance of the bucket, nevertheless, the aspect of wear needs to be studied more thoroughly.

Rope shovel position
The position of the bucket (and rope shovel) relative to the muck-pile was investigated. 6 different positions were simulated and analyzed. The position change is only in one  Fig. 11c. The baseline case, described as position 0 m from the muck-pile was the same as for the simplified model in Sect. 4.2. In Fig. 17, the average force and mass in the bucket from the simulations are shown. In Fig. 18, the material inside the bucket is shown for the different positions.
There is a clear relationship between the position, mass of ore in the bucket, and the average dig force, The position from the muck-pile will affect the depth of cut, i.e., how far the bucket reaches into the muck-pile. The effect of position seems to be relatively sensitive, a change of 0.2 m gives a large change of mass in the bucket. Since the change in position gives such a huge impact, it might be one of the easiest ways of receiving a higher filling degree in the bucket. These results might sound intuitive, nevertheless, it requires an experienced operator to

Dig velocity
The bucket velocity as it travel through the muck-pile was varied. Velocity from the baseline case was defined to have the velocity ratio 1. 4 other velocities were varied from 0.7 to 1.3 times the baseline velocity. Scaling of the velocity was performed uniformly on both the hoist and crowd to maintain identical bucket trajectories. In Fig. 19 mass of material in the bucket, and mean dig force is presented as a function of the velocity ratio. The average dig force from simulations was taken as the average during approximately the same period, however, the time during that period varied due to the velocity difference. The period was from the start of the dig until the bucket was located above the dump truck. The dig veloc-  Fig. 19 Mass in bucket and average dig force as a function of the rake angle ity showed no clear trends on the mass in the bucket or the average dig force. Smaller fluctuations in mass and dig force were observed, and in general, a slightly better performance was observed for higher velocities, due to the larger mass in the bucket and smaller average dig force. Nonetheless, the presented results are of high value for the mining industry, since the mass in the bucket was kept approximately constant when an increased dig velocity was applied. Dig velocity has been studied and the result from this work seems consistent with the previous work, [9,18,40]. The somewhat constant mass in the bucket, with respect to digging velocity indicates that the dig cycle can be reduced, and in turn cause an increased dig capacity for the rope shovel.

General discussion
A full-scale system simulation model of the loading process for a heavy-duty mining rope shovel is developed. The granular material is modeled using a unique combination of mixing particles with realistic shapes through a rigid Fig. 18 Snapshot from simulations after the bucket was filled for the studied positions from the muck pile x. Bucket is transparent in order to make the material inside visible FE-discretization and traditional spherical discrete element method (DEM) particles. The system model and specifically the granular ore model is calibrated and validated with experimental measurements of mass, general material flow behavior, and static repose angles. The presented numerical model is used to obtain a greater understanding and characterize the bucket filling. The presented numerical model shows the possibility of modeling a large-scale industrial process. In this study, the main objective is to capture the material bulk behavior. Therefore, some of the simplifications and assumptions presented previously were performed.
Since most of the particles (in terms of mass percentage) in the Aitik mine are large > 200 mm, the focus was mostly on capturing the behavior of the large particles. The larger particles, due to rock fragmentation contain very sharp corners and aspect ratios, and modeling these types of rocks with traditional spherical DEM would be unrealistic. The trade-off between computational time, available methods, and considering the objective with this study was the choice for the granular model in this work. The result from the model shows that the presented modeling technique can accurately capture the main mass flow behavior. However, to include the smaller particles (< 200 mm) would be an improvement action for the model, and will extend the model to other applications as well. For example, for conveyor transportation after the ore is crushed. In most mines or other processes where rocks can be found, a method to utilize complex rock shapes has been a research topic for some years. Though some methods exist, see Sect. 1, the utilization of such methods requires some specific software or writing an in-house solver. Since the application of this work is aimed directly to the mining industry, the demand is to use a commercial code. The results presented in Sect. 4.4, shows three different parameters influencing the loading process. These three parameters are of high interest for the mining industry, used in this study to show the capability of the model to capture these relatively small changes. However, the use of the model in the mine is very broad. For example, the utilizing of a new bucket or tooth system investment, different bucket trajectories, or different muck-pile formations can be studied and evaluated in order to have a tool to make decisions with.
In mines, there are a lot of other mining machines and equipment involved in the mining process than rope shovels, for example, wheel loaders, hydraulic mining shovels, bulldozers, and haul dump trucks that are all dealing with the same particles as modeled in this work. This signifies that the proposed model can easily be utilized for those applications.
Since only one type of experiment (pile angle) was used, the set of parameters for friction coefficients, rolling friction coefficient, stiffness scale factors and damping coefficients are not unique. Previous research has found that in most of the experimental tests used for calibrating DEM parameters, the same type of bulk behavior for a specific test may be estab-lished with the use of different DEM parameters, [41]. For example, two sets of rolling and sliding friction coefficients may lead to the same angle of repose.
Since the structure is modeled with rigid finite elements, the possibility of performing stress and fatigue analysis by changing the rigid elements to e.g elasto-plastic, with a very good load case is another model utilization. As mentioned in the introduction, wear and damage predictions on the bucket and other equipment are subjects that are of high interest to the mining industry. Utilizing this numerical model for studying the interaction between bucket and ore in order to minimize the wear behavior is a topic for future work.

Conclusion
-The numerical model can capture the main phenomenons of the granular ore and the interaction with the bucket for the loading process with a rope shovel. Validation of the model shows that the model is representing the experimentally observed data excellent. -New large scale experimental methods for model calibration and validation of granular material behavior for mining industry is presented. -The model is used to evaluate the performance of the digging process when under certain conditions, e.g. during different rake angles, bucket positions, and dig velocities with interesting and valuable results as a product. -The results from the presented paper agrees with previous experimental studies on rope shovel bucket filling, and also extends the knowledge and gives suggestions on how a more efficient digging can be obtained. -This numerical model, developed in the present work, is not limited to the use of rope shovels, and may easily be applied where granular material with complex shapes interacts with equipment.
tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.