# A Gibbs-potential-based framework for ideal plasticity of crystalline solids treated as a material flow through an adjustable crystal lattice space and its application to three-dimensional micropillar compression

AA Gibbs-potential-based framework for ideal plasticityof crystalline solids treated as a material ﬂow throughan adjustable crystal lattice space and its applicationto three-dimensional micropillar compression

Jan Kratochv´ıl a,b , Josef M´alek b , Piotr Minakowski c,d, ∗ a Czech Technical University, Faculty of Civil Engineering, Th´akurova 7, 166 29 Prague 6,Czech Republic b Charles University in Prague, Faculty of Mathematics and Physics, MathematicalInstitute, Sokolovsk´a 83, 186 75 Prague 8, Czech Republic c Heidelberg University, Interdisciplinary Center for Scientiﬁc Computing, Im NeuenheimerFeld 205, 69120 Heidelberg, Germany d University of Warsaw, Institute of Applied Mathematics and Mechanics, Banacha 2,02-097 Warsaw, Poland

Abstract

We propose an Eulerian thermodynamically compatible model for ideal plastic-ity of crystalline solids treated as a material ﬂow through an adjustable crystallattice space. The model is based on the additive splitting of the velocity gra-dient into the crystal lattice part and the plastic part. The approach extends aGibbs-potential-based formulation developed in [21] for obtaining the responsefunctions for elasto-visco-plastic crystals. The framework makes constitutiveassumptions for two scalar functions: the Gibbs potential and the rate of dissi-pation. The constitutive equations relating the stress to kinematical quantitiesis then determined using the condition that the rate of dissipation is maximalproviding that the relevant constraints are met. The proposed model is appliedto three-dimensional micropillar compression, and its features, both on the levelof modelling and computer simulations, are discussed and compared to relevantstudies.

Keywords: crystal plasticity, Gibbs potential, constitutive behaviour,micropillar compression, ﬁnite elements

1. Introduction

Severe plastic deformation experiments [30, 31] reveal that a crystalline ma-terial at yield can be seen as an anisotropic, highly viscous ﬂuid. A structural ∗ Corresponding author

Email address: [email protected] (Piotr Minakowski)

Preprint submitted to International Journal of Plasticity October 17, 2018 a r X i v : . [ c ond - m a t . m t r l - s c i ] M a y djustment of the crystal lattice space to the material ﬂow is seen as a defor-mation substructure. The ﬂow space is restricted to preferred crystallographicplanes and the directions causing anisotropy. Let us note that traditionallyﬁnite crystal plasticity is based on a Lagrangian convected coordinate represen-tation. The presented form of crystal plasticity provides a possibility to treatthe ﬂow-adjustment boundary-value problem alternatively as a ﬂuid ﬂow in theEulerian representation [14].The traditional multiplicative split of the deformation gradient F into anelastic part F e and a plastic part F p in the form F = F e F p does not have asound physical meaning. The sequential split means that the plastic action F p occurs ﬁrst, followed by the elastic action F e . However, that is not the way plas-ticity actually evolves; rather, elastic and plastic actions occur simultaneously.The formula F = F e F p is good for visualization purpose only, see [27].The additive splitting of the velocity gradient L = ∇ v , where v is thevelocity ﬁeld, into the elastic part L e and plastic part L p , i.e., L = L e + L p . (1)is far more natural and corresponds to the concept of plastic ﬂow. The de-formation gradient F can be found simply by solving the diﬀerential equation˙ F = LF , where the superposed dot denotes material time-diﬀerentiation. Incrystal plasticity, the parts L e and L p are speciﬁed in the following way: theelastic part L e is identiﬁed with the evolution of the lattice space vectors, while L p is determined constitutively; the Gibbs approach provides a suitable frame-work to achieve this, as will be shown in the next section.For modelling of crystal plasticity as a ﬂuid we adopt the Gibbs-potential-based framework proposed by Rajagopal and Shrinivasa [21]. They demon-strated the use of the Gibbs-potential-based formulation as a suitable tool fordeveloping a thermodynamically consistent model for a wide class of elastic ma-terials and viscoelastic ﬂuids that are characterized by an implicit constitutiveequation of the form f ( T , ˙ T , L ) = , (2)where T is the Cauchy stress.For the class of models characterized through (2), Rajagopal and Srinivasa[21] developed a thermodynamical framework that stems from the assumptionsthat the Gibbs potential is a function of the stress and the rate of dissipationis non-negative. In fact, the latter is strengthened by requiring that the rateof dissipation is maximum possible. The approach thus leads to models thatsatisfy the second law of thermodynamics, automatically. This thermodynam-ical framework, that is fully Eulerian (using the current state as a referencestate), yields a constitutive equation of the form (2) from specifying the con-stitutive equations for two scalar functions: the Gibbs potential and the rate The term “space” is used deliberately, as the crystal lattice is understood as a space ofpreferred positions in a crystalline phase considered regardless of particles staying or ﬂowingthrough them.

