A nodal integration scheme for meshfree Galerkin methods using the virtual element decomposition

TL;DR

This paper introduces a new nodal integration scheme for meshfree Galerkin methods based on virtual element decomposition, improving accuracy and robustness in linear and nonlinear solid mechanics simulations.

Contribution

The paper develops a novel nodal integration approach using virtual element decomposition, applicable to any linear meshfree approximant, enhancing stability and accuracy.

Findings

Achieves optimal convergence in benchmark problems

Demonstrates improved reliability over standard methods

Effective in linear and nonlinear solid mechanics simulations

Abstract

In this paper, we present a novel nodal integration scheme for meshfree Galerkin methods that draws on the mathematical framework of the virtual element method. We adopt linear maximum-entropy basis functions for the discretization of field variables, although the proposed scheme is applicable to any linear meshfree approximant. In our approach, the weak form integrals are nodally integrated using nodal representative cells that carry the nodal displacements and state variables such as strains and stresses. The nodal integration is performed using the virtual element decomposition, wherein the bilinear form is decomposed into a consistency part and a stability part that ensure consistency and stability of the method. The performance of the proposed nodal integration scheme is assessed through benchmark problems in linear and nonlinear analyses of solids for small displacements and small-strain kinematics. Numerical results are presented for linear elastostatics and linear elastodynamics, and viscoelasticity. We demonstrate that the proposed nodally integrated meshfree method is accurate, converges optimally, and is more reliable and robust than a standard cell-based Gauss integrated meshfree method.

Tables

Table 1. Table 1: Relative error in the L 2 superscript 𝐿 2 L^{2} norm for the patch test.

Method	Gauss rule	Regular	Distorted	Unstructured
MEM	1-pt (interior of the cell)	$1.0\times {10}^{-2}$	$2.0\times {10}^{-2}$	$1.7\times {10}^{-2}$
MEM	3-pt (interior of the cell)	$2.3\times {10}^{-3}$	$1.6\times {10}^{-3}$	$1.6\times {10}^{-3}$
MEM	6-pt (interior of the cell)	$5.0\times {10}^{-5}$	$8.0\times {10}^{-4}$	$1.2\times {10}^{-3}$
MEM	12-pt (interior of the cell)	$2.2\times {10}^{-7}$	$3.0\times {10}^{-4}$	$5.0\times {10}^{-4}$
NIVED	1-pt (per edge of the cell)	$3.6\times {10}^{-15}$	$4.1\times {10}^{-15}$	$2.5\times {10}^{-15}$

Table 2. Table 2: Relative error in the H 1 superscript 𝐻 1 H^{1} seminorm for the patch test.

Method	Gauss rule	Regular	Distorted	Unstructured
MEM	1-pt (interior of the cell)	$1.4\times {10}^{-2}$	$5.4\times {10}^{-2}$	$2.5\times {10}^{-2}$
MEM	3-pt (interior of the cell)	$2.6\times {10}^{-3}$	$5.3\times {10}^{-3}$	$4.8\times {10}^{-3}$
MEM	6-pt (interior of the cell)	$5.4\times {10}^{-5}$	$1.9\times {10}^{-3}$	$1.3\times {10}^{-3}$
MEM	12-pt (interior of the cell)	$2.3\times {10}^{-7}$	$7.7\times {10}^{-4}$	$4.5\times {10}^{-4}$
NIVED	1-pt (per edge of the cell)	$3.6\times {10}^{-15}$	$5.2\times {10}^{-15}$	$7.8\times {10}^{-15}$

Table 3. Table 3: NIVED radial displacements at points A 𝐴 A and B 𝐵 B for the thick-walled viscoelastic cylinder problem at time t = 20 𝑡 20 t=20 .

$({\mu }_0,{\mu }_1)=$		$(0.7,0.3)$		$(0.3,0.7)$		$(0.01,0.99)$
mesh	# dofs	$u_A$	$u_B$	$u_A$	$u_B$	$u_A$	$u_B$
Figure 24	882	0.053012	0.031904	0.118830	0.065142	0.646690	0.330050
Figure 24	1922	0.053077	0.031950	0.118970	0.065219	0.647820	0.330170
Figure 24	3362	0.053101	0.031966	0.119020	0.065246	0.648160	0.330210
Figure 24	7442	0.053118	0.031977	0.119050	0.065264	0.648370	0.330240
Figure 24	26082	0.053131	0.031983	0.119070	0.065272	0.648440	0.330210
Ref. (FEM-Q9)	79242	0.053584	0.032425	0.119970	0.066177	0.652880	0.334810

Equations

$${\mathbf{v}}^h(\mathbf{x})=\sum _{a=1}^m{\varphi }_a(\mathbf{x}){\mathbf{v}}_a,$$

$$\underset{\mathbf{\varphi }\in {\mathrm{I}\mathrm{R}}_{+}^m}{\min }\sum _{a=1}^m$$

$$\sum _{a=1}^m{\varphi }_a(\mathbf{x})=1$$

$$w_a(\mathbf{x})=\exp \left(-\frac{\gamma }{h_a^2}\Vert {\mathbf{c}}_a{\Vert }^2\right),$$

