Instability analysis of shear bands using the instantaneous growth-rate method

International Journal of Impact Engineering xxx (2015) 1e13

Contents lists available at ScienceDirect

International Journal of Impact Engineering journal homepage: www.elsevier.com/locate/ijimpeng

Instability analysis of shear bands using the instantaneous growth-rate method Miguel Arriaga, Colin McAuliffe, Haim Waisman* Department of Civil Engineering and Engineering Mechanics, Columbia University, 610 Seeley, W. Mudd Building, 500 West 120th Street, Mail Code 4709, New York, NY 10027, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Available online xxx

Shear banding, as an unstable process of localization, is a common precursor to fracture in materials under high strain rate loadings, making the detection of the instability point after which localization will occur of significant importance. Stability analysis based on the perturbation method or the acoustic tensor have been employed. However, these methodologies are limited to certain class of problems and are difficult to generalize. In this work we propose an alternative for identifying the instability point by employing the concept of generalized stability analysis. In this framework, a stability measure is obtained by computing the instantaneous growth rate of the vector tangent to the solution. Such an approach is more appropriate for non-orthogonal problems and is easier to generalize to difficult dynamic fracture problems. Under conditions where the local instability triggers the non-homogeneous solution growth, i.e. problems that are homogeneous until the appearance of local instability, the generalized stability analysis and the modal stability analysis will closely match. Therefore, the non-homogeneous growth can be approximated by the Rayleigh Quotient of the vector tangent to the solution, which is easier to compute. We show that for a particular class of problems that respect the aforementioned conditions, in 1D and 2D examples, both quantities successfully find the instability point predicted analytically and validated experimentally in past literature results. This methodology is general and can be applied to a wide array of dynamic fracture problems, for which instability that leads to localization is important. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Shearband Instantaneous growth-rate Instability Failure Rayleigh Quotient

1. Introduction Shear bands are material instabilities associated with highly localized intense plastic deformation zones. In materials such as metals they typically form under high strain rate loadings and are a strong precursor to fracture [57]. These localized bands have also been observed in materials with non-associative flow laws like porous or granular materials, polymers, rocks and soils [17,27]. The flow law of many structural materials show strain rate dependency. In a material like steel a high strain rate loading hardens the material. On the other hand, the heat produced by inelastic deformation generates thermal softening that will lower the material's flow stress potentially causing instability. These two effects added to the strain hardening of the material lead to an intricate

* Corresponding author. E-mail address: [email protected] (H. Waisman).

equilibrium that is highly sensitive to changes in material parameters and loading conditions. Given the very short time intervals associated with the formation of shear bands, adiabatic conditions might be adopted to simplify the resolution of the problem. However, this assumption has been shown to make the numerical solution mesh-dependent [14,37,60,61]. By considering the effect of thermal diffusion, the problem is given an intrinsic length scale that functions as a localization limiter, which leads to mesh insensitive results [9,10,56]. Other regularization techniques reported in the literature include strain gradient theories [2e4], which have been used in the context of shear bands in Refs. [1,33,53] and nonlocal methods [7,23,45]. Mesh alignment is another form of sensitivity that can be improved using mesh-free formulations [35,36] or Isogeometric analysis [16]. Local stability analysis of shear bands has been an active area of research as it provides valuable information that can be used to model shear bands through phenomenological models. Some

http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004 0734-743X/© 2015 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

2

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

examples of the importance of stability detection can be found in the work of Rabczuk et al. [46] that set strong displacement discontinuities with cohesive surfaces where local instability is detected. Song et al. [49] introduce additional tangential degrees of freedom at the discontinuity using phantom nodes. Belytschko et al. [15] inject coarse scale discontinuities based on isolating the unstable regions at the fine scale into a localization domain. Tabarraei et al. [51] use material instability based on the perturbation method to inject a phantom node discontinuity on a macro level model based on the cohesive laws obtained from a microstructural model. Finally, Rabczuk and Samaniego [47] use instability in an adaptive mesh-free method for discontinuous modeling of shear bands. Among others, two main approaches on the study of stability involving shear bands are the linear perturbation method [5,8,10,12,13,22,26,31,38,44,55] and the eigenvalue analysis of the acoustic tensor [18,19,50,59]. This paper focuses on the linear perturbation method and not on the loss of ellipticity detected by the eigenvalues of the acoustic tensor since for a rate-dependent material the equations remain elliptical as was shown by Lemonds and Needleman [34]. The stability of a system is traditionally studied from the perspective of the modal spectrum of that system. The underlying reasoning is that whichever mode has the largest growth rate corresponding to an exponential growth will dominate the solution. However, for non-normal systems (i.e. systems whose eigenvalues are non-orthogonal) the short time growth-rate of the solution might not correspond to the growth rate of the most unstable mode. On the other hand, non-autonomous systems do not have a constant spectrum because of the time dependency of their operators. Since the problem being studied is neither normal nor autonomous, a generalized stability analysis is of interest [24,25]. In this work we use the Finite Element Method to obtain a monolithic semi-discrete Jacobian for the shear band problem following the derivation in Refs. [41], on which a generalized stability analysis is performed. The local stability criteria proposed by Bai [8], Ling and Belytschko [38] and others, are usually derived from a perturbation analysis, which under the same assumptions is equivalent to a local eigenvalue analysis [6]. These works have shown that this instability condition predicts well the localized region and, for a certain special class of problems, it will mark the point after which the solution will rapidly localize. We focus on these problems and compare the results of the proposed generalized stability analysis with this local instability criterion. There are several different concepts that are important for a generalized stability analysis as shown in Farrell and Ioannou [24,25]. The study of the growth rate of a perturbation tangent to the solution is of particular interest since it is computationally cheap and can easily be implemented on complex multiphysics systems. Additionally, an analysis based on the concept of the Rayleigh Quotient of an eigenvector is employed and is shown to also produce a good approximation of the instability point. It has been shown that the local spectral analysis based on element eigenvalues recover the analytical local instability criterion [6]. Albeit considerably more expensive, the local spectral analysis posses the same implementation simplicity that the instantaneous growth rate and the Rayleigh Quotient have, i.e. they only require the knowledge of the tangent stiffness matrix. However, a global modal analysis based on the eigenvalues of the full system is computationally prohibitive as the problem grows larger. This fact makes the proposed approach attractive for the study of the global instability behavior since the computational cost is of the same magnitude of the analytical instability criterion.