2f dissipation. The importance of choosing the Gibbs potential from the set ofthermodynamical potentials (including further the Helmholtz free energy, theinternal energy and the enthalpy) stems from the fact that it allows one to usethe stress as a primitive quantity instead of a kinematical quantity, as a measureof strain. As pointed out by Rajagopal and Srinivasa [21], from a causal pointof view, the traction (stress) is the cause, and the kinematics (motion or defor-mation) being the eﬀect, which motivates to use of the stress as the primitivequantity. Requiring that the Gibbs potential is a function of a suitably chosenrotated stress tensor, Rajagopal and Srinivasa [21] extend their framework tomodel the response of anisotropic media. They obtain rate type ﬂuid modelsfor both anisotropic elastic materials and for anisotropic viscoelastic materialsof Maxwell type.The objective of this paper is to reﬁne and extend the Gibbs framework inorder to derive models that are capable of describing ideally plastic deforma-tions of crystalline solids considered as a ﬂow of material through an adjustablecrystal lattice and that are furthermore compatible with the basic principlesof continuum thermodynamics. The reason for focusing our attention to idealplasticity, i.e. for ignoring work hardening or softening eﬀects, is that an exten-sion of the proposed model should not modify essential features of the Gibbsapproach. To incorporate work hardening or softening would mean to enrich themodel by suitable internal (history) variables governed by evolution equations(for a review of the hardening problem see e.g. [27]).An Eulerian approach to dynamic crystal plasticity has been recently pro-posed, analysed and tested by Cazacu and Ionescu [2, 4, 3]. Their goal wasto develop an Eulerian rate-dependent model that is suitable for high strainrates, large strains and rotations of the incompressible material and the crystallattice. Particular attention is devoted to the description of the kinematics ofthe crystal lattice in the Eulerian coordinates. The viscoplastic constitutivelaw and the equations for the evolution of the crystal lattice are expressed interms of quantities attached to the current conﬁguration. In applications involv-ing large deformations and high strain rates the elastic component of strain issmall with respect to the inelastic one, and it can be therefore neglected. Usingthis argument a rigid-viscoplastic approach has been adopted. Modelling hasbeen focused on in-plane deformation and on the role of inertia in the materialresponse.The novelty of the approach described herein is in ensuring that the consid-ered models are compatible with the second law of thermodynamics, thanks tothe fact that their derivation is based on the speciﬁcation of constitutive assump-tions for the Gibbs potential and the rate of dissipation, and by an application ofthe maximal rate of entropy production principle. In contrast with the commonrestriction to incompressible rigid body motions, the proposed model covers thecase of fully elastic response. We also refrain from making the assumption thatthe material is incompressible since in processes such as high pressure torsionincompressibility seems to be a restrictive approximation. Compressibility, onthe other hand, may relax non-physically high stresses in a 2-turn equal channelangular extrusion, cf. [15]. 3e further illustrate the capabilities and the eﬃciency of the proposed modelby performing numerical simulations for solving a three-dimensional micropil-lar compression problem. A material with face-centred cubic symmetry (FCCmaterial) with 12 slip systems is considered. A numerical method is brieﬂyoutlined; more details and further numerical results conﬁrming the eﬃciency ofthe numerical method will be given in a forthcoming paper.The organization of the paper is as follows. In Subsection 2.1, we ﬁrst recallthe main ideas of the Gibbs-potential-based approach. Then, in Subsection2.2, we recall basic concepts of crystal plasticity, specify the forms of boththe Gibbs potential and the rate of dissipation, and show that, together withbalance equations, these constitutive equations for two scalar quantitities suﬃcefor the derivation of the proposed elasto-visco-plastic model of crystal plasticity.The derived model is then compared with the Eulerian approach of Cazacu andIonescu [4, 3] in Section 2.3. Section 3 brieﬂy introduces a numerical methodand provides the results of numerical simulations for a micropillar compressionansatz. We conclude with a short overview of our main results in Section 4.

2. A ﬂow model of crystal plasticity

