In this PhD some key computational and algorithmic aspects of the Combined Finite Discrete Element Method (FDEM) are critically evaluated and either alternative novel or improved solutions have been proposed, developed and tested. In particular, two novel algorithms for contact detection have been developed. Also a comparative study of different contact detection algorithms has been made. The scope of this work also included large and grand scale FDEM problems that require intensive use of CPU; thus, novel parallelization solutions for grand scale FDEM problems have been developed and implemented using the MPI (Message Passing Interface) based domain decomposition. In this context a special attention is paid to the rapidly developing multi-core desktop architectures. The proposed novel solutions have been intensively validated and verified and demonstrated using various problems from literature.

Many powders and particulate solids are stored and handled in large quantities across various industries. These solids often encounter handling and storage difficulties that are caused by the material cohesion. The cohesive strength of a bulk material is a function of its past consolidation stress. For example, high material cohesive strength as a result from high storage stresses in a silo can cause ratholing problems during discharge. Therefore, it is essential to consider the stress-history dependence when evaluating such handling behaviour. In recent years the Discrete Element Method (DEM) has been used extensively to study the complex behaviour of granular materials. Whilst extensive DEM studies have been performed on cohesionless solids, much less work exists on modelling of cohesive solids. The commonly used DEM models to model adhesion such as the JKR, DMT and linear cohesion models have been shown to have difficulty in predicting the stress-history dependent behaviour for cohesive solids. DEM modelling of cohesive solid at individual particle level is very challenging. To apply the model at single particle level accurately would require one to determine the model parameters at particle level and consider the enormous complexity of interfacial interaction. Additionally it is computationally prohibitive to model each and every individual particle and cohesion arising from several different phenomena. In this study an adhesive elasto-plastic contact model for the mesoscopic discrete element method (DEM) with three dimensional non-spherical particles is proposed with the aim of achieving quantitative predictions of cohesive powder flowability. Simulations have been performed for uniaxial consolidation followed by unconfined compression to failure using this model. Additionally, the scaling laws necessary to produce scale independent predictions for cohesionless and cohesive solids was also investigated. The influence of DEM input parameters and model implementation have been explored to study the effect of particle (meso-scale) properties on the bulk behaviour in uniaxial test simulation. The DEM model calibration was achieved using the Edinburgh Powder Tester (EPT) – an extended uniaxial tester to measure flowability of bulk solids. The EPT produced highly repeatable flowability measurements and was shown to be a good candidate for DEM model calibration. The implemented contact model has been shown to be capable of predicting the experimental flow function (unconfined compressive strength versus the prior consolidation stress) for a limestone powder which has been selected as a reference solid in the Europe wide PARDEM research network. Contact plasticity in the model is shown to affect the flowability significantly and is thus essential for producing satisfactory computations of the behaviour of a cohesive granular material. The model predicted a linear relationship between a normalized unconfined compressive strength and the product of coordination number and solid fraction. Significantly, it has been found that contribution of adhesive force to the limiting friction has a significant effect on bulk unconfined strength. Failure to include the adhesive contribution in the calculation of the frictional resistance may lead to under-prediction of unconfined strength and incorrect failure mode. The results provide new insights and propose a micromechanical based measure for characterising the strength and flowability of cohesive granular materials. Scaling of DEM input parameters in a 3D simulation of the loading regimes in a uniaxial test indicated that whilst both normal and tangential contact stiffness (loading, unloading, and load dependent) scales linearly with radius of the particle, the adhesive forces scales with the square of the radius of the particles. This is a first step towards a mesoscopic representation of a cohesive powder that is phenomenological based to produce the key bulk characteristics of a granular solid and the results indicate that it has potential to gain considerable computational advantage for large scale DEM simulations. The contact model parameters explored include particle contact normal loading stiffness, tangential stiffness, and contact friction coefficient. The DEM model implementation parameters included numerical time step, strain rate, and boundary condition. Many useful observations have been made with significant implications for the relative importance of the DEM input parameters. Finally the calibration procedure was applied to a spray dried detergent powder and the simulation results are compared to whole spectrum of loading regime in a uniaxial experiment. The experimental and simulation results were found to be in reasonable agreement for the flow function and compression behaviour.

