Elastoplastic Analysis of Plates with Radial Point Interpolation Meshless Methods

Abstract

For both linear and nonlinear analysis, finite element method (FEM) software packages, whether commercial or in-house, have contributed significantly to ease the analysis of simple and complex structures with various working conditions. However, the literature offers other discretization techniques equally accurate, which show a higher meshing flexibility, such as meshless methods. Thus, in this work, the radial point interpolation meshless method (RPIM) is used to obtain the required variable fields for a nonlinear elastostatic analysis. This work focuses its attention on the nonlinear analysis of two benchmark plate-bending problems. The plate is analysed as a 3D solid and, in order to obtain the nonlinear solution, modified versions of the Newton–Raphson method are revisited and applied. The material elastoplastic behaviour is predicted assuming the von Mises yield surface and isotropic hardening. The nonlinear algorithm is discussed in detail. The analysis of the two benchmark plate examples allows us to understand that the RPIM version explored is accurate and allows to achieve smooth variable fields, being a solid alternative to FEM.

1. Introduction

Computational mechanics is the discipline concerned with the use of computational numerical methods to study phenomena governed by the principles of solid and fluid mechanics. It is used every day to help predict, develop, and improve not only structural and fluid applications in the industrial sense but also several other science areas, such as medicine, biology, chemistry, etc. [1]. With the constant evolution of engineering, the challenges for computational mechanics are numerous. For example, when analysing the simulation of manufacturing processes, such as extrusion or moulding, it is necessary to deal with large displacements of mesh, viscoplastic effects and contact conditions, whereas in simulations of failure processes, it is required to model the fracture initiation and its propagation path, leading to the development of demanding remeshing algorithms and fracture-damage criteria.

Problems dealing with large mesh distortions or localised remeshing requirements are not efficiently solved by mesh-dependent discretization techniques, such as the finite element method (FEM) [2]. Meshless methods are an alternative to classic mesh-dependent discretization methods. In meshless methods, the problem domain is discretized by a random set of nodes, without any preestablished connectivity, rather than a structured element mesh. To define the nodal connectivity, meshless methods use the influence-domain concept, allowing one to establish automatically the nodal connectivity in a later stage of the preprocessing phase [3,4].

Meshless methods and their origin can be traced to the smooth particle hydrodynamics method (SPH) [5]. It was applied to polytropic stellar models and used statistical techniques to recover analytical expressions for the physical variables from a known distribution of fluid elements, starting with a nonaxisymmetric distribution of approximately 80 particles in three dimensions. This method was applied to solid mechanics by Larry D. Libersky and A.G. Petschek [6] in 1991, in which an elastic-perfectly-plastic constitutive model was formulated. It was found that excellent features of the SPH, such as robustness, simplicity, and ease of adding new physics, could be transferred to the analysis of solid mechanics problems [7].

SPH is a particle method, in which the set of nodes (the particles) do not form a solid continuum in the same sense as in FEM, being more suited for free surface fluid analyses. Thus, other meshless technique started to be developed, such as the diffuse element method (DEM) [8], the first mature meshless method using a weak formulation and shape functions constructed by using the moving least-squares (MLS) approximants approach. Then, Belytschko and coworkers improved the DEM and proposed in 1994 the element free Galerkin method (EFGM) [9], one of the most popular meshless methods. At the same period, other approximation meshless techniques were developed, such as the reproducing kernel particle method (RKPM) [10] and the meshless local Petrov–Galerkin method (MLPG) [11], a truly meshless method (in the sense that it does not require a background nodal independent integration mesh) based on a local symmetric weak form (LSWF). MLPG does not need an element mesh for the interpolation of the solution variables. Integrals are solved over a regular-shaped domain and its boundaries. The generalized finite difference method (GFDM) is another popular approximation method [12,13]. Because it can be applied to solve problems discretized by using only an irregular nodal distribution, GFDM is also considered truly meshless. GFDM uses the Taylor series expansions and the MLS approximation to derive explicit expressions for the required partial derivatives of unknown variables. Collocation meshless methods are another class of meshless methods available in the literature, such as the local radial basis function collocation method (LRBFCM) [14,15], which uses the strong-form formulation to establish the system of equations and obtain the variable field. LRBFCM uses radial basis functions to construct trial functions, and it is highly dependent on the shape parameters of the corresponding radial basis function. The accuracy and efficiency of LRBFCM strongly depend on the values of such shape parameters, which must be analysed and optimised for distinct problems [14,15]. The literature offers several variants of collocation methods, such as the Trefftz collocation method [16,17].

Approximation methods, such as EFGM, RKPM, and MLPG, produce highly accurate solutions. However, their shape functions lack the Kronecker delta property, hindering the imposition of the essential and natural boundary conditions, which can only be enforced with exact technique, such Lagrange multipliers [1]. Therefore, efforts were made to develop interpolating meshless techniques. The natural element method (NEM), was one of the first to be developed [18]. Using the natural neighbour and Voronoï diagram concepts, the nodal connectivity, the background integration mesh and the shape functions (and corresponding spatial derivatives), are constructed automatically, producing a true meshless method approach with interpolating properties. In 2001, G.R. Liu and his research team developed the point interpolator method (PIM) [19], in which the shape functions are constructed by using a polynomial basis. Although simple and easy to apply, the process is not capable of producing shape functions when nodes are perfectly aligned (most common in regular nodal distributions). Therefore, a radial basis function was added to the formulation, allowing for the development of one of the most popular interpolating meshless methods: the radial point interpolation method (RPIM) [20]. The use of the radial-based functions eliminated the singularities occurring from the nodal alignment. Later, Belinha and coworkers included in the formulation the natural neighbour concept and developed the natural neighbour radial point interpolation method (NNRPIM) [1,21]. With the NNRPIM, both the nodal connectivity and the node-dependent integration background mesh are constructed resorting to the Voronoï tessellation and to the Delaunay triangulation. This procedure guarantees that the shape functions developed have the Kronecker delta property and a formulation that is truly meshless. The literature shows that both RPIM and NNRPIM are highly accurate and present high convergence rates [1], being solid alternatives for structural analysis.

It is possible to find in the literature several works successfully combining meshless methods with elastoplastic formulations. In 1999, Barry and Saigal [22] extended the EFGM to elastoplasticity considering a 3D formulation with small strains, analysing a thick-walled pressure vessel with internal pressure and an extension of a strip with a circular. Compared with reference solutions [23], the method was capable of producing results very close to the ones obtained in the FEM. In 2006, Belinha and Dinis [24,25,26] used the EFGM in the analysis of two-dimensional problems, plates, and laminates in the elastoplastic regime. Various modified Newton–Raphson algorithms were used to obtain the nonlinear solution. A Cook membrane and a sandwich material, under distributed load, were studied. The obtained results were in accordance with FEM solutions. Several modified Newton–Raphson algorithm versions were used to obtain the nonlinear solution and their performance and computational cost was also analysed. The results showed that the EFGM is a solid alternative to such nonlinear problems, being capable of obtaining similar results to the FEM analysis. More recently, in 2014, Cheng et al. [27], developed a novel interpolating EFGM (IEFGM) for two-dimensional elastoplasticity. By using the interpolating moving least-squares method (IMLS) to construct the shape functions, the Kronecker delta property is satisfied, and the essential boundary conditions can be applied directly. The method has some advantages, such as simpler formulae and the direct application of the boundary conditions. The results showed that the IEFGM is in good agreement with FEM, and it was also found that the IEFGM possesses a higher precision than the classic EFGM.

Regarding the RPIM, and its versios, in 2011, Dinis and coworkers [28] extended the NNRPIM to the large-deformation elastoplastic analysis, using the Newton–Rapson nonlinear solution algorithm and an efficient “forward-Euler” procedure to return the stress state to the yield surface. Then, in 2015, Zhang et al. [29] were capable of combining the node-based smoothed radial point interpolation method (NS-RPIM) with gradient-dependent plasticity to analyse the elastoplastic material behaviour of 2D solids. More recently, in 2019, Belinha and coworkers [30] extended both RPIM and NNRPIM to the elastoplastic analysis of aluminium alloys. In order to obtain comparable results, due to the lack of a suitable elastoplastic solution, an experimental test was performed in order to obtain the force/displacement curve and the corresponding mechanical properties. A dog-bone-shaped specimen was analysed with elastoplastic assumptions, under uniaxial tensile state conditions. Two different influence-cell degrees were considered—first degree and second degree—with the results showing that the second-degree influence cell returned a smooth result, whereas the first-degree cell had some discontinuities. Moreover, the results obtained with either the RPIM and NNRPIM showed a strong correlation with the experimental data. The main disadvantage of NNRPIM is clearly addressed: its higher computational cost when compared to the FEM and RPIM, which is due to the time spent establishing the natural neighbours and constructing the background integration points. Nevertheless, being a truly meshless method is also pointed out as its main advantage. Therefore, to fully discretize the problem domain, only a nodal discretization is required. The technique is then capable of automatically generating all the other required numerical structures (nodal connectivity, integration point, and shape functions).

