Computational modelling of a multifield single-crystal gradient plasticity formulation

852 COMPUTATIONAL MODELLING OF A MULTIFIELD SINGLE-CRYSTAL GRADIENT PLASTICITY FORMULATION C. B. HIRSCHBERGER∗, B. D. REDDY† ∗Institute of Continuum Mechanics Leibniz Universit¨at Hannover Appelstr. 11, 30167 Hannover, Germany e-mail: [email protected], †Centre of Research in Applied Mechanics (CERECAM) University of Cape Town, South Africa 7701 Rondebosch, South Africa e-mail: [email protected], www.cerecam.uct.ac.za Key words: Crystal plasticity, dislocations, finite deformation, size effects, finite elements. Abstract. A model of higher-order single crystal plasticity is presented and reviewed in order to develop a corresponding finite-element framework. Contrary to the underlying model of Gurtin [Int. J. Plast. 24:702-725, 2008], here rather than the slip rate, the slip and its gradient constitute primary micro state variables. The resulting rate-dependent formulation accounts for size effects through the free energy depending on density of geometrically necessary dislocations. The relationship to multifield theories of continua with microstructure is pointed out. With the presented finite-element approach, the corresponding fully coupled initial-boundary value problem is solved monolithically, and features of the model are illustrated in two preliminary numerical examples. 1 INTRODUCTION The size-dependent behaviour of polycrystalline materials such as metals at grain sizes of the order of tens to hundreds of microns is well documented. Such behaviour stems from heterogeneities in crystallites and arise, for example, due to the existence of grain boundaries, as well as due to impurities, inclusions, and other imperfections in the crystal lattices. The appropriate modelling of behaviour at the microstructural level requires a knowledge of the underlying dynamics of dislocations, and proper incorporation of such dynamics and associated length scales into the model. Dislocation-based crystal plasticity formulations capture size effects largely by including dislocation behaviour through an 1 XI International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XI E. Oñate, D.R.J. Owen, D. Peric and B. Suárez (Eds) 853 C. B. Hirschberger, B. D. Reddy averaged field description of dislocation populations in crystals. Some representative works in an extensive literature include [17, 9, 13, 6, 10, 12, 18]. Typically, the mutual interaction of dislocations are captured in a back-stress term that counteracts the resolved shear stress driving the flow of dislocations on the lattice glide planes. The objective of this contribution is to develop finite element approximations of a model of single-crystal higher-order plasticity due to Gurtin [10]. The model is based on the use of geometrically necessary dislocation (gnd) densities as a field variable. However, in contrast to the treatment in [10], instead of the slip rate, the slip constitutes a micro state variable. The relationship between gradient of slip and gnd density is nevertheless constructed in a way that is consistent with the dissipation inequality. The slip together with macroscopic displacements are the primary unknown variables of the problem. It is shown that the formulation has a relationship to multifield theories [3, 4]. Conforming finite element approximations are based on a fully coupled, monolithic approach to solving the governing equations of the rate-dependent initial-boundary value problems. The article closes with two preliminary numerical examples. 2 A MULTIFIELD CRYSTAL PLASTICITY FRAMEWORK The model of higher-order crystal plasticity largely follows that due to Gurtin [10]. 2.1 Kinematics The starting point for the description of kinematical relations is the standard multiplicative decomposition F=Fe·Fp(1) of the deformation gradient F=∇ Xϕinto elastic and plastic parts Feand Fprespectively. Here x=ϕ(X,t) describes the motion from the material to the spatial configuration and ∇ Xis the material gradient. The plastic deformation is assumed to be isochoric implying det Fp= 1. It follows that J:= det(F) = det Fe>0. The spatial velocity gradient l=∇ xvmay be decomposed additively according to l=˙ F·F−1=le+Fe· Lp·F−1 e(2) in which  Lp=˙ FpF−1 p.(3) The elastic Cauchy–Green tensor  Ce=Ft e·Fe=F−t p·C·F−1 p(4) characterizes the deformation of the intermediate configuration, quantities in which are denoted here and henceforth by  . 2 854 C. B. Hirschberger, B. D. Reddy The motion of dislocations in a single crystal takes place on a set of nSdefined slip systems, whereby an orthonormal pair comprising a slip direction sαand slip-plane normal vector  mα(α=1,...,N) in the intermediate configuration precisely defines the α-th system. It is useful also to introduce the Schmid (projection) tensor  Zα=sα⊗ mα, which is trace-free. For edge dislocations, the slip line direction lαis defined by lα= mα×sα, so that { mα,sα,lα}form a local orthonormal basis. In case of screw dislocations, the slip line and the slip direction coincide to sα. The plastic distortion-rate tensor  Lpis determined by the slip rates acting on each of the slip systems according to  Lp= nS  α=1 ˙γαsα⊗ mα=: nS  α=1 ˙γα Zα.(5) The slip direction sα, slip plane normal  mαand dislocation line direction lαmay be mapped to their counterparts sα,mαand lαin the current configuration by setting sα=Fe·sα,mα=(Fe)−t· mα,lα=Fe·lα.(6) 2.1.1 Dislocation densities The dislocations and their interactions are accounted for by fields of spatial densities of dislocations. In the spirit of Gurtin [10] we define both the total density of dislocations, ρα , and the – polar – density of geometrically necessary dislocations (gnd), κα , per unit length, i. e.,, normalized by the Burgers vector length. According to Nye [14], only the geometrically necessary dislocations are relevant to the occurrence of size effects. Following [2], the gnd relates to the gradient of slip as ˙κα =∇ x˙γα·pα=∇ X˙γα·F−1 p· pαwith  pα=−sαfor edge dislocations (=⊥) lαfor screw dislocations (=⊙). (7) Here a pullback from the spatial form to the intermediate configuration has been carried out. The respective subscript ∈ {⊥,⊙} identifies either edge or screw dislocations. Remark 1. The formulation by Gurtin [10] bases on the slip rate ναas the primary micro state variable, the gradient of which corresponds to the gnd rate. Contrarily, we assume the slip itself and its gradient as primary micro quantities, for which we develop the governing equations, the variational formulation and computational algorithm. Our choice is particularly beneficial as it reduces the complexity of the computational treatment. Remark 2. In [7] the relationship between GND density and slip gradient is assumed to take the form κα =  ∇γα·  pα 3 855 C. B. Hirschberger, B. D. Reddy in which  ∇denotes the gradient with respect to the intermediate configuration. In the absence of the notion of a placement vector in the intermediate configuration, the nature of the gradient term is not clear. Likewise, a spatial relation of the form κα =∇ xγα·pα, which is used by some authors, does not imply nor is implied by (7). 2.2 Free energy, dissipation inequality, stresses and microstresses Dissipation inequality. The spatial form of the local dissipation inequality is D=σ:le+ nS  α=1 (ξα·∇ x˙γα+πα˙γα)−J−1˙ Ψ≥0 in Bt.(8) In the material configuration this inequality reads D0=1 2 Se:Lp v( Ce)+ nS  α=1 (ξα 0·∇ X˙γα+πα 0˙γα)−˙ Ψ≥0 in B0(9) with the Lie derivative Lp v( Ce)=Fp−t·˙ C·F−1 p=˙  Ce+2[  Ce· Lp]sym [16] and ξα 0=F−1·ξα. The free energy  Ψis assumed to depend on the elastic tensor  Ceand the set �κ of dislocation densities. This motivates the definition of a macrostress  Seand vector microstress ξα en according to  Se=2∂ Ψ ∂ Ce ,ξα en := J−1∂ Ψ ∂κα  pα.(10) By assuming that the microstress is purely energetic (that is, ξα=ξα en: see also [15]), the dissipation inequality (9) reduces to Dred = nS  α=1 πα˙γα≥0,(11) A flow rule for παwill be specified later. Macro- and micro-force balances. Balance equations are derived from a principle of virtual power [10]. The macroforce balance or equilibrium equation is expressed in terms of the first Piola-Kirchhoff stress and the referential body force 0as Div P+f0=0in B0.(12) This equation is supplemented by the boundary conditions u=upon ∂Bu tand P·N=tp 0 on ∂Bp 0, in which ∂Bu tand ∂Bp 0are non-overlapping parts that cover the boundary ∂B0. 4 856 C. B. Hirschberger, B. D. Reddy The microforce balance is given for each slip system αin the spatial configuration by div ξα−πα+σα= 0 on αin Bt(13) In the material configuration this reads Div ξα 0−πα 0+σα 0= 0 on αin B0.(14) Here σα 0is the resolved shear stress defined by σα 0= Ce· Se:Zα(15) and Div = F−1:∇ X. Dirichlet and Neumann boundary conditions are assumed to be γα=γpon ∂Bν 0and ξα 0·N=tξα 0on ∂Bξ 0. For an overview on the choice of micro-hard and micro-free boundary conditions, see e. g., [7] and references cited in this work. Remark 3. The macro- and microforce balance equations may also be derived by adopting a micromorphic or microfield continuum approach [11, 3], in which the energy functional corresponding to the incremental problem is written as the sum of the free energy  Ψand the dissipative flow potential Υ, which generates the dissipative microforce through the incremental form of the relation [15] πα=∂Υ ∂˙γα.(16) 2.3 Constitutive relations Based on an underlying hyperelastic material, the crystal plasticity constitutive framework relies on the choice of a free energy and the definition of a flow rule or dissipation function. 2.3.1 Free energy The free energy per reference volume is assumed to comprise an elastic or macro-energy  Ψmacro and a defect or micro-energy  Ψmicro:  Ψ= Ψmacro( Ce)+ nS  α=1  Ψmicro(κα ).(17) This yields the stress and energetic microstress via the definitions (10). For the sake of simplicity, latent hardening effects are neglected in the energetic part, and will instead be captured in the dissipative microforce via the flow rule . 5 857 C. B. Hirschberger, B. D. Reddy Due to the relatively small elastic deformations in crystal plasticity, it suffices to choose for the macro-energy a St. Venant–Kirchhoff relation  Ψmacro =λ 8tr2( Ce−I)+1 4µ[ Ce−I]2.(18) However, following Gurtin [10], the micro-energy is formulated for both edge and screw dislocations as  Ψmicro(κα )=1 2 Ns  α=1 (C1κα ⊥)2+C2κα ⊙2.(19) With the choices C1=µR2 8[1−ν]and C2=µR2 4, the relation of [8] is retrieved (see [7] for the relationship between the two approaches). 2.3.2 Flow rule and slip resistance Due to the physically well-established assumption that all dislocations are mobile all the time, hence we choose a viscoplastic flow rule. Various options are available: for example, πα=Sα|˙ γα| ˙γα 0m sgn(πα).(20) The slip resistance Sαis given in [10] by the evolution equation ˙ Sα= β hαβ(S)|˙γβ|(21) in which hαβ is a matrix of hardening moduli and Sdenotes the array [S1,...,S nS]tfor nSslip systems. See also [1] for further details. 3 NUMERICAL FRAMEWORK FOR CRYSTAL PLASTICITY A conforming finite element approximation is used to solve the initial-boundary value problem, with the displacement uand slips γαbeing the unknown variables. Finite-element approximations. In a standard Bubnov–Galerkin approximation, the test and trial functions of both unknowns are discretized with the same shape functions: thus uh= nu en  J=1 Nu JuJ,δuh= nu en  I=1 Nu IδuI,γh= nγ en  L=1 Nγ LγL,δγh= nγ en  K=1 Nγ KδγK,(22) with provision for different orders of interpolation of the displacement and slips. 6 858 C. B. Hirschberger, B. D. Reddy Table 1: Algorithm for the crystal plasticity material routine 1. The approximate exponential map with a backward Euler approximations gives Fn+1 p=[I−Λn+1]−1Fn pwhere Λn+1 = nS  α=1 ∆γα Zα 2. Obtain elastic right Cauchy–Green deformation tensor  Cefrom Fn+1 e=Fn+1 ·(Fn+1 p)−1 Cn+1 e=(Fn+1 e)t·Fn+1 e 3. Obtain the GND density from καn+1 =καn+[∇ X∆γα·F−1 p·sα]n+1 4. Obtain the stress, microstress and resolved shear stress from  Sn+1 e=2 ∂Ψ ∂ Cn+1 e ,ξα 0 n+1 =(Jn+1)−1(Fn+1)−1∂Ψ ∂καn+1 ,σ α 0 n+1 = Cn+1 e· Sn+1 e: Zα 3. Obtain the microforce πα 0from the evolution παn+1 =Sα|∆γα| ∆t˙γα 0m sgn πα Finite-element residual and iterative solution. A spatial discretization of the weak form of the macro and micro balance (12) and (14) yields the element residuals ruh 0I(uh)=fuh 0intI−fuh 0surI−fuh 0volI . =0(23) rγα h 0K(γαh)=fγα h 0intK−fγα h 0surK . =0(24) that have to vanish at equilibrium. The internal and external macro forces at nodes Iare assembled from the element contributions, which from the weak form of (12) are fuh intI= nel A e=1  B0 P·∇ XNu IdV, fuh surI= nel A e=1  ∂B0 t0Nu IdA, fuh volI= nel A e=1  B0 f0Nu IdV(25) Likewise, the internal and external micro forces stemming from (14) for each slip system αare determined at nodes Kas fγαh intK= nel A e=1  B0 ξα 0·∇ XNγ K+[πα 0−σα 0]Nγ KdVα=1,...,n S(26) fγαh surK= nel A e=1  ∂B0 tξ 0Nγ KdAα=1,...,n S(27) The energetic stresses Pand ξα 0are obtained from the free energy (10). The Schmid stress (15) and the dissipative micro force (20) follow from an Euler-backward time integration for the update on the plastic deformation, as summarized in Table 1. 7 859 C. B. Hirschberger, B. D. Reddy With these ingredients, the time-dependent problem is solved incrementally for the unknown uand γαusing an iterative global Newton-Raphson iterative solution procedure DuJruh 0IDγLruh 0I DuJrγαh 0KDγLrγαh 0K·∆uJ ∆γL=−ru I rγ K(28) with tangent stiffness matrices D•r◦ 0quantifying the sensitivity of the nodal residua ◦∈ {u, γα}with respect to the nodal unknowns •∈{u,γα}. 4 NUMERICAL EXAMPLES The numerical algorithm for the present crystal plasticity framework is demontrated in two benchmark-type problems. At this stage, the simulations are restricted to single slip, with a extension to multiple slip as part of future work. 4.1 Single slip in a shear layer We first study a shear layer with one slip system under an angle of θ=π/3, similar to the problem studied for example in [12]. Omitting periodic boundary conditions, a slip profile with little boundary influences is produced by choosing a relative broad shear layer. Unit values are used for material parameters in these preliminary computations. Macroscopically, a lateral displacement u1is prescribed at top and bottom in opposite directions, and homogeneous Neumann boundary conditions are prescribed on the other two sides. To mimic the dislocation distribution within the shear layer, homogeneous micro Dirichlet, i. e.,”micro-hard” boundary conditions, γ= 0, are chosen at the top and bottom boundary. On the other hand the left and right boundaries obey homogenous micro Neumann boundary conditions, tξα = 0, often referred to as ”micro free”. Neglecting the lateral boundary region, we concentrate on the central region of the boundary value problem. Within the shear layer with approximately homogeneous macro deformation, the shear problem correctly reflects the zero slip at the boundaries and a non zero slip over the height profile of the shear layer. Moreover, the gnd density is larger near the boundaries which reflects the pile-up of positive and negative dislocations against a dislocation-impenetrable boundary. 4.2 Single slip in a micro composite The typical example of a model composite presented in [5] comprises a composite material with rectangular elastic particles embedded in a plastically deforming matrix. Such a double-symmetric problem can be reduced to the symmetric unit cell shown in 2(a). The elastic inclusions are modelled by prescribing zero slip within these subdomains, and the cell is subjected to simple shear loading. While the macroscopic displacement exhibits a slightly heterogeneous deformation, the slip is, naturally, zero within the elastic inclusions and strongly heterogenous in the plastic matrix material. The slip particularly localizes horizontally just above and below 8 860 C. B. Hirschberger, B. D. Reddy the inclusion. Due to the limitation to single slip and the modelling of gliding mechanisms only, dislocations from the left bottom cannot propagate to the top right region (Fig. 2(e)). Instead dislocations pile up against the elastic inclusions with different signs, as shown in 2(f). 5 CONCLUSION In this short contribution, we present a multifield-type single crystal plasticity theory at finite strain, similar to the formulation of [10]. The presented governing equations stem from the choice of the the displacement, the plastic slip (rather than its rate [10]) and its gradient as the primary macro and micro state variables. With a thermodynamically consistent relationship between dislocation density rate and slip rate, the free energy, however, is formulated in terms of the gnd density, which is directly related to the gradient of slip. For this framework we provide a corresponding finite-element framework that employs both the displacement and the plastic slip as primary nodal degrees of freedom and hence is a strongly coupled problem. Results based on initial computations indicate that the fully coupled algorithm based on a conforming finite element framework is robust. Current work is concerned with the extension to multiple slip, and the use of alternative forms of the hardening law based on total dislocation densities. References [1] R. J. Asaro and A. Needleman. Texture development and strain hardening in rate dependent polycrystals. Acta Metall., 33:923–953, 1985. [2] M. F. Ashby. The deformation of plastically non-homogeneous materials. Phil. Mag., 21:399–424, 1970. [3] G. Capriz. Continua with latent microstructure. Arch. Rat. Mech. Anal., 90:43–56, 1985. [4] G. Capriz. Continua with Microstructure. Springer, 1989. [5] H. H. M. Cleveringa, E. van der Giessen, and A. Needleman. Comparison of discrete dislocation and continuum plasticity predictions for a composite material. Acta Mater., 45:3163–3179, 1997. [6] M. Ekh, M. Grymer, K. Runesson, and T. Svedberg. Gradient crystal plasticity as part of the computational modelling of polycrystals. Int. J. Numer. Meth. Engng, 72:197–220, 2007. [7] I. Ert¨urk, J. A. W. Dommelen, and M. G. D. Geers. Energetic dislocation interactions and thermodynamical aspects of strain gradient crystal plasticity theories. J. Mech. Phys. Solids, 57:1801–1814, 2009. 9 C. B. Hirschberger, B. D. Reddy AceFEM 0.1 Min. 0.1 Max. u1 0.75e1 0.50e1 0.25e1 0 0.249e1 0.499e1 0.750e1 (a) displacement u1 AceFEM 0.416e1 Min. 0.4165e1 Max. u2 0.15e2 0.10e2 0.51e3 0 0.516e3 0.103e2 0.154e2 (b) displacement u2 AceFEM 0.526e1 Min. 0.2688e2 Max. gam 0.45e1 0.38e1 0.32e1 0.25e1 0.19e1 0.12e1 0.64e2 (c) slip γ AceFEM 0.112 Min. 0.1127 Max. kap 0.72e1 0.48e1 0.24e1 0 0.240e1 0.481e1 0.722e1 (d) GND density κ AceFEM 0.609e1 Min. 0.1015 Max. sig11 0.104e1 0.151e1 0.197e1 0.243e1 0.289e1 0.335e1 0.381e1 (e) macro stress σ11 AceFEM 0.150 Min. 0.2677 Max. sig22 0.696e3 0.865e2 0.166e1 0.245e1 0.325e1 0.404e1 0.484e1 (f) macro stress σ22 AceFEM 0.606e4 Min. 0.1173 Max. sig12 0.422e1 0.478e1 0.533e1 0.589e1 0.645e1 0.701e1 0.757e1 (g) macro stress σ12 AceFEM 0.606e4 Min. 0.1173 Max. sig21 0.422e1 0.478e1 0.533e1 0.589e1 0.645e1 0.701e1 0.757e1 (h) macro stress σ21 Figure 1: Shear layer with single slip at θ= 60 and micro-hard boundary conditions for the slip. 1 2 3 4 5 6 X 1 2 3 4 Y (a) FE mesh and b.c. AceFEM 0 Min. 0.2 Max. u1 0.243e1 0.486e1 0.730e1 0.973e1 0.121 0.146 0.170 (b) Displacement u1 AceFEM 0.4805e3 Min. 0.7892e1 Max. E12 0.107e1 0.165e1 0.223e1 0.281e1 0.339e1 0.397e1 0.455e1 (c) Shear strain E12 AceFEM 0.784e3 Min. 0.1353e1 Max. Se12 0.865e3 0.214e2 0.342e2 0.470e2 0.599e2 0.727e2 0.855e2 (d) Elastic stress Se 12 AceFEM 0 Min. 0.4080 Max. gam 0.273e1 0.546e1 0.819e1 0.109 0.136 0.163 0.191 (e) Slip γ AceFEM 0.121 Min. 0.1219 Max. kap 0.44e1 0.29e1 0.14e1 0 0.149e1 0.298e1 0.447e1 (f) Dislocation density κ⊥ Figure 2: Micro composite unit cell with elastic inclusions and single slip at θ= 0. 10