This work is structured as follows. Section 1 contains the current introduction. In Section 2 the shear band problem is described mathematically, followed by its implementation using the Finite Element Method. Section 3 describes the method of local instability analysis that has been employed in the past and compared with experimental results and the Generalized stability method used in this work. In Section 4 the application of the proposed methodology on numerical examples is showed for both 1D and 2D cases. Finally, the conclusions are presented in Section 5 after which the Appendixes and the References are located.

2. Problem statement 2.1. Shear bands The shear band problem is described by the equations of conservation of momentum and energy, and also the equations for elastic and inelastic constitutive relations. All of these and boundary conditions are formalized as a set of coupled PDEs as follows. The Momentum Equation, which includes inertial effects but ignores body forces, is given by

€ ¼ V$s ru

(1)

with u as the displacement field, t as the time, s as the stress tensor and r as the material density. A super imposed dot corresponds to a time derivative. The Elastic Constitutive Relation, in total form, is

s ¼ Celas : ðε  εp Þ

(2)

where Celas is the fourth-order elastic constitutive tensor. The strain tensor has been additively decomposed into elastic and inelastic parts [48] as

ε¼

i 1h Vu þ ðVuÞT ¼ εe þ εp 2

(3)

with εe and εp as the elastic and inelastic parts of the strain, respectively. The energy equation is written as

rcp T_ ¼ lV2 T þ cs : Dp

(4)

where T is the temperature, and l, cp, and c are the conductivity, specific heat and Taylor-Quinney coefficient, respectively. This equation accounts for diffusion as well as heat production in proportion to the plastic work, s:Dp [52], where Dp ¼ ε_ p is the rate of plastic deformation. The rate of change of the equivalent plastic strain is defined as

g_ p ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 p D : Dp 3

(5)

The inelastic part of the rate of deformation tensor, assuming J2 plasticity, is defined as

Dp ¼

3 g_ p S 2 s

(6)

where the deviatoric stress tensor S is

1 S ¼ s  trðsÞI 3

(7)

and the effective stress is given by

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

s¼

rffiffiffiffiffiffiffiffiffiffiffiffiffi 3 S:S 2

Specifically, the residual of the momentum equation (14) reads

(8)

The inelastic constitutive relation, given in terms of the effective stress, may be written as

    s ¼ PðTÞQ gp R g_ p

(9)

where gp is the equivalent plastic strain, which serves as a strain hardening parameter and the functions P(T), Q ðgp Þ and Rðg_ p Þ depend on specific material models. Any flow law that can be decomposed as in (9) is valid [29] and hereafter the JohnsoneCook material model [30] will be used. The four equations (1), (2), (4) and (9) along with the definition of the flow law describe the evolution of the four unknown fields of displacement, stress, temperature and plastic strain. Lastly, the boundary conditions are

u¼u

on

Gu

(10)

T ¼T

on

GT

(11)

n$s ¼ t n$q ¼ q

3

on on

Gt G

(12)

q

Z U

2.2. Finite element formulation

Z Vwv sdU 

U

wv s$ndG

(19)

G

The weak form of the residual of the elastic constitutive equation (15) is

Z

Z _  ws sdU

rs ¼ U

ws C

elas

:

U

  ! 3 g s; T; gp S dU V v 2 s s

(20)

The residual of the energy equation (16) reads

Z

_ wT rcp TdU þ

rT ¼ U

Z

Z VwT lVTdU  U

Z 

wT c* htgdU

U

wT lVT$ndG

(21)

G

and finally, the weak form of the residual of the inelastic constitutive equation (17) is

Z rgp ¼

wgp g_ p dU 

U

(13)

where t and q are the prescribed traction and flux, respectively. This model considers small strains, and additionally neglects thermal strains and the thermoelastic contribution to the system energy. At t¼0 the system is considered to be unstressed, undeformed and at room temperature.

Z _ wv rvdU þ

rv ¼

Z wgp gdU

(22)

U

These equations can be grouped into a residual vector r and a solution vector x with displacements, temperature, stress and plastic strain fields as

3 v 6 T 7 7 x¼6 4s 5 gp 2

2

3 rv 6 rT 7 7 r¼6 4 rs 5 rgp

(23)

The coupled nonlinear problem can then be stated as A mixed finite element formulation is used to discretize the governing equations for the shear band problem. The matrices that result from that discretization are used for both solving the PDEs and analyzing the stability of the problem. The system of equations (1)e(9) is first reformulated in the following form

rv_ ¼ V$s s_ ¼ Celas :

(14) Vs v 

  ! 3 g s; T; gp S 2 s

  rcp T_ ¼ lV2 T þ csg s; T; gp   g_ p ¼ g s; T; gp

(15)

i 1h Vv þ ðVvÞT 2

(24)

where uibv, dropped henceforth for simplification, represents the initial and the boundary conditions. To solve (24), the Newton method is applied, assuming a Taylor series expansion around x0 and x_ 0 , such that

_ ¼ rðx0 ; x_ 0 Þ þ ddx rðx; xÞj _ x ;x_ þ ddx_ rðx; xÞj _ x ;x_ rðx; xÞ 0 0 0 0

(25)

^teaux differentials in the dx where ddx($) and ddx_ ð$Þ represent the Ga and dx_ directions, respectively, which are defined as

(16)

  d _  rðx þ εdx; xÞ ¼ K$dx dε ε¼0   d _  _ ¼ rðx; x_ þ εdxÞ ddx_ rðx; xÞ ¼ M$dx_ dε ε¼0

_ ¼ ddx rðx; xÞ (17)

with the rate of deformation tensor given by the symmetric part of the velocity gradient, or directly defined by taking the time derivative of Eq. (3), and reads

Vs v ¼

_ uibv Þ ¼ 0 rðx; x;

(18)

To obtain the weak form in mixed formulation, each equation is pre-multiplied by an appropriate shape function wv, ws, wT and wgp for the momentum equation, the elastic constitutive equation, the energy equation and the inelastic constitutive equation, respectively. Then, integrating over the domain and applying the divergence theorem where appropriate, the weak form of the residual equations is obtained.

(26)

_ ¼ 0, gives Considering that rðx; xÞ

rðx0 ; x_ 0 Þ ¼ M$dx_  K$dx

(27)