In the present work, a radial point interpolation method version combined with a nonlinear solution algorithm was originally written and coded to study the elastoplastic bending behaviour of plates by using a classic 3D formulation. Thus, in the following sections, the RPI formulation and its extension to elastoplasticity are presented, and then 3D plate examples are solved by considering elastoplastic materials.

2. Meshless Methods

The majority of meshless methods follow a generic procedure, well discribed in the literature [1,2]. After the problem geometry and the natural and essential boundary conditions are described and established, the problem domain and the respective boundary are discretized by using a nodal set, which can follow a regular or irregular distribution. Similar to mesh-dependent methods combined with weak formulations and applied to elasticity problems, the accuracy of the solution provided by meshless methods depends on the nodal density of the discretization. Thus, a higher density of nodes generally returns more accurate results. Moreover, in specific problems wherein particular geometrical or loading details can have a considerable effect on the stress distribution (material discontinuities or concentrated loads), a higher nodal density in those locations is required to obtain better results. In meshless methods, the nodal distribution does not form a mesh because there is no preestablished connectivity between the nodes.

After the nodal discretization of the problem domain, the background integration points are defined. The complete set of integration points is known as integration mesh. Background integration meshes can be obtained by using background integration cells (fitted to the problem domain, as elements in the FEM) or by using nodal integration schemes (as the ones used in NEM or NNRPIM). This last procedure uses the mathematical concepts of natural neighbours, Voronoï diagrams and Delaunay tessellation to construct the integration points uniquely based on the nodal distribution, and thus producing truly meshless formulations.

After discretizing the problem domain with nodes and integration points, it is possible to establish the nodal connectivity. In FEM, as input, it is necessary to declare the nodes forming each element. Thus, integration points inside each element interact with the nodes forming that element, and elements interact with other elements possessing common nodes. In meshless methods, because there is no preestablished declared relation between the nodes, it is necessary to enforce such connectivity. Therefore, it is necessary to apply the concept of influence domain. Each integration point ${\mathbf{x}}_I$ will search for nodes inside a certain radius. The nodes inside such radius will form the influence domain of ${\mathbf{x}}_I$. Due to its simplicity, several meshless formulations apply such concepts (RPIM, EFGM, MLPG, RKPM, etc.). There are several variants to this concept [1], some authors use fixed radial searches and other use variable radial searches. Fixing a radius will produce uneven influence domains. Integration points near a physical boundary will possess a lower number of nodes inside their influence domains when compared with the influence domains of the nodes inside the problem domain. Alternatively, fixing the number of nodes inside each influence domain (variable radius searches), will produce influence domains of the same size for all integration points. Alternatively, meshless formulations using the natural neighbour concept, such as NEM or NNRPIM, already possess the nodal connectivity established due to the Voronoï diagram [31] of the nodal discretization. Each integration point ${\mathbf{x}}_I$ has an associated node ${\mathbf{x}}_i$, and each node possesses a set of natural neighbours: ${\mathbf{V}}_i$. Thus, the influence domain of ${\mathbf{x}}_I$ is the set of natural neighbours of ${\mathbf{x}}_i$: ${\mathbf{V}}_i$, and ${\mathbf{x}}_i$ itself.

Afterward, the shape functions are constructed by using approximation or interpolation functions according to the technique adopted. Regardless of the technique, the field variable can be approximated (or interpolated) at an integration point ${\mathbf{x}}_I$ with

$$u\left({\mathbf{x}}_I\right)=\sum _{i=1}^n{\phi }_i\left({\mathbf{x}}_I\right)·u\left({\mathbf{x}}_i\right)=\mathrm{\Phi }{\left({\mathbf{x}}_I\right)}^T·{\mathbf{u}}_s,$$

(1)

where n is the number of nodes inside the influence domain of the integration point ${\mathbf{x}}_I$, ${\mathbf{u}}_s$ is a vector containing the field variable components of each node ${\mathbf{x}}_i$ included in the influence domain of ${\mathbf{x}}_I$, and ${\phi }_i\left({\mathbf{x}}_I\right)$ is the ith component of the shape function constructed for ${\mathbf{x}}_I$. Instead, the sum operation, the vectorized operation of Equation (1) can be used, in which $\mathrm{\Phi }\left({\mathbf{x}}_I\right)={\{{\phi }_1\left({\mathbf{x}}_I\right),{\phi }_2\left({\mathbf{x}}_I\right),\cdots ,{\phi }_n\left({\mathbf{x}}_I\right)\}}^T$ and ${\mathbf{u}}_s={\{u\left({\mathbf{x}}_1\right),u\left({\mathbf{x}}_2\right),\cdots ,u\left({\mathbf{x}}_n\right)\}}^T$.

After constructing the shape function of each integration point, and its spatial partial derivatives, it is possible to apply then to a strong or weak form governing the problem domain and establish the system of equations, the resolution of which provides the problem’s solution variable fields (displacements, velocities, pressures, etc.)

2.1. The Radial Point Interpolation Method

In order to obtain the background integration points and to allow the numerical integration of the integro-differential equations governing the studied phenomenon, the RPIM classic formulation uses the Gauss–Legendre quadrature integration scheme. Therefore, the problem domain can be divided with regular integration cells, Figure 1, and then each integration cell is transformed into a isoparametric cube, in which the integration points are included following the Gauss–Legendre integration scheme (Figure 1). In the end, by using linear polynomial interpolating functions $\mathbf{N}(\xi ,\eta ,\zeta )={\{N_1(\xi ,\eta ,\zeta ),N_2(\xi ,\eta ,\zeta ),\cdots ,N_8(\xi ,\eta ,\zeta )\}}^T$, whose generic equation is given by

$$N_i(\xi ,\eta ,\zeta )=\frac{1}{8}·(1+{\widehat{\xi }}_i·\xi )·(1+{\widehat{\eta }}_i·\eta )·(1+{\widehat{\zeta }}_i·\zeta ),$$

(2)

considering

$$\left\{\begin{array}{c}\mathrm{\Xi }=\{{\widehat{\xi }}_1,{\widehat{\xi }}_2,\cdots ,{\widehat{\xi }}_8\}=\{-1,1,1,-1,-1,1,1,-1\}\hfill \\ \mathbf{H}=\{{\widehat{\eta }}_1,{\widehat{\eta }}_2,\cdots ,{\widehat{\eta }}_8\}=\{-1,-1,1,1,-1,-1,1,1\}\hfill \\ \mathbf{Z}=\{{\widehat{\zeta }}_1,{\widehat{\zeta }}_2,\cdots ,{\widehat{\zeta }}_8\}=\{-1,-1,-1,-1,1,1,1,1\},\hfill \end{array}\right.$$

(3)

the isoparametric coordinates can be transformed back into Cartesian coordinates (Figure 1), with

$${\mathbf{x}}_I=\left\{\begin{array}{c}{x}_I\\ y_I\\ z_I\end{array}\right\}=\left[\begin{array}{c}\mathbf{N}{({\xi }_I,{\eta }_I,{\zeta }_I)}^T\\ \mathbf{N}{({\xi }_I,{\eta }_I,{\zeta }_I)}^T\\ \mathbf{N}{({\xi }_I,{\eta }_I,{\zeta }_I)}^T\end{array}\right]·\left[\begin{array}{c}{\mathbf{x}}_1^T\\ {\mathbf{x}}_2^T\\ ⋮\\ {\mathbf{x}}_8^T\end{array}\right]=\left[\begin{array}{cccc}{N}_1& N_2& \cdots & N_8\\ N_1& N_2& \cdots & N_8\\ N_1& N_2& \cdots & N_8\end{array}\right]·\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ x_2& y_2& z_2\\ ⋮& ⋮& ⋮\\ x_8& y_8& z_8\end{array}\right].$$

(4)

For the RPIM formulation, inside each integration cell, it can be inserted in $2\times 2\times 2$ integration points, whose isoparametric coordinates and weights are, correspondingly,