$$w_{a}(q)=\left\{\begin{array}[]{lr}1-6q^{2}+8q^{3}-3q^{4}&\quad 0\leq q\leq 1\\ 0&\quad q>1\end{array}\right.,$$

$${\varphi }_a(\mathbf{x},\mathbf{\lambda })=\frac{w_a(\mathbf{x})\exp (-\mathbf{\lambda }(\mathbf{x})\cdot {\mathbf{c}}_a(\mathbf{x}))}{Z(\mathbf{x},\mathbf{\lambda }(\mathbf{x}))},Z(\mathbf{x},\mathbf{\lambda }(\mathbf{x}))=\sum _{b=1}^m{w}_b(\mathbf{x})\exp (-\mathbf{\lambda }(\mathbf{x})\cdot {\mathbf{c}}_b(\mathbf{x})),$$

$${\mathbf{\lambda }}^{\ast }(\mathbf{x})=\arg \underset{\mathbf{\lambda }\in {\mathrm{I}\mathrm{R}}^d}{\min }\ln Z(\mathbf{x},\mathbf{\lambda }),$$

$$\nabla \cdot \mathbf{\sigma }+\mathbf{b}$$

$$\mathbf{u}$$

$$\mathbf{\sigma }\cdot {\mathbf{n}}_{\Gamma }$$

$$\begin{array}{rl}a(\mathbf{u},\mathbf{v})& =\ell (\mathbf{v})\forall \mathbf{v}(\mathbf{x})\in \mathcal{W},\\ a(\mathbf{u},\mathbf{v})={\int }_{\Omega }\mathbf{\sigma }(\mathbf{u}):\mathbf{\epsilon }(\mathbf{v})d\mathbf{x},& \ell (\mathbf{v})={\int }_{\Omega }\mathbf{b}\cdot \mathbf{v}d\mathbf{x}+{\int }_{{\Gamma }_t}\bar{\mathbf{t}}\cdot \mathbf{v}ds,\end{array}$$

$$\mathcal{V}:=\left\{\mathbf{u}(\mathbf{x}):\mathbf{u}\in [H^1(\Omega ){]}^2,\mathbf{u}=\bar{\mathbf{u}}\text{on }{\Gamma }_u\right\},\mathcal{W}:=\left\{\mathbf{v}(\mathbf{x}):\mathbf{v}\in [H^1(\Omega ){]}^2,\mathbf{v}=0\text{on }{\Gamma }_u\right\},$$

$$\mathbf{\epsilon }(\mathbf{v})=\frac{1}{2}\left(\nabla \mathbf{v}+{\nabla }^{\mathsf{T}}\mathbf{v}\right).$$

$${\mathcal{V}}^h{|}_E:=\left\{{\mathbf{u}}^h(\mathbf{x}):{\mathbf{u}}^h\in [H^1(E){]}^2\right\},{\mathcal{W}}^h{|}_E:={\mathcal{V}}^h{|}_E,$$

$${\mathcal{V}}^h$$

$${\mathcal{W}}^h$$

$$a({\mathbf{u}}^h,{\mathbf{v}}^h)=\sum _{E\in {\mathcal{T}}^h}{a}_E({\mathbf{u}}^h,{\mathbf{v}}^h)=\sum _{E\in {\mathcal{T}}^h}|E|{\mathbf{\epsilon }}_E({\mathbf{v}}^h):\mathcal{D}:{\mathbf{\epsilon }}_E({\mathbf{u}}^h)\forall {\mathbf{u}}^h,\forall {\mathbf{v}}^h\in {\mathcal{V}}^h,$$

$$\ell ({\mathbf{v}}^h)=\sum _{E\in {\mathcal{T}}^h}{\ell }_{b,E}({\mathbf{v}}^h)+\sum _{S\in {\mathcal{E}}^h}{\ell }_{t,S}({\mathbf{v}}^h)=\sum _{E\in {\mathcal{T}}^h}|E|{\mathbf{b}}_E\cdot {\mathbf{v}}^h+\sum _{S\in {\mathcal{E}}^h}|S|{\bar{\mathbf{t}}}_S\cdot {\mathbf{v}}^h\forall {\mathbf{v}}^h\in {\mathcal{V}}^h,$$

$$\Pi :{\mathcal{V}}^h{|}_E\to [\mathcal{P}(E){]}^2,\Pi \mathbf{p}=\mathbf{p}\forall \mathbf{p}\in [\mathcal{P}(E){]}^2,$$

$${\mathbf{u}}^h$$

$${\mathbf{v}}^h$$

$$a_E(\mathbf{p},{\mathbf{v}}^h-\Pi {\mathbf{v}}^h)=0\forall \mathbf{p}\in [\mathcal{P}(E){]}^2,{\mathbf{v}}^h\in {\mathcal{V}}^h{|}_E.$$

$$a_E({\mathbf{u}}^h,{\mathbf{v}}^h)=a_E(\Pi {\mathbf{u}}^h,\Pi {\mathbf{v}}^h)+a_E({\mathbf{u}}^h-\Pi {\mathbf{u}}^h,{\mathbf{v}}^h-\Pi {\mathbf{v}}^h).$$

$$a_E^h({\mathbf{u}}^h,{\mathbf{v}}^h)=a_E(\Pi {\mathbf{u}}^h,\Pi {\mathbf{v}}^h)+s_E({\mathbf{u}}^h-\Pi {\mathbf{u}}^h,{\mathbf{v}}^h-\Pi {\mathbf{v}}^h).$$

$$a_E^h(\mathbf{p},{\mathbf{v}}^h)=a_E(\mathbf{p},{\mathbf{v}}^h).$$

$$\forall {\mathbf{v}}^h\in {\mathcal{V}}^h{|}_E,{\alpha }_{\ast }{a}_E({\mathbf{v}}^h,{\mathbf{v}}^h)\le a_E^h({\mathbf{v}}^h,{\mathbf{v}}^h)\le {\alpha }^{\ast }{a}_E({\mathbf{v}}^h,{\mathbf{v}}^h).$$

$$\mathbf{\epsilon }(\mathbf{v})=\nabla \mathbf{v}-\mathbf{\omega }(\mathbf{v}),$$

$$\mathbf{\omega }(\mathbf{v})=\frac{1}{2}\left(\nabla \mathbf{v}-{\nabla }^{\mathsf{T}}\mathbf{v}\right).$$

$$a_E(\mathbf{p},{\mathbf{v}}^h-\Pi {\mathbf{v}}^h)$$

$$=\mathbf{\sigma }(\mathbf{p}):\left[{\int }_E\nabla {\mathbf{v}}^{h}d\mathbf{x}-\nabla \Pi {\mathbf{v}}^h{\int }_{E}d\mathbf{x}\right]$$

$$=\mathbf{\sigma }(\mathbf{p}):\left[{\int }_E\nabla {\mathbf{v}}^{h}d\mathbf{x}-|E|\nabla \Pi {\mathbf{v}}^h\right]$$

$$=0,$$

$$\nabla \Pi {\mathbf{v}}^h=\frac{1}{|E|}{\int }_E\nabla {\mathbf{v}}^{h}d\mathbf{x}.$$

$$\nabla \Pi {\mathbf{v}}^h=\mathbf{\epsilon }(\Pi {\mathbf{v}}^h)+\mathbf{\omega }(\Pi {\mathbf{v}}^h)=\frac{1}{|E|}{\int }_E\mathbf{\epsilon }({\mathbf{v}}^h)d\mathbf{x}+\frac{1}{|E|}{\int }_E\mathbf{\omega }({\mathbf{v}}^h)d\mathbf{x},$$

$$\Pi {\mathbf{v}}^h=\left(\frac{1}{|E|}{\int }_E\mathbf{\epsilon }({\mathbf{v}}^h)d\mathbf{x}\right)\cdot \mathbf{x}+\left(\frac{1}{|E|}{\int }_E\mathbf{\omega }({\mathbf{v}}^h)d\mathbf{x}\right)\cdot \mathbf{x}+{\mathbf{a}}_0.$$

Full text

\NME

10000000

\runningheads

Silva-Valenzuela et al.A nodal integration scheme for meshfree Galerkin methods

\noreceived\norevised\noaccepted

\corraddr

A. Ortiz-Bernardin, Department of Mechanical Engineering, Universidad de Chile, Av. Beauchef 851, Santiago 8370456, Chile. E-mail: aortizb@uchile.cl

1 INTRODUCTION

Meshfree Galerkin methods are based on the weak form where the field variables are discretized using basis functions associated with a set of scattered nodes that partition the domain of analysis. Since the inception of the element-free Galerkin method [1], efficient and accurate numerical integration of the weak form integrals has attracted broad attention in meshfree Galerkin methods. The different numerical integration techniques that exist for meshfree Galerkin methods can be grouped into two main approaches: cell integration and nodal integration techniques. Nodal integration in meshfree methods is attractive since state variables (strains, stresses, and internal variables) can be stored at the nodes, thereby avoiding the need for remapping algorithms at integration points that arise in traditional Lagrangian finite element large deformation simulations with remeshing. In the optimal transportation meshfree method [2], remeshing is avoided by using the integration points as the carriers of material state information in tandem with changes in the support of the basis functions to enable Lagrangian finite deformation simulations. In the material point method [3], integration points of an initial grid (usually, a Cartesian grid) carry the material state variables and the solution stage occurs in three phases: (a) information is mapped from particles to grid nodes; (b) equations of motion are solved on the grid nodes, and the updated information is mapped back to the particles to update their positions and velocities; and (c) the grid is reset. In this process, the mesh remains fixed throughout the simulation, and the deformation causes the material points to move within other cells in the grid. In this paper, we use the virtual element decomposition [4] and follow Reference [5] to devise an accurate and stable nodally integrated meshfree method for linear elastostatics and linear elastodynamics, as well as nonlinear viscoelasticity.

The most simple cell integration technique requires the construction of nonoverlapping cells on which standard Gauss integration is performed. However, this integration scheme is inexact due to the following two properties of meshfree basis functions [6]: (1) They are nonpolynomial functions; and (2) in general, the region that is defined by the intersecting supports of two overlapping nodal basis functions does not coincide with the integration cell. These are two issues that often lead to consistency errors (patch test is not passed) and stability problems due to underintegration. Various approaches have been put forth in meshfree Galerkin methods to address integration errors on cells. For instance, higher-order tensor-product Gauss quadrature is adopted in the element-free Galerkin method [1, 6], whereas the support of the nodal basis functions is used as the domain of integration in the meshless local Petrov-Galerkin method [7] and also in the method of finite spheres [8]. More details on cell-based integration schemes in meshfree Galerkin methods can be found in Ortiz-Bernardin et al. [5].

Nodal integration techniques perform the integration of the weak form integrals by sampling them at the nodes. This approach requires the construction of nodal cells that represent the volume of the integrals being sampled at the nodes. Nodal integration techniques are also prone to integration errors and require to be stabilized. A direct nodal integration (1-point) scheme leads to rank instabilities because meshfree basis functions have zero or nearly zero derivatives at the nodes. In an effort to stabilize the rank instability, Beissel and Belytschko [9] introduced a least-squares residual-based method where the second-order derivatives stabilize the rank instability. With the aim of computing nodal derivatives away from the nodes to avoid rank instabilities, Chen et al. [10] devised a strain smoothing procedure known as stabilized conforming nodal integration (SCNI) approach, which is the cornerstone of various meshfree nodal-based [11, 12, 13] and cell-based [14, 15, 16] integration methods, and even smoothed finite element methods [17]. Another method to compute derivatives away from the nodes is the stress-point method, which was first introduced in the smoothed particle hydrodynamics (SPH) meshfree method [18]. Although the SCNI approach [10] suppresses the rank instability and provides patch test satisfaction, instabilities due to nonzero low-energy modes are still encountered [12]. Puso et al. [12] proposed a penalty stabilization to stabilize the SCNI method. Hillman and Chen [19] combined the arbitrary order variationally consistent meshfree nodal integration framework [13] with a stabilization method devised from an implicit gradient expansion of the strains at the nodes. The resulting method is first-order variationally consistent and stable, and unlike the stabilization method of Puso et al. [12], it is devoid of tunable parameters.

Recently, the virtual element method [4, 20] (VEM) has been proposed, where an exact algebraic construction of the stiffness matrix is realized without the explicit construction of basis functions (basis functions are deemed as virtual). In the VEM, the stiffness matrix is decomposed into two parts: a consistent term that reproduces a given polynomial space and a correction term that provides stability. Such a decomposition (herein referred to as the virtual element decomposition) is formulated in the spirit of the Lax equivalence theorem (consistency $+$ stability $\to$ convergence) for finite-difference schemes and is sufficient for the method to pass the patch test. In polygonal and polyhedral finite elements, Talischi and Paulino [21], Gain et al. [22] and Manzini et al. [23] have used the virtual element decomposition to pass the patch test. The VEM can be viewed as a stabilized Galerkin method on polytopal meshes [24]. For meshfree Galerkin methods with cell-based integration, Ortiz-Bernardin et al. [5] used the virtual element decomposition to develop a method for linear elastostatics that is consistent and stable.

In this paper, on using the virtual element decomposition, a novel meshfree nodal integration technique for linear and nonlinear analyses of solids for small displacements and small-strain kinematics is presented. In particular, linear elastostatics, linear elastodynamics and viscoelasticity are considered. We use the acronym NIVED to refer to this method. The distinctions and contributions of the present work vis-à-vis our prior work on meshfree Galerkin methods using the virtual element decomposition [5] are as follows:

• •

focus herein is on nodal integration, whereas Gauss-like integration scheme is used in Ortiz-Bernardin et al. [5];

• •

we develop the method for linear (static and dynamic) and nonlinear analyses of solids for small displacements and small-strain kinematics, whereas only linear elastostatics is treated in the prior work; and

• •

unlike the prior work, the stabilization of the stiffness matrix in the nodal integration method does not require a tunable (scaling) parameter.

We consider maximum-entropy basis functions (Section 2), although the formulation is applicable to any linear meshfree approximant. The governing equations for linear elastostatics are described in Section 3. The main ingredients of the virtual element framework and the development of our proposed nodal integration technique are presented in Section 4 for linear elasticity and in Section 6 for nonlinear analysis of viscoelastic solids. The demonstration of the patch test satisfaction is provided in Section 5. The accuracy and convergence of the devised nodal integration method are assessed in Section 7 through several examples in linear elastostatics, linear elastodynamics and viscoelasticity. A summary of our main findings with some final remarks is provided in Section 8.

2 MAXIMUM-ENTROPY BASIS FUNCTIONS

Consider a convex domain represented by a set of $N$ scattered nodes and a prior (weight) function $w_{a}(\bm{x})$ associated with each node $a$. We can write down the approximation for a vector-valued function $\bm{v}(\bm{x})$ in the form:

[TABLE]

where $\bm{v}_{a}:=\bm{v}(\bm{x}_{a})$ are nodal coefficients, $\phi_{a}(\bm{x})$ is the meshfree basis function associated with node $a$ and $m\leq N$ represents the number of nodes whose basis functions take a nonzero value at the point $\bm{x}$.

On using the Shannon-Jaynes (or relative) entropy functional, the maximum-entropy basis functions $\{\phi_{a}(\bm{x})\geq 0\}_{a=1}^{m}$ are obtained via the solution of the following convex optimization problem [25]:

[TABLE]

where $\bm{c}_{a}=\bm{x}_{a}-\bm{x}$ are shifted nodal coordinates and ${\rm{I\!R}}_{+}^{m}$ is the nonnegative orthant. In this paper, we use as the prior weight function either the Gaussian radial basis function given by [26]

[TABLE]

where $\gamma$ is a parameter that controls the support size of the basis function and $h_{a}$ is a characteristic nodal spacing associated with node $a$, or the $C^{2}$ quartic polynomials given by [27]

[TABLE]

where $q=\|\bm{c}_{a}\|/(\gamma h_{a})$.

On using the method of Lagrange multipliers, the solution to (2) is given by [25]

[TABLE]

where the Lagrange multiplier vector $\bm{\lambda}(\bm{x})$ is obtained as the minimizer of the dual optimization problem ($\bm{x}$ is fixed):

[TABLE]

where $\bm{\lambda}^{*}$ is the converged solution that gives the basis functions as $\phi_{a}(\bm{x})=\phi_{a}(\bm{x},\bm{\lambda}^{*})$ for $a=1,\ldots,m$.

3 GOVERNING EQUATIONS

3.1 Strong form

Consider an elastic body that occupies the open domain $\Omega\subset{\rm{I\!R}}^{2}$ and is bounded by the one-dimensional surface $\Gamma$ whose unit outward normal is $\bm{n}_{\Gamma}$. The boundary is assumed to admit decompositions $\Gamma=\Gamma_{u}\cup\Gamma_{t}$ and $\emptyset=\Gamma_{u}\cap\Gamma_{t}$, where $\Gamma_{u}$ is the Dirichlet boundary and $\Gamma_{t}$ is the Neumann boundary. The closure of the domain is $\overline{\Omega}=\Omega\cup\Gamma$. Let $\bm{u}(\bm{x}):\overline{\Omega}\rightarrow{\rm{I\!R}}^{2}$ be the displacement field at a point $\bm{x}$ of the elastic body when the body is subjected to external tractions $\bar{\bm{t}}(\bm{x}):\Gamma_{t}\rightarrow{\rm{I\!R}}^{2}$ and body forces $\bm{b}(\bm{x}):\Omega\rightarrow{\rm{I\!R}}^{2}$. The imposed Dirichlet (essential) boundary conditions are $\bar{\bm{u}}(\bm{x}):\Gamma_{u}\rightarrow{\rm{I\!R}}^{2}$. The boundary-value problem for linear elastostatics is: find $\bm{u}(\bm{x}):\overline{\Omega}\rightarrow{\rm{I\!R}}^{2}$ such that

[TABLE]

where $\bm{\sigma}$ is the Cauchy stress tensor.

3.2 Weak form

The Galerkin weak formulation reads: find $\bm{u}(\bm{x})\in\mathscr{V}$ such that

[TABLE]

where $\mathscr{V}$ and $\mathscr{W}$ are the continuous displacement trial and test spaces:

[TABLE]

and $\bm{\varepsilon}(\bm{v})$ is the small strain tensor that is given by

[TABLE]

We consider a set of scattered nodes that discretizes $\bar{\Omega}$. Using these nodes, the domain $\Omega$ is partitioned into nonoverlapping nodal representative polygonal cells. This can be achieved using a Voronoi diagram (Figure 1) or a triangular mesh where the centroids of triangles surrounding a given node are connected to form a polygon (Figure 1). We denote by $E$ a node and its associated cell. The node $E$ has coordinates $\bm{x}_{E}$ in the Cartesian coordinate system and the area of its associated cell is $|E|$.

A node on $\Gamma_{t}$ is denoted by $S$ and its coordinates in the Cartesian coordinate system by $\bm{x}_{S}$. The length of influence of the node $S$ is represented by $|S|$. Figure 2 presents two typical cases of a node located on $\Gamma_{t}$.

The partition formed by all the nodes and their associated cells lying on $\Omega$ is denoted by $\mathcal{T}^{h}$, where $h$ is the maximum diameter of any cell in the partition. The one-dimensional partition formed by all the nodes and their associated lengths of influence lying on $\Gamma_{t}$ is denoted by $\mathcal{E}^{h}$.

Following a standard Galerkin procedure, we define the following discrete local spaces at the nodal cell level:

[TABLE]

where $\bm{u}^{h}(\bm{x})$ is given in the form (1) with $\bm{x}\in E$. These discrete local spaces are assembled to form the following discrete global spaces:

[TABLE]

Owing to the definition of the discrete global spaces, we can evaluate the weak form (6) by sampling it locally at each node in $\mathcal{T}^{h}$ and summing through all of them, as follows:

[TABLE]

where $\mathcal{D}$ is the material moduli tensor that defines the constitutive relation $\bm{\sigma}=\mathcal{D}:\bm{\varepsilon}(\bm{u})$; $\bm{\varepsilon}_{E}:=\bm{\varepsilon}(\bm{x}_{E})$, $\bm{b}_{E}:=\bm{b}(\bm{x}_{E})$ and $\bar{\bm{t}}_{S}:=\bar{\bm{t}}(\bm{x}_{S})$ are nodal quantities. Equation (8) is the result of direct nodal integration (1-point) of the weak form integrals.

4 NODAL INTEGRATION USING THE VIRTUAL ELEMENT DECOMPOSITION

As discussed in Section 1, the direct nodal integration of the bilinear form as given in (8a) is not viable due to stability issues. As a remedy, nodal integration techniques add a stabilization term to the bilinear form [12]. Here, we take a different route and develop a nodal integration scheme for meshfree methods using the mathematical framework of the virtual element method. In this approach, the consistency and stability of the nodally integrated meshfree Galerkin method is ensured by construction.

4.1 Nodal integration cell

The sampling of the weak form at an integration node in the partition $\mathcal{T}^{h}$, where the partition is constructed by means of one of the methods shown in Figure 1, is performed using 1-point Gauss rule over the edges of its nodal cell. It is stressed that in the nodal integration scheme, the integration node and its associated nodal cell are interchangeable. This means that any quantity that is evaluated on the nodal cell is also a quantity evaluated at the node. Considering the Cartesian coordinate system, $\bm{x}_{E}$ are the coordinates of the integration node and $\bm{n}$ is the unit outward normal to the nodal cell’s boundary. A sample nodal integration cell is shown in Figure 3.

4.2 Nodal contribution

Since the NIVED method uses meshfree basis functions for the discretization of the field variables, we have to consider the contributions from the neighboring nodes when performing Gauss integration over the edges of the nodal cell. The global nodal contribution list at an integration point is defined as the global indices of the nodes whose basis functions take a nonzero value at the integration point. We label these global indices from 1 to $m$ and construct a local nodal contribution list with them. We always keep track of the correspondence between the local and global nodal contribution lists as this correspondence is used later in the assembly of the global stiffness matrix and global force vector. The construction of the local nodal contribution list is schematically shown in Figure 4 for the evaluation of meshfree basis functions at an integration point on the edge of a nodal cell.

4.3 Virtual element decomposition

Due to the nonpolynomial nature of the linearly precise maximum-entropy meshfree basis functions, the approximation of the displacement field using these functions contains a linear polynomial part plus some additional nonpolynomial terms. Let $[\mathcal{P}(E)]^{2}$ represent the space of linear displacements over the nodal cell $E$.

Following the standard VEM literature (see for instance, Reference [28]) the following projection operator onto the linear displacement space is defined:

[TABLE]

which allows the splitting of the meshfree approximation of the fields into their linear polynomial part and their nonpolynomial terms, respectively, as follows:

[TABLE]

The projection $\Pi$ is required to satisfy the following orthogonality condition:

[TABLE]

Using (10) and (11), we obtain the following decomposition of the local bilinear form:

[TABLE]

The first term on the right-hand side of (12) is computable as it depends on the linear fields. However, the second term is noncomputable as it depends on the nonpolynomial terms. In the framework of the VEM, the second term is approximated by a bilinear form that can be conveniently computed adopting the form of a stabilization term. We denote this stability bilinear form by $s_{E}$ and rewrite (12) as follows:

[TABLE]

We refer to (13) as the virtual element decomposition.

The virtual element decomposition is endowed with the following crucial properties for establishing the convergence of the VEM [4, 20]:

For all $h$ and for all $E$ in $\mathcal{T}^{h}$

• •

Consistency: For all $\bm{p}\in[\mathcal{P}(E)]^{2}$ and for all $\bm{v}^{h}\in\mathscr{V}^{h}|_{E}$

[TABLE]

• •

Stability: There exists two constants $\alpha_{*}>0$ and $\alpha^{*}>0$, independent of $h$ and of $E$, such that

[TABLE]

In view of the preceding properties, it is straightforward to recognize that the first term on the right-hand side of (13) provides consistency (i.e., ensures patch test satisfaction) and the second term lends stability. The stability property (15) reveals the necessary conditions that $s_{E}$ must possess: it must be symmetric and positive definite on the kernel of $\Pi$ so that property (15) holds without violating (14).

4.4 Projection operator

The explicit form of the projection operator $\Pi$ is obtained from the orthogonality condition (11). Let the strain tensor (7) be written as

[TABLE]

where $\bm{\omega}(\bm{v})$ is the skew-symmetric tensor given by

[TABLE]

Using (16), noting that $\bm{\sigma}(\bm{p})$ and $\bm{\nabla}\Pi\bm{v}^{h}$ are constant fields over $E$ since $\bm{p},\,\Pi\bm{v}^{h}\in[\mathcal{P}(E)]^{2}$, and that $\bm{\sigma}:\bm{\omega}=0$, we can write the orthogonality condition (11) as

[TABLE]

which leads to

[TABLE]

Note that (20) defines $\bm{\nabla}\Pi\bm{v}^{h}$ as the average value of $\bm{\nabla}\bm{v}^{h}$ over the cell $E$. On using (16), (20) can be rewritten as

[TABLE]

which on integrating yields

[TABLE]

To determine $\bm{a}_{0}$ we need a projection operator onto constants $P_{0}:\mathscr{V}^{h}|_{E}\to{\rm{I\!R}}^{2}$ such that

[TABLE]

Given that the field variables computed at the integration node become the representative field variables of the cell, we define the projection operator onto constant as

[TABLE]

where $\bm{x}_{E}$ are the coordinates of the integration node (see Figure 3).

Applying (23) to (22) gives

[TABLE]

On solving for $\bm{a}_{0}$ in (4.4), we obtain

[TABLE]

We introduce the definitions $\bm{\varepsilon}_{E}(\bm{v}^{h}):=\frac{1}{|E|}\int_{E}\bm{\varepsilon}(\bm{v}^{h})\,d\bm{x}$ and $\bm{\omega}_{E}(\bm{v}^{h}):=\frac{1}{|E|}\int_{E}\bm{\omega}(\bm{v}^{h})\,d\bm{x}$ since both are constant tensors over the cell and thus they can be associated with the integration node. Finally, substituting (26) into (22) yields the projection operator onto the linear displacements, as follows:

[TABLE]

In (27), the cell averages $\bm{\varepsilon}_{E}(\bm{v}^{h})$ and $\bm{\omega}_{E}(\bm{v}^{h})$ are evaluated on the boundary of $E$ by invoking the divergence theorem, which gives the following nodal measures:

[TABLE]

and

[TABLE]

respectively.

4.5 Projection matrix

After some algebraic manipulations, (27) can be written as

[TABLE]

where

[TABLE]

[TABLE]

Using meshfree basis functions leads to the following discrete representation of (31):

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and explicitly replacing $\bm{v}^{h}$ in the form (1) into (32) gives

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

where $\phi_{a}^{E}:=\phi_{a}(\bm{x}_{E})$ is the value of the $a$-th basis function at the integration node $E$. To evaluate $q_{ia}$ in (39), a 1-point Gauss rule is used on each edge of the nodal cell (see Figure 3).

Finally, substituting (33) and (37) into (30) yields the following discrete version of the projection operator:

[TABLE]

which defines the projection matrix (the matrix form of $\Pi$) as

[TABLE]

4.6 Nodally integrated stiffness matrix

The nodally integrated local stiffness matrix is obtained by discretizing the virtual element decomposition (13). We start by working on the consistency term (i.e., the first term on the right-hand side of (13)). On using (27) to obtain $\bm{\nabla}\Pi\bm{v}^{h}=\bm{\varepsilon}_{E}(\bm{v}^{h})+\bm{\omega}_{E}(\bm{v}^{h})$ and considering the constitutive relation $\bm{\sigma}(\Pi\bm{u}^{h})=\mathcal{D}:\bm{\varepsilon}_{E}(\bm{u}^{h})$, we can write

[TABLE]

Using the symmetry of tensors $\bm{\varepsilon}_{E}$ and $\mathcal{D}$, (4.6) can be written in Voigt notation as

[TABLE]

where $\bm{D}$ is the constitutive matrix for an isotropic linear elastic material given by

[TABLE]

for plane strain and plane stress conditions, respectively, where $E_{\mathrm{Y}}$ is the Young’s modulus and $\nu$ is the Poisson’s ratio; $\bm{\upvarepsilon}_{E}(\bm{v}^{h})$ is defined in (32).

Now, substituting the discrete version of $\bm{\upvarepsilon}_{E}(\bm{v}^{h})$ as given in (37) into (44) leads to the following discrete local consistency bilinear form:

[TABLE]

where $\bm{d}$ is a column vector similar to $\bm{q}$ given in (38) that contains the nodal coefficients that are associated with the displacement trial functions.

The discrete local stability bilinear form is obtained by replacing the field discretizations $\bm{u}^{h}=\bm{N}\bm{d}$ and $\bm{v}^{h}=\bm{N}\bm{q}$, where $\bm{N}$ is defined in (36), along with (41) into the second term on the right-hand side of (13). This yields

[TABLE]

where $\bm{I}_{2m}$ is the identity ($2m\times 2m$) matrix and $\bm{S}_{E}=s_{E}(\bm{N}^{\mathsf{T}},\bm{N})$. Thus, substituting (46) and (4.6) into the first and second terms on the right-hand side of (13), respectively, gives

[TABLE]

which defines the nodally integrated local stiffness matrix as the sum of the consistency stiffness matrix $\bm{K}_{E}^{\mathrm{c}}$ and the stability stiffness matrix $\bm{K}_{E}^{\mathrm{s}}$, as follows:

[TABLE]

Regarding the stability stiffness, we make the following choice for $\bm{S}_{E}$:

[TABLE]

where $\bm{1}_{2m}$ is the unit ($2m\times 1$) column vector, $\bm{\gamma}=\mathrm{diag}(\bm{K}_{E}^{\mathrm{c}})$ is a column vector containing the diagonal of the matrix $\bm{K}_{E}^{\mathrm{c}}$ and $\odot$ is the element-wise product. This means that $\bm{S}_{E}$ is a diagonal matrix whose diagonal entries are those of the diagonal of the nodally integrated local consistency stiffness matrix, which is in the spirit of the “$D$-recipe” studied in References [29, 30].

4.7 Nodally integrated force vector

For linear fields, the simplest approximation for the body force vector is constructed by projecting both the body force vector $\bm{b}$ and the test displacements $\bm{v}^{h}$ onto constants, as follows [4, 20, 31]:

[TABLE]

where we have used $P_{0}$ as defined in (24), and

[TABLE]

where $\phi_{a}^{E}$ is defined below (40). Hence, the nodal body force vector is given by

[TABLE]

which coincides with the direct integration of the body force vector at the node with coordinates $\bm{x}_{E}$.

The integral that defines the traction force vector is similar to the integral that defines the body force vector, but it is one dimension lower. This means that we can simply apply a direct nodal integration on the Neumann edges. Proceeding likewise leads to the following nodal traction force vector:

[TABLE]

where

[TABLE]

with $\phi_{a}^{S}:=\phi_{a}(\bm{x}_{S})$.

4.8 Nodally integrated mass matrix

With the aim of testing the NIVED approach for an elastodynamic problem, we show how to compute the approximation of the mass matrix using the virtual element decomposition. The mass matrix is constructed along the same lines as the stiffness matrix, that is, the mass matrix is split into a consistency part and a stability part. However, this construction is done using an $L^{2}$-projection on $[\mathcal{P}(E)]^{2}$ of any function $\bm{v}^{h}\in[\mathcal{W}(E)]^{2}$. For linear and quadratic fields, the $L^{2}$-projection coincides with the elliptic projection $\Pi$ [32, 28]. This means,

[TABLE]

For testing our nodal integration approach, we use Newmark’s constant average acceleration time integration method, which is implicit and unconditionally stable. For such integrators, the stability part of the mass matrix is not needed [33]. Thus, generally we only require the consistency term $\int_{E}\Pi\phi_{a}\Pi\phi_{b}\,d\bm{x}$. Hence, the expression for the nodally integrated mass matrix for a solid material of density $\rho$, is:

[TABLE]

where the integral has been nodally integrated at $\bm{x}_{E}$ leading to $\bm{N}_{E}$ as defined in (52).

For a detailed treatment of time-dependent problems using the virtual element framework, the interested reader is referred to References [33, 34, 35].

5 PATCH TEST SATISFACTION

In the patch test, the exact linear field $\bm{u}(\bm{x})=\bm{p}(\bm{x})\in[\mathcal{P}(E)]^{2}$ is imposed on the entire boundary of the domain. Hence, the numerical solution for the patch test must be $\bm{u}^{h}(\bm{x})=\bm{p}(\bm{x})$ on the whole domain. Using this in (13), we get

[TABLE]

The second term on the right-handside of (57) is zero because $\bm{p}-\Pi\bm{p}=\bm{p}-\bm{p}=0$. Then,

[TABLE]

Note that $\bm{\sigma}(\bm{p})$ is a constant field. Let $\bm{\sigma}^{c}$ denote this constant stress field and observe that $\bm{\varepsilon}(\Pi\bm{v}^{h})=\bm{\varepsilon}_{E}(\bm{v}^{h})$ is the nodal strain (cell average) defined in (28) and thus a constant field as well. Hence, (58) becomes

[TABLE]

We now substitute $\bm{v}^{h}(\bm{x})=\phi_{a}(\bm{x})\bm{v}_{a}$ into (59) and sum over all the nodal cells to obtain

[TABLE]

where

[TABLE]

and $\bm{f}_{a}$ is an interior nodal force. Since the patch test produces a state of constant strains and stresses, all the interior nodal forces must be identically equal to zero. Therefore, for the patch test to be satisfied it suffices to show that

[TABLE]

where the assembly is over all the nodal cells that have a non-zero intersection with the support of $\phi_{a}$. To compute (62), we choose a 1-point Gauss rule over the edges of the integration cells, and since the evaluation of $\bm{B}_{a}$ at a given interior edge will arise from two adjacent cells in the assembly, the two contributions cancel each other. Thus, the net contribution to $\bm{f}_{a}$ from all the interior edges vanishes, and hence (62) is satisfied.

6 EXTENSION TO NONLINEAR ANALYSIS: THE VISCOELASTIC CASE

In this section, the proposed NIVED scheme is extended to solid mechanics problems in which the constitutive law for the solid material depends on the history of deformation. In particular, we consider the case of Generalized Maxwell viscoelastic model describing a standard linear isotropic solid [36], wherein the Cauchy stress in vector form and as a function of time $t$ is split as

[TABLE]

where $K$ is the material’s Bulk modulus, $\bm{\upvarepsilon}$ is the small strain tensor in vector form, $\bm{m}=[1\quad 1\quad 0]^{\mathsf{T}}$ and $\bm{s}(t)$ is the deviatoric stress given by

[TABLE]

where $\bm{e}=\bm{\upvarepsilon}-1/3\,\bm{m}\bm{m}^{\mathsf{T}}\bm{\upvarepsilon}$ is the deviatoric strain and $G(t)$ is the shear relaxation modulus given by a one-term Prony series that is of the form

[TABLE]

where $G=E/(2(1+\nu))$, $\mu_{0}>0$ and $\mu_{1}>0$ are parameters that satisfy $\mu_{0}+\mu_{1}=1$, and $\lambda_{1}$ is a relaxation time. Substituting (65) into (64) leads to

[TABLE]

where $\bm{q}^{(1)}(t)$ is the dimensionless partial deviatoric strain given by

[TABLE]

The term $\bm{q}^{(1)}(t)$ is solved using the recursion formula [36]

[TABLE]

where

[TABLE]

Thus, using the recursion formula, the stresses can be written as

[TABLE]

and

[TABLE]

On considering (71) along with (70) and (68), it is realized that $\bm{\upsigma}_{n+1}$ is determined by the knowledge of the time increment $\Delta t$ and the strain $\bm{\upvarepsilon}$ at times $t_{n}$ and $t_{n+1}$.

The numerical solution of the equilibrium equation for the standard linear isotropic solid described by the Generalized Maxwell viscoelastic model is sought incrementally using a Newton-Raphson scheme. Under the NIVED framework, the constitutive equation (71) is computed using the vector form of the nodal strain (28) and the base equilibrium equation for the Newton-Raphson scheme is obtained by linearizing the following VED representation of the virtual work:

[TABLE]

The linearization of (6) yields

[TABLE]

where $\bm{D}_{\mathsf{T}}$ is the tangent moduli that is obtained by deriving (71), as follows:

[TABLE]

Adopting the same argument used to obtain the VED decomposition (13), we approximate the second term on both sides of (6) by computable quantities, thus giving

[TABLE]

We remark that the linearized virtual work (6) is supported by the virtual element framework for inelastic analysis presented in Reference [37]. After discretizing (6), relying on the arbitrariness of nodal variations and summing through all the cells, we finally obtain the following Newton-Raphson iterations 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}$:

