# Numerical analysis of small-strain elasto-plastic deformation using local
Radial Basis Function approximation with Picard iteration

Filip Strniša<EMAIL_ADDRESS>Mitja Jančič<EMAIL_ADDRESS>Gregor
Kosec<EMAIL_ADDRESS>“Jožef Stefan” Institute, Parallel and Distributed
Systems Laboratory, Jamova cesta 39, 1000 Ljubljana, Slovenia

###### Abstract

This paper deals with a numerical analysis of plastic deformation under
various conditions, utilizing Radial Basis Function (RBF) approximation. The
focus is on the elasto-plastic von Mises problem under plane-strain
assumption. Elastic deformation is modelled using the Navier-Cauchy equation.
In regions where the von Mises stress surpasses the yield stress, corrections
are applied locally through a return mapping algorithm. The non-linear
deformation problem in the plastic domain is solved using the Picard
iteration. The solutions for the Navier-Cauchy equation are computed using the
Radial Basis Function-Generated Finite Differences (RBF-FD) meshless method
using only scattered nodes in a strong form. Verification of the method is
performed through the analysis of an internally pressurized thick-walled
cylinder subjected to varying loading conditions. These conditions induce
states of elastic expansion, perfectly-plastic yielding, and plastic yielding
with linear hardening. The results are benchmarked against analytical
solutions and traditional Finite Element Method (FEM) solutions. The paper
also showcases the robustness of this approach by solving case of thick-walled
cylinder with cut-outs. The results affirm that the RBF-FD method produces
results comparable to those obtained through FEM, while offering substantial
benefits in managing complex geometries without the necessity for conventional
meshing, along with other benefits of meshless methods.

###### keywords:

meshless , plasticity , non-linear , isotropic hardening , von Mises model ,
Picard iteration , RBF-FD , plane strain , thick-walled cylinder

††journal: Computer Methods in Applied Mechanics and Engineering

tcb@breakable

The elasto-plastic von Mises problem is solved in strong form on scattered
nodes using direct RBF-FD and Picard iteration.

The presented method provides results regardless the underlying geometry
without any special treatment of differential operators, boundary conditions
or spatial discretisation.

The presented method shows good agreement with analytical and Finite Element
Method solutions, as well as convergence behaviour.

The method effectively handles different stages of material behavior, from
elastic deformation through perfect plastic yielding to linear hardening
conditions.

##  1 Introduction

The study of plasto-elastic behaviour is a field of research concerned with
the challenges associated with modelling the deformation of materials that
exhibit both elastic and plastic behaviour under applied loads [1, 2]. These
materials, which are widely used in engineering and manufacturing, initially
exhibit linear elastic behaviour when subjected to mechanical stress, but
begin to deform plastically as soon as a threshold yield stress is exceeded.
Accurate prediction of this transition is important for the design and
analysis of structural components subjected to dynamic loads such as
earthquakes, for modelling the behaviour of soils and rocks under construction
loads or natural forces, and for additive manufacturing processes where
materials are often subjected to complex loads, to name but a few
applications.

Elastic deformations are usually modelled using the Navier-Cauchy equation. As
soon as the material begins yielding (i.e. the stress locally exceeds the
yield stress), the Navier-Cauchy solution is considered as an initial guess in
an iterative algorithm that adjusts the solution locally. Traditionally, the
finite element method (FEM) is used to solve elasto-plastic problems [3, 4, 5,
6]. Although FEM is a powerful tool that offers a mature and versatile
solution approach including all types of adaptivities [7], well-understood
error indicators [8] and an established toolset for isogeometric analysis,
i.e. coupling with computer-aided design (CAD) [9], the complexity of mesh
generation has led researchers to explore alternative methods. In FEM, the
domain is discretised by dividing a complex geometry into smaller, simpler
parts, called elements, that cover it enitely. Despite significant
developments in the field of mesh generation, the process of meshing often
remains the most time-consuming part of the overall solution procedure, while
the quality of the meshing has a direct impact on the accuracy, convergence
and computational cost of the method [10].

In contrast, meshless methods work with node clouds, which are often referred
to as scattered nodes and do not require a topological relationship between
them. Although some authors argued that arbitrary nodes could be used [11,
12], nowadays it is generally accepted that for a stable meshless
approximation the domain should be discretised with quasi-uniformly [13, 14]
distributed scattered nodes. Modern node positioning algorithms can
automatically populate complex geometries, are dimension-independent, support
variable node densities and can handle domains whose boundaries are defined by
CAD objects [15, 16, 17, 18, 19, 20], and yet are significantly less complex
than meshing [15].

The meshless methods have already been applied in elasto-plastic problems. In
[21] authors demonstrated the Reproducing Kernel Particle Method weak solution
of Drucker-Prager and Mohr-Coulomb models. Another weak form solution was
presented in [22], where authors demonstrated the Element free Galerkin method
in elasto-plastic stress analysis around the tip of a crack. In [23] a Radial
Point Interpolation method was demonstrated in elasto-plastic analysis of
frame structures.

In weak form methods, the governing equations of plasticity are formulated in
an integral form rather than a differential form. The key difference from
strong form meshless methods, which directly discretize the differential
equations at nodes, is that weak form methods satisfy these equations in an
average sense over the domain, hence potentially increasing the solution’s
smoothness and numerical stability. However, a recent study revealed that
strong form solution might perform better in capturing the peaks in the stress
comparing it to the weak from solutions [24].

In a strong form, elasto-plastic deformation was approached with the meshless
method of fundamental solutions (MFS) in the simulation of material behaviour
that hardens with plastic deformation characterized by the Chakrabarty model
[25]. In the RBF-FD context, the problem was recently addressed in the
analysis of cantilever beam and the Reinforced Concrete column-Steel [26]. In
[26] authors used a multiquadratic basis that requires shape tuning, which can
have a notable effect on the accuracy of the numerical solution due to
stagnation errors and on the stability of the solution procedure [27, 28].
Solution presented in [26] also stabilises the calculation of the first
derivatives using Finite difference method (FDM) on bilinear local nodes.
Since the study is limited to regular nodes, authors only discussed how to
extend the approach to scattered nodes by using fictitious nodes at which the
field nodes should be interpolated. This has been properly implemented in the
recent paper [29], where authors argued that direct RBF-FD, i.e. without using
FDM to compute the divergence operator, is unsuitable for solving elasto-
plastic problems. Furthermore, authors introduced a stabilisation of the
boundary conditions by adding the first two layers of nodes inside the domain
on the opposite side of the outward-facing normal vector instead of using
purely scattered nodes. They have shown that such stabilisation is essential
for the convergence of a method. Although these stabilisations improve the
stability of the method, they also reduce to some extent the elegance of the
meshless discretisation of the domain and introduce new free parameters that
need to be tuned.

In this paper, we present an alternative approach that does not require any
kind of stabilisation, i.e. a direct RBF-FD can be applied to purely scattered
nodes. In [29], the authors solve the system of nonlinear equations arising
from plastic deformation using the Newton–Raphson iteration. Here, we propose
to solve the problem using the Picard iteration, a much simpler approach where
there is no need to assembly the Jacobian matrix and only elastic material
parameters are used, resulting in a constant stiffness matrix during the
iteration. Although the Picard iteration converges significantly slower than
the Newton–Raphson iteration [30], it is less computationally expensive per
iteration and, as shown in this paper, converges without stabilisation of the
divergence operator or boundary conditions. The application of Picard
iteration in field of elasto-plastic materials was initially proposed by
Kołodziej et al. [31], where authors considered elasto–plastic properties of
prismatic bars using the MFS. The use of Picard iteration was also recently
reported in elasto-plastic torsion problem [32] using the MFS as well as the
Taylor series expansion with moving least squares in [33].

In this paper, we discuss a von Mises plasticity model with non-linear
isotropic hardening under small-strains assumption in a plane stress example
of internally pressurized thick-walled cylinder subjected to varying loading
conditions. The solution procedure is built on the direct RBF-FD approximation
with Polyharmonic Splines (PHS) basis functions and polynomial augmentation
and Picard iteration. The implementation was done using our in-house Medusa
C++ library [34] supporting all the required meshless procedures.

##  2 Iterative algorithm for solving Elasto-Plastic deformation problems

Plastic deformation, also known as _plasticity_ , is the phenomenon when a
shape of a solid material undergo non-reversible changes as a response to the
applied forces, i.e. material’s internal structure is changed during the
deformation. The transition from elastic to plastic behaviour is known as
_yielding_ and the ability to predict it is important for the engineering of
structures and components that are subjected to various loads and exposed to
diverse environmental conditions. Typically, uniaxial tension tests are
conducted to gain insights into the response characteristics of a material
when exposed to external loads. Such tests, when performed with a ductile
material, yield a _stress-strain_ relationship curve as conceptually
illustrated in Figure 1.

Figure 1: Sketch of stress and strain relation as a result of an uniaxial
tension test.

In the first segment, the material undergoes an _elastic deformation_. Since
the stress-strain relation in that segment is linear, it is often referred to
as _linear elastic deformation_ , where the slope of the curve is the elastic
Young’s modulus $E$. If the applied load is released before the _yield stress_
is achieved, the material returns to its original state without permanent
structural changes. On the other hand, if applied loads introduce internal
stresses exceeding the yield stress, the total observed deformation becomes a
sum of two contributions — the plastic part and the elastic part.

The stress-strain relation beyond the _yielding point_ is commonly referred to
as the _hardening curve_. When the material under consideration is unloaded in
that segment, the object no longer returns to its original shape. Instead, at
least parts of it, stay _plastically deformed_. Different materials exhibits
different hardening curves, marked with red in Figure 1. A commonly considered
simplification of the hardening curve is _ideal plasticity_ , where the yield
stress beyond the yielding point is constant (see dashed red line in Figure
1). Another common modelling approach is to use linear hardening curve marked
with red dotted line. Such modelling approach results in a _bilinear_ elastic-
plastic material behaviour, because both elastic and plastic regimes are
governed by a linear stress-strain relation with the Young’s modulus $E$
denoting the slope in the elastic regime and hardening rate $H$ denoting the
slope in the plastic regime. In practice, the hardening is commonly
represented as curve fitted to experimental data obtained from uniaxial
tension tests, in simplest case a piecewise linear fit between the
experimental data resulting in _multilinear material models_.

