3. Material Point Method

Similar to other numerical methods, the Material Point Method (MPM) [Sulsky1994] transforms a set of partial differential equations (PDE) into a system of linear algebraic equations. We explain our reasoning for choosing this method and our reasoning for avoiding others. Note that the derivation of traditional MPM presented here is similar to previous descriptions used by MPM researchers in natural hazards engineering (Carter Mast [Carter2013Dissertation], Wen-Chia Yang [Yang2016Dissertation], and Wookuen Shin [Shin2009Dissertation]) to maintain nomenclature consistency.

For the problem of large-deformation debris-fluid-structure interaction (DFSI), found in fluid-driven debris-fields (e.g. tsunamis, avalanches, landslides), the shortlist benefits of MPM appear ideal (large-deformation, nonlinear, multi-material, multi-phase). In theory, it is, although limitations are found. Some of these are alleviated algorithmically while others are a property of the hardware that computation is performed on. The former is to be addressed in part during discussion / derivation of modern numerical improvements in the methodology, such as the Affine Particle-in-Cell (APIC) and Moving Least Squares Material Point Method (MLS-MPM).

Loosely, a Material Point Method cycle can be described graphically by Fig. 3.1 and textually as follows: Space is discretized by a grid, where grid-nodes hold momentum. This introduces an Eulerian scheme. Material bodies are discretized by points (interchangeably referred to as particles), which hold positions and deformation gradients, as well as material parameters and velocities when necessary. This introduces a Lagrangian scheme. The equation of motion and boundary conditions are solved on the grid nodes. Material laws and advection are carried out on particles. Results of these computations pass between the particles and grid through a transfer-scheme (via shape-functions). MPM is thereby defined as a hybrid Euler-Lagrangian numerical method.

Fig. 3.1 Material Point Method numerical steps as a flowchart (Source: CB-Geo MPM)

In this context, the Material Point Method (MPM) is a numerical technique that is best suited for modeling history dependent materials in a dynamic, large deformation setting. The formulation tracks moving points relative to stationary nodes, and can be used to capture the behavior of both fluids and solids in a unified framework. The standard, or traditional implementation solves the governing equation of motion at fixed nodes that collectively form a grid.

Each body or phase in the analysis is represented by a collection of discrete points known as material points or particles. This general concept is shown in the upper portion of Fig. 3.1. Here the different components that make up an MPM simulation are classified as either Body/Phase-Based or Domain/Grid-Based, and properly understanding the role and interaction of these two categories will prove beneficial.

The Body/Phase-Based group is comprised of the continuum body itself and the computational points that, collectively, describe the object. Each particle represents a portion of the total mass, and thus caries an implied volume as well as various state variables depending on the application. For example, in solid mechanics each material point is assigned initial values for position, velocity, stress, strain, and any other state variable needed for the constitutive relationship.

That being the case, these particles form a Lagrangian frame of reference from which the state of the body is determined at any instant in time. A crucial and fundamental characteristic of the Material Point Method is the following: these objects serve only as integration points for the governing equation in space (more on this shortly). As such, these entities are not physical particles, per se, and casting the MPM as a particle-based method is arguably incorrect or, at a minimum, misleading. This problem is compounded by two additional factors. The first is visualization, where viewing dynamic MPM simulations as moving material points is questionable if the true nature of the computational particles is not understood. Even simple effects like adjusting the plotting point size can severely distort the interpretation of what a material point is, can give false impressions of the implicit volume being occupied by the point, and lead to erroneous conclusions regarding the particle nature of the method.

The second factor is the application of the method to discretely-based matter, such as granular materials, where the very nature of the medium lends itself to a particle description. In these applications the computational particles serve as integration points in space that contain the state of the continuum representation at that particular location, and not the individual constituent grains that make up the material. Those readers familiar with Gaussian quadrature in the FEM can relate to the computational particle’s role in terms of Gauss points—the spatial location(s) in which definite integrals arising from the governing equations are approximated over each element. From this perspective particles may be viewed as moving Gauss points that can travel from element to element.

