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

TL;DR

This paper introduces a thermodynamically compatible Eulerian model for ideal plasticity in crystalline solids, treating deformation as a flow through an adjustable lattice, and applies it to 3D micropillar compression simulations.

Contribution

It extends Gibbs-potential-based formulations to model crystalline plasticity as a material flow with an adjustable lattice, incorporating a new computational framework.

Findings

Model successfully simulates 3D micropillar compression.

Framework aligns with thermodynamic principles.

Results compare favorably with existing studies.

Abstract

We propose an Eulerian thermodynamically compatible model for ideal plasticity of crystalline solids treated as a material flow through an adjustable crystal lattice space. The model is based on the additive splitting of the velocity gradient into the crystal lattice part and the plastic part. The approach extends a Gibbs-potential-based formulation developed by Rajagopal and Srinivasa for obtaining the response functions for elasto-visco-plastic crystals. The framework makes constitutive assumptions for two scalar functions: the Gibbs potential and the rate of dissipation. The constitutive equations relating the stress to kinematical quantities is then determined using the condition that the rate of dissipation is maximal providing that the relevant constraints are met. The proposed model is applied to three-dimensional micropillar compression, and its features, both on the level of modelling and computer simulations, are discussed and compared to relevant studies.

Tables

Table 1. Table 1: The slip directions and normals for a FCC crystal in Miller notation.

$\alpha$	$\sqrt{2}{𝒔}^{\left(\alpha \right)}$	$\sqrt{3}{𝒎}^{\left(\alpha \right)}$	$\alpha$	$\sqrt{2}{𝒔}^{\left(\alpha \right)}$	$\sqrt{3}{𝒎}^{\left(\alpha \right)}$	$\alpha$	$\sqrt{2}{𝒔}^{\left(\alpha \right)}$	$\sqrt{3}{𝒎}^{\left(\alpha \right)}$
$1$	$\left[\overline{1}10\right]$	$\left(111\right)$	$5$	$\left[101\right]$	$\left(1\overline{1}\overline{1}\right)$	$9$	$\left[0\overline{1}\overline{1}\right]$	$\left(\overline{1}1\overline{1}\right)$
$2$	$\left[10\overline{1}\right]$	$\left(111\right)$	$6$	$\left[01\overline{1}\right]$	$\left(1\overline{1}\overline{1}\right)$	$10$	$\left[1\overline{1}0\right]$	$\left(\overline{1}\overline{1}1\right)$
$3$	$\left[0\overline{1}1\right]$	$\left(111\right)$	$7$	$\left[110\right]$	$\left(\overline{1}1\overline{1}\right)$	$11$	$\left[\overline{1}0\overline{1}\right]$	$\left(\overline{1}\overline{1}1\right)$
$4$	$\left[\overline{1}\overline{1}0\right]$	$\left(1\overline{1}\overline{1}\right)$	$8$	$\left[\overline{1}01\right]$	$\left(\overline{1}1\overline{1}\right)$	$12$	$\left[011\right]$	$\left(\overline{1}\overline{1}1\right)$

Equations

$$\mathbf{L}={\mathbf{L}}_{\mathbf{e}}+{\mathbf{L}}_{\mathbf{p}}.$$

$$\mathbf{f}(\mathbf{T},\dot{\mathbf{T}},\mathbf{L})=0,$$

$$\mathbf{D}:=(\mathbf{L}+{\mathbf{L}}^T)/2\text{ and }\mathbf{W}:=(\mathbf{L}-{\mathbf{L}}^T)/2.$$

$$\dot{\varrho }+\varrho \text{div}\mathbf{v}$$

$$\varrho \dot{\mathbf{v}}$$

$$\varrho \dot{ϵ}=\mathbf{T}:\mathbf{D}-\text{div}\mathbf{q},$$

$$\varrho \zeta :=\varrho \dot{\eta }+\text{div}\left(\frac{\mathbf{q}}{\theta }\right)\text{ and }\zeta \ge 0,$$

$$\psi :=ϵ-\theta \eta ,\mathbf{S}:=\mathbf{T}/\varrho \text{ and }\xi :=\theta \zeta ,$$

$$\xi =\mathbf{S}:\mathbf{D}-\dot{\psi }-\eta \dot{\theta }-\mathbf{q}\cdot \frac{\nabla \theta }{\varrho \theta }\text{ and }\xi \ge 0.$$

$$G=G(\theta ,\mathbf{S})\text{ or }G(t,x)=G(\theta (t,x),\mathbf{S}(t,x)).$$

$$\begin{array}{rl}\psi (\theta ,\mathbf{S})& =G(\theta ,\mathbf{S})-\frac{\partial G(\theta ,\mathbf{S})}{\partial \mathbf{S}}:\mathbf{S},\\ ϵ(\theta ,\mathbf{S})& =G(\theta ,\mathbf{S})-\frac{\partial G(\theta ,\mathbf{S})}{\partial \mathbf{S}}:\mathbf{S}-\frac{\partial G(\theta ,\mathbf{S})}{\partial \theta }\theta ,\\ \eta (\theta ,\mathbf{S})& =-\frac{\partial G(\theta ,\mathbf{S})}{\partial \theta }\theta ,\end{array}$$

$$\xi =\mathbf{S}:\left\{\mathbf{D}+\frac{{\partial }^{2}G}{\partial {\mathbf{S}}^2}\dot{\mathbf{S}}+\frac{{\partial }^{2}G}{\partial {\theta }^2}\dot{\theta }\right\}-\mathbf{q}\cdot \frac{\nabla \theta }{\varrho \theta }\text{ and }\xi \ge 0.$$

$$\xi =\mathbf{S}:\mathbf{D}+\mathbf{S}:\frac{{\partial }^{2}G}{\partial {\mathbf{S}}^2}\dot{\mathbf{S}}\text{ and }\xi \ge 0.$$

$$G=G(\overline{\mathbf{S}})=G({\mathbf{R}}^T\mathbf{S}\mathbf{R}),$$

$$\xi =\mathbf{S}:\mathbf{D}+\overline{\mathbf{S}}:\frac{{\partial }^{2}G}{\partial {\overline{\mathbf{S}}}^2}\dot{\overline{\mathbf{S}}}=\mathbf{S}:\left\{\mathbf{D}-\mathbf{R}\overline{\mathcal{A}}\dot{\overline{\mathbf{S}}}{\mathbf{R}}^T\right\}\text{ and }\xi \ge 0,$$

$$G(\overline{\mathbf{S}})=-\frac{1}{2}\overline{\mathbf{S}}:\overline{\mathcal{A}}\overline{\mathbf{S}}\text{ and }\overline{\mathcal{A}}=-({\partial }^{2}G)/(\partial {\overline{\mathbf{S}}}^2).$$

