A node-based uniform strain virtual element method for elastoplastic solids

The weak form (1) is the continuous problem. The discrete problem is formulated on a partition of the domain $\Omega$ into nonoverlapping elements with arbitrary number of edges (convex or non-convex polygons). This partition is denoted by $\mathcal{T}_{h}$, where $h$ is the maximum diameter of any element in the partition. An element in the partition is denoted by $E$ and its boundary by $\partial E$. $|E|$ is the area of the element and $N_{E}^{V}$ its number of edges/nodes. The unit outward normal to the element boundary in the Cartesian coordinate system is denoted by $\bm{n}=[n_{1}\quad n_{2}]^{\mathsf{T}}$. Fig. 1 depicts an element with seven edges ($N_{E}^{V}=7$), where the edge $e_{a}$ of length $|e_{a}|$ and the edge $e_{a-1}$ of length $|e_{a-1}|$ are the element edges incident to node $a$, and $\bm{n}_{a}$ and $\bm{n}_{a-1}$ are the unit outward normals to these edges, respectively.

Following a standard Galerkin approach, we assume approximations of $\bm{u}$ and $\bm{v}$ on the element, as follows:

where $\{\phi_{a}(\bm{x})\}_{a=1}^{N_{E}^{V}}$ are basis functions that form a partition of unity. A peculiarity of the VEM is that the basis functions are never computed, which is why they are considered virtual. For the method to work, we only need to assume their behavior on the element boundary. For linearly precise VEM approximations, the basis functions on the element boundary are assumed to be

• $\bullet$

  piecewise linear (edge by edge),

• $\bullet$

  continuous on the element edges,

which means that the basis functions possess the Kronecker-delta property on the element edges, and hence they behave like the one-dimensional hat function.

At the element level, the following discrete local spaces are defined:

The discrete local spaces are assembled to form the following discrete global spaces:

Using the preceding definitions, the discrete version of the weak form (1) reads: find $\bm{u}_{h}\in\mathscr{V}_{h}$ such that

To obtain the discrete weak form for the VEM, we follow the standard VEM literature (see for instance, Ref. [61]). We first define a projection operator $\Pi$ onto the space of polynomials of degree 1. To this end, let $[\mathcal{P}(E)]^{2}$ represent the space of polynomials of degree 1 over the element $E$. The projection operator $\Pi$ is defined as:

$\Pi$ is then used to split the displacement approximation on the element, as follows:

where $\Pi\bm{u}_{h}$ is the polynomial part of $\bm{u}_{h}$ (of degree 1) and $\bm{u}_{h}-\Pi\bm{u}_{h}$ contains its remainder terms. The remainder terms can contain polynomials of order greater than 1 or even nonpolynomial terms. The actual form of the projection $\Pi\bm{u}_{h}$ is obtained from the orthogonality condition: $\forall\bm{p}\in[\mathcal{P}(E)]^{2}$,

which, at the element level, gives [22, 62, 63, 64]

where $\bar{x}_{1}$ and $\bar{x}_{2}$ are the components of the mean of the values that the position vector $\bm{x}=\left[x_{1}\quad x_{2}\right]^{\mathsf{T}}$ takes over the vertices of the element; i.e.,

where $\bm{x}_{a}=\left[x_{1a}\quad x_{2a}\right]^{\mathsf{T}}$ are the coordinates of node $a$; $\bar{u}_{1}$ and $\bar{u}_{2}$ are the components of the mean of the values that the displacement approximation $\bm{u}_{h}=\left[u_{1h}\quad u_{2h}\right]^{\mathsf{T}}$ takes over the vertices of the element; i.e.,

In other words, $\bar{\bm{x}}$ and $\bar{\bm{u}}$ represent the geometric center of the element and its associated displacement vector, respectively; the terms $\widehat{\varepsilon}_{ij}$ are components of the element average

and $\widehat{\omega}_{12}$ is the component of the element average

where $\bm{\omega}(\bm{u}_{h})$ is the skew-symmetric tensor that represents rotations.

On substituting (6) into (4), and using the orthogonality condition (7) and noting that $\bm{u}_{h}$ and $\bm{v}_{h}\in[\mathcal{P}(E)]^{2}$, leads to the following VEM representation of the discrete weak form: find $\bm{u}_{h}\in\mathscr{V}_{h}$ such that