In this section we summarise the result achieved in [21]. These results willbe essential for our further consideration.We consider a body that occupies a conﬁguration Ω t at the current instant t .The current position of any particle is denoted by x , its velocity by v . Moreover,we introduce the symmetric and antisymmetric parts of the velocity gradient L = ∇ v through D := ( L + L T ) / W := ( L − L T ) / . (3)The mass density of the material is denoted by (cid:37) and the Cauchy stress by T .The governing balance equations for mass, momentum and angular momentum(in the absence of external forces) are given in their local forms, as follows:˙ (cid:37) + (cid:37) div v = 0 , (4a) (cid:37) ˙ v = div T , T = T T , (4b)where the superscript notation A T means the transposed tensor to any tensor A and the material time derivative of a scalar function z is given by ˙ z = z ,t + ∇ z · v (for a vector and tensor-valued function, the same relation is applied to eachcomponent). Furthermore, we will state the balance of energy (in the absenceof the body heat supply) in the form (cid:37) ˙ (cid:15) = T : D − div q , (5)where (cid:15) is the speciﬁc internal energy, q the heat ﬂux vector and D the sym-metric part of the velocity gradient introduced in (3). Finally, we express the4econd law of thermodynamics through (cid:37) ζ := (cid:37) ˙ η + div (cid:16) q θ (cid:17) and ζ ≥ , (6)where η is the speciﬁc entropy, θ the temperature and ζ the speciﬁc rate ofentropy production; here we tacitly assume that the entropy ﬂux is of the form q /θ . Introducing the (speciﬁc) Helmholtz free energy ψ , the Kirchhoﬀ stresstensor S and the (speciﬁc) rate of dissipation ξ through ψ := (cid:15) − θ η, S := T /(cid:37) and ξ := θ ζ, and using (5) and (6), we arrive at the equation for the rate of dissipation ξ = S : D − ˙ ψ − η ˙ θ − q · ∇ θ(cid:37)θ and ξ ≥ . (7)The starting point of the Gibbs-potential-based framework as developed in[21] is the assumption that the (speciﬁc ) Gibbs potential, denoted by G , is afunction of the temperature θ and the Kirchhoﬀ stress S , i.e., G = G ( θ, S ) or G ( t, x ) = G ( θ ( t, x ) , S ( t, x )) . (8)Stipulating further the Helmoltz free energy, the internal energy and the entropy,as functions of θ and S , through ψ ( θ, S ) = G ( θ, S ) − ∂G ( θ, S ) ∂ S : S ,(cid:15) ( θ, S ) = G ( θ, S ) − ∂G ( θ, S ) ∂ S : S − ∂G ( θ, S ) ∂θ θ,η ( θ, S ) = − ∂G ( θ, S ) ∂θ θ, (9)and inserting the ﬁrst and third of these relations into (7), we obtain ξ = S : (cid:26) D + ∂ G∂ S ˙ S + ∂ G∂θ ˙ θ (cid:27) − q · ∇ θ(cid:37)θ and ξ ≥ . (10)In what follows, we restrict ourselves to isothermal processes. Then, theequation (10) reduces to ξ = S : D + S : ∂ G∂ S ˙ S and ξ ≥ . (11)We have thus arrived at a representation of thermodynamics associated withthe speciﬁcation of the Gibbs potential (as given in (8)). There is however We suppress the use of the word speciﬁc for relevant quantities in what follows althoughthe quantities are taken per unit mass. D and S are both objective tensors, ˙ S and consequently( ∂ G ) / ( ∂ S ) ˙ S are not objective tensors.To overcome this diﬃculty and, moreover, to open the possibility to includeanisotropic responses, Rajagopal and Srinivasa (see [21]) propose to consider,instead of (8), the Gibbs potential depending on a rotated stress S = R T SR ,i.e., G = G ( S ) = G ( R T SR ) , (12)where R is any rotation tensor that is objective (which means that R , relatedto a motion χ , and R ∗ related to the motion χ ∗ , satisfy R ∗ = QR whenever χ and χ ∗ diﬀer by a rigid body rotation Q , all quantities being functions ofposition and time).Thus, assuming (12) (instead of (8)) and following the same scheme to theone that yields (11) from (8), we get ξ = S : D + S : ∂ G∂ S ˙ S = S : (cid:110) D − R A ˙ SR T (cid:111) and ξ ≥ , (13)where G ( S ) = − S : A S and A = − ( ∂ G ) / ( ∂ S ) . (14)Substituting S = R T SR into the expression on the right-hand side of (13)and using the fact that R T R = I we obtain S : (cid:110) D − R A ˙ SR T (cid:111) = S : (cid:110) D − R A R T (cid:16) ˙ S + R ˙ R T S + S ˙ RR T (cid:17) RR T (cid:111) = S : (cid:110) D − A ( ˙ S + S ˙ RR T + ( ˙ RR T ) T S ) (cid:111) = S : (cid:110) D − A ( ˙ S + S Ω − Ω S ) (cid:111) , where A is a fourth-order tensor such that A ijkl = R iA R jB R kC R lD A ABCD and Ω = ˙ RR T . Note that RR T = I implies that ˙ RR T = − R ˙ R T = − ( ˙ RR T ) T and Ω is antisymmetric. Next, introducing the notation (cid:79) S := ˙ S + S Ω − Ω S , (15)it is not diﬃcult to observe that (cid:79) S is an objective time derivative. With thisnotation, (13) takes the form ξ = S : ( D − A (cid:79) S ) and ξ ≥ , (16)where S , D and A (cid:79) S are objective and symmetric second-order tensors. In ad-dition, this step provides a simple and natural way for introducing, through A ,an anisotropic response as there is no restriction imposed by frame indiﬀerenceon the dependence of G on S . 6ajagopal and Srinivasa [21] use (16) to make, among others, the followingobservation: the dissipation rate ξ vanishes for arbitrary S if D = A (cid:79) S i.e. D = A { ˙ S + S Ω − Ω S } . (17)This constitutive equation, which is of the form (2), characterizes elastic (nondis-sipative) materials . In the next Section, the Gibbs potential (12) and the con-stitutive assumption (17) modiﬁed for elastic properties of the crystal latticespace will be used for modeling of crystal plasticity. A formulation of crystal plasticity in the framework of the Gibbs-potentialapproach is derived in two steps. First, we recall the standard description ofa crystal lattice space as well as slip systems and identify the lattice velocitygradient L e with the evolution of the lattice space. As a second step we specifythe constitutive equation for the rate of dissipation and by applying the assump-tion that the response of the material is such that it maximizes this dissipationfunction (and other relevant constraints) we obtain the constitutive equationfor S .A crystal lattice space can be characterized by three non-planar vectors a , a and a , which form a basis in R and span the so-called Bravais lattice (asan example, a face-centred cubic lattice is sketched in Fig. 1a). The reciprocaldual lattice vectors are deﬁned by a = a × a a · ( a × a ) , (18)with a and a deﬁned by cyclic permutation; note that a i · a j = δ ij . a) b) a a a s s s m Figure 1: Face-centred cubic (FCC) crystal structure is sketched in the left image a), whilean example of the slip system, i.e. the plane along which the crystal slides, is shown in greyin the right image b). As is well known such elastic response cannot be achieved via the Helmholtz free energyapproach, see [21] for further details and references.