It is interesting to note that (27) represents both the equation for spectral instability analysis by the Indirect method of Lyapunov (explained later) and the equations for Newton's method for solving the nonlinear problem with an implicit scheme. To emphasize this point, note that if the initial point ðx0 ; x_ 0 Þ is a solution, then its residual will vanish and (27) will translate into the Lyapunov stability equation in (40). On the other hand, if the residual does not vanish, then a typical Newton correction is obtained, as follows

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

4

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

Jdx ¼ r

(28)

where J corresponds to the Jacobian matrix (often called Tangent ^teaux derivStiffness matrix) and is computed by applying the Ga ative to the appropriate residual. Depending on the time discretization chosen, the Jacobian matrix can be expressed as a linear combination of the matrices M and K. If the process is assumed to be quasi-static, the Mass matrix M will simply become zero. Furthermore, for a Backward Euler discretization in time, the Jacobian matrix is given by

J¼

1 MK Dt

(29)

where Dt is the time step used to advance the solution. Note that the time integration scheme will affect the Jacobian matrix used to obtain a solution at a particular time step but not the matrices M and K used for the stability analysis. This happens because M and K are taken after a converged solution at a given time is reached and are therefore independent from the integration scheme. The block structure of the matrix M is given by

2

Mv 6 0 M¼6 4 0 0

0 Ms 0 0

0 0 MT 0

3 0 0 7 7 0 5 Mgp

(30)

0 6 KLsv K¼6 4 0 0

KLvs KLss KNLTs KNLgp s

0 KNLsT KLTT þ KNLTT KNLgp T

e Stage II: Inhomogenous plastic deformation e Stage III: Highly localized deformation and stress collapse It is worth noting that the concept of homogeneity in this context refers to a homogeneous (as opposed to heterogeneous) distribution of the physical quantities that describe the problem (strain, stress, temperature, etc …) where changes occur smoothly and gradually over the domain. Since it is at stage II that inhomogeneous plastic deformation begins, it is expected that the instability criterion will predict the initiation of this stage. In fact, as shown by Bai [8], the analytical local instability condition for this problem is given by:

1

and the matrix K by

2

Fig. 1. Stages of Shear band deformation. I-homogeneous solution, II-onset of localization, III-stress collapse [40].

0

3

KNLsgp 7 7 KNLTgp 5 KNLgp gp

(32)

where the subscript 0 indicates the value of that quantity at the current time and

(31)

Here M is the mass matrix and K ¼ KL þ KNL represents the sum of the stiffness matrices associated with linear material behavior such as elasticity and thermal diffusion (KL) and the tangent stiffness matrices associated with material non-linear behavior (KNL). The inner blocks in M and K are presented in Arriaga et al. [6]. The subscripts on each block indicate the fields which are coupled, thus the pair (i,j) represents the variation in the j direction of the residual associated with field i. 3. Theoretical background This section focus on the stability analysis of the problem. In section 3.1 the local analytical instability criterion is recalled, evidentiating the correspondence between this result and experimental observations. It is also shown why this analytical condition, commonly used in the literature, is equivalent to the spectral stability analysis in Arriaga et al. [6]. The presentation of these topics is essential for validation of the proposed criterion. In section 3.2 the method of generalized stability analysis is used to derive a set of additional conditions for which instability is expected. There, the relationship between this method and the spectral stability analysis is shown, motivating the prediction of the results presented in section 4. 3.1. Local instability According to Marchand and Duffy [40], three stages in a shear band deformation process can be observed in experiments and are identified in Fig. 1. These stages are. e Stage I: Homogenous plastic deformation

cs0 P0 <0 rcp Q0

 vs P0 ¼   ; vT x¼x0

 vs  Q0 ¼  vgp 

(33) x¼x0

It has been shown in Arriaga et al. [6] that the local instability can be computed by doing an eigenvalue analysis of the element stiffness matrix. This was more expensive than the computation of the analytical local instability criterion but was much easier to implement in more complex situations. Also, this result shows that the local criterion is in fact a spectral stability analysis. The technique used to prove this result is the indirect method of Lyapunov [43] which is equivalent to the linear perturbation method. This technique is briefly explained here to provide insight on the assumptions and simplifications being made and to show the correspondence with spectral stability analysis. Firstly, consider the system of equations expressed in its autonomous form as

_ ¼ FðxÞ L ðxÞ

(34)

where x represents all non-constant variables of the problem, L represents a linear operator and F(x) is a non-linear function of x. The variables x can be expressed as a sum of the equilibrium solution with a perturbation dx as:

x ¼ x0 þ dx

(35)

Linearizing F(x) by a Taylor Series Expansion around x0 and considering only the first order term gives

FðxÞ ¼ Fðx0 Þ þ F 0 ðx0 Þdx

(36)

Thus, (34) can be written as

_ ¼ Fðx0 Þ þ F 0 ðx0 Þdx L ðx_ 0 þ dxÞ

(37)

Since x0 is a solution point, L ðx_ 0 Þ ¼ Fðx0 Þ, which means that

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

_ ¼ K$dx L ðdxÞ

(38)

0

where K¼F (x0) is the linearized Jacobian matrix of F(x) at x0. In our model, this matrix will be obtained from mixed the finite element formulation described earlier. This process centers the stability analysis problem on the equilibrium point such that it becomes the origin, permitting the use of a Lyapunov indirect stability analysis [43]. Representing the operator L ($) by a matrix M and considering a linear perturbation of the form

b eut dx ¼ x

(39)

b its where u2ℂ is an eigenvalue at the complex space ℂ and x corresponding eigenvector representing the mode shape of the perturbation, gives the following generalized eigenvalue problem

uMdx ¼ Kdx

(40)

This method of stability analysis was proposed by Lyapunov [39] and has been referred to as the first method of Lyapunov or the indirect method of Lyapunov [20,43]. Here we showed that the linear perturbation method leads to an eigenvalue problem. On the other hand, by applying the linear perturbation method [8] the local instability criterion is obtained. Therefore, it is clear that the local instability criterion methodology (perturbation) is equivalent to a spectral stability analysis (eigenvalues).

3.2. Generalized stability of a perturbation tangent to the solution 3.2.1. Tangent perturbation An important quantity of interest for a general stability analysis is the growth rate of a perturbation tangent to the current solution. To understand the importance of this quantity it is useful to notice that a typical spectral analysis defines an instability by the existence of a mode with positive exponential growth rate in the system, in which case that mode will completely dominate the solution at t / ∞. However, this rationale assumes that the operators of the system do not change, or change very slowly in time. If the operators change very quickly, as in the shear band problem, then the growth rate of the unstable perturbation may not be significantly larger than the actual rate of change of the solution. Therefore, the growth of the perturbation might not be fast enough to dominate the solution before the system changes considerably. On the other hand, if the solution is already progressing in an unstable direction, it is expected that this unstable mode will dominate the progress of the solution. One way to obtain the vector that is tangent to the solution is by taking the increment of the solution on a first order integration scheme given by