The second category of objects, the Domain/Grid-Based entities, are responsible for defining the physical space a body moves in. The primary object is a node, and a collection of nodes forms a grid. Typically nodes are arranged in a repeating and regular pattern, forming a line in 1-dimensional (1D) applications, a rectangular pattern in 2D simulations, or a rectangular cuboid in fully 3D environments. This repeatable structure is not a requirement of the method but is the most common scheme to date. The nodal arrangement also defines cells, or the region contained between adjacent nodes, as well as nodal supports, defined by piece-wise continuous shape functions residing at each node in the domain. The interplay between the grid, nodes, cells, and nodal support (assuming linear shape functions, which will not be true later on) is shown in Fig. 3.1. Strictly speaking the nodal positions are arbitrary and can potentially change without penalty at any point in time. However, nodes are most commonly assumed to reside in a single location effectively creating a static grid. This facilitates an Eulerian frame of reference when viewed relative to the particle motion. The sharing of information between particles and nodes is governed solely by the shape function that serve as an effective weight for determining the importance of a given particle to any node in the domain. In general this process is referred to as mapping, and can occur from particle-to-node, or in the opposite direction, from node-to-particle.

The primary goal of any analyses is to track the system in time while monitoring the evolution of key quantities in both the Body and Domain categories. This is accomplished by splitting a finite time increment into many smaller time intervals, \(\Delta t=t_{n+1}^{}-t_{n}^{}\), and approximating key equations over each \(\Delta t\). When the governing equation is conservation of linear momentum, material point quantities of mass, momentum, and force are mapped to the appropriate nodes as indicated in Fig. 3.1. After collecting contributions from all particles in the support, the nodal acceleration and velocity vectors are determined over \(\Delta t\) as observed in Fig. 3.1 (b). The velocity gradient and the corresponding strain increment are mapped to the particle location using the updated nodal velocity. The particle stress and material state variables are computed from the desired constitutive model as part of the third step highlighted in Fig. 3.1 (c). Finally, the incremental changes in nodal velocity and position are mapped from the nodes to the particles, resulting in a fully updated system at the particle level. After Fig. 3.1 (d) the procedure begins again and the computational cycle is repeated for a prescribed time duration.

3.1. Derivation

The traditional approach [Sulsky1993] is built around conservation of linear momentum, which when expressed in differential form appears as follows:

(3.1.1)\[\rho\,\dot{\mathbf{v}} = \text{div}{\mathbf{\sigma}} + \mathbf{b} \,\,,\]

with the mass density \(\rho(\mathbf{x},t)\) at position \(\mathbf{x}\) and time \(t\), \(\dot{\mathbf{v}}(\mathbf{x},t)\) as the material time derivative of the velocity field—also known as the acceleration field. Stress divergence \(\text{div}{\mathbf{\sigma}} = \text{div} \mathbf{\sigma}\) , where \(\mathbf{\nabla}\) the gradient operator, \(\mathbf{\sigma}(\mathbf{x},t)\) is the Cauchy stress tensor. \(\mathbf{b}(\mathbf{x},t)\) is the body force per unit volume. In the present derivation the end goal is to obtain an expression for \(\dot{\mathbf{v}}(\mathbf{x},t)\) consistent with the description of the MPM given in the previous section. Thus, it is necessary to build an approximation for the acceleration field in terms of the nodes and particles that make up a given analysis. For this purpose a weighted integral statement is constructed from (3.1.1) as

(3.1.2)\[\int_{V_{\mathcal{B}}}\, \left( \rho\,\dot{\mathbf{v}} - \text{div}{\mathbf{\sigma}} - \mathbf{b} \right) \cdot\mathbf{\eta} \,dV = 0\,\,,\]

effectively transferring the strict, or strong, requirements of (3.1.1) to a weighted statement known as a weak form. Here the integration domain is over the spatial volume \(V_{\mathcal{B}}\) of a continuous body, \(\mathcal{B}\). The vector field \(\mathbf{\eta}(\mathbf{x},t)\) is an arbitrary vector-valued spatial function that is kinematically consistent with the desired boundary conditions. Separating each term according to

