An implicit Material Point Method applied to granular flows

Procedia Engineering 175 ( 2017 ) 226 – 232 1877-7058 © 2017 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of the 1st International Conference on the Material Point Method doi: 10.1016/j.proeng.2017.01.017 ScienceDirect Available online at www.sciencedirect.com 1st International Conference on the Material Point Method, MPM 2017 An implicit material point method applied to granular flows Ilaria Iaconetaa,∗, Antonia Laresea, Riccardo Rossia, Eugenio O˜ natea aInternational Center for Numerical Methods in Engineering (CIMNE), Technical University of Catalonia (UPC), Edificio C1, Campus Norte, Jordi Girona 1-3, 08034 Barcelona, Spain Abstract The main objective of this work lies in the development of a variational implicit Material Point Method (MPM), implemented in the open source Kratos Multiphysics framework. The ability of the MPM technique to solve large displacement and large deformation problems is widely recognised and its use ranges over many problems in industrial and civil engineering. In the current work the continuum based implicit MPM is applied to engineering applications, where granular material flow is involved. For the resolution of the length and time scale of these particular problems, both continuum and discrete models are typically used. Even if discrete techniques predict more feasible results, nowadays, their use is limited to the investigation of element tests of particles, or to the simulation of reduced systems, not allowing to make important decisions in the analysis and design of granular processes. Some advantages of MPM over discrete methods are tested, such as, the ability to simulate granular flow at the large scale with acceptable computational cost and the capability to get information of stress and strain state in a more straightforward way. The focus of this paper is a comparative study between an irreducible and a mixed formulation, both implemented in the MPM code, to assess the improvement in accuracy and reliability of the numerical results when the latter formulation is adopted. c 2016 The Authors. Published by Elsevier B.V. Peer-review under responsibility of the organizing committee of the 1 st International Conference on the Material Point Method. Keywords: material point method; nonlinear finite element method; implicit MPM; granular flows. 1. Introduction The numerical simulation of solid mechanics problems involving history dependent materials and large deformations, has historically represented one of the most important topics in computational mechanics. Among the Lagrangian techniques, in the last decades, the Material Point Method (MPM) [1,2] has experienced an increasing popularity due to its capabilities of solving several complex engineering problems. This method has its origin in the particle-in-cell (PIC) method, formulated by Harlow [3] for the resolution of fluid flow problems. Some decades after, Sulsky presented its extension to solid mechanics [1,2]. The Material Point Method combines a Lagrangian description of the body of interest, which is represented by a set of particles, the so-called material points, or MPM particles, with the use of a computational mesh used to solve numerically the system of governing equations. ∗Corresponding author. E-mail address: [email protected] © 2017 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the organizing committee of the 1st International Conference on the Material Point Method 227 Ilaria Iaconeta et al. / Procedia Engineering 175 ( 2017 ) 226 – 232 At each time step, the governing equations are solved on the nodes of the computational grid, while history dependent variables and material information are saved on the particles for the entire deformation process. With few notable exceptions [4–8], the majority of the MPM algorithms are written in an explicit formulation [9– 11]. This approach is generally preferable when simulating impacts at high velocities, or fast transient problems. In other cases, for example when the driving force is the gravity, or when the rate of deformation is small, the adoption of an implicit time scheme is the best choice. In implicit scheme, in fact, the stability of the method (for properly chosen dissipative methods) does not depend on the wave propagation speed within the media, which provides the typical time step limitation for explicit approaches. It is worth mentioning the work of [8] where a comparison between an explicit and implicit time scheme is performed. In this work the advantages of using an implicit scheme are discussed, as for instance, the lower limitation on the time step size and the improvement of the algorithmic accuracy for elastoplastic constitutive law. The MPM code employed in the current work has been developed by the authors within the Kratos Multiphysics open source platform [12,13]; the details of the formulation and the algorithm can be found in [14]. The code proposed is designed for an easy implementation of a wide range of pressure-dependent constitutive laws, chosen in function of the problem to be solved. A preliminary study is performed to test the capability of an irreducible formulation (i.e., with displacement as only primary variable) and a mixed formulation (i.e., with displacement and pressure1as primary variables) to give reliable results when one wants to solve problems which involve granular flow under the assumption of large displacements and large deformations. In fact, it is well known that some issues might appear when an irreducible formulation is adopted [15]. For instance, when elastoplastic constitutive laws in the framework of J2-incompressible plasticity are considered, the numerical results can be affected by locking when the volumetric strain contribution is higher than the deviatoric one. The only way to guarantee a locking free behaviour for incompressible material and mesh independence, is to develop a mixed finite element, able to evaluate the results avoiding this drawback. In the present work a granular column collapse is analysed, using an elastoplastic constitutive law with Drucker- Prager yield criterion. The authors demonstrate how using an irreducible formulation the numerical solution is affected by volumetric locking. Therefore, the adoption of a mixed formulation becomes fundamental; this latter formulation, which accounts for a second primary variable, commonly represented by the pressure, allows to correctly evaluate the strain field when the body experiences only an isochoric deformation and, as a consequence, to avoid some issues such as the volumetric locking. 2. Description of the problem In this section the equations, which describe the physical problem under investigation, are presented. Firstly, the governing equations in strong form are illustrated. Secondly, the procedure to derive the linearised solving system of algebraic equations is briefly described and referenced. Finally, the Drucker-Prager plastic model is presented. 2.1. Governing equation in strong form Let us consider the body Bwhich occupies a region Ωof the three-dimensional Euclidean space Ewith a regular boundary ∂Ωin its reference configuration. A deformation of Bis defined by a one-to-one mapping ϕ:Ω→E (1) that maps each point pof the body Binto a spatial point x x=ϕ(p)(2) which represents the location of pin the deformed configuration of B. The region of Eoccupied by Bin its deformed configuration is denoted as ϕ(Ω). 1In the present work the pressure variable refers to the volumetric stress field. 228 Ilaria Iaconeta et al. / Procedia Engineering 175 ( 2017 ) 226 – 232 The problem is governed by mass and linear momentum balance equations Dρ Dt +ρ∇·v=0inϕ(Ω) (3a) ρa−∇·σ=ρbin ϕ(Ω) (3b) where ρis the mass density, ais the acceleration, vis the velocity, σis the symmetric Cauchy stress tensor and bis the body force. Thermal effects are not considered in the present work, so the energy balance is considered implicitly fulfilled. The balance equations are solved numerically in a three-dimensional region Ω⊆R 3, in the time range t∈[0,T], given the following boundary conditions on the Dirichlet (ϕ(∂ΩD)) and Neumann boundaries (ϕ(∂ΩN)), respectively u=uon ϕ(∂ΩD) (4a) σ·n=ton ϕ(∂ΩN) (4b) where nis the unit outward normal. The constitutive model, which is needed to fully define the boundary value problem is described in detail in Section 2.3. 2.2. Weak form and linearisation of the weak form in spatial form Following the standard FEM procedure, the weak form of the momentum balance equation is obtained by employing the Galerkin method. The L2inner product of Equation 3b is derived using an arbitrary test function w, such that w={w∈V|w=0on ϕ(∂ΩD)}, where Vis the space of virtual displacements. By using the divergence theorem the weak form of momentum balance can be obtained and expressed as G(u,w)=ϕ(Ω) σ:(∇Sw)dv −ϕ(Ω) ρ(b−a)·wdv −ϕ(∂ΩN) t·wda =0,∀w∈V (5) The previous equation has the same expression of a weak form under the assumption of infinitesimal strains and displacements. However, in this work material and geometric non-linearities are considered; thus, a linearisation of the weak form is needed. An expansion in Taylor’s series of Equation 5, evaluated at the last known equilibrium configuration u∗, is performed and expressed as L(δu,w)≃G(u∗,w)+DG(u∗,w)[δu]=0,∀w∈V (6) where Lis the linearised virtual work function and DG(u∗,w)[δu] is the directional derivative of Gat u∗in the direction of δu. The detailed procedure to derive the final expression of the solving system of linearised equations (Equation 6) in integral and discrete form can be found in [16] and [14]. Concerning the mixed formulation, linear finite elements formulated in a mixed displacement-pressure field are considered. In this regard, Equation 5 can be rewritten as G(u,w)=ϕ(Ω) (dev(σ)+p1):(∇Sw)dv −ϕ(Ω) ρ(b−a)·wdv −ϕ(∂ΩN) t·wda =0 (7) where the Cauchy stress tensor σis decomposed in its deviatoric dev(σ) and volumetric component p. The pressure field pin Equation 7 represents the second primary variable, determined by the volumetric part of the material model. As in the present work a Neo-Hookean material is adopted, the resultant continuity equation is given by p−U(J)=0 (8) where U(J) is the first derivative of the volumetric term of the free energy function ψ, which will be defined in the following section once the constitutive model is characterized. By performing a L2inner product of Equation 8 with 229 Ilaria Iaconeta et al. / Procedia Engineering 175 ( 2017 ) 226 – 232 an arbitrary test function q, such that q={q∈Q|q=0on ϕ(∂ΩD)}, where Qis the space of virtual pressures, the weak form of the pressure constitutive equation can be obtained and expressed as G(p,q)=ϕ(Ω) qp−U(J)dv =0,∀q∈Q (9) For the treatment of the incompressibility constraint a pressure stabilization is adopted. In this work the Polynomial Pressure Projection (PPP), introduced by Dohrmann and Bochev [17], is used. The determination of an additional term, which has to be added to Equation 9 to stabilized the mixed finite element, is explained in detail in [17] and [18]. 2.3. Drucker-Prager law in finite strain plasticity To complete the non-linear boundary valued problem a Drucker-Prager plastic model is considered. For the elastic response a hyperelastic Neo-Hookean model is considered, while the plastic behaviour is modelled following the theory of J2-plasticity at finite strains developed by Simo and Hughes [19]. The algorithm adopted in the current work is the same one presented by Hofstetter and Taylor [20]. The present formulation uses a non-associative flow rule, assuming that the material is plastically incompressible, which means that the dilatation of the plastic strain is equal to zero. Under this assumption, an extension of the return mapping algorithm for J2plasticity [19] is required. As the Neo-Hookean material is represented by the following free energy function ψ ψ(b)=U(J)+W(b)=1 2Kln2(J)+1 2Gdev(¯ b) (10) written here as the sum of volumetric (U) and deviatoric (W) response, the Kirchhoffstress τcan be defined as τ=Kln(J)+1 2Gdev(¯ b) (11) where Kand Gare the bulk and shear moduli, respectively, J=det(F) is the determinant of the total deformation gradient and ¯ bis the volume preserving part of left Cauchy-Green tensor, defined as ¯ b=J−2 3b(12) Let Ψdenote a convex set of plastically admissible stresses Ψ={τ:f(τ)≤0}(13) the Drucker-Prager yield condition in terms of Kirchhoffstress is f(τ)=|dev(τ)|+μ √3tr(τ)−√2κ(14) where μ=2√2sinφ (3+sinφ)is a coefficient depending on the internal friction angle φ, and κthe yield strength in simple shear. The constitutive relation, based on the decomposition of the total deformation gradient F=FeFpinto an elastic and a plastic part , denoted as Feand Fp, respectively, can be written as dτn+1= Cn+1:dn+1−p n+1(15) where p n+1=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ˙ λ∂G(τ) ∂τ,if f(τ)=0 0,if f(τ)<0(16) where  Cis the fourth order incremental constitutive tensor, ˙ λis the plastic multiplier and G is the plastic potential, which in the present work can be expressed as G(τ)=|dev(τ)|. 230 Ilaria Iaconeta et al. / Procedia Engineering 175 ( 2017 ) 226 – 232 3. Numerical example In the present work a preliminary study is performed to test the capability of an irreducible and a mixed formulation, which the implicit MPM code is based on, to obtain reliable results in problems where granular flow is involved under the assumption of large displacements and large deformations. In particular the authors are interested in demonstrating if the results are affected by volumetric locking, a consequence of the inability of the element to exactly represent the isochoric strain field, when an irreducible formulation is adopted and if the adoption of a mixed formulation can fix this issue. In this section a planar granular column collapse test case, under the force gravity, g, is considered. The geometry of the sample can be observed in Fig. 1, characterized by an aspect ratio a=h0 l0=1.5 and a bottom boundary where no-penetration and no-slip conditions are applied. The material and problem data used to run the simulations can be found in Table 1 and Table 2.    Fig. 1. Planar granular column collapse. Geometry. Table 1. Planar granular column collapse. Material data. a=h0/l0Density kg/mc Young’s Modulus kPa Poisson’s ratio Friction coefficient 1.5 2200 200 0.29 0.633 Table 2. Planar granular column collapse. Problem data. MP particles per cell Mesh size [m] Time step [s] Type of geometrical element 12 0.004 0.00001 Triangle In Figure 2 the configurations at the dimensionaless time T=tho l2 0 gare shown. As one can observe, the results obtained by adopting an irreducible formulation (Figure 2(a)) are affected by volumetric locking in all the cases considered; the most visible effect of this issue is represented by the left and right upper corner, which doesn’t disappear during the simulation (Figure 2(c)). On the other hand by using a mixed formulation (Figure 2(b)) the results look to be less affected by locking, mainly visible at T =0.8, when the upper corners are completely diseappered. In the presentation it is shown that by adopting a mixed formulation this issue can be fixed and the quality of the results can be improved, demonstrated by performing a qualitatively comparison of the configuration of the sample at different dimensionless time. 4. Conclusion An implicit MPM is presented and used to solve problems, which involve granular flow, under the assumption of large displacements and large deformations. The code proposed is suitable to implement several constitutive laws; in this work a Drucker-Prager plastic model is used as first attempt to predict the behaviour of granular flow problem. 231 Ilaria Iaconeta et al. / Procedia Engineering 175 ( 2017 ) 226 – 232     (a) Irreducible formulation.     (b) Mixed Formulation. (c) Detail of the side edges of the irreducible (left) and mixed (right) formulations. Fig. 2. Planar granular column collapse. Configuration of the deformed sample at different dimensionless time T. Comparison between the use of an irreducible formulation (Figure 2(a)) and a mixed one (Figure 2(b)); detail of the deformed sample side edges at different dimensionless time T in Figure 2(c). Nevertheless, in the future further pressure-dependent constitutive laws will be considered and tested to improve the accuracy and quality of the numerical results. A preliminary study is performed; a granular column collapse test case is considered and the results obtained by adopting an irreducible and a mixed formulation are compared. It is demonstrated that by using an irreducible formulation it is not possible to avoid the volumetric locking issue under the assumption of plastic incompressibility. In the presentation it is shown that the adoption of a mixed formulation is strictly recommended to avoid locking and to obtain more reliable results. Acknowledgements The research was supported by the Research Executive Agency through the T-MAPPP project (FP7 PEOPLE 2013 ITN-G.A.n607453). References [1] D. Sulsky, Z. Chen, H. Schreyer, A particle method for history-dependent materials, Computer Methods in Applied Mechanics and Engineering 118 (1994) 179–196. [2] D. Sulsky, S.-J. Zhou, H. L. Schreyer, Application of a particle-in-cell method to solid mechanics, Computer Physics Communications 87 (1995) 236–252. [3] F. Harlow, The particle-in-cell computing method for fluid dynamics, Methods for Computational Physics 3 (1964) 319–343. 232 Ilaria Iaconeta et al. / Procedia Engineering 175 ( 2017 ) 226 – 232 [4] J. E. Guilkey, J. A. Weiss, Implicit time integration for the material point method: Quantitative and algorithmic comparisons with the finite element method, International Journal for Numerical Methods in Engineering 57 (2003) 1323–1338. [5] D. Sulsky, A. Kaul, Implicit dynamics in the material-point method, Computer Methods in Applied Mechanics and Engineering 193 (2004) 1137–1170. [6] L. Beuth, Z. Wieckowski, P. Vermeer, Solution of quasi-static large-strain problems by the material point method, International Journal for Numerical and Analytical Methods in Geomechanics 35 (2011) 1451–1465. [7] J. Sanchez, H. Schreyer, D. Sulsky, P. Wallstedt, Solving quasi-static equations with the material-point method, International Journal for Numerical Methods in Engineering 103 (2015) 60–78. [8] B. Wang, P. Vardon, M. Hicks, Z. Chen, Development of an implicit material point method for geotechnical applications, Computers and Geotechnics 71 (2016) 159–167. [9] Z. Wieckowski, The material point method in large strain engineering problems, Computer Methods in Applied Mechanics and Engineering 193 (2004) 4417–4438. [10] X. Zhang, K. Y. Sze, S. Ma, An explicit material point finite element method for hyper-velocity impact, International Journal for Numerical Methods in Engineering 66 (2006) 689–706. [11] P. C. Wallstedt, J. E. Guilkey, An evaluation of explicit time integration schemes for use with the generalized interpolation material point method, Journal of Computational Physics 227 (2008) 9628–9642. [12] P. Dadvand, A framework for developing finite element codes for multi-disciplinary applications., PhD thesis: Universidad Polit´ ecnica de Catalu˜ na, 2007. [13] P. Dadvand, R. Rossi, E. O˜ nate, An object-oriented environment for developing finite element codes for multi-disciplinary applications, Archives of Computational Methods in Engineering 17 (2010) 253–297. [14] I. Iaconeta, A. Larese, R. Rossi, Z. Guo, An implicit grid-based and a meshless mpm formulation for problems in solid mechanics, Computational Particle Mechanics (2016). Submitted to International Journal for Numerical Methods in Engineering. [15] M. Cervera, M. Chiumenti, L. Benedetti, R. Codina, Mixed stabilized finite element methods in nonlinear solid mechanics. part iii: Compressible and incompressible plasticity, Computer Methods in Applied Mechanics and Engineering 285 (2015) 752 – 775. [16] P. Wriggers, Computational Contact Mechanics, Springer, 2006. [17] C. Dohrmann, P. Bochev, A stabilized finite element method for the Stokes problem based on polynomial pressure projections, IJNME 46 (2004) 183–201. [18] O. Zienkiewicz, R. Taylor, J. Zhu, The finite element method: Its basis and fundamentals, in: The Finite Element Method: its Basis and Fundamentals (Seventh Edition), seventh edition ed., Butterworth-Heinemann, Oxford, 2013. doi:http://dx.doi.org/10.1016/B978-1-85617- 633-0.00019-8. [19] J. Simo, T. Hughes, Computational Inelasticity, Springer-Verlag New York, 1998. [20] G. Hofstetter, R. Taylor, Non-associative drucker-prager plasticity at finite strains, Communications in Applied Numerical Methods 6 (1990) 583–589.