xnþ1 ¼ xn þ xs

where x is the solution at step n and xs is the increment of the solution step. The Backward Euler integration scheme assumes that the rate of the increment is the same as the rate at the converged time step (i.e. the tangent of the solution at the converged time step). Hence, the incremental step is given by

xs ¼ Dt x_ nþ1

(where K and M were also computed for that time step) is in fact the tangent to the solution at xnþ1 and the desired vector for our study. One can immediately recognize that no extra effort was put into the calculation of this vector for the Backward Euler integration scheme. If other integration schemes were used, then the rate vector would have to be computed explicitly at the converged time step. Naturally, the discrete solution step might still be an acceptable approximation of the tangent if the time step is not excessively large. On the other hand, even in the Backward Euler scheme, if the time step is large then the operators might change significantly between time steps, rendering the tangent direction less accurate. 3.2.2. Instantaneous growth rate We define the instantaneous growth rate of a general perturbation vector. For consistency with the literature, let us first define the linearized dynamic operator as

A ¼ M1 K

(42)

with Dt is the time increment and x_ nþ1 is the rate of change of the solution at step nþ1. This means that the vector xs that represents the step from xn to the converged solution xnþ1 obtained with the Newton Method

(43)

Therefore, the perturbed problem for a converged solution, as seen in (27), may be expressed as

x_ ¼ A$x

(44)

Note that, generally speaking, A is not constant and is only valid for small perturbations with respect to the current solution. Assuming that A is autonomous (linearity is already implied when looking at the perturbed system), the solution at time t is given by

xðtÞ ¼ eAt x0

(45)

with t being the time in reference to the current time and x0¼x(0) the perturbation vector which we want to study. The growth of that perturbation vector is therefore given by:

s2 ðx0 ; tÞ ¼

xðtÞT xðtÞ xT0 x0

(46)

plugging in (45) T

2

s ðx0 ; tÞ ¼

xT0 eA t eAt x0

(47)

xT0 x0

The instantaneous growth rate will therefore be given by

GI ðx0 Þ ¼ lim

t/0

vs2 ðx0 ; tÞ vt

(48) T

Since x0 is not a function of t we can expand the matrix eA t eAt in its Taylor Series and compute its limit when t/0, which leads to:

(41)

n

5

GI ðx0 Þ ¼

xT0 HðAÞx0 xT0 x0

¼

xT0 Ax0 xT0 x0

(49)

where H(A)¼(ATþA)/2 is the Hermitian part of A. If we revert to using K and M, then

GI ðx0 Þ ¼

xT0 M1 Kx0 xT0 x0

(50)

Notice that, if x0 is an eigenvector, then is it clear that GI(x0) will be the eigenvalue of that eigenvector. This is an important result because it shows why the spectral analysis is usually considered a good way of detecting instability. The idea is that, if any mode has a positive real part, then that mode

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

6

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

will grow exponentially at a rate given by the eigenvalue of that mode. If that is the case, then the solution will be mostly changing in that direction and therefore the tangent of the solution should be very close to an eigenvalue. As can be observed in Appendix A, in a system with orthogonal operators, the instantaneous growth rate behaves as a weighted average of the eigenvalues scaled by the modal participation factors. Therefore, in this type of problem GI(x0)max(u) with u being the eigenvalues of the system. However, it is known that nonnormal systems might have growth rates larger than the eigenvalues of that system. Regardless, if the tangent to the solution is close to an eigenvector, then the difference between GI(x0) and max(u) should be small. 3.2.3. Rayleigh Quotient The previous procedure is used to find the instantaneous growth rate of any vector. However, as mentioned before, the tangent direction is expected to be close to an eigenvector with an unstable eigenvalue. Therefore, an approximate way of computing the growth rate is by using the well known Rayleigh Quotient given as:

GRQ ðx0 Þ ¼

xT0 Kx0 xT0 Mx0

(51)

Once again, if x0 is an eigenvector, then GRQ(x0) will be the eigenvalue of that eigenvector. As can be observed in Appendix B, in a system with orthogonal operators, the Rayleigh Quotient is also a weighted average of the eigenvalues scaled by the modal participation factors, but now also scaled by the modal masses. The immediate advantage of the Rayleigh Quotient over the Instantaneous Growth rate is the fact that it does not require the inversion of M. Since M is positive definite and well conditioned, this operation will be cheap using iterative methods. However, for large problems and many time steps, this cost might become relevant. Another advantage is that static modes1 will make M singular and the inversion of M impossible or very complex (pseudo inverse). This of course should not happen if the Mass matrix is positive definite. Nevertheless, referring to the example for a system with orthogonal operators presented in Appendix B, we notice that if we have a mode that is very close to being a static mode, i.e. with a small modal mass, it is likely that its eigenvalue will be very large. In the instantaneous growth grate of a perturbation formulation in equation (50) nothing is done to counter the influence of this mode, but the Rayleigh Quotient scales down this mode based on its low mass. This becomes a significant concern in non-normal systems since the eigenmodes of the operators are not orthogonal, so that almost insignificant subspaces might imply very large modal components (see Appendix C for a pictorial example).

Table 1 Material properties for 4340 Steel. Property name

Symbol

Value

Unit

Mass density Specific heat Thermal conductivity Young's modulus Poisson's ratio Shear modulus Taylor-Quinney coefficient

r

7830 477 38 200 0.29 77.5 0.9

kg/m3 J/kgK W/mK GPa e GPa e

In this section we choose numerical examples where the local instability triggers the solution to start evolving mainly in the direction of localization. In other words, when local instability happens, the direction of the solution change will almost immediately be dominated by the localization direction. These type of problems usually model very small scales where the elastic wave propagation plays an insignificant role. Problems with larger scales will have other effects included in the direction of the solution change, like elastic wave propagation and reflection, stress redistribution, shear band propagation, etc. These effects counter the property of the

1

Static modes: Modes with modal mass equal to zero.

l E

m G

c

Table 2 JohnsoneCook parameters for 4340 Steel. Property name