[TABLE]

where $\bm{\gamma}=\mathrm{diag}(\bm{K}_{E,\mathsf{T}}^{\mathrm{c}})$, and $\bm{W}_{E}$, $\bm{P}_{E}$, $\bm{I}_{2m}$, $\bm{1}_{2m}$, $\bm{N}_{E}$, $\bm{N}_{S}$, $\bm{b}_{E}$ and $\bar{\bm{t}}_{S}$ are the same as those used in the linear NIVED formulation (see Section 4).

7 NUMERICAL EXPERIMENTS

In this section, numerical experiments are presented to demonstrate the accuracy and convergence of the NIVED approach. Results are compared with a meshfree Galerkin method that uses maximum-entropy basis functions and standard Gauss quadrature with integration points defined in the interior of the cells of a background mesh of 3-node triangles. We denote this method by the acronym MEM. When available, we also provide comparisons of our results to reference solutions obtained with the finite element method (FEM). The construction of the nodal representative polygonal cells for the NIVED method is carried out using the same background mesh of 3-node triangles used for the MEM method, that is, the construction method depicted in Figure 1(b) is considered. This allows a direct comparison of the NIVED and MEM methods since the number of degrees of freedom remains the same in both methods.

For numerical integration in the NIVED approach, we use a 1-point Gauss rule per edge of the polygonal cell and at these points only the evaluation of basis functions is required. In the MEM approach, we consider 1-point, 3-point, 6-point and 12-point Gauss rules in the interior of the 3-node triangular cell. In contrast to the NIVED method, at these interior integration points the evaluation of both the basis functions and their derivatives is required. The maximum-entropy basis functions in the NIVED and MEM methods are performed using $\gamma=2.0$ in all the numerical experiments that follow.