(3.1.3)\[\int_{V_{\mathcal{B}}}\, \rho\, \dot{\mathbf{v}}\cdot\mathbf{\eta} \,dV -\int_{V_{\mathcal{B}}}\, \text{div}{\mathbf{\sigma}} \cdot\mathbf{\eta} \,dV -\int_{V_{\mathcal{B}}}\, \mathbf{b} \cdot\mathbf{\eta} \,dV = 0\,\,,\]

and using the product rule of differentiation yields

(3.1.4)\[-\int_{V_{\mathcal{B}}}\, \text{div}{\mathbf{\sigma}} \cdot\mathbf{\eta} \,dV = -\int_{V_{\mathcal{B}}}\, \text{div}{\left(\mathbf{\sigma}\cdot\mathbf{\eta}\right)} \,dV +\int_{V_{\mathcal{B}}}\, \mathbf{\sigma}:\mathbf{\nabla}_{}^{s}\mathbf{\eta} \,dV\,\,.\]

The modified form produces a term that can readily be transformed via the divergence theorem as

(3.1.5)\[\int_{V_{\mathcal{B}}}\, \text{div}{\left(\mathbf{\sigma}\cdot\mathbf{\eta}\right)} \,dV = \int_{\mathcal{S}}^{}\, \left(\mathbf{\sigma}\cdot\mathbf{n}\right) \cdot\mathbf{\eta} \,d\mathcal{S} = \int_{\mathcal{S}_{\mathbf{\sigma}}}^{}\, \tilde{\mathbf{t}} \cdot\mathbf{\eta} \,d\mathcal{S} + \int_{\mathcal{S}_\mathbf{u}}^{}\, \left(\mathbf{\sigma}\cdot\mathbf{n}\right) \cdot\mathbf{\eta} \,d\mathcal{S} \,\,,\]

where \(\mathcal{S}\) is the surface of the body \(\mathcal{B}\) (sometimes written as \(\mathcal{S}=\partial V_{\mathcal{B}}\) in the literature) and \(\mathbf{n}\) is the outward normal defined on \(\mathcal{S}\). The terms \(\mathcal{S}_{\mathbf{\sigma}}\) and \(\mathcal{S}_\mathbf{u}\) correspond to the portions of the surface \(\mathcal{S}\) where loads and displacements are prescribed, respectively. These subsets collectively form the entire surface and do not overlap. The latter statement is summarized concisely as \(\mathcal{S} = \mathcal{S}_{\mathbf{\sigma}}\cup\mathcal{S}_\mathbf{u}\) and \(\mathcal{S}_{\mathbf{\sigma}}\cap\mathcal{S}_\mathbf{u}=0\). The term \(\tilde{\mathbf{t}}=\mathbf{\sigma}\cdot\mathbf{n}\) is a prescribed traction vector residing on the surface \(\mathcal{S}_{\mathbf{\sigma}}\). Requiring that \(\mathbf{\eta} = \mathbf{0}\) on \(\mathcal{S}_\mathbf{u}\) removes the last integral and the remaining terms are collected to form

(3.1.6)\[\int_{V_{\mathcal{B}}} \dot{\mathbf{v}} \cdot \mathbf{\eta}\,\rho\,dV =-\int_{V_{\mathcal{B}}} \mathbf{\sigma}:\mathbf{\nabla}_{}^{s}\mathbf{\eta}\, dV +\int_{V_{\mathcal{B}}} \mathbf{b} \cdot \mathbf{\eta}\, dV +\int_{\mathcal{S}_{\mathbf{\sigma}}} \, \tilde{\mathbf{t}} \cdot \mathbf{\eta}\, d\mathcal{S} \,\,,\]

the very foundation of the MPM approximation scheme—not to mention several other numerical techniques (e.g. Finite Element Method). In the current form two key items need to be addressed: the arbitrary vector-valued spatial function, \(\mathbf{\eta}(\mathbf{x},t)\), and the integration procedure for each term arising in (3.1.6). These items are discussed sequentially in what follows.

The governing equations are solved at nodal points in the domain. That being the case it makes sense to build the unknown field quantities \(\dot{\mathbf{v}}(\mathbf{x},t)\) and \(\mathbf{\eta}(\mathbf{x},t)\) using the nodes themselves. These approximations are constructed as

