# A Finite Points Method Approach for Strain Localization Using the Gradient Plasticity Formulation

Authors: Luis Pérez Pozo; Andy Campos; Sheila Lascano; Sergio Oller; Antonio Rodríguez-Ferran
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The softening elastoplastic models present an unsuitable behavior after reaching
the yield strength: unbounded strain localization. Because of the material
instability, which is reflected in the loss of ellipticity of the governing partial
differential equations, the solution depends on the discretization.
The present work proposes to solve this dependency using the meshless Finite
Points Method. This meshfree spatial discretization technique allows enriching
the governing equations using gradient’s plasticity and introducing an internal
length scale parameter at the material model in order to objectify the solution.


---

## Body

## 1. Introduction

The strain localization phenomenon is usually a precursor of material failure. In geomaterials such as soils, concrete, and rocks, strain localization has been observed in many situations such as triaxial tests in the lab and earth excavation processes in the field [1].If the classical continuum theory is directly applied to elastoplastic softening materials, a loss of ellipticity of the governing partial differential equations occurs, and the boundary value cannot correctly describe the physical problem [2]. The strains are localized in zones that have a width related to the characteristic length of the subdomains that are used in the discretization. When these subdomains become smaller, the width of the strain band also decreases without a lower bound. If the discretization length tends to zero, there is no energy dissipation, which does not make physical sense [3]. To overcome the loss of ellipticity and the related problems, an implicit or explicit length scale must be incorporated into the material description or the formulation of the boundary value problem [4]. These kinds of regularization techniques, called localization limiters, are usually based on various forms of enriched continuum theories; see [5–9].Several proposed strategies consider an internal material length scale. One strategy introduces nonlocal theories and gradient methods into the model [10]. These concepts can be related by a Taylor series expansion, as shown by Peerlings et al. [11].The nonlocal constitutive theories hold that the local state of the material at a given point may not be sufficient to evaluate the stress at that point. This statement can be physically justified by the fact that no real material is an ideal continuous medium; on a sufficiently small scale, the effects of heterogeneity and discontinuous microstructure become nonnegligible [7].The gradient plasticity formulation [12–15] introduces higher-order gradient terms of a nonlocal variable into the differential equation that governs its evolution. The origins of these gradient terms are discussed in Groma and Vörös [16] for different cases that are related to atomic distance, dislocations, and finite grain size. Here, the internal material length scale is implicitly introduced through the regularization parameter. The choice of the variable to be represented nonlocally depends on the described material behavior and the chosen approach [17–21].Furthermore,meshless ormeshfree methods possess intrinsic nonlocal properties. The approximation functions are not locally constructed because the support size of the weight function is greater than the nodal spacing; therefore, the approximation is inherently nonlocal. These nonlocal properties are used to incorporate an intrinsic length scale that regularizes the problems of material instabilities [4].Meshless methods are a family of numerical techniques that do not require a mesh. In these methods, the body or the domain is discretized using a collection of points. It is divided into the local interpolation subdomains, which are also calledclouds. The clouds consist of one central point (thestar node) and several neighboring points. Generally, these methods are computationally efficient and easy to implement, and they have been successfully used in several applications. The general characteristics, classifications, advantages, and disadvantages of these methods have been described in [22–25].This work is focused on a meshless Finite Points Method (FPM) for the strain localization analysis through a gradient plasticity formulation. This method uses an implicit nonlocal approach to represent the global material behavior. The present contribution solves the solution-discretization dependency using the meshless FPM. This meshfree spatial discretization technique allows enriching the governing equations using the gradients’ plasticity and introducing an internal length scale parameter at the material model. The same shape functions that conform to the meshless approximation are used in the gradient additional terms, which reduces the computational cost.The paper is organized as follows. Section2 introduces the basic formulation of the meshless FPM and provides a brief overview of its main features. Section 3 focuses on the basic formulation of strain gradient plasticity. Section 4 introduces some aspects related to the numerical implementation. In Section 5, we give several numerical examples. Conclusions are drawn in Section 6.


## 2. The Finite Point Method (FPM)

