Plastic deformation treated as material flow through adjustable crystal lattice

Looking at severe plastic deformation experiments, it seems that crystalline materials at yield behave as a special kind of anisotropic, highly viscous fluids flowing through an adjustable crystal lattice space. High viscosity provides a possibility to describe the flow as a quasi-static process, where inertial and other body forces can be neglected. The flow through the lattice space is restricted to preferred crystallographic planes and directions causing anisotropy. In the deformation process the lattice is strained and rotated. The proposed model is based on the rate form of the decomposition rule: the velocity gradient consists of the lattice velocity gradient and the sum of the velocity gradients corresponding to the slip rates of individual slip systems. The proposed crystal plasticity model allowing for large deformations is treated as the flow-adjusted boundary value problem. As a test example we analyze a plastic flow of an single crystal compressed in a channel die. We propose three step algorithm of finite element discretization for a numerical solution in the Arbitrary Lagrangian Eulerian (ALE) configuration.


Introduction
Considering severe plastic deformation experiments as a motivation [1,2,3,4,5], plastic behaviour of crystalline solids is treated as a material flow through an adjustable crystal lattice space 1 . This point of view has been regarded by Asaro [6] as crystal plasticity "basic tenet". It seems that crystalline materials at yield behave as a special kind of incompressible, anisotropic, highly viscous fluids 2 . High viscosity provides a possibility to describe the flow as a quasi-static process, where inertial and other body forces can be neglected. The flow through the lattice space is restricted to preferred crystallographic planes and directions causing anisotropy. In the deformation process the lattice space is adjusted to the material flow by a lattice distortion, i.e. by rotation and stretching. The lattice space is treated as a solid, its distortion is measured with respect to a lattice reference configuration. 1 The term "space" is used deliberately, as the crystal lattice is understood as a space of preferred positions in a crystalline phase considered regardless of particles staying or flowing through them. 2 We mean a fluid in a generalized sense proposed e.g. by Rajagopal et al. [7]. Let us note that traditional solution of finite crystal plasticity is based on Lagrangian convected coordinate representation, see e.g. [8]. The Eulerian approach has been employed for a rigid-plastic model in a recent paper [9] with emphasis on dynamical problems. While the Eulerian formulation is very well suited for flow-like problems, as for example an equal channel angular extrusion experiment [10], in the case of a simple compression we deal with a free boundary problem. For this reason we employ the Arbitrary Lagrangian Eulerian (ALE) approach in the sense that we use the Eulerian formulation on a moving mesh which captures the free boundary.
The organization of the paper is as follows. In Section 2, we introduce the kinematics and formulate the flow model. The model of the plane strain compression and the finite element formulation are presented in Section 3. The simulation results are given and compared with the compression tests and computation of Harren et al. [11] in Section 4. Moreover, in Appendix we provide detailed derivation of the evolution equation for the Cauchy stress.

Flow model
Kinematics of the model is characterized by the velocity field v(x, t) and the lattice deformation gradient F e (x, t); x is a position at the current configuration, t means time, v describes the material flow in the current configuration, and F e controls the adjustment of the lattice space in the current configuration with respect to the lattice reference configuration.
The velocity gradient decomposes as follows L = ∇v = D + W , where D = (∇v + (∇v) T )/2 is the stretching and W = (∇v − (∇v) T )/2 is the spin; ∇ denotes the spatial gradient with respect to the coordinates x in the current configuration. The velocity gradient L consists of the lattice rate L e =Ḟ e (F e ) −1 and the plastic rate L p (a superposed dot means time derivative with respect to lattice space) where D e and D p are the symmetric parts of L e and L p , respectively. The lattice deformation gradient F e involves small elastic stretches and rigid body rotations of the lattice space as expressed in the polar decomposition where R e is the lattice rotation field. The right stretch tensor U e is related to the Green tensor E e = [(U e ) 2 − I]/2 = [(F e ) T F e − I]/2, where I stands for the identity matrix. The rate of E e can be expressed aṡ The material flow takes place on prescribed slip systems (i), i = 1, 2, . . . , N ; N is the number of slip systems. In the current configuration (i) slip system is defined by the unit vector s (i) in the direction of slip and by the unit normal to the glide plane m (i) . The lattice vectors s (i) and m (i) are transformed from the lattice reference configuration into the current configuration, where s  By taking the time derivative of the relations (4) and substitutingḞ e from (1), we obtain the equation that describes the evolution of the slip and normal directionṡ The deformation process is governed by the Cauchy stress T (x, t). The fields v and T have to satisfy the stress equilibrium where stands for the density subjected to the mass balancė The evolution of the Cauchy stress is given bẏ For detailed derivation of (9) and definition of the fourth order elastic tensor C, we refer to Appendix. The second law of thermodynamics expressed in a form convenient for the present consideration reads Ψ represents a Helmholtz potential. The Cauchy stress T controls the slip rates ν (i) through the resolved shear stresses τ (i) Considering rate dependent (visco-plastic) version of the flow model the resolved shear stresses τ (i) are assumed to be coupled with the slip rates ν (i) through a power law constitutive equation (viscous-plastic yield condition assumed by Kuroda [12]) where ν 0 is the characteristic slip rate, g (i) > 0 is the slip resistance, and m > 0 is a material parameter which controls the rate sensitivity. The slip resistances g (i) , i = 1, 2, . . . , N , are governed by the evolution equations, where the hardening matrix components H ij = H 0 sech 2 (H 0 ν acc /(g s − g 0 )), with the initial hardening rate H 0 = 8.9g 0 , the saturation strength g s = 1.8g 0 and the accumulated slip ν acc (t) = i t 0 |ν (i) |dt, see [6]. We conclude this section with formulation of an initial boundary value problem for the unknown density , the velocity v, the Cauchy stress tensor T , the slip directions s (i) , and the normal to slip planes m (i) , which is given by the system of equations where ν (i) is given by (12), (13). The system (14) holds for all x ∈ Ω and all t ∈ [0, ∞) and is endowed with the initial and boundary conditions where ∂Ω = Γ D ∪ Γ N and 0 , v 0 , φ i 0 are given functions.

Plane strain compression of single crystal
The model considered in this section is inspired by the detailed study of structure and micromechanisms of strain localization process during plane strain compression of f.c.c. single crystals in a channel die [11]. X-ray measurements were carried out to determine the lattice reorientation of the single crystals of various initial crystallographic orientations prior macroscopic shear bands appeared, i.e. up to the engineering compression strain 0.3-0.9. It was observed that except crystal symmetrically oriented with respect the compression and extension axes, the crystals had an overall common behaviour. After yield all these crystals start to reorient in such a way that their crystallographic direction

Channel-die compression
To test the present approach we consider the two-dimensional single crystal plane strain compression simulated in Harren et al. [11]. The model is represented by a rectangular-shaped single crystal with two or three active slip systems, that is subjected to channel-die compression as shown in Figures 1 and 2.
The initial boundary conditions are taken as in (15), with Γ D = Γ 1 ∪ Γ 3 and Γ N = Γ 2 ∪ Γ 4 . Namely, for the velocity (v = (v 1 , v 2 )) and the Cauchy stress are assumed to be Characteristic values are set to satisfy following ratios where E stands for Young's modulus and V , L, ν 0 are the reference velocity, length and slip rate, respectively. Note, that g 0 was defined in (13).

Finite element formulation
In the employed Arbitrary Lagrangian Eulerian (ALE) approach the time is discretized by a one step finite difference and a mixed finite element discretization is used in space. The choice of the approximations of each variable consists of P 2 elements for the velocity (v = (v x , v y )), P 1 for the density and the lattice rotations (all continuous) and P 1 -discontinuous for the Cauchy stress and the slip rates. This leads to the following definitions of the finite-dimensional spaces Here, T h is the triangulation of our domain with h representing the discretization parametr, for example the characteristic size of the triangles. The time interval (0, T ) is divided into steps of length ∆t and by superscript k we denote the variables at the time level t k .
For the solution algorithm of the system (14) at each timestep the set of the equations is decomposed into two parts which are solved subsequently in three sub-steps. Given v k−1 , k−1 ,  Step 1: Solve problem for v k , k , T k , and ν (i)k Step 2: Move the mesh by u k = v k mesh ∆t. The velocity of the mesh motion v k mesh can be taken as the material velocity v k which in fact leads to a Lagrangian description or it can be choosen arbitrary with restriction that v k mesh |∂Ω = v k |∂Ω which leads to the ALE description.
Step 3: Compute the new critical stress, the slip directions, and the new normal to the slip planes (17) and (18) are formulated in standard weak sense, discretized by means of the mixed FE spaces mentioned above. To solve the systems we employ a non-linear Newton-Raphson solver with an analytic evaluation of the Jacobian and use sparse direct solver MUMPS for solving the linear systems. Implementation is done in the software package FEniCS [13].
Simulations were conducted on a structured triangular mesh consists of 37688 vortices (74800 cells).