where $s_{E}(\bm{u}_{h}-\Pi\bm{u}_{h},\bm{v}_{h}-\Pi\bm{v}_{h})$ is a computable approximation to $a_{E}(\bm{u}_{h}-\Pi\bm{u}_{h},\bm{v}_{h}-\Pi\bm{v}_{h})$ and is meant to provide stability.

The VEM as described above is prone to volumetric locking in the limit $\nu\to 1/2$. Using the virtual element mesh and considering a typical nodal vertex $I$, the NVEM that was proposed in Ref. [60] applies a nodal averaging operator $\pi_{I}$ to (13) that precludes volumetric locking without introducing additional degrees of freedom¹¹1The nodal averaging operator permits to integrate the weak form integrals directly at the nodes. This results in a total number of incompressiblity constraints equal to the number of nodes in the mesh. If this number divides the total number of displacement equations, two degrees of freedom every one constraint is obtained in two dimensions, which is the optimal ratio to perform well in incompressible and nearly incompressible settings [65].. This leads to a nodal version of (13), as follows: find $\bm{u}_{h}\in\mathscr{V}_{h}$ such that

where the notations $a_{I}$, $s_{I}$ and $\ell_{I}$ are introduced as the nodal counterparts of $a_{E}$, $s_{E}$ and $\ell_{E}$, respectively. The construction of the nodal averaging operator is described next.

Each node of the mesh is associated with their own patch of virtual elements. The patch for node $I$ is denoted by $\mathcal{T}_{I}$ and is defined as the set of virtual elements connected to node $I$ (see Fig. 2). Each node of a virtual element $E$ in the patch is assigned the area $\frac{1}{N_{E}^{V}}|E|$; that is, the area of an element is uniformly distributed among its nodes. The representative area of node $I$ is denoted by $|I|$ and is computed by addition of all the areas that are assigned to node $I$ from the elements in $\mathcal{T}_{I}$; that is,

Similarly, each node of a virtual element $E$ is uniformly assigned the strain $\frac{1}{N_{E}^{V}}\widehat{\bm{\varepsilon}}(\bm{u}_{h})$. On considering each strain assigned to node $I$ from the elements in $\mathcal{T}_{I}$, the node-based uniform strain is defined as follows:

Since $\widehat{\bm{\varepsilon}}(\bm{u}_{h})$ is by definition given at the element level, then from (16) the following nodal averaging operator is proposed:

where $[\,\cdot\,]_{{}_{E}}$ denotes evaluation over the element $E$.

The NVEM nodal stiffness matrix is developed by substituting the discretizations (3), the projection operator (8) and the nodal averaging operator (17) into the left-hand side of (14) for a node $I$. This gives

where $\bm{d}$ and $\bm{q}$ are column vectors of element nodal displacements and element nodal values associated with $\bm{v}_{h}$, respectively; $\bm{K}_{I}$ is the NVEM nodal stiffness matrix given by

where $\bm{D}$ is the constitutive matrix and $\bm{S}$ is the stability matrix both defined in Section 3; $\bm{B}_{I}=\pi_{I}[\,\bm{B}\,]$ and $(\bm{I}-\bm{P})_{I}=\pi_{I}[\,\bm{I}-\bm{P}\,]$, where $\bm{I}$ is the identity $(2N_{E}^{V}\times 2N_{E}^{V})$ matrix, and $\bm{B}$ and $\bm{P}$ are defined as

where $q_{ia}=\frac{1}{|E|}\int_{\partial E}\phi_{a}(\bm{x})n_{i}\,ds$ and can be exactly computed on the element boundary using a trapezoidal rule giving the following algebraic expression:

where $n_{ia}$ is the $i$-th component of $\bm{n}_{a}$ and $|e_{a}|$ is the length of the edge incident to node $a$ as defined in Fig. 1;

where

and

Similarly, the NVEM nodal force vector is developed using the right-hand side of (14) for a node $I$, which leads to

where $\bm{q}$ is the column vector of element nodal values associated with $\bm{v}_{h}$, and $\bm{f}_{I}$ is the NVEM nodal force vector associated with the body force and external tractions defined as