$$\left\{\begin{array}{c}{\xi }_I\\ {\eta }_I\\ {\zeta }_I\\ {\omega }_I\end{array}\right\}=\left\{\begin{array}{cccccccc}\frac{-1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}\\ \frac{-1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\\ \frac{-1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{-1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\\ 1& 1& 1& 1& 1& 1& 1& 1\end{array}\right\}.$$

(5)

The integration weight of each integration point can be calculated with

$${\widehat{\omega }}_I=\frac{V_{cel}}{V_{iso}}·{\omega }_I$$

(6)

with $V_{cel}$ being the total volume of the integration cell and $V_{iso}$ the total volume of the isoparametric cell, which is $V_{iso}=2\times 2\times 2=8$.

Thus, following the Gauss–Legendre integration scheme, a given function $f(x,y,z)$ defined inside a hexagonal domain ${\mathrm{\Omega }}_h$ discretized with $2\times 2\times 2$ integration points, can be integrated with

$$F={\int }_{{\mathrm{\Omega }}_h}f(x,y,z)d{\mathrm{\Omega }}_h=\sum _{I=1}^{8}f(x_I,y_I,z_I)·{\widehat{\omega }}_I,$$

(7)

for which each ${\mathbf{x}}_I=\{x_I,y_I,z_I\}$ is obtained with Equation (4) and the corresponding ${\widehat{\omega }}_I$ is indicated in Equation (5).

This procedure is repeated for each integration cell, allowing us to discretize the complete solid domain with a set of $n_Q$ integration points (known as background integration mesh). Summing all integration points’ weight gives the solid domain total volume, $V_{\mathrm{\Omega }}$:

$$\sum _{I=1}^{n_Q}{\widehat{\omega }}_I=V_{\mathrm{\Omega }}.$$

(8)

As already mentioned, the RPIM enforces the nodal connectivity by establishing overlapped influence domains. In this work, a spherical influence domain will be used, imposing a fixed number of nodes inside each influence domain. Thus, each integration point ${\mathbf{x}}_I$ will search for the n closest nodes and those nodes will set the influence domain of ${\mathbf{x}}_I$. According to Wang et al. [20], for a 2D problem, each influence domain should possess between $n=9$ and $n=16$ nodes. Thus, respectively, and in accordance with [32], for a 3D problem, each influence domain should possess between $n=27$ and $n=64$ nodes. Respecting the same number of nodes inside each influence domain allows the construction of shape functions with the same level of complexity for each integration point ${\mathbf{x}}_I$.

Regarding the construction of the shape functions, RPIM uses the radial point interpolating (RPI) technique. Thus, consider a function $\mathbf{u}\left(\mathbf{x}\right)$ defined in the domain $\mathrm{\Omega }\in {\mathbb{R}}^3$ discretized by a set of N nodes: $\mathit{X}=\{{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N\}\in {\mathbb{R}}^3$. Assuming an integration point, ${\mathbf{x}}_I$, possessing an influence domain with n nodes inside, it is possible to obtain the function’s interpolated value for ${\mathbf{x}}_I$ by using the following construction:

$$\begin{array}{c}\hfill u\left({\mathbf{x}}_I\right)=\sum _{i=1}^n{r}_i\left({\mathbf{x}}_I\right)·a_i\left({\mathbf{x}}_I\right)+\sum _{j=1}^m{p}_j\left({\mathbf{x}}_I\right)·b_j\left({\mathbf{x}}_I\right)=\\ \hfill =\mathbf{r}{\left({\mathbf{x}}_I\right)}^T·\mathbf{a}\left({\mathbf{x}}_I\right)+\mathbf{p}{\left({\mathbf{x}}_I\right)}^T·{\mathbf{b}}^T\left({\mathbf{x}}_I\right)=\{\mathbf{r}{\left({\mathbf{x}}_I\right)}^T,\mathbf{p}{\left({\mathbf{x}}_I\right)}^T\}\left\{\begin{array}{c}\mathbf{a}\\ \mathbf{b}\end{array}\right\},\end{array}$$

(9)

where $\mathbf{r}\left({\mathbf{x}}_I\right)$ is a radial basis function, $\mathbf{p}\left({\mathbf{x}}_I\right)$ is a polynomial basis function, with m monomials, and $a_i\left({\mathbf{x}}_I\right)$ and $b_j\left({\mathbf{x}}_I\right)$ are nonconstant coefficients of $\mathbf{r}\left({\mathbf{x}}_I\right)$ and $\mathbf{p}\left({\mathbf{x}}_I\right)$, respectively.

The vectors present in Equation (9) are defined by

$$\begin{array}{cc}\hfill \mathbf{r}\left({\mathbf{x}}_I\right)& ={\{r_1\left({\mathbf{x}}_I\right),r_2\left({\mathbf{x}}_I\right),\cdots ,r_n\left({\mathbf{x}}_I\right)\}}^T\hfill \\ \hfill \mathbf{p}\left({\mathbf{x}}_I\right)& ={\{p_1\left({\mathbf{x}}_I\right),p_2\left({\mathbf{x}}_I\right),\cdots ,p_m\left({\mathbf{x}}_I\right)\}}^T\hfill \\ \hfill \mathbf{a}\left({\mathbf{x}}_I\right)& ={\{a_1\left({\mathbf{x}}_I\right),a_2\left({\mathbf{x}}_I\right),\cdots ,a_n\left({\mathbf{x}}_I\right)\}}^T\hfill \\ \hfill \mathbf{b}\left({\mathbf{x}}_I\right)& ={\{b_1\left({\mathbf{x}}_I\right),b_2\left({\mathbf{x}}_I\right),\cdots ,b_m\left({\mathbf{x}}_I\right)\}}^T.\hfill \end{array}.$$

(10)

Complete polynomial basis functions are added to the basis function to ensure the consistency of the RPI functions, e.g., adding to the RBF a linear polynomial ensures $C^1$ consistency and helps the RPI pass the standard patch test. Generally, low-order polynomial basis are considered, such as

$$\begin{array}{c}\hfill \mathbf{p}\left({\mathbf{x}}_I\right)=\left\{1\right\},m=1\\ \hfill \mathbf{p}\left({\mathbf{x}}_I\right)={\{1,x_I,y_I,z_I\}}^T,m=4\\ \hfill \mathbf{p}\left({\mathbf{x}}_I\right)={\{1,x_I,y_I,z_I,x_I^2,x_I·y_I,y_I^2,y_I·z_I,z_I^2,z_I·x_I\}}^T,m=10.\end{array}.$$

(11)

As recommended in the literature [1], for efficiency purposes, in this work, a constant polynomial basis function is applied, $m=1$.

Regarding the radial basis function, this work assumes the multiquadratic radial basis function (MQ-RBF), proposed initially by Hardy [33],

$$r_j\left({\mathbf{x}}_i\right)=r_{ij}={(d_{ij}^2+c^2)}^p={\left({\left(\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+{(z_j-z_i)}^2}\right)}^2+c^2\right)}^p,$$

(12)

in which c and p are the MQ-RBF shape parameters. For 3D analyses, the literature [1] shows that c should be close to zero, but not be zero, and p should be close to one, but not be one. Therefore, following the literature recommendation, in this work the following MQ-RBF shape parameters are considered: $c={10}^{-4}$ and $p=1-{10}^{-4}$. The MQ-RBF shape parameters were obtained by optimizing the Kronecker delta function of the 3D RPI shape functions. The results shown in [1] proved to be accurate and stable, allowing us to reproduce highly complex variable fields.

The polynomial basis has to satisfy an extra requirement in order to obtain a unique solution,

$$\sum _{i=1}^n{p}_j\left({\mathbf{x}}_i\right)a_i\left({\mathbf{x}}_i\right)=0$$

(13)

with $j=\{1,2,\cdots ,m\}$. Therefore, combining Equation (9) with Equation (13), a new equation matrix can be written as

$$\left[\begin{array}{cc}\mathbf{R}& \mathbf{P}\\ {\mathbf{P}}^T& 0\end{array}\right]\left\{\begin{array}{c}\mathbf{a}\\ \mathbf{b}\end{array}\right\}=\mathbf{G}\left\{\begin{array}{c}\mathbf{a}\\ \mathbf{b}\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{u}}_s\\ 0\end{array}\right\}$$

(14)

where ${\mathbf{u}}_s$ is given by

$${\mathbf{u}}_s={\{u_1,u_2,\cdots ,u_n\}}^T$$

(15)

with the symmetric matrix $\mathbf{R}$ obtained with

$$\mathbf{R}=\left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1n}\\ r_{21}& r_{22}& \cdots & r_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ r_{n1}& r_{n2}& \cdots & r_{nn}\end{array}\right]$$

(16)

and matrix $\mathbf{P}$, of size $[n\times m]$, of the polynomial basis being given by