$$\begin{array}{rl}\mathbf{S}:\left\{\mathbf{D}-\mathbf{R}\overline{\mathcal{A}}\dot{\overline{\mathbf{S}}}{\mathbf{R}}^T\right\}& =\mathbf{S}:\left\{\mathbf{D}-\mathbf{R}\overline{\mathcal{A}}{\mathbf{R}}^T\left(\dot{\mathbf{S}}+\mathbf{R}{\dot{\mathbf{R}}}^T\mathbf{S}+\mathbf{S}\dot{\mathbf{R}}{\mathbf{R}}^T\right)\mathbf{R}{\mathbf{R}}^T\right\}\\ & =\mathbf{S}:\left\{\mathbf{D}-\mathcal{A}(\dot{\mathbf{S}}+\mathbf{S}\dot{\mathbf{R}}{\mathbf{R}}^T+(\dot{\mathbf{R}}{\mathbf{R}}^T{)}^T\mathbf{S})\right\}\\ & =\mathbf{S}:\left\{\mathbf{D}-\mathcal{A}(\dot{\mathbf{S}}+\mathbf{S}\mathbf{\Omega }-\mathbf{\Omega }\mathbf{S})\right\},\end{array}$$

$$\stackrel{▽}{S}:=\dot{\mathbf{S}}+\mathbf{S}\mathbf{\Omega }-\mathbf{\Omega }\mathbf{S},$$

$$\xi =\mathbf{S}:(\mathbf{D}-\mathcal{A}\stackrel{▽}{S})\text{ and }\xi \ge 0,$$

$$\mathbf{D}=\mathcal{A}\stackrel{▽}{S}\text{ i.e. }\mathbf{D}=\mathcal{A}\{\dot{\mathbf{S}}+\mathbf{S}\mathbf{\Omega }-\mathbf{\Omega }\mathbf{S}\}.$$

$${\mathbf{a}}^1=\frac{{\mathbf{a}}_2\times {\mathbf{a}}_3}{{\mathbf{a}}_1\cdot ({\mathbf{a}}_2\times {\mathbf{a}}_3)},$$

$${\mathbf{s}}^{(\alpha )}=\sum _{i=1}^3{s}_i^{\alpha }{\mathbf{a}}_{\mathbf{i}},{\mathbf{m}}^{(\alpha )}=\sum _{i=1}^3{m}_i^{\alpha }{\mathbf{a}}^i,$$

$${\mathbf{a}}_{\mathbf{i}}(\mathbf{x},t)={\mathbf{F}}_{\mathbf{e}}(\mathbf{x},t){\mathbf{a}}_{\mathbf{i}}(\mathbf{x},0),$$

$${\dot{\mathbf{a}}}_i(\mathbf{x},t)={\mathbf{L}}_{\mathbf{e}}{\mathbf{a}}_{\mathbf{i}}(\mathbf{x},t).$$

$${\mathbf{L}}_{\mathbf{e}}=\sum _{i=1}^3{\dot{\mathbf{a}}}_{\mathbf{i}}\otimes {\mathbf{a}}^{\mathbf{i}}⟺{\dot{\mathbf{a}}}_{\mathbf{j}}={\mathbf{L}}_{\mathbf{e}}{\mathbf{a}}_{\mathbf{j}}.$$

$$\mathbf{L}={\mathbf{L}}_{\mathbf{e}}+{\mathbf{L}}_{\mathbf{p}},$$

$${\mathbf{D}}_{\mathbf{e}}$$

$${\mathbf{D}}_{\mathbf{p}}$$

$$\mathbf{D}={\mathbf{D}}_{\mathbf{e}}+{\mathbf{D}}_{\mathbf{p}}\text{ and }\mathbf{W}={\mathbf{W}}_{\mathbf{e}}+{\mathbf{W}}_{\mathbf{p}}.$$

$${\mathbf{D}}_{\mathbf{e}}=\mathcal{A}\stackrel{▽}{S}=\text{sym}\left(\sum _{i=1}^3{\dot{\mathbf{a}}}_{\mathbf{i}}\otimes {\mathbf{a}}^{\mathbf{i}}\right)\text{ and }{\mathbf{L}}_{\mathbf{e}}=\sum _{i=1}^3{\dot{\mathbf{a}}}_{\mathbf{i}}\otimes {\mathbf{a}}^{\mathbf{i}}.$$

$${\tau }^{(\alpha )}={\mathbf{s}}^{(\alpha )}\cdot \mathbf{S}{\mathbf{m}}^{(\alpha )}.$$

Full text

††journal: International Journal of Plasticity

1 Introduction

Severe plastic deformation experiments [30, 31] reveal that a crystalline material at yield can be seen as an anisotropic, highly viscous fluid. A structural adjustment of the crystal lattice space111The 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. to the material flow is seen as a deformation substructure. The flow space is restricted to preferred crystallographic planes and the directions causing anisotropy. Let us note that traditionally finite crystal plasticity is based on a Lagrangian convected coordinate representation. The presented form of crystal plasticity provides a possibility to treat the flow-adjustment boundary-value problem alternatively as a fluid flow in the Eulerian representation [14].

The traditional multiplicative split of the deformation gradient $\boldsymbol{F}$ into an elastic part $\boldsymbol{F_{e}}$ and a plastic part $\boldsymbol{F_{p}}$ in the form $\boldsymbol{F}=\boldsymbol{F_{e}}\boldsymbol{F_{p}}$ does not have a sound physical meaning. The sequential split means that the plastic action $\boldsymbol{F_{p}}$ occurs first, followed by the elastic action $\boldsymbol{F_{e}}$. However, that is not the way plasticity actually evolves; rather, elastic and plastic actions occur simultaneously. The formula $\boldsymbol{F}=\boldsymbol{F_{e}}\boldsymbol{F_{p}}$ is good for visualization purpose only, see [27].

The additive splitting of the velocity gradient $\boldsymbol{L}=\nabla\boldsymbol{v}$, where $\boldsymbol{v}$ is the velocity field, into the elastic part $\boldsymbol{L_{e}}$ and plastic part $\boldsymbol{L_{p}}$, i.e.,

[TABLE]

is far more natural and corresponds to the concept of plastic flow. The deformation gradient $\boldsymbol{F}$ can be found simply by solving the differential equation $\dot{\boldsymbol{F}}=\boldsymbol{L}\boldsymbol{F}$, where the superposed dot denotes material time-differentiation. In crystal plasticity, the parts $\boldsymbol{L_{e}}$ and $\boldsymbol{L_{p}}$ are specified in the following way: the elastic part $\boldsymbol{L_{e}}$ is identified with the evolution of the lattice space vectors, while $\boldsymbol{L_{p}}$ is determined constitutively; the Gibbs approach provides a suitable framework to achieve this, as will be shown in the next section.

For modelling of crystal plasticity as a fluid we adopt the Gibbs-potential-based framework proposed by Rajagopal and Shrinivasa [21]. They demonstrated the use of the Gibbs-potential-based formulation as a suitable tool for developing a thermodynamically consistent model for a wide class of elastic materials and viscoelastic fluids that are characterized by an implicit constitutive equation of the form

[TABLE]