In this section, we review the basic formulation of the meshless FPM and provide a brief overview of its main features.This method was proposed by Oñate et al. [26, 27] to solve convective-transport and fluid-flow problems. Its application has been extended to advection-diffusion transport [28], incompressible-flow problems [29], elasticity [30, 31], solid dynamics [32], solidification modeling [33], nonlinear material behavior problems [34–39], adaptive refinement [40, 41], and the large deflection analysis of flexible plates [42]. The lack of dependence on a mesh or an integration procedure is an important feature that makes the FPM a truly meshless method.To obtain the final system of discrete equations, the FPM approximates the local solution of a partial differential equation at each point of the discretized domain using a weighted-least-square technique and a point-collocation procedure. Because the approximation procedure used by this method is local, it is necessary to define a subdomainΩk for each node that contains the neighboring nodes, which are selected using a suitable criterion [43, 44]. This collection of points is called a cloud, and its referential central point is the star node. For example, a relevant aspect in the definition of clouds is that their superposition must produce the entire domain Ω,
(1)⋃k=1NpΩk=Ω,
where Np is the total number of nodes. Note that the definition of clouds is the basic, initial step in implementing the FPM approximation using fixed weighted least squares.


### 2.1. Fixed Weighted-Least-Square Approximation

With the defined discretized domain, let us define a functionu(x), which is approximated by u^(x); u^(x) is only valid in the cloud Ωk associated with the star node xk. The function u(x) is a linear combination of known functions p(x),
(2)ux≅u^x=pTxαk,∀x∈Ωk,
where p(x) is the vector that represents the basis of m linearly independent functions and αk is a vector of constant parameters that are only valid in Ωk. The elements of the interpolation basis may belong to any function family. Nevertheless, for simplicity, the first m monomials polynomials are used. Examples of this basis in 2D are pT(x)=[1,x,y] for m=3 and pT(x)=[1,x,y,x2,xy,y2] for m=6.Because (2) is valid for all Nc points of the kth subdomain, the approximations u^(X) conform to a Vandermonde system that is given by the relation
(3)u(Xk)≅u^(Xk)=P(Xk)·αk,
where
(4)Xk=xk,1⋯xk,NcT,u(Xk)=uxk,1⋯uxk,NcT,u^(Xk)=u^xk,1⋯u^xk,NcT,αk=αk,1⋯αk,NcT,P(Xk)=p(xk,1)⋮p(xk,Nc).In general, the number of pointsNc that conform to the cloud is greater than the number of functions m that define the basis; hence, the matrix P(Xk) is usually rectangular. Thus, the interpolation property is lost, and the problem must be addressed with a numerical approximation. The coefficients of the vector αk must be determined so that the weighted sums of the squared differences between the exact values u(x) and the approximated values u^(x) of each point are minimized according to
(5)min⁡∑j=1Nc‍w(xj)·u(xj)-u^(xj)2,
where w(xj) is a fixed weighting function that is defined in Ωk and evaluated for the node xj; see [26, 27]. The minimization process described by (6) leads to the expression for the vector αk:
(6)αk=A-1(Xk)B(Xk)λ(Xk),
where λ(Xk) represents the unknown parameters that are sought on the cloud. Ωk is defined as
(7)λ(Xk)=λxk,1⋯λxk,NcT.
Additionally, the matrices A(Xk), B(Xk), and W(Xk) are given as
(8)A(Xk)=P(Xk)W(Xk)PT(Xk),B(Xk)=PT(Xk)W(Xk),
and W(Xk) is an Nc×Nc diagonal matrix defined by
(9)W(Xk)=diag⁡wxk,1⋯w(xk,Nc),
where the weighting functions w(xj) are derived to have unit values near the star node and zero values outside the Ωk subdomains. Under the FPM, the common selection is the normalized Gaussian,
(10)w(xj)=exp⁡-hj/ζ-exp⁡-r/ζ1-exp⁡-r/ζ,ifhj≤r,0,ifhj>r,
where hj is the distance between the star node and the point xj, r=q·hmax⁡ (max. of hj) is a reference distance, and ζ=β·r. A detailed description of the effects of the constant parameters q and β on the numerical approximation and the guidelines to set their values are presented in [45]. Other considerations in selecting the function w(xj) can be found in [26, 27, 46]. Finally, replacing (7) in (2) leads to
(11)u^(x)=NT(x)λ(Xk),
where N(x) is theshape function matrix, which is defined as
(12)N(x)=pT(x)C(Xk),
where C(Xk)=A-1(Xk)B(Xk). Note that according to the least-square nature of the approximation, u(x)≅u^(x)≠λ(x); see Figure 1. In other words, the local values of the approximating function do not fit the nodal unknown values. Indeed, u^(x) is the true approximation, which we use to satisfy the differential equation and the boundary conditions. In this context, λ(x) are simply the unknown parameters that we aim to determine. According to these concepts and (9), it is possible to obtain
(13)u^xx=NxTxλXk,u^xx(x)=NxxT(x)λ(Xk),
where (·)x and (·)xx denote the first and second space derivatives, respectively. Note that these derivatives are computed by taking the derivative of the basis functions p(x) in (2).Figure 1
Weighted least squares procedure.