Symbol

Value

Unit

Yield shear stress Shear stress hardening parameter Strain hardening parameter Reference temperature Melting temperature Thermal softening exponent Strain-rate hardening parameter Reference strain-rate

A B N T0 Tm m c g_ pr

792 510 0.26 298 1793 1.03 0.014 1.0

MPa MPa e K K e e 1/s

small scale examples where the local instability triggers the global localization. We study the behavior on 1D and 2D shearband benchmark problems. The material considered in both examples is a 4340 Steel with a JohnsoneCook constitutive law [30] given in (52). The parameters associated with this material are presented in Tables 1 and 2.

!#    " h  N i g_ p T  T0 m 1 1 þ c ln s ¼ A þ B gp Tm  T0 g_ pr

(52)

Neumann boundary conditions of zero flux for the Temperature field are applied on all edges and plane strain conditions are assumed when applicable. The small strain simplification is used without loss of generality, since the focus of this work is on the behavior of the instability measures and not the quantitative results obtained from the modeling. An hyperbolic secant type of imperfection is used in order to trigger the shearband. Hence, in the analysis we scale by himp(rn) the parameters A and B, with himp(rn) given by

 himp ðrn Þ ¼ 1  ared

4. Numerical results

cp

2 ern þ ern

 (53)

where ared is the reduction of the material parameters in percentage at the center of the imperfection and rn is the normalized distance given by

Fig. 2. Velocity profile used as a loading function to generate the shearband.

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

7

Fig. 3. Representation of a 1D right half-rod. The black dot marks b.

rn ðx; yÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx  x0 Þ2 þ ðy  y0 Þ2

(54)

r0

with x0 and y0 representing the center of the imperfection and r0 its radius. A velocity-type loading is applied as dirichlet boundary conditions. Velocity control is defined as shown in Fig. 2 with vr being the steady velocity and tr the time to reach that velocity. The S-curve is a polynomial that guarantees that the derivative of the Jerk2 is continuous. We choose this to avoid discontinuities and kinks in the profile of the Jerk since the thermo-mechanical coupled system is sensitive to the applied velocity and may generate jumps in the instantaneous stability measure.

the half-rod, which means that the final length is 103 m. The total number of nodes is 401 since the central node is repeated. Fig. 3a shows the relative position of the nodes in the half-rod as a function of the relative node number and Fig. 3b shows the relative element size also as a function of the relative node number. As defined before, the material considered for the 1D problem is a 4340 Steel with a JohnsoneCook constitutive law [30] and an hyperbolic secant type of imperfection is considered at the center of the rod with ared¼0.01, r0¼10m m and x0¼L/2¼500m m. The analysis is performed on a range of nominal strain rates between 101 s1 and 104 s1 by applying the velocity profile shown in Fig. 2 to both edges of the rod. Our goal is to detect the onset of shearband localization using the proposed methods: Instantaneous growth rate method and Rayleigh Quotient method and compare those with the analytical criteria proposed by Bai [8].

4.1. 1D rod under uniform shear 4.1.1. Problem description A steel rod is studied using a 1D formulation often employed as an idealized version of the Torsional Kolsky Bar experiment [21,28,32,54], which is commonly used to study a pure shear deformation without necking effects [11,58]. To capture the shear band localization more accurately, the rod is discretized with a non-constant element size. Considering symmetry, the relative position of the elements (p) from the center to one edge (half-rod) can be defined by the following function

pðxÞ ¼

8 > < S0 x



xb > : S0 x þ ð1  S0 Þ 1b

P

for x  b for x > b

(55)

where x the relative node number, P is the exponent used for varying the element size, b is the percentage of nodes concentrated at the center of the rod with constant element size and S0 is the slope for the central elements, defined as:

S0 ¼

de N L

(56)

where L and N are respectively the length and the number of nodes of the half-rod and de defines the constant element size in the center of the rod. In the following sections, the results are presented for N¼201, b¼0.25, P¼5, de/L¼1103 and L¼0.5103 m. All values refer to

2 The Jerk is the derivative of the Acceleration with respect to time and is the quantity responsible for generating elastic waves.

4.1.2. Results Consider Fig. 4 that represents the shear stress at the edge of the rod as a function of the nominal strain (Fig. 4a) and the equivalent plastic strain along the rod at different moments in time (Fig. 4b). The various deformation states are illustrated by the different colors. It can clearly be observed that the plastic strain grows uniformly along the rod during the initial stage of the deformation (disregarding the slight inhomogeneity due to the imperfection). This is consistent with the experimental results of Marchand and Duffy [40] as it corresponds to Stage I of Fig. 1. Only after the onset of instability and the initiation of Stage II do we start to observe an inhomogeneous distribution of the solution with increasing growth rate. The special characteristic of this problem is precisely the fact that plastic deformation starts to develop homogeneously and that localization is globally initiated when Stage II starts. These types of problems are only realistic when considering a small portion of material at the vicinity of the shear band, and therefore are not big enough to model shear band propagation. However, they still provide a good way to model the region of the shear band and its immediate vicinity as a comparison with the experimental work of Marchand and Duffy [40] suggests. Fig. 5 shows the comparison in the detection of onset of shear band localization between the analytical stability criteria for the central point of the 1D rod (Eq. (32)), the Rayleigh Quotient (Eq. (51)) and the instantaneous growth rate (Eq. (50)). Five different nominal strain rates are tested. The gray line represents the stress at the center of the rod and the blue, green and red lines correspond to the analytical criterion, the Rayleigh Quotient and the instantaneous growth rate, respectively. It can be seen that all three measures agree and approximate well the local instability point.

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

8

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

Fig. 4. Behavior of the rod with progression of time. 4a shows the stress-strain behavior of a rod being sheared at a strain-rate of 104 s1. Several points are chosen for which the entire distribution of the plastic strain along the rod is represented in 4b. The red indicator marks the onset of instability. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The results in Fig. 5 can be summarized in a single plot as shown in Fig. 6. There, the nominal strain at the points where the local stability criterion, the Rayleigh Quotient and the instantaneous growth rate change sign and the peak stress at the center of the rod are represented as a function of the strain-rate. As expected, the instability has a trend similar to the peak stress and all measures for instability agree. For a simple problem like this, the peak stress is occurring roughly at a given critical value of accumulated inelastic work since plastic strain is the single source of the softening and

