The particle finite element method for transient granular material flow: modelling and validation

The prediction of transient granular material flow is of fundamental industrial importance. The potential of using numerical methods in system design for increasing the operating efficiency of industrial processes involving granular material flow is huge. In the present study, a numerical tool for modelling dense transient granular material flow is presented and validated against experiments. The granular materials are modelled as continuous materials using two different constitutive models. The choice of constitutive models is made with the aim to predict the mechanical behaviour of a granular material during the transition from stationary to flowing and back to stationary state. The particle finite element method (PFEM) is employed as a numerical tool to simulate the transient granular material flow. Use of the PFEM enables a robust treatment of large deformations and free surfaces. The fundamental problem of collapsing rectangular columns of granular material is studied experimentally employing a novel approach for in-plane velocity measurements by digital image correlation. The proposed numerical model is used to simulate the experimentally studied column collapses. The model prediction of the in-plane velocity field during the collapse agrees well with experiments.


Introduction
A common aspect of various industrial processes and natural phenomena is the flow behaviour of dense granular materials. The lack of comprehensive theoretical models results in a low operating efficiency of industrial processes including dense granular material flow. A granular material is composed of a large number of individual particles of arbitrary size and shape. Although the individual particles may be of relatively simple geometrical shape, granular materials features a wide range of complex behaviours. The mechanical behaviour of a granular material is strongly dependent on the loading conditions. For quasi-static loading conditions, the behaviour is solid-like, while the behaviour of a flowing granular material typically is liquid-like [35]. The study of B Simon Larsson simon.larsson@ltu.se 1 Division of Mechanics of Solid Materials, Luleå University of Technology, SE-97187 Luleå, Sweden 2 Mechanical Engineering Department, EAFIT University, Medellín, Colombia granular material flow is of importance in many industries, such as the mining industry, the pharmaceutical industry and the agricultural industry. Numerical modelling and simulation provide insight into mechanisms of granular material flow that are difficult or impossible to study experimentally. High-quality numerical simulations of granular material flow are of great industrial interest, and such simulations require an adequate constitutive model and numerical method, but also high-accuracy experimental data.
Typically, granular material flow is modelled either at the particle scale, or at the continuum scale. The discrete element method (DEM), originally formulated by Cundall and Strack [15], is a method that has been widely used to model granular material flow in various industrial processes. In the DEM, each particle in the granular material mass is represented using a discrete particle. The motion of the discrete particles is determined by Newton's second law of motion, and the motion of the granular material mass is governed by the motion and interactions between the individual discrete particles. In the DEM, a small overlap is allowed at the contact between particles. The overlap is related to contact forces via a force-displacement law. The time integration of Newton's equations typically proceeds in an explicit manner. Explicit integration requires that the time step size is kept low to ensure numerical stability. However, in some recent studies [66,71,72], the DEM has been implemented using implicit time integration schemes, allowing for much larger time steps. Implementations of the DEM commonly use spherical particles to represent non-spherical real granular particles. To accurately represent a non-spherical particle with a spherical particle requires careful selection of contact parameters. Attempts has also been made to use a multisphere approach to model non-spherical particles with the DEM [7]. A multi-sphere approach can be appropriate, but the authors still emphasize the difficulty and limitation of the method: the difficulty in calibrating micro-scale parameters for the DEM models. The DEM has traditionally been applied to model mainly frictional granular materials, but recently a growing interest in wet granular materials has stimulated the development of DEM models including capillary forces between particles [36]. When using the DEM, the computational cost increases with an increasing number of particles. Today, the availability of increasingly powerful computational resources has enabled simulation of large systems of granular materials, containing millions of particles [27]. However, the DEM is still impractical for the simulation of industrial size-scaled granular material flows, typically involving several billions of particles.
In a continuum approach, the granular mass is modelled as a continuum and its behaviour is predicted by fundamental laws of physics, namely the conservation of mass, momentum and energy. By this, the modelling of individual particles is avoided. Traditionally, when the continuum approach has been used for modelling of solid-like granular material behaviour, strain-rate-independent plasticity models, originating from Mohr-Coulomb plasticity, have been used, see e.g. [2,17,54,64]. In the literature [56], it has been shown that the mechanical behaviour of granular materials is strain rate independent in the quasi-static regime and strain rate dependent in the flow regime. In Andrade et al. [3], a strain-rate-dependent constitutive model for granular materials was formulated. The strain rate dependency was included by postulating a material strength that evolves with the strain rate.
When modelling fluid-like granular material behaviour, the focus has mainly been on the prediction of the steadystate flow regime. A visco-plastic rheology model, based on the use of a dimensionless inertia number, was used in [40,65] to model dense granular material flow. Promising results were obtained for flows on inclined planes and for the flow when a granular material was poured on top of a pile. However, the visco-plastic rheology approach is limited to flows where the inertial number is low, corresponding to relatively slow granular material flows. Furthermore, the solid-like granular material behaviour at quasi-static load-ing cannot be captured, and no hysteresis is included. No inclusion of hysteresis means that an event where part of the granular material is static, while some other part of it is flowing, cannot be predicted correctly.
Modelling and simulation of the behaviour of granular materials with a continuum approach require, besides the selection of an adequate constitutive model, the choice of a robust and efficient numerical method. The finite element method (FEM) is a numerical method with a long tradition that has been used for numerical modelling in a wide range of technical fields [88]. When used with a Lagrangian description of motion, large deformations tend to severely distort the FE mesh, resulting in numerical difficulties. Thus, to use the FEM to model large deformation problems, some remedy for the mesh distortion is required. The FEM used with an Eulerian description of motion has been used to model granular material flow because it avoids mesh distortions at large deformation [23,42,83]. However, the FEM with an Eulerian description suffers from difficulties in predicting free surfaces and moving boundaries. The arbitrary Lagrangian-Eulerian (ALE) method attempts to overcome the inherent drawbacks of both Lagrangian and Eulerian descriptions. Drawing on the advantages of pure Lagrangian and pure Eulerian descriptions, the ALE method was used to model granular material flow in [13,14,78]. Advantages and disadvantages of using the FEM for the numerical simulation of forming processes involving large strains are discussed thoroughly in Rodríguez et al. [68].
There exist a number of particle methods within the continuum approach, and they provide an attractive alternative to the above-mentioned numerical methods, for the modelling of granular materials. Particle methods are commonly classified into two categories, particle methods that use a background mesh and particle methods that do not use a background mesh. One example of the latter is the smoothed particle hydrodynamics (SPH). The SPH was originally developed for the simulation of astrophysical problems [28,53]. The SPH is a Lagrangian mesh-free method, where the computational domain is represented by a set of particles. The particles also serve as the frame over which the field equations are approximated. In the SPH, no direct connectivity between particles exists. Thus, it can be used to treat problems involving large deformation, without suffering from the numerical difficulties inherent in mesh-based methods. The original SPH suffers from a number of drawbacks such as tensile instability, a lack of interpolation consistency, zero-energy modes, difficulties in handling essential boundary conditions and non-physical pressure oscillations. Furthermore, the SPH requires a homogeneous and smooth particle distribution to obtain stable and reliable results. This becomes particularly important in the evaluation of the pressure field. In recent versions of SPH, many of the inherent drawbacks of the original version have been addressed and solved. Today, SPH is used for modelling in a wide range of engineering applications, including simulation of slope stability analysis and failure [6,58], and for the modelling of granular material flow [5,30,[37][38][39]48,57,63].
The material point method (MPM) is a particle method that uses a background mesh. The MPM was developed by Sulsky et al. [75,76], and it is based on a combined Lagrangian-Eulerian description of motion. In the MPM, the state variables are traced on Lagrangian material points, while the equations of motion are integrated on a background computational mesh. The material points can be chosen independently of the mesh, and the connectivity between the material points changes dynamically during the simulation. In the MPM, each particle is assigned a fixed mass, which ensures mass conservation as long as the number of particles is kept constant throughout the simulation. Mass conservation is a strong advantage of the MPM, but since it is not easy to add or remove particles inside a given mass, the concentration of mass in localized areas might create problems. Since its original formulation, the MPM has become widely used for modelling in computational mechanics. Noteworthy publications include [10] where the MPM was used to treat localized large deformations for brittle failure and in [1] to model landslides using a Mohr-Coulomb yield criterion. Furthermore, the literature contains a number of studies where the MPM has been used to model granular material flow, for instance to model the granular material column collapse of frictional materials [20,24,55,74], cohesive-frictional materials [31] and silo discharge problems [80]. In a recent study by [32], a stabilized mixed implicit MPM was developed and used to model incompressible and compressible materials with a variety of plasticity laws. The MPM, combined with appropriate constitutive models, has been shown to adequately predict granular material flow at varying flow conditions.
The particle finite element method (PFEM) is another mesh-based particle method, and it is a Lagrangian particle method based on the FEM. The PFEM was initially developed for solving fluid dynamics problems in the context of fluid-structure interaction and free-surface flow [33,34,60]. The first extension of the PFEM to solid mechanics applications was made by Oliver et al. [59]. Recently, the flexibility and robustness of the PFEM have been demonstrated in a variety of engineering applications such as modelling of water waves generated by landslides [11,70], modelling rapid landslide run-outs with a Drucker-Prager model cast as a Bingham, non-Newtonian fluid model [12] and modelling of failure of rockfill dams with a Mohr-Coulomb model in a non-Newtonian Bingham form [46]. Furthermore, the PFEM has been used to model a variety of granular material flow problems, see e.g. [9,16,45,86]. In the PFEM, a Lagrangian description of motion is used for the nodes in a finite element mesh. The nodes are considered as free particles that are allowed to separate from the domain they originally are a part of. A cloud of particles is used to identify the computational domain, and a finite element discretization is utilized to advance the solution by a time increment. The particles contain all properties and variables, and the values of those are projected onto the mesh at each time increment, where the necessary equations are solved. The PFEM predicts a smooth pressure field, and the non-physical pressure oscillations typical for the SPH method are avoided. One advantage of the PFEM compared to the MPM is that in the PFEM it is possible to add or remove particles during the simulation; thus, concentration of mass in localized areas can be avoided. If a numerical model is to be used for industrial decisionmaking, it needs to produce trustworthy results. Here, validation against experimental observations is of major importance. The continued development of numerical methods for modelling granular materials requires improved experimental methods that can provide high-accuracy results for model calibration and validation. Experimental insight into the flow dynamics can be obtained through in-plane velocity measurement, and such measurements are very useful for validation of numerical models. The introduction of digital photography has led to the development of optical experimental techniques such as digital particle image velocimetry (DPIV) [81]. In DPIV, a cross-correlation method is applied to a series of digital images to obtain the in-plane velocity field. DPIV was applied by [73] for field measurements of granular material flow during the discharge from plane hoppers. Digital image correlation (DIC) is an optical experimental technique that has been used extensively for the displacement and strain field measurement of materials subjected to large strains [41,62]. The DIC technique is based on the comparison of a series of digital photographs of a specimen surface recorded during deformation. Similar to the DPIV technique, a crosscorrelation procedure is applied to determine the in-plane displacement field.
The collapse of granular material columns is an experimental set-up that has received much attention in recent years. The simplicity of the set-up and the ability to use it to study complex granular material flow phenomena have made the column collapse popular in the field of granular material flow. The set-up was popularized by [43] and [51] where cylindrical columns producing axisymmetric collapses were used. Their work revealed that the flow dynamics and deposit morphology were primarily dependent on the initial aspect ratio between the height and radius of the columns. Nonintrusive measurements of the deposit morphology using a laser scanner equipment were presented in [79]. In [85], particle tracking velocimetry was used to obtain in-plane velocity measurements of collapsing granular columns. In [47], the DIC technique was used to study and characterize granular material flow through field measurements.
The aim of the present study was to obtain a robust numerical tool for the simulation of granular material flow at dissimilar flow conditions. Transient granular material flow was modelled by the PFEM and a novel constitutive model where the internal friction of the granular material evolves with strain rate. To assess the capability of the proposed modelling framework, a set of numerical simulations of the fundamental problem of collapsing rectangular columns of granular material were performed. An extensive experimental investigation of the column collapse problem was performed where a novel experimental technique based on DIC was applied to quantify the flow dynamics in the form of in-plane velocity measurements. The proposed numerical model was validated by comparing results from the simulations with experimental measurements.