### 2.2. Discrete-Equation System and the Point Collocation Scheme

Let us assume a boundary value problem governed by the set of differential equations with boundary conditions,(14)Au=binΩ,Bu=ton⁡Γt,u-up=0on⁡Γu,
where A and B are differential operators and u is the unknown function of the problem. In solid mechanics problems, A and B correspond to the equilibrium equations, u is the displacement field, b denotes the body forces that act on the domain Ω, t denotes the external tractions on Γt, and up represents the prescribed values of the unknown function on Γu.Using the point collocation scheme [47], the differential equations (14), and the final FPM approximation, which is defined in (11), one obtains the system of discrete equations
(15)(u^)i-bi=0inΩ,i=1,Np,B(u^)j-tj=0on⁡Γt,j=1,Nt,u^-upk=0on⁡Γu,k=1,Nu,
where Np is the number of points in Ω and Nt and Nu are the numbers of points on the boundaries Γt and Γu, respectively.


### 2.3. Stabilization Procedure

The collocation procedure usually leads to an ill-conditioned system of equations with unstable and inaccurate results mainly because the point collocation method cannot satisfy the equilibrium over a cloud of finite-sized points by simply sampling the equation at the star node in the cloud. These deficiencies are more pronounced in clouds near the boundaryΓt because the clouds in these regions are not symmetric [30].A stabilized form of the governing equation is used in this work. This formulation is derived from the finite calculus (FIC) procedure that was described in [28, 48]. The FIC method is based on imposing the typical balance laws of mechanics over a finite-sized domain. Then, the unknown fields are approximated within the finite domain using a Taylor series expansion and retaining the higher-order terms over those used in the standard infinitesimal approach [26–28, 30]. This method naturally introduces new terms in the governing differential equation, and these terms have stabilization features. The stabilized form of (15), which uses the FIC method, is
(16)A(u^)-12hk∂∂xkA(u^)k-bk=0inΩ,B(u^)-12hknk∂∂xkA(u^)k-tk=0on⁡Γt,u^k-up=0on⁡Γu,
where nk denotes the components of the unit normal to the boundary Γt and hk denotes the stabilization-length parameters that are obtained from
(17)hT=dxmin⁡dymin⁡;
see Figure 2.Definition of stabilization length for (a) the interior and (b) the boundary nodes.
(a)
(b)The efficiency of the FIC stabilization procedure in the context of FPM and solid mechanics has been presented in [30–32, 35, 41].


## 2.1. Fixed Weighted-Least-Square Approximation