The study presented in this work is limited to rate-independent materials
where the deformation of the material under consideration does not depend on
the rate of the applied loads. The research is also limited to _isotropic
hardening_ , where the evolution of yield stress surface corresponds to a
uniform (isotropic) expansion of the initial yield surface. A detailed
discussion on other stress-strain relationships and different types of
hardening can be found in [3].

Initially, the deformation is elastic, as depicted in Figure 1. The principal
assumption of linear elasticity is that, when subjected to external forces, a
body deforms linearly – obeying Hooke’s law. Thus, the stress tensor
$\bm{\sigma}$ is defined as the product of the elasticity tensor
$\mathrm{D}^{\mathrm{e}}$ and the strain tensor $\bm{\varepsilon}$,

$\bm{\sigma}=\mathrm{D}^{\mathrm{e}}\bm{\varepsilon}.$ (1)

commonly expressed as
$\mathrm{D}^{\mathrm{e}}=2\mu\mathrm{I}_{S}+\lambda\mathbf{I}\otimes\mathbf{I}$,
where $\mathrm{I}_{S}$ is the symmetric identity tensor, $\mathbf{I}$ is the
identity matrix and $\lambda$ and $\mu$ are the first and second Lamé
parameters, respectively, expressed in terms of Young’s modulus $E$, and
Poisson’s ratio $\nu$ as

$\lambda=\frac{E\nu}{\left(1-2\nu\right)\left(1+\nu\right)}\quad\text{and}\quad\mu=\frac{e}{2\left(1+\nu\right)}.$
(2)

Under small-strain assumption the strain tensor $\bm{\varepsilon}$ and
displacement vector $\vec{u}$ are related as

$\bm{\varepsilon}=\frac{\nabla\vec{u}+(\nabla\vec{u})^{\intercal}}{2}.$ (3)

Considering (1) and (3) and the stationary Cauchy momentum equation in the
absence of the external force ($\nabla\cdot\bm{\sigma}=0$), one gets a well-
known Navier-Cauchy equation

$\left(\lambda+\mu\right)\nabla(\nabla\cdot\vec{u})+\mu\nabla^{2}\vec{u}=0.$
(4)

describing the linear-elastic deformation. Two main types of boundary
conditions are frequently employed: essential boundary conditions and traction
boundary conditions, also referred to as natural boundary conditions.
Essential boundary conditions are used to define displacements $\vec{u}_{0}$
along certain parts of the domain’s boundary, denoted as
$\vec{u}=\vec{u}_{0}$. On the other hand, traction boundary conditions are
concerned with defining surface traction $\vec{t}_{0}$, expressed as
$\bm{\sigma}\vec{n}=\vec{t}_{0}$, were, $\vec{n}$ represents the outward unit
normal to the domain’s boundary.

Upon surpassing a critical stress threshold (_yield stress_ $\sigma_{y}$), the
material experiences localized irreversible plastic deformation. In the
plastic regime, the total local deformation has two contributions. One due to
the elastic strain $\bm{\varepsilon}^{\mathrm{e}}$ and another due to the
plastic strain $\bm{\varepsilon}^{p}$. The sum of both is referred to as the
total strain $\bm{\varepsilon}=\bm{\varepsilon}^{e}+\bm{\varepsilon}^{p}$. In
plastic regime the Navier-Cauchy (4) is therefore not valid anymore. To check
its validity different criteria as the Tresca, the Mohr-Coulomb, the von Mises
yield or the Drucker-Prager criteria [3] can be used. In this paper we focus
on the von Mises yield criterion assuming that the yield stress ($\sigma_{y}$)
is a function of the scalar equivalent plastic strain
$\overline{\varepsilon}^{p}$ (see Figure 14 for example of such a stress-
strain relationship). The criterion is based on the the yield function $\Phi$
defined as

$\Phi={\sigma_{\mathrm{VMS}}}_{n}-\sigma_{y}(\overline{\varepsilon}^{p}_{n}),$
(5)

where von Mises stress $\sigma_{\mathrm{VMS}}$ stands for

$\begin{split}\sigma_{\mathrm{VMS}}&=\frac{1}{\sqrt{2}}\left[\left(\sigma_{xx}-\sigma_{yy}\right)^{2}+\left(\sigma_{yy}-\sigma_{zz}\right)^{2}+\left(\sigma_{zz}-\sigma_{xx}\right)^{2}\right.\\\
&\left.\phantom{\left(\sigma_{xx}-\sigma_{yy}\right)^{2}}\quad+6\left(\sigma_{xy}^{2}+\sigma_{yz}^{2}+\sigma_{xz}^{2}\right)\right]^{\frac{1}{2}}.\end{split}$
(6)

If $\Phi>0$, the material is locally under plastic deformation, otherwise the
deformation is elastic

