Extension of B-spline Material Point Method for unstructured triangular grids using Powell-Sabin splines

TL;DR

This paper extends the B-spline Material Point Method to unstructured triangular grids by integrating Powell-Sabin splines, reducing grid-crossing errors and achieving higher-order convergence in simulations.

Contribution

It introduces a novel extension of BSMPM to unstructured meshes using Powell-Sabin splines, enhancing applicability and accuracy.

Findings

Significant reduction in grid-crossing errors.

Higher-order convergence demonstrated in numerical tests.

Effective application on unstructured triangular grids.

Abstract

The Material Point Method (MPM) is a numerical technique that combines a fixed Eulerian background grid and Lagrangian point masses to simulate materials which undergo large deformations. Within the original MPM, discontinuous gradients of the piecewise-linear basis functions lead to so-called `grid-crossing errors' when particles cross element boundaries. Previous research has shown that B-spline MPM (BSMPM) is a viable alternative not only to MPM, but also to more advanced versions of the method that are designed to reduce the grid-crossing errors. In contrast to many other MPM-related methods, BSMPM has been used exclusively on structured rectangular domains, considerably limiting its range of applicability. In this paper, we present an extension of BSMPM to unstructured triangulations. The proposed approach combines MPM with $C^1$-continuous high-order Powell-Sabin (PS) spline basis functions. Numerical results demonstrate the potential of these basis functions within MPM in terms of grid-crossing-error elimination and higher-order convergence.

Figures

Equations

$$\frac{\partial u_k}{\partial t}$$

$$\rho \frac{\partial v_k}{\partial t}$$

$$D_{kl}$$

$${\sigma }_{kl}$$

$$u_k(\mathbf{x},t)$$

$${\sigma }_{kl}(\mathbf{x},t)n_l$$

$${\int }_{\Omega }\varphi \rho a_k\mathrm{d}\Omega ={\int }_{\Omega }\varphi \frac{\partial {\sigma }_{kl}}{\partial x_l}\mathrm{d}\Omega +{\int }_{\Omega }\varphi \rho g_k\mathrm{d}\Omega ,$$

$$\begin{gathered} {\int }_{\Omega }\varphi \rho a_k\mathrm{d}\Omega ={\int }_{\partial {\Omega }_{\tau }}\varphi {\tau }_k\mathrm{d}S \\ -{\int }_{\Omega }\frac{\partial \varphi }{\partial x_l}{\sigma }_{kl}\mathrm{d}\Omega +{\int }_{\Omega }\varphi \rho g_k\mathrm{d}\Omega , \end{gathered}$$

$$a_k(\mathbf{x},t)=\sum _{j=1}^{n_{bf}}{\hat{a}}_{k,j}{\varphi }_j(\mathbf{x}),$$

$$\begin{gathered} {\int }_{\Omega }{\varphi }_i\rho \left(\sum _{j=1}^{n_{bf}}{\hat{a}}_{k,j}{\varphi }_j\right)\mathrm{d}\Omega ={\int }_{\partial {\Omega }_{\tau }}{\varphi }_i{\tau }_k\mathrm{d}S \\ -{\int }_{\Omega }\frac{\partial {\varphi }_i}{\partial x_l}{\sigma }_{kl}\mathrm{d}\Omega +{\int }_{\Omega }{\varphi }_i\rho g_k\mathrm{d}\Omega , \end{gathered}$$

$$\begin{gathered} \sum _{j=1}^{n_{bf}}\underset{M_{ij}}{\underbrace{\left({\int }_{\Omega }{\varphi }_i\rho {\varphi }_j\mathrm{d}\Omega \right)}}{\hat{a}}_{k,j}=\underset{F_{k,i}^{trac}}{\underbrace{{\int }_{\partial {\Omega }_{\tau }}{\varphi }_i{\tau }_k\mathrm{d}S}} \\ -\underset{F_{k,i}^{int}}{\underbrace{{\int }_{\Omega }\frac{\partial {\varphi }_i}{\partial x_l}{\sigma }_{kl}\mathrm{d}\Omega }}+\underset{F_{k,i}^{body}}{\underbrace{{\int }_{\Omega }{\varphi }_i\rho g_k\mathrm{d}\Omega }}, \end{gathered}$$

$$\Rightarrow \mathbf{M}{\hat{\mathbf{a}}}_k={\mathbf{F}}_k^{trac}-{\mathbf{F}}_k^{int}+{\mathbf{F}}_k^{body},$$

$$M_{ij}^t=$$

$$F_{k,i}^{trac,t}=$$

$$F_{k,i}^{int,t}=$$

$$F_{k,i}^{body,t}=$$

$${\mathbf{M}}^t{\hat{\mathbf{a}}}_k^t=\left({\mathbf{F}}_k^{trac,t}-{\mathbf{F}}_k^{int,t}+{\mathbf{F}}_k^{body,t}\right).$$

$$v_{k,p}^{t+\Delta t}=v_{k,p}^t+\Delta t{a}_k^t({\mathbf{x}}_p^t)=v_{k,p}^t+\Delta t\sum _{j=1}^{n_{bf}}{\hat{a}}_{k,j}^t{\varphi }_j({\mathbf{x}}_p^t).$$

$$v_k^{t+\Delta t}(\mathbf{x})=\sum _{j=1}^{n_{bf}}{\hat{v}}_{k,j}^{t+\Delta t}{\varphi }_j(\mathbf{x}),$$

$${\mathbf{M}}^t{\hat{\mathbf{v}}}_k^{t+\Delta t}={\mathbf{P}}_k^{t+\Delta t},$$

$$P_{k,i}^{t+\Delta t}=\sum _{p=1}^{n_p}\underset{m_p}{\underbrace{V_p^t{\rho }_p^t}}{v}_{k,p}^{t+\Delta t}{\varphi }_i({\mathbf{x}}_p).$$

$$D_{kl,p}^{t+\Delta t}=$$

$$J_p^{t+\Delta t}=$$

$${\epsilon }_{kl,p}^{t+\Delta t}=\frac{1}{2}\sum _{j=1}^{n_{bf}}\left({\hat{v}}_{k,j}^{t+\Delta t}\frac{\partial {\varphi }_j}{\partial x_l}({\mathbf{x}}_p^t)+{\hat{v}}_{l,j}^{t+\Delta t}\frac{\partial {\varphi }_j}{\partial x_k}({\mathbf{x}}_p^t)\right).$$

$$E_{kl,p}^{t+\Delta t}=\frac{1}{2}\left(D_{kl,p}^{t+\Delta t}+D_{lk,p}^{t+\Delta t}\right)-{\delta }_{kl},$$

$${\sigma }_{kl,p}^{t+\Delta t}=\lambda {\delta }_{kl}\mathrm{tr}\left(E_{kl,p}^{t+\Delta t}\right)+2\mu E_{kl,p}^{t+\Delta t}.$$

$${\sigma }_{kl,p}^{t+\Delta t}=\frac{\lambda \ln (J_p^{t+\Delta t})}{J_p^{t+\Delta t}}{\delta }_{kl}+\frac{\mu }{J_p^{t+\Delta t}}\left(D_{km,p}^{t+\Delta t}{D}_{lm,p}^{t+\Delta t}-{\delta }_{kl}\right).$$

$$V_p^{t+\Delta t}=$$

$${\rho }_p^{t+\Delta t}=$$

$$x_{k,p}^{t+\Delta t}=$$

$$u_{k,p}^{t+\Delta t}=$$

$$u_x$$

$$u_y$$

$$g_x$$

$$g_y$$

$$\lambda$$

Full text

1 Introduction

The Material Point Method (MPM) has proven to be successful in solving complex engineering problems that involve large deformations, multi-phase interactions and history-dependent material behaviour. Over the years, MPM has been applied to a wide range of applications, including modelling of failure phenomena in single- and multi-phase media andersen2010land ; Sulsky2004 ; Zabala2011 , crack growth Nairn2003 ; guo2017simulation and snow and ice dynamics Stomakhin2013 ; sulsky2007ice ; Gaume2018 . Within MPM, a continuum is discretised by defining a set of Lagrangian particles, called material points, which store all relevant material properties. An Eulerian background grid is adopted on which the equations of motion are solved in every time step. The solution on the background grid is then used to update all material-point properties such as displacement, velocity, and stress. In this way, MPM incorporates both Eulerian and Lagrangian descriptions. Similarly to other combined Eulerian-Lagrangian techniques, MPM attempts to avoid the numerical difficulties arising from nonlinear convective terms associated with an Eulerian problem formulation, while at the same time preventing grid distortion, typically encountered within mesh-based Lagrangian formulations (e.g., Donea2004 ; sulsky1994particle ).

Classically, MPM uses piecewise-linear Lagrange basis functions, also known as ‘tent’ functions. However, the gradients of these basis functions are discontinuous at element boundaries. This leads to so-called grid-crossing errors bardenhagen2004generalized when material points cross this discontinuity. Grid-crossing errors can significantly influence the quality of the numerical solution and may eventually lead to a lack of convergence (e.g., steffen2008analysis ). Different methods have been developed to mitigate the effect of grid-crossings. For example, the Generalised Interpolation Material Point (GIMP) bardenhagen2004generalized and Convected Particle Domain Interpolation (CPDI) Sadeghirad2011 methods eliminate grid-crossing inaccuracies by introducing an alternative particle representation. The GIMP method represents material points by particle-characteristic functions and reduces to standard MPM, when the Dirac delta function centered at the material-point position is selected as the characteristic function. For multivariate cases, a number of versions of the GIMP method are available such as the uniform GIMP (uGIMP) and contiguous-particle GIMP (cpGIMP) methods. The CPDI method extends GIMP in order to accurately capture shear distortion. Much research has been performed to further improve the accuracy of the CPDI approach and increase its range of applicability Sadeghirad2013 ; Nguyen2017 ; Leavy2019 . On the other hand, the Dual Domain Material Point (DDMP) method zhang2011material preserves the original point-mass representation of the material points, but adjusts the gradients of the basis functions to avoid grid-crossing errors. The DDMP method replaces the gradients of the piecewise-linear Lagrange basis functions in standard MPM by smoother ones. The DDMP method with sub-points Dhakal2016 proposes an alternative manner for numerical integration within the DDMP algorithm.

The B-spline Material Point Method (BSMPM) steffen2008analysis ; Steffen2010 solves the problem of grid-crossing errors completely by replacing piecewise-linear Lagrange basis functions with higher-order B-spline basis functions. B-spline and piecewise-linear Lagrange basis functions possess many common properties. For instance, they both satisfy the partition of unity, are non-negative and have a compact support. The main advantage of higher-order B-spline basis functions over piecewise-linear Lagrange basis function is, however, that they have at least $C^{0}$-continuous gradients which preclude grid-crossing errors from the outset. Moreover, spline basis functions are known to provide higher accuracy per degree of freedom as compared to $C^{0}$-finite elements Hughes2014 . On structured rectangular grids, adopting B-spline basis functions within MPM not only eliminates grid-crossing errors but also yields higher-order spatial convergence steffen2008analysis ; Steffen2010 ; tielen2017high ; Gan2018mpm ; Bazilevs2017 . Previous research also demonstrates that BSMPM is a viable alternative to the GIMP, CPDI and DDMP methods steffen2008examination ; Motlagh2017 ; Gan2018mpm ; wobbestaylor . While the CPDI and DDMP methods can be used on unstructured grids zhang2011material ; Nguyen2017 ; Leavy2019 , to the best of our knowledge, BSMPM for unstructured grids does not yet exist. This implies that its applicability to real-world problems is limited compared to the CPDI and DDMP methods.

In this paper, we propose an extension of BSMPM to unstructured triangulations to combine the benefits of B-splines with the geometric flexibility of triangular grids. The proposed method employs quadratic Powell-Sabin (PS) splines powell1977piecewise . These splines are piecewise higher-order polynomials defined on a particular refinement of any given triangulation and are typically used in computer aided geometric design and approximation theory speleers2012isogeometric ; dierckx1992algorithms ; Manni2007 ; Sablonniere1987 . PS splines are $C^{1}$-continuous and hence overcome the grid-crossing issue within MPM by design. We would like to remark that although this paper focuses on PS splines, other options such as refinable $C^{1}$ splines Nguyen2016 can be used to extend MPM to unstructured triangular grids. The proposed PS-MPM approach is analysed based on several benchmark problems.

The paper is organised as follows. In Section 2, the governing equations are provided together with the MPM solution strategy. In Section 3, the construction of PS-spline basis functions and their application within MPM is described. In Section 4, the obtained numerical results are presented. In Section 5, conclusions and recommendation are given.

2 Material Point Method

We will summarise the MPM as introduced by Sulsky et al. sulsky1994particle to keep this work self-contained. First, the governing equations are presented, afterwards the MPM strategy to solve these equations is introduced.

2.1 Governing equations

The deformation of a continuum is modelled using the conservation of linear momentum and a material model. It should be noted that MPM can be implemented with a variety of material models that either use the rate of deformation (i.e., the symmetric part of the velocity gradient) or the deformation gradient. However, for this study it is sufficient to consider the simple linear-elastic and neo-Hookean models that are based on the deformation gradient. Using the Einstein summation convention, the system of equations in a Lagrangian frame of reference for each direction $x_{k}$ is given by

[TABLE]

where $u_{k}$ is the displacement, $t$ is the time, $v_{k}$ is the velocity, $\rho$ is the density, $g_{k}$ is the body force, $D_{kl}$ is the deformation gradient, $\delta_{kl}$ is the Kronecker delta, $\sigma_{kl}$ is the stress tensor, $x^{0}$ is the position in the reference configuration, $\lambda$ is the Lamé constant, $E_{kl}$ is the strain tensor defined as $\frac{1}{2}\left(D_{kl}+D_{lk}\right)-\delta_{kl}$, $J$ is the determinant of the deformation gradient and $\mu$ is the shear modulus. Equation 2 and Equation 4 describe the conservation of linear momentum and the material model, respectively.

Initial conditions are required for the displacement, velocity and stress tensor. The boundary of the domain $\Omega$ can be divided into a part with a Dirichlet boundary condition for the displacement and a part with a Neumann boundary condition for the traction:

[TABLE]

where $\boldsymbol{x}=[\>x_{1}\>\>x_{2}\>]^{T}$ and $\boldsymbol{n}$ is the unit vector normal to the boundary $\partial\Omega$ and pointing outwards. We remark that the domain $\Omega$, its boundary $\partial\Omega=\partial\Omega_{u}\cup\partial\Omega_{\tau}$ and the normal unit vector $\boldsymbol{n}$ may all depend on time.

For solving the equations of motion in the Material Point Method, the conservation of linear momentum in Equation 2 is needed in its weak form. For the weak form, Equation 2 is first multiplied by a continuous test function $\phi$ that vanishes on $\partial\Omega_{u}$, and is subsequently integrated over the domain $\Omega$:

[TABLE]

in which $a_{k}=\frac{\partial v_{k}}{\partial t}$ is the acceleration. After applying integration by parts, the Gauss integration theorem and splitting the boundary into Dirichlet and Neumann part, the weak formulation becomes

[TABLE]

whereby the contribution on $\partial\Omega_{u}$ equals zero, as the test function vanishes on this part of the boundary.

2.2 Discretised equations

Equation 8 can be solved using a finite element approach by defining $n_{bf}$ basis functions $\phi_{i}$ $(i=1,\dots,n_{bf})$. The acceleration field $a_{k}$ is then discretised as a linear combination of these basis functions:

[TABLE]

in which $\hat{a}_{k,j}$ is the time-dependent $j^{th}$ acceleration coefficient corresponding to basis function $\phi_{j}$. Substituting Equation 9 into Equation 8 and expanding the test function in terms of the basis functions $\phi_{i}$ leads to

[TABLE]

which holds for $i=1,\dots,n_{bf}$. By exchanging summation and integration, this can be rewritten in matrix-vector form as follows:

[TABLE]

[TABLE]

where $\boldsymbol{M}$ denotes the mass matrix, $\boldsymbol{\hat{a}}_{k}$ the coefficient vector for the acceleration, while $\boldsymbol{F}_{k}^{trac}$, $\boldsymbol{F}_{k}^{int}$ and $\boldsymbol{F}_{k}^{body}$ denote respectively the traction, internal and body force vector in the $x_{k}$-direction. When density, stress and body force are known at time $t$, the coefficient vector $\boldsymbol{\hat{a}}_{k}$ is found from Equation 12. The properties of the continuum at time $t+\Delta t$ can then be determined using the acceleration field at time $t$.

2.3 Solution procedure

Within MPM, the continuum is discretised by a set of particles that store all its physical properties. At each time step, the particle information is projected onto a background grid, on which the momentum equation is solved. Particle properties and positions are updated according to this solution, as illustrated in Figure 1.

To solve Equation 12 at every time step, the integrals in Equation 11 have to be evaluated. Since the material properties are only known at the particle positions, these positions are used as integration points for the numerical quadrature, and the particle volumes $V$ are adopted as integration weights. Assume that the total number of particles is equal to $n_{p}$. Here, subscript $p$ corresponds to particle properties. A superscript $t$ is assigned to particle properties that change over time. The mass matrix and force vectors are then defined by

[TABLE]

in which $m_{p}$ denotes the particle mass, which is set to remain constant over time, guaranteeing conservation of the total mass. The coefficient vector $\boldsymbol{\hat{a}}_{k}^{t}$ at time $t$ for the $x_{k}$-direction is determined by solving

[TABLE]

Unless otherwise stated, Equation 17 is solved adopting the consistent mass matrix. Next, the reconstructed acceleration field is used to update the particle velocity:

[TABLE]

The semi-implicit Euler-Cromer scheme cromer1981stable is adopted to update the remaining particle properties. First, the velocity field $v_{k}^{t+\Delta t}$ is discretised as a linear combination of the same basis functions, similar to the acceleration field (Eq. 9):

[TABLE]

in which $\boldsymbol{\hat{v}}_{k}^{t+\Delta t}$ is the velocity coefficient vector for the velocity field at time $t+\Delta t$ in the $x_{k}$-direction. The coefficient vector $\boldsymbol{\hat{v}}_{k}^{t+\Delta t}$ is then determined from a density weighted $L_{2}$-projection onto the basis functions $\phi_{i}$. This results in the following system of equations:

[TABLE]

where $\boldsymbol{M}^{t}$ is the same mass matrix as defined in Equation 13. $\boldsymbol{P}_{k}^{t+\Delta t}$ denotes the momentum vector with the coefficients given by

[TABLE]

Particle properties are subsequently updated in correspondence with Equations 1-4. First, the deformation gradient and its determinant are obtained:

[TABLE]

Here, $\varepsilon_{kl,p}^{t+\Delta t}$ denotes the symmetric part of the velocity gradient:

[TABLE]

For linear-elastic materials, the particle stresses are computed as follows:

[TABLE]

[TABLE]

For neo-Hookean materials, the stresses are obtained from

[TABLE]

The determinant of the deformation gradient is used to update the volume of each particle. Based on this volume, the density is updated in such way that the mass of each particle remains constant:

[TABLE]

Finally, particle positions and displacements are updated from the velocity field:

[TABLE]

The described MPM can be numerically implemented by performing the steps in Equations 13-31 in each time step in the shown order. Note that all steps can be applied with a variety of basis functions, without any essential difference in the algorithm. In this paper, we investigate the use of Powell-Sabin (PS) spline basis functions, which are described in the next section.

3 Powell-Sabin spline basis functions

In this section, the essential properties of Powell-Sabin (PS) spline basis functions are shortly summarised, for the full construction process and the proofs of the properties we refer to the publication by Dierckx dierckx1997calculating ; dierckx1992algorithms . The considered PS-splines are piecewise quadratic splines with global $C^{1}$-continuity. They are defined on arbitrary triangulations, have local support and possess the properties of non-negativity and partition of unity. Several examples of PS-splines are shown in Figure 3.

To construct PS-splines on an arbitrary triangulation, a Powell-Sabin refinement is required, dividing each of the original main triangles into six sub-triangles as follows (see also Figure 2).

For each triangle $t_{j}$, choose an interior point $Z_{j}$ (e.g., the incenter), such that if triangles $t_{j}$ and $t_{m}$ have a common edge, the line between $Z_{j}$ and $Z_{m}$ intersects the edge. The intersection point is called $Z_{jm}$. 2. 2.

Connect each $Z_{j}$ to the vertices of its triangle $t_{j}$ with straight lines. 3. 3.

Connect each $Z_{j}$ to all edge points $Z_{jm}$ with straight lines. In case $t_{j}$ is a boundary element, connect $Z_{j}$ to an arbitrary point on each boundary edge (e.g., the edge middle).

Each PS-spline is associated with a main vertex $V_{i}$, and will be locally defined on its molecule $\Omega_{i}$, which is defined as the area consisting of the main triangles around $V_{i}$, as shown in Figure 2. Next, for each $V_{i}$, a set of PS-points is defined as the union of $V_{i}$ and the sub-triangle edge midpoints directly around it (see Figure 3). A control triangle for $V_{i}$ is then defined as a triangle that contains all its PS-points, and preferably has a small surface dierckx1997calculating . The set of control triangles uniquely defines the set of PS-splines, in which three PS-splines are associated with each control triangle or main vertex. The construction process from the control triangles ensures the properties of non-negativity and partition of unity (the details are omitted here, but can be found in dierckx1997calculating ; dierckx1992algorithms ). Figure 3 shows a PS-refinement, a main vertex and its associated molecule, a possible control triangle and its three associated PS-spline basis functions. Note that for PS-splines, each main vertex is associated with three PS-spline basis functions, in contrast to a piecewise-linear Lagrange basis, for which each main vertex is associated with only one basis function.

4 Numerical Results

To validate the proposed PS-MPM, several benchmarks involving large deformations are considered. The first benchmark describes a thin vibrating bar, where the displacement is caused by an initial velocity. For this benchmark, PS-MPM on a structured triangular grid is compared with a reference solution. After that, a vibrating plate with known analytical solution is considered. The spatial convergence of PS-MPM on an unstructured grid is determined for this benchmark. The third benchmark describes a soil-column under self-weight. For this benchmark, we propose the strategy of partial lumping to mitigate oscillations in the numerical solution that occur when the mass matrix is lumped.

4.1 Vibrating bar

In this section, a thin linear-elastic vibrating bar is considered with both ends fixed. A UL-FEM ten2007large solution on a very fine grid serves as reference. The grid used for PS-MPM and the initial particle positions are shown in Figure 4.

The bar is modelled with density $\rho=25\,$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, Young’s modulus $E=50\,$\mathrm{P}\mathrm{a}, Poisson ratio $\nu=0$, length $L=1\,$\mathrm{m}, width $W=2\,$\mathrm{m} and initial maximum velocity $v_{0}=0.1\,$\mathrm{m}\mathrm{/}\mathrm{s}$$. The chosen parameters result in a maximum normal strain in the $x$-direction of approximately $7\%$. At the left and right boundary, homogeneous Dirichlet boundary conditions are imposed for both $x$- and $y$-displacement, whereas at the top and bottom boundary, a homogeneous Dirichlet boundary condition is imposed only for the $y$-displacement, and a free-slip boundary condition for the $x$-direction. The initial displacement equals zero, but an initial $x$-velocity profile is set to $v_{x}(x^{0},y^{0},t)=0.1\sin(x^{0}\pi/L)$. The initial $y$-velocity is equal to zero.

The time step size for the simulation is $\Delta t=5\cdot 10^{-3}$ s. The corresponding Courant number is defined as $C=\tfrac{\Delta t}{h}\sqrt{E/\rho}$, in which $\sqrt{E/\rho}$ is the characteristic wave speed of the material and $h$ is the typical element length. Due to the ambiguity of $h$ for a PS-refinement on an unstructured triangular grid, we estimate the average edge length of the sub-triangles by $h\approx 0.025$ m. In this case, the Courant number is $C\approx 0.57<1$, satisfying the CFL condition courant1928partiellen .

Figure 5 shows the displacement in the middle of the bar over time and the stress profile through the entire bar at the end of the simulation. Although relatively coarse PS-MPM grid and particle configurations are adopted, the method yields accurate results and shows a smooth stress profile.

4.2 Vibrating plate undergoing axis-aligned displacement

In this section, a neo-Hookean vibrating plate on the unit square ($L\times W$, with $L=W=1$ m) is considered. The plate undergoes axis-aligned displacement. An analytic solution for this problem is constructed using the method of manufactured solutions (MMS): the analytic solution is assumed a priori, and the corresponding body force is calculated accordingly. This benchmark has been adopted from sadeghirad2011convected . The analytical solution for the displacement in terms of the reference configuration is given by

[TABLE]

Here, $\rho_{0}=10^{3}\,$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $u_{0}=0.05\,$\mathrm{m} and $E=10^{7}\,$\mathrm{P}\mathrm{a}$$. The corresponding body forces sadeghirad2011convected are

[TABLE]

in which the Lamé constant $\lambda$, the shear modulus $\mu$, and the normal components of the deformation gradient $D_{xx}$ and $D_{yy}$ are defined as

[TABLE]

Here, $\nu=0.3$ and all solutions are again given with respect to the reference configuration.

This benchmark was simulated with standard MPM and PS-MPM, using an unstructured triangular grid with material points initialised uniformly over the domain, as shown in Figure 6. A time step size of $\Delta t=2.25\cdot 10^{-4}$ s was chosen, which results in a Courant number of approximately $0.36$. For the adopted parameters, the imposed solution in Equation 32-33 has period of $T=0.02$ s.

4.2.1 Grid-crossing error

First, it is shown that the grid-crossing error typical for standard MPM does not occur in PS-MPM, by comparing the normal horizontal stress resulting from these methods (see Figure 7). The configurations for PS-MPM and standard MPM contain the same number of particles ($5120$, 4 times as many as in Figure 6) and a comparable number of basis functions, $289$ for standard MPM and $243$ basis functions for PS-MPM. As expected, the stress field in standard MPM suffers severely from grid-crossing, whereas with PS-MPM a smooth stress field is observed. The same conclusion can also be drawn when investigating individual particle stresses over time, as is shown in Figure 8. For this figure, the particle positioned at $(x^{0},y^{0})\approx(0.25$ m,$0.47$ m$)$ has been traced throughout the simulation. The figure depicts the stress and trajectory of the particle. For standard MPM, it is observed that the stress profile is highly oscillatory, resulting in a disturbed trajectory. PS-MPM shows no oscillations, but only a small offset, which was found to disappear upon further refinement of the grid.

4.2.2 Spatial convergence

PS-MPM with quadratic basis functions is expected to show third-order spatial convergence. To determine the spatial convergence, the time averaged root-mean-squared (RMS) error is used:

[TABLE]

where $n_{t}$ denotes the total number of time steps of the simulation. Under the assumption that the spatial error is much larger than the error produced by the numerical integration and time-stepping scheme, the spatial convergence of standard MPM and PS-MPM is determined by varying the typical element length $h$. This length which was defined as the average edge length for standard MPM and the average sub-triangle edge length for PS-MPM. It has been observed that the time-integration error is indeed sufficiently small, but the number of particles required for an adequately accurate integration increases rapidly as $h$ decreases.

Figure 9 shows the spatial convergence of both standard MPM and PS-MPM for different numbers of particles per element on an unstructured grid. Provided that a sufficient number of particles is adopted, standard MPM shows second-order spatial convergence, whereas PS-MPM converges with third order. Furthermore, the RMS-error is smaller associated with PS-MPM for all considered configurations. To achieve the optimal convergence order for both standard MPM and PS-MPM, the number of integration particles required increases rapidly, as the integration error otherwise dominates the total error, as shown in Figure 9. Note that the convergence rate in the number of particles for both standard MPM and PS-MPM is measured to be of first order. However, for a fixed typical element length $h$ and number of particles per element, the use of PS-MPM leads to a lower RMS-error, illustrating the higher accuracy per degree of freedom. The inaccurate integration in MPM due to the use of particles as integration points is known to limit the spatial convergence. Over the years, different measures have been proposed to decrease the quadrature error in MPM based on function reconstruction techniques like MLS gong2015improving , cubic splines tielen2017high and Taylor least squares wobbestaylor . The combination of PS-MPM with function reconstruction techniques to obtain optimal convergence rates with a moderate number of particles is subject to future research to make PS-MPM even more suited for practical applications.

4.3 Column under self-weight

In the previous examples, we considered PS-MPM adopting a consistent mass matrix. However, lumping of the mass matrix is common practice in many applications of MPM, as it speeds up the simulation and increases the numerical stability. In case elements become (almost) empty, the consistent mass matrix in Equation 13 has very small entries and becomes ill-conditioned, leading to stability issues, which is a well-known phenomena in MPM sulsky1994particle . Using a lumped version of the mass matrix to solve for the acceleration and velocity fields, is known to overcome the ill-conditioning sulsky1994particle . However, we have observed that lumping within PS-MPM also causes oscillations as a side-effect, which is not the case in standard MPM.

To illustrate this phenomenon, a soil-column benchmark under self-weight is considered, as shown in Figure 10. The soil column in this exemplary benchmark is modelled as a linear-elastic material. It is fixed in all directions at the bottom, and completely free at the top. At the left and right boundary, we impose a free-slip boundary condition for movement in the $y$-direction, but the displacement in the $x$-direction is fixed to be zero. The column is compressed by gravity due to its self-weight. The soil column is modelled with density $\rho=1\cdot 10^{3}\,$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, Young’s modulus $E=1\cdot 10^{5}\,$\mathrm{P}\mathrm{a}, Poisson ratio $\nu=0$, gravitational acceleration $g=-9.81\,$\mathrm{m}\mathrm{/}\mathrm{s}^{2}, height $H=1\,$\mathrm{m} and width $W=0.1\,$\mathrm{m}$$. The maximum strain when adopting these parameters is approximately $18\%$.

Figure 10 illustrates the grid considered for this benchmark as well as the particle positions when the deformation is maximal.

Due to the fact that part of the grid becomes empty, the use of a consistent mass matrix leads to stability issues with unlumped PS-MPM, which makes a PS-MPM simulation impossible. Therefore, a lumped mass matrix was adopted for this benchmark to restore stability. The use of a lumped mass matrix within PS-MPM indeed restores the stability, but spatial oscillations occur in the stress profile, degrading the solution of the displacement and velocity as well, as shown in Figure 11 (left column).

4.3.1 Partial lumping to mitigate spatial oscillations

A solution to mitigate the oscillations when lumping the mass matrix in PS-MPM, is the application of partial lumping. Instead of applying full lumping on the mass matrix, partial lumping only lumps those rows responsible for the ill-conditioning. These rows typically correspond to the basis functions in the part of the grid where very few particles are left. For this benchmark, all rows corresponding to a basis function with at least one empty main-triangle in its molecule are lumped. Additional empty elements are added at the top of the column, to ensure that the top-most basis functions are always lumped. Figure 11 (right column) shows the displacement and velocity of a particle over time, as well as the stress profile at $t=2.5$ s obtained with partial lumped PS-MPM. Compared to the results obtained with lumped PS-MPM, shown in Figure 11 (left), the stress profile, the velocity and displacement over time significantly improve.

Hence, the use of partial lumping within PS-MPM combines the advantages of adopting a consistent and lumped mass matrix, while minimising their drawbacks, improving the overall performance of the method. Future research will focus on techniques to further mitigate or fully overcome the oscillations caused by lumping, in particular when PS-MPM is applied for practical problems.

5 Conclusion

In this paper, we presented an alternative for B-spline MPM suited for unstructured triangulations using piecewise quadratic $C^{1}$-continuous Powell-Sabin spline basis functions. The method combines the benefits of smooth, higher-order basis functions with the geometric flexibility of an unstructured triangular grid. PS-MPM yields a mathematically sound approach to eliminate grid-crossing errors, due to the smooth gradients of the basis functions. As a first validation, a vibrating bar was considered, for which the PS-MPM solution yields accurate results on a relatively coarse grid. Numerical simulations, obtained for a vibrating plate undergoing axis-aligned displacement, have shown higher-order convergence for the particle stresses and displacements.

The use of a lumped mass matrix to increase stability within PS-MPM, leads to spatial oscillations in the stress profile. A partial lumping strategy was proposed, to combine the advantages of adopting a consistent and lumped mass matrix, which successfully mitigates this issue. Investigation of alternative unstructured spline technologies, in particular, refinable $C^{1}$ splines on irregular quadrilateral grids Nguyen2016 is underway.

Conflict of interest

Pascal de Koster has received funding support from Dutch research institute Deltares.