7.1 Patch test

We solve the boundary-value problem (5) for the patch test that is schematically depicted in Figure 5. Plane stress condition is assumed with the following material parameters: $E_{\mathrm{Y}}=3\times 10^{7}$ psi and $\nu=0.3$. The exact solution for this problem is the linear field \bm{u}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}\nu(1-x_{1})/E_{\mathrm{Y}}&x_{2}/E_{\mathrm{Y}}\\ \end{array}\right]^{\mathsf{T}}. The Dirichlet boundary conditions are imposed as follows: on the bottom boundary, the vertical displacement is fixed, and at its right corner, the displacement is additionally fixed in the horizontal direction. The Neumann boundary conditions are \bar{\bm{t}}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}0&0\\ \end{array}\right]^{\mathsf{T}} on the left and right boundaries, and \bar{\bm{t}}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}0&\sigma\\ \end{array}\right]^{\mathsf{T}} on the top boundary. The Gaussian prior weight function is used for the maximum-entropy basis functions. The background meshes used in this study are depicted in Figure 6. Numerical results for the relative error in the $L^{2}$ norm and the $H^{1}$ seminorm are presented in Tables 1 and 2, respectively, for the MEM and NIVED approaches. Numerical results confirm that the patch test is met to machine precision only for the NIVED method.