$$\mathbf{p}=\left[\begin{array}{cccc}{p}_1\left({\mathbf{x}}_1\right)& p_2\left({\mathbf{x}}_1\right)& \cdots & p_m\left({\mathbf{x}}_1\right)\\ p_1\left({\mathbf{x}}_2\right)& p_2\left({\mathbf{x}}_2\right)& \cdots & p_m\left({\mathbf{x}}_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ p_1\left({\mathbf{x}}_n\right)& p_2\left({\mathbf{x}}_n\right)& \cdots & p_m\left({\mathbf{x}}_n\right)\end{array}\right].$$

(17)

The moment matrix $\mathbf{G}$ is symmetric because the distance is directionally independent. Solving Equation (14),

$$\left\{\begin{array}{c}\mathbf{a}\\ \mathbf{b}\end{array}\right\}={\mathbf{G}}^{-1}·\left\{\begin{array}{c}{\mathbf{u}}_s\\ 0\end{array}\right\}$$

(18)

and substituting Equation (18) into Equation (9), the interpolation is obtained:

$$u\left({\mathbf{x}}_I\right)=\{\mathbf{r}{\left({\mathbf{x}}_I\right)}^T,\mathbf{p}{\left({\mathbf{x}}_I\right)}^T\}·{\mathbf{G}}^{-1}·\left\{\begin{array}{c}{\mathbf{u}}_s\\ 0\end{array}\right\}=\{\mathrm{\Phi }\left({\mathbf{x}}_I\right),\mathrm{\Psi }\left({\mathbf{x}}_I\right)\}·\left\{\begin{array}{c}{\mathbf{u}}_s\\ 0\end{array}\right\},$$

(19)

where $\mathrm{\Phi }\left({\mathbf{x}}_I\right)$ is the interpolation of the integration point ${\mathbf{x}}_I$,

$$\{\mathrm{\Phi }\left({\mathbf{x}}_I\right),\mathrm{\Psi }\left({\mathbf{x}}_I\right)\}=\{\mathbf{r}{\left({\mathbf{x}}_I\right)}^T,\mathbf{p}{\left({\mathbf{x}}_I\right)}^T\}·{\mathbf{G}}^{-1}=\{{\varphi }_1\left({\mathbf{x}}_I\right),\cdots ,{\varphi }_n\left({\mathbf{x}}_I\right),{\psi }_1\left({\mathbf{x}}_I\right),\cdots ,{\psi }_m\left({\mathbf{x}}_I\right)\}$$

(20)

In order to achieve a well-conditioned moment matrix, $\mathbf{G}$, the number of monomials, m, of the polynomial basis function should be much lower than the number of nodes, n, inside the influence domain of ${\mathbf{x}}_I$: $m\ll n$. Such a relation is assured because in 3D analyses, the polynomials have $m=1$, $m=4$ and $m=10$ for constant, linear, and quadratic polynomial basis functions, respectively, and the number of nodes inside the influence domain is typically $n\ge 27$.

The Galerkin weak form depends on the partial derivative of $\mathrm{\Phi }\left({\mathbf{x}}_I\right)$. Thus, in regard to a general direction $\xi$, the partial derivatives can be written as

$${\mathrm{\Phi }}_{\xi }\left({\mathbf{x}}_I\right)=\{\mathbf{r}{\left({\mathbf{x}}_I\right)}_{\xi }^T,\mathbf{p}{\left({\mathbf{x}}_I\right)}_{\xi }^T\}·{\mathbf{G}}^{-1},$$

(21)

with the partial derivatives of the MQ-RBF in order to $\xi$ being defined as

$${\left(r_{ij}\right)}_{\xi }=2·p·{(d_{ij}^2+c^2)}^{p-1}·({\xi }_j-{\xi }_i).$$

(22)

As the literature shows, RPI functions possess several useful properties, such as the Kronecker delta property and satisfy the partition of unity. Being interpolating functions, they allow one to impose directly the essential and natural boundary conditions, reducing the overall computational cost of the analysis when compared with approximating meshless methods [1].

2.2. Discrete System of Equations

Being a discrete numerical technique, RPIM uses the Galerkin weak formulation to construct the global system of equations. Thus, considering the solid domain $\mathrm{\Omega }$, with a body force $\mathbf{b}$, bounded by a surface boundary $\mathrm{\Gamma }$, submitted to external forces $\overline{t}$, on the natural boundary ${\mathrm{\Gamma }}_t$, and displacement constrains at the essential boundary ${\mathrm{\Gamma }}_u$, the Galerkin weak form can be expressed as

$${\int }_{\mathrm{\Omega }}\delta {\mathsf{\epsilon }}^T·\mathsf{\sigma }·d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}\delta \mathbf{u}{\left({\mathbf{x}}_I\right)}^T·\mathbf{b}·d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_t}\delta \mathbf{u}{\left({\mathbf{x}}_I\right)}^T·\overline{\mathbf{t}}·d\mathrm{\Gamma }=0,$$

(23)

recalling that because $u\left({\mathbf{x}}_I\right)=\mathrm{\Phi }{\left({\mathbf{x}}_I\right)}^T·\mathbf{u}$, it is possible to interpolate all the components of the field variable with one operation:

$$\begin{array}{c}\hfill \mathbf{u}\left({\mathbf{x}}_I\right)=\left[\begin{array}{c}u\left({\mathbf{x}}_I\right)\\ v\left({\mathbf{x}}_I\right)\\ w\left({\mathbf{x}}_I\right)\end{array}\right]=\mathbf{H}\left({\mathbf{x}}_I\right)·\mathbf{u}=\\ \hfill =\left[\begin{array}{ccccccc}{\varphi }_1\left({\mathbf{x}}_I\right)& 0& 0& \cdots & {\varphi }_n\left({\mathbf{x}}_I\right)& 0& 0\\ 0& {\varphi }_1\left({\mathbf{x}}_I\right)& 0& \cdots & 0& {\varphi }_n\left({\mathbf{x}}_I\right)& 0\\ 0& 0& {\varphi }_1\left({\mathbf{x}}_I\right)& \cdots & 0& 0& {\varphi }_n\left({\mathbf{x}}_I\right)\end{array}\right]·\left[\begin{array}{c}{u}_1\\ v_1\\ w_1\\ ⋮\\ u_n\\ v_n\\ w_n\end{array}\right]\end{array}$$

(24)

and express the deformation vector as

$$\mathsf{\epsilon }\left({\mathbf{x}}_I\right)=\mathbf{L}·\mathbf{u}\left({\mathbf{x}}_I\right)=\left[\begin{array}{ccc}\frac{\partial }{\partial x}& 0& 0\\ 0& \frac{\partial }{\partial y}& 0\\ 0& 0& \frac{\partial }{\partial z}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial z}& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial x}\end{array}\right]·\mathrm{\Phi }\left({\mathbf{x}}_I\right)·\mathbf{u}=\mathbf{B}\left({\mathbf{x}}_I\right)·\mathbf{u}=\mathbf{B}\left({\mathbf{x}}_I\right)·\left[\begin{array}{c}{u}_1\\ v_1\\ w_1\\ ⋮\\ u_n\\ v_n\\ w_n\end{array}\right]$$

(25)

being $\mathbf{B}\left({\mathbf{x}}_I\right)$ the deformation matrix of integration point ${\mathbf{x}}_I$, which can be defined as

$$\mathbf{B}\left({\mathbf{x}}_I\right)=\left[\begin{array}{ccccccc}\frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial x}& 0& 0& \cdots & \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial x}& 0& 0\\ 0& \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial y}& 0& \cdots & 0& \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial y}& 0\\ 0& 0& \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial z}& \cdots & 0& 0& \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial z}\\ \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial y}& \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial x}& 0& \cdots & \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial y}& \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial x}& 0\\ 0& \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial z}& \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial y}& \cdots & 0& \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial z}& \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial y}\\ \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial z}& 0& \frac{\partial {\varphi }_1\left({\mathbf{x}}_I\right)}{\partial x}& \cdots & \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial z}& 0& \frac{\partial {\varphi }_n\left({\mathbf{x}}_I\right)}{\partial x},\end{array}\right]$$

(26)

where n is the number of nodes inside the influence domain of integration point ${\mathbf{x}}_I$.

Then, with the linear strain–stress relation $\mathsf{\sigma }\left({\mathbf{x}}_I\right)=\mathbf{c}\left({\mathbf{x}}_I\right)·\mathsf{\epsilon }\left({\mathbf{x}}_I\right)$, stress at an integration point ${\mathbf{x}}_I$ can be defined by

$$\mathsf{\sigma }\left({\mathbf{x}}_I\right)=\mathbf{c}\left({\mathbf{x}}_I\right)·\mathsf{\epsilon }\left({\mathbf{x}}_I\right)=\mathbf{c}\left({\mathbf{x}}_I\right)·\mathbf{B}\left({\mathbf{x}}_I\right)·\mathbf{u}$$

(27)

where $\mathbf{c}\left({\mathbf{x}}_I\right)$ is the material constitutive matrix of integration point ${\mathbf{x}}_I$. Substituting the stress and the strain in Equation (23) allows us to obtain

$${\int }_{\mathrm{\Omega }}\delta {\left(\mathbf{B}\left({\mathbf{x}}_I\right)·\mathbf{u}\right)}^T·\left(\mathbf{c}\left({\mathbf{x}}_I\right)·\mathbf{B}\left({\mathbf{x}}_I\right)·\mathbf{u}\right)·d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}\delta {\left(\mathbf{H}\left({\mathbf{x}}_I\right)·\mathbf{u}\right)}^T·\mathbf{b}·d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_t}\delta {\left(\mathbf{H}\left({\mathbf{x}}_I\right)·\mathbf{u}\right)}^T·\overline{\mathbf{t}}·d\mathrm{\Gamma }=0$$

(28)

because only the small strain will be considered, $\delta \mathbf{B}\left({\mathbf{x}}_I\right)=0$ and $\delta \mathbf{H}\left({\mathbf{x}}_I\right)=0$, the expression can be presented as