(3.1.7)\[\mathbf{\eta}(\mathbf{x},t)\approx\mathbf{\eta}^h_{}(\mathbf{x},t) := \sum_{i}^{}\,N_{i}^{}(\mathbf{x})\,\mathbf{\eta}_{i}^{}(t) \quad \text{ and } \quad \dot{\mathbf{v}}(\mathbf{x},t)\approx\dot{\mathbf{v}}^h_{}(\mathbf{x},t) := \sum_{j}^{}\,N_{j}^{}(\mathbf{x})\,\dot{\mathbf{v}}_{j}^{}(t)\]

where \(N_{i}^{}(\mathbf{x})\) and \(N_{j}^{}(\mathbf{x})\) are the shape functions associated with nodes \(i\) and \(j\), respectively. \(\mathbf{\eta}_i^{}(t)\) is an arbitrary, time-dependent nodal vector at a node \(i\), and \(\dot{\mathbf{v}}_{j}(t)\) is the time-dependent nodal acceleration vector at a node \(j\). In this work the superscript \(h\) indicates a grid-based approximation. Closer inspection of the second integral term in Equation (3.1.6) reveals that \(\mathbf{\eta}^h_{}(\mathbf{x},t)\) must be sufficiently smooth in order to accommodate non-zero action of the differential operator, \(\mathbf{\nabla}\). Thus, at a very minimum, the shape functions \(N_{}^{}(\mathbf{x})\) must be linear in \(\mathbf{x}\) (at least \(C^0_{}\) continuous).

The next task is to identify an approximation scheme for the integral terms in (3.1.6) Representing the total body as a collection of particles of fixed mass \(m_p\) not only satisfies conservation of mass, but also allows integrals to be computed as sums over particles as follows:

(3.1.8)\[\begin{split}\begin{eqnarray} \int_{V_{\mathcal{B}}} \left( \bullet \right)\,\rho\,dV &= \sum_{p}^{}\,\int_{V_{p}} \left( \bullet \right)\,\rho_p\,dV_p \\ &= \sum_{p}^{}\,\int_{m_{p}} \left( \bullet \right)\,dm_p \\ &\approx \sum_{p}^{}\,\left( \bullet \right)_p\,m_{p} \,\,. \end{eqnarray}\end{split}\]

The symbol \(\sum_{p}^{}\) indicates a summation over all particles while the subscript \(p\) refers to a particle quantity. The approximation leading to the last term in (3.1.8) may be viewed either as a direct application of the mean value theorem of integration or as a single point numeric integration over the particle domain. This form is contingent upon the transformation to a mass element, defined as \(dm=\rho\,dV\). Comparing to (3.1.6), the proper mass element exists only for the first term and the other terms must be modified appropriately. The notion of a mass-specific term is introduced using the notation \((\bar{\bullet})\), which indicates the transformation of a volume-specific quantity, i.e., \(\left(\bullet\right) = \rho\,(\bullet/\rho) = \rho\,(\bar{\bullet})\) to its mass-specific counterpart. In the present example this transforms the weak form equation to

(3.1.9)\[\int_{m_{\mathcal{B}}} \dot{\mathbf{v}} \cdot \mathbf{\eta} \,dm = -\int_{m_{\mathcal{B}}} \bar{\mathbf{\sigma}}:\mathbf{\nabla}_{}^{s}\mathbf{\eta}\, dm +\int_{m_{\mathcal{B}}} \bar{\mathbf{b}}\cdot\mathbf{\eta} \, dm +\int_{\mathcal{S}_{\mathbf{\sigma}}} \, \tilde{\mathbf{t}}\cdot\mathbf{\eta} \, d\mathcal{S} \,\,,\]

where \(\bar{\mathbf{\sigma}}\) and \(\bar{\mathbf{b}}\) are the mass-specific Cauchy stress and body force, respectively. The primary integration domain has been transformed from the body volume \(V_{\mathcal{B}}\), to the body mass \(m_{\mathcal{B}}\).

3.2. Constructing the System of Equations

The discrete set of equations

(3.2.1)\[\sum_j\,m_{ij}^{}\, \dot{\mathbf{v}}_j^{} = \mathbf{f}_{i}^{int} + \mathbf{f}_{i}^{ext} \,,\]