the solution is spatially homogeneous up to the peak (ignoring the effect of the imperfection). Given that diffusion does not play a role when the solution is homogeneous, we would expect that the integral of the stress-strain curve from zero to the peak stress to be independent of the strain rate. Therefore, since the yield stress increases with strain rate, the strain for the peak stress would have to reduce slightly to maintain the value of accumulated inelastic work. This can be observed in Fig. 5 from the downward slope of the maximum stress curve with respect to the strain rate.

Fig. 5. Stability analysis of 1D rod (center point) at different strain rates. The local instability condition given by Eq. (32) (where a negative value corresponds to instability) is compared to the Rayleigh Quotient as given in (51) and the Instantaneous Growth rate as given in (50).

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

Fig. 6. Instability points as a function of the strain rate.

It becomes clear that both methods proposed in this paper for detection of the instability produce a good approximation of the instability point and agree well with the analytical criteria of Bai [8]. Consequently, it is possible to use these methods to provide a cheap and convenient alternative to detect the instability point when it is difficult to formulate the analytical criterion, which is the case in general general shear band problems. 4.2. 45 shear band under uniform compression 4.2.1. Problem description In this section a square plate is compressed under uniform  impact loading generating a shear band oriented at a 45 angle. The geometry of the square plate is given in Fig. 7 with the edge dimension of D ¼ 100 mm. The mesh consists of 4040 quad elements, modeled by an irreducible mixed-finite element formulation with B-bar to reduce shear-locking effects. The velocity and temperature fields are discretized with bilinear C0 continuous functions and the stress and equivalent plastic strain are piecewise C1 continuous. Details on the numerical implementation can be found in Arriaga et al. [6], Berger-Vergiat et al. [16], McAuliffe and Waisman [41,42]. As defined before, the material considered is a 4340 Steel with a JohnsoneCook constitutive law [30] and an hyperbolic secant type of imperfection (defined earlier in equation (53) and graphically represented in Fig. 7) is considered with ared¼0.1, r0¼D/10 and x0¼y0¼0m m. Note that this location (indicated in Fig. 7 as BL) will be the first to become unstable and for that reason it is used as the reference location for the local instability behavior.

Fig. 7. A square plate under uniform impact loading. The colors represent the factor himp defined in equation (53). Also indicated are two corner points that are used for plotting purposes: the bottom-left corner (BL) and the top-right corner (TR).

9

Fig. 8. Evolution of Von-Mises stress at the maximum imperfection point (BL-bottom left corner), the shear band imperfection free point (TR - top right corner) and the average Von-Mises stress on the plate as function of the nominal strain.

The plate is compressed uniaxially with a velocity control loading following the profile shown in Fig. 2 with vr ¼ 10 m/s and tr ¼ 0.4 ms, which corresponds to a nominal strain rate of 105 s1 for this plate. 4.2.2. Results Similarly to the 1D case, this model is part of a particular set of examples where the inhomogeneous deformation triggered locally develops almost simultaneously in the entire plate. This effect comes as a consequence of the small size of the plate width and is realistic only for studying the shear band region and its immediate vicinity. Fig. 8 shows the comparison between the Von-Mises stress at the bottom left corner of plate where the imperfection is maximum (BL), the stress at top right corner where the shear band passes but the imperfection is practically zero (TR) and the average Von-Mises stress on the plate. This comparison shows that the response in the early stage of yielding is uniform (apart from the imperfection) and that the nature of the differences between stresses are only

Fig. 9. Stability analysis of the 2D plate at the reference position (BL) for a strain rate of 105 s1. The local instability condition given by Eq. (32) (where a negative value corresponds to instability) is compared to the Rayleigh Quotient as given in (51) and the Instantaneous Growth rate as given in (50).

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

10

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

Fig. 10. Representation of the equivalent plastic strain field for the three different stages mentioned in section 3.1.

geometrical and not temporal. This illustrates well that both the elastic wave propagation speed and the shear band propagation speed are too large to be causing any significant effect on a plate of the size being modeled. Once again, this is only realistic when studying the vicinity of the shear band since for larger scales elastic wave propagation and shear band propagation play significant roles. Fig. 9 shows the values of the analytical stability criterion (Eq. (32)), the Rayleigh Quotient (Eq. (51)) and the instantaneous growth rate (Eq. (50)) for the reference position (BL). The gray line represents the Von-Mises stress and the blue, green and red

lines correspond to the analytical criterion, the Rayleigh Quotient and the instantaneous growth rate, respectively. It is shown that both quantities approximate well the local instability point. Additionally, notice that the growth rate stabilizes in the long term since the solution is now evolving in a steady-state like path, where the heat generated at the shear band is balanced by the thermal diffusion. In Fig. 10 the equivalent plastic strain field is represented at four different moments in time showing the three stages mentioned earlier in section 3.1. Fig. 10b shows a moment before the instability point, where the deformation, even though it has

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

considerable plasticity, is still homogeneous (besides the imperfection). The second figure (Fig. 10c) shows a moment immediately after the instability point and once again, apart from the effect of the imperfection, the velocity field looks homogeneous. The next figure (Fig. 10d) shows an intermediate time in stage II and mild localization can be observed. Finally, the last figure (Fig. 10e) shows a time well into stage III where extreme localization is visible. Notice that the shear band appears to be relatively wide in this representation due to the small size of the model. However, taking the dimensions of the plate into consideration, the shear band is of the order of 10e30 mm. As expected, the non-homogeneous deformation associated with a shear band only occurs after the instability point that initiates stage II. This confirms the predictive capability of the proposed method when compared to experimental results.

5. Concluding remarks

x ¼ Qv

(A.2)

where Q is the matrix that defines the eigenspace of the operators and v the components on the eigenspace resulting from the decomposition of x. Therefore

GI ¼

vT Q T M1 Q Q 1 KQ v vT Q T Q v

1

GI ¼

vT M K v vT v

i

v2i ui

GI ¼ P i

Appendix B. Rayleigh Quotient on an system with orthogonal operators The Rayleigh Quotient of a vector x is given by

GRQ ¼

xT Kx xT Mx

(B.1)

a modal decomposition of x yields

x ¼ Qv

(B.2)

where Q is the matrix that defines the eigenspace of the operators and v the components on the eigenspace resulting from the decomposition of x. Therefore

GRQ ¼

vT Q T KQ v vT Q T MQ v

vT K v v T M v

¼

(B.3)

with M*¼QTMQ and K*¼QTKQ. Since M* and K* are diagonal because of orthogonality:

