MPM simulations of granular column collapse with a new constitutive model for the solid-fluid transition

MPM simulations of granular column collapse with a new constitutive model for the solid-fluid transition V International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2017 P. Wriggers, M. Bischoff, E. Oñate, D.R.J. Owen, & T. Zohdi (Eds) MPM SIMULATIONS OF GRANULAR COLUMN COLLAPSE WITH A NEW CONSTITUTIVE MODEL FOR THE SOLID-FLUID TRANSITION IRENE REDAELLI1, FRANCESCA CECCATO2, CLAUDIO G. DI PRISCO3 AND PAOLO SIMONINI4 ¹ Politecnico di Milano (DICA), Piazza Leonardo da Vinci 32, Milano 20133, Italy irene.redaell[email protected] ² Università degli Studi di Padova (DICEA), Via Ognissanti 39, Padova 23129 francesca.c[email protected]nipd.it ³ Politecnico di Milano (DICA), Piazza Leonardo da Vinci 32, Milano 20133, Italy claudio.dipris[email protected] 4 Università degli Studi di Padova (DICEA), Via Ognissanti 39, Padova 23129 paolo.sim[email protected] Key words: dense granular flows, MPM, column collapse, constitutive modelling. Abstract. Dry granular flows are very common both in many natural phenomena (flow-like landslides) as well as in industrial processes. To model granular flows in the framework of continuum mechanics, a key issue is the formulation of a suitable constitutive model, capable of capturing the complex rheological behaviour of the granular material in a wide range of strain rates. In this contribution, a physically based model capable of reproducing the behaviour of granular materials in different flow regimes is briefly summarized. Both the stress tensor and the dissipated energy are calculated as the sum of a quasi-static and a collisional contribution: the former one is modelled by employing an elasto-plastic model incorporating the critical state concept, whereas the latter stems from the kinetic theory of granular gases. The constitutive model has been implemented in the MPM research code Anura3D and applied to the simulation of the collapse of a dry granular column. The results show that the model can deal with the phase transition between the solid- and the fluid-like behaviour. The evolution of the state variables, the state of stress and the energy contributions is investigated in details. 539 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 2 1 INTRODUCTION Flows of dry granular materials characterizes many natural and industrial processes. An important feature of granular materials is that, depending on the strain rates and the grain packing, they may behave like a solid, a liquid or a gas. Two different particle interaction mechanisms can be identified in granular materials: (i) enduring frictional contacts among grains involved in force chains and (ii) inelastic collisions. At high particle concentrations and slow strain rates (quasi-static regime), the first mechanism prevails and the material behaves like a solid, where force chains span the entire domain. At low concentrations and very high strain rates (collisional regime), the second mechanism prevails and the material behaves like a granular gas. In between these two conditions, both mechanisms are relevant, the material is in the “transition regime” and behaves like a granular liquid. The Discrete Element Methods (DEM) is commonly employed to study the flow of granular materials because it is capable of automatically take into account large displacements as well as the granular nature of the material: however, this approach becomes inefficient for high numbers of grains. Large-scale phenomena are more often studied with continuum methods that apply a macroscopic approach. A key issue in this problem is the definition of a constitutive relationship suitable for reproducing the transition form a solid- to fluid-like behaviour of the material. In the literature, quasi-static and collisional conditions are often tackled by employing different approaches. The constitutive models proposed within the soil mechanics community, consider only the quasi-static regime disregarding the effect of collisions among grains. The collisional regime has been largely studied in the context of kinetic theories of granular gases, where the granular temperature is introduced as a state variable measuring the degree of agitation of the system. To reproduce the intermediate regime, several constitutive models have been proposed in the literature [2,1]; however many of them considers only steady state conditions. This paper investigates the potentialities of a recently proposed constitutive model [3] in the simulation of dry granular column collapse with the Material Point Method (MPM). The constitutive model, summarized in Section 2, is based on kinetic theory and critical state elasto-plasticity; it can capture the mechanical response of a granular material under a wide range of strain rates. The unique state variables of this constitutive model are the granular temperature T and the void ratio e, which govern the solid-to-fluid phase transition. MPM is a point-based method specifically developed for large deformations of history dependent materials. It simulates large displacements by means of Lagrangian points moving through an Eulerian grid [4]. A brief introduction of the method is provided in Section 3, followed by a description of the numerical model used for the presented simulations. Section 4 discusses the results of the column collapse problem; in particular, we focus on the energy dissipation mechanisms and the phase transitions. 2 THE CONTINUUM APPROACH FOR GRANULAR FLOWS The constitutive model considers an assembly of identical spherical particles of diameter d and density  p. For the sake of simplicity, inherent and/or induced evolving anisotropy are disregarded. This discontinuous medium can be described, in the framework of continuum mechanics, by writing the field equations, derived from the conservation laws of mass, 540 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 3 momentum and energy as it follows: Mass balance: 0 1 1 u Dt De e (1) Momentum balance: σf u   eDt D e ss 11  (2) Energy balance: Γ Dt DE e qεσ  : 1 1 (3) where D/Dt denotes the material time derivative, e = void ratio, u = vector of the macroscopic mean velocity, f = vector of the external forces per unit volume,  = stress tensor, E = specific internal stored energy, q = vector of the energy flux per unit area and unit time, = strain rate tensor and  = dissipated energy. In this paper, for the sake of simplicity, the energy flux q is disregarded. A parallel scheme is assumed, which implies that E is calculated as the sum of the specific elastic stored energy Eq and the kinetic fluctuating energy Ec (4)  results by the sum of the energy dissipated by force chains  q and the energy dissipated by inelastic collisions  c (5)  is given by the sum of a quasi-static contribution c, associated with long elapsing frictional contacts among grains involved in force chains, and a collisional contribution q, associated with inelastic collisions among grains. (6) 2.1 Quasi static contribution In this paragraph the quasi-static contributions, Eq,  q and q are derived on the basis of standard perfect elasto-plasticity incorporating the critical state concept The main ingredients are  Strain rate additivity: ,where e and p stand for the elastic reversible and the plastic irreversible contribution, respectively  The specific hyperelastic stored energy (Gibbs energy function) [5]: )1(3 )2)(1( 1 2 0 nk p nnkp p Eq n a n q     (7) where 5.0 2 0 2 :)1( 9         g nkp p qqq ss (8) pq = mean quasi-static stress, sq = deviatoric stress tensor, k and g are two macroscopic material dimensionless constitutive parameters associated with the bulk and shear stiffness respectively, n = non-dimensional constant and pa = 541 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 4 atmospheric pressure.  A yield locus f(q)=0, of the Lade and Duncan type, together with the consistency rule . f has been defined on the basis of DEM simulations [6–8] and its size depends on the triaxial compression stress ratio Mc.  The critical state locus equation G1(pq,e)=0, function pq and e, defined on the basis of DEM simulations [9]. G1 is affected by the particle Young modulus Ep, by the critical void ratio ec (representing the critical state void ratio for pq=0) and by the material parameter aL.  The plastic potential g(q,e)=0, evolving with e and q, in order to incorporate the critical state into the model [10,11]. g is affected by Ep and by  , which is a macroparameter influencing the rapidity in achieving the critical state locus [3]. As a consequence, the flow rule, necessary to compute the plastic strain rate p ε  , is nonassociated.  An extended consistency rule, , which imposes that the critical state locus cannot be abandoned once is reached (critical regime)[3]. In fact, the basic idea of the model is that 0 represents an attractor locus for the material. Therefore, condition , implies that under the critical regime, the void ratio evolution is governed by the mean quasi static stress via the critical state locus definition. In this regime, the plastic volumetric strain rate does not obey the flow rule but is computed by depurating the total volumetric strain rate of the elastic component. For the sake of brevity, the analytical expression of f, g and G1 are here omitted, but can be found in [3,16]. According to the previous assumptions it is now possible calculate the plastic dissipated energy as p qq Eεσ  : (9) and the quasi static stress tensor p qq εDσ  (10) where Dq(q,e) is the quasi-static stiffness fourth order tensor, depending on the quasi static stress tensor q and the void ratio e. The details of the derivation of Dq are reported in [3,16]. 2.2 Collisional contribution The collisional contributions Ec,  c and c are function of both the void ratio e and the granular temperature T. In particular, according to kinetic theories of granular gases [12], the fluctuating kinetic energy is T e E sc  1 1 2 3  (11) The expression of the energy dissipated by inelastic collisions  c is based on the extended 542 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 5 kinetic theory [13], modified by Berzi and Jenkins [14], to take into account the role of particle stiffness on collisions. Its expression is reported in [14] and it is affected by p, d, the normal and effective coefficient of restitution [15] n and eff, respectively. The collisional stress tensor is modelled according to the kinetic theory described in [12], taking also into account the modification introduced by Berzi and Jenkins [14]. The analytical expression of c is given by vv c hεΦσ   : (12) In Equation 12 v(e,T) and hv(e,T) are a fourth-order and a second-order tensor, respectively. They are affected by s, d, Es, n and a material constant em and their analytical form can be found in [16]. 2.3 The visco-perfect-elasto-plastic constitutive relationship According to Equation 6, the stress tensor is the sum of the quasi static and collisional contributions. To this purpose, the quasi-static contribution is expressed in finite form by integrating in time Equation 10 qqq tσεDσˆ :  (13) where  t is the integration time step and is the quasi-static stress tensor calculated at the previous time step. From Equation 12 and Equation 13, the total stress tensor reads cεDσ  : vep (14) where Dvep (q,e,T) is the visco-elasto-plastic fourth-order tensor vqvep tΦDD  (15) and c( ,e,T) is a second-order tensor given by v q hσc ˆ (16) According to this model, different flow regimes can be defined. 1. The visco-elastic regime, when f<0 or f=0 and =0 2. The visco-ealsto-plastic regime, when f= =0 3. The critical regime: when f= =0 4. The collisional regime: when f= =0 The description of these regimes and the analytical expression of Dvep and c for each regime is reported in [3]. 3 THE NUMERICAL MODEL MPM has been developed to overcome the difficulties arising when simulating large deformation problems with Lagrangian FEM. In the last 10 years the number of applications, ranging from solid to fluid mechanics, has been significantly increasing. In the field of geomechanics it has been used to study a number of problems such as slope stability [18,17], 543 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 6 granular flow propagation [19,20], soil penetration problems [22,23,21], soil erosion and sedimentation [24,25]. MPM uses two level of discretization. On one hand, the continuum body is discretized by means of a set of material points (MP) that carry all the information of the continuum (density, velocity, stresses, state parameters etc.). On the other hand, the spatial region where the body move through is discretized by a finite element mesh that allows to solve the discretized governing equation of motion. Large deformations are simulated by MP moving through a fix mesh as shortly summarized in Figure 1. The constitutive model is implemented as a user defined model in the MPM dynamic explicit research code Anura3D (www.anura3d.eu)[26]. The constitutive equations are integrated explicitly; in particular, a Dormand-Prince method with substepping and correction for yield surface drift is applied for the quasi-static stresses [27]. The specific collisional dissipated energy and the plastic dissipated energy are integrated by means of the mid-point rule. The kinetic fluctuation energy and the elastic stored energy are computed with equations (11) and (7), respectively, at the end of the time step as function of the updated stress state and granular temperature. The simulated 2D column collapse problem considers a 0.1m-high, 0.1m-long column (aspect ratio a=l/h=1). The model is 0.005m-wide. The mesh is structured and counts 5610 tetrahedral elements (Figure 2). The bottom of the mesh is fully fixed, while roller boundary conditions are applied on the other surfaces. 10 MP are initially placed inside each active element for a total number of 14400 MP. Stresses are initialized with horizontal coefficient at rest K0=0.5. The material parameters are summarized in Table 1. Figure 1 Computation scheme of MPM: a) interpolate state variables to the grid nodes, (b) solve the governing equations of motions at the nodes, (c) update MP velocity, strain, stress etc., (d) update MP housekeeping. 544 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 7 Figure 2 Geometry and discretization of the numerical model Table 1 Material properties and initial state variables Grain diameter d [m] 0.002 Grain density  s [kg/m 3 ] 2532 Grain Young modulus E p [Pa] 3e+08 Normal restitution coefficient  n 0.879 Parameter of collisional functions e m 1.5 Critical void ratio e c 0.709 Parameter of CSL a L 0.56 elastic constant related to bulk modulus K 380 elastic constant related to shear modulus G 48 Parameter of plastic potential  0.35 Stress ratio under triaxial compression Mc 0.9 Effective restitution coefficient  eff 0.59 Initial void ratio e 0 0.71 Initial granular temperature T 0 [m2/s2] 1e-15 4 RESULTS The collapse of a granular column is a well-established experiment that consists in releasing a column of granular material by removing its lateral support on to a flat surface. The column then fails and some of its mass flows on to the flat surface before it is deposited. Figure 3 shows the deviatoric strains at the material points (MP) at different time instants, from which the development of the failure surface can be noted. The failure surface defines the boundary between a static cone with small strains (bottom left, blue MP) and a mobilised mass (top right, red MP). During the collapse, the mobilised mass slides along the failure surface and crumbles upon the base. The friction between the mobilised mass and the static region or the boundary, as well as the collisions between grains, dissipate energy and slows down the mobilised mass until static conditions are reached. Figure 4 and Figure 5 show the specific collisional dissipated energy and the specific plastic dissipated energy, respectively. It can be noted that dissipation due to inelastic collisions among particles mainly occurs inside the mobilized mass, while the plastic dissipation mainly occurs at the sliding surface. 545 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 8 Figure 3 Deviatoric strain at the MP in different time instants Figure 4 Specific collisional dissipated energy at the MP in different time instants Figure 5 Specific plastic dissipated energy at the MP in different time instants 546 I. Redaelli, F. Ceccato, C.G. di Prisco and P. Simonini 9 Integrating the specific energies at the MP, it is possible to obtain the time evolution of the collisional dissipated energy, the plastic dissipated energy, the kinetic fluctuation energy and the elastic stored energy (Fig. 6). Collisional dissipated energy is significantly higher than the plastic dissipated energy, meaning that collisions between grains is a key phenomenon in this problem. The kinetic fluctuation energy and the elastic stored energy are very small compared to the plastic and collisional dissipation. The kinetic fluctuation energy increases up to a peak and then tends to zero under static conditions. In addition to the form of internal energy shown in Figure 6, it is also interesting to monitor the evolution of the potential energy (EP) and the kinetic energy (EK) of the system. There are computed by summing up the potential energy and the kinetic energy of each MP according to Equations 17 and 18 respectively: (17) (18) where mp=mass of the MP, vp=velocity, yp=elevation, N=total number of MP, g=gravity acceleration. The potential energy of the system decreases with time (Fig. 7). This decrement of potential energy (EP) obviously coincide with the work of the external forces. The potential energy transforms in kinetic energy and internal energy of the system, the latter is partially stored and partially dissipate as previously discussed. In contrast, the kinetic energy of the system increases up to a peak and then tends to zero. The peak of kinetic energy occurs approximately at time , where corresponds to the time taken by a single particle in free fall to travel from the centre of the column to the base (Fig. 7). Figure 6 Evolution of internal energy types. 547