with

(3.2.2)\[\begin{eqnarray} \mathbf{f}_{i}^{int} = -\sum_p\,\bar{\mathbf{\sigma}}_p^{} \cdot \mathbf{\nabla} N_{ip}^{} \, m_p^{}\,. \hspace{0.375in} \textrm{and} \hspace{0.375in} \mathbf{f}_{i}^{ext}= \sum_p\, \bar{\mathbf{b}}_p^{}\,N_{ip}^{}\,m_p^{} + \int_{\mathcal{S}_{\mathbf{\sigma}}} \, \tilde{\mathbf{t}}\,N_{ip}^{} \, dS \end{eqnarray}\]

is obtained for the unknown nodal accelerations \(\dot{\mathbf{v}}_j^{}\) by substituting the grid-based definitions given in listing (3.1.7) and the integral approximation scheme outlined in (3.1.8) into the weak form (3.1.9). The resulting system utilizes \(m_{ij}^{}=\sum_p\,N_{ip}^{}\,N_{jp}^{}\,m_p^{}\), the consistent mass matrix coefficients with \(N_{ip}^{}\) as the shape function evaluated at the particle location, i.e., \(N_{ip}^{} = N_i^{}(\mathbf{x}_p^{})\). Frequently the off diagonal coupling terms in \(m_{ij}^{}\) are eliminated by approximating the mass matrix as a purely diagonal matrix: \(m_{i}^{}=\sum_p\,N_{ip}^{}\,m_p^{}\). In doing so the system in (3.2.1) is reduced to a series of \(i\) uncoupled equations for the \(i\) nodes describing the spatial domain.

The external surface force term in (3.2.2) can be problematic in the MPM. The root of the problem lies in the fact that surface tractions must be applied on the body—a.k.a. the particles—and these objects move throughout nodal supports in time. Thus, the particle area and force orientation must be tracked appropriately so these terms can be applied at the correct nodes for any given position/orientation of the particle/surface. This is in contrast to other techniques, such as the Finite Element Method (FEM), where this term is applied directly to nodal values.

The primary goal of any analyses is to track the system in time while monitoring the evolution of key quantities at both the particle and nodal levels. This is accomplished in part by splitting a finite time increment, \(T\), into many smaller time intervals, \(\Delta t=t_{}^{n+1}-t_{}^{n}\ll~T\). Over each time step \(\Delta t\) the current state is mapped to the nodes, a grid-based time integration is performed, and particle values are updated. In this section the details of each step are presented. The computational cycle is broken down and visualized as individual components in Fig. 3.1.

The first step involves the transfer of particle quantities to the nodes. This is shown in Fig. 3.1 (a) and is accomplished by way of

(3.2.3)\[\begin{split}\begin{eqnarray} \mathbf{p}_{i}^{n}=\sum_p N_{ip}^{}\,m_p^{}\,\mathbf{v}_{p}^{n} \hspace{0.125in} \textrm{,} \hspace{0.125in} \mathbf{f}_{i}^{int} = -\sum_p\,\bar{\mathbf{\sigma}}_p^{} \cdot \mathbf{\nabla} N_{ip}^{} \, m_p^{}\, \nonumber \\ \hspace{0.25in} \textrm{and} \hspace{0.25in} \mathbf{f}_{i}^{ext} = \sum_p\, \bar{\mathbf{b}}_p^{}\,N_{ip}^{}\,m_p^{} + \int_{\mathcal{S}_{\mathbf{\sigma}}} \, \tilde{\mathbf{t}}\,N_{ip}^{} \, dS \end{eqnarray}\end{split}\]

for the momentum \(\mathbf{p}_{i}^{n}\), internal force \(\mathbf{f}_{i}^{int}\), and external force \(\mathbf{f}_{i}^{ext}\) contributions, respectively. These values are used to solve the linear systems

(3.2.4)\[\begin{eqnarray} \dot{\mathbf{v}}_i^{} = \frac{1}{m_{i}}\left(\mathbf{f}_{i}^{int} + \mathbf{f}_{i}^{ext}\right) \hspace{0.25in} \textrm{and} \hspace{0.25in} \mathbf{v}_{i}^{n} = \frac{\mathbf{p}_{i}^{n}}{m_{i}} \,\,, \end{eqnarray}\]