Granular materials are in abundance in nature and are estimated to constitute over 75% of all raw materials passing through the industry. Granular or particulate solids are thus of considerable interest to many industrial sectors and research communities, where many unsolved challenges still remain. This thesis investigates the micro- and macro-phenomena in densely packed particulate systems by means of the Discrete Element Method (DEM), which is a numerical tool for analysing the internal complexities of granular material as the mechanical interactions are considered at the grain scale. It presents an alternative approach to phenomenological continuum approaches when studying localisation problems and finite deformation problems in granular materials. In order to develop a comprehensive theoretical understanding of particulate matter and to form a sound base to improve industrial processes, it is desirable to study the mechanical behaviour of granular solids subject to a variety of loading conditions. In this thesis, three loading actions were explored in detail, which are biaxial compression, rigid object penetration and progressive formation of granular piles. The roles of particle shape and contact friction in each of these loading scenarios were investigated. The resulting packing structures were compared and studied to provide a micromechanical insight into the development of contact force network which governs the collective response. The interparticle contact forces and displacements were then used to evaluate the equivalent continuum stress and strain components thus providing the link between micro- and macroscopic descriptions. The information collected from the evolution of strong contact network illustrates the underlying mechanism of force transmission and propagation. DEM simulations presented in this thesis demonstrate strong capability in predicting the bulk behaviour as well as capturing local phenomenon occurring in the system. The research first simulates a testing environment of biaxial compression in DEM, in which the phenomenon of strain localisation was investigated, with special attention given to the interpretation of underlying failure mechanism. Several key micromechanical quantities of interest were extracted to understand the bifurcation instability, such as force chains, contact orientation, particle rotation and void ratio. In the simulation of progressive formation of granular piles, a counterintuitive pressure profile with a significant pressure dip under the apex was predicted for three models under certain conditions. Both particle shape and preparation history were shown to be important in the resulting pressure distribution. During the rigid body penetration into a granular sample, the contact forces were used to evaluate the equivalent continuum stress components. Significant stress concentration was developed around the punch base which further led to successive collapse and reformation of force chains. Taking the advantage of micromechanical analysis at particle scale, two distinct bearing failure mechanisms were identified as the penetration proceeded. To further quantify the nature of strain mobilisation leading to failure, Particle Image Velocimetry (PIV) was employed to measure the deformation over small strain interval in association with shear band propagation in the biaxial test and deformation pattern in the footing test. The captured images from DEM simulation and laboratory experiments were evaluated through PIV correlation. This optical measuring technique is able to yield a significant improvement in the accuracy and spatial resolution of the displacement field over highly strained and localised regions. Finally, a series of equivalent DEM simulations were also conducted and compared with the physical footing experiments, with the objective of evaluating the capability of DEM in producing satisfactory predictions.

