Axisymmetric Virtual Elements For Problems of Elasticity and Plasticity

VEM is employed to solve 2D elasticity problems. Recall, for elasticity problems, the standard weak form: Find $\mathbf{u}\in\bm{\mathcal{V}}$ such that

where the bilinear form is

and the linear form is

An example geometric domain of interest is provided in Figure 1a. Calculating field variables such as stress, strain, and displacements are the objective. By discretizing the domain with polygons, Figure 1b, it is then possible to use VEM. For VEM the polygons can be convex or non-convex with arbitrary number of sides. The interpolation space is necessarily populated with polynomial as well as non-polynomial functions as a consequence of selecting polygonal elements for discretization. The resulting polygon elements are able to compatibly interface with adjacent elements along their sides. Essentially, the polygon edges have polynomial interpolation functions, but polynomial plus (possibly) non-polynomial functions exist on the element interior. Nevertheless, it is only necessary to know polynomial functions along the edges to implement VEM. Consequently, the edge polynomials are chosen by selecting a polynomial order, $k$. The present work chooses first order polynomials exclusively.

On element E, scalar-valued polynomials of order equal to $k$ or less are in the space $\mathcal{P}_{k}(E)$. Consequently, $\bm{\mathcal{P}}_{k}\equiv[\mathcal{P}_{k}]^{2}$ is used to represent a 2D vector space of polynomials with two variables. The expression $\mathbf{P}_{k}=\{\mathbf{p}_{\alpha}\}_{\alpha=1,...,n_{k}}$ denotes a basis for the polynomial space. Illustrated below is a polynomial basis of order $k=1$.

or

In equation (5), scaled monomials are used to create the components. The scaled monomials are defined as

where $\bar{\mathbf{x}}=(\bar{x},\bar{y})$ is the element $E$ centroid location and the diamter is $h_{E}$(i.e., diameter of smallest circle that encloses all element vertices).

Before the stiffness matrix is presented, it is necessary to define a variety of matrices and operators^{8, 9}. To begin, recall that the modular matrix for plane stress is

and for plane strain is

where $E_{Y}$ is Young’s modulus of elasticity, and $\nu$ is Poisson’s ratio.

Degree of Freedom Operator, $dof_{i}$. Define an operator, $dof_{i}(\mathbf{p}_{j})$, which is found as follows: first, calculate polynomial vector $\mathbf{p}_{j}$ evaluated at the vertex coordinates associated with dof $i$, second, take the component of the vector that is directed along dof $i$. The operator, $dof_{i}$, extracts the value of its argument, at dof $i$, in the direction of dof $i$, and the dofs (degrees of freedom) range from $1$ to $2n_{v}$.

The $\mathbf{D}$ matrix. At coordinates of the degrees of freedom of polygon $E$ polynomial vector values are computed. These quantities are used to construct matrix $\mathbf{D}$. The result is

The engineering strain operator, $\bm{\varepsilon}$. An engineering strain operator, $\bm{\varepsilon}$, is applied to the polynomial basis as follows:

Acting on the polynomial base functions, the strain operator produces the matrix,

where in (11), $p_{i,j}$ represents the $j$th (1st or 2nd) component for polynomial vector $\mathbf{p}_{i}$. For example, in 2D the polynomial vectors $\mathbf{p}_{i}$ are shown in equations (4) and (5). Depending on the context, the strain operator $\bm{\varepsilon}$ has a specific use. For single vectors, like $\mathbf{p}_{1}$ or $\mathbf{v}^{h}$, in 2D containing two components, the operator output is $3\times 1$. For a group of vectors, like $\mathbf{P}_{1}$, as shown in (11), the operator output is $3\times n_{k}$, where $n_{k}=6$ for the first order polynomial base.

The engineering stress operator, $\bm{\sigma}$. An engineering stress operator, $\bm{\sigma}(\mathbf{p}_{\alpha})$, is obtained by matrix multiplication