$\left\\{\begin{matrix}\Phi\leq 0,&\quad\text{(elastic deformation), solution
is locally valid}\\\ \Phi>0,&\quad\text{(plastic deformation), solution is
locally invalid.}\end{matrix}\right.$ (7)

In the elastic regime, the Navier-Cauchy solution is accepted, while in the
plastic regime a correction is applied that adjusts the calculated stresses
and strains to account for plastic effects.

Since the state of stress reflects the history of applied stress, the problem
can be solved in incremental way by applying partial loads to the body

that eventually add-up to the total applied load [3, 30]. The key concept
behind such approach is that every partial load applied to the body is assumed
to result in a purely elastic deformation. Thus, a first estimate is obtained
by solving the Navier-Cauchy Equation (4) resulting in a vector field of
displacement increments $\delta\bm{u}$ corresponding to the amount of applied
partial load. In the $i$-th load step, $i\leq N_{\mathrm{load}}$, one obtains
the following set of bulk state predictions

$\displaystyle\vec{u}_{i}$ $\displaystyle=\vec{u}_{i-1}+\delta\vec{u}$ (8)
$\displaystyle\delta\bm{\varepsilon}$
$\displaystyle=\frac{\nabla\left(\delta\vec{u}\right)+(\nabla\left(\delta\vec{u}\right))^{\intercal}}{2}$
$\displaystyle\bm{\varepsilon}^{\mathrm{e}}_{i}$
$\displaystyle=\bm{\varepsilon}^{\mathrm{e}}_{i-1}+\delta\bm{\varepsilon}$
$\displaystyle\bm{\varepsilon}^{p}_{i}$
$\displaystyle=\bm{\varepsilon}^{p}_{i-1}$
$\displaystyle\overline{\varepsilon}^{p}_{i}$
$\displaystyle=\overline{\varepsilon}^{p}_{i-1}$
$\displaystyle\bm{\sigma}_{i}$
$\displaystyle=\mathrm{D}^{\mathrm{e}}\bm{\varepsilon}^{\mathrm{e}}_{i}$

where a scalar value $\overline{\varepsilon}^{p}$ denotes the accumulated
plastic strain, and the stress tensor is calculated only from the elastic
strain. In next step, we check the validity of the stress by assessing the
yield criterion (7). If $\Phi\leq 0$, the solution is valid as the yield
condition has not been violated, and the assumed elastic deformation computed
using (4) holds. If $\Phi>0$, on the other hand, the prediction was incorrect,
thus, solution needs to first be locally adjusted according to the _return
mapping procedure_ , where we use the single-equation return mapping scheme
[3]. For this purpose an incremental plastic multiplier $\Delta\gamma$ is
introduced. The goal of the return mapping procedure is to find such
$\Delta\gamma$, that the following equality is satisfied [3]

$\Phi={\sigma_{\mathrm{VMS}}}_{i}-3\mu\Delta\gamma-\sigma_{y}(\overline{\varepsilon}^{p}_{i}+\Delta\gamma)=0.$
(9)

We use a Newton-Raphson approach (see Algorithm 1 for details) to solve (9)
and finally, with the computed $\Delta\gamma$, the local bulk state variables
are temporarily updated

$\displaystyle\mathbf{s}_{i}^{j+1}$
$\displaystyle=\bigg{(}1-\frac{\Delta\gamma
3\mu}{{\sigma_{\mathrm{VMS}}}_{i}^{j}}\bigg{)}\mathbf{s}_{i}^{j}$ (10)
$\displaystyle\bm{\sigma}_{i}^{j+1}$
$\displaystyle=\mathbf{s}_{i+1}^{j+1}+p_{i+1}^{j+1}\mathbf{I}$
$\displaystyle\bm{\varepsilon}^{\mathrm{e}\ j+1}_{i}$
$\displaystyle=\frac{1}{2\mu}\mathbf{s}_{i+1}^{j+1}+\frac{1}{3}\varepsilon^{\mathrm{e}\
j}_{v\ i}\mathbf{I}$ $\displaystyle\overline{\varepsilon}^{p\ j+1}_{i}$
$\displaystyle=\overline{\varepsilon}^{p\ j}_{i}+\Delta\gamma$

where $\mathbf{s}$ and $p$ denote the deviatoric and hydrostatic components of
the stress tensor and
$\varepsilon^{\mathrm{e}}_{v}=\sum\mathrm{diag}(\bm{\varepsilon}^{\mathrm{e}})$
is the volumetric component of the elastic strain and $j$ denotes the global
optimisation index (explained in following few sentences). The updated stress
$\bm{\sigma}_{i}^{j+1}$ does not fulfill Navier-Cauchy equation anymore, i.e.

$\nabla\cdot\bm{\sigma}_{i}^{j+1}=-\delta\vec{r},$ (11)

where vector $\vec{r}\neq 0$ stands for residual. To drive the state towards
correct solution a global optimisation is performed, where the whole process
repeats using (11) to calculate stress guess until

$\left\|\delta\vec{r}\right\|<\epsilon_{\mathrm{tol}}$ (12)

holds in all nodes with $\epsilon_{\mathrm{tol}}$ standing for prescribed
tolerance [30]. Note, that during global optimization, the boundary conditions
need to be updated accordingly, e.g., at the $i$-th load step, the total
prescribed stress at the boundary $\Omega$ is

$\bm{\sigma}(\Omega)=\bm{\sigma}_{\mathrm{load}\
i}=\sum_{i}\delta\bm{\sigma},$ (13)

and after $j$ global optimization steps, this boundary condition reads

$\bm{\sigma}(\Omega)=\bm{\sigma}_{\mathrm{load}\
i}-\bm{\sigma}(\Omega)_{i}^{j}.$ (14)

Once the current load step converges, a new partial load is applied, and the
whole process is repeated until the sum of partial loads equals the prescribed
total load.

Algorithm 1 Return mapping algorithm.

Input: Material properties $\Pi$, maximum return mapping iterations
$I_{\mathrm{max}}$ and return mapping tolerance $\epsilon_{\mathrm{tol}}$.
Output: Incremental plastic multiplier $\Delta\gamma$.

1:function return_mapping($\Pi,I_{\mathrm{max}},\epsilon_{\mathrm{tol}}$)

2: $k\leftarrow 0$ $\triangleright$ Initializing iteration counter.

3:
$\Phi^{\star}\leftarrow{\sigma_{\mathrm{VMS}}}_{n}-\sigma_{y}(\overline{\varepsilon}^{p}_{n})$
$\triangleright$ Initializing residual yield function value.

4: // Newton-Raphson iteration:

5: do

6:
$H\leftarrow\frac{\mathrm{d}\sigma_{y}}{\mathrm{d}\overline{\varepsilon}^{p}}\Big{\rvert}_{\overline{\varepsilon}^{p}_{n}+\Delta\gamma}$
$\triangleright$ Hardening curve slope evaluated at
$\overline{\varepsilon}^{p}_{n}+\Delta\gamma$.

7: $\Delta\gamma\leftarrow\Delta\gamma+\frac{\Phi^{\star}}{3\mu+H}$
$\triangleright$ New guess for $\Delta\gamma$.

8:
$\Phi^{\star}={\sigma_{\mathrm{VMS}}}_{n}-3\mu\Delta\gamma-\sigma_{y}(\overline{\varepsilon}^{p}_{n}+\Delta\gamma)$
$\triangleright$ Re-evaluate residual.

9: $k=k+1$ $\triangleright$ Increase number of Newton-Raphson iterations.

10: while $\lvert\Phi^{\star}\rvert>\epsilon_{\mathrm{tol}}\text{ and
}k<k_{\mathrm{max}}$

11: return $\Delta\gamma$ $\triangleright$ Return incremental plastic
multiplier $\Delta\gamma$.

12:end function

The backbone of above solution procedure is the numerical treatment of Navier-
Cauchy (4) PDE that requires appropriate domain discretisation, partial
differential operators discretisation and ultimately solution of global sparse
system resulting from the discretisation of the (4). All these steps will be
discussed in details in following section.

The entire solution procedure demonstrated in Algorithm 2.

Algorithm 2 Solution procedure algorithm.

Input: The problem, nodal density function $h$, stencil size $n$, linear
differential operators $\mathcal{L}$, approximation basis $\xi$, material
properties $\Pi$, number of load steps $N_{\mathrm{load}}$, external load
$\sigma_{\mathrm{load}}$, maximum return mapping iterations $I_{\mathrm{max}}$
and return mapping tolerance $\epsilon_{\mathrm{tol}}$.
Output: Deformation field ${\mathbf{u}}$ and corresponding stress and strain
tensors $\bm{\sigma}$ and $\bm{\varepsilon}$ respectively.

1:function
solve($\mathrm{problem},h,n,\xi,\Pi,N_{\mathrm{load}},\sigma_{\mathrm{load}},I_{\mathrm{max}},\epsilon_{\mathrm{tol}}$)

2: $\Omega\leftarrow\textsc{discretise\\_domain}(h)$ $\triangleright$
Discretize the domain.

3: $\mathcal{L}\leftarrow\textsc{aproximate}(\mathrm{problem},\xi,\Omega,n)$
$\triangleright$ Differential operator approximation.

4: for $i$ $\leftarrow$ 0 to $N_{\mathrm{load}}$ do

5: $\delta\vec{u}=\mathbf{0}$ $\triangleright$ Initalize displacement changes.

6: $\delta\vec{r}=\mathbf{0}$ $\triangleright$ Initalize residuum force
changes.

7:

8: do

9:
$\delta\vec{u}\leftarrow\textsc{solve\\_elastic}(\text{problem},\Pi,\delta\vec{r},\partial\sigma_{\mathrm{load}})$
$\triangleright$ Linear-elastic guess for partial load
$\partial\sigma_{\mathrm{load}}=\frac{\sigma_{\mathrm{load}}}{N_{\mathrm{load}}}$.

10: $\vec{u}=\vec{u}+\delta\vec{u}$

11: $\bm{\varepsilon},\bm{\sigma}\leftarrow\textsc{process}(\vec{u},\Pi)$
$\triangleright$ Compute strain and stress tensors.

12:

13: for each $\vec{p}$ in $\Omega$ do

14: if $\sigma_{\mathrm{VMS}}(\vec{p})>\sigma_{\mathrm{yield}}$ then

15:
$\bm{\sigma}(\vec{p}),\bm{\varepsilon}(\vec{p})\leftarrow\textsc{return\\_mapping}(I_{\mathrm{max}},\epsilon_{\mathrm{tol}},\Pi)$

16: end if

17: end for

18: $\delta\vec{r}\leftarrow\textsc{compute\\_residuum}(\bm{\sigma})$
$\triangleright$ Estimate residuum forces.

19: while $\textsc{max}(\left\|\delta\vec{r}\right\|)>\epsilon_{\mathrm{tol}}$
$\triangleright$ Repeat until residuum forces are sufficiently small.

20: end for

21: return $\bm{\varepsilon},\bm{\sigma},\vec{u}$ $\triangleright$ Return
strains, stresses and displacements for all computational nodes.

22:end function

##  3 Solution to Navier-Cauchy equation using RBF-FD approximation

In this section we will discuss how to numerically solve the Navier-Cauchy
equation (4) using the RBF approximation on scattered nodes. In the first
step, the domain $\Omega$ and its boundary are populated with discrete nodes
at which the solution will be computed. Afterwards, a linear differential
operators governing the equation at hand are approximated at each node
$\mathbf{x}_{i}\in\Omega$. In the third step, the Navier-Cauchy equation is
approximated with system of algebraic equations that is represented as a
global sparse system. Finally, in the last step, the sparse system is solved,
resulting in the numerical solution. In following subsections each of the
above steps are discussed.

###  3.1 Node generation

In RBF community, a well established approach towards stable and accurate RBF
approximation is to discretise domain of interest with quasi-uniformly [13,
14] scattered nodes. Although there is still no generally accepted consensus
on measure for node quality, there are two quantities that are often
considered for that purpose, namely the minimum spacing between any pair of
nodes that is often also referred to as separation distance, and the maximal
empty space without nodes commonly also referred to as fill distance. Here we
use an algorithm that generates nodes according the above guidelines [15] and
has been used several times in RBF related research [35, 36, 37, 38, 39].¸

The algorithm is of an iterative advancing-front type. It operates by starting
with seed nodes that are put in an _expansion queue_. In our particular case,
we initially distribute nodes on the boundary, and since the boundary is one-
dimensional, they can trivially be distributed uniformly (left-most plot of
Figure 4). In each iteration a node is selected from a queue and expanded by
generating new candidate (in our case $12$ candidates) nodes uniformly on a
randomly rotated annulus around the processed point. Note that for different
randomizations, i.e. different random _seeds_ , of the annulus rotations we
get different node layouts that all conform to the quasi-uniformity
guidelines. These candidates are then processed, discarding those outside the
domain or too close to existing nodes. In last iteration step, the nodes that
were not discarded are put in a queue. The iteration continues in a breadth-
first-search manner until the queue is empty, ensuring the domain is fully
covered with the nodes. The algorithm is dimension-independent and capable of
populating complex geometries. It also allows for variable internodal
distances that can be used as a cornerstone for $h$ [40] or $hp$ [36]
adaptivities. It has been also promoted to generation of nodes on parametric
surfaces [16] as well as domains whose boundaries are represented as CAD
models [20], and has been tailored for execution in parallel on shared memory
architectures [41]. Interested readers are referred to the original paper [15]
for more details on the node generation algorithm or its stand-alone C++
implementation in the _Medusa library_ [34].

###  3.2 Approximation of Differential Operators

After computational nodes have been positioned, linear differential operators
of the governing system of PDEs are approximated using a set of nearby nodes,
commonly referred to as _stencil nodes_. There have been different stencil
strategies proposed [42, 43], nevertheless, we resort to the most common
approach using cloud of closest $n$ nodes. Let us assume a linear differential
operator $\mathcal{L}$ in point $\mathbf{x}_{c}\in\Omega$ and its stencil
nodes $\left\\{\mathbf{x}_{c}\right\\}_{i=1}^{n}=\mathcal{N}$. The
approximation is then sough using ansatz

$(\mathcal{L}g)(\vec{x}_{c})\approx\sum_{i=1}^{n}w^{\mathcal{L}}_{i}g(\vec{x}_{i}),$
(15)

for any function $g$ and unknown weights $\mathbf{w}^{\mathcal{L}}$ that are
obtained by enforcing the equality of approximation (15) for a chosen set of
basis functions $\mathbf{f}(\vec{x})$.

$\sum_{i=1}^{n}w_{i}f_{j}(\vec{x_{i}})=(\mathcal{L}f_{j})(\vec{x_{c}}).$ (16)

Gathering all $n$ equations for a given stencil results in a system

$\mathbf{F}\mathbf{w}=\bm{\ell}_{f},\qquad\mathbf{F}=\begin{bmatrix}f(\left\|\vec{x_{1}}-\vec{x_{1}}\right\|)&\cdots&f(\left\|\vec{x_{n}}-\vec{x_{1}}\right\|)\\\
\vdots&\ddots&\vdots\\\
f(\left\|\vec{x_{1}}-\vec{x_{n}}\right\|)&\cdots&f(\left\|\vec{x_{n}}-\vec{x_{n}}\right\|)\\\
\end{bmatrix},\qquad\ell_{f}^{i}=\Big{(}\mathcal{L}f_{i}(\vec{x})\Big{)}\Big{|}_{\vec{x}=\vec{x_{c}}}.$
(17)

In this work, Polyharmonic Splines (PHS) are used as an approximation basis
defined as

$f(r)=\begin{cases}r^{k},&k\text{ odd}\\\ r^{k}\log r,&k\text{
even}\end{cases}$ (18)

for Euclidean distance $r$, the dependency on the shape parameter is
eliminated [44]. It has previously been observed that using only PHS as basis
functions does not guarantee convergent behaviour or solvability [13].
Convergent behaviour and conditional positive definiteness can be assured if
the approximation basis is expanded with monomials [45, 13, 46]. Therefore,
the system (17) is augmented with $N_{p}=\binom{m+d}{d}$ monomials $p$ with
orders up to (and including) degree $m$ for a $d$-dimensional domain, yielding
an overdetermined system of equations describing the constrained optimization
problem [47]. In practice, the weights $\mathbf{w}$ are expressed as a
solution to the

$\begin{bmatrix}\mathbf{F}&\mathbf{P}\\\
\mathbf{P}^{\mathsf{T}}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{w}\\\
\bm{\lambda}\end{bmatrix}=\begin{bmatrix}\bm{\ell}_{f}\\\
\bm{\ell}_{p}\end{bmatrix},$ (19)

with the matrix $\mathbf{P}$ of the evaluated monomials and Lagrangian
multipliers $\mathbf{\lambda}$ that are discarded once the weights are
computed.

The discussed approach is also known as a RBF generated finite differences
method (RBF-FD) that is a popular choice in the meshless community due to its
stability [48, 49].

###  3.3 Navier-Cauchy equation Discretization

In previous sections we discussed how to numerically approximate differential
operators on scattered nodes, in this subsection, we discuss the
discretisation of the Navier-Cauchy equation (4). Since we deal with a two
dimensional domain, we first express the governing problem component-wise
$\vec{u}=(u,v)$ as

$\displaystyle(\lambda+\mu)\frac{\partial}{\partial x}\left(\frac{\partial
u}{\partial x}+\frac{\partial v}{\partial
y}\right)+\mu\left(\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)$ $\displaystyle=0$ (20)
$\displaystyle(\lambda+\mu)\frac{\partial}{\partial y}\left(\frac{\partial
u}{\partial x}+\frac{\partial v}{\partial
y}\right)+\mu\left(\frac{\partial^{2}v}{\partial
x^{2}}+\frac{\partial^{2}v}{\partial y^{2}}\right)$ $\displaystyle=0$ (21)

and natural boundary conditions

$\displaystyle{t_{0}}_{x}$ $\displaystyle=\mu n_{y}\frac{\partial u}{\partial
y}+\lambda n_{x}\frac{\partial v}{\partial
y}+(2\mu+\lambda)n_{x}\frac{\partial u}{\partial x}+\mu n_{y}\frac{\partial
v}{\partial x}$ (22) $\displaystyle{t_{0}}_{y}$ $\displaystyle=\mu
n_{x}\frac{\partial u}{\partial y}+(2\mu+\lambda)n_{y}\frac{\partial
v}{\partial y}+\lambda n_{2}\frac{\partial u}{\partial x}+\mu
n_{x}\frac{\partial v}{\partial x}$ (23)

with surface traction $\vec{t_{0}}=({t_{0}}_{x},{t_{0}}_{x})$ and outward unit
normal to the domain’s boundary $\vec{n}=(n_{x},n_{y})$.

The goal here is to transform the above PDEs into a system of algebraic
equations. To do so, equations (21) and (23) are considered in $N$ nodes,
where corresponding partial differential operators are discretised using (15),
resulting in a global system schematically represented as [50]

$\begin{bmatrix}\mathbf{W}_{11}&\mathbf{W}_{12}\\\
\mathbf{W}_{21}&\mathbf{W}_{22}\end{bmatrix}\begin{bmatrix}\mathbf{u}\\\
\mathbf{v}\end{bmatrix}=\begin{bmatrix}\mathbf{b}_{1}\\\
\mathbf{b}_{2}\end{bmatrix},$ (24)

where $\mathbf{u}$ and $\mathbf{v}$ stand for vector of size $N$ containing
the unknown discrete displacements components. Vectors $\bm{b_{1}}$ and
$\bm{b_{2}}$ hold boundary conditions values and in blocks
$\mathbf{W}_{11},\mathbf{W}_{12},\mathbf{W}_{21},\mathbf{W}_{22}$ stencil
weights (19) are stored. To discretise Navier-Cauchy equation, stencil weights
for approximation of first and second derivatives $\mathbf{w}_{i}^{\partial
x},\mathbf{w}_{i}^{\partial
y},\mathbf{w}_{i}^{\partial^{2}x},\mathbf{w}_{i}^{\partial^{2}y},\mathbf{w}_{i}^{\partial
xy}$ are required in each node $i$ (how to compute those is discussed in
section 3.2). To assemble the system (24), we need also a vector of stencil
nodes of a $i$-th discretisation node $\mathcal{N}(i)$ that is in our case
vector containing $n$ closest nodes. The system is finally expressed as

$\displaystyle\left.\begin{aligned}
W_{11}^{i,\mathcal{N}(i)_{j}}&=\left[(\lambda+2\mu)\mathbf{w}^{\partial^{2}}_{i}+\mu\mathbf{w}^{\partial^{2}y}_{i}\right]_{j}\\\
W_{12}^{i,\mathcal{N}(i)_{j}}&=\left[(\lambda+\mu)\mathbf{w}^{\partial
x\partial y}_{i}\right]_{j}\\\ {\bm{b}_{1}}_{i}&=0\end{aligned}\right\\},$
(25) $\displaystyle\left.\begin{aligned}
W_{21}^{i,\mathcal{N}(i)_{j}}&=\left[(\lambda+\mu)\mathbf{w}^{\partial
x\partial y}_{i}\right]_{j}\\\
W_{22}^{i,\mathcal{N}(i)_{j}}&=\left[\mu\mathbf{w}^{\partial^{2}x}_{i}+(\lambda+2\mu)\bm{\chi}^{\partial^{2}y}_{i}\right]_{j}\\\
{\bm{b}_{2}}_{i}&=0\\\ \end{aligned}\right\\},$ (26)