$$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}\delta {\mathbf{u}}^T·\mathbf{B}{\left({\mathbf{x}}_I\right)}^T·\mathbf{c}\left({\mathbf{x}}_I\right)·\mathbf{B}\left({\mathbf{x}}_I\right)·\mathbf{u}·d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}\delta {\mathbf{u}}^T·\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\mathbf{b}·d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_t}\delta {\mathbf{u}}^T·\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\overline{\mathbf{t}}·d\mathrm{\Gamma }=\\ \hfill =\delta {\mathbf{u}}^T{\int }_{\mathrm{\Omega }}\mathbf{B}{\left({\mathbf{x}}_I\right)}^T·\mathbf{c}\left({\mathbf{x}}_I\right)·\mathbf{B}\left({\mathbf{x}}_I\right)·d\mathrm{\Omega }·\mathbf{u}-\delta {\mathbf{u}}^T{\int }_{\mathrm{\Omega }}\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\mathbf{b}·d\mathrm{\Omega }-\delta {\mathbf{u}}^T{\int }_{{\mathrm{\Gamma }}_t}\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\overline{\mathbf{t}}·d\mathrm{\Gamma }=0\end{array}$$

(29)

leading to its final form, representing the system of equations

$$\mathbf{K}·\mathbf{u}-{\mathbf{f}}_b-{\mathbf{f}}_t=0\Rightarrow \mathbf{u}={\mathbf{K}}^{-1}·\left({\mathbf{f}}_b+{\mathbf{f}}_t\right),$$

(30)

in which $\mathbf{K}$, $\mathbf{u}$, ${\mathbf{f}}_b$ and ${\mathbf{f}}_t$ are the global stiffness matrix, the global displacement field, and the global body and external force vectors, respectively. The discretization procedure allows us to obtain these entities by using a numerical integration:

$$\mathbf{K}={\int }_{\mathrm{\Omega }}\mathbf{B}{\left({\mathbf{x}}_I\right)}^T·\mathbf{c}\left({\mathbf{x}}_I\right)·\mathbf{B}\left({\mathbf{x}}_I\right)·d\mathrm{\Omega }=\sum _{I=1}^{N_Q}\mathbf{B}{\left({\mathbf{x}}_I\right)}^T·\mathbf{c}\left({\mathbf{x}}_I\right)·\mathbf{B}\left({\mathbf{x}}_I\right)·{\widehat{\omega }}_I$$

(31)

$${\mathbf{f}}_b={\int }_{\mathrm{\Omega }}\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\mathbf{b}·d\mathrm{\Omega }=\sum _{I=1}^{N_Q}\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\mathbf{b}·{\widehat{\omega }}_I$$

(32)

$${\mathbf{f}}_t={\int }_{{\mathrm{\Gamma }}_t}\mathbf{H}{\left({\mathbf{x}}_I\right)}^T·\overline{\mathbf{t}}·d\mathrm{\Gamma }=0=\sum _{J=1}^{n_Q}\mathbf{H}{\left({\mathbf{x}}_J\right)}^T·\overline{\mathbf{t}}·{\widehat{\omega }}_J.$$

(33)

Because the RPI function possesses the Kronecker delta property, the essential and natural boundary conditions can be directly imposed by using the direct imposition method or penalty methods, such as in other discrete numerical methods using interpolation functions.

3. Elastoplastic Analysis

The mathematical theory of plasticity provides a theoretical description for the relationship between stress and strain for a material exhibiting an elastoplastic behaviour [23]. Material plasticity is characterized by an irreversible deformation. Before occuring, plastic flow, the material shows a linear elastic behaviour, for which a linear relation following the generalized Hooke’s law is valid. After the yield stress ${\sigma }_Y$ is reached (established by the yield criterion), plastic behaviour occurs and the stress field starts following a flow rule. In order to implement material plasticity, it is necessary to establish three fundamental concepts:

• a yield criterion, indicating the stress level at which plastic flow initiates, and an initial yield surface, defining the elastic limit in a multiaxial stress state;
• a hardening law, describing the evolution of the subsequent yield surfaces; and
• a flow rule, relating the plastic deformation to the current stresses in the material and the stress increments outside the yield surface.

Concerning the yield criterion, in this work, the von Mises yield criterion is considered:

$$F\left(\mathsf{\sigma }\right)=\sqrt{\frac{1}{2}\left[{({\sigma }_{xx}-{\sigma }_{yy})}^2+{({\sigma }_{yy}-{\sigma }_{zz})}^2+{({\sigma }_{zz}-{\sigma }_{xx})}^2+6({\sigma }_{yz}^2+{\sigma }_{zx}^2+{\sigma }_{xy}^2)\right]}-{\sigma }_Y=0.$$

(34)

Assuming small strains, the total strain $\mathsf{\epsilon }$ can be decomposed into two components [34,35],

$$\mathsf{\epsilon }={\mathsf{\epsilon }}^p+{\mathsf{\epsilon }}^e,$$

(35)

where ${\mathsf{\epsilon }}^p$ is the plastic component of the total strain and ${\mathsf{\epsilon }}^e$ the elastic component. Concerning the hardening law, this work assumes that the yield surface grows linearly from the initial yield surface, only changing its size and expanding uniformly with the increase of the effective plastic strain, ${\overline{\mathsf{\epsilon }}}^{\mathit{p}}$:

$${\overline{\mathsf{\epsilon }}}^p=\int d{\overline{\mathsf{\epsilon }}}^p=\int {\left[\frac{2}{3}d{\overline{\mathsf{\epsilon }}}_{ij}^p·d{\overline{\mathsf{\epsilon }}}_{ij}^p\right]}^{\frac{1}{2}}$$

(36)

where $d{\mathsf{\epsilon }}^p$ is the plastic component of strain occurring during a strain increment.

The flow rule serves as a mathematical description of the evolution of the infinitesimal increments of the plastic strain $d{\mathsf{\epsilon }}^p$ in the course of the load history of the body [36]. As the incremental load increases, and the yield function f reaches a value equal to the yield stress of the material, ${\sigma }_Y$, the frontier of elasticity starts to end and plasticity begins. The variation of the yield function can be written as a function of the variation of the stress field

$$df={\left(\frac{\partial f}{\partial \mathsf{\sigma }}\right)}^T·d\mathsf{\sigma },$$

(37)

where $\partial f/\partial \mathsf{\sigma }$ is the gradient of f, and therefore an orthogonal vector to the yield surface for a certain stress increment $d\mathsf{\sigma }$. Applying the Prandtl–Reuss flow rule, it is possible to write

$$d{\mathsf{\epsilon }}^p=d\lambda ·\frac{\partial f}{\partial \mathsf{\sigma }},$$

(38)

where $d\lambda$ is the plastic multiplier (a scalar). The vector $\partial f/\partial \mathsf{\sigma }$ is a vector normal to the yield surface at the current stress state, as shown in Figure 2a [23,36], also known as a flow vector: $\mathbf{a}=\partial f/\partial \mathsf{\sigma }$.

3.1. Elastoplastic Constitutive Model

For its implementation in computational mechanics, it is convenient to present the previous elastoplastic theory in matrix formulation [23]. Thus, generically, the yield criterion of Equation (34) can be defined as

$$F(\mathsf{\sigma },\kappa )=f(\mathsf{\sigma },\kappa )-{\sigma }_Y\left(\kappa \right)=0,$$

(39)

where $\kappa$ is a generic hardening parameter (which in this work will be assumed as the effective plastic deformation, Equation (36)). By differentiating Equation (39), it is possible to write

$$dF=\left(\frac{\partial F}{\partial \mathsf{\sigma }}\right)·d\mathsf{\sigma }+\frac{\partial F}{\partial \kappa }·d\kappa =0,$$

(40)

which is equivalent to Equation (37) but for a material with hardening. Thus, Equation (40) can be written as function of the plastic multiplier $d\lambda$,

$${\mathbf{a}}^T·d\mathsf{\sigma }-A·d\lambda =0,$$

(41)

where the vector $\mathbf{a}$ is termed the flow vector, $\mathbf{a}=\partial f/\partial \mathsf{\sigma }$. The hardening parameter A depends on the hardening rule [23], and it can be written as

$$A=-\frac{1}{d\lambda }·\frac{\partial F}{\partial \kappa }·d\kappa .$$

(42)

From decomposition equation, Equation (35), it possible to write ${\mathsf{\epsilon }}^e=\mathsf{\epsilon }-{\mathsf{\epsilon }}^p$. Thus, with Equation (38), $d{\mathsf{\epsilon }}^p=d\lambda ·\mathbf{a}$, and using the Hooke law,

$$d\mathsf{\sigma }=\mathbf{c}·d{\mathsf{\epsilon }}^e=\mathbf{c}·\left(d\mathsf{\epsilon }-d{\mathsf{\epsilon }}^p\right)=\mathbf{c}·\left(d\mathsf{\epsilon }-d\lambda ·\mathbf{a}\right),$$

(43)

it is possible to develop Equation (41) as

$${\mathbf{a}}^T·\left(\mathbf{c}·\left(d\mathsf{\epsilon }-d\lambda ·\mathbf{a}\right)\right)-A·d\lambda =0⟹d\lambda =\frac{{\mathbf{a}}^T·\mathbf{c}·d\mathsf{\epsilon }}{{\mathbf{a}}^T·\mathbf{c}·\mathbf{a}+A}.$$

