Numerical solutions of strain localization with nonlocal softening plasticity

Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Numerical solutions of strain localization with nonlocal softening plasticity Xilin Lu a,b,c, Jean-Pierre Bardet b,*, Maosong Huang a,c a

Department of Geotechnical Engineering, Tongji University, Siping 1239, Shanghai 200092, China Sonny Astani Department of Civil and Environmental Engineering, University of Southern California, Los Angeles, CA 90089-2531, USA c Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai 200092, China b

a r t i c l e

i n f o

Article history: Received 13 September 2008 Received in revised form 22 July 2009 Accepted 1 August 2009 Available online 8 August 2009 Keywords: Nonlocal Softening plasticity Characteristic length Strain localization Spectral analysis

a b s t r a c t Softening plasticity, when it is deprived of a characteristic length, often leads to boundary value problems that are ill-posed and produce unreliable numerical solutions. As meshes are refined, plastic strains localize into bands that become narrower and narrower, and generate load–displacement responses with excessive softening, which may even degenerate into unrealistic snapbacks. In the case of discrete systems, the underlying reasons for this poor numerical performance can be examined using a spectral analysis of the tangential stiffness matrix derived from incremental equilibrium equations. In one-dimensional boundary-valued problems, the introduction of a weak softening element into a mesh of elastoplastic elements produces a dominant eigenvector that clearly exhibits a strain localization width directly related to the weak element size. The dominant eigenvector that controls the incremental solution varies spuriously when the mesh is refined, which results in mesh dependency. Over-nonlocal softening plasticity introduces a length scale, forces more elements to become plastic, and preserves the width of the plastic zone when the mesh is refined. When one-dimensional meshes are refined, spectral analysis shows that the dominant eigenvector does not vary erratically, and that its deformation patterns are smooth and display a constant localization width. This explains why the spatial distribution of plastic strains and the global load–displacement responses converge to analytical solutions in one dimension. The generalization to higher dimensions will be the object of future studies. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Numerous studies have shown that numerical solutions suffer from mesh dependency in case of softening plasticity [1], mainly because the governing equations become hyperbolic in some unknown area and make the boundary value problem (BVP) ill-posed [2]. Ill-posed problems can be regularized using various strategies, including visco-plasticity [1,3], Cosserat or micro-polar models [4– 6], higher-order gradient models [7,8] and integral type nonlocal models [9–14]. All these methods explicitly or implicitly introduce a material characteristic length to control the width of the localization band, thus prevent strain from localizing into infinitely narrow zones. Nonlocal integral formulation modifies softening plasticity [10,15,16] and damage models [9] by introducing nonlocal variables, which are space-weighted averages of local variables. Earlier versions of nonlocal integral formulation have been shown to be efficient in regularizing ill-posed problems associated with

* Corresponding author. Tel.: +1 213 740 0609; fax: +1 213 744 1426. E-mail address: [email protected] (J.-P. Bardet). 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.08.002

damage, but to require an over-nonlocal formulation [12,17] for softening plasticity. An over-nonlocal variable is defined by linearly combining a local variable with a negative weight and a nonlocal variable with a weight greater than one [18]. The overnonlocal approach was shown to be necessary for wave propagation [17,19] in softening plasticity [12]. In this paper, the spectral analysis [20–22], which is used to study the global stability of solids and structures discretized into finite elements, is applied to examine the mesh dependency of local and over-nonlocal softening plasticity. To our knowledge, it is the first time that spectral analysis is applied to investigate strain localization of discrete systems made of nonlocal plasticity, and to compare the numerical performance of local and nonlocal plasticity. Following the introduction, the second section examines the mesh dependency of local softening plasticity in a one-dimensional example. Following the third section that reviews the over-nonlocal formulation, the fourth section regularizes and derives analytical solutions for plastic strain distributions and load–displacement responses, which are useful for validating numerical solutions. The fourth section also examines the convergence of numerical solutions in the case over-nonlocal approach.

3703

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

2. Mesh dependency in local plasticity 2.1. Rate form of local plasticity Most rate-independent elasto-plastic materials have the following yield function:

Fðr; jÞ ¼ f ðrÞ  j;

ð1Þ

where r is the stress tensor; and j is the yield limit which can be hardening (or softening). The rate-type stress–strain relationship is

r_ ¼ De : ðe_  e_ p Þ;

ð2Þ

where De is the elastic moduli tensor. The plastic strain rate e_ p is assumed to follow an associative rule, i.e.,