with $i\in[1,N]$ and $i\in[1,n]$. In the same manner the lines representing
the normal boundary conditions are assembled

$\displaystyle\left.\begin{aligned} W_{11}^{i,\mathcal{N}(i)_{j}}&=\left[\mu
n_{2}\mathbf{w}^{\partial y}_{i}+(2\mu+\lambda)n_{1}\mathbf{w}^{\partial
x}_{i}\right]_{j}\\\ W_{12}^{i,\mathcal{N}(i)_{j}}&=\left[\lambda
n_{1}\mathbf{w}^{\partial y}_{i}+\mu n_{2}\mathbf{w}^{\partial
x}_{i}\right]_{j}\\\ {\bm{b}_{1}}_{i}&=t_{0x,1}\end{aligned}\right\\},$ (27)
$\displaystyle\left.\begin{aligned} W_{21}^{i,\mathcal{N}(i)_{j}}&=\left[\mu
n_{1}\mathbf{w}^{\partial y}_{i}+\lambda n_{2}\mathbf{w}^{\partial
x}_{i}\right]_{j}\\\ W_{22}^{i,\mathcal{N}(i)_{j}}&=\left[\mu
n_{1}\mathbf{w}^{\partial x}_{i}+(2\mu+\lambda)n_{2}\mathbf{w}^{\partial
y}_{i}\right]_{j}\\\ {\bm{b}_{2}}_{i}&=t_{0y,i}\end{aligned}\right\\},$ (28)

while in essential boundary conditions nodes, we simply express

$\displaystyle\begin{aligned} W11_{i,i}&=1\\\
{\bm{b}_{1}}_{i}&=u_{0x,i}\end{aligned}$ and $\displaystyle\begin{aligned}
W22_{i,i}&=1\\\ {\bm{b}_{2}}_{i}&=u_{0y,i}\end{aligned}\quad.$ (29)

The example of assembled system is depicted in Figure 5. Solving this system
provides the numerical solution to the Navier-Cauchy problem as formulated in
Equation (4).

###  3.4 Note on implementation and used parameters

The entire solution procedure is implemented in a C++ 111Source code is
available at https://gitlab.com/e62Lab/public/p_2d_plasticity under tag
_v1.0_. and it is strongly dependent on our in-house a open-source mesh-free
project _Medusa_ [34]. The code was compiled using g++ (GCC) 11.3.0 for Linux
with -O3 -DNDEBUG flags on AMD EPYC 7713P 64-Core Processor computer. In cases
where we show wall-clock times, the CPU frequency was fixed to 2.27 GHz, boost
functionality disabled and CPU affinity assured using the taskset command. The
highest RAM usage during the computation was approximately $4\,\mathrm{GB}$
out of available $12\,\mathrm{GB}$, thus, no RAM-related issues were
encountered due to using the direct SparseLU solver in the process.

In all following analyses, a RBF-FD aproximation with PHS basis of order $k=3$
and augmenting monomials of order $m=3$ on stencil size of $n=20$ is used. The
Picard iteration tolerance is set to $\epsilon_{\mathrm{tol}}=10^{-6}$ The
final sparse linear system, is solved with the direct SparseLU using Eigen C++
library [51]. Post-processing was done with Python 3.10.6 and Jupyter
notebooks, also available in the provided git repository1.

##  4 Numerical Examples

The proposed solution procedure is applied to three different cases:

  1. 1.

the linear-elastic case in Section 4.1,

  2. 2.

the perfectly plastic elasto-plastic case in Section 4.2 and

  3. 3.

the elasto-plastic case with linear isotropic hardening in Section 4.3.

The first two problems have been extensively studied by other researchers. De
Souza et al. [3] used FEM to obtain numerical solutions while Hill [2] derived
the analytic solutions. Both analytical and FEM solutions are used here to
evaluate the accuracy of present solution procedure.

The computational domain is the same for all three cases: a thick-walled pipe
with inner and outer radii $a=100\,\mathrm{mm}$ and $b=200\,\mathrm{mm}$
respectively. As schematically presented in Figure 2, only a quarter of the
pipe from the first quadrant is considered. The material properties were
constant for all considered cases: Young’s modulus $E=210\,\mathrm{GPa}$,
Poisson’s ratio $\nu=0.3$ and yield stress $\sigma_{y}=0.24\,\mathrm{GPa}$.

Figure 2: Schematic representation of a pressurized thick-walled pipe
expansion.

The left side of the pipe is constrained to move freely along the
$y$-direction, while the bottom side of the pipe is constrained to move only
in the $x$-direction. The outer wall of the pipe is traction free, while the
inner wall is subjected to internal pressure $p$, which, as explained in
Section 2, is gradually increased until the full magnitude is applied to the
pipe’s inner surface, i,e, at each load step $p/N_{\mathrm{load}}$ is applied.
One can find all boundary conditions also on a problem scheme in Figure 2.

The amount of internal pressure depends on the case being considered:
$0.05\,\mathrm{GPa}$ for the linear-elastic case, $0.19\,\mathrm{GPa}$ for
perfectly plastic case and $0.175\,\mathrm{GPa}$ for the case with linear
isotropic hardening. Note that all of the following analyses are carried out
assuming plane strain conditions, which means that the out-of-plane or
$z$-direction strains are zero, i.e.
$\varepsilon_{zz}=\varepsilon_{xz}=\varepsilon_{yz}=0$.

###  4.1 The Linear Elastic Case

To establish the basic confidence level in the proposed solution procedure,
the internal pressure $p$ is set to a value of $0.05\,\mathrm{GPa}$, which is
intentionally low, to assure that all deformations are well within the elastic
regime. The case can be solved in a closed form [3, 2]

$\begin{aligned}
\sigma_{r}&=-\left(\frac{\frac{b^{2}}{r^{2}}-1}{\frac{b^{2}}{a^{2}}-1}\right)p,\\\
\sigma_{\theta}&=\left(\frac{\frac{b^{2}}{r^{2}}+1}{\frac{b^{2}}{a^{2}}-1}\right)p,\end{aligned}\qquad\qquad\begin{aligned}
\sigma_{z}&=\frac{2\nu}{\frac{b^{2}}{a^{2}}-1}p,\\\
u_{r}&=\frac{p}{e}\frac{\left(1+\nu\right)\left(1-2\nu\right)r+\left(1+\nu\right)\frac{b^{2}}{r}}{\frac{b^{2}}{a^{2}}-1}\end{aligned}$
(30)

in cylindrical coordinates, which we will use for a proper verification of the
Navier-Cauchy solver. To verify the performance of the proposed solution
procedure using scattered nodes, we numerically solve the problem in Cartesian
coordinate system (left plot of Figure 4). In Figure 3) the RBF-FD solution is
compared with the (30).

Figure 3: Assessment of the present solution of the purely elastic case. The
scattered plot is displaying the relative error in all nodes.

We can see that the present method is in good agreement with the closed-form
solution, with the relative error of the numerical solution well below 1
$\permil$ using $N\approx 95\,000$ nodes. In right plot of Figure 4, we vary
the internodal distance $h$, examining its effect on the numerical solution’s
accuracy as the local field description improves. The error
$\left\|{\mathbf{u}}(\mathbf{x}_{i})-u_{r}(\mathbf{x}_{i})\mathbf{e}_{r}(\mathbf{x}_{i})\right\|$
is computed for each computational point $\mathbf{x}_{i}\in\Omega$ and the
$\ell^{2}$-norm of the numerical solution error is plotted against the number
of discretization nodes $N$. For clarity, an example spatial distribution of
the displacement magnitude error is shown in the middle of Figure 4. Since
there is practically an infinite number of possible node layouts for a given
$h$, depending on the number and position of the seed nodes, the number of
expansion candidates and different randomizations (random _seed_) in candidate
generation, we also test the spread of the solutions in relation to this
variation. The left plot of Figure 4 shows two different node layouts with
different initial node generation settings. For each data point in the right
plot of Figure 4, we varied the node generation setting and obtained an
ensemble of $10$ solutions as a result. The ensemble converges in both
respects, i.e. the finer the domain discretizations the higher the accuracy
and the lower the spread.

Figure 4: An example of scattered nodes using two different random seeds
(left), spatial distribution of the displacement magnitude error (middle) and
convergence of the $\ell^{2}\mathrm{-norm}$ of the displacement vector for the
linear-elastic case (right).

Note that all solutions in this section have been obtained by immediately
enforcing the full amount of internal pressure $p$, omitting the gradual load
increase over $N_{\mathrm{load}}$ load steps for computation by setting
$N_{\mathrm{load}}=1$. This is justified as only elastic deformations are
expected and the solution procedure never encountered a local violation of the
von Mises yield criterion from Equation (5), i.e. the return mapping Algorithm
1 was never used.

The sparsity patterns of the final large sparse systems are shown in Figure 5.
The top row shows the sparsity of the systems for three different domain
discretization qualities. The bottom row of Figure 5 shows the spectra of the
three matrices. The ratios between the real and imaginary parts of the
eigenvalues are in good agreement with previous studies [35, 27, 44]. Most of
the eigenvalues have an order of magnitude larger real-valued part compared to
the relatively small imaginary parts.

Figure 5: Sparsity of final matrices (top row) and spectrum (bottom row) for
three different discretization qualities.

###  4.2 Perfectly-Plastic Yielding of an Internally Pressurized Thick-Walled
Cylinder

With the established confidence in the implementation of elastic deformation,
the internal pressure $p$ is increased to $0.19\,\mathrm{GPa}$, which is
enough to introduce the irreversible plastic deformations. The pressure $p$ is
gradually increased until a prescribed load is reached, as proposed in
Algorithm 2.

For the present cylindrically symmetrical problem, the plastic yielding starts
at the inner surface ($r=a$) and gradually develops towards the outer surface
($r=b$) in the form of a circular front located at $r=c$. The material
collapses when the plastic front reaches the outer surface ($c=b$). At that
point, the entire domain is plastically deformed and the external load
required to achieve that state is referred to as the _limit load_ 222For the
present dimensions and material properties, the limit load is
$p_{\mathrm{lim}}\approx 0.19209\,\mathrm{GPa}$ as discussed by de Souza [3]..

Analytical solutions for the normal stress tensor components and for the
radial displacement are given in [2, 1]. Hill, in his derivation, assumes the
Tresca yielding criterion, but argues that this can make a fair approximation
for the von Mises yielding criterion [2]. Provided that the cylinder has
started yielding, and the position of the plastic front within the cylinder is
located at $c$, the analytical approximations read

$\displaystyle\begin{aligned} \left.\begin{aligned}
\sigma_{r}&=-\frac{2c^{2}}{\sqrt{3}b^{2}}\left(\frac{b^{2}}{r^{2}}-1\right)\sigma_{y}\\\
\sigma_{\theta}&=\frac{2c^{2}}{\sqrt{3}b^{2}}\left(\frac{b^{2}}{r^{2}}+1\right)\sigma_{y}\\\
\sigma_{z}&=\frac{4\nu
c^{2}}{\sqrt{3}b^{2}}\sigma_{y}\end{aligned}\right\\}c\leq r\leq b,\\\
\end{aligned}\qquad\qquad\begin{aligned} \left.\begin{aligned}
\sigma_{r}&=-\frac{2}{\sqrt{3}}\left(1-\frac{c^{2}}{b^{2}}+\ln\frac{c^{2}}{r^{2}}\right)\sigma_{y}\\\
\sigma_{\theta}&=\frac{2}{\sqrt{3}}\left(1+\frac{c^{2}}{b^{2}}+\ln\frac{c^{2}}{r^{2}}\right)\sigma_{y}\\\
\sigma_{z}&=\frac{4\nu
c^{2}}{\sqrt{3}}\left(\frac{c^{2}}{b^{2}}-\ln\frac{c^{2}}{r^{2}}\right)\sigma_{y}\end{aligned}\right\\}a\leq
r\leq c,\end{aligned}$ (31) $\displaystyle u_{r}$
$\displaystyle=\left(1-\nu\right)\frac{2\sigma_{y}c^{2}}{\sqrt{3}\mu
r}+\left(1-2\nu\right)\frac{\sigma_{r}r}{2\mu},\quad a\leq r\leq b,$

where $c$ is determined from its relation to the cylinder’s internal pressure
$p$ as

$p=\frac{2\sigma_{y}}{\sqrt{3}}\left(1-\frac{c^{2}}{b^{2}}+\ln\frac{c^{2}}{a^{2}}\right).$
(32)

As in the previous case, we use Cartesian coordinates, even though the
solution is given in cylindrical coordinates. In this case, the internal
pressure has been gradually increased over $N_{\mathrm{load}}=25$ load steps.
Example solution for the last load step corresponding to internal pressure
$p=0.19\,\mathrm{GPa}$ is shown in Figure 6. Figure shows the normal stresses
$\sigma_{r}$, $\sigma_{\theta}$ and $\sigma_{z}$ and the von Mises stress
$\sigma_{\mathrm{VMS}}$. Additionally, the plastic front at $c=186.67\
\mathrm{mm}$ is shown with a dashed pink line.

Figure 6: Example of numerical solution for perfectly plastic case at
$p=0.19\,\mathrm{GPa}$, $N=95\,300$ nodes at the final load step
$N_{\mathrm{load}}=25$. The dashed line represents plastic front position
$c=186.67\ \mathrm{mm}$. Figure 7: Displacement magnitude
$\left\|\mathbf{u}\right\|$ for the perfectly plastic case. For clarity, only
every second point is shown for the FEM solution.

The accuracy of the numerical solution is further evaluated in Figures 7 and
8, where the present solution procedure is compared to the analytic solution
(see Equations (31)) and to the FEM solution. To obtain the FEM solution,
commercial software package Abaqus was used. The FEM domain was discretized
with a structured mesh consisting of 23 600 8-node biquadratic plane strain
quadrilaterals (CPE8) with $71\,473$ nodes, corresponding to element side
length of approximately $1\,\mathrm{mm}$.

In Figure 7 the displacement magnitude with respect to the radial coordinate
$r$ is shown for all $N=95\,300$ nodes of the present method and for nodes
along the $y=0$ line for the FEM approximation. We are pleased to observe that
the proposed solution procedure is in good agreement with the FEM solution
throughout the entire domain. Both FEM and RBF-FD solutions, however, are
slightly offset from the analytical solution. The most likely reason for the
mismatch compared to the analytical solution is the different yield criteria
used.

Figure 8: Stress analysis in case of perfectly elasto-plastic problem.

Moreover, in Figure 8, the performance of the present method is evaluated in
terms of stresses. Specifically, the figure shows the normal stress components
and the von Mises stress $\sigma_{\mathrm{VMS}}$ with respect to the radial
coordinate $r$. Good agreement of the present method and the FEM solution is
observed, while both again show some differences compared to the analytic
solution, particularly near the elastic-plastic interface.

####  4.2.1 The effect of load stepping

This section briefly discusses the impact of $N_{\mathrm{load}}$ on solution
accuracy and the convergence of the Picard iteration. In the left plot of
Figure 9, an effect of $N_{\mathrm{load}}$ on the final displacement
accompanied with the average number of Picard iterations per step is
presented. We see that in this particular case $N_{\mathrm{load}}$ has no
effect on the final solution, i.e. the maximum displacement is almost constant
with respect to the $N_{\mathrm{load}}$ (note the scale of the left $y$-axis).
On the other hand, increasing the $N_{\mathrm{load}}$ leads to additional
computational effort. Although the average number of required Picard
iterations per step decreases with increasing $N_{\mathrm{load}}$, the
cumulative number of iterations increases. The right plot of Figure 9 shows
the Picard convergence for a setup with $N_{\mathrm{load}}=10$ in the sixth
step, which is similar to the one presented in [30].

Figure 9: The effect of number of load steps $N_{\mathrm{load}}$ on the
accuracy of the numerical solution in terms of maximum displacement magnitude
$\left\|\mathbf{u}\right\|$ and number of Picard iteration per step (left).
Example of residual convergence at sixth load step in simulation with ten load
steps (right).

However, it is important to consider that, for some problems, the degree of
non-linearity of the return-mapping equations is so high that the convergence
radius of iteration scheme becomes reduced [3]. In such cases, improved
initial guesses or a increased number of load steps is recommended [3, 52].
Some researchers also proposed more sophisticated algorithms [3, 52, 53] in
order to improve the convergence.

####  4.2.2 The Plastic Front Position

The following set of analyses is devoted to study the behaviour of hoop
($\sigma_{\theta}$) and axial ($\sigma_{z}$) stresses in Figure 8. Neither of
the two numerical solutions is able to accurately follow the analytic
solution. Nevertheless, a zoomed in section of the hoop stress shows the
breaking point, positioned at the elastic-plastic interface, hints that the
two conceptually different numerical solutions are in good agreement. It
appears as if both numerical approaches overestimate the position of the
plastic front by approximately $4\,\mathrm{mm}$ compared to the analytic
solution (see also Table 2). A possible reason for this observed discrepancy
between analytic and numerical solutions may be in the derivation of the
analytic solution, originally derived for the Tresca yielding criterion. The
Tresca yielding theory333From the design point of view, the Tresca yield
condition is safer, but it also leads the engineer to take unnecessary
measures in attempt to prevent an unlikely failure. Experiments suggest that
the von Mises yield criterion is the one which provides better agreement with
observed behaviour than the Tresca yield criterion. However, the Tresca yield
criterion is still used because of its mathematical simplicity. is known to be
more conservative than the von Mises theory [54].

The position of plastic front $c$ is estimated from the numerical solution
$\widehat{u}$. Example demonstrating how the elastic-plastic interface has
been determined is shown in Figure 10. Red and blue points are used to mark
the plastic and elastic regime, respectively. To capture the front position we
fitted following function to the computed data

$f_{\mathrm{fit}}(r)=C_{1}r^{-1}+C_{2}r^{-2}+C_{3}\log(r^{-1})+C_{4}$ (33)

with constants $C_{i}$ (provided in Table 1) to fit the data and a full black
line to denote them in Figure 10. The intersection of the two fitted lines,
gives us an estimate on the position of the plastic front $c$. For clarity,
only analysis for the present solution procedure is shown in Figure 10 while
Table 2 also shows the numerical data for the FEM solution.

Table 1: Constant values to reproduce the fit lines from Figure 10. regime | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}$
---|---|---|---|---
plastic | $-9.015\cdot 10^{1}$ | $2.751\cdot 10^{3}$ | $8.877\cdot 10^{-2}$ | $1.123$
elastic | $3.036\cdot 10^{2}$ | $-9.989\cdot 10^{3}$ | $-7.823\cdot 10^{-1}$ | $-5.176$
Figure 10: Example of plastic-elastic interface determination. Blue dots
represent the points where the accumulated plastic strain
$\bar{\varepsilon}^{p}=0$, and red dots $\bar{\varepsilon}^{p}>0$.

Table 2 clearly shows that the two numerical solutions are in good agreement
in terms of all studied variables. The von Mises stress in both extremal parts
of the domain, i.e. $r=a$ and $r=b$ are in good agreement. The displacement
magnitude $\left\|\mathbf{u}\right\|$ shows a marginal difference but the
plastic front position $c$ again differs for only one tenth of a millimetre.
As previously observed in Figure 8, the closed-form solution predicts the
elastic-plastic interface closer to the origin by approximately
$4\,\mathrm{mm}$.

Table 2: Results of the perfect plastic case - von Mises stress, and
displacement magnitude at inner and outer surfaces, and position of the
elastic-plastic interface. Under “Present” results we present the solution
with $N\approx 95\,000$ vs. FEM vs. analytical approximation.

[b]  $\sigma_{\mathrm{VMS}}$ [GPa] $\left\|\mathbf{u}\right\|$ [mm] $c$ [mm]
$N$ $r=a$ $r=b$ $r=a$ $r=b$ Present 0.2400✳ 0.2106✳ 0.3623✳ 0.2053✳ 186.67✩
$95\,300$ FEM 0.2400 0.2109 0.3593 0.2057 186.79✩ $71\,473$ Analytical
0.2400$\triangle$ 0.2060 0.3546 0.2008 182.89 -

  * ✳

Average for all nodes on the edge.

  * ✩

Determined by fitting lines through points that are in the plastic zone
($\varepsilon^{p}>0$), and points that are outside it ($\varepsilon^{p}=0$),
and finding the intersection between the two fits.

  * $\triangle$

Presumed von Mises stress criterion.

####  4.2.3 Shape of the Elastic-Plastic Front

Due to the cylindrical symmetry of the governing problem, the obtained
numerical solution is expected to show symmetrical properties. For instance,
the elastic-plastic front shape is expected to preserve a circular shape.
Therefore, in this section, the front positioning algorithm from Section 4.2.2
is extended and applied to $N_{\mathrm{seg}}$ segments, which equally divide
the computational domain. A schematic presentation of the domain division is
shown in Figure 11 on the right, where a coarse domain discretization with
$N=6078$ is shown for clarity.

To determine the front position within each segment, the front position is
estimated using the same procedure as in Section 4.2.2, except that only nodes
within a given segment are used compared to the algorithm in Section 4.2.2
working on all computational nodes.

Figure 11: The elastic-plastic front distance to the origin computed for each
of the $N_{\mathrm{seg}}=20$ segments.

The shape of the front has been computed for the finest domain discretization.
Results for $N_{\mathrm{seg}}=20$ segments with approximately $5000$ nodes
within each segment are shown in Figure 11 in the centre. For reference, the
fully-symmetrical FEM solution is added with a full black line. The numerical
solution obtained with the present approach indeed preserves the circular
shape of the elastic-plastic front position. The circular shape is further
evaluated in Figure 11 on the right, where the distance of the plastic front
from the origin $r_{\mathrm{present}}$ is compared to distance
$r_{\mathrm{FEM}}$ of the FEM-based solution.

###  4.3 Linearly-Hardening Yielding of an Internally Pressurized Thick-Walled
Cylinder

In this section, we consider the case of isotropic linear hardening in the
plastic regime (see dotted red line in Figure 1). Although pure linear
hardening is hardly observed in real materials, it forms strong foundations
for the simulation of realistic materials. The hardening curve for real
materials is typically measured experimentally during an uniaxial tension
test, where the experimental data is interpolated with a piecewise-linear
function to create a model that can be used in numerical simulation.

Figure 12: Displacement magnitude $\left\|\mathbf{u}\right\|$ vs radial
coordinate $r$. Present solution compared with FEM. Note that the FEM solution
is only displaying every second point for clarity.

Specifically, in our case, the hardening curve is defined as follows

$\sigma_{y}(\overline{\varepsilon}^{p})=0.24\,\mathrm{MPa}+10\,\mathrm{MPa}\cdot\overline{\varepsilon}^{p},$
(34)

resulting in a constant hardening slope $H=10\,\mathrm{MPa}$ throughout the
entire plastic regime.

Example solution fields obtained with the present solution procedure are
similar to the ones shown in Figure 6 and are, thus, not shown again. Instead,
the numerical solution obtained with the present method is compared to the
community standard – the FEM solution from Abaqus. The Abaqus setup was the
same as in the previous subsection. The two solutions are compared in terms of
the displacement field magnitude with respect to the radial coordinate in
Figure 12. Good agreement between the two conceptually different approaches is
observed throughout the entire domain in radial coordinates. These
observations are further confirmed with numerical data in Table 3.

In the next step, we evaluate the accuracy of the numerical solution in terms
of the secondary variables in Figure 13. Note that, the inside plot of Hoop
stress $\sigma_{\theta}$ in Figure 13 is a close-up to the breaking point,
which is positioned at the plastic-elastic interface. All secondary variables
are also in good agreement with the FEM-based solution procedure.

Figure 13: Stress analysis in case of isotropic linear hardening.

A summary of the performance evaluation can be found in Table 3. We show the
von Mises stresses and displacement magnitudes at both extremal values of the
domain, i.e. $r=a$ and $r=b$. Additionally, the plastic front position $c$ is
shown. Note that the same procedure as explained in the previous Section 4.2.2
has been used to determine the plastic-elastic interface position $c$, with
constants $C_{i}$ gathered in Table 4.

Table 3: Results of the linear-hardening plastic case - von Mises stress, and
displacement magnitude at inner and outer surfaces, and position of the
elastic-plastic interface.

[b]  $\sigma_{\mathrm{VMS}}$ [GPa] $\left\|\mathbf{u}\right\|$ [mm] $c$ [mm]
$N$ $r=a$ $r=b$ $r=a$ $r=b$ Present 0.25589✳ 0.13635✳ 0.22212✳ 0.13295✳
149.42✩ $95\,346$ FEM 0.25592 0.13634 0.22224 0.13294 149.40✩ $71\,473$
$\frac{\widehat{u}-\widehat{u}_{\mathrm{FEM}}}{\widehat{u}_{\mathrm{FEM}}}\cdot
10^{3}$ -0.117 0.073 -0.540 0.075 0.134 -

  * ✳

Average for all nodes on the edge.

  * ✩

Determined numerically as explaneed in Section 4.2.2.

Table 4: Constant values $C_{i}$ used to determine the plastic front position $c$ in case of linear hardening. regime | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}$
---|---|---|---|---
plastic | $-8.120\cdot 10^{1}$ | $2.426\cdot 10^{3}$ | $1.016\cdot 10^{-1}$ | $1.157$
elastic | $2.993\cdot 10^{-1}$ | $3.055\cdot 10^{3}$ | $-8.508\cdot 10^{-4}$ | $7.102\cdot 10^{-2}$

##  5 Additional Example

In addition to already solved cases, we demonstrate a solution for irregular
domain with a more realistic stress-strain relation. Additionally, the domain
under consideration is cut with the following four spheres:

$\displaystyle B_{0}$ $\displaystyle=\left\\{\mathbf{x}\in\mathbb{R}^{2},\
\left\|\mathbf{x}-\frac{a+b}{2}\Biggl{(}\cos\biggl{(}\frac{\pi}{16}\biggl{)},\sin\biggl{(}\frac{\pi}{16}\biggl{)}\Biggl{)}\right\|<20\right\\},$
(35) $\displaystyle B_{1}$
$\displaystyle=\left\\{\mathbf{x}\in\mathbb{R}^{2},\
\left\|\mathbf{x}-\frac{a+b}{2}\Biggl{(}\cos\biggl{(}\frac{7\pi}{16}\biggl{)},\sin\biggl{(}\frac{7\pi}{16}\biggl{)}\Biggl{)}\right\|<10\right\\},$
(36) $\displaystyle B_{2}$
$\displaystyle=\left\\{\mathbf{x}\in\mathbb{R}^{2},\
\left\|\mathbf{x}-0.8\cdot
a\Biggl{(}\cos\biggl{(}\frac{\pi}{4}\biggl{)},\sin\biggl{(}\frac{\pi}{4}\biggl{)}\Biggl{)}\right\|<30\right\\},$
(37) $\displaystyle B_{3}$
$\displaystyle=\left\\{\mathbf{x}\in\mathbb{R}^{2},\
\left\|\mathbf{x}-1.1\cdot
b\Biggl{(}\cos\biggl{(}\frac{\pi}{4}\biggl{)},\sin\biggl{(}\frac{\pi}{4}\biggl{)}\Biggl{)}\right\|<50\right\\},$
(38)

resulting in a domain as shown in Figure 15. The edge of $B_{2}$ is subjected
to internal pressure $p=0.13\,\mathrm{MPa}$, while the remaining three sphere
edges are traction free, i.e. $\bm{\sigma}\vec{n}=0$.

For the material properties, we simulate a set of uniaxial test results and
use a piecewise linear interpolation in the solution procedure. The stress-
strain relation is given as provided in Table 5 and illustrated in Figure 14.

$\sigma$ [MPa] | $\overline{\varepsilon}^{p}$
---|---
0.240 | 0.000
0.290 | 0.001
0.330 | 0.002
0.370 | 0.003
0.400 | 0.004
0.430 | 0.005
0.450 | 0.006
0.470 | 0.007
0.483 | 0.008
0.495 | 0.009
0.500 | 0.010
Table 5: Synthetic stress-strain relation data.

Figure 14: Plotted synthetic stress-strain relationship.

The domain is discretized with $20\,292$ scattered nodes corresponding to
approximately 1 mm internodal spacing. Example solution is shown in Figure 15
with displacement magnitudes on the left and accumulated plastic strain
$\overline{\varepsilon}^{p}$ on the right.

Figure 15: Example solution on irregular domain. Displacement magnitude (left)
and accumulated plastic strain (right).

The proposed solution procedure is clearly able to cope with irregular domains
and real-like material properties. The largest plastic deformations are
observed in the bottom right part of the domain, while the largest
deformations are observed in the top right part of the domain, effectively
yielding an ellipse-like deformed shape. Note that no special technique or
algorithm adjustments were required to solve such problem. The solution
procedure used is exactly as described in Algorithm 2.

##  6 Conclusions

This paper discusses the use of meshless methods, in particular RBF-FD, in
solving an elasto-plastic problem using the Picard iteration. The study
confirms that the presented solution method effectively accounts for different
stages of material behaviour, from elastic deformation to perfect plastic
yielding and linear hardening conditions. The method shows good agreement with
analytical and FEM solutions (obtained with the commercial Abaqus solver),
confirming its accuracy and robustness. It shows stable behaviour and provides
results that are independent of the underlying geometry, without the need for
special treatment in the approximation of differential operators, boundary
conditions or spatial discretisation.

The iterative numerical algorithm initially assumes elastic deformation, i.e.
in the first step, displacement and strain fields are calculated by solving
the Navier-Cauchy equation. If the stress at any node exceeds the yield
stress, indicating plastic deformation, a correction is made using a return
mapping algorithm that adjusts the calculated stresses and strains to account
for plastic effects, followed by a Picard iteration to solve for non-linear
deformations.

Although the Picard iteration converges monotonically and slowly, it avoids
the divergence issues common with the Newton-Raphson method, which otherwise
converges abruptly, especially when using the consistent tangent matrix [30].
This was also reported in a recent RBF-FD analysis of elasto-plastic benchmark
problems [29], where authors had to introduce additional stabilisations in the
calculation of the divergence operator and boundary conditions in order to
achieve stable convergence with Newton-Raphson iterations.

In this paper, we show that a direct RBF-FD on pure scattered nodes can be
used with the help of Picard iteration, which simplifies both the
implementation and the discretisation of the underlying domain considerably,
but at the price of slow convergence.

Here, the Navier-Cauchy equation is solved in a strong form using a direct
RBF-FD approximation on pure scattered nodes generated by advancing front
algorithm [15] that can be generalised to higher dimension [35] and supports
$h$-adaptive discretization [40] via spatially variable node spacing. At each
node, the linear differential operators relevant to the problem are
approximated using the RBF-FD approximation. Nearby nodes, known as stencil
nodes, are used to approximate the differential operators based on RBFs
augmented with polynomials. In addition to increased stability and solving
problems due to the Gibss phenomenon [55], polynomial augmentation also gives
us direct control over the order of the method, which can be dynamically
adjusted node by node. This property enables a relatively straightforward
application of $hp$-adaptivity [36].

While the local meshless approach offers considerable potential advantages, it
is not without drawbacks. Despite extensive research, the RBF-FD method still
has some unsolved problems. Node quality in a meshless context is not nearly
as well understood as in FEM, where mesh quality measures such as aspect
ratio, skewness, Jacobian ratio, etc. are well-defined concepts [56]. There
are also no clear guidelines for the selection of stencils. In this work, we
simply used the nearest $n$ nodes as a stencil, which leads to a relatively
computationally intensive method due to the rather high stencil size
reuqirements [45].

Future work could focus on extending this method to three-dimensional problems
and investigating the effects of dynamic loading conditions in more complex
hardening models. As for the approximation method, symmetric stencils [57]
could be used in addition to simple closest neighbourhood stencils, to improve
the execution performance. The presented method is also compatible with the
recently introduced $h$ [40] and $hp$ [36] adaptivities, which could be used
to improve the performance of elasto-plastic deformation simulation. As for
the iteration, an improved Picard iteration could be considered, which has
been recently applied to flow in porous media [58]. The presented solution
could also be coupled with a more sophisticated Newton-Raphson iteration to
alleviate its stability problems [30].

## Acknowledgments

Authors acknowledge the financial support from the Slovenian Research And
Innovation Agency (ARIS) research core funding No. P2-0095.

## References

  * [1] J. Chakrabarty, Theory of Plasticity, McGraw-Hill international editions. Engineering mechanics series, McGraw-Hill, 1987.
  * [2] R. Hill, The Mathematical Theory of Plasticity, Oxford classic texts in the physical sciences, Clarendon Press, 1998.
  * [3] E. A. de Souza Neto, D. Peric, D. R. J. Owen, Computational Methods for Plasticity: Theory and Applications, Wiley, 2008.
  * [4] B. Schröder, D. Kuhl, Small strain plasticity: classical versus multifield formulation, Archive of Applied Mechanics 85 (2015) 1127–1145.
  * [5] A. A. Roostaei, H. Jahed, A cyclic small-strain plasticity model for wrought mg alloys under multiaxial loading: Numerical implementation and validation, International Journal of Mechanical Sciences.
  * [6] G. Y. Amouzou, A. Soulaïmani, Numerical algorithms for elastoplacity: Finite elements code development and implementation of the mohr–coulomb law, Applied Sciences 11 (10).  doi:10.3390/app11104637.
  * [7] W. F. Mitchell, M. A. McClain, A comparison of hp-adaptive strategies for elliptic partial differential equations, ACM Transactions on Mathematical Software (TOMS) 41 (1) (2014) 1–39.
  * [8] K. Segeth, A review of some a posteriori error estimates for adaptive finite element methods, Mathematics and Computers in Simulation 80 (8) (2010) 1589–1600.
  * [9] J. A. Cottrell, T. J. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [10] G.-R. Liu, Y.-T. Gu, An introduction to meshfree methods and their programming, Springer Science & Business Media, 2005.
  * [11] G.-R. Liu, Mesh free methods: moving beyond the finite element method, CRC press, 2002.  doi:10.1201/9781420040586.
  * [12] K. Reuther, B. Sarler, M. Rettenmayr, Solving diffusion problems on an unstructured, amorphous grid by a meshless method, Int. J. Therm. Sci. 51 (2012) 16–22.  doi:10.1016/j.ijthermalsci.2011.08.017.
  * [13] H. Wendland, Scattered data approximation, Vol. 17, Cambridge university press, 2004\.
  * [14] G.-R. Liu, Meshfree methods: moving beyond the finite element method, CRC press, 2009.
  * [15] J. Slak, G. Kosec, On generation of node distributions for meshless pde discretizations, SIAM Journal on Scientific Computing 41 (5) (2019) A3202–A3229.  doi:10.1137/18M1231456.
  * [16] U. Duh, G. Kosec, J. Slak, Fast variable density node generation on parametric surfaces with application to mesh-free methods, SIAM Journal on Scientific Computing 43 (2) (2021) A980–A1000.
  * [17] B. Fornberg, N. Flyer, Fast generation of 2-D node distributions for mesh-free PDE discretizations, Computers & Mathematics with Applications 69 (7) (2015) 531–544.  doi:10.1016/j.camwa.2015.01.009.
  * [18] K. van der Sande, B. Fornberg, Fast variable density 3-d node generation, SIAM Journal on Scientific Computing 43 (1) (2021) A242–A257.
  * [19] V. Shankar, R. M. Kirby, A. L. Fogelson, Robust node generation for meshfree discretizations on irregular domains and surfaces, SIAM J. Sci. Comput. 40 (4) (2018) 2584–2608.  doi:10.1137/17m114090x.
  * [20] U. Duh, V. Shankar, G. Kosec, Discretization of non-uniform rational b-spline (nurbs) models for meshless isogeometric analysis, arXiv preprint arXiv:2303.02638.
  * [21] Z. Ji-fa, Z. Wen-pu, Z. Yao, A meshfree method and its applications to elasto-plastic problems, Journal of Zhejiang University-SCIENCE A 6 (2) (2005) 148–154.
  * [22] M. Kargarnovin, H. E. Toussi, S. Fariborz, Elasto-plastic element-free galerkin method, Computational Mechanics 33 (2004) 206–214.
  * [23] J. Belinha, M. Aires, D. E. Rodrigues, Elastoplastic analysis of frame structures using radial point interpolation meshless methods, Applied Sciences 13 (23).  doi:10.3390/app132312591.
URL https://www.mdpi.com/2076-3417/13/23/12591

  * [24] G. Kosec, J. Slak, M. Depolli, R. Trobec, K. Pereira, S. Tomar, T. Jacquemin, S. P. Bordas, M. A. Wahab, Weak and strong from meshless methods for linear elastic problem under fretting contact conditions, Tribology InternationalPublisher: Elsevier.  doi:10.1016/j.triboint.2019.05.041.
  * [25] M. A. Jankowska, On elastoplastic analysis of some plane stress problems with meshless methods and successive approximations method, Engineering Analysis with Boundary Elements 95 (2018) 12–24.  doi:https://doi.org/10.1016/j.enganabound.2018.06.013.
  * [26] P. Jiang, H. Zheng, J. Xiong, P. Wen, Nonlinear elastic-plastic analysis of reinforced concrete column-steel beam connection by rbf-fd method, Engineering Analysis with Boundary Elements 128 (2021) 188–194.
  * [27] V. Bayona, N. Flyer, B. Fornberg, G. A. Barnett, On the role of polynomials in rbf-fd approximations: Ii. numerical solution of elliptic pdes, Journal of Computational Physics 332 (2017) 257–273.
  * [28] J. Wang, G. Liu, On the optimal shape parameters of radial basis functions used for 2-d meshless methods, Computer methods in applied mechanics and engineering 191 (23-24) (2002) 2611–2630.
  * [29] G. Vuga, B. Mavrič, B. Šarler, An improved local radial basis function method for solving small-strain elasto-plasticity, Computer Methods in Applied Mechanics and Engineering 418 (2024) 116501.  doi:10.1016/j.cma.2023.116501.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782523006254

  * [30] V. M. Yarushina, M. Dabrowski, Y. Y. Podladchikov, An analytical benchmark with combined pressure and shear loading for elastoplastic numerical models, Geochemistry, Geophysics, Geosystems 11 (8).  doi:https://doi.org/10.1029/2010GC003130.
  * [31] J. A. Kolodziej, M. A. Jankowska, M. Mierzwiczak, Meshless methods for the inverse problem related to the determination of elastoplastic properties from the torsional experiment, International Journal of Solids and Structures 50 (25) (2013) 4217–4225.  doi:https://doi.org/10.1016/j.ijsolstr.2013.08.025.
URL https://www.sciencedirect.com/science/article/pii/S0020768313003417

  * [32] F. Moayyedian, J. K. Grabski, Elastic–plastic torsion problem with non-linear hardenings using the method of fundamental solution, Archives of Civil and Mechanical Engineering 21 (4) (2021) 155.
  * [33] B. Xu, R. Zhang, K. Yang, G. Yu, Y. Chen, Application of generalized finite difference method for elastoplastic torsion analysis of prismatic bars, Engineering Analysis with Boundary Elements 146 (2023) 939–950.  doi:https://doi.org/10.1016/j.enganabound.2022.11.028.
URL https://www.sciencedirect.com/science/article/pii/S0955799722004295

  * [34] J. Slak, G. Kosec, Medusa: A c++ library for solving pdes using strong form mesh-free methods, ACM Transactions on Mathematical Software (TOMS) 47 (3) (2021) 1–25.
  * [35] M. Jančič, J. Slak, G. Kosec, Monomial augmentation guidelines for rbf-fd from accuracy versus computational time perspective, Journal of Scientific Computing 87 (1) (2021) 9.
  * [36] M. Jančič, G. Kosec, Strong form mesh-free hp-adaptive solution of linear elasticity problem, Engineering with Computers (2023) 1–21.
  * [37] M. Najafi, M. Dehghan, B. Šarler, G. Kosec, B. Mavrič, Divergence-free meshless local Petrov–Galerkin method for Stokes flow, Engineering with Computers 38 (6) (2022) 5359–5377.  doi:10.1007/s00366-022-01621-w.
URL https://link.springer.com/10.1007/s00366-022-01621-w

  * [38] J. M. Berljavac, P. K. Mishra, J. Slak, G. Kosec, RBF-FD analysis of 2D time-domain acoustic wave propagation in heterogeneous media, Computers & Geosciences 153 (2021) 104796, publisher: Elsevier.
  * [39] M. Jančič, M. Založnik, G. Kosec, Meshless interface tracking for the simulation of dendrite envelope growth, Journal of Computational Physics 507 (2024) 112973.  doi:10.1016/j.jcp.2024.112973.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999124002225

  * [40] J. Slak, G. Kosec, Adaptive radial basis function-generated finite differences method for contact problems, International Journal for Numerical Methods in Engineering 119 (7) (2019) 661–686.  doi:https://doi.org/10.1002/nme.6067.
  * [41] M. Depolli, J. Slak, G. Kosec, Parallel domain discretization algorithm for RBF-FD and other meshless numerical methods for solving PDEs, Computers & Structures 264 (2022) 106773.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045794922000335

  * [42] O. Davydov, D. T. Oanh, Adaptive meshless centres and rbf stencils for poisson equation, Journal of Computational Physics 230 (2) (2011) 287–304.
  * [43] O. Davydov, D. T. Oanh, N. M. Tuong, Improved stencil selection for meshless finite difference methods in 3d, Journal of Computational and Applied Mathematics (2023) 115031.
  * [44] V. Bayona, Comparison of moving least squares and rbf+ poly for interpolation and derivative approximation, Journal of Scientific Computing 81 (2019) 486–512.
  * [45] V. Bayona, An insight into rbf-fd approximations augmented with polynomials, Computers & Mathematics with Applications 77 (9) (2019) 2337–2353.
  * [46] N. Flyer, B. Fornberg, V. Bayona, G. A. Barnett, On the role of polynomials in rbf-fd approximations: I. interpolation and accuracy, Journal of Computational Physics 321 (2016) 21–38.
  * [47] N. Flyer, B. FornXberg, V. Bayona, G. A. Barnett, On the role of polynomials in rbf-fd approximations: I. interpolation and accuracy, Journal of Computational Physics 321 (2016) 21–38.  doi:https://doi.org/10.1016/j.jcp.2016.05.026.
  * [48] S. Le Borne, W. Leinen, Guidelines for rbf-fd discretization: Numerical experiments on the interplay of a multitude of parameter choices, Journal of Scientific Computing 95 (1) (2023) 8.
  * [49] M. Jančič, G. Kosec, Stability analysis of rbf-fd and wls based local strong form meshless methods on scattered nodes, in: 2022 45th Jubilee International Convention on Information, Communication and Electronic Technology (MIPRO), 2022, pp. 275–280.  doi:10.23919/MIPRO55190.2022.9803334.
  * [50] J. Slak, G. Kosec, Refined meshless local strong form solution of cauchy-navier equation on an irregular domain, Engineering Analysis with Boundary Elements 100 (2019) 3–13, improved Localized and Hybrid Meshless Methods - Part 1.  doi:https://doi.org/10.1016/j.enganabound.2018.01.001.
  * [51] G. Guennebaud, B. Jacob, et al., Eigen v3, http://eigen.tuxfamily.org (2010).
  * [52] K. Krabbenhøft, Basic computational plasticity, University of Denmark.
  * [53] A. Perez-Foguet, A. Rodriguez-Ferran, A. Huerta, Consistent tangent matrices for substepping schemes, Computer Methods in Applied Mechanics and Engineering 190 (35) (2001) 4627–4647.  doi:https://doi.org/10.1016/S0045-7825(00)00336-4.
URL https://www.sciencedirect.com/science/article/pii/S0045782500003364

  * [54] A. A. Roostaei, H. Jahed, 2 - fundamentals of cyclic plasticity models, in: H. Jahed, A. A. Roostaei (Eds.), Cyclic Plasticity of Metals, Elsevier Series on Plasticity of Materials, Elsevier, 2022, pp. 23–51.  doi:https://doi.org/10.1016/B978-0-12-819293-1.00011-5.
  * [55] B. Fornberg, N. Flyer, Solving PDEs with radial basis functions, Acta Numerica 24 (2015) 215–258.  doi:10.1017/S0962492914000130.
URL
https://www.cambridge.org/core/product/identifier/S0962492914000130/type/journal_article

  * [56] N. S. Gokhale, Practical finite element analysis, Finite to infinite, 2008.
  * [57] O. Davydov, D. T. Oanh, N. M. Tuong, Improved stencil selection for meshless finite difference methods in 3D, Journal of Computational and Applied Mathematics 425 (2023) 115031.  doi:10.1016/j.cam.2022.115031.
URL https://linkinghub.elsevier.com/retrieve/pii/S037704272200629X

  * [58] S. Zhu, L. Zhang, L. Wu, L. Tan, H. Chen, Application of improved picard iteration method to simulate unsaturated flow and deformation in deformable porous media, Computers and Geotechnics 166 (2024) 105979\.  doi:https://doi.org/10.1016/j.compgeo.2023.105979.
URL https://www.sciencedirect.com/science/article/pii/S0266352X2300736X