7n crystalline solids the material ﬂow in a plastic regime is carried by amotion of dislocations. Crystal plasticity theory is based on the observation thatthe motion of dislocations takes place in preferred crystallographic directions(so-called slip directions) on preferred crystallographic planes (so-called slipplanes); in our consideration climb of dislocations is excluded. The unit vectorin the slip direction s ( α ) and the unit normal to the slip plane m ( α ) form the α -slip system. The set of N potentially active slip systems is one of the basicmicroscopic characteristics of crystal plasticity, hence α = 1 , . . . , N . The vectors s ( α ) and m ( α ) of the slip system α can be expressed through the lattice vectors a , a and a as s ( α ) = (cid:88) i =1 s αi a i , m ( α ) = (cid:88) i =1 m αi a i , (19)where s αi and m αi are constant coeﬃcients, called Miller indexes. For example,for s = 1, s = 1 and s = 0, and m = 1, m = 1 and m = 1 the slip directionvector s = a + a and the slip plane normal m (1) = a + a + a , expressed ina standard way through the Miller indexes, read [110] and (111), respectively;[110](111) is a typical slip system of face-centred cubic crystals as seen in Fig.1b). (Negative values are denoted by bars above the (positive) numbers.)In a deformation process the crystal lattice space adjusts to the stress andthe material ﬂow. In the proposed model, the crystal lattice is treated as asolid described by the lattice deformation gradient F e from the natural (latticereference) conﬁguration to the current conﬁguration. According to Srinivasaand Srinivasan [27, Section 13.4, p. 467] this evolution equation for the latticevectors a i is characterised in the following way L e = (cid:88) i =1 ˙ a i ⊗ a i ⇐⇒ ˙ a j = L e a j . (22)Indeed, to see that (22) holds, we ﬁrst multiply (cid:80) i =1 ˙ a i ⊗ a i by a j ( a i · a j = δ ij )and observe that L e = (cid:80) i =1 ˙ a i ⊗ a i ⇒ ˙ a j = L e a j . On the other hand, theopposite implication, i.e. L e = (cid:80) i =1 ˙ a i ⊗ a i ⇐ ˙ a j = L e a j , follows fromthe deﬁnition of the reciprocal basis (18). Note that the lattice vectors remain The concept of the natural conﬁguration has been developed and applied in a number ofareas beyond the area of plasticity by Rajagopal and his co-workers, see for example [18, 19,20]. In our setting, the natural conﬁguration is the lattice reference conﬁguration and thus a i ( x , t ) = F e ( x , t ) a i ( x , , (20)where a i ( x ,

0) are the lattice vectors in the natural (lattice reference) conﬁguration.Since the lattice velocity gradient L e is supposed to be linked with the rate of F e throughthe relation L e = ˙ F e F e − (that is analogous to the kinematical relationship between L and˙ F ), the rate form of (20) reads ˙ a i ( x , t ) = L e a i ( x , t ) . (21) F e .Our basic kinematical assumption is represented by the additive splitting ofthe velocity gradient in the form (1), i.e. L = L e + L p , (23)where L e is linked with the a i via the equations (22) and the form of L p will bespeciﬁed by maximization of the rate of dissipation based on the identiﬁcationof the constitutive equation for the rate of dissipation.Following the decomposition (23), we set D e := ( L e + L e T ) / W e := ( L e − L e T ) / , (24) D p := ( L p + L p T ) / W p := ( L p − L p T ) / , (25)and from (3) we have that D = D e + D p and W = W e + W p . (26)Considering (12) and (16) with the constitutive equation (17) modiﬁed tothe crystal lattice space D e = A (cid:79) S we get from (22) that D e = A (cid:79) S = sym (cid:32) (cid:88) i =1 ˙ a i ⊗ a i (cid:33) and L e = (cid:88) i =1 ˙ a i ⊗ a i . (27)Next, we employ the procedure of maximization of the rate of dissipation ξ to get D p ; this, together with the second equation in (22), also allows us todetermine uniquely the rate of material ﬂow L p for a crystalline solid deformedby slip. To do so, we need to specify the constitutive equation for the rate ofdissipation ξ . The driving force of slip is the resolved shear stress τ ( α ) controlledby the critical value τ ( α ) c ([27]), τ ( α ) = s ( α ) · S m ( α ) . (28)For the sake of the following procedure the critical values τ ( α ) c are assumed to bea given positive material parameter in agreement with our restriction to idealplasticity. If τ ( α ) < τ ( α ) c the slip rate and the corresponding dissipation rate ofthe α -slip system is very small. On the other hand, if τ ( α ) > τ ( α ) c the dissipationrate is high. One of the suitable possibilities is to assume that the dissipationrate is a sum of contributions of the slip systems proportional to a power of thequotients τ ( α ) /τ ( α ) c , i.e.,ˆ ξ ( S ) = N (cid:88) α =1 c (cid:18) | τ ( α ) | τ ( α ) c (cid:19) β +1 = N (cid:88) α =1 c (cid:18) | s ( α ) · S m ( α ) | τ ( α ) c (cid:19) β +1 , (29) We decompose the lattice deformation gradient F e into the elastic stretch U e and therotation R e , F e = R e U e . The angle is preserved by the rotation R e , and the loss oforthonormality is caused by the elastic stretching. c > β > ν ( α ) := c sgn( τ ( α ) ) τ ( α ) c (cid:18) | τ ( α ) | τ ( α ) c (cid:19) β , (30)we observe that ∂ ˆ ξ ( S ) ∂ S = ( β + 1) N (cid:88) α =1 ν ( α ) sym( s ( α ) ⊗ m ( α ) ) (31)and consequently ∂ ˆ ξ ( S ) ∂ S : S = ( β + 1) ˆ ξ ( S ) . (32)In order to determine the constitutive equation for the evolution of S , wemaximize ˆ ξ given in (29) with respect to S under the constraint (16) that canbe rewritten in a more compact form, namely ξ = S : { D − D e } = S : D p and ξ ≥ . (33)Hence, setting L ( S ) := ˆ ξ ( S ) + (cid:96) ( ˆ ξ ( S ) − S : D p ) , the constrained optimality condition reads0 = ∂L ( S ) ∂ S = (1 + (cid:96) ) ∂ ˆ ξ ( S ) ∂ S − (cid:96) D p ⇐⇒ D p = 1 + (cid:96)(cid:96) ∂ ˆ ξ ( S ) ∂ S . (34)Multiplying the resulting equation scalarly by S , and using (33) and (32), weconclude that 1 + (cid:96)(cid:96) = 1 β + 1 . Consequently, the second condition in (34) together with (27) and (31) implythat D − A (cid:79) S = D p = N (cid:88) α =1 ν ( α ) sym( s ( α ) ⊗ m ( α ) ) . (35)Note that ν ( α ) can be interpreted as a slip rate of the α -slip system. Thus, wehave identiﬁed D p . Since (22) and (24) imply that W = W e + W p uniquelydetermines W p , we stipulate that L p = N (cid:88) α =1 ν ( α ) ( s ( α ) ⊗ m ( α ) ) . (36)Note that this formula is consistent with the results obtained so far for thesymmetric parts D , D e and D p since we have (see (35)) D = D e + D p = A (cid:79) S + N (cid:88) α =1 ν ( α ) sym( s ( α ) ⊗ m ( α ) ) . (37)10n addition, the choice (36) and the decomposition rule (23) lead to ∇ v = (cid:88) i =1 ˙ a i ⊗ a i + N (cid:88) α =1 ν ( α ) ( s ( α ) ⊗ m ( α ) ) , (38)which we view (and use) as the evolutionary equation for the vectors a i .Recalling deﬁnition (15) and combining it with the obtained constitutiverelation (35) we obtain D − A { ˙ S + S Ω − Ω S } = N (cid:88) α =1 ν ( α ) sym( s ( α ) ⊗ m ( α ) ) , where Ω = ˙ RR T for any rotation tensor R that is objective. By choosing dif-ferent rotations R , one ends up with diﬀerent objective time derivatives, and,consequently, with diﬀerent models, with the restriction for Ω being antisym-metric. For example, by choosing R as the rotational part R e of the polardecomposition of the lattice deformation gradient F e = R e U e will yield theGreen-McInnis-Naghdi stress rate. The second example, namely the Zaremba-Jaumann derivative, extensively used in the crystal plasticity literature, resultsfrom setting Ω = ˙ RR T = W e .From now on we focus on the Zaremba-Jaumann rate (cid:79) S = ˙ S + SW e − W e S . (39)At this point we compare the constitutive relation (35) with relations thatare typically used in crystal plasticity, see e.g. [16, 17]. In these studies, basedon the Helmholtz potential approach, the resulting constitutive equation takesthe form ˙ S − W e S − D e S + SW e − SD e = C ( D e ) , (40)where the rate derivative on the left-hand side is the Oldroyd derivative. Toobtain the Zaremba-Jaumann rate, the terms D e S + SD e are often neglectedor transferred to the right-hand side of (40) by deﬁning the fourth-order tensor K ( D e ) = C ( D e ) + D e S + SD e . In this sense the Gibbs and the Helmholtzpotential based approaches are comparable. In summary, the primary variables of the proposed crystal plasticity modelbased on the Gibbs potential are: the density (cid:37) , the velocity v , the Kirchhoﬀstress S , and the vectors a , a and a being functions of the current position x and time t . These variables are governed by the system of equations (4a),(4b), (35) and (38), i.e., 11 (cid:37) + (cid:37) div v = 0 , (41a) (cid:37) ˙ v = div ( (cid:37) S ) , (41b) A (cid:79) S = (cid:88) i =1 sym( ˙ a i ⊗ a i ) , (41c) (cid:88) i =1 ˙ a i ⊗ a i = ∇ v − N (cid:88) α =1 ν ( α ) ( s ( α ) ⊗ m ( α ) ) . (41d)The slip system vectors s ( α ) and m ( α ) in (41d) are expressed in terms ofthe lattice vectors a , a and a through (18) and (19); the slip rates ν ( α ) aregoverned by the equations (30), (28).We also multiply (41c) by C := A − and obtain (cid:79) S = C (cid:88) i =1 sym( ˙ a i ⊗ a i ) = C D e . (42)The elastic tensor C is anisotropic; in what follows we assume that C is of cubicsymmetry characterized by three independent elements, c , c , and c , whichcorrespond to C , C , and C respectively, through the standard Voigtnotation.In comparison with (41) the primary variables of the rigid-plastic modelproposed by Cazacu and Ionescu [2, 3, 4] are: the velocity v , which is divergence-free, the Cauchy stress T = (cid:37) S , and the slip system vectors s ( α ) and m ( α ) , α = 1 , . . . , N as functions of the current position x and time t . These variablesare governed by the system of the equations:div v = 0 , (43a) (cid:37) ˙ v = div T , (43b) D = N (cid:88) α =1 ν ( α ) sym( s ( α ) ⊗ m ( α ) ) , (43c)˙ s ( α ) = ( W e − N (cid:88) α =1 ν ( α ) skew( s ( α ) ⊗ m ( α ) )) s ( α ) , (43d)˙ m ( α ) = ( W e − N (cid:88) α =1 ν ( α ) skew( s ( α ) ⊗ m ( α ) )) m ( α ) . (43e)The expression in the brackets on the right-hand side of (43d) and (43e) is therate of the lattice rotation ˙ R e R − e (for a rigid lattice F e = R e ); (43d) and(43e) are the rigid plasticity versions of (41d) of our model evaluated directlyin terms of the vectors of the slip system s ( α ) and m ( α ) , α = 1 , . . . , N insteadof the lattice vectors a i , i = 1 , , N can be large (for e.g.for FCC crystal structure N = 12), working with the lattice vectors a i makes12he system (41) advantageous with respect to the system (43). It reduces thenumber of unknowns in the system and thus plays an important role in theeﬃciency of three-dimensional numerical simulations.Due to the assumed rigidity, the density (cid:37) in (43b) is a material constant.The equation (43c) is an implicit constitutive equation for the stress T as couldbe seen if one used for the slip rates ν ( α ) the power law (30) and for the resolvedshear stress τ ( α ) = s ( α ) · T m ( α ) . Note that in this case all N slip systems areactive, and the power law determines their relative activity. However, Cazacuand Ionescu preferred, mainly from the computational point of view, to employinstead the power law (30) the visco-plastic extension of the classical rigid-plasticSchmid law using the overstress approach ν ( α ) = 1 η ( α ) [ | τ ( α ) | − τ ( α ) c ] + sign( τ ( α ) ) , (44)where η α is the viscosity; [ z ] + = ( z + | z | ) / z . Using this overstress approach, it is important to mention that onlyup to 5 resolved shear stresses τ ( α ) can be independent (there are 5 componentsof the stress deviator controlling slip). The slip rates ν ( α ) given by (44) are alsonot independent as they have to satisfy the kinematic constrain (43c). Giventhe deformation rate D , the slip rate ν ( α ) of the active slip systems can bedetermined by minimizing the internal power under the constrain (43c). Thedeviatoric part of the stress T corresponding to a given deformation rate D isobtained as the Lagrange multiplier of this minimization problem.One of the main diﬀerences is that in the rigid-plastic model (43) the elasticproperties of the lattice are neglected as they are assumed to be small in com-parison with the large plastic deformations. This assumption is reasonable inmodelling e.g. metal production processes and generally it is suitable from amacroscopic point of view. If the interest is focused on a meso-scale, where thedeformation is spontaneously heterogeneous, not only the local misorientationsof the crystal lattice but also lattice strains connected with a stress redistri-bution may play a role. It is worth mentioning that classical crystal plasticitymodels of shear bands consider anisotropic elasticity as an important ingredi-ent. Moreover, elasticity and compressibility have a stabilizing eﬀect, helpingthe computational strategy as indicated in Section 3, see also [15].In some aspects the models (41) and (43) are similar. Both use the ratedecomposition rule (23) in the current conﬁguration as the basic kinematic as-sumption (the reference to decomposition rule (23) in Cazacu and Ionescu [4] isonly indicative) and their Eulerian descriptions are conceptually equivalent (ex-cept the acceptance of elasticity and the diﬀerence in the power vs. overstressrules mentioned above). A substantial diﬀerence is in the formulations of themodels. The rigid-plastic model (43) is based on the slip rate composition of theplastic velocity gradient (36) and the resolved shear stress (28) as the primaryassumptions. On the other hand, the proposed constitutive modelling basedon the Gibbs-potential formulation is constructed from two scalar functions ofthe Kirchhoﬀ stress S : the Gibbs potential G ( S ) given by (14), S = R Te SR e ,and the rate of dissipation ξ ( S ) given by (29). The form of the material ﬂow13radient L p given by (36) is the consequence of the dissipation rate maximiza-tion. Moreover, the Gibbs-potential framework guarantees that the model isthermodynamically consistent.

3. Numerical example: three-dimensional micropillar compression

The algorithm presented in this section is based on the Eulerian frameworkdescribed above and thus contributes to the Crystal Plasticity Finite ElementMethod (CPFEM). A recent overview by Roters et al. [23, 22] summarizesapplications of the CPFEM. Most of the methods that have been developed sofar adopt a Lagrangian description of the continuum problem. In the context ofthe present paper we speciﬁcally highlight the following references focused onthe application of CPFEM to micropillar compression: [8, 9, 12, 10, 26].In [10], Jung et al. compare diﬀerent primary slip plane inclination anglesand use a hardening rule that accounts for anisotropic slip system interactions.Currently, it seems that the main interest in the ﬁeld concentrates on ﬁnite el-ement analysis of geometrically necessary dislocations, see e.g. [8, 9, 12]. Sincethe primary purpose of this study is to computationally analyze the Eulerianapproach based on the Gibbs-potential-based formulation, we focus on the caseof ideal plasticity here and postpone the study of the hardening (and possiblynon-local) eﬀects to our further research. Note, however that in Minakowskiet al. [16] the authors presented a crystal plasticity model including isotropichardening and proposed an Arbitrary Lagrangian Eulerian (ALE) based numer-ical method for a two-dimensional plastic ﬂow of a single crystal compressed ina channel die.To examine the potential of the derived model we consider three-dimensionalmicropillar compression. For comparison we use the same setup as Kuroda in[12], where the behaviour of single crystal micropillars subjected to compressiveloading is described. A three-dimensional ﬁnite element method incorporatinghigher-order gradient crystal plasticity is analyzed. The author ﬁrst simulateda compression test for a model without higher strain gradients, which we usefor comparison.The compression testing methodology is shown schematically in Figure 2(like Kuroda [12] we consider samples with a conical base). The setting mimicsthe compression experiment commonly performed on macroscopic samples [28,29, 25]. The compressive stress is deﬁned as the sum of the nodal forces in the − x direction at the top surface divided by the initial cross-section area. Thenominal compressive strain is the ratio of the displacement of the top surface inthe − x direction over the original height of the sample.We use the same set of values of parameters of the compressed materialas in [12], which corresponds to the experimental and numerical works [5, 25].The ratio between the height L and the diameter D is selected to be 2 . C are the same ones as for the simplest anisotropic solid, the three values c , c and c are independent. The elastic constants are chosen to be c =1480 . c = 188 . c = 133 . l = 10 µ m , the velocity v = 1 µm/s , thereference slip rate c = 0 . s , the critical resolved shear stress τ c = 420MPa.The rate sensitivity parameter m is taken to be 0 .

05, for β = m , compare(30). A review of computational strategies in an eﬀort to overcome numericalinstabilities connected with the power law in (30) can be found e.g. in [6, 11].In accordance with the assumptions of our approach, the hardening eﬀects areneglected. These material parameters are those for a single-crystal Ni-basesuperalloy, which has the nominal composition of NiCoCrTa.2Al WReMo.2Hf,see [25].We solve the set of equations derived in Section 2, supplemented with properboundary and initial conditions, see also (41). Let T >

0, Ω ⊂ R open andbounded, (cid:37) : Ω → (cid:60) , v : Ω → R , S : Ω → (cid:60) × , a i : Ω → R , i = 1 , , v D : Γ D → R , t : Γ N → R . The system is satisﬁed inside the domain Ω ⊂ R ,the boundary ∂ Ω consists of two non-intersecting parts Γ D and Γ N correspond-ing to the Dirichlet and the Neumann boundary conditions, respectively. Welook for the density (cid:37) , the velocity v , the Kirchhoﬀ stress S , and the latticevectors a , a , a ( s ( α ) and m ( α ) are given through lattice vectors, equations(18) and (19)) satisfying˙ (cid:37) + (cid:37) div v = 0 , (45a) (cid:37) ˙ v = div ( (cid:37) S ) , (45b) (cid:79) S = C (cid:32) D − N (cid:88) α =1 ν ( α ) sym (cid:16) s ( α ) ⊗ m ( α ) (cid:17)(cid:33) , (45c) (cid:88) i =1 ˙ a i ⊗ a i = ∇ v − N (cid:88) α =1 ν ( α ) (cid:16) s ( α ) ⊗ m ( α ) (cid:17) , (45d) (cid:37) ( x,

0) = (cid:37) ( x ) ∀ x ∈ Ω , (45e) v ( x,

0) = v ( x ) ∀ x ∈ Ω , (45f) S ( x,

0) = S ( x ) ∀ x ∈ Ω , (45g) a i ( x,

0) = a i ( x ) ∀ x ∈ Ω , (45h) v = v D on Γ D × (0 , T ) , (45i) Sn = t on Γ N × (0 , T ) , (45j)where (cid:37) ( x ) = 1, S ( x ) n = for x ∈ Γ N × [0 , ∞ ), S ( x ) = for x ∈ Ω, a ( x ) = (1 , , a ( x ) = (0 , , a ( x ) = (0 , ,

1) and n is the outernormal vector. Moreover, as already explained the slip rates ν ( α ) are governedby the equations (30), (28).Two kinds of top velocity boundary condition are considered, see Figure 2(b-c). In the ﬁrst case, corresponding to Kuroda’s condition BC1 [12], we ﬁx15 D L ( a ) ( b ) ( c )Γ Γ Γ Figure 2: Scheme of micropillar compression (a), compressed sample with constrained (b) andunconstrained (c) boundary conditions. the velocity on the top boundary Γ and require v ( x, t ) | Γ = (0 , , − . (46)This condition is referred below to as the constrained case and the setting mimicsthe situation with a very stiﬀ loading axis and the friction between the topsurface and the plunger is extremely high. The second case (correspondingto Kuroda’s condition BC2 [12]) is called below the unconstrained case and itallows the top surface to move in the plane perpendicular to the z -axis, whichmeans that we require that v · n = v z ( x, t ) | Γ = − Sn · τ i = 0 , i = 1 , , (47)where τ , τ are two vectors generating the basis of the tangent plane orthogonalto n at the boundary. This simply mimics the case where there is no frictionbetween the top surface and the plunger. In the rest of the Dirichlet boundary,namely the bottom of the specimen, the velocity is ﬁxed v ( x, t ) | Γ = (0 , , S ( x, t ) n = on the lateralsurface Γ .The crystal orientation is relative to the ﬁxed sample coordinates in thereference conﬁguration. The [¯123] direction is set to be parallel to the sampleaxis, i.e. the X direction, the [1¯11] direction is chosen to coincide with the X direction, and the X direction coincides with [54¯1] (the overbar denotes thenegative value). The speciﬁcation of the slip systems is given in Table 1. Theprimary slip system is number 2, (111)[10¯1], the angle between s (2) and X − X plane reads 49 . ◦ . 16 √ s ( α ) √ m ( α ) α √ s ( α ) √ m ( α ) α √ s ( α ) √ m ( α ) Table 1: The slip directions and normals for a FCC crystal in Miller notation.

In the solution of the micropillar boundary-value problem without higher-order gradients, Kuroda [12] used a generally adopted ﬁnite element Lagrangianframework. In his numerical experiments, he evaluated deformed ﬁnite elementmeshes, contour maps of slip, maps of lattice rotations, and stress-strain di-agrams for nominal compressive strain up to 0 .

2. In the present example, weemploy the Arbitrary Lagrangian Eulerian (ALE) approach in the sense that weuse an Eulerian formulation on a moving mesh, which captures the free bound-ary of the micropillar. In this way, we are able to extend the tested range ofnominal compressive strain up to 0 .

5. For a brief description of our numericalscheme we refer to the Appendix. The detailed description concerning accuracyand eﬃciency of the numerical method and further numerical results of the ap-plication of the ALE method to three-dimensional crystal plasticity problemswill be described in detail in a forthcoming work.( a )( b ) A B A B

Figure 3: The velocity ﬁeld (on the sample surface) for X − X and X − X plane views(denoted by A and B, respectively) with constrained (left) and unconstrained (right) boundaryconditions for nominal strains 0.2 (a) and 0.4 (b). a )( b )( c ) A B A B

Figure 4: Contours of slip, on the specimen surface, on the primary system ( s (2) ) for X − X and X − X plane views (denoted by A and B, respectively) with constrained (left) andunconstrained (right) boundary conditions for nominal strains 0.2 (a). 0.3 (b) and 0.4 (c). Since the computations are conducted in the Eulerian coordinates our pri-mary variable is the velocity. In Figure 3 we present the velocity ﬁeld at variousstrain levels in the whole sample for the X − X and X − X plane views.Due to anisotropy of the acting slip systems the deformation proﬁles are notthe same for diﬀerent plane views, Fig. 3 of [12] exhibits the same features.We focus on the eﬀect of the applied boundary conditions in comparison with ascanning electron microscope compression experiments, see [25] and the numer-ical studies [12, 25]. The presented approach allows for nominal strains up to0 . .

2. Nevertheless, for nominalstrains higher than 0 . .3. Activity of slip systems In Figure 4 we present the contour maps of slip on the primary system ( s (2) )and compare the results for two kinds of boundary conditions. The sample de-formation and the primary slip activity exhibit the same features as the oneshown in Figure 3 of [12], in the nominal strain region up to 0.2. There isa measurable diﬀerence between the alignment of slip activity regions (shearbands) between diﬀerent types of boundary conditions, see Figure 4. The ma-jority of deformation has been localized to an intense slip band. Slips are highlyconcentrated in a shear zone.In Figure 5 we present the average slip rate for each of twelve individualslip systems, and compare two kinds of boundary conditions. The slip rates arenot qualitatively similar for diﬀerent types of boundary conditions. While theconstrained boundary conditions give rise to stable constant-in-time values, theunconstrained case changes the behaviour for strains over 0 .

05. In both caseswe observe instabilities for high strains over 0.3, see Figure 4. The inactive slipsystems α ∈ { , , } are oriented in such a way that the slip direction s ( α ) isperpendicular to the compression axis, namely m ( α ) · ( − , ,

3) = 0, see Table 1. − . − . − . . . . . ν (1) ν (2) ν (3) ν (4) ν (5) ν (6) ν (7) ν (8) ν (9) ν (10) ν (11) ν (12) − . − . − . . . . . ν (1) ν (2) ν (3) ν (4) ν (5) ν (6) ν (7) ν (8) ν (9) ν (10) ν (11) ν (12) s li p r a t e s ( a ) s li p r a t e s ( b ) strain [%] Figure 5: The average slip rate for individual slip systems for a constrained (a) and anunconstrained domain (b).

Figure 6 shows the nominal compressive stress-nominal compressive straincurves. Despite of the considered zero hardening, the constrained condition ex-hibited signiﬁcant stress increase while the unconstrained condition showed a19 .00 0.05 0.10 0.15 0.20 0.25 0.3002004006008001000 S − BC S − BC N o m i n a l c o m p r e ss i v e s t r e ss [ M P a ] Nominal compressive strain [%]

Figure 6: The nominal compressive stress for a constrained ( BC

1) and an unconstraineddomain ( BC decrease of the calculated nominal compressive stress. The stress increase inthe constrained case is caused by the Dirichlet velocity boundary condition onthe top and bottom boundaries. There, plastic deformation is suppressed, andpart of the sample deformed elastically causing the stress to increase. The eﬀectis associated with barrelling of the sample. On the other hand, the decrease ofthe nominal stress with the increasing nominal strain in the case of the uncon-strained boundary condition is associated with the inclination of the sample. Asis seen from the Figure 6 the results in the strain range up to 0 .

4. Summary • We presented a new approach to the derivation of a thermodynamicallycompatible rate-type ﬂuid model for crystal plastic materials. We em-ployed a Gibbs-potential-based approach and derived a rate-type stress-strain constitutive equation. Instead of the standard crystal plasticityhypothesis that the plastic part of the velocity gradient equals a sum ofcontributions of the individual slip systems L p = (cid:80) Nα =1 ν ( α ) s ( α ) ⊗ m ( α ) we have employed the procedure of maximization of the rate of dissipationto derive this relation from a single scalar function, (29). It is worth em-phasizing that the Gibbs-potential-based approach, as developed in [21]and presented here, automatically provides an objective rate derivative20or the stress as well as the possibility to incorporate anisotropy of thematerial response. • We have incorporated an Eulerian form for the evolution of the lattice basisvectors instead of computing directly the evolution of the slip directionsand the normals to the slip plane. This is one of the novel ingredients ofthe developed model. It allows us to perform eﬃcient three-dimensionalcomputations of the material with a full set of slip systems (e.g. 12 slipsystems in the case of FCC structure), as illustrated on two test problemsconcerning three-dimensional micropillar compression. • We have formulated a ﬁnite element discretization scheme, which has beenused to solve the fully coupled problem. Our solver is monolithic. TheArbitrary Lagrangian Eulerian (ALE) approach has been applied in orderto use the Eulerian formulation on a moving mesh, that captures thefree boundary. We developed a new algorithm for a three-dimensionalcompression and achieved qualitative correspondence between our resultsand the experimental and numerical results concerning the same problem,as presented in [12]. • The fact that our method is based on the Eulerian formulation, and con-sequently the primal role is given to the velocity of the material and thedistortion of the crystal lattice space, enabled us to demonstrate that thewhole approach is capable of simulating deformation processes with largestrains for suitably deﬁned problems.

Appendix

The ALE approach is a ﬁnite element formulation in which the computa-tional system is not a priori ﬁxed in space (Eulerian) or attached to the mate-rial (Lagrangian) [7, 24]. In fact we can formulate the problem on an arbitrarymoving, time-dependent domain and deﬁne the ALE mapping between a ﬁxedcomputational domain and an arbitrary domain. The ALE mapping is the keyto transforming our system to a ﬁxed computational domain. To this end theproper deﬁnition of the ALE mapping is of a great importance. The class ofALE methods is widely applied, e.g. to ﬂuid-structure interaction problems. Incrystal plasticity, the ALE method has been successfully tested in the solutionof a two-dimensional problem of a channel-die compression in [2, 4, 3], wherethe nominal compressive strain of 0 .

45 has been reached.In the present paper, the ALE method is used to solve the system (45).The implementation was done using FEniCS [1, 13]. A mixed ﬁnite elementdiscretization is used in space, and time is discretized by the backward Eulerscheme. The choice of the conforming simplical tetrahedral elements for eachvariable consists of P elements for the velocity, P for the density and thelattice basis vectors, and P − disc (discontinuous elements) for the Kirchhoﬀstress. We employ the ALE approach in the sense that we use our Eulerianformulation on a moving mesh, which captures the free boundary.21o solve (45) at each time step, we perform two actions. Given v k , (cid:37) k , S k , a i k from the previous time level we proceed as follows: Step 1:

Solve problem for v k +1 h , (cid:37) k +1 h , S k +1 h , and a i k +1 h ( i ∈ { , , } ) (cid:32) (cid:37) k +1 h − (cid:37) kh ∆ t + div ( (cid:37) k +1 h v k +1 h ) , z (cid:37) (cid:33) = 0 ,(cid:37) k +1 h (cid:32) v k +1 h − v kh ∆ t + (cid:0) v k +1 h − v kh − mesh (cid:1) ∇ v k +1 h , z v (cid:33) + (cid:16) (cid:37) k +1 h S k +1 h , ∇ z v (cid:17) = 0 , (cid:32) (cid:79) S k +1 h − C D e k +1 h , Z S (cid:33) = 0 , (cid:32) a i k +1 h − a i kh ∆ t − L e k +1 h a i k +1 h , z a i (cid:33) = 0 , where L e k +1 h = ∇ v k +1 h − N (cid:88) α =1 ν ( α ) k +1 h (cid:16) s ( α ) k +1 h ⊗ m ( α ) k +1 h (cid:17) and D e h = sym L e h . Step 2:

Move the mesh by u kh = v kh − mesh ∆ t . The velocity of the meshmotion v kh − mesh is taken as the material velocity v kh , which in fact leads to anupdated Lagrangian description.The mesh velocity in the second step can be chosen arbitrarily, so that v kh − mesh | ∂ Ω = v kh | ∂ Ω, which leads to the ALE description.

Acknowledgements

The authors acknowledge the support of GACR [grant number P107/12/0121].P. Minakowski also acknowledges the support of NSC (Poland) [grant number2012/07/N/ST1/03369]. The authors thank Martin Kruˇz´ık for general discus-sions on this topic, Jaroslav Hron for speciﬁc discussions on numerics-relatedissues and Endre S¨uli for several suggestions improving the ﬁnal form of thetext.