_ ;r : e_ p ¼ kF

ð3Þ

The loading–unloading condition is:

F k_ ¼ 0;

k_ P 0;

F 6 0:

ð4Þ

When j depends on the accumulated plastic strain eep ðeep ¼ ffi R qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p p 2deij deij =3dtÞ, the differential consistency condition is:

F ;r : r_ þ F ;j j;eep eep;ep : F ;r k_ ¼ 0;

Fig. 1. Stress–strain relationship of strain-softening plasticity.

ð5Þ

where ‘‘:” denotes a tensor product contracted on two inner indices. The plastic multiplier is:

k_ ¼

F ;r : De : e_ ; F ;r : De : F ;r þ H

ð6Þ

where H is plastic tangential modulus. In case of J2 plasticity [23], pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi f ðrÞ ¼ re ¼ 3sij sij =2 where the deviatoric stress is sij ¼ rij  dij rkk =3. H becomes

H ¼ F ;j j;eep eep;ep : F ;r ¼ j;eep ;



De 

De : F ;r  F ;r : De F ;r : De : F ;r þ H

w L Fig. 2. One-dimensional tension bar with a weak part in the center.

ð7Þ

H is positive for hardening and negative for softening. Using Eqs. (3), (6) and (7), (2) becomes

r_ ¼ Dep : e_ ¼

Δδ



: e_ :

ð8Þ

Eq. (8) is generally used to obtain the tangential stiffness of incremental equilibrium equations.

nite steps. At the beginning time t, all quantities are known, and the solution must be computed at the end of time t þ Dt, corresponding to time increment Dt. The bar is divided into a finite number of elements with a linear displacement interpolation, i.e., the shape function is N ¼ ððx  x2 Þ=h; ðx  x1 Þ=hÞ, where h is the element size; x1 and x2 are the nodal coordinates. The bar is divided into 7, 15 and 31 elements, which correspond to coarse, medium and fine mesh, respectively (Fig. 3). The elements and nodes are numbered from left to right, and strain localization is as-

2.2. Numerical solutions In one-dimensional cases, the internal variable eep is equal to the plastic strain ep and His a negative constant in linear softening plasticity. By integrating the rate-type stress–strain relationship (Eq. (8)), the stress–strain relationship becomes