(44)

Substituting Equation (44) in Equation (43) permits us to obtain

$$d\mathsf{\sigma }=\mathbf{c}·d\mathsf{\epsilon }-\frac{{\mathbf{a}}^T·\mathbf{c}·d\mathsf{\epsilon }}{{\mathbf{a}}^T·\mathbf{c}·\mathbf{a}+A}·\mathbf{c}·\mathbf{a}=\left(\mathbf{c}-\frac{\mathbf{c}·\mathbf{a}·{\mathbf{a}}^T\mathbf{c}}{{\mathbf{a}}^T·\mathbf{c}·\mathbf{a}+A}\right)·d\mathsf{\epsilon }={\mathbf{c}}^p·d\mathsf{\epsilon },$$

(45)

with ${\mathbf{c}}^p$ being the elastoplastic tangential constitutive matrix. The elastic constitutive matrix can be defined as

$$\mathbf{c}={\left[\begin{array}{cccccc}\frac{1}{E}& \frac{-\nu }{E}& \frac{-\nu }{E}& 0& 0& 0\\ \frac{-\nu }{E}& \frac{1}{E}& \frac{-\nu }{E}& 0& 0& 0\\ \frac{-\nu }{E}& \frac{-\nu }{E}& \frac{1}{E}& 0& 0& 0\\ 0& 0& 0& \frac{2·(\nu +1)}{E}& 0& 0\\ 0& 0& 0& 0& \frac{2·(\nu +1)}{E}& 0\\ 0& 0& 0& 0& 0& \frac{2·(\nu +1)}{E}\end{array}\right]}^{-1},$$

(46)

where E is the elastic modulus of the material (observed in the linear elastic phase) and $\nu$ the Poisson’s coefficient.

The hardening parameter A can be obtained by means of the hardening hypothesis [23] considering the associated flow rule. In this work, the analysed nonlinear materials exhibit a “linear-elastic/linear plastic” hardening behaviour. Thus, as shown in the literature [23], the hardening parameter A can be written as

$$A=\frac{E_T}{1-\frac{E_T}{E}},$$

(47)

where $E_T$ is the tangent modulus (tangent to the stress–strain curve in the plastic regime).

3.2. Stress-Returning Algorithm

In this work, an incremental procedure is adopted, in which an incremental stress leads to an incremental strain. Thus, trial stresses are a sum of a previous stress state with an incremental stress state. Consequently, it is possible to obtain trial stress states beyond the yield surface. For such case, the trial stress state must be forced to return to the yield surface. In this work, a “backward-Euler” procedure is adopted [23], Figure 2b.

Considering a trial stress state ${\mathsf{\sigma }}^{trial}\left({\mathbf{x}}_I\right)$ at a given integration point ${\mathbf{x}}_I$, the trial stress state is constituted by the previous stress state on the integration point, $\mathsf{\sigma }\left({\mathbf{x}}_I\right)$, and the incremental stress state of the present load increment (or iteration), $d\mathsf{\sigma }\left({\mathbf{x}}_I\right)$,

$${\mathsf{\sigma }}^{trial}\left({\mathbf{x}}_I\right)=\mathsf{\sigma }\left({\mathbf{x}}_I\right)+d\mathsf{\sigma }\left({\mathbf{x}}_I\right),$$

(48)

if the trial stress state does not respect the yield criterion, Equation (39), being observed:

$$f({\mathsf{\sigma }}^{trial}\left({\mathbf{x}}_I\right),\kappa )>{\sigma }_Y\left(\kappa \right),$$

(49)

then, it is necessary to return the stress state to the yield surface. As Figure 2b shows, consider a two-dimensional Westgaard stress space in which a stress state represented by point B has passed beyond the yield surface. In order to avoid extra, unnecessary computations, instead calculating the intersection stress state represented with point A, it is assumed a “forward” yield function $f_B$, containing the stress state B,

$$f_B=\overline{{\sigma }_B}-{\sigma }_Y^{\ast },$$

(50)

in which ${\sigma }_Y^{\ast }$ is the updated yield stress, accounting for the hardening parameter

$${\sigma }_Y^{\ast }=\sigma +[E_T·\overline{{\epsilon }_{i-1}^p}],$$

(51)

where $\overline{{\epsilon }_{i-1}^p}$ is the accumulated effective plastic strain from the previous increment, $(i-1)$.

Thus, the flow vector is calculated from stress state B, ${\mathbf{a}}_B$, and the plastic multiplier $d\lambda$ is given by

$$d\lambda =\frac{f_B}{{\mathbf{a}}_B^T·\mathbf{c}·{\mathbf{a}}_B+H^{\prime }}.$$

(52)

Now it is possible to estimate the stress in point C:

$${\mathsf{\sigma }}_C={\mathsf{\sigma }}_B-d\lambda ·\mathbf{c}·{\mathbf{a}}_B.$$

(53)

As can be observed, stress state ${\mathsf{\sigma }}_C$ is an estimated stress state. In order to achieve accurately stress state ${\mathsf{\sigma }}_C$, it would be required to use the flow vector calculated directly in point C,

$${\mathsf{\sigma }}_C={\mathsf{\sigma }}_B-d\lambda ·\mathbf{c}·{\mathbf{a}}_C.$$

(54)

However, stress state ${\mathsf{\sigma }}_C$ is unknown, naturally hindering the calculation of ${\mathbf{a}}_C$. Therefore, an iterative procedure must be implemented to reach point C from point B. For each “stress return iteration”, the effective plastic strain is updated in order to apply Equation (51) and pass to the next “stress return iteration”, until point C is satisfactorily approximated, concluding the iterative process.

3.3. Nonlinear Solution Algorithm

The elastoplastic analysis of structures involves a set of nonlinear equations that need to be solved in order to obtain the problem’s corresponding solution. Generally, such nonlinear equations are solved by using iterative-incremental methods. In this work, the Newton–Raphson method is used as a nonlinear solution technique. Therefore, for a given increment i and iteration j, the the displacement, $\delta u_i^j$ for each load increment $\delta f_i^j$ is given by

$$\delta u_i^j={\left(\frac{df}{du}\right)}^{-1}·\delta f_i^j=K_T^{-1}·\delta f_i^j$$

(55)

and the total displacements are given by

$$u_i^j=u_i^{j-1}+\partial u_i^{j-1}.$$

(56)

During the incremental-iterative process, the increment of load acts as a predictor. It provides a starting solution, $u_i^0$, in order for the iterative process to have a starting point to begin the convergence process. In addition, the complete version of the Newton–Raphson method, in which the stiffness matrix is calculated in each iteration (of each increment), there are several other versions documented in the literature [23]. For instance, the KT0 approach, in which only the initial stiffness matrix, calculated in the first iteration of the first increment (the elastic stiffness matrix), is considered for the complete analysis. Another version is the KT1, in which the stiffness matrix is calculated only at the first iteration of each increment, being then used in all the iterations of such increments. Both KT0 and KT1 versions allow us to reduce the number of computations in slightly nonlinear problems; however, for highly nonlinear problems the best approach is to use the complete version of the Newton–Raphson method, which allows us to significantly reduce the number of iterations in each increment.

The complete version of the Newton–Raphson used in this work can be summarised as follows. First, the problem domain is discretized with a nodal mesh, the background integration mesh is constructed and the nodal connectivity is imposed. Then, the RPI shape functions are constructed for each integration point. In this work, an incremental-iterative procedure was implemented for the elastoplastic deformation problem. Thus, the discretized system of equation represented in Equation (30) can be written as

$${\mathbf{K}}_T·\Delta \mathbf{u}-\Delta \mathbf{f}={\mathbf{f}}_{res},$$

(57)

where ${\mathbf{K}}_T$ is the tangent stiffness matrix, $\Delta \mathbf{u}$ the incremental displacement field, $\Delta \mathbf{f}$ is the incremental load vector, and ${\mathbf{f}}_{res}$ is the residual force vector. Thus, the process can start for the first increment ($i=1$).

1.

Set the incremental load, $\Delta \mathbf{f}={\mathbf{f}}^i$.

2.

Loop over all integration points, ${\mathbf{x}}_I$.

(a)

Construct $\mathbf{c}{\left({\mathbf{x}}_I\right)}_p$, Equation (45).

(b)

Construct $\mathbf{K}{\left({\mathbf{x}}_I\right)}_T$, Equation (31).

3.

Assemble the global tangent stiffness matrix: ${\mathbf{K}}_T=A_{I=1}^{n_Q}\mathbf{K}{\left({\mathbf{x}}_I\right)}_T$.

4.

Solve $\Delta {\mathbf{u}}^i={\mathbf{K}}_T^{-1}·\Delta \mathbf{f}$.

5.

Update the displacement field: ${\mathbf{u}}^i={\mathbf{u}}^{i-1}+\Delta {\mathbf{u}}^i$.

6.

Loop over all integration points, ${\mathbf{x}}_I$.

(a)

Evaluate the strain increment: $\Delta \mathsf{\epsilon }{\left({\mathbf{x}}_I\right)}^i=\mathbf{B}\left({\mathbf{x}}_I\right)·\Delta {\mathbf{u}}^i$.