7.2 Numerical stability

To assess the stability of the NIVED method, eigenvalue analyses are performed on a unit square domain. As a reference for comparison, we include the MEM approach in the analyses. The quartic polynomials and the Gaussian radial basis function are considered as the prior weight function in the evaluation of the maximum-entropy basis functions in these analyses. Three zero eigenvalues, which correspond to the three zero-energy rigid-body modes, are obtained for both methods. The three mode shapes that follow the three rigid-body modes are depicted in Figure 7 for the MEM and NIVED methods using the quartic polynomials as the prior weight function. For the quartic prior, instabilities are observed for the MEM with a 1-point Gauss rule since it exhibits nonsmooth mode shapes (Figures 7–7) and even a 12-point Gauss rule is not sufficient for removing the instabilities (Figures 7–7). In stark contrast, Figures 7–7 show the smooth mode shapes that are obtained in the NIVED approach when using the quartic prior.

Similarly, the three mode shapes that follow the three rigid-body modes are depicted in Figure 8 for the MEM and NIVED methods using the Gaussian radial basis function as the prior weight function. For the Gaussian prior, instabilities are observed for the MEM with a 1-point Gauss rule since it exhibits nonsmooth mode shapes (Figures 8–8), but a 6-point Gauss rule can effectively remove the instabilities (Figures 8–8). The NIVED approach using the Gaussian prior is also free of instabilities as revealed by the smooth mode shapes depicted in Figures 8–8.