Non-symmetric double slip
The method presented in Section 3.2 was employed first to solve initial boundary value problem (14) with the non-symmetric double slip. We address the problem of the compression of the single crystal with the non-zero initial angle α 0 = 0 of the lattice rotation with respect to the compression axis. In Figure 3 we present result of simulations with the different initial orientations for the initially square and rectangular shaped domains. The initial orientations were set at 0.2 rad = 11.46 • and 0.4 rad = 22.9 • .
The simulations indicate that the angle of lattice rotation has a tendency to decrease; during die-compression the orientation of the slip systems approaches the symmetric orientation. Such reorientation effect was observed in most the single crystals measured by X-ray in [11]. Figure 4 shows the behavior of individual slip systems. Notice the difference between the alignment of the slip activity regions (shear bands) between the square and rectangle shaped domains. It is also worth of noticing that if we add third slip system perpendicular to the compression axis (α 3 0 = α 0 + π/2), the third slip system is almost inactive. For similar studies we refer to [9].
In the investigation of the influence of the initial angle of the lattice rotation we took the rate sensitivity parameter 1/m = 20, i.e. the value significantly smaller then 1/m = 200 employed in the symmetric case. The reason was that for high 1/m the suggested method is not yet able to capture the case presented in Section 4.2 for initial angle different then zero.
For similar computational studies we refer to [9] where the authors considered a rigid-plastic model with 3 slip systems and the initial orientation of 0.417 rad.

Symmetric double slip
As a special case we consider the symmetric double slip for high rate sensitivity parameter m = 0.005. The values of the magnitude of the velocity, the density, and the slip rates for three different strains are presented in Figure 5. By strain we understand the nominal strain e = (h − H 0 )/H 0 , where h and H 0 stand for current and initial height of the specimen, respectively. All the results where computed up to nominal strain 0.3.
The shear banding seen in Figure 5 to occur, some sort of a material or geometric imperfection must be present in the specimen. Following [11] the imperfection of the specimen was introduce by the perturbation of the right-hand side of the boundary, namely

Conclusions
• Plastic behavior of crystalline solids is treated as a material flow through an adjustable crystal lattice space. The material flow is characterized by the velocity field and the lattice space is treated as a solid, its distortion is measured with respect to a lattice reference configuration. • While the Eulerian formulation is very well suited for flow-like problems, as for example an equal channel angular extrusion experiment [10], in the case of a simple compression simulated here we deal with a free boundary problem. For this reason we employ the Arbitrary Lagrangian Eulerian (ALE) approach in the sense that we use the Eulerian formulation on a moving mesh which captures the free boundary.  • To test the proposed approach based on ALE method we consider the two-dimensional single crystal plane strain compression as in Harren et al. [11]. In agreement with Harren's et al observations, our model predicts that the lattice reorientation of the single crystals of various initial crystallographic orientations tends to bring the crystals toward the geometry coincided with a stable state of symmetric slip. The simulations in the case of symmetric double slip lead to the formation of the shear bands which corresponds to Harren's et al results.
• Since our method is based on Eulerian formulation it should be capable of solving deformations with very large strains for suitably defined problems, i.e. in the sense that we treat the material as a fluid. The reason is that the priority is not the displacement and the strain but the velocity of the material flow and the distortion of the crystal lattice space.
Here, when combined with mesh moving to capture the free boundary motion the method is able to reproduce a solutions obtained by Lagrangian based methods, as e.g. in [11].   where C is a fourth order material elasticity tensor (as a derivative of tensorial function w.r.t. second order tensor) and C is a fourth order spatial elasticity tensor. They are related by On the left hand side of the formula (24) we obtained Truesdell rate of the Kirchhoff stress. Equation (24) is equivalent to the one containing the Jaumann rate of the Kirchhoff stresṡ S − W e S + SW e = C(D − D p ) + D e S + SD e (26) It is common to neglect term D e S + SD e on the right hand side of (26), see eg. [14]. Without this assumption the Helmholtz potential approach does not lead to the Jaumann rate of the Kirchhoff stress, see [15].
To reformulate (24) in terms of the Cauchy stress we substitute S = T / and deducė Moreover we substitute the Helmholtz potential to the dissipation formula ( If the dissipation is caused by the material flow only, it follows from eqs. (10), (5) and (11) that The validity of the dissipation inequality (10) is guaranteed by the requirement g (i) ≥ 0 and the yield conditions (12).