(b)

Considering the elastic constitutive matrix, Equation (46), obtain the incremental stress state: $\Delta \sigma {\left({\mathbf{x}}_I\right)}^i=\mathbf{c}\left({\mathbf{x}}_I\right)·\Delta \mathsf{\epsilon }{\left({\mathbf{x}}_I\right)}^i$.

(c)

Calculate the trial stress state: $\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{trial}^i=\sigma {\left({\mathbf{x}}_I\right)}^{i-1}+\Delta \sigma {\left({\mathbf{x}}_I\right)}^i$.

(d)

Verify if: $f\left(\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{trial}^i\right)\le {\sigma }_Y\left(\kappa \right)$. If so, $\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{new}^i=\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{trial}^i$. Otherwise, the algorithm presented in Section 3.2 is applied, which will lead to $\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{new}^i=\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{returned}^i\ne \mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{trial}^i$.

(e)

Update the stress field: $\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}^i=\mathsf{\sigma }{\left({\mathbf{x}}_I\right)}_{new}^i$.

7.

Calculate the residual force vector: ${\mathbf{f}}_{res}^i=\Delta \mathbf{f}-{\int }_{\mathrm{\Omega }}{\mathbf{B}}^T\left({\mathsf{\sigma }}^i-{\mathsf{\sigma }}^{i-1}\right)d\mathrm{\Omega }$.

8.

Apply a residual force convergence criteria: $\mathrm{\Lambda }={\left({\left({\mathbf{f}}_{res}^i\right)}^T·{\mathbf{f}}_{res}^i\right)}^{1/2}·{\left({\left({\mathbf{f}}^i\right)}^T·{\mathbf{f}}^i\right)}^{-1/2}$.

(a)

if $\mathrm{\Lambda }\le \vartheta$, being $\vartheta$ a small threshold value (such as $\vartheta ={10}^{-3}$), then, the process moves to point 2, initialising a new increment: $i=i+1$, and assuming a new force increment: $\Delta \mathbf{f}={\mathbf{f}}^{i+1}$.

(b)

If $\mathrm{\Lambda }>\vartheta$, then the iterative process begins, attempting to achieve $\mathrm{\Lambda }\le \vartheta$. Thus, the process moves to point 2, initialising iteration j, assuming a new force increment: $\Delta \mathbf{f}={\mathbf{f}}_{res}^i$.

4. Numerical Results

In this section, simply supported and clamped 3D square plates are analysed. The three nonlinear solution algorithms (KT-0, KT-1, and KT-ALL) are used and compared by using the RPIM and the FEM for both plate configurations. The results obtained are compared with solutions documented in the literature [26], which were obtained by using the element free Galerkin method (EFGM) and an equivalent single-layer deformation theory, the first order shear deformation theory, and matched/validated with advanced FEM formulations, equally assuming an equivalent single layered deformation theory.

Regarding the geometric properties of the plates, following Figure 3 representation, the length L and thickness h of all analysed cases are: $L=6.0$ m and $h=0.2$ m. Concerning the plate’s material, the following mechanical properties were considered: elasticity (Young’s) modulus: $E=206.7$ GPa, Poisson’s coefficient: $\nu =0.3$, tangent elasticity modulus: $E_T=2.067$ GPa, and yield stress: ${\sigma }_Y=206.7$ MPa.

Regarding the loading conditions, all plates are submitted to a uniform distributed load covering the complete surface of the plate, as shown in Figure 3. In order to reduce the computational cost of the nonlinear analyses, instead computing the complete plate, only one quarter of the plate is analysed, assuming the corresponding boundary conditions. Thus, in surface $A{A}^{\prime }{B}^{\prime }B$ the displacement on $Ox$ direction is assumed constrained: $\overline{u}=0$, in surface $A{A}^{\prime }{D}^{\prime }D$ the displacement on $Oy$ direction is assumed constrained: $\overline{v}=0$, and in line $A{A}^{\prime }$ both $Ox$ and $Oy$ directions are constrained: $\overline{u}=0\wedge \overline{v}=0$. These symmetry boundary conditions are common to both simply supported and clamped plate examples. In addition, for the simply supported plate, the line of nodes on the middle layer of the plate (with $z=0$) have all degrees of freedom constrained: $\overline{u}=0\wedge \overline{v}=0\wedge \overline{w}=0$. For the clamped plate, in addition to the symmetry boundary conditions, all the points belonging to surfaces $B{B}^{\prime }{C}^{\prime }C$ and $D{D}^{\prime }{C}^{\prime }C$ have all degrees of freedom constrained: $\overline{u}=0\wedge \overline{v}=0\wedge \overline{w}=0$.

All examples are analysed by using the FEM and the RPIM. Concerning FEM, eight-node hexagonal elements are considered, assuming full integration: $2\times 2\times 2$ integration points following the Gauss–Legendre quadrature scheme. Regarding the RPIM, the background integration mesh coincides with the FEM element mesh, being formed by hexagonal integration cells filled with $2\times 2\times 2$ integration points respecting the Gauss–Legendre quadrature scheme, and each influence domain is composed of 27 nodes. The RPI shape functions are constructed with a constant polynomial basis, $\mathbf{p}\left({\mathbf{x}}_I\right)=1$, and the MQ-RPI shape parameters: $c={10}^{-4}$ and $p=1-{10}^{-4}$.

The problem domain is dicretized with a 3D regular nodal of $16\times 16\times 7=1792$ nodes, which corresponds to $15\times 15\times 6=1350$ hexagonal elements in FEM, or hexagonal integration cells in RPIM.

4.1. Normalized Central Displacement

In this section, the w displacements of point P (obtained along $Oz$) as a function of the applied load will be presented. Thus, in order to represent normalized load-displacement curves, the obtained w displacement are normalized with

$${\overline{w}}_P=\frac{w_P·D}{M_p·L^2},$$

(58)

and the force is normalized with

$$\overline{f}=\frac{f·L^2}{M_p},$$

(59)

being,

$$D=\frac{100·E·h^3}{12(1-{\nu }^2)},$$

(60)

and

$$M_p={\sigma }_Y\frac{h^2}{4},$$

(61)

where $w_P$ is the obtained transversal displacement (along $Oz$) on point P, h is the thickness of the plate, E the elasticity (Young’s) Modulus, $\nu$ the Poisson’s coefficient, $M_P$ the plastic moment, and ${\sigma }_Y$ the yield stress of the material being used.

The normalized load-transversal displacement curves are presented in Figure 4 for both simply supported and clamped plates, and for both numerical methods (FEM and RPIM).

Regarding the simply supported plate, it can be observed that by using the KT1 and KT-All nonlinear solution algorithms allows similar solutions for FEM and RPIM. Comparing the results with the literature, it is perceptible that both FEM and RPIM (using the KT-All nonlinear solution algorithm) are capable of achieving close results with the EFGM formulation by using an equivalent single layer theory [26].

For the clamped plate, Figure 4b shows that the FEM with KT1 and KT-All is capable of delivering a solution close to the reference solution. However, because the clamped plate is a more demanding benchmark, the solution is not as close as for the simply supported plate. FEM requires more nodes along the thickness to achieve a closer solution. For the RPIM (Figure 4d), both KT1 and KT-All algorithms produce solutions very close with the literature.

Notice that in the present work, a full 3D-deformation theory is considered, which is fundamentally different from the first-order shear deformation theory used by the reference solution. The full 3D-deformation theory strongly depends on the number of nodes along the thickness. Thus, the similarity between the obtained results and the solution documented in the literature indicates that both FEM and RPIM nonlinear solution algorithms were successfully implemented.

It would be expected that we would obtain similar results regardless of the nonlinear solution algorithm used (KT0, KT1 or KT-All). However, due to the total computational time of the analysis, per increment only a maximum of 100 iterations is allowed. The KT-All never requires more than roughly 10 iterations; however, when a large portion of integration point yield, KT1 starts to require more than 100 iterations per increment, and from third increment (very early), KT0 already requires more than 100 iterations per increment. This is because in each increment the solution is not fully corrected (the number of iterations is insufficient), the solutions start to accumulate convergence errors and, consequently, to diverge. This, observation highlights the advantage of using the full Newton–Raphson method (KT-ALL).

4.2. Stress Diagrams through the Plate’s Thickness

In this section, the stress distributions along with the plates’ thickness for distinct increment levels, and for the three numerical methods, will be presented. These results were obtained only with the KT-All nonsolution algorithm. As for units, the load increments are in [kN/m²], the stress value (horizontal axis) in MPa, and the thickness coordinate (vertical axis) is in meters.

First, the stress distribution along line A-${\mathrm{A}}^{\prime }$ is shown for the simply supported plate. Thus, in Figure 5, the FEM and RPIM results are shown side by side.