Materials and experimental study
In the present study, the flow dynamics of two granular materials was investigated experimentally. The main purpose of the experimental study was to obtain qualitative and quantitative measures of the flow of the granular materials, to be used for the calibration and validation of the proposed numerical model.

Materials
In the present study, two granular materials with very different properties were investigated. The first granular material was a potassium chloride (KCl) fertilizer, commonly known as muriate of potash (MOP). MOP denotes mixtures of KCl, at 95% or greater purity, and NaCl, which are adequate for agricultural use [77]. The particle size for the granular MOP was in the range of 2.0-4.0 mm, the particle shape was angular, the particle density was 1.99 g/cm 3 , and the bulk density was 1.00 g/cm 3 . The second granular material was a sintered aluminium oxide (Al 2 O 3 ) used as a grinding media in liquid fine grinding in stirred media mills. The particle shape of the Al 2 O 3 was spherical, and the particle size was in the range 1.2-2.0 mm. The particle density was 3.41 g/cm 3 , and the bulk density was 2.13 g/cm 3 . Optical light microscope images of the two granular materials are shown in Fig. 1, and it is observed that the granular materials have very different particle shapes. The particle shape of a granular material affects the internal angle of friction, which in turn affects its flowability. It has been shown [4,82] that particles with fairly spherical shape have a significantly lower internal angle of friction compared to that of particles with an angular and rough shape. In the present study, two granular materials with different particle shapes were chosen. The Al 2 O 3 consists of spherical particles and the KCl with angular and rough particles. This choice was made deliberately to investigate the ability of the proposed numerical model to represent granular materials with very different properties and bulk flow characteristics. Furthermore, an objective of the present study was to evaluate the performance of the proposed optical experimental technique for granular materials with different properties. Thus, this choice of granular materials was considered adequate by the authors.