With the defined discretized domain, let us define a functionu(x), which is approximated by u^(x); u^(x) is only valid in the cloud Ωk associated with the star node xk. The function u(x) is a linear combination of known functions p(x),
(2)ux≅u^x=pTxαk,∀x∈Ωk,
where p(x) is the vector that represents the basis of m linearly independent functions and αk is a vector of constant parameters that are only valid in Ωk. The elements of the interpolation basis may belong to any function family. Nevertheless, for simplicity, the first m monomials polynomials are used. Examples of this basis in 2D are pT(x)=[1,x,y] for m=3 and pT(x)=[1,x,y,x2,xy,y2] for m=6.Because (2) is valid for all Nc points of the kth subdomain, the approximations u^(X) conform to a Vandermonde system that is given by the relation
(3)u(Xk)≅u^(Xk)=P(Xk)·αk,
where
(4)Xk=xk,1⋯xk,NcT,u(Xk)=uxk,1⋯uxk,NcT,u^(Xk)=u^xk,1⋯u^xk,NcT,αk=αk,1⋯αk,NcT,P(Xk)=p(xk,1)⋮p(xk,Nc).In general, the number of pointsNc that conform to the cloud is greater than the number of functions m that define the basis; hence, the matrix P(Xk) is usually rectangular. Thus, the interpolation property is lost, and the problem must be addressed with a numerical approximation. The coefficients of the vector αk must be determined so that the weighted sums of the squared differences between the exact values u(x) and the approximated values u^(x) of each point are minimized according to
(5)min⁡∑j=1Nc‍w(xj)·u(xj)-u^(xj)2,
where w(xj) is a fixed weighting function that is defined in Ωk and evaluated for the node xj; see [26, 27]. The minimization process described by (6) leads to the expression for the vector αk:
(6)αk=A-1(Xk)B(Xk)λ(Xk),
where λ(Xk) represents the unknown parameters that are sought on the cloud. Ωk is defined as
(7)λ(Xk)=λxk,1⋯λxk,NcT.
Additionally, the matrices A(Xk), B(Xk), and W(Xk) are given as
(8)A(Xk)=P(Xk)W(Xk)PT(Xk),B(Xk)=PT(Xk)W(Xk),
and W(Xk) is an Nc×Nc diagonal matrix defined by
(9)W(Xk)=diag⁡wxk,1⋯w(xk,Nc),
where the weighting functions w(xj) are derived to have unit values near the star node and zero values outside the Ωk subdomains. Under the FPM, the common selection is the normalized Gaussian,
(10)w(xj)=exp⁡-hj/ζ-exp⁡-r/ζ1-exp⁡-r/ζ,ifhj≤r,0,ifhj>r,
where hj is the distance between the star node and the point xj, r=q·hmax⁡ (max. of hj) is a reference distance, and ζ=β·r. A detailed description of the effects of the constant parameters q and β on the numerical approximation and the guidelines to set their values are presented in [45]. Other considerations in selecting the function w(xj) can be found in [26, 27, 46]. Finally, replacing (7) in (2) leads to
(11)u^(x)=NT(x)λ(Xk),
where N(x) is theshape function matrix, which is defined as
(12)N(x)=pT(x)C(Xk),
where C(Xk)=A-1(Xk)B(Xk). Note that according to the least-square nature of the approximation, u(x)≅u^(x)≠λ(x); see Figure 1. In other words, the local values of the approximating function do not fit the nodal unknown values. Indeed, u^(x) is the true approximation, which we use to satisfy the differential equation and the boundary conditions. In this context, λ(x) are simply the unknown parameters that we aim to determine. According to these concepts and (9), it is possible to obtain
(13)u^xx=NxTxλXk,u^xx(x)=NxxT(x)λ(Xk),
where (·)x and (·)xx denote the first and second space derivatives, respectively. Note that these derivatives are computed by taking the derivative of the basis functions p(x) in (2).Figure 1
Weighted least squares procedure.


## 2.2. Discrete-Equation System and the Point Collocation Scheme

Let us assume a boundary value problem governed by the set of differential equations with boundary conditions,(14)Au=binΩ,Bu=ton⁡Γt,u-up=0on⁡Γu,
where A and B are differential operators and u is the unknown function of the problem. In solid mechanics problems, A and B correspond to the equilibrium equations, u is the displacement field, b denotes the body forces that act on the domain Ω, t denotes the external tractions on Γt, and up represents the prescribed values of the unknown function on Γu.Using the point collocation scheme [47], the differential equations (14), and the final FPM approximation, which is defined in (11), one obtains the system of discrete equations
(15)(u^)i-bi=0inΩ,i=1,Np,B(u^)j-tj=0on⁡Γt,j=1,Nt,u^-upk=0on⁡Γu,k=1,Nu,
where Np is the number of points in Ω and Nt and Nu are the numbers of points on the boundaries Γt and Γu, respectively.


## 2.3. Stabilization Procedure

The collocation procedure usually leads to an ill-conditioned system of equations with unstable and inaccurate results mainly because the point collocation method cannot satisfy the equilibrium over a cloud of finite-sized points by simply sampling the equation at the star node in the cloud. These deficiencies are more pronounced in clouds near the boundaryΓt because the clouds in these regions are not symmetric [30].A stabilized form of the governing equation is used in this work. This formulation is derived from the finite calculus (FIC) procedure that was described in [28, 48]. The FIC method is based on imposing the typical balance laws of mechanics over a finite-sized domain. Then, the unknown fields are approximated within the finite domain using a Taylor series expansion and retaining the higher-order terms over those used in the standard infinitesimal approach [26–28, 30]. This method naturally introduces new terms in the governing differential equation, and these terms have stabilization features. The stabilized form of (15), which uses the FIC method, is
(16)A(u^)-12hk∂∂xkA(u^)k-bk=0inΩ,B(u^)-12hknk∂∂xkA(u^)k-tk=0on⁡Γt,u^k-up=0on⁡Γu,
where nk denotes the components of the unit normal to the boundary Γt and hk denotes the stabilization-length parameters that are obtained from
(17)hT=dxmin⁡dymin⁡;
see Figure 2.Definition of stabilization length for (a) the interior and (b) the boundary nodes.
(a)
(b)The efficiency of the FIC stabilization procedure in the context of FPM and solid mechanics has been presented in [30–32, 35, 41].


## 3. Strain Gradient Plasticity Formulation

In this section, we present the basic formulation of strain gradient plasticity. We start from the definition of the nonlocal approach and obtain the implicit gradient form of the accumulated plastic strain.The nonlocal variableg¯(x) is defined from the internal variable g(x) as the weighted average over a domain V in the vicinity of point x [49]:
(18)g¯(x)=∫V‍ϑ(x,ξ)g(ξ)dξ,
where the weight function ϑ depends only on the distance |x-ξ|; hence, ϑ≡ϑ(x,ξ). A possible choice for ϑ is the Gaussian weight function, as proposed in [50].The weight functionϑ introduces the internal length scale c, which is considered a material parameter with the dimension of length. This parameter defines the dimension of the neighborhood that affects the nonlocal function. Note that the average is extended to the whole domain V, but because of the shape of the weight function the internal length scale c defines the region of the body that surrounds point x, which significantly influences the behavior at that point [51].In this work, the local accumulated plastic strainκ(x) is considered the evolutionary internal variable. Then, the explicit nonlocal gradient formulation κ¯(x) of this plastic strain becomes [11]
(19)κ¯(x)=κ(x)-c2∇2κ(x).To obtain the implicit gradient formulation of the accumulated plastic strainκ, (19) is rewritten as
(20)κ(x)=κ¯(x)+c2∇2κ¯(x).
Several advantages of using an implicit formulation instead of an explicit formulation have been reported in the literature. For example, the continuity requirements for the shape functions in numerical simulations are lower, the dynamical behavior is more realistic, a better approximation of the underlying nonlocal integral model is provided, and the singularities in the strain field are more adequately damped (see [11, 52–54]).The implicit formulation of (20) introduces the gradient effects and is accompanied by the boundary condition
(21)∇κ¯(x)·n=0,
where n is the outward unit normal to the boundary. In the literature, there is no agreement respect to the physical sense of this boundary condition. Nevertheless, some thermodynamical aspects related to the nonlocal residual energy have been interpreted by Polizzotto [55].In the nonlocal gradient plasticity formulation, the scalar yield function is defined by(22)f(σ,κ¯,∇2κ¯)=σ(κ¯)-σth(κ¯,∇2κ¯),
where σ is the stress tensor and the threshold stress in softening is given by
(23)σth(κ¯,∇2κ¯)=σ0th+h·(κ¯-c2∇2κ¯),
and parameter h is the linear softening modulus.More details about strain gradient plasticity and some thermodynamics aspects can be found in [1–4, 21, 56].


## 4. Numerical Implementation