As a first observation, it is possible to visualise that the stress distribution along the thickness for the normal stresses ${\sigma }_{xx}$ and ${\sigma }_{yy}$ are virtually identical. Such an observation indicates that the symmetry was properly implemented. Additionally, the normal stresses obtained with RPIM are very similar to the ones obtained with the FEM. It is possible to visualise that, as the material starts to yield, the stress distribution along the thickness leaves the linear behaviour and starts showing a pronounced curve at the top and bottom layers, consistent with the maximum (yield) stress allowed.

The most significant differences can be found in the shear components of the stress tensor: ${\tau }_{xz}$ and ${\tau }_{yz}$. However, the differences become smaller as the load increment grows. At the centre of the plate, along the thickness, the distribution of shear stresses ${\tau }_{xz}$ and ${\tau }_{yz}$ should be zero. This is true for the first increment (around 807–815 kN/m²). Figure 5e–h does not show such a null line for the first increment because the stress on the nodes is extrapolated from the integration points in its vicinity (which do not belong exactly to the centre of the plate, and therefore possess non-null shear stresses). Afterward, yielding starts to occur at the centre of the plate and the stress field of the integration points along the thickness starts to increase due to the reconfiguration of the stress, a consequence of the stress-returning algorithm.

For the clamped plate, four interest lines along the plate thickness are analysed. Thus, Figure 6, Figure 7, Figure 8 and Figure 9 represent the stress distributions obtained along lines A-${\mathrm{A}}^{\prime }$, B-${\mathrm{B}}^{\prime }$, C-${\mathrm{C}}^{\prime }$ and D-${\mathrm{D}}^{\prime }$, as indicated in Figure 3, respectively.

Observing Figure 6, Figure 7, Figure 8 and Figure 9 it is noticeable that the load increments of FEM and RPIM are different. Such differences come from the elastoplastic algorithm implemented by the authors. The algorithm first identifies the load increment corresponding to the problem’s elastic limit, $f_0=f_{el}$; then such a load increment is multiplied by a factor, which in this work was 2.5, and an ultimate load is achieved: $f_u=2.5\times f_{el}+f_{el}$. Then, the number of increments is decided. In this work, 10 increments were used for each elastoplastic analysis. Thus, after the elastic limit, each incremental load will have the following magnitude: $f_i=2.5\times f_{el}/10$. Hence, if two distinct numerical methods predict two distinct elastic limits (because distinct numerical methods produce distinct numerical results), the ultimate load, and consequently the incremental load magnitude, of each numerical method will also be different.

Regarding line A-${\mathrm{A}}^{\prime }$, Figure 6, it is possible to visualise that FEM and RPIM produce similar results for similar load increment levels. Additionally, as already observed in the simply supported plate, the normal stresses ${\sigma }_{xx}$ and ${\sigma }_{yy}$ evolve from a linear elastic distribution along the thickness to a curved yield configuration at the top and bottom layers of the plate. The shear stresses ${\tau }_{xz}$ and ${\tau }_{yz}$ start with a value very close to zero (it is not zero due to the same reason already discussed for the simply supported plate) and both evolve to a bilinear configuration, with lower shear stress values in the top and bottom layers of the plate and maximum stress in the plate’s middle layer. Such stress evolution is due to the material yielding and the consequent required stress returning procedure.

Along line B-${\mathrm{B}}^{\prime }$, Figure 7 shows that the normal stresses ${\sigma }_{yy}$ obtained with FEM and RPIM are close for similar incremental load levels. However, the shear stress results are very different. The FEM obtains an expected shear stress ${\tau }_{yz}$ distribution, showing an unsymmetrical irregular path. In opposition, the RPIM produces the expected result, a parabolic distribution, symmetric with respect to the plate’s middle layer, that becomes gradually shaped (with a central peak at the plate’s middle layer) as the material yielding evolves.

Line C-${\mathrm{C}}^{\prime }$ was selected to study the variation of the shear stress ${\tau }_{xy}$. Being the plate’s corner, shear stresses ${\tau }_{xy}$ are at maximum at this location. Hence, Figure 8 presents, along the plate’s thickness, the evolution of shear stress ${\tau }_{xy}$ with the plate’s yielding. At first sight, the most evident observation is the difference between FEM and RPIM results. The shear stress ${\tau }_{xy}$ produced with FEM maintains its almost linear distribution along the thickness as the plate yields. On the other hand, with RPIM for the first load increments, ${\tau }_{xy}$ maintains a relatively low magnitude; however, as the material starts to yield at this location, the stress level rises rapidly.

Concerning the stress distribution along the line D-${\mathrm{D}}^{\prime }$, Figure 9, the FEM and RPIM results are consensual. The normal stress ${\sigma }_{xx}$ and the shear stress ${\tau }_{xz}$, obtain with FEM and RPIM, are very close. Furthermore, along this line D-${\mathrm{D}}^{\prime }$, the shear stress ${\tau }_{xz}$ obtained with FEM follows the same evolution as the one observed with RPIM as the plate yields.

4.3. Effective Stress Maps

In this section, the effective stress distribution maps for both plate configurations, simply supported and clamped plate, and both numerical method, FEM and RPIM, will be presented. The FEM and RPIM results were once again obtained assuming the full Newton–Raphson nonlinear solution algorithm (KT-All). After extrapolating the stress state from the integration points near the plate’s top layer to its closest nodes (belonging to the plate’s top surface), the (extrapolated) von Mises equivalent stress on each of those nodes was calculated:

$${\overline{\sigma }}_{eq}=\sqrt{\frac{1}{2}\left[{({\sigma }_{xx}-{\sigma }_{yy})}^2+{({\sigma }_{yy}-{\sigma }_{zz})}^2+{({\sigma }_{zz}-{\sigma }_{xx})}^2+6({\sigma }_{yz}^2+{\sigma }_{zx}^2+{\sigma }_{xy}^2)\right]}.$$

(62)

The obtained FEM von Mises equivalent stress maps for the simply supported plate are presented in Figure 10, and the RPIM results are shown in Figure 11. The load increment is in [kN/m²], and the maximum and minimum stress values of each incremental load (subfigure) are the ones indicated in the caption of the corresponding subfigure. It is possible to observe that both FEM and RPIM produce very similar results. As expected, the simply supported plate forms a rupture line along the plate diagonal, connecting the plate’s centre and the plate’s corner. Notice that the material enters the plastic regime along that diagonal line first, and then expands from that line to the remaining area as the load increment increases.

Regarding the clamped plate results, the corresponding von Mises equivalent stress maps are presented in Figure 12 and Figure 13 for FEM and RPIM, respectively.

Both FEM and RPIM solutions are again very similar. It is clear that the initiation of material yields at the boundaries of the plate, and then it starts to yield at the centre of the plate. Afterward, the yielding advances to a diagonal line merging the boundary yield material to the yield material of the plate’s centre.

Both the rupture paths observed in the simply supported and clamped plates are the ones described in the literature for these two benchmarks [26].

5. Conclusions

In this work, the radial point interpolation method (RPIM) was combined with an elastoplastic formulation to analyse 3D plates. Original nonlinear RPIM and FEM codes were written in order to have a more detailed comprehension of the differences and similarities between both numerical methods (at the formulation and results levels). This process allowed us to directly compare performances (global computational times and efficiency) and the results accuracy.

Concerning the nonlinear solutions algorithms tested, it was found that using the KT-All algorithm is more efficient. Although it recalculates the stiffness matrix in every iteration, it required many fewer iterations than KT0 or KT1 when a large part of the material starts to yield. Thus, the overall computation time is less than the one required by KT0 or KT1 algorithms to achieve the same solution.

Regarding the stress distribution along the thickness of the plate, the results show that generally FEM and RPIM produce similar normal stresses distributions. However, concerning shear stresses, the results are not as close. In opposition to FEM, RPIM is capable of producing the expected parabolic shear stress distribution along the thickness, demonstrating an higher accuracy. The stress distributions shown were acquired at the nodes. Therefore, such results are an extrapolation of the stress values on the integration points. This explains why shear stress obtained with FEM and RPIM is not null at the top and bottom layers of the plate (as it should be).

The presented stress distribution maps allowed us to visualise the expected rupture lines of both benchmarks. In the simply supported plate, the diagonal rupture line can identified, as well as in the clamped plate the boundary rupture lines at first, and then the central rupture and, in the end, the diagonal lines. Such results indicate that the nonlinear solution algorithm was successfully combined with the linear RPIM formulation.

The main advantage of RPIM is its potential to be directly integrated in existing FEM software. For instance, the RPIM can assume the finite elements as integration cells, distributing inside them the required integration points (following the most appropriated scheme). Moreover, the element nodes can be utilised to form the global nodal discretization. Then, by using the influence-domain concept, the nodal connectivity can be established. Afterward, the process moves on, constructing the RPI functions, establishing the system of equations, obtaining the displacement, strain and stress fields and following the elastoplastic analysis. The final results can then be depicted in existent FEM software. Moreover, because RPIM is an interpolating numerical method, it can be combined with FEM. Allowing to use meshless methods in parts of the domain in which higher accuracy is required (for instance, the vicinity of cracks or other locations in which stress concentration is expected) and the FEM can be applied in the remaining part of the solid (to speed-up the analysis).