7.3 Cantilever beam

We conduct a convergence study for the problem of a cantilever beam of unit thickness subjected to a parabolic end load $P=-1000$ lbf. A schematic representation of the problem is shown in Figure 9. Plane strain condition is assumed with material parameters given by $E_{\mathrm{Y}}=10^{7}$ psi and $\nu=0.3$. The essential boundary conditions on the clamped edge are applied according to the analytical solution given by Timoshenko and Goodier [38]:

[TABLE]

where $\bar{E}_{\mathrm{Y}}=E_{\mathrm{Y}}/\left(1-\nu^{2}\right)$ and $\bar{\nu}=\nu/\left(1-\nu\right)$; $L=8$ inch is the length of the beam, $D=4$ inch is the height of the beam, and $I=D^{3}/12$ is the second-area moment of the beam section. The exact stress field is:

[TABLE]

The Neumann boundary conditions are applied using the exact stress field, which gives \bar{\bm{t}}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}0&0\\ \end{array}\right]^{\mathsf{T}} on the top and bottom boundaries, and \bar{\bm{t}}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}0&\sigma_{12}\\ \end{array}\right]^{\mathsf{T}} on the right boundary. For the evaluation of the maximum-entropy basis functions, the Gaussian prior is used. The sequence of background integration meshes used in the study are shown in Figure 10.

The convergence of the MEM and NIVED methods in the $L^{2}$ norm and the $H^{1}$ seminorm of the error are shown in Figure 11. The MEM approach needs a 3-point Gauss rule to deliver the optimal rate in the $L^{2}$ norm. Regarding the $H^{1}$ seminorm, the convergence of the MEM method behaves erratically for 1-point and 3-point Gauss rules due to integration errors, and a 6-point Gauss rule is needed to recover the optimal rate. In contrast, the proposed NIVED method delivers the optimal rate of convergence in both the $L^{2}$ norm and the $H^{1}$ seminorm of the error.