From pharmaceutical powders to agricultural grains, a great proportion of the materials handled in industrial situations are granular or particulate in nature. The variety of stesses that the matierals may experience and the resulting bulk behaviours may be complex. In agricultural engineering, a better understanding into agricultural processes such as seeding, harvesting, transporting and storing will help to improve the handling of agricultural grains with optimised solutions. A detailed understanding of a granular system is crucial when attempting to model a system, whether it is on a micro (particle) or macro (bulk) scale. As numerical capabilities are ever increasing, the Discrete Element Method (DEM) is becoming an increasingly popular numerical technique for computing the behaviour of discrete particels for both industrial and scientific applications. A look into the literature shows a lack of validation of what DEM can predict, specifically with respect to bulk behaviour. In addition, when validation studies are conducted, discrepancies between bulk responses in physical tests and numerical predictions using measured particles properties may arise. The aire of this research is to develop a methodology to calibrate DEM models for agricultural grains using data meaured in bulk physical tests. The methodology will have a wider application to granular solids in general and will advance understanding in the area of DEM model calibration. A contrasting set of granular materials were used to develop the methodology including 3 inorganic solids (single and paired glass beads, and polyethylene terephthalate pellets) and two organic materials (black eyes beans and black kidney beans). The developed methodology consists of three steps: 1. The development of bulk physical tests to measure the bulk responses that will be used to calibrate the DEM models, 2. The creation of the numerical dataset that will describe how the DEM input parameters influence the bulk responses , and 3. The optimisation of the DEM parameters using a searching algorithm and the results from Step 1 and 2. Two laboratory devices were developed to provide calibration data for the proposed methodology: a rotating drum and an confined compression test. These devices were chosen as they can produce bulk responses that are repeatable and easy to quantify, as well as generate discriminating results in numerical simulations when DEM parameters are varied. The bulk response determined from the rotating drum device was the dynamic angle of repose Ør formed when the granular material in a 40% filled drum is rotating at a speed of 7 rpm. the confined compression apparatus was used to determine the bulk stiffness of a system by monitoring the change in void ratio from the stress applied during a loading and unloading cycle. The gradient of the loading and unloadng curves termed λ and κ respectively were chosen as the bulk responses to calibrate the DEM models. The experimental results revealed that the dynamic Ør was significantly influences by the particle aspect ration and boundary conditions. The stiffness parameters were found to be predominantly influences by the initial packing arrangement. The numerical dataset describing how the DEM input parameters influence the numerical bulk responses was created by simulating the bulk physical tests, varying selected DEM parameters and monitoring the effects on bulk parameters. To limit the number of simulations required, design of experiment (DOE) methods were used to determine a reduced factorial matrix of simulations. In additions, an extensive parametric investigation on the non-optimised parameters as well as a scaling sensitivity study was carried out. The final step in determining the optimised parameters is to use a searching algorithm to infer the DEM parameters based on the numerical dataset and used the experimental results as calibration data. To perform a comparative study, tow searching algorithms were explored: the first was a simple method based on Microsoft Excel's Solver algorithm coupled with a weighted inverse distance method. The second made used of the statistical analysis program Statistica. It was shown that the Excel Solver algorithm is simpler and quicker to use but for the present first implementation, could only perform an optimisation based on two bulk responses. Statistica required the creation of a staistical model based on the numerical dataset before using the profiling and desirability searching technique, but was able to optimise the parameter using all three bulk responses. A verification and validation of the optimisation methodology was conducted using the optimised parameters for the black eyed beans. A verification was cnducted by simulating the two calibration experiments using the optimsed parameters and comparing these with the experiments. In addition, a validation was peformed by predicting the response of ta shallow footing penetration on a bed of black eyed beans. It was found that DEM simulations using optimised parameters predicted vertical stress on the footing during penetration to an acceptable degree of accuracy for industrial applications (<10%) at penetration depths up to 30mm.

The main objective of this research was to model the airflow paths through grain bulks and predict the resistance to airflow. The discrete element method (DEM) was used to simulate the pore structures of grain bulks. A commercial software package PFC3D (Particle Flow Code in Three Dimension) was used to develop the DEM model. In the model, soybeans kernels were considered as spherical particles. Based on simulated positions (coordinates) and radii of individual particles, the characteristics of airflow paths (path width, tortuosity, turning angles, etc.) in the vertical and horizontal directions of the grain bed were calculated and compared. The discrete element method was also used to simulate particle packing in porous beds subjected to vertical vibration. Based on the simulated spatial arrangement of particles, the effect of vibration on critical pore structure parameters (porosity, tortuosity, pore throat width) was quantified. A pore-scale flow branching model was developed to predict the resistance to airflow through the grain bulks. Delaunay tessellation was also used to develop a pore network model to predict airflow resistance. Experiments were conducted to measure the resistance to airflow to validate the models. It was found that the discrete element models developed using PFC3D was capable of predicting the pore structures of grain bulks, which provided a base for geometrically constructing airflow paths through the pore space between particles. The tortuosity for the widest and narrowest airflow paths predicted based on the discrete element model was in good agreement with the experimental data reported in the literature. Both pore-scale models (branched path and network) proposed in this study for predicting airflow resistance (pressure drop) through grain bulks appeared promising. The predicted pressure drop by the branched path model was slightly (<12%) lower than the experimental value, but almost identical to that recommended by ASABE Standard. The predicted pressure drop by the network model was also lower than the measured value (2.20 vs. 2.44 Pa), but very close to that recommended by ASABE Standard (2.20 vs. 2.28Pa). / February 2017

Soil disturbance and cutting force are two of the most common performance indicators for soil-engaging tools. In this study the interaction of two soil-engaging tools (a disc opener for fertilizer banding and a hoe opener from an air drill) with soil were modeled using Particle Flow Code in Three Dimensions (PFC3D), a discrete element modeling software. When comparing the disc model to the experiment results, the relative error was 11% for the average soil throw, 1.9% for the average draft force, and 51% for the average vertical force. Results from the soil-hoe model showed a relative error of 15% between the simulated soil throw and the measured one. In conclusion, both the soil-disc and soil-hoe models could simulate the selected soil dynamic properties (except for the vertical forces of the disc opener) with a reasonably good accuracy, considering the highly variable nature of the soil. / October 2016

Aggregate strength, gradation, and shape play a vital role in controlling asphalt mixture performance. Many studies have demonstrated the effects of these factors on asphalt mixture performance in terms of resistance to fatigue cracking and rutting. This study introduces numerical and analytical approaches supported with imaging techniques for studying the interrelated effects of aggregate strength, gradation, and shape on resistance of asphalt mixtures to fracture. The numerical approach relies on the discrete element method (DEM). The main advantage of this approach is the ability to account for the interaction between the internal structure distribution and aggregate properties in the analysis of asphalt mixture response and performance. The analytical approach combines aggregate strength variability and internal force distribution in an asphalt mixture to predict the probability of aggregate fracture. The numerical and analytical approaches were calibrated and verified using laboratory tests on various aggregate types and mixtures. Consequently these approaches were used to: (1) determine the resistance of various mixture types with different aggregate properties to fracture, (2) study the effects of aggregate strength variability on fracture, (3) quantify the influence of blending different types of aggregate on mixture strength, (4) develop a mathematical expression for calculating the probability of aggregate fracture within asphalt mixture, and (5) relate cracking patterns (cohesive: aggregate - aggregate and matrix - matrix, and adhesive: aggregate - matrix) in an asphalt mixture to internal structure distribution and aggregate properties. The results of this dissertation established numerical and analytical techniques that are useful for developing a virtual testing environment of asphalt mixtures. Such a virtual testing environment would be capable of relating the microscopic response of asphalt mixtures to the properties of the mixture constituents and internal structure distribution. The virtual testing environment would be an inexpensive mean to evaluate the influence of changing different material and design factors on the mixture response.

Aggregate strength, gradation, and shape play a vital role in controlling asphalt mixture performance. Many studies have demonstrated the effects of these factors on asphalt mixture performance in terms of resistance to fatigue cracking and rutting. This study introduces numerical and analytical approaches supported with imaging techniques for studying the interrelated effects of aggregate strength, gradation, and shape on resistance of asphalt mixtures to fracture. The numerical approach relies on the discrete element method (DEM). The main advantage of this approach is the ability to account for the interaction between the internal structure distribution and aggregate properties in the analysis of asphalt mixture response and performance. The analytical approach combines aggregate strength variability and internal force distribution in an asphalt mixture to predict the probability of aggregate fracture. The numerical and analytical approaches were calibrated and verified using laboratory tests on various aggregate types and mixtures. Consequently these approaches were used to: (1) determine the resistance of various mixture types with different aggregate properties to fracture, (2) study the effects of aggregate strength variability on fracture, (3) quantify the influence of blending different types of aggregate on mixture strength, (4) develop a mathematical expression for calculating the probability of aggregate fracture within asphalt mixture, and (5) relate cracking patterns (cohesive: aggregate - aggregate and matrix - matrix, and adhesive: aggregate - matrix) in an asphalt mixture to internal structure distribution and aggregate properties. The results of this dissertation established numerical and analytical techniques that are useful for developing a virtual testing environment of asphalt mixtures. Such a virtual testing environment would be capable of relating the microscopic response of asphalt mixtures to the properties of the mixture constituents and internal structure distribution. The virtual testing environment would be an inexpensive mean to evaluate the influence of changing different material and design factors on the mixture response.

The main objective of the project was to develop a model for the motion of granular media under vibration in a tub vibrator. For such a system, it was decided that a discrete element method (DEM) was the most appropriate tool to model bulk velocity and circulation of media. In the first phase of the work, a vibratory finisher was modified to introduce planar vibration into a single layer of particles. The motion of the tub was measured using accelerometers and the corresponding granular media behavior was determined by video recording. A discrete element model, based on Cundall’s approach to contact, was developed to model granular flow in different vibratory beds, and the results were compared with experimental measurements of bulk flow velocity and bed expansion for the tub finisher. The sensitivity of the model predictions to the contact parameters was considered and the parameters were optimized with respect to the experimental results. After optimization, the difference between the model predictions of the bulk flow velocity and the measurements was less than 20% at four locations in media beds of two depths. The average bulk density of the vibrating beds was also predicted to be within 20% of the measured values. In the next phase, a two-dimensional discrete element model was developed to model single-cell circulation in vibratory beds that had both vertical and horizontal components of motion. The model predictions were compared with experimental measurements of the onset and growth of circulation in beds of steel and glass spheres as a function of bed depth, inter-particle and wall friction coefficients, and the amplitude of vibration. While the values from the DEM showed an error of up to 50% in the predicted circulation strength, depending on the type of the media and system conditions, the trends predicted by the model closely matched those in the experiments. Finally, a physical model was developed to describe the relationship between the onset and direction of circulation with the vibration of the container. A similar model was used to describe the experimental results as well as the transition in circulation patterns in terms of the resultant shear forces at the vibrating container walls and the interlocking of media close to the container walls. It was also demonstrated that a two-dimensional DEM could model a granular flow in which the media had three-dimensional contact and freedom of movement, but that was driven by vibrations in a plane. In summary, it was found that the linear optimization procedure for the contact parameters is an efficient way to improve the results from DEM. Additionally, the circulation in a tub-vibrator increased with the depth of the particulate media in the container, and with the magnitude of the wall-particle and particle-particle friction coefficients. The strength of circulation also increased with the amplitude of vibration. A strong correlation existed between the total shear force along the vibrating container walls and the circulation behavior. Bulk circulation increased sharply when increasing bed depth increased the pressure and the shear forces at the walls and between particle layers. It was also concluded that dimensionless bed depth (the ratio of bed depth to particle diameter) was not a proper dimensionless group when discussing the circulation behavior and it should act in conjunction with other parameters.

The main objective of the project was to develop a model for the motion of granular media under vibration in a tub vibrator. For such a system, it was decided that a discrete element method (DEM) was the most appropriate tool to model bulk velocity and circulation of media. In the first phase of the work, a vibratory finisher was modified to introduce planar vibration into a single layer of particles. The motion of the tub was measured using accelerometers and the corresponding granular media behavior was determined by video recording. A discrete element model, based on Cundall’s approach to contact, was developed to model granular flow in different vibratory beds, and the results were compared with experimental measurements of bulk flow velocity and bed expansion for the tub finisher. The sensitivity of the model predictions to the contact parameters was considered and the parameters were optimized with respect to the experimental results. After optimization, the difference between the model predictions of the bulk flow velocity and the measurements was less than 20% at four locations in media beds of two depths. The average bulk density of the vibrating beds was also predicted to be within 20% of the measured values. In the next phase, a two-dimensional discrete element model was developed to model single-cell circulation in vibratory beds that had both vertical and horizontal components of motion. The model predictions were compared with experimental measurements of the onset and growth of circulation in beds of steel and glass spheres as a function of bed depth, inter-particle and wall friction coefficients, and the amplitude of vibration. While the values from the DEM showed an error of up to 50% in the predicted circulation strength, depending on the type of the media and system conditions, the trends predicted by the model closely matched those in the experiments. Finally, a physical model was developed to describe the relationship between the onset and direction of circulation with the vibration of the container. A similar model was used to describe the experimental results as well as the transition in circulation patterns in terms of the resultant shear forces at the vibrating container walls and the interlocking of media close to the container walls. It was also demonstrated that a two-dimensional DEM could model a granular flow in which the media had three-dimensional contact and freedom of movement, but that was driven by vibrations in a plane. In summary, it was found that the linear optimization procedure for the contact parameters is an efficient way to improve the results from DEM. Additionally, the circulation in a tub-vibrator increased with the depth of the particulate media in the container, and with the magnitude of the wall-particle and particle-particle friction coefficients. The strength of circulation also increased with the amplitude of vibration. A strong correlation existed between the total shear force along the vibrating container walls and the circulation behavior. Bulk circulation increased sharply when increasing bed depth increased the pressure and the shear forces at the walls and between particle layers. It was also concluded that dimensionless bed depth (the ratio of bed depth to particle diameter) was not a proper dimensionless group when discussing the circulation behavior and it should act in conjunction with other parameters.