A variational discrete element method for quasi-static and dynamic elasto-plasticity

The domain $\Omega$ is discretized with a mesh $\mathcal{T}_{h}$ of size $h$ made of polyhedra with planar facets in three space dimensions or polygons with straight edges in two space dimensions. We assume that $\Omega$ is itself a polyhedron or a polygon so that the mesh covers $\Omega$ exactly, and we also assume that the mesh is compatible with the partition of the boundary $\partial\Omega$ into the Dirichlet and Neumann parts.

Let $\mathcal{C}$ denote the set of mesh cells and $\mathcal{Z}^{\partial}$ the set of mesh vertices sitting on the boundary of $\Omega$. Vector-valued volumetric degrees of freedom (dofs) for a generic displacement field $v_{\mathcal{C}}:=(v_{c})_{c\in\mathcal{C}}\in\mathbb{R}^{d\#(\mathcal{C})}$ are placed at the barycentre of every mesh cell $c\in\mathcal{C}$, where $\#(S)$ denotes the cardinality of any set $S$. Additional vector-valued boundary degrees of freedom $v_{\mathcal{Z}^{\partial}}:=(v_{z})_{z\in\mathcal{Z}^{\partial}}\in\mathbb{R}^{d\#(\mathcal{Z}^{\partial})}$ for the displacement are added at every boundary vertex $z\in\mathcal{Z}^{\partial}$. The reason why we introduce boundary vertex dofs is motivated in Section 3.3. These dofs are also used to enforce the Dirichlet condition on $\partial\Omega_{D}$. We use the compact notation $v_{h}:=(v_{\mathcal{C}},v_{\mathcal{Z}^{\partial}})$ for the collection of all the cell dofs and all the boundary vertex dofs. Figure 1 illustrates the position of the displacement dofs. In addition a (trace-free) symmetric tensor-valued dof representing the internal plasticity variable $\eta_{p,c}\in\mathbb{R}^{d\times d}_{\mathrm{sym}}$ is attached to every mesh cell $c\in\mathcal{C}$, as well as a scalar dof $p_{c}$ representing the cumulated plastic deformation. We write $\eta_{p,\mathcal{C}}:=(\eta_{p,c})_{c\in\mathcal{C}}\in Q_{h}:=\left(\mathbb{R}^{d\times d}_{\mathrm{sym}}\right)^{\#(\mathcal{C})}$ and $p_{\mathcal{C}}:=(p_{c})_{c\in\mathcal{C}}\in P_{h}:=\mathbb{R}^{\#(\mathcal{C})}$.

Let $\mathcal{F}$ denote the set of mesh facets. We partition this set as $\mathcal{F}=\mathcal{F}^{i}\cup\mathcal{F}^{b}$, where $\mathcal{F}^{i}$ is the collection of the internal facets shared by two mesh cells and $\mathcal{F}^{b}$ is the collection of the boundary facets sitting on the boundary $\partial\Omega$ (such facets belong to the boundary of only one mesh cell). Using the cell and boundary-vertex dofs introduced above, we reconstruct a collection of displacements $v_{\mathcal{F}}:=(v_{F})_{F\in\mathcal{F}}\in\mathbb{R}^{d\#(\mathcal{F})}$ on the mesh facets. The facet reconstruction operator is denoted $\mathcal{R}$ and we write

The precise definition of the facet reconstruction operator is given in Section 3.3. Using the reconstructed facet displacements and a discrete Stokes formula, it is possible to devise a discrete $\mathbb{R}^{d\times d}$-valued piecewise-constant gradient field for the displacement that we write $G_{\mathcal{C}}(v_{\mathcal{F}}):=(G_{c}(v_{\mathcal{F}}))_{c\in\mathcal{C}}\in\mathbb{R}^{d^{2}\#(\mathcal{C})}$. Specifically we set in every mesh cell $c\in\mathcal{C}$,

where the summation is over the facets $F$ of $c$ and $n_{F,c}$ is the outward normal to $c$ on $F$. Note that for all $v_{h}\in V_{h}$, we have

since $\sum_{F\in\partial c}|F|n_{F,c}=0$. We define a constant linearized strain tensor in every mesh cell $c\in\mathcal{C}$ such that

Finally, we define two additional reconstructions. The first is a cellwise nonconforming $P^{1}$ reconstruction $\mathfrak{R}$ defined for all $c\in\mathcal{C}$ by

where $\mathbf{x}\in c$ and $\mathbf{x}_{c}$ is the barycentre of the cell $c$. The second is a $P^{1}$ reconstruction on boundary facets, written $\mathfrak{R}^{\partial}(v_{\mathcal{Z}^{\partial}})$, and computed using the vertex dofs of the boundary facets. In case of simplicial facets, it reduces to classical $P^{1}$ Lagrange interpolation. For non-simplicial facets, generalised barycentric coordibnates need to be used, see [5], for instance.

Let us set $V_{h}:=\mathbb{R}^{d\#(\mathcal{C})}\times\mathbb{R}^{d\#(\mathcal{Z}^{\partial})}$ and (recall that $\partial\Omega_{D}$ is a closed set)

(Note that the spaces $V_{hD}$ and $W_{hD}$ are actually time-dependent if the Dirichlet data is time-dependent.) The discrete stiffness bilinear form is parameterized by a member $\eta_{p,\mathcal{C}}\in Q_{h}$ and is such that, for all $(v_{h},\tilde{v}_{h})\in V_{hD}\times V_{h0}$ (compare with (8)),

Here $s_{h}$ is a weakly consistent stabilization bilinear form intended to render $a_{h}$ coercive and which is defined on $V_{h}\times V_{h}$ as follows:

where $h_{F}$ is the diameter of the facet $F\in\mathcal{F}$. For an interior facet $F\in\mathcal{F}^{i}$, writing $c_{-}$ and $c_{+}$ the two mesh cells sharing $F$, i.e., $F=\partial c_{-}\cap\partial c_{+}$, and orienting $F$ by the unit normal vector $n_{F}$ pointing from $c_{-}$ to $c_{+}$, one has

The sign of the jump is irrelevant in what follows. The role of the summation over the interior facets in (19) is to penalize the jumps of the cell reconstruction $\mathfrak{R}$ across the interior facets. For a boundary facet $F\in\mathcal{F}^{b}$, we denote $c_{-}$ the unique mesh cell containing $F$, we orient $F$ by the unit normal vector $n_{F}:=n_{c_{-}}$ which points outward $\Omega$, and we define

The role of the summation over the boundary facets in (19) is to penalize the jumps between the cell reconstruction $\mathfrak{R}$ and the boundary facet reconstruction $\mathfrak{R}^{\partial}$. The parameter $\eta>0$ in (19) is user-defined with the only requirement that $\eta>0$. The bilinear form $s_{h}$ is classical in the context of discontinuous Galerkin methods (see [3, 10] for instance, see also [9] for cell-centred Galerkin methods). In practice, the penalty parameter $\eta$ scales as $\eta=\beta\mu$ where $\mu$ is the second Lamé coefficient of the material and $\beta$ is a dimensionless factor that remains user-dependent. Notice that this choice is robust with respect to $\nu\to 0.5$. We present a numerical test in Section 3.5.1 illustrating the moderate impact of $\beta$ on the numerical computations.

We can next define a discrete mass bilinear form $m_{h}$ similar to (7) and a discrete load linear form $l_{h}(t)$ similar to (10); details are given below. Then the space semi-discrete version of the evolution problem (9) amounts to seeking $(u_{h},\varepsilon_{p,\mathcal{C}},p_{\mathcal{C}}):(0,T)\rightarrow V_{hD}\times Q_{h}\times P_{h}$ such that, for all $t\in(0,T)$,

where the time-dependency is left implicit in the second and third equations and where we introduced in every mesh cell $c\in\mathcal{C}$ the local stress tensor

Note that the plasticity relations in (22) are enforced cellwise, i.e., a mesh cell $c\in\mathcal{C}$ is either in the elastic state or in the plastic state depending on the value of $\varphi(\Sigma_{c}(u_{h}),p_{c})$. This is due to the fact that stresses are cell-wise constant and thus the second relation of (22) can only be enforced cell-wise.

The definition of the discrete mass bilinear form $m_{h}$ hinges on subdomains to condense the mass associated with the dofs. Figure 2 represents our choice for the subdomains. For all the interior cells, the subdomain $\omega_{c}$ is chosen as the whole cell, i.e., $\omega_{c}=c$. For the boundary vertices and for the cells having a boundary face, a dual barycentric subdomain is constructed, leading to subdomains denoted by $\omega_{z}$ and $\omega_{c}$, respectively (see Figure 2). For the discrete load linear form, we compute averages of the external loads $f$ and $g_{N}$ in the mesh cells and on the Neumann boundary facets, respectively. As a consequence, $m_{h}(\cdot,\cdot)$ and $l_{h}(t;\cdot)$ can be written as follows for all $(v_{h},\tilde{v}_{h})\in V_{h}\times V_{h}$ (compare with (7) and (10)):

with $m_{z}:=\int_{\omega_{z}}\rho$, $m_{c}:=\int_{\omega_{c}}\rho$, $f_{c}(t):=\int_{c}f(t)$ and $g_{F}(t):=\int_{F}g_{N}(t)$.

The reconstruction operator $\mathcal{R}$ is constructed in the same way as in the finite volume methods studied in [13, Sec. 2.2] and in the cell-centered Galerkin methods from [9]. For a given facet $F\in\mathcal{F}$, we select neighbouring boundary vertices collected in a subset denoted $\mathcal{Z}_{F}^{\partial}$ and neighbouring cells collected in a subset denoted $\mathcal{C}_{F}$, as well as coefficients $(\alpha_{F}^{z})_{z\in\mathcal{Z}_{F}^{\partial}}$ and $(\alpha_{F}^{c})_{c\in\mathcal{C}_{F}}$, and we set

The neighbouring dofs should stay $\mathcal{O}(h)$ close to the facet $F$. An algorithm is detailed thereafter to explain the selection of the neighbouring dofs in $\mathcal{Z}^{\partial}_{F}$ and $\mathcal{C}_{F}$. The coefficients $\alpha_{F}^{z}$ and $\alpha_{F}^{c}$ are chosen as the barycentric coordinates of the facet barycentre $\mathbf{x}_{F}$ in terms of the location of the boundary vertices in $\mathcal{Z}_{F}^{\partial}$ and the barycentres of the cells in $\mathcal{C}_{F}$. For any interior facet $F\in\mathcal{F}^{i}$, the set $\mathcal{Z}_{F}^{\partial}\cup\mathcal{C}_{F}$ is constructed so as to contain exactly $(d+1)$ points forming the vertices of a non-degenerate simplex. Thus, the barycentric coefficients $\alpha^{\mathbf{z}}_{F}$ and $\alpha^{c}_{F}$ are computed by solving the linear system:

where the vertex and its position are written $\mathbf{z}$ and $\mathbf{x}_{c}$ is the position of the barycentre of the cell $c$.

The main rationale for choosing the neighboring dofs is to ensure as much as possible that all the coefficients $\alpha_{F}^{z}$ or $\alpha_{F}^{c}$ lie in the interval $(0,1)$, so that the definition of $\mathcal{R}(v_{h})_{F}$ in (26) is based on an interpolation formula (rather than an extrapolation formula if some coefficients lie outside the interval $(0,1)$.) For most internal facets $F\in\mathcal{F}^{i}$, far from the boundary $\partial\Omega$, it is possible to choose an interpolation-based reconstruction operator using only cell dofs, i.e., we usually have $\mathcal{Z}_{F}^{\partial}:=\emptyset$. Figure 3 presents an example for an interior facet $F$ using three neighbouring cell dofs located at the cell barycentres $\mathbf{x}^{i}$, $\mathbf{x}^{j}$ and $\mathbf{x}^{k}$. Close to the boundary $\partial\Omega$, the use of boundary vertex dofs helps to prevent extrapolation. In all the cases we considered, interpolation was always possible using the algorithm described below.

On the boundary facets, the reconstruction operator only uses interpolation from the boundary vertices of the facet, i.e., we always set $\mathcal{C}_{F}:=\emptyset$ for all $F\in\mathcal{F}^{b}$. In three space dimensions, the facet can be polygonal and the barycentric coordinates are generalized barycentric coordinates. This is achieved using [5] and the package 2D Triangulation of the geometric library CGAL. In the case of simplicial facets (triangles in three space dimensions and segments in two space dimensions), we recover the classical barycentric interpolation as described above.

The advantage of using interpolation rather than extrapolation is relevant in the context of explicit time-marching schemes where the time step is restricted by a CFL condition depending on the largest eigenvalue $\lambda_{\max{}}$ of the stiffness matrix associated with the discrete bilinear form $a_{h}(0;\cdot,\cdot)$ (see, e.g., (48)). It turns out that using extrapolation can have an adverse effect on the maximal eigenvalue of the stiffness matrix, thereby placing a severe restriction on the time step, and this restriction is significantly alleviated if enough neighboring dofs are used in (26) to ensure that interpolation is being used. We refer the reader to Table 1 below for an illustration.

Let us now briefly outline the algorithm used for selecting the reconstruction dofs associated with a given mesh inner facet $F\in\mathcal{F}^{i}$. This algorithm has to be viewed more as a proof-of-concept than as an optimized algorithm, and we observe that this algorithm is only used in a pre-processing stage of the computations. Fix an integer $I\geq d+1$.

1. 1.

   Compute a list of points $(\mathbf{x}_{i})_{1\leq i\leq I}$ ordered by increasing distance to $\mathbf{x}_{F}$; each point can be either the barycentre of a mesh cell or a boundary vertex. To this purpose, we use the KDTree structure of the scipy.spatial module of Python.

2. 2.

   Check if $\mathbf{x}_{F}$ lies in the convex hull of the set $(\mathbf{x}_{i})_{1\leq i\leq I}$. To this purpose we use the ConvexHull structure of scipy.spatial. If that is not true, then extrapolation must be used. Otherwise find a subset of $(\mathbf{x}_{i})_{1\leq i\leq I}$ containing $(d+1)$ points and use the barycentric coordinates of the resulting simplex to evaluate the interpolation coefficients to be used in (26).

Note that $I$ must be large enough to enable the recovery of at least one simplex per inner facet that is not too flat, independently of the use of extrapolation or interpolation. In our computations, we generally took $I=10$ if $d=2$ and $I=25$ if $d=3$.

To illustrate the performance of the above algorithm in alleviating time step restrictions based on a CFL condition, we report in Table 1 the largest eigenvalue $\lambda_{\max{}}$ of the stiffness matrix and the percentage of the mesh facets where extrapolation has to be used as a function of the integer parameter $I$ from the above algorithm. The results are obtained on the two three-dimensional meshes (called "coarse" and "fine") considered in Section 5.3.1 together with the DEM discretization for the dynamic elasto-plastic evolution of a beam undergoing flexion. Note that the minimal value is $I=7$ for the coarse mesh and $I=9$ for the fine mesh.

In this section we rewrite the first equation in (22) as a particle method by introducing the dofs of the discrete displacement $u_{h}(t)\in V_{hD}$ attached to the mesh cells and to the boundary vertices lying on the Neumann boundary, which we write $U_{\mathrm{DEM}}:=(U_{p}(t))_{p\in\mathcal{P}}$ with $\mathcal{P}:=\mathcal{C}\cup\mathcal{Z}_{N}^{\partial}$ and $\mathcal{Z}_{N}^{\partial}:=\{z\in\mathcal{Z}^{\partial}\ |\ z\in\partial\Omega_{N}\}$. Here $\mathcal{P}$ can be viewed as the indexing set for the set of particles. Note that $U_{\mathrm{DEM}}$ is a collection of dof values forming a column vector, whereas $u_{h}$ is a piecewise-constant function. Recalling the definition of the discrete mass bilinear form, we set $m_{p}:=\int_{\omega_{c}}\rho$ if $p=c\in\mathcal{C}$ and $m_{p}:=\int_{\omega_{\mathbf{z}}}\rho$ if $p=\mathbf{z}\in\mathcal{Z}_{N}^{\partial}$. Concerning the external loads, we set $F_{\mathrm{DEM}}(t):=(F_{p}(t))_{p\in\mathcal{P}}$ with $F_{p}(t):=f_{c}(t)=\int_{c}f(t)$ if $p=c$ and $F_{p}(t):=\sum_{F\in\mathcal{F}_{\mathbf{z}}}\alpha_{F}^{\mathbf{z}}g_{F}(t)=\sum_{F\in\mathcal{F}_{\mathbf{z}}}\alpha_{F}^{\mathbf{z}}\int_{F}g_{N}(t)$ if $p=\mathbf{z}$, where $\mathcal{F}_{\mathbf{z}}\subset\mathcal{F}^{b}$ is the collection of the boundary facets to which $\mathbf{z}$ belongs and the coefficients $\alpha_{F}^{\mathbf{z}}$ are those used in (26) for the facet reconstruction. Since the Neumann boundary $\partial\Omega_{N}$ is relatively open in $\partial\Omega$, all the facets in $\mathcal{F}_{\mathbf{z}}$ belong to $\partial\Omega_{N}$ if $\mathbf{z}\in\mathcal{Z}_{N}^{\partial}$.

Recall that $\varepsilon_{p,\mathcal{C}}:(0,T)\to Q_{h}$ is the discrete tensor of remanent plastic strain. Let us use the shorthand notation $\Sigma_{c}(t):=\Sigma_{c}(u_{h}(t))$ as defined in (23), as well as $\Sigma_{\mathcal{C}}(t):=(\Sigma_{c}(t))_{c\in\mathcal{C}}$. For a piecewise-constant function defined on the mesh cells, say $w_{\mathcal{C}}=(w_{c})_{c\in\mathcal{C}}$, we define the mean-value $\{w_{\mathcal{C}}\}_{F}:=\frac{1}{2}(w_{c_{-}}+w_{c_{+}})$ for all $F=\partial c_{-}\cap\partial c_{+}\in\mathcal{F}^{i}$. Recall that the interior facet $F$ is oriented by the unit normal vector $n_{F}$ pointing from $c_{-}$ to $c_{+}$ and that the jump across $F\in\mathcal{F}^{i}$ is defined such that $[w_{\mathcal{C}}]_{F}:=w_{c_{-}}-w_{c_{+}}$. Recall also that for a boundary facet $F\in\mathcal{F}^{b}$, $c_{-}$ denotes the mesh cell to which $F$ belongs and that $n_{F}$ is the unit normal vector to $F$ pointing outward $\Omega$. Then a direct calculation shows that for all $\tilde{v}_{h}\in V_{h0}$,

To simplify some expressions, we are going to neglect the second term on the above right-hand side since this term is of higher-order (it is essentially the product of two jumps). Recall that, by definition, the reconstruction operator $\mathcal{R}$ on a boundary facet $F\in\mathcal{F}^{b}$ makes use only of the vertex dofs of that facet. Then, letting $(\tilde{V}_{p})_{p\in\mathcal{P}}$ be the collection of the dofs of the discrete test function $\tilde{v}_{h}$, we infer that

where $\Phi_{p}^{\mathrm{ep}}(t)$ is the elasto-plastic force acting on the particle $p\in\mathcal{P}$ and $\Phi_{p}^{\mathrm{pen}}(t)$ is the force induced by the penalty and acting on the same particle. For all $c\in\mathcal{C}$, we have $\Phi_{c}^{\mathrm{ep}}(t):=\sum_{F\in\mathcal{F}_{c}^{i}\cup\mathcal{F}_{c}^{N}}\Phi_{c,F}^{\mathrm{ep}}(t)$ with $\mathcal{F}_{c}^{i}:=\{F\in\mathcal{F}^{i}\ |\ F\subset\partial c\}$, $\mathcal{F}_{c}^{N}:=\{F\in\mathcal{F}^{b}\ |\ F\subset\partial c\cap\partial\Omega_{N}\}$, and

with $\iota_{c,F}:=n_{c}\cdot n_{F}$, and for all $z\in\mathcal{Z}_{N}^{\partial}$, we have

with $\Phi^{\mathrm{ep}}_{c_{-},F}$ defined in (30) (recall that $c_{-}$ denotes the unique mesh cell containing the boundary facet $F\in\mathcal{F}^{b}$). Note that the principle of action and reaction is encoded in the fact that $\iota_{c_{-},F}+\iota_{c_{+},F}=0$ for all $F=\partial c_{-}\cap\partial c_{+}\in\mathcal{F}^{i}$.

The penalty force is composed of two terms, that is, for all $p\in\mathcal{P}$, $\Phi_{p}^{\mathrm{pen}}(t):=\Phi_{p}^{\mathrm{pen,1}}(t)+\Phi_{p}^{\mathrm{pen,2}}(t)$. We define the first term for every cell $c\in\mathcal{C}$ and every facet $F\in\mathcal{F}$ such that $F\subset\partial c$ as

and we define it for every Neumann boundary vertex $z\in\mathcal{Z}_{N}^{\partial}$ as

The second term is more intricate since it is a byproduct of the extended stencil of the method. The set of facets using the dof of the particle $p\in\mathcal{P}$ (whether cell or boundary vertex) is defined as $\mathfrak{F}_{p}$. This means that if $F\in\mathfrak{F}_{p}$, then one has either $p\in\mathcal{C}_{F}$ or $p\in\mathcal{Z}_{F}^{\partial}$. Let us recall that for a facet $F\in\mathcal{F}$, the cell dofs used for the reconstruction $\mathcal{R}(u_{h}(t))_{F}$ are collected in $\mathcal{C}_{F}$ and the boundary vertex dofs in $\mathcal{Z}_{F}^{\partial}$. Then the second penalty term writes for all $p\in\mathcal{P}$ and all $F\in\mathfrak{F}_{p}$ as

where $\alpha_{F}^{p}$ is the barycentric coordinate associated with the particle $p\in\mathcal{P}$ in the reconstruction over the facet $F$. As a consequence, the total second penalty term writes for all $p\in\mathcal{P}$ as

Finally, putting everything together, we infer that the the first equation in (22) becomes