7.4 Infinite plate with a circular hole

The rates of convergence are studied for the problem of an infinite plate with a circular hole that is loaded at infinity according to the tractions $\sigma_{11}=T$ and $\sigma_{22}=\sigma_{12}=0$ (Figure 12). Due to the symmetry of the geometry and boundary conditions, we consider the domain of analysis shown in Figure 12. Plane stress condition is assumed with $E_{\mathrm{Y}}=10^{3}$ psi and $\nu=0.3$. The exact solution on the domain of analysis is given by [38]

[TABLE]

where $G=E_{\mathrm{Y}}/(2(1+\nu))$ and $\kappa=(3-\nu)/(1+\nu)$. The exact stress field is:

[TABLE]

where $r$ is the radial distance from the center $(x_{1}=0,x_{2}=0)$ to a point $(x_{1},x_{2})$ in the domain of analysis. In the computations, the following data are used: $T=100$ psi, $r_{0}=1$ inch and $a=5$ inch. The Dirichlet boundary conditions on the domain of analysis are imposed as follows: $\bar{u}_{1}=0$ on the left side and $\bar{u}_{2}=0$ on the bottom side. The Neumann boundary conditions are prescribed using the exact stresses, as follows: \bar{\bm{t}}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}\sigma_{12}&\sigma_{22}\\ \end{array}\right]^{\mathsf{T}} on the top side and \bar{\bm{t}}=\left[\begin{array}[]{ccccccccccccccccccccccccccccccccccccccccccccccccccc}\sigma_{11}&\sigma_{12}\\ \end{array}\right]^{\mathsf{T}} on the right side. The sequence of background integration meshes used in the study are shown in Figure 13. The Gaussian prior is used for the evaluation of the maximum-entropy basis functions.

The convergence rates that are delivered by the MEM and NIVED approaches in the $L^{2}$ norm and the $H^{1}$ seminorm of the error are compared in Figure 14. We observe that the convergence of the MEM using a 1-point Gauss rule behaves erratically in the $L^{2}$ norm and that the optimal rate of 2 is recovered using a 3-point Gauss rule. Regarding its convergence in the $H^{1}$ seminorm, we observe that it is optimal when using at least a 6-point Gauss rule. On the other hand, the NIVED approach delivers the optimal rates of 2 and 1 in the $L^{2}$ norm and the $H^{1}$ seminorm, respectively.

7.5 $L$-shaped domain under traction load

An $L$-shaped domain under traction load is considered. The geometry and boundary conditions are shown in Figure 15, where $H=100$ inch and $p=1$ psi. The thickness is 1 inch. A stress singularity occurs at the re-entrant corner and the exact solution is not known. However, an estimation of the exact strain energy using Richardson’s extrapolation is known and has the value 15566.46 [39]. We use this value to assess the convergence of the strain energy. Plane stress condition is assumed with $E_{\mathrm{Y}}=1$ psi and $\nu=0.3$. The sequence of background meshes used in this study are depicted in Figure 16. The convergence of the strain energy for the MEM with several Gauss integration rules and for the NIVED approach is shown in Figure 17. We observe that even with a 12-point Gauss rule, the MEM does not uniformly approach the reference value of the strain energy. In contrast, the proposed NIVED scheme uniformly approaches this reference value.

7.6 Manufactured elastostatic problem

In this example, the manufactured elastostatic problem found in Reference [15] is used to assess the performance of the NIVED method. The analysis is conducted on a $2\times 2$ square domain. Plane stress condition is considered with $E_{\mathrm{Y}}=10^{5}$ psi and $\nu=0.3$. The entire domain boundary is prescribed with the following Dirichlet boundary conditions:

[TABLE]

which correspond to the exact solutions of the linear elastostatic problem (5) manufactured with a body force given by

[TABLE]

The exact stress field is

[TABLE]

Figure 18 depicts the background integration meshes used in the study. The Gaussian prior is used for the evaluation of the maximum-entropy basis functions.

The convergence in the $L^{2}$ norm and the $H^{1}$ seminorm for the MEM and NIVED methods are compared in Figure 19. The $L^{2}$ norm is not convergent for the MEM approach using a 1-point Gauss rule, but the optimal rate of 2 is recovered using a 3-point Gauss rule. The convergence of the MEM method in the $H^{1}$ seminorm is not good for 1-, 3- and 6-point Gauss rules, but the optimal rate of 1 is recovered when using a 12-point Gauss rule. The plots show that the NIVED approach delivers the optimal rates of 2 and 1 in the $L^{2}$ norm and the $H^{1}$ seminorm, respectively.

The computational cost of the MEM and NIVED methods are compared in Figure 20. The methods are assessed in terms of accuracy and cell refinement using the normalized CPU time, which is defined as the ratio of the CPU time of a particular model analyzed to the maximum CPU time of any of the models analyzed. In terms of accuracy, we observe that the computational cost of the NIVED scheme is similar to the computational cost of the MEM method with 12-point Gauss rule (Figure 20). This occurs because for a given level of accuracy, the MEM requires a coarser mesh than the NIVED scheme, as shown in Figure 19. In terms of cell refinement, we observe that the computational cost of the NIVED method is similar to the computational cost of the MEM method with a 3-point Gauss rule (Figure 20). However, the MEM with a 3-point Gauss rule does not converge in the $H^{1}$ seminorm, which means that the computational cost of a convergent MEM approach is greater than the computational cost of the NIVED method. In fact, from Figure 19 we know that a 12-point Gauss rule is needed for optimal $H^{1}$-convergence of the MEM method, but this rule is the most expensive as shown in Figure 20. Combining both the cost in terms of accuracy and cell refinement, we conclude that in general the NIVED scheme can be expected to be faster but less accurate than a convergent MEM approach with the same number of degrees of freedom.

7.7 Manufactured elastodynamic problem

A manufactured elastodynamic problem that is found in Reference [40] is considered. The domain of analysis is a $2\times 2$ square domain and the background meshes considered are the same used for the manufactured elastostatic problem (Figure 18). The material parameters are $E_{\mathrm{Y}}=1\times 10^{5}$ psi, $\nu=0.3$ and $\rho=800$ lb/$\mathrm{in}^{3}$, and plane stress condition is assumed. The Gaussian prior is used for the evaluation of the maximum-entropy basis functions.

The problem is manufactured with the following body force:

[TABLE]

where $L$ is the side length of the domain of analysis and

[TABLE]

The exact displacement field solution is

[TABLE]

and the exact stress field is

[TABLE]

At the initial condition $t=0$ seconds, the body is at rest. The exact displacement field solution is used to impose the Dirichlet boundary conditions along the entire boundary of the domain. The following values for the parameters are used in the computations: $\alpha=0.001$, $\beta=0.001$. The Newmark’s constant average acceleration method with a time increment $\Delta t=0.01$ seconds is used as the time integration algorithm.

The convergence of the MEM and NIVED methods in the $L^{2}$ norm and the $H^{1}$ seminorm after 100 time steps (i.e., $t=1$ s) are compared in Figure 21. The MEM approach exhibits suboptimal convergence in the $L^{2}$ norm with a 1-point Gauss rule, but the optimal rate of 2 is recovered using a 3-point Gauss rule. Its convergence in the $H^{1}$ seminorm is optimal when using at least a 12-point Gauss rule. The optimal rates of 2 and 1 are obtained in the $L^{2}$ norm and the $H^{1}$ seminorm, respectively, for the NIVED method.

Similar to the static case, Figure 22 (also obtained at $t=1$ s) reveals that in terms of cell refinement, the computational cost of a convergent MEM approach is greater than the NIVED counterpart thereby making the NIVED approach superior in performance. Moreover, in the dynamic regime the outperformance of the NIVED over the MEM, in terms of cell refinement, is more pronounced as the computational cost of the NIVED approach is about the same as the computational cost of the MEM using a 1-point Gauss rule (Figure 22).

7.8 Thick-walled viscoelastic cylinder subjected to internal pressure