GRQ ¼ P

Appendix A. Instantaneous Growth rate on an system with orthogonal operators

(A.5)

v2i

with vi as the ith component of the modal decomposition of x, ui ¼ Ki as the eigenvalue associated with the ith mode and Ki and Mi as Mi the ith diagonal entry of K* and M*, respectively.

i

The financial support of the Department of Energy through the Early Career Research Program, No. DE-SC0008196, is gratefully acknowledged.

(A.4)

with M*¼QTMQ and K*¼QTKQ. Since M* and K* are diagonal because of orthogonality:

P Acknowledgment

(A.3)

If the operators are normal, then Q1¼QT, which gives

P

In this work we have shown that the instantaneous growth rate of the vector tangent to the solution, derived from general stability analysis, predicts reasonably well the onset of shear band localization. This technique can easily be implemented within existing implicit FEM codes since the matrices and vector required to obtain the growth rates are the same as those computed for the Newton iterations. Moreover, we have also demonstrated that the instantaneous growth rate approach can be simplified into a special choice of the Rayleigh Quotient approximation which requires a few matrix vector products and therefore leads to a more efficient implementation. The Rayleigh Quotient approximation of the instantaneous growth rate works well when the global localization of the model is triggered by the appearance of local instability. Additionally, the proposed methodology is particularly useful to study problems for which an analytical criterion depends on functions that are fitted to experimental results, for example when different constitutive laws are used, thus making the analytical criterion case specific and hence more cumbersome to derive. Another advantage of this methodology is the possibility to modify the equations to include additional phenomena like nonlocal techniques, phase field or non-constant parameters like the Taylor Quinney coefficient, thus avoiding having to derive the analytical criterion for it. This is of significant importance because, even though the analytical criterion might convey more information about the location of the unstable regions of the material, it is usually difficult to derive and simplifying assumptions might be necessary to obtain it.

11

i

v2i Ki v2i Mi

P i

v2i ui Mi

¼ P i

v2i Mi

(B.4)

with vi as the ith component of the modal decomposition of x, ui ¼ Ki th mode and Ki and Mi > 0  as the eigenvalue associated with the i Mi th * * as the i diagonal entry of K and M , respectively. Appendix C. Non-orthogonality amplification

The instantaneous growth rate of a vector x is given by

xT M1 Kx GI ¼ xT x a modal decomposition of x yields

(A.1)

In Figure C1 is shown the effect of doing a decomposition of a vector in a non-orthogonal basis. We can observe that the components of a vector in this basis might have a disproportionate weights when compared to the actual magnitude that they represent.

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

12

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13

Figure C1: Amplification of the relative weight of the components on a non-orthogonal decomposition. References [1] Abu Al-Rub RK, Voyiadjis GZ. A finite strain plastic-damage model for high velocity impact using combined viscosity and gradient localization limiters: Part I - theoretical formulation. Int J Damage Mech 2006;15(4):293e334. [2] Abu Al-Rub RK, Voyiadjis GZ. A physically based gradient plasticity theory. Int J Plast 2006;22(4):654e84. [3] Aifantis EC. On the microstructural origin of certain inelastic models. J Eng Mater Technol 1984;106(4):326e30. [4] Aifantis EC. The physics of plastic deformation. Int J Plast 1987;3(3):211e47. [5] Anand L, Kim K, Shawki T. Onset of shear localization in viscoplastic solids. J Mech Phys Solids 1987;35(4):407e29. [6] Arriaga M, McAuliffe C, Waisman H. Onset of shear band localization by a local generalized eigenvalue analysis. Comput Methods Appl Mech Eng 2015;289: 179e208. sek M. Nonlocal integral formulations of plasticity and damage: [7] Ba zant ZP, Jira survey of progress. J Eng Mech 2002;128(11):1119e49. [8] Bai Y. Thermo-plastic instability in simple shear. J Mech Phys Solids 1982;30(4):195e207. [9] Batra R. The initiation and growth of, and the interaction among, adiabatic shear bands in simple and dipolar materials. Int J Plast 1987;3(1):75e89. [10] Batra R, Chen L. Effect of viscoplastic relations on the instability strain, shear band initiation strain, the strain corresponding to the minimum shear band spacing, and the band width in a thermoviscoplastic material. Int J Plast 2001;17(11):1465e89. [11] Batra R, Kim C. Analysis of shear banding in twelve materials. Int J Plasticity 1992;8(4):425e52. [12] Batra R, Wei Z. Instability strain and shear band spacing in simple tensile/ compressive deformations of thermoviscoplastic materials. Int J Impact Eng 2007;34(3):448e63. [13] Baucom JN, Zikry MA. Perturbation analysis of high strain-rate shear localization in B.C.C. crystalline materials. Acta Mech 1999;137(1e2):109e29. [14] Belytschko T, Chiang H-Y, Plaskacz E. High resolution two-dimensional shear band computations: imperfections and mesh dependence. Comput Methods Appl Mech Eng 1994;119(1e2):1e15. [15] Belytschko T, Loehnert S, Song J-H. Multiscale aggregating discontinuities: a method for circumventing loss of material stability. Int J Numer Methods Eng 2008;73(6):869e94. [16] Berger-Vergiat L, McAuliffe C, Waisman H. Isogeometric analysis of shear bands. Comput Mech 2014;54(2):503e21. http://dx.doi.org/10.1007/s00466014-1002-8.

