Temporal and null‐space filter for the material point method

This paper presents algorithm improvements that reduce the numerical noise and increase the numerical stability of the material point method formulation. Because of the linear mapping required in each time step of the material point method algorithm, a possible mismatch between the number of material points and grid nodes leads to a loss of information. These null‐space related errors may accumulate and affect the numerical solution. To remove the null‐space errors, the presented algorithm utilizes a null‐space filter. The null‐space filter shown removes the null‐space errors, resulting from the rank deficient mapping and the difference between the number of material points and the number of nodes. The presented algorithm enhancements also include the use of the explicit generalized‐α integration method, which helps optimizing the numerical algorithmic dissipation. This paper demonstrates the performance of the improved algorithms in several numerical examples.

PIC and attempt to preserve the information, somewhat reducing energy dissipation. In this paper, a generalized-α dissipative scheme 22 is adopted in the framework of the MPM. The generalized-α scheme simply damps out the high-frequency noises only. Therefore, it avoids the dissipation present in the PIC schemes without cumbersome coupling with other numerical methods.
Furthermore, the MPM formulation requires transfer of the state variables between the material points and the background grid. This frequent mapping between material points and grid nodes can cause significant numerical errors and lead to instabilities, such as ringing instability (presented in detail by Brackbill 3 for PIC methods), which relates to a null-space problem. 23 Similar to the PIC methods, the null-space in the mapping is also presented in MPM, which leads to spurious numerical oscillations in the solution. In general, the null-space related errors may grow continuously over multiple time steps, and, consequently, the numerical solutions may diverge. 23, 24 Tran et al 24 have shown that, if there are more material points than nodes, null-space errors exist despite the choice of the shape function. Therefore, as the null-space existence is due to the larger number of the material points than the grid nodes, the null-space error is present in all the variants of the MPM, including GIMP, 8 DDMP, 9 and B-spline MPM. [10][11][12][13] To remove the null-space errors, Gritton et al 23,25 proposed a null-space filter using the single-value-decomposition (SVD) operator in the 1D formulation of the original MPM. This paper shows how to apply a null-space filter to higher-order interpolation function, for example, in the GIMP or DDMP formulations of the MPM. The presented general formulation also removes the null-space instability at lower computational cost, when compared with the SVD. 23,25 Furthermore, the paper extends the null-space filter to a two-phase hydro-mechanically coupled formulation of MPM. The potential of the proposed formulation is demonstrated in several numerical examples in both one-phase and two-phase simulations.