We conclude this section with a demonstration of the nonlinear NIVED formulation for small-displacement and small-strain analysis. The test problem consists of a thick-walled cylinder subjected to internal pressure whose material constitution is a linear isotropic solid described by the Generalized Maxwell viscoelastic model considered in Section 6. This example has been used to assess the virtual element method for inelastic analysis in Reference [37]. Due to symmetry, only a quarter of the cross section of the cylinder is considered, as shown in Figure 23, where $r_{i}=2$ inch, $r_{0}=4$ inch and $p=10$ psi. The material parameters are set to $E=1000$ psi, $\nu=0.3$ and $\lambda_{1}=1$. The pressure is applied suddenly and the structural response is computed through 20 unit time steps.

In the first part of this study, the convergence of the radial displacement at control points $A$ and $B$ (Figure 23) is assessed. For comparison purposes, the numerical solution obtained with a highly refined mesh of 9-node FEM quadrilateral elements (FEM-Q9) is used. Three sets of choices are considered for the material parameters $\mu_{0}$ and $\mu_{1}$: $(\mu_{0},\mu_{1})=\{(0.7,0.3),\,(0.3,0.7),\,(0.01,0.99)\}$. Using $E$ and $\nu$ to calculate the material’s bulk modulus $K$ and (65) to calculate $G(t=20)$, the effective Poisson’s ratio is calculated as $\nu_{\mathrm{eff}}=\frac{3K-2G(20)}{2(3K+G(20))}$. For the three preceding sets of parameters this gives at time $t=20$ the values $\nu_{\mathrm{eff}}=\{0.3542,\,0.4338,\,0.4977\}$, respectively. Thus, for the last set, the material is near-incompressible.

The convergence of the NIVED radial displacement at control points $A$ and $B$ is summarized in Table 3 for the sequence of background integration meshes depicted in Figure 24. Good agreement with the FEM-Q9 reference solution is observed for the three sets of material parameters.

Curves depicting the time history of the radial displacement measured at control points $A$ and $B$ are shown in Figure 25 for the three sets of material parameters. The response curves for the NIVED scheme is plotted for the mesh shown in Figure 24. Good agreement with the FEM-Q9 reference solutions is observed.

Lastly, the performance of the NIVED method for near-incompressible material behavior is studied. The NIVED approach is tested on the background mesh shown in Figure 26 and its radial stress field solution is compared with the corresponding FEM-Q9 solution on the mesh shown in Figure 26. Also on this mesh, the radial stress field is computed with 4-node FEM quadrilateral elements (FEM-Q4). The comparisons are summarized in Figure 27. The radial stress field corresponding to the compressible case $(\mu_{0},\mu_{1})=(0.7,0.3)$ ($\nu_{\mathrm{eff}}=0.3542$) is shown in Figures 27–27 for FEM-Q4, FEM-Q9 and NIVED methods, respectively. It is observed that the three methods perform well in this case. Of course, because of its higher approximation order, the FEM-Q9 radial stress field is smoother than the FEM-Q4 and NIVED radial stress fields.

The radial stress field corresponding to the near-incompressible case $(\mu_{0},\mu_{1})=(0.01,0.99)$ ($\nu_{\mathrm{eff}}=0.4977$) is depicted in Figures 27–27 for FEM-Q4, FEM-Q9 and NIVED methods, respectively. In this case, the FEM-Q9 radial stress field looks smooth and no sign of locking is observed. However, locking behavior is exhibited in the FEM-Q4 case as evidenced by the oscillatory radial stress field. The NIVED scheme performs better than the FEM-Q4 in the incompressible limit, but some small oscillations do occur in the radial stress field and thus it is not completely free of locking. This is not surprising since use of the virtual element framework in the NIVED approach ensures consistency (patch test satisfaction) and stability for a displacement-based formulation. To render the NIVED scheme locking-free in the near-incompressible limit would require introducing, for instance, a form of a nodal averaging technique [41, 42].

8 CONCLUDING REMARKS

In this paper, we proposed a novel nodal integration scheme for meshfree Galerkin methods that is devised from the virtual element decomposition [4], wherein the bilinear form is decomposed into a consistency part and a stability part that ensure that the method is consistent and stable. Linear maximum-entropy meshfree basis functions were adopted, but the formulation is applicable to any other linear meshfree approximant. We referred to this new nodal integration scheme as NIVED. As in any nodal integration method, a nodal representative cell is needed in the NIVED approach. The nodal cell can be constructed from a Voronoi diagram or a Delaunay triangulation by connecting the centroids of triangles surrounding a node.

Numerical tests in linear elastostatic and linear elastodynamic boundary-value problems were conducted, where the focus was on comparing the performance of the NIVED method and the maximum-entropy meshfree method (MEM) using standard Gauss integration. We also provided the extension of the NIVED approach to nonlinear analysis where the material constitution is a linear isotropic solid described by the Generalized Maxwell viscoelastic model. Comparisons with reference solutions obtained with 4-node and 9-node FEM quadrilateral elements were conducted. In the NIVED approach, the weak form integrals are integrated by sampling them at the nodes and the integration at a node is performed using a 1-point Gauss rule over the edges of its nodal representative cell. This task only involves the evaluation of basis functions (no derivatives are needed). In the MEM approach, the integration of the weak form integrals is performed using standard Gauss integration on the interior of triangular cells, which requires the evaluation of basis functions derivatives.

Our main findings from this work are as follows. The NIVED scheme passes the patch test to machine precision, whereas the MEM does not achieve this level of accuracy due to integration errors. The numerical stability test showed that both the NIVED and MEM methods deliver the three zero-energy rigid-body modes, but instability is exhibited only by the MEM method as evidenced by the presence of the nonsmooth eigenmodes that follow the three rigid-body modes. Convergence was assessed through several numerical experiments, which included a cantilever beam subjected to a parabolic end load, an infinite plate with a circular hole and manufactured (elastostatic and elastodynamic) problems. The convergence studies revealed that the NIVED method delivers the optimal rate of convergence in the $L^{2}$ norm and the $H^{1}$ seminorm. On the other hand, the convergence of the MEM was dependent on the number of Gauss points used inside the triangular integration cell. In the tests that were conducted, the MEM required a 3-point Gauss rule for optimal convergence in the $L^{2}$ norm and, depending on the problem, a 6-point or a 12-point Gauss rule for optimal convergence in the $H^{1}$ seminorm. In terms of computational cost, it was shown that for the same number of degrees of freedom the proposed NIVED scheme is faster but less accurate than a convergent MEM approach. We also tested the NIVED and MEM methods in an $L$-shaped domain under traction load, where the presence of a re-entrant corner introduces a stress singularity that results in a nonsmooth solution. In this case, a study to assess the convergence to a reference value of the strain energy was conducted. The results showed that only the NIVED scheme uniformly approaches the reference strain energy. The performance of the nonlinear NIVED formulation was assessed by solving the problem of a thick-walled viscoelastic cylinder subjected to internal pressure. The radial displacement solution was found to be in good agreement with a reference solution obtained using quadratic (9-node) quadrilateral finite elements. We also showed that for a near-incompressible viscoelastic material, the NIVED scheme is not completely devoid of locking.

In closing, we mention that a desirable feature that can be offered by the (nodally integrated) NIVED approach is that state variables such as strains, stresses and other internal variables in nonlinear computations can be stored at the nodes, which is attractive for Lagrangian large deformation meshfree simulations. This extension of the NIVED method to large deformations is appealing, and is planned as part of future work.

ACKNOWLEDGEMENTS

RSV and AOB acknowledge the support provided by the Chilean National Fund for Scientific and Technological Development (FONDECYT) through grant CONICYT/FONDECYT No. 1181192. NS acknowledges the research support of Sandia National Laboratories. EA gratefully acknowledges the partial financial support of the University of Rome Tor Vergata Mission Sustainability Programme through project SPY-E81I18000540005; and the partial financial support of PRIN 2017 project “3D PRINTING: A BRIDGE TO THE FUTURE (3DP_Future). Computational methods, innovative applications, experimental validations of new materials and technologies,” grant 2017L7X3CS_004. NHK is grateful for the support provided by the Chilean National Fund for Scientific and Technological Development (FONDECYT) through grant CONICYT/FONDECYT No. 1181506.