For computing the nodal body force vector $\bm{f}_{I}^{b}$ in (27), $\bar{\bm{N}}_{I}=\pi_{I}[\,\bar{\bm{N}}\,]$ and $\widehat{\bm{b}}_{I}=\pi_{I}[\,\widehat{\bm{b}}\,]$, where

and

Regarding the nodal traction force vector $\bm{f}_{I}^{t}$ in (27), the nodal components are now computed with respect to the one-dimensional domain on the Neumann boundary; that is, the representative nodal area reduces to a representative nodal length $|I_{\Gamma}|=\sum_{e\in\mathcal{T}_{I}}\frac{1}{2}|e|$, where $e$ is an element’s edge located on the Neumann boundary and $|e|$ its length; $\mathcal{T}_{I}$ now represents the set of edges connected to node $I$ on the Neumann boundary. Using the preceding definitions, the nodal averaging operator on the Neumann boundary is defined as

where $[\,\cdot\,]_{e}$ denotes evaluation over the edge $e$. The remainder nodal matrices are then obtained as $\bar{\bm{N}}_{\Gamma,I}=\pi_{I,\Gamma}\bigl{[}\bar{\bm{N}}_{\Gamma}\bigr{]}$ and $\widehat{\bm{t}}_{N,I}=\pi_{I,\Gamma}\bigl{[}\,\widehat{\bm{t}}_{N}\bigr{]}$, where

and

For linear elastostatics as developed in Ref. [60], the discrete equilibrium equation is obtained by the discretization of (14). This is accomplished by summing (18) and (26) through all the nodes in the domain and invoking the arbitrariness of $\bm{q}$. This results in the following system of equations:

Eq. (33) is also the discrete equilibrium equation for elastoplasticity when $\bm{D}$ and $\bm{S}$ are nonlinear functions obtained from the elastoplastic constitutive law. As usual, the solution for this problem requires the linearization of (33). Doing this gives

where $\breve{\bm{D}}^{ep}$, $\breve{\bm{S}}$, and $\breve{\bm{\sigma}}$, which are given in Section 3, are the elastoplastic consistent tangent operator, the stability matrix, and the nonlinear stress, respectively, evaluated at node $I$ (by means of the elastoplastic constitutive law) using the node-based uniform strain $\widehat{\bm{\varepsilon}}_{I}$.

The linearized equilibrium equation (34) is used to solve the equilibrium state $\bm{d}_{n+1}^{(k)}=\bm{d}_{n+1}^{(k-1)}+\Delta\bm{d}^{(k)}$ at time $t_{n+1}$ with a time increment $\Delta t=t_{n+1}-t_{n}$ via Newton-Raphson iterations, as follows:

where

The strain tensor $\bm{\varepsilon}$ is split into elastic ($\bm{\varepsilon}^{e}$) and plastic ($\bm{\varepsilon}^{p}$) parts; that is,

The elastic part is governed by the standard linear elastic law, and the plastic part by the von Mises model with mixed linear hardening (for details on this model, see for instance Ref. [68]). The yield function for this model is

where $\bm{\beta}$ is the backstress tensor, $\sigma_{y}$ is a function of the accumulated plastic strain $\bar{\varepsilon}^{p}$ and defines the radius of the yield surface, and $\bm{\eta}$ is the relative stress given by

where $\bm{s}$ is the deviatoric stress. The plastic flow is described by the following associative law:

where $\dot{\bm{\varepsilon}}^{p}$ is the rate of the plastic strain tensor and $\dot{\gamma}$ is the plastic multiplier.

The mixed linear hardening combines linear isotropic and linear kinematic hardening models. The linear isotropic hardening model is defined by the linear function

where $\sigma_{y0}$ is the initial yield stress and $H_{i}$ is the linear isotropic hardening modulus. The linear kinematic hardening model describes the evolution law for the backstress as the following linear function:

where $H_{k}$ is the linear kinematic hardening modulus.

Let $n$ and $n+1$ be the subindices that denote the previous and current states, respectively. These subindices are used for labeling the time at which the variables that are involved in the constitutive law are evaluated. Using the preceding notation, the elastoplastic consistent tangent operator (in Voigt notation) for the above model is

where $G$ and $K$ are the shear modulus and bulk modulus of the material, respectively. Here, we use the three-dimensional constitutive law from where the plain strain state components are extracted. The remainder quantities that appear in (43) are defined as follows: