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.

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 [59].

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 [1215] 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 [1721].

Furthermore, meshless or meshfree 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 called clouds. The clouds consist of one central point (the star 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 [2225].

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. Section 2 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 [3439], 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 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 , where 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 function , which is approximated by ; is only valid in the cloud associated with the star node . The function is a linear combination of known functions , where is the vector that represents the basis of linearly independent functions and is a vector of constant parameters that are only valid in . The elements of the interpolation basis may belong to any function family. Nevertheless, for simplicity, the first monomials polynomials are used. Examples of this basis in 2D are for and for .

Because (2) is valid for all points of the th subdomain, the approximations conform to a Vandermonde system that is given by the relation where

In general, the number of points that conform to the cloud is greater than the number of functions that define the basis; hence, the matrix is usually rectangular. Thus, the interpolation property is lost, and the problem must be addressed with a numerical approximation. The coefficients of the vector must be determined so that the weighted sums of the squared differences between the exact values and the approximated values of each point are minimized according to where is a fixed weighting function that is defined in and evaluated for the node ; see [26, 27]. The minimization process described by (6) leads to the expression for the vector : where represents the unknown parameters that are sought on the cloud. is defined as Additionally, the matrices , , and are given as and is an diagonal matrix defined by where the weighting functions are derived to have unit values near the star node and zero values outside the subdomains. Under the FPM, the common selection is the normalized Gaussian, where is the distance between the star node and the point , (max. of ) is a reference distance, and . A detailed description of the effects of the constant parameters and on the numerical approximation and the guidelines to set their values are presented in [45]. Other considerations in selecting the function can be found in [26, 27, 46]. Finally, replacing (7) in (2) leads to where is the shape function matrix, which is defined as where . Note that according to the least-square nature of the approximation, ; see Figure 1. In other words, the local values of the approximating function do not fit the nodal unknown values. Indeed, is the true approximation, which we use to satisfy the differential equation and the boundary conditions. In this context, are simply the unknown parameters that we aim to determine. According to these concepts and (9), it is possible to obtain where and denote the first and second space derivatives, respectively. Note that these derivatives are computed by taking the derivative of the basis functions in (2).

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, where and are differential operators and is the unknown function of the problem. In solid mechanics problems, and correspond to the equilibrium equations, is the displacement field, denotes the body forces that act on the domain , denotes the external tractions on , and represents the prescribed values of the unknown function on .

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 where is the number of points in and and are the numbers of points on the boundaries and , 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 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 [2628, 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 where denotes the components of the unit normal to the boundary and denotes the stabilization-length parameters that are obtained from see Figure 2.

The efficiency of the FIC stabilization procedure in the context of FPM and solid mechanics has been presented in [3032, 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 variable is defined from the internal variable as the weighted average over a domain in the vicinity of point [49]: where the weight function depends only on the distance ; hence, . A possible choice for is the Gaussian weight function, as proposed in [50].

The weight function introduces the internal length scale , 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 , but because of the shape of the weight function the internal length scale defines the region of the body that surrounds point , which significantly influences the behavior at that point [51].

In this work, the local accumulated plastic strain is considered the evolutionary internal variable. Then, the explicit nonlocal gradient formulation of this plastic strain becomes [11]

To obtain the implicit gradient formulation of the accumulated plastic strain , (19) is rewritten as 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, 5254]).

The implicit formulation of (20) introduces the gradient effects and is accompanied by the boundary condition where 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 where is the stress tensor and the threshold stress in softening is given by and parameter is the linear softening modulus.

More details about strain gradient plasticity and some thermodynamics aspects can be found in [14, 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 where and are the total and the plastic strain tensors, respectively, and is the elastic moduli tensor.

The rate-type stress-strain relationship obtained from (24) is given by which can be rewritten as where is the plastic multiplier, which is given by and is the gradient to the yield surface, which is defined in (22).

However, the yield function (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:

During the plastic flow, the stress point must remain on the yield surface according to the equation which is known as the Prager’s consistency condition.

Considering the above dependencies, (29) can be rewritten as which corresponds to the consistency condition equation under plastic deformation gradients.

The plastic consistency equation for a time step can be approximated using a Taylor series to yield where the increments correspond to

If after a time step , the plastic flow condition has been reached at any point; that is, ; (31) can be reworked to give 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)): which reduces the computational cost.

More details about the numerical implementation issues can be found in [2, 10, 5759]. 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.

For an internal length scale of , 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.

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.

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.

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 (Figure 10). 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.

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.

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 Figure 11. 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.

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.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

Luis Pérez Pozo acknowledges the financial support from the Chilean Agency CONICYT (FONDECYT project 1140583).