GENERALIZED-INTEGRATION SCHEME FOR MPM
The MPM algorithm has been described in detail by many researchers 1,4 and summarized in Figure 1. In the MPM framework, the velocity update scheme is based on the acceleration, stemming from the FLIP scheme 26 as follows: The other approach is to overwrite the velocity at material points by the nodal velocity using the PIC scheme 2 v t+1 Jiang et al 19,20 showed that the PIC scheme can filter the spurious modes but exhibits an excessive energy dissipation. In contrast, the FLIP scheme used in MPM is energy conservative, but the spurious modes are preserved and accumulated. Stomakhin et al 27 and Nairn 28 updated the velocity by combining the PIC and FLIP scheme velocity update, with PIC component as a damping factor. In this paper, we apply the explicit generalized-α time integration scheme proposed by Hulbert and Chung, 22 a pure FLIP scheme, but the acceleration is evaluated both in the beginning and the end of the time step. This approach is shown to be able to damp out the high frequency modes only without causing excessive numerical dissipation for structural engineering. 29 Kontoe 30 has shown that the generalized-α algorithm is more accurate and has better numerical dissipation characteristics than other dissipated time integration schemes. 31,32 Moreover, the generalized-α algorithm can mitigate the spurious high frequency mode but does not affect the lower modes strongly. Additionally, it allows controlling the numerical dissipation with a single parameter, which is more appropriate than damping related to the numerical time step size. 29 In the generalized-α algorithm, 33 the acceleration is updated as follows: m t a t+1−α m = f t+1−α f (5) a t+1−α m = (1 − α m ) a t+1 + α m a t (6) where α m and α f are time integration parameters, a t+1−α m is the nodal acceleration evaluated at t = (1 − α m )t t+1 + α m t t , and f = f int + f ext is the sum of internal and external forces and evaluated at t = (1 − α f )t t+1 + α f t t . For the update of displacements and velocities, the original Newmark's equations are where γ and β are time integration parameters. Chung and Hulbert 33 showed that the generalized-α method is secondorder accurate, provided Moreover, the generalized-α method is unconditionally stable, provided To derive the explicit generalized-α method, Hulbert and Chung 22 selected α f = 1 and rederived the numerical algorithmic parameters. Because α f = 1, the explicit scheme becomes conditionally stable and Equation (5) becomes Hilber and Hughes 29 showed that the stability of a time integration scheme depends on the natural frequency ω = √ K∕M, where K and M are stiffness and mass, respectively, and time step is dt. Therefore, Ω = ω * dt is considered to evaluate the stability of the scheme. The explicit generalized-α integration 22 is stable with Ω ∈ [0, Ω s ], where Ω s is the critical limit. The maximum dissipation is obtained at a bifurcation limit Ω b . The spectral radius at the bifurcation limit ρ b is a user-specified parameter to control the algorithmic dissipation. If ρ b is equal to 1, there is no dissipation. When the spectral radius ρ b reduces, the numerical damping increases and gets to maximum values at ρ b = 0 (see Figure 2). The bifurcation limit Ω b and the critical limit Ω s are the function of ρ b as follows: FIGURE 2 Evolution of spectral radius for explicit generalized-α method 22 Other algorithm parameters could be linked to ρ b as follows: Overall, the generalized-α method requires to compute the nodal accelerations twice, first at the beginning a t i and then at the end of the time step a t+1 i . This scheme doubles the amount of computations needed, hence we apply another approach for MPM based on tracking the acceleration a t p at material points. Therefore, each material points have two kinetic variables: velocity v t p and acceleration a t p , which are mapped to the grid as follows: The nodal accelerations a t+1−α m i and a t+1 i and nodal velocities v t+1 i are now computed explicitly as follows: Finally, the velocity, position, and acceleration of material points are updated, ie, The algorithm is controlled by only a single parameter ρ b , as all the other parameters α m , β, γ are automatically updated (see

Choice of shape functions gradient ∇S ip
The original MPM 1 selects the gradient of the shape function ∇S ip such that the sign of the internal force changes when material points cross to a new cell. Consequently, cell-crossing creates unphysical noise and possibly leads to numerical instability. The cell-crossing can be mitigated by using characteristic functions in the GIMP. 8 The characteristic functions in GIMP give the material points domains (red square in Figure 3) and continuous gradients. Another approach, which reduces the cell-crossing errors, is DDMP. 9 The DDMP replaces the gradient ∇S ip with a continuous function between the grid cells while leaving the shape functions intact. The shape function gradient in the cell, which contains the material points, is the shape function gradient of the original MPM, denoted ∇S i , whereas the shape function gradient in the neighbor-node cells is the node-based function, denoted∇S i . The DDMP enlarges the influence domain of material points through dual domains, a cell contains the material points (gray cell in Figure 3) and neighbor-node cells (blue dash cell in Figure 3), without using the characteristic functions to track the domain of the material points.
In DDMP, the gradient of the shape function is modified to with the node-based function∇S i being∇ where V j = ∫ Ω S j dv is the nodal volume at node j. In Equation (27), the integral ∫ Ω S j ∇S i dv is calculated as a convolution of the gradient of the shape function (∇S i ) and the shape function of neighbor nodes j (S j ). Therefore, it can be solved analytically, as the background grid remains unchanged during the calculation. Figure 4 shows the gradient of the shape function in a 1D case. On the right-hand side in Equation (26), the first term represents the mapping of material points to particle-based nodes (defined as particle-based mapping) and the second term represents the mapping of material points to node-based nodes (defined as node-based mapping). Both mappings are weighted by α, a weighting function. The gradient of the shape function is continuous only if α is continuous and α = 0 at cell boundaries. The choice of α can retrieve the gradient of the shape function of other versions of MPM as follows: The Appendix gives more details related to function α. In this paper, the dual domain technique is used to switch the shape function gradient between MPM, cpGIMP, and DDMP by the choice of weighted function α in Equation (28).