In (12), the correct modular matrix, $\mathbf{C}$, for plane strain or plane stress is inserted. The engineering stress matrix is formed by using the results of (12)

This prepares the way for using (13) in (14).

The $\tilde{\mathbf{B}}$ matrix. The $\tilde{\mathbf{B}}$ matrix is formed, for rows $\alpha=1$ to $n_{k}$ two columns at a time, for $j$ values ranging over vertex numbers $1$ to $n_{v}$, according to the following formulas:

For equation (14), with reference to Figure 2, edge length $j$ is written as $|e_{j}|$ and the $j$th edge outward unit normal vector components are $[n_{e1}\ n_{e2}]^{T}_{j}$. For vertex $j=1$, edge length $|e_{j-1}|$ is the length from vertex $j=n_{v}$ to $j=1$, where the polygon element has $n_{v}$ vertices, and normal $n_{e_{j-1}}$ is the outward unit normal to the edge from $j=n_{v}$ to $j=1$. The terms $[1\ \ 0]^{T}$ or $[0\ \ 1]^{T}$ are dotted with the vector that arises from the right side of (14), as shown.

The $\breve{\mathbf{B}}$ matrix. The $\breve{\mathbf{B}}$ matrix is computed as

where

The $\bar{\mathbf{B}}$ matrix. A $\bar{\mathbf{B}}$ matrix is constructed by creating the $\tilde{\mathbf{B}}$ matrix and then replacing the first three rows with the $\breve{\mathbf{B}}$ matrix. The infinitesimal strain, $\bm{\varepsilon}(\mathbf{p}_{\alpha})$, is zero for $\alpha=1,2,3$ of the polynomial basis, $\mathbf{P}_{1}$. It is for this reason that the beginning three rows of matrix $\tilde{\mathbf{B}}$ given in (14) are populated with $\breve{\mathbf{B}}$ to create the $\bar{\mathbf{B}}$ matrix. If this was not done the first three rows would be zero (see ^{8, 9} for more details). Note, the $\bar{\mathbf{B}}$ matrix described herein should not be confused with the $\bar{\mathbf{B}}$ matrix ¹⁴ used with finite elements and incompressibility problems.

The projector matrix, $\tilde{\mathbf{\Pi}}$. The projector matrix, $\tilde{\mathbf{\Pi}}$, is constructed as

where

Engineering stress calculation, $\bm{\sigma}(\mathbf{v}^{h})$. The stress terms needed to calculate the $\tilde{\mathbf{B}}$ matrix are found by using equation (12). Yet, actual engineering stresses, calculated during post processing of elasticity problems, are found by using the applicable modular matrix $\mathbf{C}$ and equation (20). The resulting formula for engineering stresses is

Element stiffness matrix, $\mathbf{k}_{E}$. Finally, analogous to finite element analysis, the element stiffness matrix for VEM is found found for each polygon element. The element stiffness for VEM is calculated from two contributions: a consistency part (23) and a stability part (24).

The consistency part is written as

where $\mathbf{C}$ is the 2D plane stress or plane strain elastic modular matrix (see equations (7) and (8)), $A_{E}$ is the area of polygon element $E$, and the polygon element has thickness $t$.

The element stability stiffness part, given by Sukumar and Tupek ¹⁰, is given as

where the diagonal matrix $\mathbf{S}^{d}_{E}$ is $2n_{v}\times 2n_{v}$ in size and is scaled as needed. The matrix has diagonal terms: $(\mathbf{S}^{d}_{E})_{ii}=\mathsf{max}(\alpha_{0}\ \text{tr}(\mathbf{C})/m,(\mathbf{k}^{c}_{E})_{ii})$, where $m=3$ in 2D, the modular matrix $\mathbf{C}$ in 2D is the applicable case for plane strain or plane stress, and since elements are constructed from scaled monomials which create diameters on the order of 1, $\alpha_{0}=1$.

The above stability part of the element stiffness matrix is found to be effective for compressible materials in linear elasticity problems. It is also appropriate for plasticity problems. Later, in section 5, similar to Park et al.¹⁵, a stabilization matrix like (24) is provided for axisymmetric problems of near incompressibility.