which yields the acceleration and velocity (\(\dot{\mathbf{v}}_i^{}\) and \(\mathbf{v}_{i}^{n}\)) at time \(t^{n}\) for each node in the domain. Here it is assumed that the consistent mass matrix is approximated using a diagonal, or lumped mass matrix as explained in the previous section. For the explicit integration scheme the nodal acceleration is assumed constant over the time step, resulting in the updated velocity field

(3.2.5)\[\mathbf{v}_{i}^{n+1} =\mathbf{v}_{i,n}^{} + \Delta\mathbf{v}_i^{} \hspace{0.25in} \textrm{with velocity increment} \hspace{0.25in} \Delta\mathbf{v}_i^{}= \Delta t\,\dot{\mathbf{v}}_{i}^{}\]

describing the total field at the end of the time step \(\Delta t\). The velocity field at the beginning and end of each time step are used to define the effective nodal velocity

(3.2.6)\[\mathbf{v}_{i}^{n+\theta} = (1 - \theta)\,\mathbf{v}_{i}^{n} + \theta\,\mathbf{v}_{i}^{n+1} \,\,,\]

where \(\theta\in[0,1]\) is an integration parameter that extracts the field at an arbitrary time, \(t_{n+\theta} = t_n^{} + \theta\Delta t\), between or at \(t_{n}^{}\) and \(t_{n+1}^{}\). The effective velocity gives way to the position increment according to

(3.2.7)\[\Delta\mathbf{x}_i^{} = \mathbf{v}_{i}^{n+\theta}\,\Delta t \,\,.\]

The series of computations outlined in (3.2.3) – (3.2.7) are depicted in Fig. 3.1 (b) and together form the grid-based time integration portion of the MPM analysis. This series of nodal equations implies the nodes themselves are moving. Strictly speaking this statement is true. However, as noted previously, the nodal position is arbitrary at the beginning of each time step. It is common practice to continuously assume nodal positions coincide with their original position at \(t=t^{0}\) for the start of each new time step. This may be interpreted as discarding the old grid and creating a new series of nodes each time step.

At this stage in the computational cycle the motion at the nodes is well defined over the time step and will not change. Therefore the resulting deformation, incurred in an incremental fashion as a result of the change in motion, is determined based on the last known state. This stage is represented in Fig. 3.1 (c). The velocity gradient is computed at the particle level according to

(3.2.8)\[\mathbf{\nabla}\mathbf{v}_{p}^{h,n+\theta} = \sum_i \mathbf{v}_{i}^{n+\theta}\otimes\mathbf{\nabla}\,N_{ip} \,\,.\]

Multiple deformation, or strain, measures exist depending on the type of analysis. A single presentation cannot possibly accommodate all the options in this regard. This section will focus on a large deformation measure obtainable from the incremental deformation gradient:

(3.2.9)\[\mathbf{f}_{p} = \mathbf{1} + \Delta t\,\mathbf{\nabla}\mathbf{v}_{p}^{h,n+\theta} \,\,.\]

The particle strain is updated according to

(3.2.10)\[\mathbf{\varepsilon}_{p}^{n+1} = \tilde{\mathbf{\varepsilon}}(\mathbf{\varepsilon}_{p}^{n}, \mathbf{f}_{p}^{})\]

where \(\tilde{\mathbf{\varepsilon}}\) is a general strain function of the known deformation state at \(t^{n}\) and the incremental change over the time step. The strain function can take many forms and be dependent on other quantities, but this is beyond the scope of the current section.

The particle stress is determined from

(3.2.11)\[\bar{\mathbf{\sigma}}_{p}^{n+1} = \frac{\partial \bar{\varPsi}(\mathbf{\varepsilon}_{p}^{n+1})}{\partial\mathbf{\varepsilon}_{p}^{n+1}} \,\,,\]