(

r¼

e < rE0 elastic; cEe þ H r0 ; e P rE0 plastic; Ee;

cE

ð9Þ

where E is the elastic modulus; r0 is tensile strength; and H is plastic modulus, c ¼ H=ðE þ HÞ. As shown in Fig. 1, in the elastic regime, Young’s modulus E governs the relation between stress and strain. After r0 is attained, the stress softens linearly with strain. Fig. 2 shows the one-dimensional problem used hereafter to investigate strain localization. For illustration purposes, the geometrical and material parameters of the bar are arbitrarily selected to be L ¼ 10 cm; A ¼ 10 cm2 ; E ¼ 21 GPa; H ¼ 2:1 GPa; r0 ¼ 10 Mpa. Without loss of generality, we assume the strain localization initiates from the center, i.e., x ¼ L=2. The bar is fixed in one side, i.e., x ¼ 0 and are subjected to an imposed displacement Dd monotonically increasing in other side, i.e., x ¼ L. The problem is solved incrementally using the finite element method. The imposed displacement Dd is divided into smaller fi-

Fig. 3. Three kinds of meshes, course, medium and fine in one-dimensional tension bar, black points denote nodes.

3704

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

sumed to initiate at the center element. The incremental equilibrium equations lead to the conventional finite- element matrix equation:

K Du ¼ Df ;

ð10Þ

where Du is the incremental nodal displacement; Df is the out-ofbalance force, i.e., difference between external and internal nodal forces; K is the global tangential stiffness matrix, which is asseme bled from the element stiffness matrix k : e

k ¼

Z

BT CB dX;

ð11Þ

X

where X denotes the whole domain; B is the strain–displacement matrix, i.e., B ¼ @N=@x ¼ f1=h; 1=hg; C is the modulus depending on the state of the bar. C ¼ E for elastic states, while C ¼ EH=ðE þ HÞ ¼ cE for plastic states. The integral in Eqs. (11) and (12) are calculated using Gaussian quadrature [24]. When there is one integration point per element, the element stiffness matrix is: e

k ¼2

 T     h 11 11 AC 1 1 ¼ AC   : 2 hh hh h 1 1

ð12Þ

Similarly to K; Df can be assembled using element unbalanced force vectors. Here, body forces are neglected and the load is displacement-controlled. The element unbalanced force vector is e

fDfp g ¼ kpq g q  f1; 1gT Artn ;

Fig. 4. Plastic strain distribution after the initiation of strain localization in local plasticity.

ð13Þ

where the subscript q corresponding to the node with a prescribed displacement; rtn denotes the stress at the last equilibrium state. Using Eqs. (10), (12) and (13), the incremental nodal displacements fDug can be found, as well as the incremental strain:

De ¼ BfDug:

ð14Þ

For elastic states, De is totally elastic; the stress is obtained using Hooke’s law. For plastic states, De is the sum of elastic strain Dee and plastic strain Dep , where Dep can be found by linearizing the rate stress–strain relationship Eq. (2) as follows:

DrðxÞ ¼ EDee ¼ EðDe  Dep Þ:

ð15Þ

By applying the consistency condition at and at time F t ¼ rt  ðr0 þ Hept Þ ¼ 0 F tþDt ¼ rtþDt  ðr0 þ HeptþDt Þ ¼ 0, we get

Dr  HDk ¼ 0:

time t, t þ Dt,

i.e., i.e.,

ð16Þ e

p

Using Eqs. (15) and (16), De and De are:

(

H Dee ¼ EþH De ¼ cDe; E Dep ¼ Dk ¼ EþH De ¼ cHE De:

ð17Þ Fig. 5. Load–displacement responses from local plasticity.

Using Eqs. (16) and (17), the tangent modulus C is:

C¼

Dr H ¼ ¼ cE: De ð1 þ H=EÞ

ð18Þ

Eq. (18) implies that the tangent modulus is spatial independent and constant within a plastic zone. The stress-equilibrium condition implies that the stress must be constant across the entire bar. The stress in the bar initially increases linearly with applied displacement. After the stress reaches its peak, strain localization is initiated from the bar center; the bar is in a loading state inside the band, and an unloading state outside the band. Due to softening inside the band, the stress in the bar decreases. When the stress is reduced to 0:5r0 , the plastic strain distributions along the bar (Fig. 4) depend on the mesh refinement; the band width becomes narrower as the mesh is refined. Influenced by the distribution of plastic strain, the global load–displacement curve (Fig. 5) is also mesh sensitive. The coarse mesh

produces negative slope, which is negative as for the stress–strain relationship, while medium and fine meshes revert to a positive slope (snapback). The underlying reason for this mesh dependency will now be investigated using a spectral approach. 2.3. Spectral analysis All the unknown variables, e.g., stresses, total strains and plastic strains are obtained from the nodal displacements Du, which themselves can be expressed in terms of the eigenvectors of stiffness matrix K in Eq. (10). Hereafter we review the spectral analysis for elastic and elasto-plastic states. First in the case of elastic states, K is obtained by assembling the element stiffness matrix ke (Eq. (11)) and enforcing the boundary conditions. Eq. (10) can be written explicitly as

3705

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

38 9 8 2 1 0 0 0 0 0 > u2 > > 0 > > > > > 7> 6 > > > > 6 1 2 1 0 > >0 > 0 0 0 7> u3 > > > > > 7> 6 > > > > > 7> > 6 > >0 > 0 0 7> 6 0 1 2 1 0 > > > > u4 > = < < 7 6 EA 6 7 0 1 2 1 0 0 7 u5 ¼ 0 6 0 > 7> h 6 > > > > > 6 0 > > >0 > u6 > 0 0 1 2 1 0 7 > > 7> 6 > > > > > 7> 6 > > > > > > 7 6 0 > > 0 0 0 1 2 1 5> u > > > 4 7 > >0 > > > : : ; > 0 0 0 0 0 1 1 u8 Df 2

9 > > > > > > > > > > > > > = : > > > > > > > > > > > > > ;

ð19Þ

ð20Þ

where the coefficient ai can be obtained from the inner product of nodal displacement and eigenvector, accounting Tfor the fact that v v eigenvector are orthogonal and normalized, i.e., kv iT k kv jj k ¼ dij : i

Du vi : kDuT k kv i k T

ai ¼

2 1 0 0 0 0 0 6 1 2 1 0 0 0 0 7 7 6 7 6 7 6 0 1 1 þ c c 0 0 0 7 6 EA 6 7 K¼ 0 7: 0 c 1 þ c 1 0 6 0 7 h 6 6 0 0 0 1 2 1 0 7 7 6 7 6 0 0 0 1 2 1 5 4 0 0

All the eigenvalues ki and eigenvectors v i of K can be calculated. Here, ki ’s are numbered according to their absolute value (1 corresponding to the smallest absolute values). In a spectral analysis, the nodal displacements can be expressed as a linear combination of the fundamental eigenvectors.

Du ¼ ai v i ;

3

2

ð21Þ

Table 1 lists the values of coefficient ai corresponding to the first eigenvectors v 1 ; v 2 ; v 3 , and v 4 . As shown in the particular case of Table 2, the values of ai decrease rapidly with i, which implies that Du can be expressed using only the first eigenvectors. The other eigenvectors ði > 4Þ contribute very little to Du and can be neglected. Table 1 list also the four smallest absolute eigenvalues k1 ; k2 ; k3 , and k4 (in addition to coefficients a1 ; a2 ; a3 , and a4 Þ. Fig. 6a–c plots the corresponding eigenvectors ðv 1 ; v 2 ; v 3 ; v 4 Þ in the case of coarse, medium and fine meshes. The first eigenvector v 1 , which contributes the most to the solution, are referred to as the dominant eigenvector. The different meshes in Fig. 7a give almost the same dominant eigenvector, which implies that the nodal displacement of Fig. 7b is mesh-independent for elastic solutions. When the center element becomes plastic, the global tangential stiffness K, which is obtained by introducing one elasto-plastic elee ment stiffness k4 into the global elastic stiffness matrices, becomes:

0

0

0

0

1

ð22Þ

1

Fig. 8a–c shows the eigenvectors v 1 ; v 2 ; v 3 ; v 4 for coarse, medium and fine mesh, whereas Table 1 list the corresponding eigenvalues k1 ; k2 ; k3 , and k4 and coefficients a1 ; a2 ; a3 , and a4 . Compared to the elastic eigenvectors of Fig. 6, the elasto-plastic eigenvectors of Fig. 8a–c display abrupt variations corresponding to inhomogeneous and localized deformation. These abrupt jumps in the eigenvectors become narrower when the mesh is refined, which corresponds to a narrower zone of strain localization. As listed in Table 1, the first four coefficients ai and eigenvalues ki indicate that the dominant eigenvector always corresponds to the smallest absolute eigenvalue, which may be positive or negative. As shown in Fig. 9a, the dominant eigenvectors obtained from three meshes are quire different. In contrast to the coarse mesh, the medium and fine meshes have dominant eigenvectors that deform in the same direction on each side of the localization zone; these modes correspond to a snapback. When the mesh is refined, the eigenvalues shift from positive to negative values and the dominant eigenvectors vary in shape, which imply that the incremental solution (Fig. 9b) becomes mesh-dependent. In the absence of a characteristic length, local plasticity confines plastic zones that are controlled by the size of weak element, which induces boundary value problems that are mathematically ill-posed and need to be regularized. 3. Integral nonlocal averaging Integral nonlocal plasticity is a kind of regularization which introduces a characteristic length. Nonlocal variables are obtained through a spatial-weighted averaging of either scalar or vector variables. The nonlocal average of a local field is defined as

^f ðxÞ ¼

Z

aðx; nÞf ðnÞ dn;

ð23Þ

V

Table 1 Coefficients ai and corresponding eigenvalues ki (in parenthesis) of first four eigenvectors for different meshes using elasticity, local plasticity and nonlocal plasticity. Material type

Number of elements

a1 ðk1 =ðEAÞÞ

a2 ðk2 =ðEAÞÞ

a3 ðk3 =ðEAÞÞ

a4 ðk4 =ðEAÞÞ

Elasticity

7 15 31

0.9932 (0.0306) 0.9928 (0.0154) 0.9927 (0.0077)

0.1080 (0.2674) 0.1100 (0.1376) 0.1103 (0.0692)

0.0384 (0.7000) 0.0394 (0.3770) 0.0400 (0.1917)

0.0179 (1.2537) 0.0199 (0.7237) 0.0205 (0.3739)

Local plasticity

7 15 31

0.9821 (0.0308) 0.9004 (0.0313) 0.9838 (0.0109)

0.1408 (0.0981) 0.3360 (0.0600) 0.1729 (0.0897)

0.1222 (0.3814) 0.2681 (0.1954) 0.0135 (0.1252)

0.0031 (1.0625) 0.0159 (0.5354) 0.0207 (0.2474)

Over-nonlocal plasticity

7 15 31

0.9997 (0.0053) 0.9996 (0.0026) 0.9996 (0.0013)

0.0215 (0.0869) 0.0264 (0.0420) 0.0279 (0.0209)

0.0061 (0.1576) 0.0056 (0.1207) 0.0058 (0.0628)

0.00014 (0.2283) 0.0020 (0.2208) 0.0020 (0.1245)

Table 2 Coefficients ai of all eigenvectors for 8-node mesh using elasticity, local plasticity and nonlocal plasticity. Material

a1

a2

a3

a4

a5

a6

a7

Elasticity Local Nonlocal

0.9932 0.9821 0.9997

0.1080 0.1408 0.0215

0.0384 0.1222 0.0061

0.0179 0.0031 0.00014

0.0105 0.0256 0.0111

0.0053 2.9889e4 4.6319e5

0.0042 0.0078 0.0021

3706

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

Fig. 7. Dominant eigenvector and normalized incremental nodal displacement in coarse, medium and fine meshes.

where V is the volume of the entire body; aðx; nÞis the averaging function, it should be positive, symmetric about x, and ensure consistency of order 0 (i.e., reproducibility of constant functions). In order to satisfy these conditions, the averaging function is defined as

R a ðx; nÞf ðnÞ dn aðx; nÞ ¼ RV 1 : a ðx  nÞ dn V 1

ð24Þ

Eq. (24) guarantees that a constant field is not modified due to nonlocal averaging. a1 can be Gaussian function, bi-linear exponential function and bell-shaped function. Here we choose bi-linear exponential function



Fig. 6. First eigenvectors for coarse, medium and fine meshes (elasticity).

a1 ¼ exp 

 jx  nj : l

ð25Þ

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

3707

Fig. 9. Dominant eigenvector and normalized incremental nodal displacement in coarse, medium and fine meshes (local plasticity).

Because the averaging zone is limited to an area of finite extent, the computational efficiency can be improved by using a truncated averaging function.

^f ðxÞ ¼

Z

xþR

aðx; nÞf ðnÞ dn;

ð26Þ

xR

where R denotes the effective influence area of nonlocal averaging, which is related to a characteristic length. Over-nonlocal variables are obtained by combing local and nonlocal variables as follows

Fig. 8. First eigenvectors in coarse, medium and fine meshes (local plasticity).

^f m ðxÞ ¼ ð1  mÞf ðxÞ þ m^f ðxÞ Z xþR ¼ ½ð1  mÞdðx; nÞ þ maðx; nÞf ðnÞ dn; xR

ð27Þ

3708

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

where dðn  xÞ is dirac delta function; mðm > 1Þ is an over-nonlocal weighting factor. For m ¼ 0, an over-nonlocal variable degenerates into a local one. For 0 < m 6 1, the formulation is the same as standard nonlocal formulation. 4. Strain localization in over-nonlocal plasticity 4.1. Analytical solutions We use over-nonlocal model to re-analyze the bar problem in Section 3. After replacing the internal variable by its over-nonlocal counterpart, Eq. (1) becomes:

F ¼ r  rð^epm Þ ¼ r  r0  H^epm ¼ r  ðr0 þ ð1  mÞHep þ mH^ep Þ ¼ 0; ^pm

ð28Þ p

^p

where e is the over-nonlocal plastic strain; e and e are the local and nonlocal plastic strain, respectively. For uniform strain, local, nonlocal and over-nonlocal formulations generate identical stress–strain relationships (i.e., Eq. (9)). When the strain localizes far away from the bar extremities, boundary effects can be ignored. Using Eqs. 9, 8, 7, 6, 5, 4, 3, 2 and 28, the consistency condition becomes:

r  r0 Hð1  mÞ

¼ ep þ

m 2lð1  mÞ

Z

xL=2þR

jnðxL=2Þj l

e

ep ðnÞ dn:

ð29Þ

xL=2R

u¼

Eq. (29) can be transformed into the following second order ordinary differential equation [25]:

@ 2 ep ðxÞ 1 r  r0 þ 2 ep ðxÞ ¼ 2 : @x2 Hl ð1  mÞ l ðm  1Þ

ð30Þ

Plastic strain is assumed to develop symmetrically about the bar center, i.e., x ¼ L=2. Considering only the symmetric terms, the solution of Eq. (30) is:



x  L=2



ep ¼ Q cos pffiffiffiffiffiffiffiffiffiffiffiffiffi þ l m1

r  r0 H

;

ð31Þ

where Q is an integration constant, which can be determined after imposing the following condition at the elastic–plastic transition at x ¼ L=2  w=2.



Fig. 10. The influence of m on band width of strain localization.

r E

Lþ2

Z

L=2þw=2

r ðr  r0 Þ ep dx ¼ L þ 2pd:

ð37Þ

H

E

L=2

After combing Eqs. (37), (35) becomes

ep ¼

u  rE0 L 1 ; f ðxÞ 2pd þ HL E

L=2  pd 6 x 6 L=2 þ pd;

u>

r0 E

L:

ð38Þ

The global force and displacement response is

F¼

! u þ 2pHdr0 A: 2pd þ EL H

ð39Þ

Snapbacks can be avoided and strain localization can be confined within the bar by enforcing the following conditions

HL=ð2pEÞ < d < L=ð2pÞ:

ð40Þ

p

e ðL=2  w=2Þ ¼ 0; e_ p ðL=2  w=2Þ ¼ 0;

ð32Þ

where w denotes the width of the localization zone. Eq. (32) yields:

Q ¼ ðr  rt Þ=H;

ð33Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffi w ¼ 2pl m  1 ¼ 2pd;

ð34Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffi where d ¼ l m  1. As shown in Fig. 10, the normalized bandwidth w=l is independent of stress level and plastic modulus and stays constant during the deformation process. Using Eqs. (31), (33) and (34), the plastic strain is

ðr  r0 Þ 1 ; H f ðxÞ

L=2  pd 6 x 6 L=2 þ pd;

ð35Þ



. Eq. (35) implies that the plastic strain where f ðxÞ ¼ 1 þ cos xL=2 d increases monotonically from zero at the elastic–plastic transition to a maximum value at the center of the plastic region. By adding the elastic strain, the total strain is

(

e¼

r þ ðrr0 Þ 1 E

H

r; E

f ðxÞ

;

The stress is updated as in the case of local plasticity theory [23] by enforcing the consistency condition at time t, i.e.:

F t ¼ rt  ðr0 þ ð1  mÞHept þ mH^ept Þ ¼ 0:

and

ep ¼

4.2. Numerical solutions

L=2  pd 6 x 6 L=2 þ pd; x < L=2  pd; or x > L=2 þ pd:

ð36Þ

The integration of the stress–strain response along the entire bar produces the complete force–displacement response, i.e.:

ð41Þ

And at time t þ Dt:

F tþDt ¼ rtþDt  ðr0 þ ð1  mÞHeptþDt þ mH^eptþDt Þ ¼ 0:

ð42Þ

Combining Eqs. (41) and (42), we get

Dr  ðð1  mÞHDkðxÞ þ mHD^kðxÞÞ ¼ 0:

ð43Þ

After adopting the bi-linear exponential averaging function, Eq. (43) turns into

Dr m ¼ DkðxÞ þ 2lð1  mÞ ð1  mÞH

Z

xL=2þR

jnðxL=2Þj l

e

DkðnÞ dn:

ð44Þ

xL=2R

The plastic multiplier Dk can be obtained by solving a Fredholm equation of second kind [25]. The same as Eqs. (29), (44) can be transformed into the following second order ODE:

@ 2 DkðxÞ 1 Dr þ 2 DkðxÞ ¼ 2 : @x2 Hl ð1  mÞ l ðm  1Þ

ð45Þ

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

3709

After imposing a symmetry about the bar center x ¼ L=2, the solution of Eq. (45) becomes:

Dk ¼

Dr 1 ; H f ðxÞ

L=2  pd 6 x 6 L=2 þ pd:

ð46Þ

Using Eqs. (17) and (46), the tangent modulus CðxÞ becomes:

CðxÞ ¼

DrðxÞ ¼ Ecf ðxÞ: DeðxÞ

ð47Þ

In contrast to the local plasticity solution, the nonlocal plasticity solution has a tangent modulus C that varies spatially in the plastic area. This spatial variation is denoted using f ðxÞ in Eq. (47). The

Fig. 11. Plastic strain distributions in nonlocal plasticity (band width is 6.28 cm).

Fig. 12. Load–displacement responses from nonlocal plasticity.

Fig. 13. First eigenvectors in coarse, medium and fine meshes (nonlocal plasticity).

3710

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711

nodal displacement is obtained using Eqs. (11) and (47), which yields the spatial distribution of plastic strain as well as the global load–displacement response. As shown in the numerical results of Fig. 11, the bandwidth and plastic strain distribution remain almost identical for different meshes. As shown in Fig. 12, the load–displacement responses displays in case of coarse mesh a strong softening whereas it shows very slight difference in case of medium and fine meshes; the numerical solution converges to the analytical solution when the mesh is refined.

2

2

1

0

0

0

0

0

0

0

3

6 1 1 þ f c f2 c 0 0 0 0 7 7 6 2 7 6 7 6 0 c ðf þ f Þ c f c 0 0 0 f 2 1 2 1 7 6 EA 6 7 K¼ ðf1 þ f2 Þc f2 c 0 0 7; 0 f1 c 6 0 7 h 6 6 0 1 þ f2 c 1 0 7 0 0 f2 c 7 6 7 6 0 0 0 1 2 1 5 4 0 0

0

0

1

1 ð48Þ

4.3. Spectral analysis In the case of over-nonlocal regularization, the tangential stiffness matrix K in the plastic area is different from that in Eq. (22). For a coarse mesh (8 nodes), K becomes:

where the coefficients f1 and f2 are calculated at the centers of elements 3–5, respectively, i.e.:

Fig. 15. The number of negative eigenvalue in different meshes: local and nonlocal plasticity.

Fig. 14. Dominant eigenvector and normalized incremental nodal displacement in coarse, medium and fine meshes (nonlocal plasticity).

Fig. 16. The band width of different meshes: local and nonlocal plasticity.

X. Lu et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3702–3711



f1 ¼ f ðL=2Þ ¼ 0:5; f2 ¼ f ðL=2  hÞ ¼ 0:8651:

ð49Þ

Fig. 13 shows the eigenvectors corresponding to eigenvalues k1 ; k2 ; k3 and k4 for three different meshes containing 8, 16, 32 nodes, respectively. When the mesh is refined, the eigenvectors become smoother while displaying a steady transition area of constant width, which is associated with strain localization. Table 1 lists the coefficients ai and corresponding eigenvalues ki of the first four eigenvectors v i . Higher-order modes have small coefficients ai and can be neglected in expanding the nodal displacements Du. As shown in the three meshes of Fig. 14a, the dominant eigenvectors always correspond to the least negative eigenvalues. Fig. 14b shows that the incremental nodal displacements are almost identical to the dominant eigenvectors. The number of negative eigenvalue increases with mesh refinement (Fig. 15), while it stays equal to one in local plasticity for any mesh size. In local plasticity, the size of the plastic zone depends on the size of the weak element; the size of the plastic area becomes smaller for finer mesh, and tends to zero when mesh becomes infinitely small. In nonlocal plasticity, when the mesh is refined, more elements are involved in defining a plastic zone of constant size. As shown in Fig. 16, the width of strain localization can be obtained from the eigenvectors; it converges to the analytical solution when the mesh becomes finer. By introducing a characteristic length, over-nonlocal softening plasticity confines the localized deformation into a finite band, which is constant and independent to the mesh size. 5. Conclusions The analytical and numerical solutions of one-dimensional strain localization problems have been obtained for strain-softening local and nonlocal plasticity models. When the mesh is refined, strain-softening local plasticity generates plastic strains that localize into narrower and narrower band, and produce load–displacement response that exhibit excessive strain-softening, even snapback. Such numerical shortcoming was investigated using the spectral analysis of the incremental equilibrium equations. In case of a weak element, the tangential stiffness matrix was found to have only one negative eigenvalue, and the corresponding eigenvector to display an abrupt discontinuity across the weak element; the dominant eigenvector, which controls the deformation mode, varies substantially with the mesh size. This finding implies that the solution is mesh-dependent and the boundary-valued problem is ill posed, as one would expect in the absence of characteristic length. This ill-posed problem was regularized using over-nonlocal softening plasticity, and produced numerical solutions that are mesh-independent. In the case of nonlocal plasticity, the spectral analysis revealed that the number of negative eigenvalue increases with mesh refinement, which means more elements become plastic and can accommodate strain localization. All the eigenvectors, including those corresponding to negative eigenvalues, display smooth patterns. The few eigenvectors that display localized patterns, reminiscent of strain localization, have a width constant and independent of the mesh size. Numerical solutions agree with the conclusion of spectral analysis; the plastic strain distribution and load–displacement response were both found to converge to analytical

3711

solutions when the mesh was refined. This one-dimensional analysis demonstrates that over-nonlocal plasticity can regularize illposed BVPs and alleviate the problems associated with strain softening and mesh dependency in finite element analysis. Additional research is needed to extend our conclusion to two- and threedimensional conditions. 6. Acknowledgments This work was financially supported by National Science Foundation of China (NSFC through Grant No. 50825803) and Chinese Scholarship Council of Oversea Studies. These supports are gratefully acknowledged. We are grateful to the anonymous reviewers for their helpful comments and suggestions. References [1] A. Needleman, Material rate dependence and mesh sensitivity in localization problems, Comput. Methods Appl. Mech. Engrg. 67 (1) (1988) 69–85. [2] H.L. Schreyer, M.K. Neilsen, Analytical and numerical tests for loss of material stability, Int. J. Numer. Methods Engrg. 39 (10) (1996) 1721–1736. [3] L.J. Sluys, R.d. Borst, Wave propagation and localization in rate dependent crack medium, model formulation and one-dimensional examples, Int. J. Solids Struct. 29 (1992) 2945–2958. [4] H.B. Mühlhaus, I. Vardoulakis, The thickness of shear bands in granular materials, Géotechnique 37 (1987) 271–283. [5] J.P. Bardet, J. Proubet, A numerical investigation of the structure of persistent shear bands in granular media, Géotechnique 41 (4) (1991) 559–613. [6] M.I. Alsaleh, G.Z. Voyiadjis, K.A. Alshibli, Modelling strain localization in granular materials using micropolar theory: mathematical formulations, Int. J. Numer. Anal. Methods Geomech. 30 (2006) 1501–1524. [7] R. de Borst, H.B. Mühlhaus, Gradient-dependent plasticity: formulation and algorithmic aspects, Int. J. Numer. Methods Engrg. 35 (1992) 521–539. [8] C. Comi, Computational modelling of gradient-enhanced damage in quasibrittle materials, Mech. Cohesive Frict. Mater. 4 (1) (1999) 17–36. [9] Z.P. Bazˇant, T. Belytschko, T.P. Chang, Continuum theory for strain softening, J. Engrg. Mech., ASCE 110 (1984) 1666–1692. [10] Z.P. Bazˇant, F.B. Lin, Nonlocal yield-limit degradation, Int. J. Numer. MethodsEngrg. 26 (1988) 1805–1823. [11] G. Pijaudier-Cabot, Z.P. Bazˇant, Non-local damage theory, J. Engrg. Mech. ASCE 113 (10) (1987) 1512–1533. [12] L. Strömberg, M. Ristinmaa, FE-formulation of a nonlocal plasticity theory, Comput. Methods Appl. Mech. Engrg. 136 (1996) 127–144. [13] C. Comi, A non-local model with tension and compression damage mechanisms, Eur. J. Mech. A/Solids 20 (2001) 1–22. [14] Z.P. Bazˇant, M. Jirásek, Nonlocal integral formulations of plasticity and damage: survey of progress, J. Engrg. Mech., ASCE 128 (11) (2002) 1119–1149. [15] A.C. Eringen, On nonlocal plasticity, Int. J. Engrg. Sci. 19 (1981) 1461–1474. [16] A.C. Eringen, Theories of nonlocal plasticity, Int. J. Engrg. Sci. 21 (1983) 741– 751. [17] G.D. Luzio, Z.P. Bazant, Spectral analysis of localization in nonlocal and overnonlocal materials with softening plasticity or damage, Int. J. Solids Struct. 42 (23) (2005) 6071–6100. [18] P.A. Vermeer, R.B.J. Brinkgreve, A new effective non-local strain measure for softening plasticity, in: Localization and Bifurcation Theory for Soil and Rocks, Balkema, Rotterdam, 1994. [19] M. Jirásek, S. Rolshoven, Comparison of integral-type nonlocal plasticity models for strain-softening materials, Int. J. Engrg. Sci. 41 (2003) 1553–1602. [20] J.P. Bardet, Finite element analysis of surface instability in hypo-elastic solids, Comput. Methods Appl. Mech. Engrg. 78 (1990) 273–296. [21] J.P. Bardet, Finite element analysis of plane strain bifurcation within compressible solids, Comput. Struct. 36 (6) (1990) 993–1007. [22] R. de Borst, Computation of post-bifurcation and post-failure behavior of strain softening solids, Comput. Struct. 25 (1987) 211–224. [23] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer-Verlag, New York, 1998. [24] T.B. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, New York, 2000. [25] A.D. Polyanim, A.V. Manzhirov, Hand Book of Integral Equations, CRC Press, Boca Raton, FL, 1998.