Null-space instability in the gradient mapping matrix ∇S
Considering a mapping from a material point stress p to a grid node force i, the mapping can be written as follows: For a global mapping from all N p material point's stress to all N n grid node's internal force, the global mapping can be written as follows: where N n and N p are the total number of the nodes and particles, respectively. Equation (30) can be written in the matrix form where F int is the vector of all nodal internal forces and ( p V p ) is the product of the material point stress tensor and the corresponding material point volume, for all material points. The gradient mapping matrix ∇S depends on the value of the shape function gradient ∇S i (x p ), with the size of the matrix ∇S being N n x N p in 1D case. In reverse, the velocity gradient vector at material points is computed by the transpose of ∇S as follows: where ∇v p is the velocity gradient vector of all material points and v i is the nodal velocity vectors of all grid nodes. The null-space of the gradient mapping matrix ∇S, written as null(∇S), is the set of all nonzero independent vectors null p , which satisfy that ∇S · null p = 0, ie, The presence of the null-space may lead to infinite solutions. Hence, the mapping in Equation (31) may result in a nonunique solution. In other words, there possibly exists a set of stress vectors null p at material points, which leads to zero internal forces after mapping with Equation (31). Because the balance equations are solved at the grid, the null p in the null-space of ∇S cannot be detected on the background grid and represents the null-space noise. This noise can accumulate and leads to instability in the solution.
The rank of the gradient mapping matrix ∇S, written as rank(∇S), is the set of linearly independent column vector spaces. Equation (31) and Equation (32) give the gradient mapping based on the ∇S and its transpose ∇S T , which transfer the data between material points and grid nodes and vice versa. Therefore, the following four fundamental subspaces are considered: • the row space, rank(∇S T ), a subspace of R N n ; • the column space, rank(∇S), a subspace of R N p ; • the right null-space, null(∇S), a subspace of R N p ; • the left null-space, null(∇S T ), a subspace of R N n .
Since the column rank is equal to the row rank, 34 we have rank(∇S) = rank (∇S T ). According to the rank-nullity theorem, 34 the sum of column space and null-space is equal to the column numbers of the matrix as follows: In case there are more grid nodes than material points, the linear equations may lead to rank-deficiency problem. A similar problem may be encountered in the finite element method using a low order of quadrature rules (see rank-deficiency, page 239, in the work of Hughes 35 ). This issue can be mitigated in finite element method (and, in theory, in MPM) by using a high-order interpolation. Commonly, there are more material points than grid nodes, therefore, we only consider that N p > N n in the analysis. Considering N p > N n , as can be seen from Equation (34) and Equation (35), when ∇S is full rank, rank(∇S) is maximum and is equal to the number of grid nodes N n . In contrast, if ∇S is rank-deficient, that is, rank(∇S) = r < N n , the left null-space exists and null(∇S T ) = N n − r.
In the original MPM, the gradient mapping is rank-deficient 24 when there are more material points than nodes N p > N n . In such case, ∇S could have rank(∇S) = N n (ie, full rank) by using higher-order shape function (for example cubic B-spline) so that the left null-space is equal to zero null(∇S T ) = 0. For example, MPM, DDMP, and GIMP contain rank-deficient mappings, whereas cubic B-spline MPM 24 is full-rank mapping. However, right null-space still exists null(∇S) > 0 despite the choice of the shape function even for cubic B-spline MPM because it depends only on the difference between the number of material points and the number of nodes rather than on the order of the interpolating function. As such, a higher number of material points than the number of grid nodes in a cell leads to a right null-space and associated instability. That may be the reason why using one material point per cell could give less errors than a higher number of material points per cell 8 because N p < N n . As the null-space instability exists when N p > N n , one can enhance the stability of the solution with an application of a null-space filter, which removes the null-space error. The next section presents how to detect the null-space components of velocity gradients.