where \(\bar{\varPsi}\) is the mass specific free energy function. Much like the strain function, \(\tilde{\mathbf{\varepsilon}}\), the free energy function is typically cast in terms of several additional variables, including state dependent quantities required for elastoplastic constitutive laws. For the time being these details are omitted. The key point to take from this presentation is the stress is a function of the updated strain. This implies a hierarchical structure that may be exploited. For the special case of a linear elastic material, the particle stress is obtained as

(3.2.12)\[\bar{\mathbf{\sigma}}_{p}^{n+1} = \frac{K}{\rho_{0}} (\text{tr}{\mathbf{\varepsilon}_{p}^{n+1}})\, \mathbf{1} + \frac{2\,G}{\rho_{0}}\text{dev}{\mathbf{\varepsilon}_{p}^{n+1}}\]

with initial mass density \(\rho_{0}\), bulk modulus \(K\), and shear modulus \(G\). The terms \(\text{tr}(\bullet)\) and \(\text{dev} (\bullet)\) are the standard trace and deviatoric operators on a second order tensor.

Depicted in Fig. 3.1, the final portion of the computational cycle is the particle velocity and position update

(3.2.13)\[\mathbf{v}_{p}^{n+1} = \mathbf{v}_{p}^{n} + \sum_i\,N_{ip}^{}\,\Delta \mathbf{v}_i^{} \hspace{0.25in} \textrm{and} \hspace{0.25in} \mathbf{x}_{p}^{n+1} = \mathbf{x}_{p}^{n} + \sum_i\,N_{ip}^{}\,\Delta \mathbf{x}_i^{} \,\,,\]

obtained from the incremental change in nodal velocity and position over the time step. Upon completion of this last step the cycle repeats until the analysis time reaches a user prescribed value.

The details provided here highlight the very basics of the Material Point Method. Of course, this traditional form is subject to change depending on the implementation strategy or details arising due to the extension of the traditional framework.

3.3. Numerical integration

The conservation equation for any given general scalar variable \(\phi\) can be given as

(3.3.1)\[\underbrace{\frac{\partial \left( \rho \phi \right)}{\partial t}}_{\text{Transient term}} + \underbrace{\nabla \cdot \left( \rho \mathbf{v} \phi \right)}_{\text{Convective term}} = \underbrace{\nabla \cdot \left( \Gamma^{\phi} \nabla\phi \right)}_{\text{Diffusion term}} + \underbrace{Q^{\phi}}_{\text{Source term}}\]

where \(\rho\) is the density of the fluid, \(\mathbf{v}\) is the velocity vector. Considering the quasi-static flow, by dropping the transient term, and integrating over the volume of an element, we have

(3.3.2)\[\int_{\Omega} {\nabla \cdot \left(\rho\mathbf{v}\phi\right) \ d\Omega} = \int_{\Omega} {\nabla \cdot \left( \Gamma^{\phi} \nabla\phi \right) \ d\Omega} + \int_{\Omega} {Q^{\phi} \ d\Omega}\]

Using the divergence theorem, the above can be re-written in terms of surface integrals as

(3.3.3)\[\int_{\Gamma} {\left(\rho\mathbf{v}\phi\right) \ d\Gamma} = \int_{\Gamma} {\left( \Gamma^{\phi} \nabla\phi \right) \ d\Gamma} + \int_{\Omega} {Q^{\phi} \ d\Omega}\]

The above integral form requires the evaluation of the flux integration over the elemental faces for the convective and diffusion terms. This can be alternatively written as the summation of the flux terms over each of the individual faces of the element, i.e.

(3.3.4)\[\begin{split}\begin{split} \int_{\Gamma} {\left(\rho\mathbf{v}\phi\right) \ d\Gamma} &= \sum_{i = 1}^{n\left(\Omega\right)} \left( \int_{\Gamma_i} {\left(\rho\mathbf{v}\phi\right) \ d\Gamma} \right) \\ \int_{\Gamma} {\left( \Gamma^{\phi} \nabla\phi \right) \ d\Gamma} &= \sum_{i = 1}^{n\left(\Omega\right)} \left( \int_{\Gamma_i} {\left( \Gamma^{\phi} \nabla\phi \right) \ d\Gamma} \right) \end{split}\end{split}\]