Some aspects related to the numerical implementation of this proposal are shown in this section.The elastoplastic constitutive model is expressed by the stress-strain relationship(24)σ=D:(ε-εp),
where ε and εp are the total and the plastic strain tensors, respectively, and D is the elastic moduli tensor.The rate-type stress-strain relationship obtained from (24) is given by
(25)σ˙=D:(ε˙-ε˙p),
which can be rewritten as
(26)σ˙=D:(ε˙-γ˙m),
where γ˙ is the plastic multiplier, which is given by
(27)γ˙=κ¯˙,
and m is the gradient to the yield surface, which is defined in (22).However, the yield functionf(σ,κ¯,∇2κ¯) (Equation (22)) and the plastic multiplier γ˙ (Equation (27)) evolve according to the loading-unloading or the Kuhn-Tucker conditions and the consequent consistency (or persistency) condition, which establishes that the stress state should persist on the yield surface for the plastic flow to occur:
(28)f(σ,κ¯,∇2κ¯)⩽0,γ˙⩾0,f(σ,κ¯,∇2κ¯)·γ˙=0.During the plastic flow, the stress point must remain on the yield surface according to the equation(29)f˙(σ,κ¯,∇2κ¯)=ddt(σ,κ¯,∇2κ¯)=0,
which is known as the Prager’s consistency condition.Considering the above dependencies, (29) can be rewritten as
(30)∂f∂σ:σ˙+∂f∂κ¯κ¯˙+∂f∂∇2κ¯∇2κ¯˙=0,
which corresponds to the consistency condition equation under plastic deformation gradients.The plastic consistency equation for a time stepΔt can be approximated using a Taylor series to yield(31)dfdtσ,κ¯,∇2κ¯t+Δt≈fσ,κ¯,∇2κ¯t+∂f∂σt:Δσ+∂f∂κ¯tΔκ¯+∂f∂∇2κ¯tΔ∇2κ¯,
where the increments correspond to(32a)Δσ=gt+Δtσ-vtσ,(32b)Δκ=κt+Δt-κt,(32c)Δ∇2κ¯=∇2t+Δtκ¯-∇2tκ¯.If after a time stepΔt, the plastic flow condition has been reached at any point; that is, (df/dt)(σ,κ¯,∇2κ¯)|t+Δt=0; (31) can be reworked to give
(33)∂f∂σt:Δσ+∂f∂κ¯tΔκ¯+∂f∂∇2κ¯tΔ∇2κ¯=-fσ,κ¯,∇2κ¯t.
It should be noted that, in this work, all gradient terms are obtained using the same shape functions that conform to the FPM meshless approximation (see (13)):
(34)∇2κ¯≈NxxT(x)κ¯(Xk),
which reduces the computational cost.More details about the numerical implementation issues can be found in [2, 10, 57–59]. On the other hand, numerical implementation aspect in the meshfree context can be found in [1, 4, 60].


## 5. Numerical Examples

### 5.1. Uniaxial Tensile Test

This numerical problem is based on the reference solution that was presented by de Borst and Muehlhaus [2]. The central tenth of the bar is weakened (10% reduction in the Young’s modulus) to induce localization, as shown in Figure 3. The numerical tests are displacement controlled. The nondimensional geometrical and material parameters are summarized in Table 1.Table 1
Uniaxial tensile test. Geometrical and material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
20000
Initial yield stress
σ
y
2
Softening modulus
h
-
2000
Length of bar
L
0
100
Length of weakened part
l
w
10
Internal length scale
c
5
Prescribed displacement
u
p
0.02Figure 3
Uniaxial tensile test statement.For an internal length scale ofc/L0=0.05, different equidistant discretizations have been studied (see Figure 4). As shown in that figure, dependence on the discretization becomes negligible for nodal densities of 161 and beyond.Uniaxial tensile test: stress-displacement curves for several discretizations.
(a)
(b)


### 5.2. Biaxial Compression Test