Detect null-space components of velocity gradients
To detect the null-space component, the transpose of the gradient mapping matrix needs to be decomposed in the QR form first, ie, where R is the upper triangular matrix with r ij are nonzero if ∇S is full rank. If ∇S is rank-deficient, the zero rows will develop from the bottom of matrix R and contribute to the left null-space (see Equation (36)). Q is the unitary matrix, in which the columns q i are orthonormal vectors that satisfy q i T q j = ij and span in the space R N p . Figure 5 shows the QR decomposition for the gradient mapping matrix. The unitary matrix Q is divided into the column space Q * N n x r and the null-space. The null-space includes the left null-space, null(∇S T ) = N n − r, and the right null-space null(∇S) = N p − N n . The velocity gradient vector of all the material points ∇v p ∈ R N p can be expressed as a linear combination of the orthonormal vector q i with the coefficients a i , ie, ∇v p = a 1 q 1 + · · · + a r q r ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟ rank(∇S)=r In Equation (37), the solution of the velocity gradient vector of the material points consists of a linear combination of column space and null-space. The left null-space, null(∇S T ), represents the loss of information due to rank-deficient mapping, whereas the remaining components, null(∇S) − null(∇S T ) = N p − N n , represent the loss of information due to the difference between the number of material points and the number of nodes. The nonnull-space component, after removing the null-space errors, is

Null-space filter algorithm for MPM
This section presents the null-space filter algorithm in greater detail. The null-space filter algorithm requires two steps. First, the algorithm computes the orthonormal basis of column space {q 1 , q 2 , … , q r }, as described in section A. Afterwards, it projects the velocity gradients to the orthonormal basis of the column space, as described in section B.

Compute the orthonormal basis of column space {q 1 , q 2 , … , q r }
Although Equation (36) is written for the global case, ie, for ∇S N n x N p , where N p and N n are the total number of material points and nodes, respectively, in practice, the null-space filter can be done locally with the use of the local QR algorithm gradient mapping for each cell is ∇S c . The size of ∇S c is n i by n p , with n i being the number of the nodes of the cell and n p the number of the material points inside this cell. Because the grid is not deformed during the calculation, the value of the gradient of the linear basis function is constant within the grid cell and independent from the material point coordinates in one dimension. Therefore, the rank of the local ∇S c can be computed analytically. As such, we can optimize the computational efficiency, with given r = rank (∇S c ), by modified householder algorithm (see Algorithm 1). Instead of computing all the n p columns of matrix Q and n i columns of matrix R, the modified householder only computes the column space of Q up to order of r. Therefore, the modified householder algorithm reduces the computational cost compared to the SVD 23 or QR decomposition 24 (see Figure 6 for numerical example 6.1). In a 3D case, with a structured grid, the rank r is computed, and then Algorithm 1 calculates three sets of orthonormal basis that are calculated separately for each dimension, leading to Q x c , Q y c , and Q z c . Note that, in this case, the values of the gradient ∇S c depend on the position of the material points.