Body forces and tractions impose external forces according to equation (3). Point loads at individual nodes (or vertices) of polygon elements are allowed also. Standard finite element methods are used to accomplish application of forces. The following linear form includes the three types of forces mentioned.

Mengolini et al. ⁸ provides additional discussion on the topic of load application. In the examples provided later, point loads or pressure loads are used. An external force vector, $\mathbf{F}_{ext}$, is assembled from the point or pressure loads. Subsequently, the solution for nodal displacements is found by using the external force vector and the global stiffness. In anticipation of nonlinear problems, such as plasticity, equilibrium between external and internal forces is accomplished by use of a Newton-Raphson scheme.

For nonlinear problems, in order to take advantage of an implicit analysis with iterations for equilibrium, local element internal forces are required. For plastic problems results from the stress integration process are used and a stabilization matrix adjustment is included. As a result, the internal forces for a single element $E$ are,

where for element $E$ the polygon area is $A_{E}$, the engineering stress vector is $\bm{\sigma}=[\sigma_{x}\ \sigma_{y}\ \sigma_{xy}]^{T}$, and $\mathrm{d}\mathbf{u}^{E}$ is the latest incremental displacemen. Note that, for a particular time step in the nonlinear analysis, the stress vector is constant across the given element. Then, analogous to FEM, the global internal force vector is assembled from the individual element internal force vectors using the assembly operator ¹⁴. That is,

Consider Figure 4 with an axisymmetric element shown. The example element shown in the figure has an axisymmetric cross-section with 4 nodes. In general, the cross-section has an arbitrary polygonal cross-section with nodes $1$ to $n_{v}$. For the cross-section the nodal displacements in the radial direction and $z$ direction are represented in terms of nodal shape functions as follows:

where $\mathbf{d}_{r}=[u_{r}^{1}\ u_{r}^{2}\ \cdots\ u_{r}^{n_{v}}]^{T}$ and $\mathbf{d}_{z}=[u_{z}^{1}\ u_{z}^{2}\ \cdots\ u_{z}^{n_{v}}]^{T}$.

The 2D coordinates are relabeled in terms of the axisymmetric variables $r\equiv x$, $z\equiv y$. As a result, for a given element, the axisymmetric cross-sectional strains are expressed as

where the vector $\mathbf{d}$ organizes $r$ and $z$ displacements as shown. Notice the strain displacement matrix, $\mathbf{B}_{2}$. Aside from coordinate direction relabeling, it is conceptually identical to the 2D elasticity VEM strain displacement matrix determined in section 2.3, equation (19). The shape functions themselves are not known, but the strain displacement matrix is computable from the VEM formulation.

For a given element in the VEM discretization, the tangential or circumferential strain is

where $C_{o}$ is the original circumference, $C_{f}$ is the final (or current) circumference, and $R_{E}$ is the distance from the axis of symmetry to the element centroid (see Figure 4), and the shape functions associated with tangential strain are evaluated at the element centroid coordinates. Note also that the strain displacement matrix $\mathbf{B}_{t}$ requires shape function values. The shape function value are not available from the VEM formulation. Hence, an alternative means of calculating shape function values is needed. The mean value coordinates (MVC) shape functions of section 3 are used.

The strain displacement matrix for the axisymmetric case is constructed by combining the results from the in-plane case, $\mathbf{B}_{2}$, which comes from the VEM formulation, and the circumferential or tangential strain case, $\mathbf{B}_{t}$, which is constructed by using MVC shape functions. The result is the axisymmetric strain displacement matrix

The axisymmetric modular matrix is

where $\zeta=\frac{E_{Y}(1-\nu)}{(1+\nu)(1-2\nu)}$.

The axisymmetric stresses are expressed as