Next, the classical biaxial compression test [61] is analyzed using the nonlocal plastic model with the proposed linear softening. The weakened imperfection at the left bottom corner of the geometry is used to illustrate the discretization insensitivity, as shown in Figure 5. The numerical results are compared with those shown by Rodríguez-Ferran et al. [58]. The initial yield stress of the weakened area is reduced by 10%, and the study is performed for three different discretizations. The nondimensional material parameters are summarized in Table 2.Table 2
Biaxial compression test: material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
11920
Poisson’s ratio
ν
0.3
Initial yield stress
σ
y
100
Softening modulus
h
-
400
Internal length scale
c
0.5
Prescribed displacement
v
p
-
1.5Figure 5
Biaxial compression test: problem statement.The results are clearly discretization-insensitive and consistent with Rodríguez-Ferran et al. [58], as shown in the force-displacement curves in Figure 6 and the deformation patterns and the equivalent plastic strain in Figure 7.Figure 6
Biaxial compression test: results for different discretizations.Biaxial compression test: deformation pattern and equivalent plastic strain for (a)37×19 nodes, (b) 61×31 nodes, and (c) 85×54 nodes.
(a)
(b)
(c)


### 5.3. Sensitivity Test for Biaxial Compression

This example is a variation of the biaxial compression test, which was presented in the previous section. The goal is to illustrate the imperfection-size insensitivity (Figure10). The initial yield stress of the weakened areas was reduced by 10%, and the study was performed for two imperfection sizes, as shown in Figure 8. The numerical results are compared with those shown by Rodríguez-Ferran et al. [58]. The nondimensional material parameters are summarized in Table 3.Table 3
Sensitivity test for biaxial compression: material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
11920
Poisson’s ratio
ν
0.3
Initial yield stress
σ
y
100
Softening modulus
h
-
400
Internal length scale
c
0.5
Prescribed displacement
v
p
-
2.4Sensitivity test for biaxial compression: (a) a small and (b) a large imperfection.
(a)
(b)The results are clearly insensitive to the size of the imperfection and consistent with Rodríguez-Ferran et al. [58], as shown in the force-displacement curves in Figure 6 in addition to the deformation patterns and the equivalent plastic strain in Figure 9.Sensitivity test for biaxial compression: the deformation patterns and the equivalent plastic strain for (a) small and (b) large imperfections.
(a)
(b)Figure 10
Sensitivity test for biaxial compression: normalized force-displacement curves for the two imperfection sizes.


### 5.4. Diagonally Weakened Biaxial Tensile Test

In this numerical example, a localization strain band is induced in a weakened area on the diagonal, as shown in Figure11. The numerical results are compared with those shown by Liebe et al. [62]. The initial yield stress of the weakened area was reduced by 20%, and the study is performed for two different discretizations. The nondimensional material parameters are summarized in Table 4.Table 4
Diagonally weakened biaxial tensile test: material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
50000
Poisson’s ratio
ν
0.3
Initial yield stress
σ
y
25
Softening modulus
h
-
15
Internal length scale
c
0.5
Prescribed displacement
v
p
0.2Figure 11
Diagonally weakened biaxial tensile test: problem statement.The results are discretization-insensitive and consistent with Liebe et al. [62], as shown in the reaction-displacement curves in Figure 12, in addition to the deformation patterns and the equivalent plastic strain in Figure 13.Figure 12
Diagonally weakened biaxial tensile test: reaction-displacement curves for two discretizations.Figure 13
Diagonally weakened biaxial tensile test: examples of the deformation pattern and the equivalent plastic strain for a discretization of33×33 nodes.


## 5.1. Uniaxial Tensile Test

This numerical problem is based on the reference solution that was presented by de Borst and Muehlhaus [2]. The central tenth of the bar is weakened (10% reduction in the Young’s modulus) to induce localization, as shown in Figure 3. The numerical tests are displacement controlled. The nondimensional geometrical and material parameters are summarized in Table 1.Table 1
Uniaxial tensile test. Geometrical and material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
20000
Initial yield stress
σ
y
2
Softening modulus
h
-
2000
Length of bar
L
0
100
Length of weakened part
l
w
10
Internal length scale
c
5
Prescribed displacement
u
p
0.02Figure 3
Uniaxial tensile test statement.For an internal length scale ofc/L0=0.05, different equidistant discretizations have been studied (see Figure 4). As shown in that figure, dependence on the discretization becomes negligible for nodal densities of 161 and beyond.Uniaxial tensile test: stress-displacement curves for several discretizations.
(a)
(b)


## 5.2. Biaxial Compression Test