Experimental set-up and procedure
The study of the transient material flows that occur during the collapse of columns of granular materials has been the focus in a number of studies. The collapse of axisymmetric granular material columns was studied experimentally in [43,51]. The rectangular channel column collapse was studied in [4,52]. The granular column collapse problem includes the kinematics of granular materials on different stages. Initially, the material is at rest in its container, and it then undergoes acceleration during the collapse and deceleration when the material comes to rest. Thus, the granular material column collapse experiment provides a good foundation for evaluating a numerical model of transient granular material flows. The experiments carried out in the present study constitute An illustrative drawing of the experimental set-up used in the present study is shown in Fig. 2. The experimental setup was designed as a 50-mm-wide, 590-mm-long and 230mm-high closed rectangular channel. The front panel was made of hardened glass with a thickness of 4 mm. The other panels of the channel were made of steel with a thickness of 6 mm. The bottom of the channel was open, allowing it to be placed on surfaces made of different materials, with different surface properties. For all experiments in the present study, a smooth and horizontal bottom surface made of steel was used. A 6-mm-thick steel door was used to confine the granular materials in a reservoir prior to the collapse. The position of the door could be varied, enabling the study of a wide range of initial shapes of the granular mass. The design of the experimental set-up used in the present study was based on the set-ups used for rectangular channel column collapse experiments in [4,44].
The experimental procedure consisted in an initial positioning of the door and thus selecting a length l i of the reservoir. The reservoir was then partly filled by carefully pouring a granular material mass to a height h i . Thus, forming a rectangular column of granular material with a length of l i , a height of h i and a width of 50 mm. The top surface of the granular material was evened out by hand. The door was then quickly removed vertically via a weight, rope and pulley system (Fig. 2). The use of a weight, rope and pulley system made it possible to remove the door in a reproducible manner, keeping the vertical speed of the door constant at approximately 1.5 m/s for all the experiments. When the door was removed, the granular mass collapsed under the influence of gravity and spread horizontally in the channel until it came to a rest, forming a deposit profile. A conceptual initial set-up and final deposit profile are shown in Fig. 3. The ratio between the initial height h i and length l i was used to express the aspect ratio a i = h i /l i of the granular mass. The aspect ratio was varied by either using different amounts of granular material for a fixed door position and thus varying h i , or by re-positioning the door and thus varying l i . This procedure enabled the investigation of the collapse of columns of different masses but with the same aspect ratio. In total, 17 experiments were carried out for each granular material and the experimental parameters are presented in Table 1. Since all possible combinations of experimental parameters would result in a huge test matrix, the choice of h i and l i was made arbitrarily with the aim to cover the span of aspect ratios given in Table 1. The experiments were recorded with a high-speed digital camera. A MATLAB script was used to process the digital images to extract the final height h ∞ , the final length l ∞ and the granular material deposit profile. For the KCl, l ∞ was defined as the horizontal position where the grains remained in contact with the rest of the granular mass. Thus, individual grains that had separated from the mass were not considered. For the Al 2 O 3 , l ∞ was defined as the horizontal position where the granular material layer ceased to have at least two grains in thickness. Furthermore, the digital images were processed using a commercial digital image correlation (DIC) software [29]. The methodology described in [47] was used to obtain the in-plane velocity field. In short, the DIC technique is based on the comparison of a series of digital images that are divided into overlapping sub-images. The in-plane velocity field is determined by applying a cross-correlation algorithm, which requires that the object to be traced is covered with a random surface pattern. The correlation algorithm is then able to trace the motion of the sub-images, and thus, the velocity field can be obtained. The granular materials that were used in the present study form a natural random surface pattern, and with a sufficient surface texture, the DIC technique could be used to obtain the in-plane velocity field. A more detailed description of the DIC and its application to quantify granular material flows can be found in [47].

Data acquisition
To record the experiments, a Redlake MotionPro X3 highspeed digital camera was used. During the recording, the granular materials were illuminated using two Dedocool floodlights equipped with 250 W lamps. The experiments were recorded with the high-speed camera set to capture 1000 images per second, at a resolution of 1280 × 720 pixels and with a shutter speed of 0.25 ms.

Numerical modelling and simulation
A granular material is a discrete media. However, in the present study the modelling of granular materials was based on the assumption that a granular material can be represented as a continuous media.
The assumption of using a continuum representation of discrete media is valid as long as the particles are much smaller than the smallest characteristic dimension of the process considered [21]. In the present study, two-dimensional computational domains were used to represent the granular materials. The PFEM, implemented in a MATLAB program, was used, and the granular materials were modelled using two different constitutive models.

Governing equations
The balance of linear momentum can be expressed in a Lagrangian description as where v i and b i are the velocity and body force components, ρ is the density, x j are the material point positions, σ i j is the Cauchy stress tensor and Dv i Dt is the material derivative of the velocity field The Cauchy stress tensor can be split into a mean stress component σ 0 = 1 3 tr(σ i j ) and a deviatoric component s i j according to where δ i j is the Kronecker delta. Furthermore, it is assumed that the mass of a continuum body is conserved and that it is a continuous function of volume. The conservation of mass can be stated as where κ is the elastic bulk modulus, Dσ 0 Dt is the material derivative of the mean stress and ε V is the volumetric strain rate. The volumetric strain rate is defined as the trace of the rate of deformation tensor d i j , which is given by

Constitutive models
Modelling a granular material as a continuum requires a constitutive model where the stresses in the material are related to some measure of deformation. Constitutive models may be dependent or independent of the strain rate, and in the present study, two strain-rate-dependent constitutive models were evaluated and compared.

Flow formulation
The first constitutive model is based on a constitutive relation for the flow of plastic and visco-plastic solids. It was originally outlined in [87] and was specialized in [8] to a Drucker-Prager yield surface [19], with a non-associated flow rule.
For large deformation, under plastic or visco-plastic conditions, elastic deformations can be neglected. A constitutive model linking the stresses and strain rates, where the viscosity is dependent on the current strain rates, can be formulated using the analogy with a viscous non-Newtonian incompressible fluid. The constitutive relation for an incompressible viscous fluid can be expressed as where μ is the viscosity andε i j is the strain rate tensor. Equation (6) can be rewritten using the split of the Cauchy stress tensor from Eq. (3) and following the definition of Perzyna non-associated viscoplasticity [69], Eq. (7) can be written aṡ where μ p is a 'pseudo-viscosity', F = F(σ i j ) = 0 is a plastic yield surface and G = G(σ i j ) is a plastic potential function. The use of Macaulay brackets in Eq. (8) means that F = F if F > 0 and F = 0 if F ≤ 0, thus ensuring no development of plastic flow if the stress state is inside the yield surface. If the viscosity parameter μ p → 0, it implies that F → 0 in order forε i j to be a finite quantity. Thus, ε i j →λ∂G/∂σ i j , whereλ is the plastic multiplier. In other words, the visco-plastic relation in Eq. (8) reduces to rate independent plasticity theory when μ p → 0. Cante et al. [8] specialized the Perzyna relationship to the Drucker-Prager yield surface (Fig. 4), which has the following functional form The parameter b 1 controls the influence of the mean stress on the yield limit, and it can be interpreted as the internal coefficient of friction of a granular material. The parameter b 2 corresponds to the yield strength of the material under pure shear, and in the context of granular materials, it can be interpreted as the granular material cohesion.
In the present study, the flow rule assigned to the Drucker-Prager yield surface is non-associated and consists of a purely deviatoric strain rate, and it can be expressed as The material flows at F ≥ 0 and from Eqs. (8) to (10) the following expression can be obtained whereε i j = 2 3 ||dev(ε i j )|| is the effective strain rate. For ideal plasticity,μ → 0 and the viscosity can be written as Thus, a relationship between the deviatoric stresses and the strain rate has been obtained. Using the above results, this expression can be written in its final form as

Flow formulation with strain-rate-dependent residual strength
The second constitutive model used in the present study is based on a strain-rate-dependent plasticity model introduced in [3]. Considering the Drucker-Prager yield surface outlined in the previous section, the parameter b 1 is interpreted as the frictional resistance of the granular material. Andrade et al. [3] proposed a frictional resistance that is dependent on the dilatancy β and on a residual resistanceμ The dilatancy is considered to be a function of the deviatoric shear strain ε s , and its evolution is given by where β * is the maximum dilatancy and ε * s is the corresponding deviatoric shear strain. Following the form outlined in [40], the evolution of the residual resistance is a function of the deviatoric shear strain rateε s , and it is given bȳ whereμ l andμ u are the lower and upper bounds for the residual resistance, respectively. The lower and upper bounds are represented byε s → 0 andε s → ∞, respectively. The parameterε * s is a characteristic deviatoric shear strain rate at which the residual resistance isμ = 1/2(μ l +μ u ). Thus, the frictional resistance b 1 is dependent on both the deviatoric shear strain and the deviatoric shear strain rate. The evolution of the frictional resistance is given by The role of the dilatancy is to couple the deviatoric and volumetric components of deformation, and it describes the volume change of a material under shear deformation. The dilatancy is important for the mechanical behaviour of granular materials at quasi-static loading. It is in contrast to other materials, such as metals which are non-dilative. Dilatancy is important for granular materials in the solid-like state, but it can be neglected in the fluid-like state. In most granular material flows, the variation of the volumetric fraction is small [25], and if the granular material is considered as incompressible, the dilatancy and frictional equations are decoupled. Thus, the incompressible assumption greatly simplifies the constitutive model. In the present study, the granular materials were modelled as quasi-incompressible. Thus, the effect of the dilatancy was not included and the evolution of the frictional resistance is given by The conceptual evolution of the residual strength as a function ofε s andε * s is shown in Fig. 5.

The particle finite element method
The PFEM is a particle-based numerical method where a background mesh is used and on which the FEM is used to solve the governing equations. The PFEM is founded upon modelling using an updated Lagrangian formulation. In the updated Lagrangian formulation, the equations are formulated in the current configuration, and the variables are assumed to be known at the last calculated configuration, at time t. The new variables are sought at the updated configuration, at time t + t. As outlined in [67], the PFEM can be divided into the following basic steps: A finite element mesh is generated, using the set of particles as nodes. The finite element mesh is obtained using a Delaunay triangulation [49]. 3. The alpha-shape method [22,84] is used to identify the external boundaries onto which the boundary conditions are imposed. 4. The nonlinear governing equations are solved for displacement, velocity and pressure at every node of the mesh. 5. Computed velocities and pressures are used to update the position of the particles. 6. Return to step 2 and repeat for the next time increment.
Thus, the PFEM can be interpreted as an updated Lagrangian approach, where the FEM is used to solve the incremental problem. In the PFEM, the mesh works as the background mesh for integration of the differential equations, and simultaneously, the mesh is used to keep track of free surfaces and contacts. Similar to the standard FEM, the accuracy of the solution in the PFEM depends on the mesh density and quality.
In a Lagrangian description of motion, the particles in the finite element mesh also represent material particles. Thus, the particles will move with the flow of the material. The motion of the particles might result in regions of increased concentration of particles and consequently regions where the particle spacing is large. The accuracy of the solution is affected if the distribution of particles becomes too irregular. In the present implementation of the PFEM, this issue is addressed by allowing the removal and addition of particles. A geometric criterion based on a characteristic element size and distance between particles governs the addition and removal of particles. In the PFEM, contact between the deforming material domain and fixed boundaries is detected Fig. 6 Conceptual illustration of the particle discretization of the domain in the two-dimensional plane deformation column collapse simulation. The initial particle disposition is regular and rectangular. In the figure, the fixed particles used to represent the bottom surface and left wall are shown in a darker shade automatically during the mesh generation, and no contact search algorithm is required. Penetration of the nodes of the deforming material into the fixed boundaries is prevented by the incompressibility condition. In the present implementation of the PFEM, frictional contact between the deforming domain and the fixed boundaries is modelled via the frictional resistance of the deforming material. More details regarding the automatic contact treatment of the PFEM can be found in [61].

Simulation procedures
The column collapse experiment was simulated using a two-dimensional plane deformation implementation of the PFEM. The granular material mass was represented using particles which were initially arranged in a regular rectangular pattern. The bottom surface and the left wall were modelled as stationary particles. The initial particle disposition and the location of the fixed particles used as boundaries are conceptually illustrated in Fig. 6. Throughout the present study, a stabilized linear triangular mixed velocity-pressure finite element formulation was used to solve the Lagrangian equations [18]. A fully implicit scheme was used for the time integration where the time step size was a function of the maximum velocity and the minimum distance between the particles. A maximum allowed value of the time step was set to t = 1.0 × 10 −4 s, and a convergence criterion of 10 −4 was used.
In the literature [44,52], the time evolution of the flow front in the column collapse is commonly described using a characteristic time scale based on the free-fall time of the granular column τ c = √ h i /g. In the present study, the simulations were terminated at the normalized timet = t/τ c = 4.0, at which the flow front propagation was assumed to have ceased for the investigated range of aspect ratios. In a comprehensive experimental study of the collapse of granular columns along a horizontal channel, Lube et al. [52] derived at = 3.3. Thus, the assumption of a ceased flow front att = 4.0 is considered adequate and conservative. Since the granular materials in the present study were assumed to be dry and cohesionless, the

Results and discussion
In the following section, the experimental and numerical results are presented and discussed. The flow dynamics of the rectangular column collapse was investigated for a range of initial aspect ratios. The PFEM was used with two different constitutive models to simulate the experiments, using a twodimensional plane deformation formulation. The numerical mesh convergence was studied, and the constitutive models were calibrated by inverse modelling. The constitutive

Experimental observations
The flow dynamics of the granular column collapse was studied experimentally for the two granular materials. A series of representative examples showing how the flow dynamics varied with the aspect ratio a i are shown in Figs. 7 and 8. When the value of a i was low (Fig. 7), the flow was contained in the top surface layer, and most of the granular mass was stationary during the collapse. The final deposit profile at low values of a i had a characteristic truncated cone shape, where a large part of the granular mass remained undisturbed during the collapse. Increasing the value of a i resulted in a larger proportion of the granular mass being disturbed during the collapse, and the final deposit profile converged to become increasingly more cone shaped.
It is noted from Figs. 7 and 8 that the two granular materials resulted in final deposit profiles of different shapes, for similar values of a i . The KCl resulted in a lower value of l ∞ compared to that of the Al 2 O 3 , for both the low and the high a i . Both granular materials resulted in similar values of h ∞ for the low a i experiment, while h ∞ was larger for the KCl than that of the Al 2 O 3 for high a i . Thus, the properties of the granular materials had an effect on the flow dynamics and on the shape of the deposit profiles. The angular shape of the grains of the KCl yields a higher internal friction compared to the spherical grains of the Al 2 O 3 . Thus, the Al 2 O 3 flows more easily than the KCl. These observations are in line with the results of previous studies [4,43], where it was reported that the flow and spreading dynamics were mainly dependent on the value of a i , but also on the internal friction of the granular materials.
It was observed that repeated experiments with the same input parameters resulted in slightly different final heights and lengths. This was expected, at least to some extent, since the nature of the granular material is very random and it is difficult to obtain exactly the same material distribution in the reservoir before each experiment. This observed randomness in the experiments is probably dominating over possible instrumentation errors.

Mesh convergence study
The mesh convergence of the current implementation of the PFEM was studied by running a test problem and varying the initial distance between the particles. For this, the rectangular granular column collapse problem with an initial height h i = 240 mm and initial length l i = 24 mm was chosen. The selected dimensions corresponded to a i = 10. A high aspect ratio was selected to produce a case with large deformations and large pressure and velocity gradients. Three initial particle distances were selected: 1.5 mm, 0.75 mm and 0.375 mm, resulting in models containing approximately 3500, 12100 and 44800 particles, respectively. For the mesh convergence study, the flow formulation constitutive model was used with b 1 = 1.0 and b 2 = 10 −6 Pa. The bulk density was set to 1.00 g/cm 3 . To evaluate the mesh size dependency, the final deposit profile was extracted and compared for the three models. In Fig. 9, the final deposit profiles are compared, showing that they are barely distinguishable. Thus, an initial particle distance of 1.5 mm was considered adequate and was used throughout the present study.