[17] Bigoni D. Bifurcation and instability of non-associative elastoplastic solids. In: Petryk H, editor. Material instabilities in elastic and plastic solids. New York: Springer; 2000. p. 1e52. [18] Bigoni D, Zaccaria D. On the eigenvalues of the acoustic tensor in elastoplasticity. Eur J Mech A/Solids 1994;13(5):621e38. [19] Borja RI. Bifurcation of elastoplastic solids to shear band mode at finite strain. Comput Methods Appl Mech Eng 2002;191(46):5287e314. [20] Chen G. Stability of nonlinear systems. In: Webster J, editor. Wiley encyclopedia of electrical and electronics engineering. New York: John Wiley & Sons, Inc; 2001. p. 4881e96. [21] Chen WW, Song B. Split hopkinson (Kolsky) bar: design, testing and applications. New York: Springer; 2010. [22] Dai L, Bai Y. Basic mechanical behaviors and mechanics of shear banding in BMGs. Int J Impact Eng 2008;35(8):704e16. [23] Eringen AC. Nonlocal continuum field theories. New York: Springer; 2002. [24] Farrell BF, Ioannou PJ. Generalized stability theory. Part I: autonomous operators. J Atmos Sci 1996;53(14):2025e40. [25] Farrell BF, Ioannou PJ. Generalized stability theory. Part II: nonautonomous operators. J Atmos Sci 1996;53(14):2041e53. [26] Fressengeas C, Molinari A. Instability and localization of plastic flow in shear at high strain rates. J Mech Phys Solids 1987;35(2):185e211. [27] Gajo A, Bigoni D, Wood D. Multiple shear band development and related in stabilities in granular materials. J Mech Phys Solids 2004;52(12):2683e724. [28] Hopkinson B. A method of measuring the pressure produced in the detonation of high explosives or by the impact of bullets. In: The Royal Society, editor. Philosophical transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character. London: Royal Society of London; 1914. p. 411e3. 213(497508). [29] Hor A, Morel F, Lebrun J-L, Germain G. Modelling, identification and application of phenomenological constitutive laws over a large strain rate and temperature range. Mech Mater 2013;64:91e110. [30] Johnson GR, Cook WH. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Eng Fract Mech 1985;21(1):31e48. [31] Kim KH. Shear localization in viscoplastic solids [Ph.D. dissertation]. Cambridge: Department of Mechanical Engineering, Massachusetts Institute of Technology; 1987. PhD thesis. [32] Kolsky H. An investigation of the mechanical properties of materials at very high rates of loading. Proc Phys Soc Sect B 1949;62(11):676e700. [33] Lages EN, Paulino GH, Menezes IFM, Silva RR. Nonlinear finite element analysis using an object-oriented philosophy e application to beam elements and to the cosserat continuum. Eng Comput 1999;15(1):73e89. [34] Lemonds J, Needleman A. Finite element analyses of shear localization in rate and temperature dependent solids. Mech Mater 1986;5(4):339e61. [35] Li S, Hao W, Liu WK. Mesh-free simulations of shear banding in large deformation. Int J Solids Struct 2000;37(48e50):7185e206. [36] Li S, Liu WK. Numerical simulations of strain localization in inelastic solids using mesh-free methods. Int J Numer Methods Eng 2000;48(9):1285e309. [37] Li S, Liu WK, Rosakis AJ, Belytschko T, Hao W. Mesh-free galerkin simulations of dynamic shear band propagation and failure mode transition. Int J Solids Struct 2002;39(5):1213e40. [38] Ling X, Belytschko T. Thermal softening induced plastic instability in ratedependent materials. J Mech Phys Solids 2009;57(4):788e802. [39] Lyapunov AM. General problem of the stability OfMotion. 1st ed. London: CRC Press; 1992. [40] Marchand A, Duffy J. An experimental study of the formation process of adiabatic shear bands in a structural steel. J Mech Phys Solids 1988;36(3): 251e83. [41] McAuliffe C, Waisman H. Mesh insensitive formulation for initiation and growth of shear bands using mixed finite elements. Comput Mech 2013;51(5):807e23. [42] McAuliffe C, Waisman H. A PianeSumihara type element for modeling shear bands at finite deformation. Comput Mech 2014;53(5):925e40. [43] Murray RM, Li Z, Sastry SS. A mathematical introduction to robotic manipulation. Boca Raton: CRC Press; 1994. [44] Pan J. Perturbation analysis of shear strain localization in rate sensitive materials. Int J Solids Struct 1983;19(2):153e64. [45] Peerlings R, Geers M, de Borst R, Brekelmans W. A critical comparison of nonlocal and gradient-enhanced softening continua. Int J Solids Struct 2001;38(44e45):7723e46. [46] Rabczuk T, Areias PMA, Belytschko T. A simplified mesh-free method for shear bands with cohesive surfaces. Int J Numer Methods Eng 2007;69(5): 993e1021. [47] Rabczuk T, Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Comput Methods Appl Mech Eng 2008;197(6e8):641e58. [48] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [49] Song J-H, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. Int J Numer Methods Eng 2006;67(6): 868e93.  L. On the eigenvalues of the fourth-order constitutive tensor and loss [50] Szabo of Strong ellipticity in elastoplasticity. Int J Plasticity 1997;13(10):809e35. [51] Tabarraei A, Song J-H, Waisman H. A two-scale Strong discontinuity approach for evolution of shear bands under dynamic impact loads. Int J Multiscale Comput Eng 2013;11(6):543e63.

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004

M. Arriaga et al. / International Journal of Impact Engineering xxx (2015) 1e13 [52] Taylor GI, Quinney H. The latent energy remaining in a metal after cold working. Proc R Soc Lond Ser A 1934;143(849):307e26. [53] Voyiadjis GZ, Abu Al-Rub RK. A finite strain plastic-damage model for high velocity impacts using combined viscosity and gradient localization limiters: Part II - numerical aspects and simulations. Int J Damage Mech 2006;15(4):335e73. [54] Wright T. Approximate analysis for the formation of adiabatic shear bands. J Mech Phys Solids 1990;38(4):515e30. [55] Wright T. Shear band susceptibility: work hardening materials. Int J Plasticity 1992;8(5):583e602. [56] Wright T, Batra R. The initiation and growth of adiabatic shear bands. Int J Plasticity 1985;1(3):205e12. [57] Wright TW. The physics and mathematics of adiabatic shear bands. Cambridge, UK: Cambridge University Press; 2002.

13

[58] Wright TW, Walter JW. On stress collapse in adiabatic shear bands. J Mech Phys Solids 1987;35(6):701e20. [59] Xue L, Belytschko T. Fast methods for determining instabilities of elasticeplastic damage models through closed-form expressions. Int J Numer Methods Eng 2010;84(12):1490e518. [60] Zhou M, Needleman A, Clifton R. Finite element simulations of shear localization in plate impact. J Mech Phys Solids 1994;42(3):423e58. [61] Zhou M, Rosakis A, Ravichandran G. Dynamically propagating shear bands in impact-loaded prenotched platesdI. Experimental investigations of temperature signatures and propagation speed. J Mech Phys Solids 1996;44(6): 981e1006.

Please cite this article in press as: Arriaga M, et al., Instability analysis of shear bands using the instantaneous growth-rate method, International Journal of Impact Engineering (2015), http://dx.doi.org/10.1016/j.ijimpeng.2015.04.004