where $\boldsymbol{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 assumptions that the Gibbs potential is a function of the stress and the rate of dissipation is non-negative. In fact, the latter is strengthened by requiring that the rate of dissipation is maximum possible. The approach thus leads to models that satisfy the second law of thermodynamics, automatically. This thermodynamical framework, that is fully Eulerian (using the current state as a reference state), yields a constitutive equation of the form (2) from specifying the constitutive equations for two scalar functions: the Gibbs potential and the rate of dissipation. The importance of choosing the Gibbs potential from the set of thermodynamical potentials (including further the Helmholtz free energy, the internal energy and the enthalpy) stems from the fact that it allows one to use the stress as a primitive quantity instead of a kinematical quantity, as a measure of strain. As pointed out by Rajagopal and Srinivasa [21], from a causal point of view, the traction (stress) is the cause, and the kinematics (motion or deformation) being the effect, which motivates to use of the stress as the primitive quantity. Requiring that the Gibbs potential is a function of a suitably chosen rotated stress tensor, Rajagopal and Srinivasa [21] extend their framework to model the response of anisotropic media. They obtain rate type fluid models for both anisotropic elastic materials and for anisotropic viscoelastic materials of Maxwell type.

The objective of this paper is to refine and extend the Gibbs framework in order to derive models that are capable of describing ideally plastic deformations of crystalline solids considered as a flow of material through an adjustable crystal lattice and that are furthermore compatible with the basic principles of continuum thermodynamics. The reason for focusing our attention to ideal plasticity, i.e. for ignoring work hardening or softening effects, is that an extension of the proposed model should not modify essential features of the Gibbs approach. To incorporate work hardening or softening would mean to enrich the model 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 proposed, analysed and tested by Cazacu and Ionescu [2, 4, 3]. Their goal was to develop an Eulerian rate-dependent model that is suitable for high strain rates, large strains and rotations of the incompressible material and the crystal lattice. Particular attention is devoted to the description of the kinematics of the crystal lattice in the Eulerian coordinates. The viscoplastic constitutive law and the equations for the evolution of the crystal lattice are expressed in terms of quantities attached to the current configuration. In applications involving large deformations and high strain rates the elastic component of strain is small with respect to the inelastic one, and it can be therefore neglected. Using this argument a rigid-viscoplastic approach has been adopted. Modelling has been focused on in-plane deformation and on the role of inertia in the material response.

The novelty of the approach described herein is in ensuring that the considered models are compatible with the second law of thermodynamics, thanks to the fact that their derivation is based on the specification of constitutive assumptions for the Gibbs potential and the rate of dissipation, and by an application of the maximal rate of entropy production principle. In contrast with the common restriction to incompressible rigid body motions, the proposed model covers the case of fully elastic response. We also refrain from making the assumption that the material is incompressible since in processes such as high pressure torsion incompressibility seems to be a restrictive approximation. Compressibility, on the other hand, may relax non-physically high stresses in a 2-turn equal channel angular extrusion, cf. [15].

We further illustrate the capabilities and the efficiency of the proposed model by performing numerical simulations for solving a three-dimensional micropillar compression problem. A material with face-centred cubic symmetry (FCC material) with 12 slip systems is considered. A numerical method is briefly outlined; more details and further numerical results confirming the efficiency of the numerical method will be given in a forthcoming paper.

The organization of the paper is as follows. In Subsection 2.1, we first recall the main ideas of the Gibbs-potential-based approach. Then, in Subsection 2.2, we recall basic concepts of crystal plasticity, specify the forms of both the Gibbs potential and the rate of dissipation, and show that, together with balance equations, these constitutive equations for two scalar quantitities suffice for the derivation of the proposed elasto-visco-plastic model of crystal plasticity. The derived model is then compared with the Eulerian approach of Cazacu and Ionescu [4, 3] in Section 2.3. Section 3 briefly introduces a numerical method and provides the results of numerical simulations for a micropillar compression ansatz. We conclude with a short overview of our main results in Section 4.

2 A flow model of crystal plasticity

2.1 The Gibbs-potential-based formulation by Rajagopal and Srinivasa [21]

In this section we summarise the result achieved in [21]. These results will be essential for our further consideration.

We consider a body that occupies a configuration $\Omega_{t}$ at the current instant $t$. The current position of any particle is denoted by $\boldsymbol{x}$, its velocity by $\boldsymbol{v}$. Moreover, we introduce the symmetric and antisymmetric parts of the velocity gradient $\boldsymbol{L}=\nabla\boldsymbol{v}$ through

[TABLE]

The mass density of the material is denoted by $\varrho$ and the Cauchy stress by $\boldsymbol{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:

[TABLE]

where the superscript notation $\boldsymbol{A}^{T}$ means the transposed tensor to any tensor $\boldsymbol{A}$ and the material time derivative of a scalar function $z$ is given by $\dot{z}=z_{,t}+\nabla z\cdot\boldsymbol{v}$ (for a vector and tensor-valued function, the same relation is applied to each component). Furthermore, we will state the balance of energy (in the absence of the body heat supply) in the form

[TABLE]

where $\epsilon$ is the specific internal energy, $\boldsymbol{q}$ the heat flux vector and $\boldsymbol{D}$ the symmetric part of the velocity gradient introduced in (3). Finally, we express the second law of thermodynamics through

[TABLE]

where $\eta$ is the specific entropy, $\theta$ the temperature and $\zeta$ the specific rate of entropy production; here we tacitly assume that the entropy flux is of the form $\boldsymbol{q}/\theta$. Introducing the (specific) Helmholtz free energy $\psi$, the Kirchhoff stress tensor $\boldsymbol{S}$ and the (specific) rate of dissipation $\xi$ through

[TABLE]

and using (5) and (6), we arrive at the equation for the rate of dissipation

[TABLE]

The starting point of the Gibbs-potential-based framework as developed in [21] is the assumption that the (specific222We suppress the use of the word specific for relevant quantities in what follows although the quantities are taken per unit mass.) Gibbs potential, denoted by $G$, is a function of the temperature $\theta$ and the Kirchhoff stress $\boldsymbol{S}$, i.e.,

[TABLE]

Stipulating further the Helmoltz free energy, the internal energy and the entropy, as functions of $\theta$ and $\boldsymbol{S}$, through

[TABLE]

and inserting the first and third of these relations into (7), we obtain

[TABLE]

In what follows, we restrict ourselves to isothermal processes. Then, the equation (10) reduces to

[TABLE]

We have thus arrived at a representation of thermodynamics associated with the specification of the Gibbs potential (as given in (8)). There is however a problem: while $\boldsymbol{D}$ and $\boldsymbol{S}$ are both objective tensors, $\dot{\boldsymbol{S}}$ and consequently $({\partial^{2}G})/{(\partial\boldsymbol{S}^{2})}\,\dot{\boldsymbol{S}}$ are not objective tensors.

To overcome this difficulty and, moreover, to open the possibility to include anisotropic responses, Rajagopal and Srinivasa (see [21]) propose to consider, instead of (8), the Gibbs potential depending on a rotated stress $\overline{\boldsymbol{S}}=\boldsymbol{R}^{T}\boldsymbol{S}\boldsymbol{R}$, i.e.,

[TABLE]

where $\boldsymbol{R}$ is any rotation tensor that is objective (which means that $\boldsymbol{R}$, related to a motion $\boldsymbol{\chi}$, and $\boldsymbol{R}^{*}$ related to the motion $\boldsymbol{\chi}^{*}$, satisfy $\boldsymbol{R}^{*}=\boldsymbol{Q}\boldsymbol{R}$ whenever $\boldsymbol{\chi}$ and $\boldsymbol{\chi}^{*}$ differ by a rigid body rotation $\boldsymbol{Q}$, all quantities being functions of position and time).

Thus, assuming (12) (instead of (8)) and following the same scheme to the one that yields (11) from (8), we get

[TABLE]

where

[TABLE]

Substituting $\overline{\boldsymbol{S}}=\boldsymbol{R}^{T}\boldsymbol{S}\boldsymbol{R}$ into the expression on the right-hand side of (13) and using the fact that $\boldsymbol{R}^{T}\boldsymbol{R}=\boldsymbol{I}$ we obtain

[TABLE]

where $\boldsymbol{\mathcal{A}}$ is a fourth-order tensor such that $\mathcal{A}_{ijkl}=\boldsymbol{R}_{iA}\boldsymbol{R}_{jB}\boldsymbol{R}_{kC}\boldsymbol{R}_{lD}\overline{\mathcal{A}}_{ABCD}$ and $\boldsymbol{\Omega}=\dot{\boldsymbol{R}}\boldsymbol{R}^{T}$. Note that $\boldsymbol{R}\boldsymbol{R}^{T}=\boldsymbol{I}$ implies that $\dot{\boldsymbol{R}}\boldsymbol{R}^{T}=-\boldsymbol{R}\dot{\boldsymbol{R}}^{T}=-(\dot{\boldsymbol{R}}\boldsymbol{R}^{T})^{T}$ and $\boldsymbol{\Omega}$ is antisymmetric. Next, introducing the notation

[TABLE]

it is not difficult to observe that $\stackrel{{\scriptstyle\triangledown}}{{\boldsymbol{S}}}$ is an objective time derivative. With this notation, (13) takes the form

[TABLE]

where $\boldsymbol{S}$, $\boldsymbol{D}$ and $\boldsymbol{\mathcal{A}}\stackrel{{\scriptstyle\triangledown}}{{\boldsymbol{S}}}$ are objective and symmetric second-order tensors. In addition, this step provides a simple and natural way for introducing, through $\boldsymbol{\mathcal{A}}$, an anisotropic response as there is no restriction imposed by frame indifference on the dependence of $G$ on $\overline{\boldsymbol{S}}$.

Rajagopal and Srinivasa [21] use (16) to make, among others, the following observation: the dissipation rate $\xi$ vanishes for arbitrary $\boldsymbol{S}$ if

[TABLE]

This constitutive equation, which is of the form (2), characterizes elastic (nondissipative) materials 333 As is well known such elastic response cannot be achieved via the Helmholtz free energy approach, see [21] for further details and references.. In the next Section, the Gibbs potential (12) and the constitutive assumption (17) modified for elastic properties of the crystal lattice space will be used for modeling of crystal plasticity.

2.2 Crystal plasticity

A formulation of crystal plasticity in the framework of the Gibbs-potential approach is derived in two steps. First, we recall the standard description of a crystal lattice space as well as slip systems and identify the lattice velocity gradient ${\boldsymbol{L}}_{e}$ with the evolution of the lattice space. As a second step we specify the constitutive equation for the rate of dissipation and by applying the assumption that the response of the material is such that it maximizes this dissipation function (and other relevant constraints) we obtain the constitutive equation for $\boldsymbol{S}$.

A crystal lattice space can be characterized by three non-planar vectors $\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$ and $\boldsymbol{a}_{3}$, which form a basis in $\mathbb{R}^{3}$ and span the so-called Bravais lattice (as an example, a face-centred cubic lattice is sketched in Fig. 1a). The reciprocal dual lattice vectors are defined by

[TABLE]

with $\boldsymbol{a}^{2}$ and $\boldsymbol{a}^{3}$ defined by cyclic permutation; note that $\boldsymbol{a}^{i}\cdot\boldsymbol{a}_{j}=\delta_{ij}$.

In crystalline solids the material flow in a plastic regime is carried by a motion of dislocations. Crystal plasticity theory is based on the observation that the motion of dislocations takes place in preferred crystallographic directions (so-called slip directions) on preferred crystallographic planes (so-called slip planes); in our consideration climb of dislocations is excluded. The unit vector in the slip direction $\boldsymbol{s}^{(\alpha)}$ and the unit normal to the slip plane $\boldsymbol{m}^{(\alpha)}$ form the $\alpha$-slip system. The set of $N$ potentially active slip systems is one of the basic microscopic characteristics of crystal plasticity, hence $\alpha=1,\dots,N$. The vectors $\boldsymbol{s}^{(\alpha)}$ and $\boldsymbol{m}^{(\alpha)}$ of the slip system $\alpha$ can be expressed through the lattice vectors $\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$ and $\boldsymbol{a}_{3}$ as

[TABLE]

where $s^{\alpha}_{i}$ and $m^{\alpha}_{i}$ are constant coefficients, called Miller indexes. For example, for $s^{1}_{1}=1$, $s^{1}_{2}=1$ and $s^{1}_{3}=0$, and $m^{1}_{1}=1$, $m^{1}_{2}=1$ and $m^{1}_{3}=1$ the slip direction vector $\boldsymbol{s}^{1}=\boldsymbol{a}_{1}+\boldsymbol{a}_{2}$ and the slip plane normal $\boldsymbol{m}^{(1)}=\boldsymbol{a}_{1}+\boldsymbol{a}_{2}+\boldsymbol{a}_{3}$, expressed in a standard way through the Miller indexes, read [110] and (111), respectively; 110 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 and the material flow. In the proposed model, the crystal lattice is treated as a solid described by the lattice deformation gradient $\boldsymbol{F_{e}}$ from the natural (lattice reference) configuration 444The concept of the natural configuration has been developed and applied in a number of areas beyond the area of plasticity by Rajagopal and his co-workers, see for example [18, 19, 20]. In our setting, the natural configuration is the lattice reference configuration and thus

$\boldsymbol{a_{i}}(\boldsymbol{x},t)=\boldsymbol{F_{e}}(\boldsymbol{x},t)\,\boldsymbol{a_{i}}(\boldsymbol{x},0)\,,$

(20)

where $\boldsymbol{a_{i}}(\boldsymbol{x},0)$ are the lattice vectors in the natural (lattice reference) configuration. Since the lattice velocity gradient $\boldsymbol{L_{e}}$ is supposed to be linked with the rate of $\boldsymbol{F_{e}}$ through the relation $\boldsymbol{L_{e}}=\dot{\boldsymbol{F_{e}}}\boldsymbol{F_{e}}^{-1}$ (that is analogous to the kinematical relationship between $\boldsymbol{L}$ and $\dot{\boldsymbol{F}}$), the rate form of (20) reads

${\dot{\boldsymbol{a}}_{i}}(\boldsymbol{x},t)=\boldsymbol{L_{e}}\,\boldsymbol{a_{i}}(\boldsymbol{x},t)\,.$

(21) to the current configuration. According to Srinivasa and Srinivasan [27, Section 13.4, p. 467] this evolution equation for the lattice vectors $\boldsymbol{a_{i}}$ is characterised in the following way

[TABLE]

Indeed, to see that (22) holds, we first multiply $\sum_{i=1}^{3}\dot{\boldsymbol{a}}_{\boldsymbol{i}}\otimes\boldsymbol{a^{i}}$ by $\boldsymbol{a_{j}}$ ($\boldsymbol{a^{i}}\cdot\boldsymbol{a_{j}}=\delta_{ij}$) and observe that $\boldsymbol{L_{e}}=\sum_{i=1}^{3}\dot{\boldsymbol{a}}_{\boldsymbol{i}}\otimes\boldsymbol{a^{i}}\Rightarrow\dot{\boldsymbol{a}}_{\boldsymbol{j}}=\boldsymbol{L_{e}}\boldsymbol{a_{j}}.$ On the other hand, the opposite implication, i.e. $\boldsymbol{L_{e}}=\sum_{i=1}^{3}\dot{\boldsymbol{a}}_{\boldsymbol{i}}\otimes\boldsymbol{a^{i}}\Leftarrow\dot{\boldsymbol{a_{j}}}=\boldsymbol{L_{e}}\boldsymbol{a_{j}}$, follows from the definition of the reciprocal basis (18). Note that the lattice vectors remain orthogonal, since they are transformed by the same lattice deformation gradient $\boldsymbol{F_{e}}$555 We decompose the lattice deformation gradient $\boldsymbol{F_{e}}$ into the elastic stretch $\boldsymbol{U_{e}}$ and the rotation $\boldsymbol{R_{e}}$, $\boldsymbol{F_{e}}=\boldsymbol{R_{e}}\boldsymbol{U_{e}}$. The angle is preserved by the rotation $\boldsymbol{R_{e}}$, and the loss of orthonormality is caused by the elastic stretching..

Our basic kinematical assumption is represented by the additive splitting of the velocity gradient in the form (1), i.e.

[TABLE]

where $\boldsymbol{L_{e}}$ is linked with the $\boldsymbol{a}_{\boldsymbol{i}}$ via the equations (22) and the form of $\boldsymbol{L_{p}}$ will be specified by maximization of the rate of dissipation based on the identification of the constitutive equation for the rate of dissipation.

Following the decomposition (23), we set

[TABLE]

and from (3) we have that

[TABLE]

Considering (12) and (16) with the constitutive equation (17) modified to the crystal lattice space $\boldsymbol{D_{e}}=\boldsymbol{\mathcal{A}}\,\stackrel{{\scriptstyle\triangledown}}{{\boldsymbol{S}}}$ we get from (22) that

[TABLE]

Next, we employ the procedure of maximization of the rate of dissipation $\xi$ to get $\boldsymbol{D}_{p}$; this, together with the second equation in (22), also allows us to determine uniquely the rate of material flow $\boldsymbol{L_{p}}$ for a crystalline solid deformed by slip. To do so, we need to specify the constitutive equation for the rate of dissipation $\xi$. The driving force of slip is the resolved shear stress $\tau^{(\alpha)}$ controlled by the critical value $\tau^{(\alpha)}_{c}$ ([27]),

[TABLE]

For the sake of the following procedure the critical values $\tau^{(\alpha)}_{c}$ are assumed to be a given positive material parameter in agreement with our restriction to ideal plasticity. If $\tau^{(\alpha)}<\tau^{(\alpha)}_{c}$ the slip rate and the corresponding dissipation rate of the $\alpha$-slip system is very small. On the other hand, if $\tau^{(\alpha)}>\tau^{(\alpha)}_{c}$ the dissipation rate is high. One of the suitable possibilities is to assume that the dissipation rate is a sum of contributions of the slip systems proportional to a power of the quotients $\tau^{(\alpha)}/\tau^{(\alpha)}_{c}$, i.e.,

[TABLE]

where $c_{0}>0$ and $\beta>0$ are material constants.

Introducing the notation

[TABLE]

we observe that

[TABLE]

and consequently

[TABLE]

In order to determine the constitutive equation for the evolution of $\boldsymbol{S}$, we maximize $\hat{\xi}$ given in (29) with respect to $\boldsymbol{S}$ under the constraint (16) that can be rewritten in a more compact form, namely

[TABLE]

Hence, setting

[TABLE]

the constrained optimality condition reads

[TABLE]

Multiplying the resulting equation scalarly by $\boldsymbol{S}$, and using (33) and (32), we conclude that

[TABLE]

Consequently, the second condition in (34) together with (27) and (31) imply that

[TABLE]

Note that $\nu^{(\alpha)}$ can be interpreted as a slip rate of the $\alpha$-slip system. Thus, we have identified $\boldsymbol{D_{p}}$. Since (22) and (24) imply that $\boldsymbol{W}=\boldsymbol{W_{e}}+\boldsymbol{W_{p}}$ uniquely determines $\boldsymbol{W_{p}}$, we stipulate that

[TABLE]

Note that this formula is consistent with the results obtained so far for the symmetric parts $\boldsymbol{D}$, $\boldsymbol{D}_{e}$ and $\boldsymbol{D}_{p}$ since we have (see (35))

[TABLE]

In addition, the choice (36) and the decomposition rule (23) lead to

[TABLE]

which we view (and use) as the evolutionary equation for the vectors $\boldsymbol{a}_{i}$.

Recalling definition (15) and combining it with the obtained constitutive relation (35) we obtain

[TABLE]

where $\boldsymbol{\Omega}=\dot{\boldsymbol{R}}\boldsymbol{R}^{T}$ for any rotation tensor $\boldsymbol{R}$ that is objective. By choosing different rotations $\boldsymbol{R}$, one ends up with different objective time derivatives, and, consequently, with different models, with the restriction for $\boldsymbol{\Omega}$ being antisymmetric. For example, by choosing $\boldsymbol{R}$ as the rotational part $\boldsymbol{R_{e}}$ of the polar decomposition of the lattice deformation gradient $\boldsymbol{F_{e}}=\boldsymbol{R_{e}}\boldsymbol{U_{e}}$ will yield the Green-McInnis-Naghdi stress rate. The second example, namely the Zaremba-Jaumann derivative, extensively used in the crystal plasticity literature, results from setting $\boldsymbol{\Omega}=\dot{\boldsymbol{R}}\boldsymbol{R}^{T}=\boldsymbol{W_{e}}$.

From now on we focus on the Zaremba-Jaumann rate

[TABLE]

At this point we compare the constitutive relation (35) with relations that are typically used in crystal plasticity, see e.g. [16, 17]. In these studies, based on the Helmholtz potential approach, the resulting constitutive equation takes the form

[TABLE]

where the rate derivative on the left-hand side is the Oldroyd derivative. To obtain the Zaremba-Jaumann rate, the terms $\boldsymbol{D_{e}}\boldsymbol{S}+\boldsymbol{S}\boldsymbol{D_{e}}$ are often neglected or transferred to the right-hand side of (40) by defining the fourth-order tensor $\mathcal{\boldsymbol{K}}(\boldsymbol{D_{e}})=\boldsymbol{\mathcal{C}}(\boldsymbol{D_{e}})+\boldsymbol{D_{e}}\boldsymbol{S}+\boldsymbol{S}\boldsymbol{D_{e}}$. In this sense the Gibbs and the Helmholtz potential based approaches are comparable.

2.3 The model and its comparison to the Eulerian approach by Cazacu and Ionescu

In summary, the primary variables of the proposed crystal plasticity model based on the Gibbs potential are: the density $\varrho$, the velocity $\boldsymbol{v}$, the Kirchhoff stress $\boldsymbol{S}$, and the vectors $\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$ and $\boldsymbol{a}_{3}$ being functions of the current position $\boldsymbol{x}$ and time $t$. These variables are governed by the system of equations (4a), (4b), (35) and (38), i.e.,

[TABLE]

The slip system vectors $\boldsymbol{s}^{(\alpha)}$ and $\boldsymbol{m}^{(\alpha)}$ in (41d) are expressed in terms of the lattice vectors $\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$ and $\boldsymbol{a}_{3}$ through (18) and (19); the slip rates $\nu^{(\alpha)}$ are governed by the equations (30), (28).

We also multiply (41c) by $\boldsymbol{\mathcal{C}}:=\boldsymbol{\mathcal{A}}^{-1}$ and obtain

[TABLE]

The elastic tensor $\boldsymbol{\mathcal{C}}$ is anisotropic; in what follows we assume that $\boldsymbol{\mathcal{C}}$ is of cubic symmetry characterized by three independent elements, $c_{11}$, $c_{12}$, and $c_{44}$, which correspond to $\mathcal{C}_{1111}$, $\mathcal{C}_{1122}$, and $\mathcal{C}_{2332}$ respectively, through the standard Voigt notation.

In comparison with (41) the primary variables of the rigid-plastic model proposed by Cazacu and Ionescu [2, 3, 4] are: the velocity $\boldsymbol{v}$, which is divergence-free, the Cauchy stress $\boldsymbol{T}=\varrho\boldsymbol{S}$, and the slip system vectors $\boldsymbol{s}^{(\alpha)}$ and $\boldsymbol{m}^{(\alpha)}$, $\alpha=1,\ldots,N$ as functions of the current position $\boldsymbol{x}$ and time $t$. These variables are governed by the system of the equations:

[TABLE]

The expression in the brackets on the right-hand side of (43d) and (43e) is the rate of the lattice rotation $\dot{\boldsymbol{R}}_{e}\boldsymbol{R}_{e}^{-1}$ (for a rigid lattice $\boldsymbol{F_{e}}=\boldsymbol{R}_{e}$); (43d) and (43e) are the rigid plasticity versions of (41d) of our model evaluated directly in terms of the vectors of the slip system $\boldsymbol{s}^{(\alpha)}$ and $\boldsymbol{m}^{(\alpha)}$, $\alpha=1,\ldots,N$ instead of the lattice vectors $\boldsymbol{a_{i}}$, $i=1,2,3$ in our case. Since $N$ can be large (for e.g. for FCC crystal structure $N=12$), working with the lattice vectors $\boldsymbol{a_{i}}$ makes the system (41) advantageous with respect to the system (43). It reduces the number of unknowns in the system and thus plays an important role in the efficiency of three-dimensional numerical simulations.

Due to the assumed rigidity, the density $\varrho$ in (43b) is a material constant. The equation (43c) is an implicit constitutive equation for the stress $\boldsymbol{T}$ as could be seen if one used for the slip rates $\nu^{(\alpha)}$ the power law (30) and for the resolved shear stress $\tau^{(\alpha)}=\boldsymbol{s}^{(\alpha)}\cdot\boldsymbol{T}\,\boldsymbol{m}^{(\alpha)}$. Note that in this case all $N$ slip systems are active, and the power law determines their relative activity. However, Cazacu and Ionescu preferred, mainly from the computational point of view, to employ instead the power law (30) the visco-plastic extension of the classical rigid-plastic Schmid law using the overstress approach

[TABLE]

where $\eta^{\alpha}$ is the viscosity; $[z]_{+}=(z+|z|)/2$ denotes the positive part of the real number $z$. Using this overstress approach, it is important to mention that only up to 5 resolved shear stresses $\tau^{(\alpha)}$ can be independent (there are 5 components of the stress deviator controlling slip). The slip rates $\nu^{(\alpha)}$ given by (44) are also not independent as they have to satisfy the kinematic constrain (43c). Given the deformation rate $\boldsymbol{D}$, the slip rate $\nu^{(\alpha)}$ of the active slip systems can be determined by minimizing the internal power under the constrain (43c). The deviatoric part of the stress $\boldsymbol{T}$ corresponding to a given deformation rate $\boldsymbol{D}$ is obtained as the Lagrange multiplier of this minimization problem.

One of the main differences is that in the rigid-plastic model (43) the elastic properties of the lattice are neglected as they are assumed to be small in comparison with the large plastic deformations. This assumption is reasonable in modelling e.g. metal production processes and generally it is suitable from a macroscopic point of view. If the interest is focused on a meso-scale, where the deformation is spontaneously heterogeneous, not only the local misorientations of the crystal lattice but also lattice strains connected with a stress redistribution may play a role. It is worth mentioning that classical crystal plasticity models of shear bands consider anisotropic elasticity as an important ingredient. Moreover, elasticity and compressibility have a stabilizing effect, helping the computational strategy as indicated in Section 3, see also [15].

In some aspects the models (41) and (43) are similar. Both use the rate decomposition rule (23) in the current configuration as the basic kinematic assumption (the reference to decomposition rule (23) in Cazacu and Ionescu [4] is only indicative) and their Eulerian descriptions are conceptually equivalent (except the acceptance of elasticity and the difference in the power vs. overstress rules mentioned above). A substantial difference is in the formulations of the models. The rigid-plastic model (43) is based on the slip rate composition of the plastic velocity gradient (36) and the resolved shear stress (28) as the primary assumptions. On the other hand, the proposed constitutive modelling based on the Gibbs-potential formulation is constructed from two scalar functions of the Kirchhoff stress $\boldsymbol{S}$: the Gibbs potential $G(\overline{\boldsymbol{S}})$ given by (14), $\overline{\boldsymbol{S}}=\boldsymbol{R}_{e}^{T}\boldsymbol{S}\boldsymbol{R}_{e}$, and the rate of dissipation $\xi(\boldsymbol{S})$ given by (29). The form of the material flow gradient $\boldsymbol{L_{p}}$ given by (36) is the consequence of the dissipation rate maximization. Moreover, the Gibbs-potential framework guarantees that the model is thermodynamically consistent.

3 Numerical example: three-dimensional micropillar compression

The algorithm presented in this section is based on the Eulerian framework described above and thus contributes to the Crystal Plasticity Finite Element Method (CPFEM). A recent overview by Roters et al. [23, 22] summarizes applications of the CPFEM. Most of the methods that have been developed so far adopt a Lagrangian description of the continuum problem. In the context of the present paper we specifically highlight the following references focused on the application of CPFEM to micropillar compression: [8, 9, 12, 10, 26].

In [10], Jung et al. compare different primary slip plane inclination angles and use a hardening rule that accounts for anisotropic slip system interactions. Currently, it seems that the main interest in the field concentrates on finite element analysis of geometrically necessary dislocations, see e.g. [8, 9, 12]. Since the primary purpose of this study is to computationally analyze the Eulerian approach based on the Gibbs-potential-based formulation, we focus on the case of ideal plasticity here and postpone the study of the hardening (and possibly non-local) effects to our further research. Note, however that in Minakowski et al. [16] the authors presented a crystal plasticity model including isotropic hardening and proposed an Arbitrary Lagrangian Eulerian (ALE) based numerical method for a two-dimensional plastic flow of a single crystal compressed in a channel die.

To examine the potential of the derived model we consider three-dimensional micropillar compression. For comparison we use the same setup as Kuroda in [12], where the behaviour of single crystal micropillars subjected to compressive loading is described. A three-dimensional finite element method incorporating higher-order gradient crystal plasticity is analyzed. The author first simulated a compression test for a model without higher strain gradients, which we use for comparison.

The compression testing methodology is shown schematically in Figure 2 (like Kuroda [12] we consider samples with a conical base). The setting mimics the compression experiment commonly performed on macroscopic samples [28, 29, 25]. The compressive stress is defined as the sum of the nodal forces in the $-x_{3}$ direction at the top surface divided by the initial cross-section area. The nominal compressive strain is the ratio of the displacement of the top surface in the $-x_{3}$ direction over the original height of the sample.

We use the same set of values of parameters of the compressed material as in [12], which corresponds to the experimental and numerical works [5, 25]. The ratio between the height $L_{0}$ and the diameter $D$ is selected to be $2.3$. In a crystal with cubic symmetry, such as face-centred cubic FCC, with the Cartesian axes oriented along the cube edges (Figure 1), the non-zero elements of $\boldsymbol{\mathcal{C}}$ are the same ones as for the simplest anisotropic solid, the three values $c_{11}$, $c_{12}$ and $c_{44}$ are independent. The elastic constants are chosen to be $c_{11}=280.5\text{GPa}$, $c_{12}=188.5\text{MPa}$, and $c_{44}=133.5\text{GPa}$, see (42). The reference values are as follows: the length $l_{0}=10\mu\,m$, the velocity $\boldsymbol{v}_{0}=1\mu m/s$, the reference slip rate $c_{0}=0.1s$, the critical resolved shear stress $\tau_{c}=420\text{MPa}$. The rate sensitivity parameter $m$ is taken to be $0.05$, for $\beta=\frac{1}{m}$, compare (30). A review of computational strategies in an effort to overcome numerical instabilities 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 effects are neglected. These material parameters are those for a single-crystal Ni-base superalloy, which has the nominal composition of NiָCoַCrַTaֶ.2Alֵ WֳReֲMoְ.2Hf, see [25].

We solve the set of equations derived in Section 2, supplemented with proper boundary and initial conditions, see also (41). Let $T>0$, $\Omega\subset\mathbb{R}^{3}$ open and bounded, $\varrho_{0}:\Omega\to\Re$, $\boldsymbol{v}_{0}:\Omega\to\mathbb{R}^{3}$, $\boldsymbol{S}_{0}:\Omega\to\Re^{3\times 3}$, $\boldsymbol{a_{i}}_{0}:\Omega\to\mathbb{R}^{3},\quad i=1,2,3$, $\boldsymbol{v}_{D}:\Gamma_{D}\to\mathbb{R}^{3}$, $\boldsymbol{t}:\Gamma_{N}\to\mathbb{R}^{3}$. The system is satisfied inside the domain $\Omega\subset\mathbb{R}^{3}$, the boundary $\partial\Omega$ consists of two non-intersecting parts $\Gamma_{D}$ and $\Gamma_{N}$ corresponding to the Dirichlet and the Neumann boundary conditions, respectively. We look for the density $\varrho$, the velocity $\boldsymbol{v}$, the Kirchhoff stress $\boldsymbol{S}$, and the lattice vectors $\boldsymbol{a_{1}}$, $\boldsymbol{a_{2}}$, $\boldsymbol{a_{3}}$ ($\boldsymbol{s}^{(\alpha)}$ and $\boldsymbol{m}^{(\alpha)}$ are given through lattice vectors, equations (18) and (19)) satisfying

[TABLE]

where $\varrho_{0}(x)=1$, $\boldsymbol{S}_{0}(x)\boldsymbol{n}=\boldsymbol{0}\text{ for }x\in\Gamma_{N}\times[0,\infty)$, $\boldsymbol{S}_{0}(x)=\boldsymbol{0}\text{ for }x\in\Omega$, $\boldsymbol{a_{1}}_{0}(x)=(1,0,0)$, $\boldsymbol{a_{2}}_{0}(x)=(0,1,0)$, $\boldsymbol{a_{3}}_{0}(x)=(0,0,1)$ and $\boldsymbol{n}$ is the outer normal vector. Moreover, as already explained the slip rates $\nu^{(\alpha)}$ are governed by the equations (30), (28).

Two kinds of top velocity boundary condition are considered, see Figure 2 (b-c). In the first case, corresponding to Kuroda’s condition BC1 [12], we fix the velocity on the top boundary $\Gamma_{1}$ and require

[TABLE]

This condition is referred below to as the constrained case and the setting mimics the situation with a very stiff loading axis and the friction between the top surface and the plunger is extremely high. The second case (corresponding to Kuroda’s condition BC2 [12]) is called below the unconstrained case and it allows the top surface to move in the plane perpendicular to the $z$-axis, which means that we require that

[TABLE]

where $\boldsymbol{\tau_{1}}$, $\boldsymbol{\tau_{2}}$ are two vectors generating the basis of the tangent plane orthogonal to $\boldsymbol{n}$ at the boundary. This simply mimics the case where there is no friction between the top surface and the plunger. In the rest of the Dirichlet boundary, namely the bottom of the specimen, the velocity is fixed $\boldsymbol{v}_{0}(x,t)|_{\Gamma_{3}}=(0,0,0)$. The degrees of freedom for the velocity on the conical base are fully constrained. Moreover, we put the Neumann boundary condition $\boldsymbol{S}(x,t)\boldsymbol{n}=\boldsymbol{0}$ on the lateral surface $\Gamma_{2}$.

The crystal orientation is relative to the fixed sample coordinates in the reference configuration. The $[\bar{1}23]$ direction is set to be parallel to the sample axis, i.e. the $X_{3}$ direction, the $[1\bar{1}1]$ direction is chosen to coincide with the $X_{1}$ direction, and the $X_{2}$ direction coincides with $[54\bar{1}]$ (the overbar denotes the negative value). The specification of the slip systems is given in Table 1. The primary slip system is number 2, $(111)[10\bar{1}]$, the angle between $\boldsymbol{s}^{(2)}$ and $X_{1}-X_{2}$ plane reads $49.1^{\circ}$.

3.1 Finite element formulation

In the solution of the micropillar boundary-value problem without higher-order gradients, Kuroda [12] used a generally adopted finite element Lagrangian framework. In his numerical experiments, he evaluated deformed finite element meshes, contour maps of slip, maps of lattice rotations, and stress-strain diagrams for nominal compressive strain up to $0.2$. In the present example, we employ the Arbitrary Lagrangian Eulerian (ALE) approach in the sense that we use an Eulerian formulation on a moving mesh, which captures the free boundary of the micropillar. In this way, we are able to extend the tested range of nominal compressive strain up to $0.5$. For a brief description of our numerical scheme we refer to the Appendix. The detailed description concerning accuracy and efficiency of the numerical method and further numerical results of the application of the ALE method to three-dimensional crystal plasticity problems will be described in detail in a forthcoming work.

3.2 Results

Since the computations are conducted in the Eulerian coordinates our primary variable is the velocity. In Figure 3 we present the velocity field at various strain levels in the whole sample for the $X_{1}-X_{3}$ and $X_{2}-X_{3}$ plane views. Due to anisotropy of the acting slip systems the deformation profiles are not the same for different plane views, Fig. 3 of [12] exhibits the same features. We focus on the effect of the applied boundary conditions in comparison with a scanning electron microscope compression experiments, see [25] and the numerical studies [12, 25]. The presented approach allows for nominal strains up to $0.5$ and overcomes numerical difficulties encountered in [12, 25], where the authors report results for the nominal strain up to $0.2$. Nevertheless, for nominal strains higher than $0.5$ a collapse of our calculations appeared. The amount of mesh distortion that can be handled by the ALE method is limited and a full remeshing of the problem is eventually required to overcome this limit.

3.3 Activity of slip systems

In Figure 4 we present the contour maps of slip on the primary system ($\boldsymbol{s}^{(2)}$) and compare the results for two kinds of boundary conditions. The sample deformation and the primary slip activity exhibit the same features as the one shown in Figure 3 of [12], in the nominal strain region up to 0.2. There is a measurable difference between the alignment of slip activity regions (shear bands) between different types of boundary conditions, see Figure 4. The majority of deformation has been localized to an intense slip band. Slips are highly concentrated in a shear zone.

In Figure 5 we present the average slip rate for each of twelve individual slip systems, and compare two kinds of boundary conditions. The slip rates are not qualitatively similar for different types of boundary conditions. While the constrained boundary conditions give rise to stable constant-in-time values, the unconstrained case changes the behaviour for strains over $0.05$. In both cases we observe instabilities for high strains over 0.3, see Figure 4. The inactive slip systems $\alpha\in\{7,8,9\}$ are oriented in such a way that the slip direction $\boldsymbol{s}^{(\alpha)}$ is perpendicular to the compression axis, namely $\boldsymbol{m}^{(\alpha)}\cdot(-1,2,3)=0$, see Table 1.

3.4 Stress-strain behaviour

Figure 6 shows the nominal compressive stress-nominal compressive strain curves. Despite of the considered zero hardening, the constrained condition exhibited significant stress increase while the unconstrained condition showed a decrease of the calculated nominal compressive stress. The stress increase in the constrained case is caused by the Dirichlet velocity boundary condition on the top and bottom boundaries. There, plastic deformation is suppressed, and part of the sample deformed elastically causing the stress to increase. The effect is associated with barrelling of the sample. On the other hand, the decrease of the nominal stress with the increasing nominal strain in the case of the unconstrained boundary condition is associated with the inclination of the sample. As is seen from the Figure 6 the results in the strain range up to $0.2$ agree qualitatively with the numerical results presented in [12, Fig.11] and [25, Fig. 3(A)], and the experimental results [25, Fig.4]. However, to obtain an exact agreement with the experiments would require the application of boundary conditions that exactly mimic the experimental setting; this is however beyond of the scope of this work.

4 Summary

We presented a new approach to the derivation of a thermodynamically compatible rate-type fluid model for crystal plastic materials. We employed a Gibbs-potential-based approach and derived a rate-type stress-strain constitutive equation. Instead of the standard crystal plasticity hypothesis that the plastic part of the velocity gradient equals a sum of contributions of the individual slip systems $\boldsymbol{L_{p}}=\sum_{\alpha=1}^{N}\nu^{(\alpha)}\boldsymbol{s}^{(\alpha)}\otimes\boldsymbol{m}^{(\alpha)}$ we have employed the procedure of maximization of the rate of dissipation to derive this relation from a single scalar function, (29). It is worth emphasizing that the Gibbs-potential-based approach, as developed in [21] and presented here, automatically provides an objective rate derivative for the stress as well as the possibility to incorporate anisotropy of the material response.

• 2.

We have incorporated an Eulerian form for the evolution of the lattice basis vectors instead of computing directly the evolution of the slip directions and the normals to the slip plane. This is one of the novel ingredients of the developed model. It allows us to perform efficient three-dimensional computations of the material with a full set of slip systems (e.g. 12 slip systems in the case of FCC structure), as illustrated on two test problems concerning three-dimensional micropillar compression.

• 3.

We have formulated a finite element discretization scheme, which has been used to solve the fully coupled problem. Our solver is monolithic. The Arbitrary Lagrangian Eulerian (ALE) approach has been applied in order to use the Eulerian formulation on a moving mesh, that captures the free boundary. We developed a new algorithm for a three-dimensional compression and achieved qualitative correspondence between our results and the experimental and numerical results concerning the same problem, as presented in [12].

• 4.

The fact that our method is based on the Eulerian formulation, and consequently the primal role is given to the velocity of the material and the distortion of the crystal lattice space, enabled us to demonstrate that the whole approach is capable of simulating deformation processes with large strains for suitably defined problems.

Appendix

The ALE approach is a finite element formulation in which the computational system is not a priori fixed in space (Eulerian) or attached to the material (Lagrangian) [7, 24]. In fact we can formulate the problem on an arbitrary moving, time-dependent domain and define the ALE mapping between a fixed computational domain and an arbitrary domain. The ALE mapping is the key to transforming our system to a fixed computational domain. To this end the proper definition of the ALE mapping is of a great importance. The class of ALE methods is widely applied, e.g. to fluid-structure interaction problems. In crystal plasticity, the ALE method has been successfully tested in the solution of a two-dimensional problem of a channel-die compression in [2, 4, 3], where the 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 finite element discretization is used in space, and time is discretized by the backward Euler scheme. The choice of the conforming simplical tetrahedral elements for each variable consists of $\mathcal{P}_{2}$ elements for the velocity, $\mathcal{P}_{1}$ for the density and the lattice basis vectors, and $\mathcal{P}_{1}-disc$ (discontinuous elements) for the Kirchhoff stress. We employ the ALE approach in the sense that we use our Eulerian formulation on a moving mesh, which captures the free boundary.

To solve (45) at each time step, we perform two actions. Given $\boldsymbol{v}^{k}$, $\varrho^{k}$, $\boldsymbol{S}^{k}$, $\boldsymbol{a_{i}}^{k}$ from the previous time level we proceed as follows:

Step 1: Solve problem for $\boldsymbol{v}_{h}^{k+1},\varrho_{h}^{k+1},\boldsymbol{S}_{h}^{k+1},$ and $\boldsymbol{a_{i}}_{h}^{k+1}$ ($i\in\{1,2,3\}$)

[TABLE]

where

[TABLE]

Step 2: Move the mesh by $\boldsymbol{u}_{h}^{k}=\boldsymbol{v}^{k}_{h-\mathrm{mesh}}\Delta t$. The velocity of the mesh motion $\boldsymbol{v}^{k}_{h-\mathrm{mesh}}$ is taken as the material velocity $\boldsymbol{v}_{h}^{k}$, which in fact leads to an updated Lagrangian description.

The mesh velocity in the second step can be chosen arbitrarily, so that $\boldsymbol{v}^{k}_{h-\mathrm{mesh}}|\partial\Omega=\boldsymbol{v}^{k}_{h}|\partial\Omega$, 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 number 2012/07/N/ST1/03369]. The authors thank Martin Kružík for general discussions on this topic, Jaroslav Hron for specific discussions on numerics-related issues and Endre Süli for several suggestions improving the final form of the text.