Next, the classical biaxial compression test [61] is analyzed using the nonlocal plastic model with the proposed linear softening. The weakened imperfection at the left bottom corner of the geometry is used to illustrate the discretization insensitivity, as shown in Figure 5. The numerical results are compared with those shown by Rodríguez-Ferran et al. [58]. The initial yield stress of the weakened area is reduced by 10%, and the study is performed for three different discretizations. The nondimensional material parameters are summarized in Table 2.Table 2
Biaxial compression test: material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
11920
Poisson’s ratio
ν
0.3
Initial yield stress
σ
y
100
Softening modulus
h
-
400
Internal length scale
c
0.5
Prescribed displacement
v
p
-
1.5Figure 5
Biaxial compression test: problem statement.The results are clearly discretization-insensitive and consistent with Rodríguez-Ferran et al. [58], as shown in the force-displacement curves in Figure 6 and the deformation patterns and the equivalent plastic strain in Figure 7.Figure 6
Biaxial compression test: results for different discretizations.Biaxial compression test: deformation pattern and equivalent plastic strain for (a)37×19 nodes, (b) 61×31 nodes, and (c) 85×54 nodes.
(a)
(b)
(c)


## 5.3. Sensitivity Test for Biaxial Compression

This example is a variation of the biaxial compression test, which was presented in the previous section. The goal is to illustrate the imperfection-size insensitivity (Figure10). The initial yield stress of the weakened areas was reduced by 10%, and the study was performed for two imperfection sizes, as shown in Figure 8. The numerical results are compared with those shown by Rodríguez-Ferran et al. [58]. The nondimensional material parameters are summarized in Table 3.Table 3
Sensitivity test for biaxial compression: material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
11920
Poisson’s ratio
ν
0.3
Initial yield stress
σ
y
100
Softening modulus
h
-
400
Internal length scale
c
0.5
Prescribed displacement
v
p
-
2.4Sensitivity test for biaxial compression: (a) a small and (b) a large imperfection.
(a)
(b)The results are clearly insensitive to the size of the imperfection and consistent with Rodríguez-Ferran et al. [58], as shown in the force-displacement curves in Figure 6 in addition to the deformation patterns and the equivalent plastic strain in Figure 9.Sensitivity test for biaxial compression: the deformation patterns and the equivalent plastic strain for (a) small and (b) large imperfections.
(a)
(b)Figure 10
Sensitivity test for biaxial compression: normalized force-displacement curves for the two imperfection sizes.


## 5.4. Diagonally Weakened Biaxial Tensile Test

In this numerical example, a localization strain band is induced in a weakened area on the diagonal, as shown in Figure11. The numerical results are compared with those shown by Liebe et al. [62]. The initial yield stress of the weakened area was reduced by 20%, and the study is performed for two different discretizations. The nondimensional material parameters are summarized in Table 4.Table 4
Diagonally weakened biaxial tensile test: material parameters (nondimensional).
Meaning
Symbol
Value
Elasticity modulus
E
50000
Poisson’s ratio
ν
0.3
Initial yield stress
σ
y
25
Softening modulus
h
-
15
Internal length scale
c
0.5
Prescribed displacement
v
p
0.2Figure 11
Diagonally weakened biaxial tensile test: problem statement.The results are discretization-insensitive and consistent with Liebe et al. [62], as shown in the reaction-displacement curves in Figure 12, in addition to the deformation patterns and the equivalent plastic strain in Figure 13.Figure 12
Diagonally weakened biaxial tensile test: reaction-displacement curves for two discretizations.Figure 13
Diagonally weakened biaxial tensile test: examples of the deformation pattern and the equivalent plastic strain for a discretization of33×33 nodes.


## 6. Conclusions

A meshless Finite Point Method for the strain localization analysis using a gradient plasticity formulation has been presented.This meshfree spatial discretization technique allows enriching the governing equations using the gradients’ plasticity and introducing an internal length scale parameter at the material model.The differentiability order of the polynomial base of approximation allows us to use the same shape function to approximate the displacements field and the internal plastic deformation variable, which reduces the computational cost.The numerical results, which were obtained using the proposed technique, demonstrate that the solution and the response are independent of the discretization density.



---

*Source: 782079-2014-12-11.xml*