where \(n\left(\Omega\right)\) represents the number of faces of the element with volume \(\Omega\), \(\Gamma_i\) represents the \(i\)-th face of the element with volume \(\Omega\). The above form of discretization ensures the conservation of quantities. It is important to note that the quantities of interest are conservative in nature. These include mass, volume, energy etc. Thus, without the convervative properties, the overall solution process can lead to unphysical solutions. In other words, the flux across the face of two shared elements need to have equal magnitudes but of opposite signs. The flux leaving the face of the first element should be equal to the flux entering, through the face, into the second element.

The resulting integral equations, shown above, include surface integral over each face of the element. These integral equations need to be converted to algebraic equations and are hence further simplified using the Gaussian quadrature as

(3.3.5)\[\int_{\Gamma_i} {\left( \alpha \right) \ d\Gamma} = \sum_{p} {\alpha_p w_p} A_{\Gamma_i}\]

where \(\alpha\) represents the quantity of interest (here the advection or diffusion term), \(A_{\Gamma_i}\) represents the area of the face, \(\alpha_p\) represents the quantity of interest at the \(p-th\) integration point, \(w_p\) represents the weight at the \(p-th\) integration point. The accuracy of the integration depends on the number of integration points used. In the case of a 2-D problem, the faces are 1-D line units and the integration points are given as

• One-integration point (or also known as Trapezoidal rule): \(\xi_{p} = 1\) and \(w_{p} = 1\)

• Two-integration points: \(\xi_{1} = \left(3-\sqrt{3}\right)/6, \ \xi_{2} = \left(3+\sqrt{3}\right)/6\) and \(w_{1} = w_{2} = 1/2\).

• Three-integration points: \(\xi_{1} = \left(5-\sqrt{15}\right)/10, \ \xi_{2} = 1/2, \ \xi_{3} = \left(5+\sqrt{15}\right)/10\) and \(w_{1} = 5/18, \ w_{2} = 4/9, \ w_{3} = 5/18\).

Similarly, the volume integration of the source term can be achieved using the Gaussian quadrature. Similar to the integration over the boundary, a volume integral need to be performed.

Once the PDE’s have been converted to a summation form above, it is necessary to express the face and volume fluxes in terms of the values of the variable at the neighboring cell centers or in otherwods, we need to linearize the fluxes.

Points to note

• One other aspect that needs to be considered is the effect of hydrostatic pressure. Thus, one also needs to consider which direction the gravity is pointed in.

3.4. References

MoMaDa2016

6. Moukalled, L. Mangani and M. Darwish, “The finite volume method in computational fluid dynamics,” Springer International Publishing Switzerland (2016)

HiNi1981

3. 23. Hirt and B. D. Nichols , “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39(1), pp. 201-225 (1981)

Beetal2009

5. Berberovic, N. P. Van Hinsberg, S. Jakirlic, I. V. Roisman, C. Tropea, “Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution,” Physical Review E, vol. 79 (2009)

Carter2013Dissertation

Mast, C. M. (2013). Modeling Landslide-Induced Flow Interactions with Structures using the Material Point Method. PhD thesis. University of Washington, Seattle, WA. https://digital.lib.washington.edu/researchworks/bitstream/handle/1773/23580/Mast_washington_0250E_11795.pdf?sequence=1

Shin2009Dissertation

Shin, W. (2009). Modeling Mixing and Separation of Solid Matter and Fluid in Landslides and Debris Flows by Representing the Multiphase Material through Distinct Phases. PhD thesis. University of Washington, Seattle, WA.

Sulsky1993

Sulsky, D., Chen, Z., & Schreyer, H. L. (1993). A particle method for history-dependent materials. Other Information: PBD: Jun 1993. https://doi.org/10.2172/10177049

Sulsky1994

Sulsky, D., Chen, Z., & Schreyer, H. L. (1994). A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering, 118(1), 179-196. https://doi.org/https://doi.org/10.1016/0045-7825(94)90112-0

Yang2016Dissertation

Yang, W.-C. (2016). Study of Tsunami-Induced Fluid and Debris Load on Bridges using the Material Point Method. PhD thesis. University of Washington, Seattle, WA.