- Email: [email protected]
A higher-order strain gradient plasticity theory with a corner-like effect
International Journal of Solids and Structures 58 (2015) 62–72
Contents lists available at ScienceDirect
International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr
A higher-order strain gradient plasticity theory with a corner-like effect Mitsutoshi Kuroda ⇑ Graduate School of Science and Engineering, Mechanical Engineering, Yamagata University, Yonezawa, Yamanata 992-8510, Japan
a r t i c l e
i n f o
Article history: Received 20 November 2014 Available online 25 December 2014 Keywords: Finite strains Size effect Shear bands Flow localization Finite element method
a b s t r a c t A corner-like plasticity model originally proposed as a size-independent theory is extended to include a size effect resulting from plastic strain gradients. A method of solving boundary value problems at finite strains is also presented. The efficiency of the new theory is demonstrated through two typical numerical examples: a constrained simple shear problem in which an infinitely long strip bounded by two hard materials is subjected to large shear under plane strain conditions, and a problem of shear band formation in plane strain tension. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Whether or not a vertex (or a corner) develops on the yield surface of a plasticity material is a very important inquiry in predictions of plastic instability that accompanies an abrupt change in the deformation mode without a large variation in the stress state (Støren and Rice, 1975; Hutchinson and Tvergaard, 1981). Plasticity theories accounting for the effects of a yield surface corner have been proposed by Christoffersen and Hutchinson (1979) and Gotoh (1985). Later, for the purpose of developing simplified and efficient numerical procedures, Hughes and Shakib (1986) proposed a pseudo-corner theory representing the reduced stiffness and increased plastic flow of the corner-theory-like response to nonproportional loading. Simo (1987) proposed a non-associative flow rule that represents a corner-like effect on an apparently smooth yield surface. An experimental investigation on the use of an abrupt strain path change to determine the shape of a subsequent yield surface in the vicinity of a current loading point was carried out by Kuwabara et al. (2000), which was based on a method proposed by Kuroda and Tvergaard (1999). The recorded trajectory of stress immediately after a strain path change, which inevitably travels close to the current yield surface, strongly indicated the existence of a yield surface vertex at the loading point. Furthermore, the direction of the plastic strain rate was clearly inclined toward the forward direction of the stress path that is regarded as part of the section of the current yield surface. That is, a large deviation of the plastic strain rate direction from the normal to the apparent ⇑ Tel.: +81 238 26 3211; fax: +81 238 26 3205. E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijsolstr.2014.12.019 0020-7683/Ó 2014 Elsevier Ltd. All rights reserved.
yield surface was observed. This is interpreted to mean that when the stress point moves along what appears to be a smooth yield surface, a vertex moves with the stress point, which explains the apparent nonnormality. These experimental observations are similar to the numerical predictions obtained using the Taylor polycrystal model (Kuroda and Tvergaard, 1999). The observed apparent nonnormality effect on a smooth yield surface is consistent with the basic idea of Simo’s corner-like plasticity model (Simo, 1987). Kuroda and Tvergaard (2001a) modified Simo’s model so that it can represent observations in the polycrystal model (Kuroda and Tvergaard, 1999) and the experiments of Kuwabara et al. (2000) as closely as possible at a macroscopic level of modeling. The modified model of Kuroda and Tvergaard (2001a) was applied to finite element analysis to predict the shear band development in plane strain tensile specimens, and it was confirmed that the results were reasonably close to crystal plasticity predictions (Kuroda and Tvergaard, 2001c). Conventional plasticity theories including the aforementioned classes of corner theories do not account for any size effect in real material behavior. A wide array of experiments on micron-size specimens have revealed significant size-dependent mechanical behaviors in plastically strained materials involving spatial gradients of strain (e.g., Fleck et al., 1994; Stölken and Evans, 1998). On the theoretical side, a considerable number of studies have been conducted with the aim of incorporating the strain gradient effects into theories of plasticity since the pioneering studies of Aifantis (1984, 1987). Part of them were carried out within a context of the crystal plasticity framework (e.g., Gurtin, 2002; Evers et al., 2004; Kuroda and Tvergaard, 2006; Borg, 2007). However, in parallel, extensions of phenomenological plasticity theories (in most cases, J2-based theories) have attracted considerable
63
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
attention owing to their simplicity and practical efficiency (e.g., Fleck et al., 1994; Fleck and Hutchinson, 2001; Gudmundson, 2004; Gurtin and Anand, 2009; Hutchinson, 2012). The former approach is advantageous for performing physically based simulations directly accounting for the effects of the geometrically necessary dislocations (GNDs), which correspond to gradients of crystallographic slips (Ashby, 1970). The latter approach is suitable in cases where one intends to model a resultant size effect in polycrystalline materials or introduce an assumption of isotropy, as a first approximation, for micron-scale plasticity. It is now widely recognized that these gradient plasticity theories must be higherorder in the sense that it should be possible to impose extra boundary conditions with respect to plastic strains or their gradients, which are outside the scope of conventional plasticity theories. In the present study, the corner-like plasticity model of Kuroda and Tvergaard (2001a), whose original version was proposed by Simo (1987), is extended to include a size effect resulting from plastic strain gradients (Aifantis, 1984, 1987). Then, a method of solving boundary value problems at finite strains using the proposed constitutive model is presented. Two typical examples are shown to demonstrate the efficiency of the new theory: one is a constrained simple shear problem in which an infinitely long strip bounded by two hard materials is subjected to large shear under plane strain conditions, and the other is a problem of shear band formation in plane strain tension. 2. Theory and solution method 2.1. Constitutive modeling Let x be the current position of a particle labeled X in the undeformed configuration of the body under consideration, F ¼ @x=@X _ _ be the deformation gradient, L ¼ F_ F1 ¼ @ x=@x ¼ u_ r ¼ @ u=@x be the velocity gradient (where u is the displacement, a superposed dot denotes a material-time derivative, r is the spatial gradient operator with respect to the current configuration, and is the tensor product operator), D ¼ 12 ðL þ LT Þ be the deformation rate tensor, and W ¼ 12 ðL LT Þ be the continuum spin tensor. An additive decomposition of the deformation rate tensor into elastic and plastic parts is postulated as De ¼ D Dp . The elastic response is assumed to be given by the following hypoelasticity relation
r ¼ r_ W r þ r W ¼ C : De
ð1Þ
with
C ¼ k1 1 þ 2lIð4sÞ ;
ð2Þ
where r is the Cauchy (true) stress, r is its Jaumann rate, C is a fourth-order elastic moduli tensor, k and l are Lamé constants, 1 is the second-order identity tensor, and I(4s) is the fourth-order symmetric identity tensor. For plastic deformation, the following form of the flow rule is adopted: p
p
_ ; D ¼ /N
ð3Þ
R(ep) is a function of the equivalent plastic strain ep, which represents a work-hardened state of the material, b is a length scale coefficient assumed to be constant in the present study, a prime, ðÞ0 , denotes the deviatoric part of the tensor ðÞ, and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jðÞj ¼ tr½ðÞT ðÞ. The introduction of the gradient term into yield conditions was mainly motivated by the desire to predict the postlocalization features of material behavior. The inclusion of the gradient term is necessary to determine the shear band width in the post-bifurcation regime as first discussed by Aifantis (1984, 1987). Although more complicated or generalized strain gradient formulations have been proposed by different researchers (e.g., Fleck and Hutchinson, 2001; Gudmundson, 2004), Eq. (4) is used here as the first choice. Eq. (4) is the simplest, but involves the primary effect of the plastic strain gradients. A physical basis for the introduction of r gp ð¼ br2 ep for a constant b) can be strengthened by an argument based on the dislocation theory, as discussed by Kuroda and Tvergaard (2006, 2010). That is, a dislocation-induced longrange internal stress arises in response to spatial gradients of the GND density (Groma et al., 2003; Evers et al., 2004), and the GND density is equated with the spatial gradient of crystallographic slip (Ashby, 1970). Thus, the dislocation-induced internal stresses correspond to the second gradients of crystallographic slips. The macroscopic yield condition with the term br2 ep is mathematically similar to microscopic yield conditions for each slip system in the gradient crystal plasticity theories (Groma et al., 2003; Evers et al., 2004). Based on this view, we can identify the term br2 ep with a resultant of the GND density-induced internal stress. On this understanding, plastic dissipation should be accounted for by Re_ p ðP 0Þ as pointed out by Kuroda and Tvergaard (2010). The direction tensor Np in Eq. (3) is taken to be
Np ¼ n þ^d m
ð7Þ
with
n ¼
@f =@ r r0 ; ¼ j@f =@ rj jr0 j
D0 ðn : D0 Þ n
m ¼
ð8Þ
jD0 ðn : D0 Þ n j
The framework of this flow rule with a corner-like effect was originally proposed by Simo (1987). In the original model, ^ d was formulated so that Dp is always coaxial to D0 when D0 lies inside a hypercone defined by semi-angle Hpcrit (measured from the direction n; Dp must lie within the hypercone). Later, Kuroda and Tvergaard ^ (2001a) modified the relation for d to
( ^d ¼ tan Hp ;
1
H ¼ cos
p
H ¼
" # n : D0 ; jD0 j
aH
H
p crit
a¼
c
for aH 6 Hpcrit for aH > Hpcrit Rðep Þ
l
;
ð9Þ
1 þ1
;
ð10Þ
p
where /_ is a scalar plastic multiplier and the tensor N represents the direction of Dp , which will be defined later. A gradient-dependent yield condition is introduced, following Aifantis (1984), as
f ¼ re þ r gp Rðep Þ ¼ 0
ð4Þ
for a plastic loading condition, and f < 0 for elastic states (i.e. before initial yielding or in an unloaded state), where
re ¼
rffiffiffi 3 0 jr j; 2
gp ¼ brep ;
1 3
r0 ¼ r r : ð1 1Þ;
ð5Þ ð6Þ
where c is a coefficient that introduces non-coaxiality between the total deviatoric strain rate D0 and the plastic strain rate Dp . A schematic diagram of the present corner-like model is shown in Fig. 1. The introduction of non-coaxiality between D0 and Dp is a necessity to reproduce the corner effect observed in crystal plasticity. Kuroda and Tvergaard (2001a) suggested that a range of c from 3 to 10 might be realistic through identification using Taylor model computations involving an abrupt strain path change. In a subsequent study on shear band simulations (Kuroda and Tvergaard, 2001c), it was found that a slightly smaller value, c = 2, was most suitable for reproducing the shear band development behavior predicted in crystal plasticity simulations. The present corner-like model
64
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
with
v ¼ n gp ;
ð15Þ
where n is a unit vector normal to the surface Sp of the plastic loading zone in the current configuration. Using the relations dV = JdV0, p J = det F, r ¼ r0 F1 , and ndSp ¼ n0 F1 JdS0 (known as Nanson’s formula) with subscript ‘‘0’’ referring to the reference configuration, Eq. (14) is rewritten as
Z p
V0
fdwðJR J re Þ þ r0 dw ðJF1 gp ÞgdV p0
¼
Z Sp0
dwn0 ðJF1 gp ÞdSp0 :
ð16Þ
For the incremental analysis of elastoplasticity problems, Eq. (16) is differentiated with respect to time, noting that _ ðF1 Þ ¼ F1 L and J ¼ J tr L, and again using JdV p0 ¼ dV p :
Fig. 1. Illustration of the corner-like effect adopted in the present plasticity model (Kuroda and Tvergaard, 2001a)).
Z hn o i R_ r_ e þ ðtr LÞðR re Þ dw þ fg_ p L gp þ ðtrLÞgp g rdw dV p Vp Z ¼ v^_ dwdSp Sp
ð17Þ
reduces to the J2-flow theory when Hpcrit ¼ 0 or c ! 1, while Simo’s original model is recovered for c = 0 (i.e. a = 1) and Hpcrit > 0 . Considering the relation p
r : D ¼ jr j/_ ¼ re 0
vb ¼ n fg_ p L gp þ ðtr LÞgp g:
rffiffiffi 2_ /; 3
ð11Þ
the equivalent plastic strain is naturally defined as
e_ p ¼
rffiffiffi 2_ /; 3
ep ¼
Z
t
e_ p dt;
ð12Þ
0
where t is time.1 In Eq. (1), the standard Jaumann stress rate has been used in the hypoelastic constitutive assumption. In the study of Kuroda and _ p with Tvergaard (2001b), the relation, W ¼ x þ Wp ¼ x þ /X Xp ¼ gðr n n rÞ, for the plastic spin for the present corner-like model was proposed and discussed through an extensive investigation of isotropic as well as anisotropic materials, followed by a comparison with experimental observations. According to this model, the plastic spin vanishes when n is coaxial to r0 , as in the present model. Thus, the spin of the substructure, x, is identified with the continuum spin, i.e. x = W. 2.2. Strategy for solving boundary value problems 2.2.1. Gradient-dependent yield condition The yield condition of Eq. (4) contains the second-order gradient term r gp , and thus this equation is viewed as a partial differential equation (PDE), unlike a standard yield condition in conventional plasticity theories. To get solutions to Eq. (4), its weak form is derived from an integral equation
Z Vp
fre þ r gp Rðep ÞgdwdV p ¼ 0;
ð13Þ
where V p is the volume of the plastic loading zone in the current (deformed) configuration and dw is an arbitrary scalar weighting function. Eq. (13) is integrated by parts to give an expression
Z Vp
1
fðR re Þdw þ gp rdwgdV p ¼
with
Z Sp
vdwdSp
ð14Þ
The t is not necessarily the ‘‘natural time’’ in the present theory since the constitutive model is time-independent.
ð18Þ
In Eq. (17), all integrations are performed in the current configuration, and thus this relation is suitable for an updated Lagrangian _ r_ e ; and g_ p ; computation scheme. The material-time derivatives, R; in Eq. (17) are calculated straightforwardly as
dR R_ ¼ p e_ p ¼ he_ p ; de
r_ e ¼
@ re @ re :r¼ :C:D @r @r
ð19Þ rffiffiffi 3 @ re : C : Np e_ p ; 2 @r
g_ p ¼ ðbrep Þ ¼ bre_ p gp L
ð20Þ ð21Þ
Eq. (14) with Eq. (15) and Eq. (17) with Eq. (18) enable a deduction of nonstandard boundary conditions, which are extra conditions compared with conventional elastoplasticity problems. The surface Sp consists of two parts: one is an external part Spext that belongs to the external boundary surface of the whole domain under consideration; the other is an internal part Spint that is an interface with the elastic (or unloaded) zone. Unlike the external boundary Spext , the internal boundary Spint moves inside the body during elastic–plastic deformation. Here, we define the plastic strain gradient flux as dF ¼ vdSp ¼ n gp dSp . When we consider a free external surface with no constraint on plastic flow, it may be acceptable to set zero plastic strain gradient flux, i.e. dF ¼ 0, _ ¼ 0. The latter is accomplished by or equivalently dF v^_ ¼ n fg_ p L gp þ ðtr LÞgp g ¼ 0 in an incremental analysis based on Eq. (17). The treatment of the conditions at the internal elastic– plastic boundary Spint may not be so simple since this boundary moves in the domain.2 We assume that there is no explicit constraint on e_ p at an elastic–plastic boundary. The corresponding condition for the plastic strain gradient flux should be considered only in an incremental manner. In the present paper, it is mimicked by _ ¼ 0), similarly to the external boundary, but ^_ ¼ 0 (i.e. dF setting v the condition dF ¼ 0 in the total form no linger holds at Spint . At a 2 If we consider a viscoplastic material in which plastic flow is potentially active everywhere, neither internal elastic–plastic boundary nor extra internal boundary condition appears in the problem (e.g., Borg et al., 2006).
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
boundary Sp interfacing with a hard or conceptually rigid material, zero plastic strain, e_ p ¼ 0, is a natural choice as a boundary condition. 2.2.2. Equilibrium condition The equilibrium condition in conventional solid mechanics remains the same in the present theory:
rr¼0
ð22Þ
with the body force effect neglected for brevity. The principle of virtual work expressed in the current configuration reads as a weak form of the PDE of Eq. (22),
Z
r : ðdu_ rÞdV ¼
V
Z
t du_ dS
ð23Þ
65
Isoparametric eight-node quadrilaterals with quadratic shape _ functions are used for the solution of the displacement rates u, while four-node quadrilaterals with linear shape functions are used for the solution of the equivalent plastic strain rate e_ p . The volume integrations are carried out by a 2 2 Gaussian integral. The structure of the finite element equations, Eq. (27), is similar to that of the equations explored by Mühlhaus and Aifantis (1991), Niordson and Hutchinson (2003) and Niordson and Redanz (2004). The relationship between the present formulation and those in the previous studies will be discussed in Section 4. 3. Numerical examples 3.1. Constrained simple shear with boundary layers
S
where t is the surface traction, V is the volume of the analysis domain or any part of it, and S is the surface of V (here, n is a unit normal to S). Differentiating the above equation through the same procedure as for Eq. (14), we obtain
A strip with a thickness H0 in the X2-direciton and an infinite length in the X1-direction is subjected to simple shear under a plane strain condition. Periodicity is assumed for all field quantities, taking an arbitrary periodic length L0 in the X1-direction. Consequently, values of the field quantities vary only in the thickness direction. The top and bottom surfaces are bounded by hard materials. The imposed boundary conditions are
Z
u_ 1 ¼ 0; u_ 2 ¼ 0; e_ p ¼ 0 on X 2 ¼ 0;
ð28Þ
~_ u_ 2 ¼ 0; e_ p ¼ 0 on X 2 ¼ H0 ; u_ 1 ¼ U;
ð29Þ
with
t ¼ n r;
ð24Þ
T
fr_ L r þ ðtrLÞrg : ðdu_ rÞdV ¼
V
Z
^t_ dudS _
ð25Þ
S
with
^t_ ¼ n fr_ L r þ ðtrLÞrg:
ð26Þ
Substituting the constitutive relations, Eqs. (1) and (3), gives the basis for finite element equations in an updated Lagrangian scheme. 2.2.3. Finite element analysis Based on Eqs. (25) and (17), finite element discretized equations can be derived in a standard manner, taking the displacement rate u_ and the equivalent plastic strain rate e_ p as independent variables, as
KðuuÞ KðpuÞ
( _ ) ( _ ) CðuÞ FðuÞ U ¼ þ : KðppÞ CðpÞ e_ p F_ ðpÞ KðupÞ
ð27Þ
_ is a vector type array of nodal displacement rates and fe_ p g Here, fUg is that of nodal equivalent plastic strain rates. The upper and lower groups of equations in Eq. (27) are derived from Eqs. (25) and (17), respectively, except for fCðuÞ g and fCðpÞ g, which are equilibrium correction terms calculated from the total form equations (23) and (14), respectively. Details of Eq. (27) are given in Appendix. Actual numerical computations are performed using a forward Euler time integration scheme in which increments of all quantities are taken to be DðÞ ¼ ðÞ _ Dt, and the time increment Dt is adaptively determined using the rmin method proposed by Yamada et al. (1968). If e_ p is computed to be negative despite the assumption of continuous plastic loading, elastic unloading is said to occur. In such a case, Dt is reduced to one-tenth of the value determined by the rmin method. As can be seen from Eq. (8), the present constitutive model is not linear in D. In the actual computations, the value of D used to construct Np is taken to be the value determined in the previous increment without iteration3 as in the study of Kuroda and Tvergaard (2001c). This is a good approximation for the small time increments used in the forward Euler method. The same kind of approximation has been employed in calculations based on the J2 corner theory (Hutchinson and Tvergaard, 1981; Tvergaard et al., 1981). 3 In Kuroda and Tvergaard (2001a), an iterative method was used for the same constitutive model.
~_ is the prescribed displacement rate in the X1-direction on where U the top surface. The macroscopic shear strain is defined by ~ 0. C ¼ U=H The strain hardening function R(ep) is chosen to be
ep n Rðep Þ ¼ r0 1 þ ;
e0
ð30Þ
where r0 is the initial yield stress, e0 is the offset strain, and n is the hardening exponent. The material parameter values used are E/ r0 = 500, m = 0.3, e0 = 0.002, and n = 0.1, where E is the Young’s modulus and m is the Poisson’s ratio. These values are the same as those used in Kuroda and Tvergaard (2001b).4 The length scale coefficient b is defined as
b¼l
2
r0 ;
ð31Þ
where l is a length-scale parameter. In this problem, plastic yielding occurs simultaneously at all material points, and no elastic unloading occurs during the subsequent increase in the shear deformation. Therefore, no internal boundary between the elastic and elastic–plastic zones emerges. Fig. 2 shows curves of stresses versus macroscopic shear strain for several typical cases. Results for uniform simple shear without a strain gradient, which are obtained by employing an extra ^_ ¼ 0 instead of e_ p ¼ 0 in Eqs. (28) and (29), boundary condition v are also shown for reference. In this case, all size effects are excluded, i.e. the results are completely identical to those for the conventional size-independent theory with b = 0. The parameter c, which governs the non-coaxiality between D0 and Dp , is taken as c = 2–3, and Hpcrit is chosen to be 20°. A remarkable feature of the present corner-like plasticity model is the generation of considerable amounts of normal stresses. In the case of the J2-flow theory (i.e. Hpcrit = 0°), the normal stresses are close to zero.5 The amounts of normal stresses predicted by the present model are consistent 4 In Kuroda and Tvergaard (2001b), a viscoplastic version of the present model was used with a very small value of the rate sensitivity parameter (i.e. m = 0.002) in a power law relation, which exhibited nearly rate-independent behavior. 5 These very small normal stresses are due to the use of the hypoelastic assumption (i.e. Eq. (1)).
66
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
2 l /H
1.5
l /H0 = 0.3
ij /
Stress
0
12
1 0.5
The deformed mesh at C = 1 in the case of l/H0 = 0.3 is illustrated in Fig. 4(b), which shows the effect of the extra boundary condition e_ p ¼ 0 on the overall deformation mode of the strip.
.5
=0 0
3.2. Shear band development in plane strain tension
Conventional
c =3,
p
c =2,
p
J2F (
crit =
20
crit = 20 p crit = 0 )
0 22
-0.5
l /H
0=
0.5
l /H0 = 0.3
A rectangular specimen under plane strain tension with X2 being the tensile axis is considered. The initial dimensions of the specimen are 2H0 along the X1-direction and 2L0 along the X2direction, and the aspect ratio L0/H0 is taken to be 3. Considering the geometrical symmetry, only one quadrant of the specimen is analyzed. The boundary conditions for the quadrant are
^_ ¼ 0 on X 1 ¼ 0; u_ 1 ¼ 0; ^t_ 2 ¼ 0; v
ð32Þ
^_ ¼ 0 on X 2 ¼ 0; u_ 2 ¼ 0; ^t_ 1 ¼ 0; v
ð33Þ
^t_ 1 ¼ 0; ^t_ 2 ¼ 0; v ^_ ¼ 0 on X 1 ¼ H0 þ DH0 ;
ð34Þ
_ ^t_ 1 ¼ 0; v ^_ ¼ 0 on X 2 ¼ L0 ; u_ 2 ¼ U;
ð35Þ
-1 0
0.5
1 1.5 2 Macroscopic shear strain
2.5
Fig. 2. Curves of stresses vs macroscopic shear strain for constrained simple shear problem of gradient-dependent corner-like plasticity solids. Stress responses in the absence of a gradient effect (designated as ‘conventional’) are also depicted for comparison. ‘‘J2F’’ stands for the J2-flow theory.
with the amounts of axial stresses observed in fixed-end torsion tests as discussed by Kuroda and Tvergaard (2001b). Next, we examine size effects in constrained simple shear. For the computations with finite length scales, a finite element mesh of 1 (in the X1-direction) 30 (in the X2-direction) quadrilaterals is used, which leads to the convergence of solutions on the plotted graph scale. It can be seen in Fig. 2 that large amounts of additional hardening are produced with increasing length scale (l/H0 = 0.3 and 0.5). Note that the stress components r22 and r12 must be uniform in the thickness direction due to the equilibrium requirement from Eq. (22). The J2-flow theory (i.e. Hpcrit = 0°) produces only insignificant normal stresses even in this size dependent setting. The stress distributions at C = 1 for the corner-like plasticity solids (Hpcrit = 20° and c = 2) with the length scales of l/H0 = 0.3 and 0.5 are plotted in Fig. 3. Although the distributions of the normal stresses, r11 and 11 and r 33 , r33, are strongly nonuniform, their volume averages, r 11 ¼ r22 and r 33 ¼ 0, respectively, simisatisfy the conditions r larly to the observation in the case of uniform simple shear. The stress distributions observed here are very similar to those for the finite strain gradient crystal plasticity model reported in Kuroda and Tvergaard (2008). The distributions of ep are plotted in Fig. 4(a). Those for l/H0 = 0.3 and 0.5 are close to each other.
1
Position X2 / H0
l /H0 = 0.3 l /H0 = 0.5
0.5
12 11
0 -1
-0.5
0
0.5 Stress
DH0 ¼ H0 ½n1 cosðpX 2 =L0 Þ þ n2 cosðmw pX 2 =L0 Þ
ð36Þ
with amplitudes of the imperfection waves, n1 and n2, and the inte^_ ¼ 0 is also valid on ger wave number mw. Note that the condition v the symmetry planes as specified in Eqs. (32) and (33). The macro 0 . The nondimensionalscopic nominal strain is defined by eN ¼ U=L ized tensile load T is defined by
T¼
P ; H 0 r0
ð37Þ
where P is the sum of the nodal reaction forces at X2 = L0 for the one-quadrant model with a unit depth in the X3-direction. Eq. (30) is employed for the strain hardening function R(ep). The material parameter values used are E/r0 = 333, m = 0.3, e0 = 0.0015, and n = 0.1, which are the same as those used in the study of Kuroda and Tvergaard (2001c). The controlling parameters for the corner-like effect are taken as Hpcrit X= 20° and c = 2 unless stated otherwise. These values were used in the study of Kuroda and Tvergaard (2001c) in the conventional size-independent case. For the length-scale coefficient b, Eq. (31) is also employed in this problem. The extra boundary conditions at internal boundaries between the elastic zone and the elastic–plastic zone enter this problem since elastic unloading occurs in a broad area at the onset of diffuse necking and the unloading region gradually develops with increasing overall elongation of the specimen. In the present simulations, no special numerical treatment is applied at the internal elastic– n o plastic boundaries, i.e. no input is given to F_ ðpÞ in Eq. (27) for the corresponding elements. This approximately satisfies the con^_ ¼ 0 on the internal boundaries. As a result, all the comdition of v n o ponents of F_ ðpÞ remain zero since the external boundary
33
22
_ is the prescribed rate of end-displacement and DH0 reprewhere U sents the initial geometrical imperfection of the form (Tvergaard et al. 1981; Kuroda and Tvergaard, 2001c)
1 ij /
1.5
2
0
Fig. 3. Distributions of stresses at a macroscopic shear strain of C = 1.0 in the constrained simple shear problem of gradient-dependent corner-like plasticity solids (Hpcrit = 20° and c = 2) with length scales of l/H0 = 0.3 and 0.5.
conditions, as well as the symmetry conditions, are all given by v^_ ¼ 0 as in Eqs. (32)–(35). Niordson and Redanz (2004) reported that in their plane strain sheet necking analysis using the Fleck– Hutchinson’s higher-order theory, the two conditions n re_ p ¼ 0 (in their formulation) and e_ p ¼ 0 made no significant difference. First, we examine the effect of the mesh density on strain localization predictions. Fig. 5 displays curves of computed tensile load versus macroscopic nominal strain for the corner-like plasticity solids (Hpcrit = 20° and c = 2) with and without the length-scale
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a) 1
67
(b) l /H0 = 0.3
Position
X2 / H0
l /H0 = 0.5
= 0.5
0.5
1.0
1.5
2.0
2.5
0 0
0.5 1 1.5 Equivalent plastic strain
2
2.5
p
Fig. 4. (a) Distributions of equivalent plastic strain ep at various stages of macroscopic shear strain C in the constrained simple shear problem of gradient-dependent cornerlike plasticity solids (Hpcrit = 20° and c = 2) with length scales of l/H0 = 0.3 and 0.5; (b) deformed mesh with contours of ep for gradient-dependent corner-like plasticity model (Hpcrit = 20° and c = 2) with l/H0 = 0.3 at C = 1.
Nondimensionalized tensile load T
2 l / H0 = 0.04
1.5 Conventional
1 30 x 90 quadrilaterals 12 x 36 quadrilaterals
0.5
0 0
0.1 0.2 Macroscopic nominal strain
0.3 N
Fig. 5. Computed curves of tensile load vs macroscopic nominal strain for plane strain tension problem employing corner-like-effect parameters of Hpcrit = 20and c = 2 with and without the length-scale effect. Two different mesh densities (12 36 and 30 90) are considered. In all the computations, the initial geometric imperfection is specified by n1 = 0.005 and n2 = 0.
effect, employing two finite element discretizations: a mesh of 12 36 quadrilaterals and a mesh of 30 90 quadrilaterals6 with a pure cosine mode imperfection (n1 = 0.005 and n2 = 0). Deformed meshes and contours of equivalent plastic strain at stages where strain localizations in the form of a shear band have been well developed are depicted in Fig. 6. In the conventional size-independent case (l = 0), significant mesh dependence appears as expected (see Fig. 6 (a) and (b)). The shear band that appears in the mesh of 30 90 quadrilaterals is much more severe than that appearing in the mesh of 12 36 quadrilaterals. The tensile load T for the fine mesh begins to decrease slightly earlier than that for the coarse mesh as shown in Fig. 5. In the finite length-scale case of l/ H0 = 0.04, the widths of the shear band for both the finite element discretizations are almost the same (Figs. 6(c) and (d)). The corre-
6 In both the finite element discretizations, the initial size of each element in the X2-direction is given by a geometric progression with a common ratio of 1.02, so that elements are more densely assigned in the potential neck region than in the anticipated unloading region.
sponding curves of T versus eN for these two computations are also almost identical (Fig. 5). Fig. 7 displays curves of computed tensile load T versus macroscopic nominal strain eN for a somewhat complex imperfection with n1 = 0.005, n2 = 0.001, and mw = 8. The computations are carried out for several different length scales using the mesh of 30 90 quadrilaterals. Deformed meshes and contours of equivalent plastic strain at stages where strain localizations have been well developed are depicted in Fig. 8. In the conventional sizeindependent case (l = 0), the most severely strained point is on the free surface at X2 = 0 and a localized band is formed from the free surface (Fig. 8(a)). With an increase in the length scale relative to the specimen size, the patterns of the flow localization zone drastically change. For l/H0 = 0.01, a severely strained spot also appears at the center of the specimen, and the localized band from the free surface ceases to develop. The deformation mode for l/ H0 = 0.04 (Fig. 8(c)) becomes almost identical to that in the case of the pure cosine imperfection (Fig. 6(d)). Thus, the gradient term introduced into the present corner-like plasticity model governs not only the ‘‘width’’ of shear bands, but also the patterns of strain localization. For l/H0 = 0.1, the shear band completely disappears and only the diffuse neck-type deformation is observed. Finally, the effect of c is illustrated in Fig. 9. Here, c is increased to 3 and the other settings are unchanged from the previous computations whose results are shown in Figs. 7 and 8. Compared with the results shown in Fig. 8 for c = 2, significantly different shear band patterns are observed for c = 3. The most severely strained point appears at the center of the specimen even in the size-independent case (Fig. 9(a)). For larger length scales, l/H0 = 0.04 and 0.1, the results for c = 3 are rather similar to those for c = 2 with the same length scales (Fig. 8 (c) and (d)).
4. Discussion 4.1. On modeling of corner plasticity The key to the successful extension from the conventional sizeindependent theory to the strain gradient-enhanced version in the present study is that the constitutive relation for Dp is given by a form similar to the conventional flow rule, i.e. the product of the qffiffi
direction tensor Np and the plastic multiplier /_ ¼ 32e_ p as is seen in Eq. (3) with Eq. (12). Owing to this simple constitutive structure,
68
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a)
(b)
(c)
(d)
Fig. 6. Results of mesh dependence tests with and without length scale effect (the corner-like-effect parameters are taken to be Hpcrit = 20and c = 2) for plane strain tension problem showing deformed meshes and contours of the equivalent plastic strain ep for (a) a size-independent (conventional) model with a mesh of 12 36 quadrilaterals at eN = 0.240; (b) a conventional model with a mesh of 30 90 quadrilaterals at eN = 0.240; (c) the gradient-dependent model of l/H0 = 0.04 with a mesh of 12 36 quadrilaterals at eN = 0.270; (d) the gradient-dependent model of l/H0 = 0.04 with a mesh of 30 90 quadrilaterals at eN = 0.270. In all the computations, an initial geometric imperfection specified by n1 = 0.005 and n2 = 0 (a pure cosine mode) is introduced.
Nondimensionalized tensile load T
2
1.5
Conventional
1
l / H0 = 0.01 l / H0 = 0.04
0.5
l / H0 = 0.1
0 0
0.1 0.2 Macroscopic nominal strain
0.3 N
Fig. 7. Computed curves of tensile load vs macroscopic nominal strain for plane strain tension problem, employing corner-like-effect parameters of Hpcrit = 20° and c = 2 with different length scales. The initial geometric imperfection is specified by n1 = 0.005, n2 = 0.001, and mw = 8, and a mesh of 30 90 quadrilaterals is used in all the computations.
the extension to the strain gradient version could be carried out in a similar manner to the case of the strain gradient theory based on the classical J2 flow theory. The other corner theories (i.e. Christoffersen and Hutchinson (1979) and Gotoh (1985)) assumes
the dependence of the plastic strain rate direction on the stress rate, which are not given in a simple form like Eq. (3) with Eq. (12). Although the actual formulations are not pursued in this paper, the extension of these corner theories to corresponding higher-order strain gradient versions would not be so straightforward unlike the case of the present theory. Under conventional size-independent postulations, these corner theories as well as the present theory have basically similar predictability for plastic flow localization phenomena. In fact, the present corner-like theory provides a prediction of shear band formation very similar to that by the Christoffersen and Hutchinson’s theory, depending on the choice of the parameter values, as illustrated in Kuroda and Tvergaard (2001c). Thus, in the context of enhancement of the theories with introduction of the size effect, the present corner-like higher-order gradient plasticity theory may have some worth as an alternative to other corner-type theories whose extension to strain gradient versions is not straightforward. 4.2. On predictability of plastic flow localization As shown in Tvergaard et al. (1981), the initial geometrical imperfection mode plays an important role in determining the pattern of shear band formation. In the present study, the same finding can be seen on comparison between the results shown in Figs. 6(b) and 8(a). The difference in these two computations was only the modes of the initial geometrical imperfection, and the both did not involve the size effect. It is clearly seen in Fig. 8 that
69
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a)
(b)
(c)
(d)
Fig. 8. Effects of the length scale on patterns of strain localization (the corner-like-effect parameters are taken to be Hpcrit = 20° and c = 2) in plane strain tension problem, introducing an initial geometric imperfection specified by n1 = 0.005, n2 = 0.001, and mw = 8 showing deformed meshes and contours of the equivalent plastic strain ep for (a) a size-independent (conventional) model at eN = 0.230; (b) l/H0 = 0.01 at eN = 0.250; (c) l/H0 = 0.04 at eN = 0.274; (d) l/H0 = 0.1 at eN = 0.300. A mesh of 30 90 quadrilaterals is used in all the computations.
the introduction of the gradient effect varies not only the width of shear band but also the fundamental patterns of flow localization. This is a new point found in the present study with respect to the shear band behavior. 4.3. On formulations and computations of strain gradient plasticity In the present paper, the terms ‘higher-order stress’ and higherorder traction’ have not appeared so far at all, which were among the central items in several previous studies on strain gradient plasticity. In this section, the connection between the present treatment of the higher-order strain gradient effects and those considered by other researchers is discussed. In the present study, the length-scale coefficient b has been taken as a constant for simplicity. In this case, the term r gp in Eq. (4) can be written as br2 ep . The introduction of the vector quantity gp would be essential for mathematical transformation to different expressions of the theory as shown below. In Eq. (14), formally replacing dw by a virtual equivalent plastic strain rate de_ p , we have
Z Vp
fðR re Þde_ p þ gp rde_ p gdV p ¼
Z Sp
vde_ p dSp :
ð38Þ
Now, the following constraint for the virtual quantities is considered:
ðdu_ rÞsym ¼ dD ¼ dDe þ dDp ¼ dDe þ
rffiffiffi 3 p p de_ N : 2
ð39Þ
Upon noting that de_ p is only taken into account within the plastic zone and considering Eqs. (11) and (12), Eq. (38) with Eq. (39) is rewritten as
Z
ðr : dDp þ Rde_ p þ gp rde_ p dV ¼
V
Z
vde_ p dS:
ð40Þ
S
Adding Eq. (40) to the virtual work principle (Eq. (23)) gives
Z V
fðr : dDe þ Rde_ p þ gp rde_ p gdV ¼
Z S
_ þ t dudS
Z
vde_ p dS: ð41Þ S
The relation (41) is the same as the modified virtual work statement in Fleck and Hutchinson’s strain gradient plasticity theory (Fleck and Hutchinson, 2001). Eq. (41) was introduced as a major assumption at the starting point of their theory. In Eq. (41), gp (in their paper a quantity equivalent to gp is denoted by s) was introduced as a work conjugate to the plastic strain gradient rep . Fleck and Hutchinson referred to it as a higher-order stress vector quantity. Similarly, v was introduced as a work conjugate to the equivalent plastic strain ep on the surface, called a higher-order traction quantity. In the study of Fleck and Hutchinson (2001), the major premise of Eq. (41) yields the conventional equilibrium condition (Eq. (22)) and the gradient dependent yield condition (Eq. (4)), as opposed to the present formulation of gradient plasticity, where all the PDEs are given initially based on physics and physically based insights. It is obvious that the present theory can be equivalent to an approach based on a relation in the form of Eq. (41), as long as the same constitutive relations are postulated for Dp and gp.
70
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a)
(b)
(c)
(d)
Fig. 9. Deformed meshes and contours of equivalent plastic strain ep for plane strain tension problem when c = 3 with different length scales (Hpcrit remains at 20°), introducing an initial geometric imperfection specified by n1 = 0.005, n2 = 0.001, and mw = 8: (a) size-independent (conventional) model at 0.250; (b) l/H0 = 0.01 at 0.250; (c) l/ H0 = 0.04 at 0.274; (d) l/H0 = 0.1 at 0.300. A mesh of 30 90 quadrilaterals is used in all the computations.
Gurtin and Anand (2009) addressed a thermodynamically consistent derivation of the constitutive relations for the higher-order quantities in the Fleck and Hutchinson’s framework (Eq. (41)): i.e. in their argument, if we assume that a defect energy, wp , of the form wp ¼ 12 bjrep j2 augments the free energy (per volume), the Clausius–Duhem inequality leads to the constitutive equation gp ¼ brep that is identical to Eq. (6). According to this view, we can claim that gp is recoverable (or energetic in terminology by Gurtin and Anand (2009)); thus, plastic dissipation must occur via Re_ p ðP 0Þ. This view is consistent with the argument presented in Section 2.1, which is based on the dislocation theory. An incremental version of Eq. (41) is derived by adding Eq. (17) to Eq. (25) with de_ p that replaces dw
Z n o ^_ r ^_ T : ðdu_ rÞ þ ðR ^_ p rde_ p dV ^_ e Þde_ p þ g P V Z Z _ _ ^_ de_ p dS þ v ¼ ^t dudS S
ð42Þ
S
with definitions
^_ ¼ r_ L r þ ðtrLÞr; P
ð43Þ
^_ ¼ R_ þ ðtr LÞR; R
ð44Þ
r^_ e ¼ r_ e þ ðtr LÞre ;
ð45Þ
^_ p ¼ g_ p L gp þ ðtr LÞgp : g
ð46Þ
_ are viewed as the nominal These quantities with superposed ‘^’ rates of the corresponding quantities at the current configuration. Eq. (42) is identified with the incremental virtual work principle of Niordson and Redanz (2004), which was derived on the basis of the Fleck and Hutchinson’s theory (a small strain theory). The treatment of the quantity gp (in their study it is a higher-order stress, s) is, however, somewhat different. In their model, a constitutive rela^_ p in an incremental form, which was tion regarding gp was given to g deduced from the Fleck and Hutchinson’s theory, while in the present theory, the relation for gp is defined in a total form (Eq. (6)) as in the original Aifantis model. It is emphasized that Eq. (42) can give exactly the same finite element equations (Eq. (27)) as long as all the constitutive relations are the same. The finite element formulation presented earlier in Niordson and Hutchinson (2003) was a small strain version of that in Niordson and Redanz (2004). The basic structure of Eq. (27) is also the same as the formulation of Mühlhaus and Aifantis (1991), which was the earliest proposed finite element procedure for the higher-order gradient plasticity, although it included a small strain assumption. Mühlhaus and Aifantis (1991) considered an incremental potential functional accounting for the gradient effects. Applying the stationary condition to the potential functional gives an equation very similar to the incremental form of Eq. (41), except for the higherorder traction term on the right-hand side. Mühlhaus and
71
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
Aifantis (1991) did not explicitly define any higher-order quantity, and they only considered the situations v_ = 0 or e_ p = 0 (under a small strain setting), assuming that the product of these quantities is always zero. Kuroda and Tvergaard (2010) discussed a different mathematical treatment of higher-order gradient plasticity theories, in which the increment of the plastic strain gradient is taken as an additional independent variable instead of the increment of the plastic strain, without introducing the concept of higher-order stress. This idea can be applied to the present theory in principle. However, in the case of large strain problems, the treatment of the extra bound^_ ¼ n fg_ p L gp þ ðtr LÞgp g ¼ 0, is not so simple ary condition, v when g_ p is chosen as an independent variable. To avoid this complexity, the solution strategy shown in section 2.2 has been adopted in the present computations. On the computational side, all the numerical simulations performed here have adopted quadratic shape functions for the displacement rates and linear shape functions for the equivalent plastic strain rates. Empirically, this combination gives high-quality numerical solutions with fast convergence with reasonable mesh densities and with no stress oscillation phenomenon as found in Fredriksson et al. (2009) and Kuroda (2011) for some higher-order gradient plasticity computations. A more extensive and systematic investigation regarding the quality of finite element solutions for the various formulations of strain gradient plasticity is needed towards further expanding the use of strain gradient plasticity. This is left for a subsequent study.
_ ¼ fu_ 1 ; u_ 2 gT at a material point inside the element, we write fug _ _ ¼ ½NfUg, the finite element discretized relations as fug _ (with fDg ¼ fD11 ; D22 ; 2D12 gT ), fLg ¼ ½BL fUg, _ fDg ¼ ½BD fUg and _ where [BD], [BL] and [divB] contain spatial derivatr L ¼ ½divBfUg, tives of the shape functions in appropriate arrays. Similarly, a shape function matrix [Ne] relates the nodal equivalent plastic strains fe_ p g to the equivalent plastic strain at a material point inside the element as e_ p ¼ ½N e fe_ p g, and its spatial gradients, fre_ p g ¼ f@ e_ p =@x1 ; @ e_ p =@x2 g, are computed as fre_ p g ¼ ½Be fe_ p g with the matrix [Be] containing spatial derivatives of the components of [Ne]. The orders of the shape functions in [N] for the displacement rates and in [Ne] for the equivalent plastic strain rates are not necessarily the same. Upon the above preparation, the finite element equations in Eq. (27) are derived from Eqs. (25) and (17). Detailed expressions for each component in Eq. (27) are shown as follows. The incremental virtual work principle (Eq. (25)) gives
KðuuÞ ¼
8 9 Z < = ½BD T ½C½BD þ ½BL T ½C r ½BL dV; |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} ; V :
ð49Þ
finite strain effect
where [C] is elastic moduli (Eq. (2)) in matrix form,
½KðupÞ ¼
rffiffiffi Z 3 2l½BD T Np ½Ne dV p 2 Vp
with the plastic strain T N p ¼ N p11 ; N p22 ; Np12 , and
rate
ð50Þ direction
of
n o Z _ F_ ðuÞ ¼ ½NT f^tgdS:
5. Concluding remarks
the
form
ð51Þ
S
The corner-like plasticity model that was originally proposed by Simo (1987) and later modified by Kuroda and Tvergaard (2001a) has been extended to include the size effect resulting from the plastic strain gradient. A method of solving boundary value problems at finite strains has also been given, and the efficiency of the model in problems including higher-order boundary conditions and plastic flow localization has been demonstrated through finite element simulations of problems involving constrained simple shear and plane strain tension. Appendix
Eq. (17) gives
KðpuÞ ¼
8 9 Z < rffiffiffi = 3 T T T T p 2l½N e fng ½BD þ ðR re Þ½N e ½divB þ ½Be ½g ½BD dV p |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}; 2 Vp : finite strain effect
ð52Þ T
with fng ¼ fn11 ; n22 ; n12 g ,
"
p
½g ¼
g p1 g p2
ð47Þ
KðppÞ ¼
Z Vp
g p2
g p1
# ;
n o ðh þ 3lÞ½Ne T ½Ne þ b½Be T ½Be dV p ;
ð53Þ
ð54Þ
Sp
In Eqs. (49) and (52), the terms with the note ‘‘finite strain effect’’ are peculiar to the present finite strain formulation, and would not appear in a small strain formulation. The total form of the virtual work principle (Eq. (23)) yield the equilibrium correction vector
Z Z fCðuÞ g ¼ a ½BD frgdV þ ½NT ftgdS V
ð55Þ
S
with {r} = {r11, r22, r12}T and {t} = {t1, t2}T, and the weak form of the yield function (Eq. (14)) gives
Z CðpÞ ¼ a
Vp
h
Z i ½Ne T ðR re Þ þ ½Be T fg p g dV p þ ½N e T vdSp ; Sp
ð56Þ
ð48Þ
Introducing a matrix array of shape functions, [N], which relates the _ to the displacement rates displacement rates fUg
nodal
g p2
n o Z ^_ dSp : F_ ðpÞ ¼ ½Ne T v
or in a symbolic expression
r þ ½C r fLg:
and
A practical form of Eq. (27) is presented for plane strain condi^_ def tions. The terms P ¼ r_ L r þ ðtrLÞr in Eq. (25) are viewed as a nominal stress rate in the current configuration. With the introduc ^_ is tion of the Jaumann rate r of the Cauchy stress (Eq. (1)), P expressed in a vector–matrix-type expression that may be useful in coding a computer program (e.g., Yamada and Sasaki, 1995): 8 9 > > r11 > > 8 9 > > > > > 2 38 9 _ > > > > > ^ > > > P11 > 0 r12 0 r11 > L11 > > > > > > > > > > > > > > > > > > 6 > > > > 7> > > > > > > r _ > > > > > 6 7 ^ 12 = < 12 = 6 0 1 ðr11 r22 Þ 1 ðr11 þ r22 Þ r12 7< L21 >
Contents lists available at ScienceDirect
International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr
A higher-order strain gradient plasticity theory with a corner-like effect Mitsutoshi Kuroda ⇑ Graduate School of Science and Engineering, Mechanical Engineering, Yamagata University, Yonezawa, Yamanata 992-8510, Japan
a r t i c l e
i n f o
Article history: Received 20 November 2014 Available online 25 December 2014 Keywords: Finite strains Size effect Shear bands Flow localization Finite element method
a b s t r a c t A corner-like plasticity model originally proposed as a size-independent theory is extended to include a size effect resulting from plastic strain gradients. A method of solving boundary value problems at finite strains is also presented. The efficiency of the new theory is demonstrated through two typical numerical examples: a constrained simple shear problem in which an infinitely long strip bounded by two hard materials is subjected to large shear under plane strain conditions, and a problem of shear band formation in plane strain tension. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Whether or not a vertex (or a corner) develops on the yield surface of a plasticity material is a very important inquiry in predictions of plastic instability that accompanies an abrupt change in the deformation mode without a large variation in the stress state (Støren and Rice, 1975; Hutchinson and Tvergaard, 1981). Plasticity theories accounting for the effects of a yield surface corner have been proposed by Christoffersen and Hutchinson (1979) and Gotoh (1985). Later, for the purpose of developing simplified and efficient numerical procedures, Hughes and Shakib (1986) proposed a pseudo-corner theory representing the reduced stiffness and increased plastic flow of the corner-theory-like response to nonproportional loading. Simo (1987) proposed a non-associative flow rule that represents a corner-like effect on an apparently smooth yield surface. An experimental investigation on the use of an abrupt strain path change to determine the shape of a subsequent yield surface in the vicinity of a current loading point was carried out by Kuwabara et al. (2000), which was based on a method proposed by Kuroda and Tvergaard (1999). The recorded trajectory of stress immediately after a strain path change, which inevitably travels close to the current yield surface, strongly indicated the existence of a yield surface vertex at the loading point. Furthermore, the direction of the plastic strain rate was clearly inclined toward the forward direction of the stress path that is regarded as part of the section of the current yield surface. That is, a large deviation of the plastic strain rate direction from the normal to the apparent ⇑ Tel.: +81 238 26 3211; fax: +81 238 26 3205. E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijsolstr.2014.12.019 0020-7683/Ó 2014 Elsevier Ltd. All rights reserved.
yield surface was observed. This is interpreted to mean that when the stress point moves along what appears to be a smooth yield surface, a vertex moves with the stress point, which explains the apparent nonnormality. These experimental observations are similar to the numerical predictions obtained using the Taylor polycrystal model (Kuroda and Tvergaard, 1999). The observed apparent nonnormality effect on a smooth yield surface is consistent with the basic idea of Simo’s corner-like plasticity model (Simo, 1987). Kuroda and Tvergaard (2001a) modified Simo’s model so that it can represent observations in the polycrystal model (Kuroda and Tvergaard, 1999) and the experiments of Kuwabara et al. (2000) as closely as possible at a macroscopic level of modeling. The modified model of Kuroda and Tvergaard (2001a) was applied to finite element analysis to predict the shear band development in plane strain tensile specimens, and it was confirmed that the results were reasonably close to crystal plasticity predictions (Kuroda and Tvergaard, 2001c). Conventional plasticity theories including the aforementioned classes of corner theories do not account for any size effect in real material behavior. A wide array of experiments on micron-size specimens have revealed significant size-dependent mechanical behaviors in plastically strained materials involving spatial gradients of strain (e.g., Fleck et al., 1994; Stölken and Evans, 1998). On the theoretical side, a considerable number of studies have been conducted with the aim of incorporating the strain gradient effects into theories of plasticity since the pioneering studies of Aifantis (1984, 1987). Part of them were carried out within a context of the crystal plasticity framework (e.g., Gurtin, 2002; Evers et al., 2004; Kuroda and Tvergaard, 2006; Borg, 2007). However, in parallel, extensions of phenomenological plasticity theories (in most cases, J2-based theories) have attracted considerable
63
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
attention owing to their simplicity and practical efficiency (e.g., Fleck et al., 1994; Fleck and Hutchinson, 2001; Gudmundson, 2004; Gurtin and Anand, 2009; Hutchinson, 2012). The former approach is advantageous for performing physically based simulations directly accounting for the effects of the geometrically necessary dislocations (GNDs), which correspond to gradients of crystallographic slips (Ashby, 1970). The latter approach is suitable in cases where one intends to model a resultant size effect in polycrystalline materials or introduce an assumption of isotropy, as a first approximation, for micron-scale plasticity. It is now widely recognized that these gradient plasticity theories must be higherorder in the sense that it should be possible to impose extra boundary conditions with respect to plastic strains or their gradients, which are outside the scope of conventional plasticity theories. In the present study, the corner-like plasticity model of Kuroda and Tvergaard (2001a), whose original version was proposed by Simo (1987), is extended to include a size effect resulting from plastic strain gradients (Aifantis, 1984, 1987). Then, a method of solving boundary value problems at finite strains using the proposed constitutive model is presented. Two typical examples are shown to demonstrate the efficiency of the new theory: one is a constrained simple shear problem in which an infinitely long strip bounded by two hard materials is subjected to large shear under plane strain conditions, and the other is a problem of shear band formation in plane strain tension. 2. Theory and solution method 2.1. Constitutive modeling Let x be the current position of a particle labeled X in the undeformed configuration of the body under consideration, F ¼ @x=@X _ _ be the deformation gradient, L ¼ F_ F1 ¼ @ x=@x ¼ u_ r ¼ @ u=@x be the velocity gradient (where u is the displacement, a superposed dot denotes a material-time derivative, r is the spatial gradient operator with respect to the current configuration, and is the tensor product operator), D ¼ 12 ðL þ LT Þ be the deformation rate tensor, and W ¼ 12 ðL LT Þ be the continuum spin tensor. An additive decomposition of the deformation rate tensor into elastic and plastic parts is postulated as De ¼ D Dp . The elastic response is assumed to be given by the following hypoelasticity relation
r ¼ r_ W r þ r W ¼ C : De
ð1Þ
with
C ¼ k1 1 þ 2lIð4sÞ ;
ð2Þ
where r is the Cauchy (true) stress, r is its Jaumann rate, C is a fourth-order elastic moduli tensor, k and l are Lamé constants, 1 is the second-order identity tensor, and I(4s) is the fourth-order symmetric identity tensor. For plastic deformation, the following form of the flow rule is adopted: p
p
_ ; D ¼ /N
ð3Þ
R(ep) is a function of the equivalent plastic strain ep, which represents a work-hardened state of the material, b is a length scale coefficient assumed to be constant in the present study, a prime, ðÞ0 , denotes the deviatoric part of the tensor ðÞ, and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jðÞj ¼ tr½ðÞT ðÞ. The introduction of the gradient term into yield conditions was mainly motivated by the desire to predict the postlocalization features of material behavior. The inclusion of the gradient term is necessary to determine the shear band width in the post-bifurcation regime as first discussed by Aifantis (1984, 1987). Although more complicated or generalized strain gradient formulations have been proposed by different researchers (e.g., Fleck and Hutchinson, 2001; Gudmundson, 2004), Eq. (4) is used here as the first choice. Eq. (4) is the simplest, but involves the primary effect of the plastic strain gradients. A physical basis for the introduction of r gp ð¼ br2 ep for a constant b) can be strengthened by an argument based on the dislocation theory, as discussed by Kuroda and Tvergaard (2006, 2010). That is, a dislocation-induced longrange internal stress arises in response to spatial gradients of the GND density (Groma et al., 2003; Evers et al., 2004), and the GND density is equated with the spatial gradient of crystallographic slip (Ashby, 1970). Thus, the dislocation-induced internal stresses correspond to the second gradients of crystallographic slips. The macroscopic yield condition with the term br2 ep is mathematically similar to microscopic yield conditions for each slip system in the gradient crystal plasticity theories (Groma et al., 2003; Evers et al., 2004). Based on this view, we can identify the term br2 ep with a resultant of the GND density-induced internal stress. On this understanding, plastic dissipation should be accounted for by Re_ p ðP 0Þ as pointed out by Kuroda and Tvergaard (2010). The direction tensor Np in Eq. (3) is taken to be
Np ¼ n þ^d m
ð7Þ
with
n ¼
@f =@ r r0 ; ¼ j@f =@ rj jr0 j
D0 ðn : D0 Þ n
m ¼
ð8Þ
jD0 ðn : D0 Þ n j
The framework of this flow rule with a corner-like effect was originally proposed by Simo (1987). In the original model, ^ d was formulated so that Dp is always coaxial to D0 when D0 lies inside a hypercone defined by semi-angle Hpcrit (measured from the direction n; Dp must lie within the hypercone). Later, Kuroda and Tvergaard ^ (2001a) modified the relation for d to
( ^d ¼ tan Hp ;
1
H ¼ cos
p
H ¼
" # n : D0 ; jD0 j
aH
H
p crit
a¼
c
for aH 6 Hpcrit for aH > Hpcrit Rðep Þ
l
;
ð9Þ
1 þ1
;
ð10Þ
p
where /_ is a scalar plastic multiplier and the tensor N represents the direction of Dp , which will be defined later. A gradient-dependent yield condition is introduced, following Aifantis (1984), as
f ¼ re þ r gp Rðep Þ ¼ 0
ð4Þ
for a plastic loading condition, and f < 0 for elastic states (i.e. before initial yielding or in an unloaded state), where
re ¼
rffiffiffi 3 0 jr j; 2
gp ¼ brep ;
1 3
r0 ¼ r r : ð1 1Þ;
ð5Þ ð6Þ
where c is a coefficient that introduces non-coaxiality between the total deviatoric strain rate D0 and the plastic strain rate Dp . A schematic diagram of the present corner-like model is shown in Fig. 1. The introduction of non-coaxiality between D0 and Dp is a necessity to reproduce the corner effect observed in crystal plasticity. Kuroda and Tvergaard (2001a) suggested that a range of c from 3 to 10 might be realistic through identification using Taylor model computations involving an abrupt strain path change. In a subsequent study on shear band simulations (Kuroda and Tvergaard, 2001c), it was found that a slightly smaller value, c = 2, was most suitable for reproducing the shear band development behavior predicted in crystal plasticity simulations. The present corner-like model
64
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
with
v ¼ n gp ;
ð15Þ
where n is a unit vector normal to the surface Sp of the plastic loading zone in the current configuration. Using the relations dV = JdV0, p J = det F, r ¼ r0 F1 , and ndSp ¼ n0 F1 JdS0 (known as Nanson’s formula) with subscript ‘‘0’’ referring to the reference configuration, Eq. (14) is rewritten as
Z p
V0
fdwðJR J re Þ þ r0 dw ðJF1 gp ÞgdV p0
¼
Z Sp0
dwn0 ðJF1 gp ÞdSp0 :
ð16Þ
For the incremental analysis of elastoplasticity problems, Eq. (16) is differentiated with respect to time, noting that _ ðF1 Þ ¼ F1 L and J ¼ J tr L, and again using JdV p0 ¼ dV p :
Fig. 1. Illustration of the corner-like effect adopted in the present plasticity model (Kuroda and Tvergaard, 2001a)).
Z hn o i R_ r_ e þ ðtr LÞðR re Þ dw þ fg_ p L gp þ ðtrLÞgp g rdw dV p Vp Z ¼ v^_ dwdSp Sp
ð17Þ
reduces to the J2-flow theory when Hpcrit ¼ 0 or c ! 1, while Simo’s original model is recovered for c = 0 (i.e. a = 1) and Hpcrit > 0 . Considering the relation p
r : D ¼ jr j/_ ¼ re 0
vb ¼ n fg_ p L gp þ ðtr LÞgp g:
rffiffiffi 2_ /; 3
ð11Þ
the equivalent plastic strain is naturally defined as
e_ p ¼
rffiffiffi 2_ /; 3
ep ¼
Z
t
e_ p dt;
ð12Þ
0
where t is time.1 In Eq. (1), the standard Jaumann stress rate has been used in the hypoelastic constitutive assumption. In the study of Kuroda and _ p with Tvergaard (2001b), the relation, W ¼ x þ Wp ¼ x þ /X Xp ¼ gðr n n rÞ, for the plastic spin for the present corner-like model was proposed and discussed through an extensive investigation of isotropic as well as anisotropic materials, followed by a comparison with experimental observations. According to this model, the plastic spin vanishes when n is coaxial to r0 , as in the present model. Thus, the spin of the substructure, x, is identified with the continuum spin, i.e. x = W. 2.2. Strategy for solving boundary value problems 2.2.1. Gradient-dependent yield condition The yield condition of Eq. (4) contains the second-order gradient term r gp , and thus this equation is viewed as a partial differential equation (PDE), unlike a standard yield condition in conventional plasticity theories. To get solutions to Eq. (4), its weak form is derived from an integral equation
Z Vp
fre þ r gp Rðep ÞgdwdV p ¼ 0;
ð13Þ
where V p is the volume of the plastic loading zone in the current (deformed) configuration and dw is an arbitrary scalar weighting function. Eq. (13) is integrated by parts to give an expression
Z Vp
1
fðR re Þdw þ gp rdwgdV p ¼
with
Z Sp
vdwdSp
ð14Þ
The t is not necessarily the ‘‘natural time’’ in the present theory since the constitutive model is time-independent.
ð18Þ
In Eq. (17), all integrations are performed in the current configuration, and thus this relation is suitable for an updated Lagrangian _ r_ e ; and g_ p ; computation scheme. The material-time derivatives, R; in Eq. (17) are calculated straightforwardly as
dR R_ ¼ p e_ p ¼ he_ p ; de
r_ e ¼
@ re @ re :r¼ :C:D @r @r
ð19Þ rffiffiffi 3 @ re : C : Np e_ p ; 2 @r
g_ p ¼ ðbrep Þ ¼ bre_ p gp L
ð20Þ ð21Þ
Eq. (14) with Eq. (15) and Eq. (17) with Eq. (18) enable a deduction of nonstandard boundary conditions, which are extra conditions compared with conventional elastoplasticity problems. The surface Sp consists of two parts: one is an external part Spext that belongs to the external boundary surface of the whole domain under consideration; the other is an internal part Spint that is an interface with the elastic (or unloaded) zone. Unlike the external boundary Spext , the internal boundary Spint moves inside the body during elastic–plastic deformation. Here, we define the plastic strain gradient flux as dF ¼ vdSp ¼ n gp dSp . When we consider a free external surface with no constraint on plastic flow, it may be acceptable to set zero plastic strain gradient flux, i.e. dF ¼ 0, _ ¼ 0. The latter is accomplished by or equivalently dF v^_ ¼ n fg_ p L gp þ ðtr LÞgp g ¼ 0 in an incremental analysis based on Eq. (17). The treatment of the conditions at the internal elastic– plastic boundary Spint may not be so simple since this boundary moves in the domain.2 We assume that there is no explicit constraint on e_ p at an elastic–plastic boundary. The corresponding condition for the plastic strain gradient flux should be considered only in an incremental manner. In the present paper, it is mimicked by _ ¼ 0), similarly to the external boundary, but ^_ ¼ 0 (i.e. dF setting v the condition dF ¼ 0 in the total form no linger holds at Spint . At a 2 If we consider a viscoplastic material in which plastic flow is potentially active everywhere, neither internal elastic–plastic boundary nor extra internal boundary condition appears in the problem (e.g., Borg et al., 2006).
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
boundary Sp interfacing with a hard or conceptually rigid material, zero plastic strain, e_ p ¼ 0, is a natural choice as a boundary condition. 2.2.2. Equilibrium condition The equilibrium condition in conventional solid mechanics remains the same in the present theory:
rr¼0
ð22Þ
with the body force effect neglected for brevity. The principle of virtual work expressed in the current configuration reads as a weak form of the PDE of Eq. (22),
Z
r : ðdu_ rÞdV ¼
V
Z
t du_ dS
ð23Þ
65
Isoparametric eight-node quadrilaterals with quadratic shape _ functions are used for the solution of the displacement rates u, while four-node quadrilaterals with linear shape functions are used for the solution of the equivalent plastic strain rate e_ p . The volume integrations are carried out by a 2 2 Gaussian integral. The structure of the finite element equations, Eq. (27), is similar to that of the equations explored by Mühlhaus and Aifantis (1991), Niordson and Hutchinson (2003) and Niordson and Redanz (2004). The relationship between the present formulation and those in the previous studies will be discussed in Section 4. 3. Numerical examples 3.1. Constrained simple shear with boundary layers
S
where t is the surface traction, V is the volume of the analysis domain or any part of it, and S is the surface of V (here, n is a unit normal to S). Differentiating the above equation through the same procedure as for Eq. (14), we obtain
A strip with a thickness H0 in the X2-direciton and an infinite length in the X1-direction is subjected to simple shear under a plane strain condition. Periodicity is assumed for all field quantities, taking an arbitrary periodic length L0 in the X1-direction. Consequently, values of the field quantities vary only in the thickness direction. The top and bottom surfaces are bounded by hard materials. The imposed boundary conditions are
Z
u_ 1 ¼ 0; u_ 2 ¼ 0; e_ p ¼ 0 on X 2 ¼ 0;
ð28Þ
~_ u_ 2 ¼ 0; e_ p ¼ 0 on X 2 ¼ H0 ; u_ 1 ¼ U;
ð29Þ
with
t ¼ n r;
ð24Þ
T
fr_ L r þ ðtrLÞrg : ðdu_ rÞdV ¼
V
Z
^t_ dudS _
ð25Þ
S
with
^t_ ¼ n fr_ L r þ ðtrLÞrg:
ð26Þ
Substituting the constitutive relations, Eqs. (1) and (3), gives the basis for finite element equations in an updated Lagrangian scheme. 2.2.3. Finite element analysis Based on Eqs. (25) and (17), finite element discretized equations can be derived in a standard manner, taking the displacement rate u_ and the equivalent plastic strain rate e_ p as independent variables, as
KðuuÞ KðpuÞ
( _ ) ( _ ) CðuÞ FðuÞ U ¼ þ : KðppÞ CðpÞ e_ p F_ ðpÞ KðupÞ
ð27Þ
_ is a vector type array of nodal displacement rates and fe_ p g Here, fUg is that of nodal equivalent plastic strain rates. The upper and lower groups of equations in Eq. (27) are derived from Eqs. (25) and (17), respectively, except for fCðuÞ g and fCðpÞ g, which are equilibrium correction terms calculated from the total form equations (23) and (14), respectively. Details of Eq. (27) are given in Appendix. Actual numerical computations are performed using a forward Euler time integration scheme in which increments of all quantities are taken to be DðÞ ¼ ðÞ _ Dt, and the time increment Dt is adaptively determined using the rmin method proposed by Yamada et al. (1968). If e_ p is computed to be negative despite the assumption of continuous plastic loading, elastic unloading is said to occur. In such a case, Dt is reduced to one-tenth of the value determined by the rmin method. As can be seen from Eq. (8), the present constitutive model is not linear in D. In the actual computations, the value of D used to construct Np is taken to be the value determined in the previous increment without iteration3 as in the study of Kuroda and Tvergaard (2001c). This is a good approximation for the small time increments used in the forward Euler method. The same kind of approximation has been employed in calculations based on the J2 corner theory (Hutchinson and Tvergaard, 1981; Tvergaard et al., 1981). 3 In Kuroda and Tvergaard (2001a), an iterative method was used for the same constitutive model.
~_ is the prescribed displacement rate in the X1-direction on where U the top surface. The macroscopic shear strain is defined by ~ 0. C ¼ U=H The strain hardening function R(ep) is chosen to be
ep n Rðep Þ ¼ r0 1 þ ;
e0
ð30Þ
where r0 is the initial yield stress, e0 is the offset strain, and n is the hardening exponent. The material parameter values used are E/ r0 = 500, m = 0.3, e0 = 0.002, and n = 0.1, where E is the Young’s modulus and m is the Poisson’s ratio. These values are the same as those used in Kuroda and Tvergaard (2001b).4 The length scale coefficient b is defined as
b¼l
2
r0 ;
ð31Þ
where l is a length-scale parameter. In this problem, plastic yielding occurs simultaneously at all material points, and no elastic unloading occurs during the subsequent increase in the shear deformation. Therefore, no internal boundary between the elastic and elastic–plastic zones emerges. Fig. 2 shows curves of stresses versus macroscopic shear strain for several typical cases. Results for uniform simple shear without a strain gradient, which are obtained by employing an extra ^_ ¼ 0 instead of e_ p ¼ 0 in Eqs. (28) and (29), boundary condition v are also shown for reference. In this case, all size effects are excluded, i.e. the results are completely identical to those for the conventional size-independent theory with b = 0. The parameter c, which governs the non-coaxiality between D0 and Dp , is taken as c = 2–3, and Hpcrit is chosen to be 20°. A remarkable feature of the present corner-like plasticity model is the generation of considerable amounts of normal stresses. In the case of the J2-flow theory (i.e. Hpcrit = 0°), the normal stresses are close to zero.5 The amounts of normal stresses predicted by the present model are consistent 4 In Kuroda and Tvergaard (2001b), a viscoplastic version of the present model was used with a very small value of the rate sensitivity parameter (i.e. m = 0.002) in a power law relation, which exhibited nearly rate-independent behavior. 5 These very small normal stresses are due to the use of the hypoelastic assumption (i.e. Eq. (1)).
66
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
2 l /H
1.5
l /H0 = 0.3
ij /
Stress
0
12
1 0.5
The deformed mesh at C = 1 in the case of l/H0 = 0.3 is illustrated in Fig. 4(b), which shows the effect of the extra boundary condition e_ p ¼ 0 on the overall deformation mode of the strip.
.5
=0 0
3.2. Shear band development in plane strain tension
Conventional
c =3,
p
c =2,
p
J2F (
crit =
20
crit = 20 p crit = 0 )
0 22
-0.5
l /H
0=
0.5
l /H0 = 0.3
A rectangular specimen under plane strain tension with X2 being the tensile axis is considered. The initial dimensions of the specimen are 2H0 along the X1-direction and 2L0 along the X2direction, and the aspect ratio L0/H0 is taken to be 3. Considering the geometrical symmetry, only one quadrant of the specimen is analyzed. The boundary conditions for the quadrant are
^_ ¼ 0 on X 1 ¼ 0; u_ 1 ¼ 0; ^t_ 2 ¼ 0; v
ð32Þ
^_ ¼ 0 on X 2 ¼ 0; u_ 2 ¼ 0; ^t_ 1 ¼ 0; v
ð33Þ
^t_ 1 ¼ 0; ^t_ 2 ¼ 0; v ^_ ¼ 0 on X 1 ¼ H0 þ DH0 ;
ð34Þ
_ ^t_ 1 ¼ 0; v ^_ ¼ 0 on X 2 ¼ L0 ; u_ 2 ¼ U;
ð35Þ
-1 0
0.5
1 1.5 2 Macroscopic shear strain
2.5
Fig. 2. Curves of stresses vs macroscopic shear strain for constrained simple shear problem of gradient-dependent corner-like plasticity solids. Stress responses in the absence of a gradient effect (designated as ‘conventional’) are also depicted for comparison. ‘‘J2F’’ stands for the J2-flow theory.
with the amounts of axial stresses observed in fixed-end torsion tests as discussed by Kuroda and Tvergaard (2001b). Next, we examine size effects in constrained simple shear. For the computations with finite length scales, a finite element mesh of 1 (in the X1-direction) 30 (in the X2-direction) quadrilaterals is used, which leads to the convergence of solutions on the plotted graph scale. It can be seen in Fig. 2 that large amounts of additional hardening are produced with increasing length scale (l/H0 = 0.3 and 0.5). Note that the stress components r22 and r12 must be uniform in the thickness direction due to the equilibrium requirement from Eq. (22). The J2-flow theory (i.e. Hpcrit = 0°) produces only insignificant normal stresses even in this size dependent setting. The stress distributions at C = 1 for the corner-like plasticity solids (Hpcrit = 20° and c = 2) with the length scales of l/H0 = 0.3 and 0.5 are plotted in Fig. 3. Although the distributions of the normal stresses, r11 and 11 and r 33 , r33, are strongly nonuniform, their volume averages, r 11 ¼ r22 and r 33 ¼ 0, respectively, simisatisfy the conditions r larly to the observation in the case of uniform simple shear. The stress distributions observed here are very similar to those for the finite strain gradient crystal plasticity model reported in Kuroda and Tvergaard (2008). The distributions of ep are plotted in Fig. 4(a). Those for l/H0 = 0.3 and 0.5 are close to each other.
1
Position X2 / H0
l /H0 = 0.3 l /H0 = 0.5
0.5
12 11
0 -1
-0.5
0
0.5 Stress
DH0 ¼ H0 ½n1 cosðpX 2 =L0 Þ þ n2 cosðmw pX 2 =L0 Þ
ð36Þ
with amplitudes of the imperfection waves, n1 and n2, and the inte^_ ¼ 0 is also valid on ger wave number mw. Note that the condition v the symmetry planes as specified in Eqs. (32) and (33). The macro 0 . The nondimensionalscopic nominal strain is defined by eN ¼ U=L ized tensile load T is defined by
T¼
P ; H 0 r0
ð37Þ
where P is the sum of the nodal reaction forces at X2 = L0 for the one-quadrant model with a unit depth in the X3-direction. Eq. (30) is employed for the strain hardening function R(ep). The material parameter values used are E/r0 = 333, m = 0.3, e0 = 0.0015, and n = 0.1, which are the same as those used in the study of Kuroda and Tvergaard (2001c). The controlling parameters for the corner-like effect are taken as Hpcrit X= 20° and c = 2 unless stated otherwise. These values were used in the study of Kuroda and Tvergaard (2001c) in the conventional size-independent case. For the length-scale coefficient b, Eq. (31) is also employed in this problem. The extra boundary conditions at internal boundaries between the elastic zone and the elastic–plastic zone enter this problem since elastic unloading occurs in a broad area at the onset of diffuse necking and the unloading region gradually develops with increasing overall elongation of the specimen. In the present simulations, no special numerical treatment is applied at the internal elastic– n o plastic boundaries, i.e. no input is given to F_ ðpÞ in Eq. (27) for the corresponding elements. This approximately satisfies the con^_ ¼ 0 on the internal boundaries. As a result, all the comdition of v n o ponents of F_ ðpÞ remain zero since the external boundary
33
22
_ is the prescribed rate of end-displacement and DH0 reprewhere U sents the initial geometrical imperfection of the form (Tvergaard et al. 1981; Kuroda and Tvergaard, 2001c)
1 ij /
1.5
2
0
Fig. 3. Distributions of stresses at a macroscopic shear strain of C = 1.0 in the constrained simple shear problem of gradient-dependent corner-like plasticity solids (Hpcrit = 20° and c = 2) with length scales of l/H0 = 0.3 and 0.5.
conditions, as well as the symmetry conditions, are all given by v^_ ¼ 0 as in Eqs. (32)–(35). Niordson and Redanz (2004) reported that in their plane strain sheet necking analysis using the Fleck– Hutchinson’s higher-order theory, the two conditions n re_ p ¼ 0 (in their formulation) and e_ p ¼ 0 made no significant difference. First, we examine the effect of the mesh density on strain localization predictions. Fig. 5 displays curves of computed tensile load versus macroscopic nominal strain for the corner-like plasticity solids (Hpcrit = 20° and c = 2) with and without the length-scale
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a) 1
67
(b) l /H0 = 0.3
Position
X2 / H0
l /H0 = 0.5
= 0.5
0.5
1.0
1.5
2.0
2.5
0 0
0.5 1 1.5 Equivalent plastic strain
2
2.5
p
Fig. 4. (a) Distributions of equivalent plastic strain ep at various stages of macroscopic shear strain C in the constrained simple shear problem of gradient-dependent cornerlike plasticity solids (Hpcrit = 20° and c = 2) with length scales of l/H0 = 0.3 and 0.5; (b) deformed mesh with contours of ep for gradient-dependent corner-like plasticity model (Hpcrit = 20° and c = 2) with l/H0 = 0.3 at C = 1.
Nondimensionalized tensile load T
2 l / H0 = 0.04
1.5 Conventional
1 30 x 90 quadrilaterals 12 x 36 quadrilaterals
0.5
0 0
0.1 0.2 Macroscopic nominal strain
0.3 N
Fig. 5. Computed curves of tensile load vs macroscopic nominal strain for plane strain tension problem employing corner-like-effect parameters of Hpcrit = 20and c = 2 with and without the length-scale effect. Two different mesh densities (12 36 and 30 90) are considered. In all the computations, the initial geometric imperfection is specified by n1 = 0.005 and n2 = 0.
effect, employing two finite element discretizations: a mesh of 12 36 quadrilaterals and a mesh of 30 90 quadrilaterals6 with a pure cosine mode imperfection (n1 = 0.005 and n2 = 0). Deformed meshes and contours of equivalent plastic strain at stages where strain localizations in the form of a shear band have been well developed are depicted in Fig. 6. In the conventional size-independent case (l = 0), significant mesh dependence appears as expected (see Fig. 6 (a) and (b)). The shear band that appears in the mesh of 30 90 quadrilaterals is much more severe than that appearing in the mesh of 12 36 quadrilaterals. The tensile load T for the fine mesh begins to decrease slightly earlier than that for the coarse mesh as shown in Fig. 5. In the finite length-scale case of l/ H0 = 0.04, the widths of the shear band for both the finite element discretizations are almost the same (Figs. 6(c) and (d)). The corre-
6 In both the finite element discretizations, the initial size of each element in the X2-direction is given by a geometric progression with a common ratio of 1.02, so that elements are more densely assigned in the potential neck region than in the anticipated unloading region.
sponding curves of T versus eN for these two computations are also almost identical (Fig. 5). Fig. 7 displays curves of computed tensile load T versus macroscopic nominal strain eN for a somewhat complex imperfection with n1 = 0.005, n2 = 0.001, and mw = 8. The computations are carried out for several different length scales using the mesh of 30 90 quadrilaterals. Deformed meshes and contours of equivalent plastic strain at stages where strain localizations have been well developed are depicted in Fig. 8. In the conventional sizeindependent case (l = 0), the most severely strained point is on the free surface at X2 = 0 and a localized band is formed from the free surface (Fig. 8(a)). With an increase in the length scale relative to the specimen size, the patterns of the flow localization zone drastically change. For l/H0 = 0.01, a severely strained spot also appears at the center of the specimen, and the localized band from the free surface ceases to develop. The deformation mode for l/ H0 = 0.04 (Fig. 8(c)) becomes almost identical to that in the case of the pure cosine imperfection (Fig. 6(d)). Thus, the gradient term introduced into the present corner-like plasticity model governs not only the ‘‘width’’ of shear bands, but also the patterns of strain localization. For l/H0 = 0.1, the shear band completely disappears and only the diffuse neck-type deformation is observed. Finally, the effect of c is illustrated in Fig. 9. Here, c is increased to 3 and the other settings are unchanged from the previous computations whose results are shown in Figs. 7 and 8. Compared with the results shown in Fig. 8 for c = 2, significantly different shear band patterns are observed for c = 3. The most severely strained point appears at the center of the specimen even in the size-independent case (Fig. 9(a)). For larger length scales, l/H0 = 0.04 and 0.1, the results for c = 3 are rather similar to those for c = 2 with the same length scales (Fig. 8 (c) and (d)).
4. Discussion 4.1. On modeling of corner plasticity The key to the successful extension from the conventional sizeindependent theory to the strain gradient-enhanced version in the present study is that the constitutive relation for Dp is given by a form similar to the conventional flow rule, i.e. the product of the qffiffi
direction tensor Np and the plastic multiplier /_ ¼ 32e_ p as is seen in Eq. (3) with Eq. (12). Owing to this simple constitutive structure,
68
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a)
(b)
(c)
(d)
Fig. 6. Results of mesh dependence tests with and without length scale effect (the corner-like-effect parameters are taken to be Hpcrit = 20and c = 2) for plane strain tension problem showing deformed meshes and contours of the equivalent plastic strain ep for (a) a size-independent (conventional) model with a mesh of 12 36 quadrilaterals at eN = 0.240; (b) a conventional model with a mesh of 30 90 quadrilaterals at eN = 0.240; (c) the gradient-dependent model of l/H0 = 0.04 with a mesh of 12 36 quadrilaterals at eN = 0.270; (d) the gradient-dependent model of l/H0 = 0.04 with a mesh of 30 90 quadrilaterals at eN = 0.270. In all the computations, an initial geometric imperfection specified by n1 = 0.005 and n2 = 0 (a pure cosine mode) is introduced.
Nondimensionalized tensile load T
2
1.5
Conventional
1
l / H0 = 0.01 l / H0 = 0.04
0.5
l / H0 = 0.1
0 0
0.1 0.2 Macroscopic nominal strain
0.3 N
Fig. 7. Computed curves of tensile load vs macroscopic nominal strain for plane strain tension problem, employing corner-like-effect parameters of Hpcrit = 20° and c = 2 with different length scales. The initial geometric imperfection is specified by n1 = 0.005, n2 = 0.001, and mw = 8, and a mesh of 30 90 quadrilaterals is used in all the computations.
the extension to the strain gradient version could be carried out in a similar manner to the case of the strain gradient theory based on the classical J2 flow theory. The other corner theories (i.e. Christoffersen and Hutchinson (1979) and Gotoh (1985)) assumes
the dependence of the plastic strain rate direction on the stress rate, which are not given in a simple form like Eq. (3) with Eq. (12). Although the actual formulations are not pursued in this paper, the extension of these corner theories to corresponding higher-order strain gradient versions would not be so straightforward unlike the case of the present theory. Under conventional size-independent postulations, these corner theories as well as the present theory have basically similar predictability for plastic flow localization phenomena. In fact, the present corner-like theory provides a prediction of shear band formation very similar to that by the Christoffersen and Hutchinson’s theory, depending on the choice of the parameter values, as illustrated in Kuroda and Tvergaard (2001c). Thus, in the context of enhancement of the theories with introduction of the size effect, the present corner-like higher-order gradient plasticity theory may have some worth as an alternative to other corner-type theories whose extension to strain gradient versions is not straightforward. 4.2. On predictability of plastic flow localization As shown in Tvergaard et al. (1981), the initial geometrical imperfection mode plays an important role in determining the pattern of shear band formation. In the present study, the same finding can be seen on comparison between the results shown in Figs. 6(b) and 8(a). The difference in these two computations was only the modes of the initial geometrical imperfection, and the both did not involve the size effect. It is clearly seen in Fig. 8 that
69
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a)
(b)
(c)
(d)
Fig. 8. Effects of the length scale on patterns of strain localization (the corner-like-effect parameters are taken to be Hpcrit = 20° and c = 2) in plane strain tension problem, introducing an initial geometric imperfection specified by n1 = 0.005, n2 = 0.001, and mw = 8 showing deformed meshes and contours of the equivalent plastic strain ep for (a) a size-independent (conventional) model at eN = 0.230; (b) l/H0 = 0.01 at eN = 0.250; (c) l/H0 = 0.04 at eN = 0.274; (d) l/H0 = 0.1 at eN = 0.300. A mesh of 30 90 quadrilaterals is used in all the computations.
the introduction of the gradient effect varies not only the width of shear band but also the fundamental patterns of flow localization. This is a new point found in the present study with respect to the shear band behavior. 4.3. On formulations and computations of strain gradient plasticity In the present paper, the terms ‘higher-order stress’ and higherorder traction’ have not appeared so far at all, which were among the central items in several previous studies on strain gradient plasticity. In this section, the connection between the present treatment of the higher-order strain gradient effects and those considered by other researchers is discussed. In the present study, the length-scale coefficient b has been taken as a constant for simplicity. In this case, the term r gp in Eq. (4) can be written as br2 ep . The introduction of the vector quantity gp would be essential for mathematical transformation to different expressions of the theory as shown below. In Eq. (14), formally replacing dw by a virtual equivalent plastic strain rate de_ p , we have
Z Vp
fðR re Þde_ p þ gp rde_ p gdV p ¼
Z Sp
vde_ p dSp :
ð38Þ
Now, the following constraint for the virtual quantities is considered:
ðdu_ rÞsym ¼ dD ¼ dDe þ dDp ¼ dDe þ
rffiffiffi 3 p p de_ N : 2
ð39Þ
Upon noting that de_ p is only taken into account within the plastic zone and considering Eqs. (11) and (12), Eq. (38) with Eq. (39) is rewritten as
Z
ðr : dDp þ Rde_ p þ gp rde_ p dV ¼
V
Z
vde_ p dS:
ð40Þ
S
Adding Eq. (40) to the virtual work principle (Eq. (23)) gives
Z V
fðr : dDe þ Rde_ p þ gp rde_ p gdV ¼
Z S
_ þ t dudS
Z
vde_ p dS: ð41Þ S
The relation (41) is the same as the modified virtual work statement in Fleck and Hutchinson’s strain gradient plasticity theory (Fleck and Hutchinson, 2001). Eq. (41) was introduced as a major assumption at the starting point of their theory. In Eq. (41), gp (in their paper a quantity equivalent to gp is denoted by s) was introduced as a work conjugate to the plastic strain gradient rep . Fleck and Hutchinson referred to it as a higher-order stress vector quantity. Similarly, v was introduced as a work conjugate to the equivalent plastic strain ep on the surface, called a higher-order traction quantity. In the study of Fleck and Hutchinson (2001), the major premise of Eq. (41) yields the conventional equilibrium condition (Eq. (22)) and the gradient dependent yield condition (Eq. (4)), as opposed to the present formulation of gradient plasticity, where all the PDEs are given initially based on physics and physically based insights. It is obvious that the present theory can be equivalent to an approach based on a relation in the form of Eq. (41), as long as the same constitutive relations are postulated for Dp and gp.
70
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
(a)
(b)
(c)
(d)
Fig. 9. Deformed meshes and contours of equivalent plastic strain ep for plane strain tension problem when c = 3 with different length scales (Hpcrit remains at 20°), introducing an initial geometric imperfection specified by n1 = 0.005, n2 = 0.001, and mw = 8: (a) size-independent (conventional) model at 0.250; (b) l/H0 = 0.01 at 0.250; (c) l/ H0 = 0.04 at 0.274; (d) l/H0 = 0.1 at 0.300. A mesh of 30 90 quadrilaterals is used in all the computations.
Gurtin and Anand (2009) addressed a thermodynamically consistent derivation of the constitutive relations for the higher-order quantities in the Fleck and Hutchinson’s framework (Eq. (41)): i.e. in their argument, if we assume that a defect energy, wp , of the form wp ¼ 12 bjrep j2 augments the free energy (per volume), the Clausius–Duhem inequality leads to the constitutive equation gp ¼ brep that is identical to Eq. (6). According to this view, we can claim that gp is recoverable (or energetic in terminology by Gurtin and Anand (2009)); thus, plastic dissipation must occur via Re_ p ðP 0Þ. This view is consistent with the argument presented in Section 2.1, which is based on the dislocation theory. An incremental version of Eq. (41) is derived by adding Eq. (17) to Eq. (25) with de_ p that replaces dw
Z n o ^_ r ^_ T : ðdu_ rÞ þ ðR ^_ p rde_ p dV ^_ e Þde_ p þ g P V Z Z _ _ ^_ de_ p dS þ v ¼ ^t dudS S
ð42Þ
S
with definitions
^_ ¼ r_ L r þ ðtrLÞr; P
ð43Þ
^_ ¼ R_ þ ðtr LÞR; R
ð44Þ
r^_ e ¼ r_ e þ ðtr LÞre ;
ð45Þ
^_ p ¼ g_ p L gp þ ðtr LÞgp : g
ð46Þ
_ are viewed as the nominal These quantities with superposed ‘^’ rates of the corresponding quantities at the current configuration. Eq. (42) is identified with the incremental virtual work principle of Niordson and Redanz (2004), which was derived on the basis of the Fleck and Hutchinson’s theory (a small strain theory). The treatment of the quantity gp (in their study it is a higher-order stress, s) is, however, somewhat different. In their model, a constitutive rela^_ p in an incremental form, which was tion regarding gp was given to g deduced from the Fleck and Hutchinson’s theory, while in the present theory, the relation for gp is defined in a total form (Eq. (6)) as in the original Aifantis model. It is emphasized that Eq. (42) can give exactly the same finite element equations (Eq. (27)) as long as all the constitutive relations are the same. The finite element formulation presented earlier in Niordson and Hutchinson (2003) was a small strain version of that in Niordson and Redanz (2004). The basic structure of Eq. (27) is also the same as the formulation of Mühlhaus and Aifantis (1991), which was the earliest proposed finite element procedure for the higher-order gradient plasticity, although it included a small strain assumption. Mühlhaus and Aifantis (1991) considered an incremental potential functional accounting for the gradient effects. Applying the stationary condition to the potential functional gives an equation very similar to the incremental form of Eq. (41), except for the higherorder traction term on the right-hand side. Mühlhaus and
71
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
Aifantis (1991) did not explicitly define any higher-order quantity, and they only considered the situations v_ = 0 or e_ p = 0 (under a small strain setting), assuming that the product of these quantities is always zero. Kuroda and Tvergaard (2010) discussed a different mathematical treatment of higher-order gradient plasticity theories, in which the increment of the plastic strain gradient is taken as an additional independent variable instead of the increment of the plastic strain, without introducing the concept of higher-order stress. This idea can be applied to the present theory in principle. However, in the case of large strain problems, the treatment of the extra bound^_ ¼ n fg_ p L gp þ ðtr LÞgp g ¼ 0, is not so simple ary condition, v when g_ p is chosen as an independent variable. To avoid this complexity, the solution strategy shown in section 2.2 has been adopted in the present computations. On the computational side, all the numerical simulations performed here have adopted quadratic shape functions for the displacement rates and linear shape functions for the equivalent plastic strain rates. Empirically, this combination gives high-quality numerical solutions with fast convergence with reasonable mesh densities and with no stress oscillation phenomenon as found in Fredriksson et al. (2009) and Kuroda (2011) for some higher-order gradient plasticity computations. A more extensive and systematic investigation regarding the quality of finite element solutions for the various formulations of strain gradient plasticity is needed towards further expanding the use of strain gradient plasticity. This is left for a subsequent study.
_ ¼ fu_ 1 ; u_ 2 gT at a material point inside the element, we write fug _ _ ¼ ½NfUg, the finite element discretized relations as fug _ (with fDg ¼ fD11 ; D22 ; 2D12 gT ), fLg ¼ ½BL fUg, _ fDg ¼ ½BD fUg and _ where [BD], [BL] and [divB] contain spatial derivatr L ¼ ½divBfUg, tives of the shape functions in appropriate arrays. Similarly, a shape function matrix [Ne] relates the nodal equivalent plastic strains fe_ p g to the equivalent plastic strain at a material point inside the element as e_ p ¼ ½N e fe_ p g, and its spatial gradients, fre_ p g ¼ f@ e_ p =@x1 ; @ e_ p =@x2 g, are computed as fre_ p g ¼ ½Be fe_ p g with the matrix [Be] containing spatial derivatives of the components of [Ne]. The orders of the shape functions in [N] for the displacement rates and in [Ne] for the equivalent plastic strain rates are not necessarily the same. Upon the above preparation, the finite element equations in Eq. (27) are derived from Eqs. (25) and (17). Detailed expressions for each component in Eq. (27) are shown as follows. The incremental virtual work principle (Eq. (25)) gives
KðuuÞ ¼
8 9 Z < = ½BD T ½C½BD þ ½BL T ½C r ½BL dV; |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} ; V :
ð49Þ
finite strain effect
where [C] is elastic moduli (Eq. (2)) in matrix form,
½KðupÞ ¼
rffiffiffi Z 3 2l½BD T Np ½Ne dV p 2 Vp
with the plastic strain T N p ¼ N p11 ; N p22 ; Np12 , and
rate
ð50Þ direction
of
n o Z _ F_ ðuÞ ¼ ½NT f^tgdS:
5. Concluding remarks
the
form
ð51Þ
S
The corner-like plasticity model that was originally proposed by Simo (1987) and later modified by Kuroda and Tvergaard (2001a) has been extended to include the size effect resulting from the plastic strain gradient. A method of solving boundary value problems at finite strains has also been given, and the efficiency of the model in problems including higher-order boundary conditions and plastic flow localization has been demonstrated through finite element simulations of problems involving constrained simple shear and plane strain tension. Appendix
Eq. (17) gives
KðpuÞ ¼
8 9 Z < rffiffiffi = 3 T T T T p 2l½N e fng ½BD þ ðR re Þ½N e ½divB þ ½Be ½g ½BD dV p |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}; 2 Vp : finite strain effect
ð52Þ T
with fng ¼ fn11 ; n22 ; n12 g ,
"
p
½g ¼
g p1 g p2
ð47Þ
KðppÞ ¼
Z Vp
g p2
g p1
# ;
n o ðh þ 3lÞ½Ne T ½Ne þ b½Be T ½Be dV p ;
ð53Þ
ð54Þ
Sp
In Eqs. (49) and (52), the terms with the note ‘‘finite strain effect’’ are peculiar to the present finite strain formulation, and would not appear in a small strain formulation. The total form of the virtual work principle (Eq. (23)) yield the equilibrium correction vector
Z Z fCðuÞ g ¼ a ½BD frgdV þ ½NT ftgdS V
ð55Þ
S
with {r} = {r11, r22, r12}T and {t} = {t1, t2}T, and the weak form of the yield function (Eq. (14)) gives
Z CðpÞ ¼ a
Vp
h
Z i ½Ne T ðR re Þ þ ½Be T fg p g dV p þ ½N e T vdSp ; Sp
ð56Þ
ð48Þ
Introducing a matrix array of shape functions, [N], which relates the _ to the displacement rates displacement rates fUg
nodal
g p2
n o Z ^_ dSp : F_ ðpÞ ¼ ½Ne T v
or in a symbolic expression
r þ ½C r fLg:
and
A practical form of Eq. (27) is presented for plane strain condi^_ def tions. The terms P ¼ r_ L r þ ðtrLÞr in Eq. (25) are viewed as a nominal stress rate in the current configuration. With the introduc ^_ is tion of the Jaumann rate r of the Cauchy stress (Eq. (1)), P expressed in a vector–matrix-type expression that may be useful in coding a computer program (e.g., Yamada and Sasaki, 1995): 8 9 > > r11 > > 8 9 > > > > > 2 38 9 _ > > > > > ^ > > > P11 > 0 r12 0 r11 > L11 > > > > > > > > > > > > > > > > > > 6 > > > > 7> > > > > > > r _ > > > > > 6 7 ^ 12 = < 12 = 6 0 1 ðr11 r22 Þ 1 ðr11 þ r22 Þ r12 7< L21 >
> > _ > > > > > > > 6 r21 1 ðr11 þ r22 Þ 12 ðr22 r11 Þ 0 7> L12 > ^ > > > > > > 2 P r 21 > > > > 21 > > > > > 4 5> > > > > > > > : > > ; > >_ > > > > > : > ; > r 0 r12 0 L22 22 > > ^ > > > P22 > > : r22 > ;
n o ^_ ¼ P
g p1
T where fg g ¼ g p1 ; g p2 . The value of gp is evaluated using the evolution equation (21) or the gradients of the total nodal equivalent plastic strain {ep} in the current configuration. The introduction of the coefficient a in Eqs. (55) and (56) is peculiar to the rmin method p
72
M. Kuroda / International Journal of Solids and Structures 58 (2015) 62–72
because Dt is adaptively determined after solving the system of equation (27). If Dt were determined to equal 1/a, the equilibrium and yield condition corrections would be fully done. But, this cannot be fulfilled in the incremental method adopted here. In the actual computations, a is approximately evaluated from the previous time increment size. In the present applications shown in Section 3, however, values of the correction terms were always small due to the use of small increment size; thus difference between results with and without these terms would not be visible on the scales plotted. During the incremental computations, the integrations for [K(up)], [K(pu)], [K(pp)] and {C(p)} are done only at Gaussian points in the plastic loading state. For [K(pp)], a matrix of 108E [1] with [1] being an identity matrix is set when the point is elastic in order to avoid an occurrence of numerical instability. The computer code used is a fully in-house program written in FORTRAN. References Aifantis, E.C., 1984. On the microstructural origin of certain inelastic models. Trans. ASME J. Eng. Mater. Technol. 106, 326–330. Aifantis, E.C., 1987. The physics of plastic deformation. Int. J. Plast. 3, 211–247. Ashby, M.F., 1970. The deformation of plastically non-homogeneous alloys. Philos. Mag. 21, 399–424. Borg, U., 2007. Strain gradient crystal plasticity effects on flow localization. Int. J. Plast. 23, 1400–1416. Borg, U., Niordson, C.F., Fleck, N.A., Tvergaard, V., 2006. A viscoplastic strain gradient analysis of materials with voids or inclusions. Int. J. Solids Struct. 43, 4906– 4916. Christoffersen, J., Hutchinson, J.W., 1979. A class of phenomenological corner theories of plasticity. J. Mech. Phys. Solids 27, 465–487. Evers, L.P., Brekelmans, W.A.M., Geers, M.G.D., 2004. Non-local crystal plasticity model with intrinsic SSD and GND effects. J. Mech. Phys. Solids 52, 2379–2401. Fleck, N.A., Hutchinson, J.W., 2001. A reformulation of strain gradient plasticity. J. Mech. Phys. Solids 49, 2245–2271. Fleck, N.A., Muller, G.M., Ashby, M.F., Hutchinson, J.W., 1994. Strain gradient plasticity: theory and experiment. Acta Matetallurgica et Materialia 42, 475– 487. Fredriksson, P., Gudmundson, P., Mikkelsen, L.P., 2009. Finite element implementation and numerical issues of strain gradient plasticity with application to metal matrix composites. Int. J. Solids Struct. 46, 3977–3987. Gotoh, M., 1985. A simple plastic constitutive equation with vertex effect. Eng. Fract. Mech. 21, 673–684. Groma, I., Csikor, F.F., Zaiser, M., 2003. Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics. Acta Mater. 51, 1271–1281. Gudmundson, P., 2004. Unified treatment of strain gradient plasticity. J. Mech. Phys. Solids 52, 1379–1406.
Gurtin, M.E., 2002. A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations. J. Mech. Phys. Solids 50, 5–32. Gurtin, M.E., Anand, L., 2009. Thermodynamics applied to gradient theories involving the accumulated plastic strain: the theories of Aifantis and Fleck and Hutchinson and their generalization. J. Mech. Phys. Solids 57, 405–421. Hughes, T.R., Shakib, F., 1986. Pseudo-corner theory: a simple enhancement of J2flow theory for applications involving non-proportional loading. Eng. Comput. 3, 116–120. Hutchinson, J.W., Tvergaard, V., 1981. Shear band formation in plane strain. Int. J. Solids Struct. 17, 451–470. Hutchinson, J.W., 2012. Generalizing J2 flow theory: fundamental issues in strain gradient plasticity. Acta Mech. Sin. 28, 1078–1086. Kuroda, M., 2011. On large-strain finite element solutions of higher-order gradient crystal plasticity. Int. J. Solids Struct. 48, 3382–3394. Kuroda, M., Tvergaard, V., 1999. Use of abrupt strain path change for determining subsequent yield surface: illustrations of basic idea. Acta Mater. 47, 3879–3890. Kuroda, M., Tvergaard, V., 2001a. A phenomenological plasticity model with nonnormality effects representing observations in crystal plasticity. J. Mech. Phys. Solids 49, 1239–1263. Kuroda, M., Tvergaard, V., 2001b. Plastic spin associated with a non-normality theory of plasticity. Eur. J. Mech. A/Solids 20, 893–905. Kuroda, M., Tvergaard, V., 2001c. Shear band development predicted by a nonnormality theory of plasticity and comparison to crystal plasticity predictions. Int. J. Solids Struct. 38, 8945–8960. Kuroda, M., Tvergaard, V., 2006. Studies of scale dependent crystal viscoplasticity models. J. Mech. Phys. Solids 54, 1789–1810. Kuroda, M., Tvergaard, V., 2008. A finite deformation theory of higher-order gradient crystal plasticity. J. Mech. Phys. Solids 56, 2573–2584. Kuroda, M., Tvergaard, V., 2010. An alternative treatment of phenomenological higher-order strain-gradient plasticity theory. Int. J. Plast. 26, 507–515. Kuwabara, T., Kuroda, M., Tvergaard, V., Nomura, K., 2000. Use of abrupt strain path change for determining subsequent yield surface: experimental study with metal sheets. Acta Mater. 48, 2071–2079. Mühlhaus, H.-B., Aifantis, E.C., 1991. A variational pronciple for gradient plasticity. Int. J. Solids Struct. 28, 845–857. Niordson, C.F., Hutchinson, J.W., 2003. Non-uniform plastic deformation of micron scale object. Int. J. Numer. Methods Eng. 56, 961–975. Niordson, C.F., Redanz, P., 2004. Size-effects in plane strain sheet-necking. J. Mech. Phys. Solids 52, 2431–2454. Simo, J.C., 1987. A J2-flow theory exhibiting a corner-like effect and suitable for large-scale computation. Comput. Methods Appl. Mech. Eng. 62, 169–194. Stölken, J.S., Evans, A.G., 1998. A microbend test method for measuring the plasticity length scale. Acta Mater. 46, 5109–5115. Støren, S., Rice, J.R., 1975. Localized necking in thin sheets. J. Mech. Phys. Solids 23, 421–441. Tvergaard, V., Needleman, A., Lo, K.K., 1981. Flow localization in the plane strain tensile test. J. Mech. Phys. Solids 29, 115–142. Yamada, Y., Yoshimura, N., Sakurai, T., 1968. Plastic stress-strain matrix and its application for the solution of elastic–plastic problems by the finite element method. Int. J. Mech. Sci. 10, 343–354. Yamada, Y., Sasaki, M., 1995. Elastic–plastic large deformation analysis problem and lamina compression test. Int. J. Mech. Sci. 37, 691–707.