For the condition of plane strain a linear elastic cylinder is subject to internal pressure, $p=10$. The cylinder has inner radius $a=4$, and outer radius, $b=10$. The modulus of elasticity $E_{Y}=1000$ and Poisson’s ratio $\nu=0.2$ or $\nu=0.49999$. Results are presented for the case of plane strain in the $r-\theta$ plane. Boundary conditions are set so that nodes at the top and bottom of the discretization are restrained in the $z$ direction causing $\varepsilon_{z}=0$ everywhere. Plots for Figures 5abcde, are accomplished with 300 VEM elements. The theoretical solution ²⁵, for radial displacements and stresses, is

For problems of near incompressibility the results for Possion’s ratio equal to 0.49999 are handled without trouble. Convex and concave results are accomplished with indistinguishable results. The convergence characteristics shown in Figure 5f are within expected ranges for both compressible and near incompressible values of Poisson’s ratio and are not affected by concave or convex elements.

In Figures 6ab, the top half of an axisymmetric sphere is modeled with 1000 polygonal elements. The applied internal pressure is $p=10$. The sphere internal radius $a=4$, the external radius $b=10$, modulus of elasticity $E_{Y}=1000$, and Poisson’s ratio $\nu=0.2$ or $\nu=0.49999$. Vertical supports are provided at all nodes for which $z=0$, and horizontal supports are provided at all nodes for which $r=0$. The theoretical solution ²⁵, for radial displacements and stresses, is

An elastic circular plate is loaded with a point load at the center, as shown in Figure 7a. The radius of the plate is $r_{o}=10$. The plate has thickness $t=0.25$ and modulus of elasticity $E_{Y}=1000$. Poisson’s ratio used in the various results are indicated in Figures 7cde. In particular, for some results Poisson’s ratios as high as 0.49999 are used to demonstrate that the formulation is not susceptible to locking for near incompressible problems. Results for displacements, radial moments, and tangential moments versus radial distance, $r$, are shown to be in good agreement with analytical solutions. Furthermore, Figures 7cde illustrate results for plates with simply supported and fixed supported conditions. These problems are accomplished with an axisymmetric model using 18000 convex polygonal elements.

A study of normalized displacements and various Poisson ratios is provided in Table 1. For all results of the table, the circular plate is modeled with 9000 convex or mostly concave elements, as indicated. The results are in good agreement with the theoretical displacements which account for both bending and shear displacements. The axisymmetric VEM results have error of at most 3.76 % and for most cases less than 2%. The axisymmetric VEM formulation using concave or convex elements provides nearly indistinguishable results.

For the circular plate loaded by a point load at the center, the center displacement (bending plus shear) is $u_{z}$, the radial moment is $M_{r}$, and the tangential moment is $M_{t}$. The theoretical solutions ^{26, 27} for an elastic circular plate with edges simply supported are

For the case of an elastic circular plate with edges fixed supported, the theoretical solutions ^{26, 27} are

For the condition of plane strain an elasto-plastic cylinder with 250 convex elements (Figure 8a) is loaded by an internal pressure, $p$. The cylinder has inner radius $a=4$, and outer radius, $b=10$. The modulus of elasticity $E_{Y}=1000$, Poisson’s ratio is $\nu=0.3$, the yield stress is $\sigma_{yield}=10$, and for the case of perfect plasticity the hardening modulus is $E_{h}=0$. When a cylinder begins to yield, the ‘plastic front’ location $c$ progresses from the inner radius $a$ toward the outer radius $b$. Hence, for certain pressures, $a<c<b$. For the case of pressure, $p=9.29$, the plastic front is located at $c=6.86$, which is in good agreement with the tangential stress results of Figures 8b and d. Radial stresses, $\sigma_{rr}$ versus radial coordinate, $r$, are in good agreement with analytical results as shown in Figure 8c.

The theoretical solution of an elasto-plastic cylinder is given by Hill ²⁸. Consider the pressure level associated with Figures 8bcd where the cylinder is yielded and the location of the plastic front, $c$, is such that $a<c<b$. For the elastic region, where $r\geq c$, the theoretical solution is

For the plastic region, where $r<c$, the theoretical solution is