Constitutive model parameter calibration
The calibration of the constitutive model parameters was conducted in two steps, using the rectangular column collapse experiments. In the first step, the parameter b 1 was calibrated using a low aspect ratio experiment, with a i = 0.73 for the KCl and a i = 0.75 for the Al 2 O 3 . Experimental and simulated final deposit profiles were compared for different values of b 1 , as shown in Fig. 10. In addition to comparing the shape of the deposit profile, the experimentally and numerically obtained values of h ∞ and l ∞ were also compared. From Fig. 10, it is observed that lowering the value of b 1 resulted in a less viscous behaviour of the granular mass, which is   Increasing the value of a i resulted in increasing strain rates during the column collapse, as shown in Fig. 11. The value of b 1 obtained through the initial calibration failed to predict h ∞ and l ∞ for the experiments with larger a i , for both the KCl and for the Al 2 O 3 . For the KCl, h ∞ was overpredicted, while l ∞ was underpredicted for increasing values of a i , indicating that a b 1 = 1.13 overpredicted the material strength at increasing strain rates. For the Al 2 O 3 , it was the other way around, h ∞ was underpredicted, while l ∞ was overpredicted for increasing values of a i , indicating that a b 1 = 0.83 underpredicted the material strength at increasing strain rates.
Thus, to accurately model the rectangular column collapse at increasing values of a i required a constitutive model able to account for the strain-rate-dependent material strength. The second constitutive model of the present study allows the internal coefficient of friction (the parameter b 1 ) of the material to be dependent on the strain rate. This model requires the definition of some additional parameters governing the evolution of the residual resistanceμ, which is equal to b 1 for a granular material in which the effect of dilatancy can be neglected. As presented in Sect. 3.2.2, these parameters are the lower and upper bounds ofμ (μ l andμ u ) and the equivalent deviatoric shear strain rateε * s . In the literature [3,40,50], the identification of model parameters corresponding toμ l andμ u has been discussed. The proposed relationships between the quasi-static material strength and the material strength at high strain rates vary depending on the granular material and the type of flow considered. It is thus difficult to know this relationship a priori for arbitrary granular materials and flow conditions. In the present study, the choice of relationship betweenμ l andμ u was based on experimental observations from the column collapses. For the KCl, the value ofμ l was set to be equal to the previously obtained value of b 1 = 1.13, while the value ofμ u was set to 0.9 × b 1 . Thus, the residual resistance was set to decrease with increasing strain rate. For the Al 2 O 3 , the value ofμ l was set equal to b 1 = 0.83 and the value ofμ u was set to 1.1 × b 1 . Thus,μ was set to increase with increasing strain rate.
The remaining model parameterε * s governs the transition betweenμ l andμ u with increasing strain rate shown conceptually in Fig. 5. The value ofε * s was obtained using three column collapse experiments, where a i was varied between 0.73 and 5.13 for the KCl and between 0.75 and 4.92 for the Al 2 O 3 . The same procedure that was used to calibrate b 1 was s . The best match to the experimentally obtained h ∞ and l ∞ was obtained withε * s = 25 for the KCl andε * s = 10 for the Al 2 O 3 . In Fig. 13, the shape of the deposits obtained experimentally and with the two constitutive models of the present work is shown for the KCl and the Al 2 O 3 . The evolution of the parameter b 1 as a function of the strain rate is shown in Fig. 12 for both granular materials.
The discrepancy between the experimental and the simulated deposit profiles in Fig. 13 is slight for the Al 2 O 3 , but more prominent for the KCl. A possible explanation is that the assumption of using a continuum representation of a discrete media might be questionable in the case of KCl. To accurately model a discrete media as a continuum requires that the particles are much smaller than the smallest characteristic dimension of the process considered [21]. The size and the angular shape of the grains of the KCl, and the length scale of the experiments of the present study cause the flowing layer of grains to be thin in comparison with the grain size in some of the experiments, thus making the continuum representation questionable at that location. A further possible explanation for the discrepancy is that the shape of the grains of the KCl might result in some dilatation during the column collapse, thus resulting in a slight volume increase. The dilatation is typically small for a granular material in the fluid-like state [25], and it is typically neglected. In the present study, the granular materials are modelled as incompressible; thus, any dilatation is not taken into account which is an additional contributing factor to the slight deviation between simulated and experimental profiles for the KCl.
It must also be noted that when using the PFEM, the remeshing method may cause a slight variation of the volume of the computational domain. To prevent this, the choice of the value of the alpha-shape parameter should be carefully considered. The issue of volume conservation and remeshing in the PFEM is discussed in detail in [26] where it is suggested that values of the alpha-shape parameter close to 1.2 keep the volume variation at acceptable levels for problems of highly unsteady flows. Thus, throughout the present study, the value of the alpha-shape parameter was set to 1.2.
In the present study, a strategy for calibration of the constitutive parameters based on a comparison between experimental and simulated column collapses is presented. Using this approach, a set of constitutive parameters were obtained for the two granular materials. To obtain a unique set of constitutive parameters is indeed a difficult task, partly because it is difficult to experimentally measure the evolution of the frictional resistance as a function of the strain rate in a granular material. Thus, the parameters obtained in this study are considered to be of use in the present application and at the investigated loading conditions. A strategy to obtain a unique set of constitutive parameters for arbitrary loading conditions and granular materials would indeed be an improvement but lies outside the scope of the present study.

Model validation
To validate the proposed numerical model, a number of column collapses with varying a i were simulated and compared to experimental results. In total, the 17 cases from Table 1 were simulated for each granular material, where a i was varied between 0.73 and 9.60 for the KCl and between 0.60 and 9.70 for the Al 2 O 3 . Experimental and simulated values of l ∞ and h ∞ are presented in Tables 2 and 3. The experimentally measured h ∞ and l ∞ were normalized with respect to the initial length l i and were plotted as a function of a i , using a logarithmic scale on both the horizontal and the vertical axes, as shown in Figs. 14 and 15. The normalized h ∞ and l ∞ obtained from the simulations were plotted together with the experimental results, as shown in Figs. 14 and 15. When comparing the experimental and numerical h ∞ and l ∞ , it is observed that the strain-rate-dependent residual strength constitutive model is able to accurately predict the column collapse at the investigated range of a i , for both the KCl and the Al 2 O 3 . An error percentage for simulated l ∞ and h ∞ compared to experimental results was calculated and the median of the error percentage was determined, excluding the tests used for calibration ofε * s , as given in Tables 2 and  3. For the KCl, comparing the median percentage errors, the most accurate prediction was obtained for the constitutive model without strain-rate-dependent residual strength. For the Al 2 O 3 , the strain-rate-dependent residual strength model resulted in a more accurate prediction of l ∞ , while for h ∞ a similar accuracy was obtained for both constitutive models. The best model predictions resulted in median percentage errors less than 5, for both granular materials.
To further compare the experimentally observed flow dynamics to the PFEM simulations, the horizontal and vertical velocity fields were extracted from the column collapse experiment using DIC, as described in Sect. 2.2. In Figs. 16 and 17, the horizontal and vertical velocity fields are compared for the Al 2 O 3 with a i = 0.60 and l i = 80 mm and for the KCl with a i = 0.83 and l i = 72 mm. Compared to the experimentally measured velocities, the proposed strainrate-dependent residual strength model was able to predict the flow dynamics of the column collapse accurately. Furthermore, the time evolution of the column height during the collapse was measured experimentally from the high-speed recording of test no. 19, for both the KCl and the Al 2 O 3 . The test was simulated, and the time evolution is compared in Fig. 18. The PFEM model resulted in a slight underprediction of the time it takes for the column to settle, with a more accurate prediction for the Al 2 O 3 compared to that of the KCl. A possible reason for the underprediction is that the vertical removal of the door is not included in the simulation. The finite time required to remove the door in the experiments might affect the flow dynamics, causing a lower vertical velocity of the top layer of the column during the collapse, compared to a collapse where the door is not included.
In the present implementation of the PFEM, friction between granular materials and surrounding structures is not treated explicitly. This is due to the use of the incompressibility condition to model the interaction between the deforming domain and the fixed boundaries. Thus, the granular material strength governs the flow at the interface between granular mass and fixed boundaries. The use of a simplified contact treatment is given some validity from an experimental study by Lube et al. [51], where column collapses of a number of different granular materials were conducted on three different surfaces: a smooth wooden surface, a smooth plastic surface and a rough surface made of sand paper. The authors found that the shape of the deposits was not significantly affected by the surface properties. It was suggested that a possible explanation for the independence of surface friction was the development of a dynamic interface a few particles from the base surface, separating the flow between stationary and moving granular material.
It should be noted that regardless of the choice of numerical model for the transient granular material flow, validation has to be performed to ensure reliable model predictions. The accurate DIC measurements of the granular material flow dynamics provided a foundation to assess the validity of the proposed PFEM model. In this work, a two-dimensional PFEM model was applied and validated for the fundamental problem of collapsing rectangular columns.

Conclusions
The particle finite element method (PFEM) is used to model the transient granular material flow of a collapsing rectangular column. A novel experimental methodology for quantification of the flow dynamics of the collapsing column of granular material is designed. The experimental results are used to calibrate and validate the proposed numerical model. A conclusion from the present study is that the flow dynamics of the column collapse can be quantified by measuring the in-plane velocity field using digital image correlation. A numerical model, where the PFEM is used with two strain-rate-dependent constitutive models, is evaluated and compared to experimental results. It is concluded that the PFEM model of the present study accurately predicts the flow dynamics of the column collapse for two granular materials with different material properties and over a range of aspect ratios. By validation, it is shown that the strain-ratedependent residual strength constitutive model is the most accurate for the Al 2 O 3 . In general, the best model prediction is obtained for the Al 2 O 3 , while some discrepancy between experimental and simulated results is observed for the KCl. One possible cause of the discrepancy is that the length scale of the granular material flow is too small for the KCl to accurately model it as a continuum. The proposed novel strain-rate-dependent residual strength constitutive model requires the calibration of only three parameters, the lower and upper bounds of the residual resistanceμ l andμ u and the equivalent deviatoric shear strain rate parameterε * s . The number of parameters of the proposed model is low compared to other numerical methods commonly used for the simulation of granular material flow, such as the DEM. In conclusion, the proposed PFEM model is a robust numerical tool that is useful for modelling transient granular material flow.