Project the velocity gradients to the orthonormal basis of column space
To implement the projection for local null-space filter, we consider a cell that contains n i nodes and n p material points with the gradient mapping matrix ∇S c . Then, the velocity gradient is computed by the particle-based mapping, ie, Afterwards, the nonnull-space components of velocity gradients, ∇v * p , could be extracted by a projection with the orthonormal basis of column spaces {q 1 , q 2 , … , q r } 34 by In the 3D case, with a structured grid, the projection of the velocity gradient ∇v is done for each dimension as follows: where Q x c , Q x c , Q x c are the sets of orthonormal basis in each dimension and ∇v * p = [ ∇v * ,x p ; ∇v * ,y is the filtered value of the velocity gradient. Finally, the velocity gradient is weighted between ∇v * p and the velocity gradient using the node-based mapping∇v p , based on Equation (26), as The selection of α in Equation (49) leads to replication of the shape function gradient of different MPM formulations. Setting (α = 1) recovers the original formulation of MPM, GIMP is obtained with (α = α 2 ), and DDMP with (α = α 1 ) with the choice of α described in the Appendix.

DOUBLE-POINT MPM FORMULATION WITH NULL-SPACE FILTER
In this paper, the two-phase formulation bases on a double-point formulation of MPM, 24,36 in which each phase is represented by a set of two material points, one for the liquid phase and one for the solid phase. The information from the material points is projected to the grid separately for the solid and liquid phases (see Figure 7). The algorithms mapping the data from the material points to the grid and back are the same for both liquid and solid phases. The solid and liquid phases are interacting with each other via dragging forces. The governing and time discretization equations are presented in detail by Bandara and Soga 36 and Tran and Sołowski 37 and summarized in Figure 8. Note that the gradient mapping matrix ∇S is used to compute the nodal internal forces and the material point's strain. Therefore, to implement the null-space filter into the two-phase double-point formulation, the calculation of the nodal internal forces and the material point's strain is modified using the dual domain technique and the null-space filter and presented in the next section.

Calculation of nodal internal forces
The vectors of the internal forces for the solid and liquid phase f int,t αi (with α = s, w) are calculated using the dual domain technique, which, when used with specific values of coefficients, can replicate DDMP, GIMP, and the original MPM formulation, ie, where t sp and p t wp = are the effective stress and pore water pressure tensor, V t sp and V t wp are the solid and liquid volume.
is a weighted function evaluated at the material point x t αp , and the choice of weighted function is described in the Appendix to switch different version of MPM including MPM, GIMP, and DDMP. The first term in the right-hand side of Equation (43) and Equation (44) is the MPM particle-based mapping for the node i and the second term, in the right-hand side of Equation (44), corresponds to the node-based mapping from the neighbor nodes j. The nodal effective stresses for the solid phase t sj or nodal pore water pressures for the liquid phase p t wj at the neighbor nodes j are calculated as follows:

Calculation of the material points strain with null-space filter
The strain increment vectors Δε k+1 αp,1 using MPM particle-based mapping at the solid and liquid material points (x k αp ) are calculated as follows: The strain increment is filtered with null-space filter, using Equation (40). The filtered components using MPM particle-based mapping Δε * k+1 αp,1 is The strain increment using node-based mapping Δ k+1 αp,2 is where nodal strain increments Δ t αj calculated as follows: The final strain increment vectors Δ k+1 αp is calculated with the weighted function α(x), ie,

Square wave in an elastic bar
The first numerical example is a simulation of an elastic bar under compression. Figure 9 presents the boundary and loading condition of the numerical model. The length of the bar L = 1 m, the Young's modulus E = 10e 6 Pa, and the density ρ = 1000 kg/m 3 . The elastic wave velocity is calculated as c = √ E∕ρ = 100 m/s. The compression force f is applied to the material points at the right boundary to replicate the square stress wave in the elastic bar as follows:  Figure 10 shows that MPM, using PIC scheme (Equation (4)), can filter the velocity oscillations, hence, filter the stress/strain oscillations, whereas MPM, GIMP, and DDMP show stress oscillations (see Figure 11, Figure 12, and Figure 13). Similar oscillations were reported in the literature. 17,38 The generalized-α scheme can remove the oscillations for all methods including MPM, GIMP, and DDMP. In this example, the null-space filter does not change the results significantly because the errors are mainly stemming from the spurious noises from the high-frequency mode. To check the energy conservation, the total energy E is the sum of the strain U and the kinetic energy K are computed in each time step as follows: Figure 14 shows the total energy after releasing the compression forces at 0.005 s. The oscillations in MPM using FLIP scheme induces the increase of the total energy, whereas both MPM using PIC scheme and generalized-α scheme demonstrate approximately constant total, energy after releasing the compression force.

Colliding disks
The second example is the colliding impact of two elastic disks, replicated the simulation by Sulsky et al 1 with given geometry and parameters in Figure 15. The background grid is discretized with the cell size of 0.05 m in plane strain condition and there are total 896 material points and no contact law is applied. The time step is 0.005 s to be equal to 10%   Figure 16, Figure 17, and Figure 18 show the visualization of the disk impact with the velocity colorbar at initial condition, 1.05 s and 2.5 s. Figure 19, Figure 20, and Figure 21 present the evolution of the kinetic, strain, and total energy for the cpGIMP using PIC scheme, FLIP scheme, and generalized-α scheme, respectively. As can be seen,  cpGIMP using PIC update scheme shows a strong energy dissipation although it can filter the velocity oscillations in the previous example. In contrast, cpGIMP (with the originally used FLIP update scheme) improves significantly the energy conservation but leads to oscillations under compression impact in the previous example. Some other PIC variants (aPIC 19 and xPIC(m) 21 ) are developed to retain the filter property and improve the energy conservation but the total energy are lower compared with MPM (with the originally used FLIP update scheme). The cpGIMP using the α-generalized scheme can filter the velocity and stress noise due to high-frequency oscillations in the previous example and preserve similar total energy conservation with the cpGIMP without the excessive dissipation (see Figure 22).

Gaussian explosion
In the next example, we simulate a 2D wave propagation in an elastic domain with Gaussian-type explosion. 39,40 An explosive source was placed at [x c ; y c ] the center of a homogeneous, isotropic elastic square (see Figure 23), replicated in the simulation by an exponential force applying to material points with the position [x; y] as follows: where    (56) were applied to the material points as external forces. To avoid any reflection from the boundary, the wave propagation is computed only for the 8e −3 s with the time step being as Δt = 10e −4 s. Figure 24 and Figure 25 show the snapshot of the velocity wave for the original MPM and MPM enriched with the generalized-α scheme at t = 8e −3 s. When the wave propagates in the elastic domain, it leaves behind some velocity oscillations near the explosive source. This noise is reduced when the simulation is repeated with the generalized-α scheme for the MPM formulation (see Figure 26 at t = 8e −3 s). In such impact problems, the generalized-α scheme can reduce considerable the noises, whereas null-space filter again does not change significantly the results, similarly as in the previous example.

1D wave propagation
The following example is the wave propagation of an elastic bar due to the nonlinear predescribed displacement. The length of the bar is 2 m, Young's modulus is E = 10e 6 Pa, and the density ρ = 1000 kg/m 3 . Assuming the magnitude of displacement is A = 0.001m, the initial displacement and strain are generated as follows:  Figure 28 and Figure 29 present the numerical model using PIC scheme and FLIP scheme (MPM), respectively. The PIC shows a diffusive solution, as reported by Hammerquist and Nairn, 21 and both PIC and MPM show numerical noises in the middle of the bar due to the null-space errors. This noise can be filtered by using the null-space filter (see Figure 30 and Figure 31). In case there is only 1 material point per cell, leading to 200 material points in total (N p = 200), the null-space errors in material points will be eliminated as N p = rank (∇S). Figure 32 demonstrates that there is no oscillation for 1 particle per cell because there are no null-space errors.

Convergence rate study
In the example earlier, cell-crossing does not occur in the simulation. Therefore, the original MPM could solve the problem accurately. This section compares spatial convergence for MPM, GIMP, DDMP, and DDMP with null-space filter in two scenarios (no cell-crossing and cell-crossing as described in Figure 33). The same example as in the previous section, with  To compare between MPM, GIMP, and DDMP, the cell length of 0.004 m is selected as the reference case. It corresponds to 500 cells and 1000 material points in the entire 2-m bar. Cell-crossing causes significant errors in the right boundary of the bar when using MPM formulation ( Figure 34). These errors increase dramatically with the increase of the grid cell number, leading to cell-crossing instability. This instability can be reduced by using GIMP or DDMP formulation ( Figure 35 and Figure 36). However, in GIMP and DDMP, the strain still slightly oscillates in the middle of the bar and this oscillation can be reduced using the null-space filter (Figure 37 and Figure 38).
To study the convergence rate of the proposed formulation, the spatial error is defined as follows: where ε exact (x, t) is the exact strain from Equation (63) and ε(x, t) is the strain of the numerical solution. N p is the total number of material points. Figure 39 and Figure 40 show the spatial convergence rate of strain at the t = 0.005 s for the "no cell-crossing" cases and the "cell-crossing" case, respectively. The MPM can converge in the "no cell-crossing" case but fail in the "cell-crossing" case. Although the cell-crossing errors is dominant, in both cases, the null-space filter can reduce the errors for GIMP and DDMP.

Two-phase gravity loading and consolidation
In the previous section, the results of the simulation of wave propagation have indicated that the null-space filter can help in removing the numerical oscillations. This section shows that, in a two-phase quasi-static problem under large deformations, the null-space errors can grow significantly in time and result in a nonconvergence.

Numerical parameters
The simulations of the gravity loading and large-strain consolidation use double-point MPM with DDMP interpolation. The calculation considers a 1D 1-m high fully saturated elastic soil column, discretized into 50 grid cells. Each cell contains 2 solid material points and 2 water material points. In total, 100 solid material points and 100 water material points are   distributed equally in the column. Soil is modeled as a linear elastic material with Young's modulus of 10 MPa, the solid density equal to 2143 kg/m 3 and the initial porosity and permeability equal to 0.3 and 10 −3 m/s, respectively. For the liquid phase, the water density is taken as 1000 kg/m 3 and the bulk modulus is equal to 2.2 GPa. The time step is set to be Δt = 10e −6 s, which satisfies the CFL condition and the permeability-dependent criterion. 41 Table 1 summarizes the parameters of the numerical model.

Numerical solutions of gravity loading
The first example simulates the application of the gravity loading. Initially, the pore pressure and the effective stress are equal to zero. Then, the gravity loading is applied immediately after the first time step. As the dragging force, proportional to gravity acceleration, acts as a natural damping to counter the dynamic effect to reach the steady state, the gravity does not need to increase gradually with time, even though no numerical damping is applied. The final time of the simulation is 2 seconds, which is high enough to reach the steady state solution. For small strain, with the gravitational acceleration of 10 m/s 2 , both DDMP (rank-deficient mapping) and B-spline (full rank mapping) replicate the analytical solution (see Figure 41 and Figure 42). To obtain large displacements, the gravity is increased to 1500 m/s 2 , leading to approximately 10% of engineering strain. Figure 43 and Figure 44 present the numerical results of DDMP and cubic B-spline MPM without the null-space filter. Clearly, the pore water pressure oscillates under the large deformation without the filter, whereas the null-space filter corrects the water pressure so it just increases linearly with depth (see Figure 45).

Analytical solutions of large strain oedometric consolidation
The second example is simulations of the large strain oedometric consolidation validated with an analytical solution given by Xie and Leo. 42 In Terzaghi small strain consolidation theory, the soil thickness H, the permeability k, and Young's modulus of soil E assume to be constant. These assumptions are not valid for large strain consolidation, and therefore Xie and Leo 42 considered the soil thickness H(t), permeability k t , and modulus of soil E t being time dependent for the derivation of the consolidation equations. The MPM naturally considers the change of soil thickness and updates the porosity n every time step, therefore, we update the permeability and Young's modulus from porosity as follows: where n o , k o , and E o are the initial porosity, permeability, and modulus, respectively. The governing equations for the excess pore water pressure is where c v is the coefficient of consolidation. The coefficient of consolidation can be written as follows: where k is the permeability and m vl is the compressibility in the large strain related to void ratio as The excess pore water pressure along the depth z = [0, H] can be written as follows: where q is the traction forces and the time factor is given by the following: The average degree of consolidation defined by strain (U s ) is given by the following: The settlement at the top surface (S t ) is given by the following: .
The final top surface settlement, when U s → 1, is equal to Considering the motion of the soil skeleton (v s ), the velocities of soil skeleton are

Numerical solutions of large strain oedometric consolidation
In the numerical model, the compression force on top of the layer is equal to 2 GPa, which is enough to create approximately 20% of engineering strain. Therefore, initial pore pressure is set to 2 GPa, leading to a zero effective stress (no gravity is considered). The top boundary layer is drained, whereas the bottom boundary layer is undrained. Therefore, the water material points can freely move outside the solid domain through the top of the solid column. Indeed, it was observed during the analysis that some water particles flew over the top soil surface. The transition from a 2-phase porous medium to a single liquid phase led to oscillations in the top boundary layer. However, as there was no gravitational loading, these water particles were not important in the analysis. As such, the water particles that moved from soil body to a free water zone were simply removed, which reduced pressure and velocity oscillations. In the further research, the interface of free water-porous media transition will be investigated in greater detail. The numerical results (see Figure  identical as the analytical solution (see Equation (74)). The dynamic motion of the soil phase is captured in the numerical model and compared with the analytical model in Figure 48. It may be that the constant time step is somewhat too large at the beginning of the analysis under a high value of velocity, leading to some velocity oscillation. However, overall, the results show an acceptable correlation between the numerical and analytical solutions for approximately 20% of engineering strain.

Influence of the interpolation in the null-space instability
Both DDMP and cubic B-spline MPM can replicate the small strain consolidation (see Figure 49 and Figure 50), but not for large strain consolidation. Figure 51, Figure 52, and Figure 53 give the comparison between the numerical model and the analytical solution for the average consolidation degree defined by strain of 0.5 for DDMP (rank-deficient mapping), MPM with cubic B-spline (full rank mapping), and DDMP with null-space filter, respectively. The result shows that null-space filter removes the pore water pressure instability and leads to a good agreement between the numerical and analytical solutions.

Influence of the grid cell size on the null-space instability
To investigate the influence of grid size, the consolidation problem is solved with three different grids, consisting of 20 cells, 50 cells, and 100 cells. Figure 54, Figure 51, and Figure 55 show the DDMP solution of the large strain consolidation with a different number of cells (20 cells, 50 cells, and 100 cells, respectively). The figures show that grid refinement reduces the magnitude of the oscillations but does not remove them fully. The null-space errors, however, can be mitigated significantly by using the null-space filter (see Figure 53, Figure 56, and Figure 57). Note that the slight discrepancy between numerical and analytical solution in Figure 57 may be influenced by the removal of the liquid particles when moving outside the solid domain.

Influence of the material points per cell in the null-space instability
Next, the problem is solved with 2, 3, and 4 material points per grid cell, with the problem discretized with the total of 20 cells. Figure 54, Figure 58, and Figure 59 show the DDMP solution of the large strain consolidation with a different number of material points per cell (2 material points per cell, 3 material points per cell, and 4 material points per cell, respectively). Increasing the number of material points per cell does not reduce the oscillations and generates more null modes, which is also observed in the PIC method. 3 The null-space errors, again, can be mitigated significantly by using the null-space filter (see Figure 56, Figure 60, and Figure 61). In calculations without the null-space filter, the null-space errors accumulate rapidly in time and lead to a nonconvergent solution (the simulation collapse when pore water pressure approaches infinity). During the large deformation consolidation, the grid remains undeformed, but the number of nodes involved in the calculation reduces in time (see Figure 62). Therefore, the number of the column spaces reduces and the number of the null-space increases, which may be the reasons for the rapid increase of the null-space errors.

CONCLUSION
This paper shows how to implement the explicit generalized-α integration scheme in the MPM. The generalized-α integration scheme suppresses the high-frequency spurious oscillations and controls the dissipation with just a single parameter. Numerical examples show that the scheme can reduce the oscillations in the impact and wave propagation problems for MPM formulation. Furthermore, this paper gives a null-space filter algorithm, which improves the simulation accuracy and stability. The algorithm is suitable for different interpolations, including original MPM, GIMP, or DDMP. The proposed algorithm is much less expensive, when compared with the original method, 23 which used the SVD or QR decomposition. Several simulations demonstrated the improvement in stability in calculations with the null-space filter, for both single phase and two-phase material.
Apart from the null-space related errors, MPM simulation results are also affected by the cell-crossing errors and quadrature errors. While the former occurs when a material point crosses to a new cell, the latter results from the first-order integration of the field quantities, 43 leading to a reduction of the convergence rate in the numerical solution. The null-space error, on the other hand, occurs due to a null-space existing due to the linear mapping between material point and grid nodes. It explains why, in the vibration example (no cell-crossing), null-space error does not exist when there is only 1 particle per cell. Moreover, in the quasi-static hydro-coupling problems, the null-space errors increase dramatically because the number of nodes involved in the calculation reduces in time, leading to the expansion of the null-space.

APPENDIX
The function α, for DDMP α 1 , can be written as follows: where n i is the number of the nodes of the cell, which stores the material point (x). S i (x) is the linear basis shape function of the node "i" and n dim = 1, 2, 3 is the number of the dimension of the problem. one dimension, the function α can also retrieve the gradient shape function of MPM with α = 1 or of cpGIMP with α 2 set as follows: , L i − l p < |Δx| ≤ L i 1 − |Δx| l p · (Δx−sign(Δx)(Li+lp)) (Δx−sign(Δx)(Li+|Δx|)) , L i < |Δx| ≤ L i + l p 1 , else, where the gradient of the shape function in cpGIMP depends on the interval length Δx = x − x i , cell spacing L i , and the current material point's domain 2l p , which evolves with the deformation of the material point, corresponding to cpGIMP. Figure A1 shows function α for both GIMP and DDMP.