Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 48, 4, pp. 973-1001, Warsaw 2010 NOWACKI’S DOUBLE SHEAR TEST IN THE FRAMEWORK OF THE ANISOTROPIC THERMO-ELASTO-VISCOPLASTIC MATERIAL MODEL1 Adam Glema Tomasz Łodygowski Wojciech Sumelka Poznan University of Technology, Faculty of Civil and Environmental Engineering, Poznań, Poland e-mail: adam.glema@put.poznan.pl; tomasz.lodygowski@put.poznan.pl; wojciech.sumelka@put.poznan.pl The paper is dedicated to the memory of Prof. Wojciech K. Nowacki In the paper, the numerical simulation ofNowacki’s double shear test in the framework of recently proposed viscoplasticity theory for anisotropic solids is presented.Thenumerical analysis comprises the full spatialmodelling and is carried out for theDH-36 steel sheet in adiabatic conditions (the analysis of anisotropic bodies canbe led only on 3Dmodels).During analyses, strain rates of order 104-107 s−1 are observed and the process time duration up to full damage(loss of continuity in the localisationzone) is around 150-300µs. The novelty of the research is focused on the formulation that includes the anisotropy of the intrinsic microdamage process. Thus, it makes possible to obtain qualitatively and quantitatively new results compared with the existing models, like tracing the softening directions and better (closer to experiment) prediction of damage paths. Key words: microdamage anisotropy, metals, constitutive relation 1. Introduction Modern constitutive models need sophisticated experimental techniques for the identification of material functions and parameters assumed in the the- 1 The support of the Polish Ministry of Science and Higher Education under grant NN519419435 ”The evolution of properties and failure criteria of materials and structu- res under fast dynamic loadings” is kindly acknowledged. 974 A. Glema et al. oretical description. To themost challenging tasks from the experimental and simultaneously the theoretical points of view, concerning the analysis ofmetals behaviour (the central point in this paper), belong to the investigation of stra- in, stress, temperature and softening fields up to full damage during extremely fast dynamic processes. Such experimental evidences include crucial proper- ties like e.g., local strain rates of the order 104-107 s−1, temperature rise in the strain localisation zones up to the melting point and anisotropy induced by evolving the intrinsically anisotropic microdamage process. One of experimental techniques which allows for the analysis of the men- tioned phenomena is the recently developed test for dynamic plane shearing of steel sheet (Gary and Nowacki, 1994; Klepaczko et al., 1999) – which we call Nowacki’s double shear test after the name of its inventor. The experimental set-up was broadly analysed both experimentally and theoretically (Klepacz- ko et al., 1999; Nowak et al., 2007; Pęcherski et al., 2009). Nevertheless, we are going to present numerical calculations up to full material degradation in the specimen (loss of continuity in the softening zones) and, moreover, to includemodelling of the intrinsic microdamage process as an anisotropic one. Notice, that the anisotropy of microdamage plays a crucial role in the proper description of damage, giving (as it will be discussed in the following sections) qualitatively andquantitatively newresults comparedwith the existingmodels (Sumelka, 2009). Thus, the fundamental aim of the paper concerns the numerical simulation of Nowacki’s double shear test in the framework of viscoplasticity theory for anisotropic solids (Perzyna, 2008; Glema et al., 2009; Sumelka, 2009). Section 2 deals with the description of Nowacki’s double shear test. The direct impact variation of this experimental technique is selected. The crucial results of this section state a base for proper numerical modelling presented in Section 4. In Section 3, the fundamental results of viscoplasticity theory for anisotro- pic solids are presented. The discussion on amicrodamage tensor as ameasu- re of directional softening, definition of kinematics of the body, constitutive postulates and constitutive relations with a complete definition of material functions for an adiabatic process are included. Section 4 comprises description of the numerical modelling. The material model implementation in a finite element software Abaqus by using the VU- MAT user subroutine is discussed. The numerical calibration of the material model forDH-36 steel is presented. Finally, the definition of the numericalmo- del for Nowacki’s experimental set-up with numerical results of the specimen shearing are shown. Nowacki’s double shear test... 975 2. Nowacki’s double shear test Let us focus the attention on the variation of the Nowacki’s double shear test for the high strain rates where loading results from the direct impact in Hopkinson transmitter tube (Klepaczko, 1994). The main advantage of the test considered is that for this type of loading path, large strains without the occurrence of a plastic instability are obtained and the application of the Hopkinson tube assures the reduction of Pochhammer-Chree vibrations what makes easier the constitutive modelling. The crucial element of the whole experimental set-up is the special devi- ce which transforms the compression into nearly homogeneous shear over the total length of the specimen, see Fig.1. The shear device, following the de- scription byKlepaczko et al. (1999), consists of two coaxial parts, the external stiff housing and the internal cylindrical element. Both parts of this device are divided into two symmetrical clamps. Between the clamps, the sheet speci- men is fixed by robust tightening of eight screws (six in the external part and two in the internal cylinder). The direct impact imposed on the cylindrical part induces the deformation of zones of the specimen between the internal and external clamps of the device producing the state of plane shear in their bounds. Fig. 1. Nowacki’s double shear experimental set-up (Nowacki et al., 2006) This experimental technique (also called ”block-bar” loading scheme in compression) enables investigation of the velocity of striker up to 200ms−1, givingglobal strain rates in thematerial tested reaching several thousands s−1. 976 A. Glema et al. For the maximal velocities of the projectile, the intensity of plastic energy dissipation through the localised band can be so high, that the temperature generated exceeds themeltingpoint leading to thephasechangeof thematerial tested –what states as a base formodernmaterial modelling including except elastic and inelastic ranges also the phase change. Nevertheless, the device can also be used for quasi-static conditions, so one can be sure of the reliability of experimental results, with a substantial reduction of scatter. Duringtheexperimental tests foradvanceddeformations, shear localisation zones are observed. Through the localised shear bands, the intensification of the intrinsic anisotropicmicrodamageprocess evolution is evidenced leading to final fracture (loss of continuity) in the specimen.Such results are fundamental for proper constitutive description of damage. It should be emphasised that the method described imposes the following restrictions on geometry of the specimen (see Fig.2): • length of the double shear specimen allowed by the shear device is up to l0 =30mm, • the shear zone should satisfy the condition: 2 < a0/t < 10 to avoid buc- kling of the sheared zone and damage of the specimen during fixation in the device (t denotes the sheet thickness) and a0/l0 ¬ 10−1 tominimise the error due to non-homogeneity of the shear stresses and strains at both ends of the specimen. The common specimen dimensions, which are also assumed in Section 4 con- cerning numerical analysis, are presented in Fig.2. Fig. 2. Geometry of the specimen Nowacki’s double shear test... 977 3. Constitutive model 3.1. Intrinsic anisotropic microdamage process 3.1.1. Experimental motivation Now, let us focus the attention on experimental observations of microda- mage anisotropy in metals being the source of the overall metal anisotropy. Other sources of anisotropy like different sizes and shapes of adjacent grains (Narayanasamy et al., 2009), presence of different phases like pearlite or ferrite (Pęcherski et al., 2009) are not discussed. So, the always existing defects inmetal structures likemicrocracks, micro- voids,mobile and immobile dislocations densities (Abu-Al-Rub andVoyiadjis, 2006; Voyiadjis andAbu-Al-Rub, 2006), see Fig.3, cause anisotropy ofmetals. The anisotropy plays the fundamental role concerning damage phenomena. It is then clear, that for proper mathematical modelling of the metal behaviour one should include the anisotropy description. The frequently used isotropic simplification for metals should be thought of as a first approximation which carries not enough information for modern applications (Glema et al., 2010); though it certainly does not disavow such an approach in many applications cf. Klepaczko (2007), Rusinek and Klepaczko (2009). Fig. 3. Anisotropy of HSLA-65 steel microstructure (Narayanasamy et al., 2009) Coming back to the experimental results, being the crux of the matter of this section, showing theanisotropyofmicrodamages,wepropose the following two statements (Sumelka, 2009): (i) intrinsic microdefects are anisotropic, (ii) evolution of microdamages is directional. 978 A. Glema et al. Statement (i) confirms the experimental results that the anisotropy inme- tals caused by the intrinsic defects, comes not only from their existence but especially from their inhomogeneous structure. As an example, let us consider the effects of flat plate impact experiment in 1145 aluminium (Seaman, 1976). The separation observed is preceded by the evolution of microdamages (microvoids), consisting for the undamagedmate- rial of three stages: nucleation, growth and coalescence. Notice in Fig. that all of the microdefects are elongated perpendicularly to the impact direction, thus to themaximal tensile stresses. In this experiment, they have approxima- tely an ellipsoidal shape. So, intrinsic defects have a directional nature. Their anisotropy influences the whole deformation process, having a considerable impact on it. Fig. 4. Cracks anisotropy in 1145 aluminium after flat plate impact experiment (Seaman, 1976) For the constitutive modelling purposes, we apply in the material model the directional measure formicrodamage, since what causes the obtainedmo- del is more reliable. Notice that commonly used isotropic damage assumption (thought as ideally spherical microdamage assumption) is no longer valid for today’s sophisticated industrial requirements. With such an approximation, we lose the information on the directions of microstructure evolution. In con- sequence, one can not predict damage directions or final failure modes with satisfactory accuracy. Statement (ii) explicitly expresses the experimental fact that the aniso- tropic properties of a body of continuum evolve anisotropically during the deformation process (cf. experimental results presented byGrebe et al., 1985). Notice, that it is a consequence of the structure rearrangement itself, but espe- cially, by directional evolution of intrinsic defects. As an example in Fig.5, the Nowacki’s double shear test... 979 evolution of microvoids in the region of forming of the shear band is presen- ted. It is clearly seen that the evolution is directional, microvoids are being elongated through the shear band. So, the existing or nucleating growth of microdamages is directional according to the imposed deformation process, inducing the anisotropic evolution of material properties. Fig. 5. Anisotropic microcracks in the shear band region in Ti-6 pct Al-4 pct V alloy (after Grebe et al., 1985) Thus, for purposes of constitutive modelling, the material model should take into consideration not only the anisotropy of microdamage but also its anisotropic evolution. Such an approach is the most natural and bases on the experimental results discussed. As a concluding remark of this section, let us recall that themicrodamage evolution mechanism in metals has generally three stages: defects nucleation, their growth and coalescence. All of themare anisotropic and shouldbedescri- bed by the material model. The concept of microdamage tensor which fulfils all those requirements is presented in the next section. 3.1.2. Concept of microdamage tensor As an introductory remark for constitutive modelling, notice that thema- terial model presented is stated in terms of continuum mechanics in the fra- mework of the thermodynamical theory of viscoplasticity together with the phenomenological approach by Perzyna (2008), Glema et al. (2009), Sumelka andGlema (2009). Formally, the constitutive structure belongs to the class of simple materials with fading memory, and due to its final form and the way of incorporating the fundamental variables, belongs to materials of the rate typewith internal state variables (Truesdell andNoll, 1965). Suchan approach 980 A. Glema et al. locates themodel inmacro (meso-macro) space scale, thus all variables in the model reflect the homogenised reaction from smaller scales of observations. So, themodelling of microstructure in terms of continuummechanics cau- ses thatwe replace realmicro geometry of its structure (its anisotropic nature) by directions of tensorial field. We introduce the microdamage tensorial field of the second order (as a state variable), denoting it by ξ, cf. Perzyna (2008), Glema et al. (2009), Sumelka (2009), which reflects the experimentally obse- rved anisotropy of metals in the mathematical (constitutive) model. The concept of microdamage tensor is constructed as follows. Let us suppose that for selected points Pi in material body B, on three perpendicular planes, the ratio of the damaged area to the assumed charac- teristic area of the representative volume element (RVE) can be measured, i.e. A p i A (3.1) where A p i is the damaged area and A denotes the assumed characteristic area of the RVE, see Fig.6. Based on the calculated ratios (A p i/A), three vectors are calculated. Theirmodules are equal to those ratios and are normal toRVE planes (see Fig.6). Fig. 6. The concept of microdamage tensor Such measurements can be repeated in any different configuration obta- ined by the rotation of those three planes through point O, Fig.6. So, for every measurement configuration, from those three vectors, we compose the resultant. Next, we choose such a configuration in which the resultantmodule Nowacki’s double shear test... 981 is the largest one. Such a resultant is called themain microdamage vector and is denoted by ξ̂ (m) (Sumelka and Glema, 2007), i.e. ξ̂ (m) = A p 1 A ê1+ A p 2 A ê2+ A p 3 A ê3 (3.2) where (̂·) denotes theprincipaldirections ofmicrodamage, and Ap1 ­ A p 2 ­ A p 3. In the following step, based on the main microdamage vector, we build a vector called the microdamage vector, denoted by ξ̂ (n) (Sumelka and Glema, 2007) ξ̂ (n) = 1 ‖ξ̂(m)‖ [(Ap1 A )2 ê1+ (Ap2 A )2 ê2+ (Ap3 A )2 ê3 ] (3.3) Finally, we postulate the existence of microdamage tensorial field ξ ξ=    ξ11 ξ12 ξ13 ξ21 ξ22 ξ23 ξ31 ξ32 ξ33    (3.4) which we define in its principal directions by applying the formula combining themicrodamage vector andmicrodamage tensor (Sumelka andGlema, 2007) ξ̂ (n) = ξ̂n (3.5) where n= √ 3‖ξ̂(m)‖−1(ξ̂(m)1 ê1+ ξ̂ (m) 2 ê2+ ξ̂ (m) 3 ê3) (3.6) and the fundamental result that ξ̂= √ 3 3   ξ̂ (m) 1 0 0 0 ξ̂ (m) 2 0 0 0 ξ̂ (m) 3   (3.7) So, we arrive at the following physical interpretation of the components of microdamage tensor:The diagonal components ξii ofmicrodamage tensor ξ in its principal directions are proportional to the components of the main micro- damage vector ξ (m) i , which defines the ratio of the damaged area to the assumed characteristic area of the RVE on the plane perpendicular to the direction i. Moreover, taking the Euclidean norm from the microdamage field ξ̂, we obtain √ ξ : ξ= √ 3 3 √ (Ap1 A )2 + (Ap2 A )2 + (Ap3 A )2 (3.8) 982 A. Glema et al. Now, assuming that the characteristic length of RVE cube is l we can rewrite Eq. (3.8) as √ ξ : ξ= √ 3 l2 3 √ (A p 1) 2+(A p 2) 2+(A p 3) 2 l3 (3.9) FromEq. (3.8), another physical interpretation formicrodamage tensor appe- ars. Namely, the Euclidean norm of the microdamage field defines the scalar quantity called the volume fraction porosity or simply porosity (Perzyna 2006, 2008) √ ξ : ξ= ξ = V −Vs V = Vp V (3.10) where ξ denotes porosity (scalar damage parameter), V is the volume of a material element and Vs is the volume of the solid constituent of thatmaterial element and Vp denotes the void volume Vp = √ 3 l2 3 √ (A p 1) 2+(A p 2) 2+(A p 3) 2 (3.11) Notice that the interpretations of the microdamage tensorial field impose mathematical bounds for microdamage evolution, namely ξ ∈< 0,1 > and ξ̂ii ∈< 0,1 > (3.12) However the physical bounds are different. The experimental evidence shows that there exists the initial porosity (denoted by ξ0) which in metals is of order ξ0 ∼=10−4-10−3 (Nemes andEftis, 1991) and also, that the porosity can not reach the theoretical full saturation, i.e. ξ = 1, during the deformation process. The real maximal fracture porosity depends then on the material tested and is of order 0.09-0.35 (Dorn and Perzyna, 2002a,b, 2006). 3.2. Constitutive model for an adiabatic process The key features of the formulation presented are: (i) the description is invariant with respect to any diffeomorphism (the obtained model in cova- riant (Marsden andHughes, 1983), (ii) well-posedness of the obtained regula- rization evolution problem, (iii) rate sensitivity, (iv) finite elasto-viscoplastic deformations, (v) plastic non-normality, (vi) dissipation effects, (vii) thermo- mechanical couplings and (viii) length scale sensitivity. Below, the fundamen- tal results for an adiabatic process are given – for detailed and more general formulation see Sumelka (2009). Nowacki’s double shear test... 983 The abstract body is a differential manifold. The kinematics of the body assumes finite elasto-viscoplastic deformations which are governed by multi- plicative decomposition of the total deformation gradient to elastic and visco- plastic parts (Lee, 1969; Perzyna, 1995) F(X, t)=Fe(X, t)Fp(X, t) (3.13) In (3.13) F = ∂φ(X, t)/∂X is the deformation gradient, φ describes motion, Xdenotesmaterial coordinates, t is timeand Fe,Fp are elastic andviscoplastic parts, respectively. From the spatial deformation gradient, denoted by l l(x, t)= ∂υ(x, t) ∂x (3.14) where υ denotes the spatial velocity and x are spatial coordinates, taking its decomposition to symmetric and antisymmetric parts, we obtain l=d+w=de+we+dp+wp (3.15) d= 1 2 (l+ l⊤) w= 1 2 (l− l⊤) Now taking the Lie derivative of the assumed strain measure (the Euler- Almansi strain), we have the relation d ♭ = Lυ(e ♭) (3.16) and simultaneously d e♭ = Lυ(e e♭) dp♭ = Lυ(e p♭) (3.17) where Lυ stands for the Lie derivative and e for the Euler-Almansi strain, showing the fundamental result that the symmetric part of spatial deforma- tion gradient is directly the Lie derivative of the Euler-Almansi strain. Notice that the following constitutive structure operates in spatial configuration. The spatial covariance of themodelmentioned is just due to consequentlymaking use of the Lie derivative. Next, assuming that the balance principles hold, namely: conservation of mass, balance of momentum, balance of moment of momentum and balance of energy, and the entropy production inequality is satisfied, we define four constitutive postulates (Perzyna, 2005): 984 A. Glema et al. (i) The existence of the free energy function ψ. Formally, we apply the following form ψ = ψ̂(e,F,ϑ;µ) (3.18) where µdenotes a set of internal state variables that govern the descrip- tion of dissipation effects, and ϑ is temperature. (ii) Theaxiomof objectivity (spatial covariance). Thematerialmodel should be invariant with respect to any superposedmotion (diffeomorphism). (iii) The axiom of the entropy production. For every regular process, the constitutive functions should satisfy the second law of thermodynamics. (iv) The evolution equation for the internal state variable vector µ should be of the form Lυµ= m̂(e,F,ϑ,µ) (3.19) where evolution function m̂ has to be determined based on the experi- mental observations. The determination of m̂ appears the most challenge problem in modern con- stitutive modelling. The explicit statement of the complete set of governing equations for an adiabatic processweprecedebyadditional assumptions, seeDorn andPerzyna (2002a,b), Perzyna (2005), Glema et al. (2008), Sumelka and Glema (2009), Sumelka (2009): (i) microdamage does not considerably influence the elastic range, (ii) in every material point of the body there exists the initial microdamage state, (iii) conductivity and thermo-elastic effects are omitted. Assuming that the above holds, the body under going deformation in an adiabatic regime is governed by the following set of equations. They state the initial boundary value problem (IBVP). Find φ, υ, ρ, τ, ξ, ϑ as functions of t and position x such that (Perzyna, 1994; Łodygowski, 1996; Łodygowski and Perzyna, 1997a,b): Nowacki’s double shear test... 985 (i) the field equations φ̇ =υ υ̇= 1 ρRef [ divτ + τ ρ · gradρ− τ 1− (ξ : ξ) 1 2 grad(ξ : ξ) 1 2 ] ρ̇ =−ρdivυ+ ρ 1− (ξ : ξ) 1 2 (Lυξ : Lυξ) 1 2 (3.20) τ̇ =Le :d+2τ ·d−Lthϑ̇− (Le+gτ +τg) :dp ξ̇=2ξ ·d+ ∂g ∗ ∂τ 1 Tm 〈 Φg [ Ig τeq(ξ,ϑ,ǫp) −1 ]〉 ϑ̇ = χ∗ ρcp τ : dp+ χ∗∗ ρcp k : Lυξ (ii) the boundary conditions (a) displacement φ is prescribed on the part Γφ of Γ(B) and tractions (τ ·n)a are prescribed on the part Γτ of Γ(B), where Γφ∩Γτ =0 and Γφ∪Γτ = Γ(B) (b) heat flux q ·n=0 is prescribed on Γ(B) (iii) the initial conditions φ, υ, ρ, τ, ξ, ϑ are given for each particle X∈B at t =0, are satisfied. In the abovewedenoted: ρRef a referential density, τ –Kirchhoff stress tensor, ρ – current density, Le – elastic constitutive tensor, Lth – ther- mal operator, g –metric tensor, ∂g∗/∂τ – evolution directions for anisotropic microdamage growth processes, Tm – relaxation time, Ig – stress intensity invariant, τeq – threshold stress, χ ∗, χ∗∗ – irreversibility coefficients and cp – specific heat. For evolution problem (3.20), we assume as follows: 1. For the elastic constitutive tensor Le L e =2µI+λ(g⊗g) (3.21) where µ,λ are Lamé constants. 2. For the thermal operator Lth L th =(2µ+3λ)θg (3.22) where θ is the thermal expansion coefficient. 986 A. Glema et al. 3. For the viscoplastic flow phenomenon dp (Perzyna, 1963, 1966) d p = Λvpp (3.23) where Λvp = 1 Tm 〈 Φvp (f κ −1 )〉 = 1 Tm 〈(f κ −1 )mpl〉 f = { J′2+[n1(ϑ)+n2(ϑ)(ξ : ξ) 1 2 ]J21 }1 2 n1(ϑ)= 0 n2(ϑ)= n = const κ = {κs(ϑ)− [κs(ϑ)−κ0(ϑ)]exp[−δ(ϑ)ǫp]} [ 1− ((ξ : ξ) 1 2 ξF )β(ϑ)] (3.24) ϑ = ϑ−ϑ0 ϑ0 κs(ϑ)= κ ∗ s −κ∗∗s ϑ κ0(ϑ)= κ∗0−κ∗∗0 ϑ δ(ϑ)= δ∗− δ∗∗ϑ β(ϑ)= β∗−β∗∗ϑ p= ∂f ∂τ ∣∣∣∣ ξ=const (∥∥∥ ∂f ∂τ ∥∥∥ )−1 = 1 [2J′2+3A 2(trτ)2] 1 2 [τ ′+Atrτδ] ξF = ξF ∗− ξF∗∗ 〈(‖Lυξ‖−‖Lυξc‖ ‖Lυξc‖ )mF〉 and f denotes potential function (Shima and Oyane, 1976; Perzyna, 1986a,b;Glema et al., 2009), κ is isotropicworkhardeningsoftening func- tion (Perzyna, 1986b; Nemes and Eftis, 1993; Glema et al., 2006), τ ′ – stress deviator, J1,J ′ 2 are the first and second invariants of theKirchhoff stress tensor and the deviatoric part of the Kirchhoff stress tensor, re- spectively, A = 2(n1 +n2(ξ : ξ) 1 2), ξF ∗ can be thought as quasi-static fracture porosity and ‖Lυξc‖ denotes the equivalent critical velocity of microdamage. Notice that Eqs. (3.24)12 reflects an experimental fact that the fracture porosity changes for fast processes. Such an approach is consistent with the so called cumulative fracture criterion (Campbell, 1953, Klepaczko, 1990), namely there exists a critical time needed for saturation of microdamage to its fracture limit. 4. For themicrodamagemechanism,assumetheadditional conditions (Dor- nowski, 1999; Glema et al., 2006, 2009): • an increment in themicrodamage state is coaxial with the principal directions of the stress state, • only the positive principal stresses induce growth of the microda- mage, Nowacki’s double shear test... 987 thus we have ∂g∗ ∂τ = 〈∂ĝ ∂τ 〉∥∥∥ 〈∂ĝ ∂τ 〉∥∥∥ −1 ĝ = 1 2 τ :G : τ (3.25) Φg ( Ig τeq(ξ,ϑ,ǫp) −1 ) = ( Ig τeq −1 )mmd where τeq = c(ϑ)[1− (ξ : ξ) 1 2 ] ln 1 (ξ : ξ) 1 2 {2κs(ϑ)− [κs(ϑ)−κ0(ϑ)]F(ξ0,ξ,ϑ)} c(ϑ)= const (3.26) F = ( ξ0 1− ξ0 1− (ξ : ξ) 1 2 (ξ : ξ) 1 2 )2 3 δ + (1− (ξ : ξ) 1 2 1− ξ0 )2 3 δ and Ig = b1J1+ b2(J ′ 2) 1 2 + b3(J ′ 3) 1 3 (3.27) bi (i = 1,2,3) are material parameters, J ′ 3 is the third invariant of de- viatoric part of the Kirchhoff stress tensor. Now, taking into account the postulates formicrodamage evolution, and assuming that the tensor G can be written as a symmetric part of the fourth order unity tensor I (Łodygowski et al., 2008), i.e. G=Is Gijkl = 1 2 (δikδjl+ δilδjk) (3.28) we can write the explicit form of the growth function ĝ as ĝ = 1 2 (τ2I + τ 2 II + τ 2 III) (3.29) The gradient of ĝ with respect to the stress field gives us the following matrix representation of the tensor describing the anisotropic evolution of microdamage ∂ĝ ∂τ =   g11τI 0 0 0 g22τII 0 0 0 g33τIII   (3.30) In (3.30), τI, τII, τIII are the principal values of the Kirchhoff stress tensor. Notice that the definition of the threshold stress for the microcrack growth function τeq indicates that the growth term in the evolution func- tion for microdamage is active only after nucleation – before nucleation we have an infinite threshold limξ→0τeq =∞. 988 A. Glema et al. 5. For temperature evolution, Eqs. (3.20), we take k= τ (3.31) To conclude this Section, notice that evolution problem(3.20) iswell-posed (Łodygowski et al., 1994; Łodygowski, 1995, 1996). The relaxation time Tm can be viewed as a regularization parameter which implicitly introduces the length scale. Thus, it can be proved (Łodygowski, 1996; Glema, 2004) that the so called Cauchy problem, defined above, has a unique and stable solution. 4. Numerical analysis 4.1. Material model implementation Thesolution to the IBVPdefinedbyEqs. (3.20) hasbeenobtainedbyusing thefinite elementmethod.Abaqus/Explicit commercial finite element codehas been adapted as a solver. Themodel has been implemented in the software by taking advantage of user subroutine VUMAT, which is coupled with Abaqus system [1]. Abaqus/Explicit uses the central-difference time integration rule together with diagonal (”lumped”) element mass matrices. Some comments are needed on the stress update inVUMATuser subrouti- ne. During computations, user subroutine VUMAT controls the evolution of stresses, viscoplastic deformation, temperature andmicrodamagefields.Recall that in the presented material model, the Lie derivative has been taken into account for all rates, including the stress rate. Hence Lυτ = τ̇ − l⊤τ −τ l (4.1) while in contrast, in Abaqus/Explicit VUMAT user subroutine, the Green- Naghdi rate is calculated by default, through the following formula [1] τ(G−N) ◦ = τ̇ +τΩ−Ωτ (4.2) where Ω = Ω(G−N) = ṘR⊤ represents the angular velocity of the material (Dienes, 1979) (or spin tensor (Xiao et al., 1997)) and R denotes the rotation tensor. The important is also that the material model in Abaqus/Explicit VUMATuser subroutine is defined in a corotational coordinate system, being described by the spin tensor Ω (see Fig.7). To keep the algorithm objective in the Lie sense inVUMATuser subrouti- ne,we have applied the following approach. In the iterative procedure,we take Nowacki’s double shear test... 989 Fig. 7. Initial (XY Z) and corotational (X̃Ỹ Z̃) coordinate systems for the material derivative of the second order tensor, the forward difference scheme. Thus, for the material derivative of the Kirchhoff stress tensor, we have τ̇ ∣∣ i = τ|i+1−τ|i ∆t (4.3) UsingEqs. (4.1) and (4.2),we canwrite in the corotational coordinate system, respectively τ̃|i+1 =R⊤ ∣∣ i+1 [τ|i+∆tLυτ|i+∆t(l⊤|iτ|i+τ|i l|i)]R|i+1 (4.4) and τ̃ ∣∣ i+1 =R⊤ ∣∣ i+1 [τ|i+∆tτ(G−N) ◦ ∣∣ i +∆t(Ω|iτ|i−τ|iΩ|i)]R|i+1 (4.5) Thus we see that the Green-Naghdi rate produces an additional term ∆t(Ω|iτ|i−τ|iΩ|i) (4.6) That is why one has to subtract this term if another rate is proposed. Hence, in the presented formulation for stress update in VUMAT, we have (Sumelka and Glema, 2009) τ̃ ∣∣ i+1 =R⊤ ∣∣ i+1 [τ|i+∆t(2τ|id|i+Lυτ|i)+Υ|i]R|i+1 (4.7) where Υ|i =−∆t(Ω|iτ|i−τ|iΩ|i) and τ|i =R|iτ̃|iR⊤|i. Such an approach is needed for stresses only because other variables are kept as scalars in VUMAT user subroutine. The detailed algorithm for the whole process can be found in Sumelka (2009). 4.2. Numerical identification To solve IBVP defined by Eqs. (3.20), one has to determine 28 unknown parameters, that characterise the analysed material (steel). In Table 1, we 990 A. Glema et al. present a complete set of parameters (identified in the sense of numerical calibration) for DH-36 steel. The identification procedure bases on the results obtained experimentally byNemat-Nasser andGuo (2003). These parameters are used in the following section concerning numerical simulation of Nowacki’s double shear test. Table 1.Material parameters for DH-36 steel λ =121.154GPa µ =80.769GPa ρRef =7800kg/m 3 mmd =1 c =0.067 b1 =0 b2 =0.5 b3 =0 ξF ∗ =0.36 ξF ∗∗ =0 mF− ‖Lυξc‖− s−1 δ∗ =6.0 δ∗∗ =1.4 Tm =2.5 µs mpl =0.14 κ∗s =605MPa κ ∗∗ s =137MPa κ ∗ 0 =395MPa κ ∗∗ 0 =90MPa β∗ =11.0 β∗∗ =2.5 n1 =0 n2 =0.25 χ∗ =0.8 χ∗∗ =0.1 θ =10−6K−1 cp =470J/kgK DH-36 steel belongs to the class of high strength structural low-alloy ste- els. This steel is commonly used in shipbuilding, thus it may be subjected to high-rate loading due to collision, impact or explosion connected with a wide range of temperatures. It is then natural that DH-36 steel distingu- ishes itself due to good weldability, toughtness, good ductility and plastici- ty (true strain > 60%) under various temperatures and loading rates (see Nemat-Nasser and Guo, 2003), where a comprehensive experimental study of DH-36 over strain rates ranging from 0.001s−1 to about 800s−1 and simultaneously with initial temperatures ranging from 77K to 1000K is presented). The steel has characteristics of the bcc structure, hence it be- longs to the so called ferritic steels. Mechanical properties of DH-36 ste- el are strongly affected by impurities in its internal structure. It is impor- tant that the processing (rolling) of DH-36 steel can induce anisotropy of its structure. Figure 8 shows the adjustment of model predictions to experimental data. Notice that the numerical solution is obtained from full 3D thermomechanical analysis accounting for the anisotropic intrinsicmicrodamage processmentio- ned – in other words, the presented numerical results take into account the whole local process. The curve fitting shows that using the presentedmaterial model one can obtain numerical simulations in a very good agreement with experimental observations. Nowacki’s double shear test... 991 Fig. 8. Comparison of experimental (Nemat-Nasser andGuo, 2003) and numerical results for strain rate 3000s−1 and initial temperature 296K 4.3. Numerical simulation of Nowacki’s double shear test 4.3.1. Computational model The full spatial model of the specimen is under consideration. Geometry of the specimen is presented in Fig.2. The thickness of the specimen is equal to 0.64mm. In zone A (cf. Fig.2) at the bottom and upper planes of the specimen, the fixed boundary conditions were imposed. They model the part of Nowacki’s device that is fixed (cf. Fig.1). In zone B (cf. Fig.2) we induced rigid motion of the bottom and upper planes of the specimen (using so called rigid body constraints inAbaqus) and, simultaneously, we added additional concentrated mass (equal 0.5kg, see Klepaczko et al. (1999)) at the centre of gravity of the zone B which models the mass of movable part of Nowacki’s device and the striker. Notice that this additional mass plays ameaningful role in proper modelling of the experiment. To the rigid surfaces in zone B, an initial velocity of 19.7ms−1 (Nowak et al., 2007) and the initial temperature of 296Kwas applied. Because of lack of the experimental data concerning the initial microdamage distribution in the specimen we assumed it to be equal in every material point in the body and, simultaneously, we assumed its initial isotropy. This simplification is crucial concerning the fact that the way of mapping of the initial microdamage state has a strong impact on the final failure mode and the global answer from the specimen (Sumelka, 2009). The components of microdamage tensor were chosen in a way to obtain the initial porosity equal to 6 ·10−4, namely 992 A. Glema et al. ξ0 =   34.64 ·10−5 0 0 0 34.64 ·10−5 0 0 0 34.64 ·10−5   The contact plays an important part in the ofmodelling. It is due to strong concentrations at the corners of shear zones (the one at the pathway from the shear zone and themovable part ofNowacki’s device, and at the pathway from the shear zone and the fixed part of Nowacki’s device). At those corners, self contact occurs and has an influence in the direction of failure path. That is why the general contact in Abaqus/Explicit was applied, which assures the self contact conditions. The properties of contact were: hard normal contact (without penetration and unlimited contact stresses) and tangential contact with Coulomb’s friction model (friction coefficient equal to 0.05). The C3D8R finite elements (8-node linear brick, reduced integration ele- ment) were chosen for spatial discretisation with total 1.4Mdof for the speci- menmodel. 4.3.2. Results Let us discuss and interpret the obtained results of the Huber-Mises- Hencky (HMH) stress, the shear stress, temperature, porosity and directions of microdamage fields in the selected time instance of a process. The time up to full damage is around 180 µs. In Figs. 9 and 10, HMH and shear stresses are presented (due to the sym- metry, only a half of themodel in shown in the following figures). One should notice that bothmentioned fields indicate a homogeneous response during the test as long as the failure (loss of continuity) occurs from the corners. The important is that neither HMH stresses nor shear stresses can not be used as a failure predictor.Due to the homogeneous response in the shear zone, before discontinuity occurs, one can not estimate the failure pattern. Nevertheless, the fundamental result coming from the discussed stress fields are that the level of HMH stresses confirms the flow stresses observed in the experiment (cf. Fig.8), while the shear stress level confirms the results fromNowacki’s test that the resultant macro force acting on the specimen is around 900-1200N (Nowacki et al., 2006; Nowak et al., 2007). The temperature field versus time is presented in Fig.11. It is clearly seen that due to viscoplastic deformation in the shear zone a substantial tempera- ture rise is observed. In the strain localisation zones, the temperature reaches Nowacki’s double shear test... 993 Fig. 9. HMH stresses evolution in Nowacki’s specimen vs. time Fig. 10. Shear stresses evolution in Nowacki’s specimen vs. time 994 A. Glema et al. around 900K. It should be noticed that for the highest velocities of the striker the temperature rise can be so large that themelting point is reached in shear bands causing that no deformation is observed in the shear zone while the movable part of Nowacki’s device moves through melted steel. This creates an interesting research area for consecutive modelling of metals including ef- fects of phase transformation. Notice that the temperature field gives a better estimation for possible failure patterns (cf. Fig.11). Fig. 11. Temperature evolution in Nowacki’s specimen vs. time In Fig.12, the time history of porosity (being the normof the intrinsic ani- sotropicmicrodamagemeasure) is presented. Similarly, as for the temperature field, this field gives a better estimation of the failure pattern in comparison to the stress one. Intensified zones of microdamage evolution in the corners of the specimen are seen. The fracture that occurs between 90 µs and 120 µs starts its evolution exactly when the intensification of porosity is observed. Notice that after initiation of fracture, the information of the process from the porosity field becomes less readable. It is due to the fact that the initial porosity is of the order 6 ·10−4, while the fracture porosity is 0.36, making it impossible to obtain detailed information from the contour map where the quantities differ nearly four times. Nowacki’s double shear test... 995 Fig. 12. Porosity evolution in Nowacki’s specimen vs. time Figure 13 presents directions of the intrinsic anisotropicmicrodamage pro- cess in time. For readability, only the information on directions in selected points are shown. The arrowsmap the principal directions of themicrodama- ge tensor ξ. Due to physical interpretation given, the damage plane is the one perpendicular to the direction ofmaximal principal value presented by arrows in the plot.Notice that the orientation of softening directions changes through the specimen length showing that the problem is truly three dimensional (see also Fig.14, where the isolines of displacement field in the direction of the spe- cimen thickness are presented). Based on those principal directions, one can predict the possible damage plane shortly after the stabilisation of the process (after several viscoplastic wave reflections from boundaries), see Fig.13 where the failure pattern is marked by the dotted line. 5. Conclusions In the paper the numerical analysis of Nowacki’s double shear test in the fra- mework of the anisotropic thermo-elasto-viscoplastic model is presented. The novel is that computations cover not only elastic and viscoplastic ranges but, 996 A. Glema et al. Fig. 13. Evolution of principal directions of the microdamage tensor through the damage zone in Nowacki’s specimen vs. time Fig. 14. Evolution of thickness in Nowacki’s specimen vs. time what is fundamental, themodel describes the anisotropy of intrinsicmicroda- mages, what enables one to obtain qualitatively andquantitatively new results comparedwith the existingmodels, namely: tracing the directions of softening and predict the damage path in time. It is crucial that the theory described gives a clear physical interpretation for all assumed variables in the model. The performed analysis shows not only the possibilities of the recently developed constitutive technique, see Sumelka (2009) but, what is crucial, can Nowacki’s double shear test... 997 state a base for improving the discussed experimental technique, namely to answer the question on microdamage evolution in the process of time. This will help us to precise the definition of microdamage state, which influences the reliability of numerical analysis (Glema et al., 2010). In the presented numerical calculations, the evolution equations are purely phenomenological. Nevertheless, themodel is able to reproduce experimental observations with a very good agreement and gives a deep insight into the local thermomechanical process, what is crucial for the requirements of modern industry. References 1. Abaqus Version 6.9ef1 Theory Manual, 2010 2. Abu Al-Rub R.K., Voyiadjis G.Z., 2006, A finite strain plastic-damage model for high velocity impact using combined viscosity and gradient localiza- tion limiters: Part I – theoretical formulation, International Journal of Damage Mechanics, 15, 4, 293-334 3. Campbell J.D., 1953, The dynamic yielding of mild steel,Acta Metallurgica, 1, 6, 706-710 4. Dienes J.K., 1979, On the analysis of rotation and stress rate in deforming bodies,Acta Mechanica, 32, 217-232 5. Dornowski W., 1999, Influence of finite deformations on the growthmecha- nism of microvoids contained in structural metals,Archives of Mechanics, 51, 1, 71-86 6. Dornowski W., Perzyna P., 2002a, Analysis of the influence of various ef- fects oncycle fatiguedamage indynamicprocess,Archive ofAppliedMechanics, 72, 418-438 7. Dornowski W., Perzyna P., 2002b, Localized fracture phenomena in thermo-viscoplastic flow process under cyclic dynamic loadings,ActaMechani- ca, 155, 233-255 8. DornowskiW., PerzynaP., 2006,Numerical investigation of localized frac- ture phenomena in inelastic solids, Foundations of Civil and Environmental Engineering, 7, 79-116 9. Gary G., Nowacki W.K., 1994, Essai de cisaillement plan appliqué à des tôles minces, Journal de Physique IV, 4, C8, 65-70 10. Glema A., 2004, Analiza natury falowej zjawiska lokalizacji odkształceń pla- stycznych w ciałach stałych, PublishingHouse of PoznanUniversity of Techno- logy,Rozprawy, 379 [in Polish] 998 A. Glema et al. 11. Glema A., Łodygowski T., Perzyna P., 2008, Numerical investigation of dynamic shear bands in inelastic solids as a problem of mesomechanics, Computational Mechanics, 41, 2, 219-229 12. Glema A., Łodygowski T., Perzyna P., Sumelka W., 2006, Constituti- ve anisotropy induced by plastic strain localization, In: 35th Solid Mechanics Conference, Kraków, Poland, September 4-8, 139-140 13. Glema A., Łodygowski T., Sumelka W., 2010, Towards the modelling of an anisotropic solids,Computational Methods in Science and Technology 14. Glema A., Łodygowski T., Sumelka W., Perzyna P., 2009, The nu- merical analysis of the intrinsic anisotropic microdamage evolution in elasto- viscoplastic solids, International Journal of DamageMechanics, 18, 3, 205-231 15. Grebe H.A., Pak H.-R., Meyers M.A., 1985, Adiabatic shear localization in titanium and Ti-6 pct Al-4 pct V alloy, Metallurgical and Materials Trans- actions A, 16, 5, 761-775 16. Klepaczko J.R., 1990, Behavior of rock like materials at high strain rates in compression, International Journal of Plasticity, 6, 415-432 17. Klepaczko J.R., 1994, An experimental technique for shear testing at high and very high strain rates. The case of mild steel, International Journal of Impact Engineering, 15, 25-39 18. Klepaczko J.R., 2007, Constitutive relations in dynamic plasticity, pureme- tals and alloys,Advances in Constitutive Relations Applied in Computer Codes, CISM, Udine, Italy, July 23-27 19. Klepaczko J.R., Nguyen H.V., Nowacki W.K., 1999, Quasi-static and dynamic shearing of sheet metals,European Journal of Mechanics – A/Solids, 18, 271-289 20. Lee E.H., 1969, Elastic-plastic deformation at finite strain,ASME Journal of Applied Mechanics, 36, 1-6 21. Łodygowski T., 1995, On avoiding of spurious mesh sensitivity in numeri- cal analysis of plastic strain localization, Computer Assisted Mechanics and Engineering Sciences, 2, 231-248 22. Łodygowski T., 1996, Theoretical and Numerical Aspects of Plastic Strain Localization, PublishingHouse of PoznanUniversity of Technology,312, D.Sc. Thesis 23. ŁodygowskiT.,GlemaA., SumelkaW., 2008,Anisotropy inducedby evo- lution ofmicrostructure in ductilematerial, In: 8thWorld Congress onCompu- tational Mechanics (WCCM8), 5th European Congress on Computational Me- thods in Applied Sciences and Engineering (ECCOMAS 2008), Venice, Italy, June 30–July 5 Nowacki’s double shear test... 999 24. Łodygowski T., Perzyna P., 1997a, Localized fracture of inelastic polycry- stalline solids under dynamic loading process, International Journal Damage Mechanics, 6, 364-407 25. Łodygowski T., Perzyna P., 1997b, Numerical modelling of localized frac- ture of inelastic solids in dynamic loading process, International Journal for Numerical Methods in Engineering, 40, 4137-4158 26. Łodygowski T., Perzyna P., Lengnick M., Stein E., 1994, Viscoplastic numerical analysis of dynamic plastic shear localization for a ductile material, Archives of Mechanics, 46, 4, 541-557 27. Marsden J.E., Hughes T.J.H., 1983, Mathematical Foundations of Elasti- city, Prentice-Hall, New Jersey 28. Narayanasamy R., Parthasarathi N.L., Narayanan C.S., 2009, Effect of microstructure on void nucleation and coalescence during forming of three different HSLA steel sheets under different stress conditions, Materials and Design, 30, 1310-1324 29. Nemat-Nasser S., Guo W.-G., 2003, Thermomechanical response of DH-65 steel plates over a wide range of strain rates and temperatures, Mechanics of Materials, 35, 1023-1047 30. Nemes J.A., Eftis J., 1991, Several features of a viscoplastic study of plate- impact spallationwithmultidimensional strain,Computers and Structures,38, 3, 317-328 31. Nemes J.A., Eftis J., 1993, Constitutive modelling of the dynamic fracture of smooth tensile bars, International Journal of Plasticity, 9, 2, 243-270 32. NowackiW.K.,GadajP., LucknerJ., 2006,Effect of strain rate onductile fracture, Report IPPT,Warsaw, January 33. Nowak Z., Nowacki W.K., Perzyna P., Pęcherski R.B., 2007, Model- ling of localized fracture phenomena in high strength steels, In: 10th European Mechanics of Materials Conference EMMC10, Kazimierz Dolny, Poland, June 11-14 34. Perzyna P., 1963, The constitutive equations for rate sensitive plastic mate- rials,Quarterly of Applied Mathematics, 20, 321-332 35. Perzyna P., 1966, Fundamental problems in viscoplasticity,Advances in Ap- plied Mechanics, 9, 243-377 36. Perzyna P., 1986a, Constitutive modelling for brittle dynamic fracture in dissipative solids,Archives of Mechanics, 38, 725-738 37. Perzyna P., 1986b, Internal state variable description of dynamic fracture of ductile solids, International Journal of Solids and Structures, 22, 797-818 38. Perzyna P., 1994, Instability phenomena and adiabatic shear band localiza- tion in thermoplastic flow process,Acta Mechanica, 106, 173-205 1000 A. Glema et al. 39. Perzyna P., 1995, Interactions of elastic-viscoplastic waves and localization phenomena in solids, In: IUTAM Symposium on Nonlinear Waves in So- lids, Wegner J.L. and Norwood F.R. (Edit.), Victoria, Canada, August 15-20, 114-121 40. Perzyna P., 2005, The thermodynamical theory of elasto-viscoplasticity,En- gineering Transactions, 53, 235-316 41. Perzyna P., 2006, The thermodynamical theory of elasto-viscoplasticity ac- counting for microshear banding and induced anisotropy effects, In: 35th Solid Mechanics Conference, Kraków, Poland, September 4-8, 35-36 42. PerzynaP., 2008,The thermodynamical theory of elasto-viscoplasticity acco- unting for microshear banding and induced anisotropy effects, Mechanics, 27, 1, 25-42 43. Pęcherski R.B., Nowacki W.K., Nowak Z., Perzyna P., 2009, Effect of strain rate on ductile fracture. A new methodology, In: Workshop, Dynamic Behaviour of Materials, In Memory of Our Friend and Colleague Prof. J.R. Klepaczko, Metz, France, May 13-15, 65-73 44. Rusinek A., Klepaczko J.R., 2009, Experiments on heat generated during plastic deformation and stored energy for TRIP steels, Materials and Design, 30, 1, 35-48 45. Seaman L., CurranD.R., ShockeyD.A., 1976,Computationalmodels for ductile and brittle fracture, Journal of Applied Physics, 47, 11, 4814-4826 46. Shima S., Oyane M., 1976, Plasticity for porous solids, International Journal of Mechanical Sciences, 18, 285-291 47. Sumelka W., 2009, The Constitutive Model of the Anisotropy Evolution for Metals with Microstructural Defects, PublishingHouse of PoznanUniversity of Technology, Poznań, Poland 48. Sumelka W., Glema A., 2007, The evolution of microvoids in elastic solids, In: 17th International Conference on Computer Methods in Mechanics CMM- 2007, Łódź-Spała, Poland, June 19-22, 347-348 49. Sumelka W., Glema A., 2009, Theoretical and computational aspects of implementation of anisotropic constitutemodel formetalswithmicrostructural defects, In: 18th Int. Conf. on Computer Methods in Mechanics CMM-2009, Zielona Góra, Poland,May 18-21, 451-452 50. TruesdellC.,NollW., 1965,TheNon-Linear Field Theories ofMechanics, In: Handbuch der Physik III/3, Springer-Verlag, S. Flügge (Edit.) 51. VoyiadjisG.Z.,AbuAl-RubR.K., 2006,Afinite strain plastic-damagemo- del for high velocity impacts using combined viscosity and gradient localization limiters: Part II – Numerical aspects and simulations, International Journal of Damage Mechanics, 15, 4, 335-373 Nowacki’s double shear test... 1001 52. Xiao H., Bruhns O.T., Meyers A., 1997, Logarithmic strain, logarithmic spin and logarithmic rare,Acta Mechanica, 124, 89-105 Analiza testu podwójnego ścinania Nowackiego w ramach anizotropowego termo-sprężysto-lepkoplastycznego modelu materiału Streszczenie Wpracy przedstawiono numeryczne symulacje testu podwójnego ścinania, zapro- ponowanego przez prof. Nowackiego, w ramach sformułowania teorii lepkoplastycz- ności dla anizotropowych ciał stałych. Analizy numeryczne obejmują modele trój- wymiarowe i są wykonane dla stali DH-36 wwarunkach adiabatycznych (analiza ciał anizotropowychmożebyćprzeprowadzonawyłącznienamodelach trójwymiarowych). W trakcie analiz obserwuje się prędkości deformacji rzędu 104-107 s−1, a czas trwa- nia procesu do całkowitego zerwania próbki (utraty ciągłości w strefie lokalizacji) jest z przedziału 150-300µs. Oryginalność badańwynika z faktu uwzględnienia w de- finicji konstytutywnego modelu anizotropowego wewnętrznego procesu mikrouszko- dzeń.W rezultacie, uzyskane wyniki dają jakościowo i ilościowo nowy obraz procesu, wszczególnościumożliwiają śledzeniekierunkówosłabieniaorazdokładniejsze (bliższe rezultatom eksperymentalnym) odwzorowanie ścieżki zniszczenia. Manuscript received April 21, 2010; accepted for print June 18, 2010