## Abstract

A micromechanical model for finite single crystal plasticity was introduced by Kochmann & Hackl (2011 *Contin. Mech. Thermodyn.* 23, 63–85 (doi:10.1007/s00161-010-0714-5)). This model is based on thermodynamic variational principles and leads to a non-convex variational problem. Based on the Lagrange functional, an incremental strategy was outlined to model the time-continuous evolution of a first-order laminate microstructure. Although this model provides interesting results on the material point level, owing to the global minimization in the evolution equations, the calculation time and numerical instabilities may cause problems when applying this model to macroscopic specimens. In this paper, a smooth transition zone between the laminates is introduced to avoid global minimization, which makes the numerical calculations cumbersome compared with the model in Kochmann & Hackl. By introducing a smooth viscous transition zone, the dissipation potential and its numerical treatment have to be adapted. We outline rate-dependent time-evolution equations for the internal variables based on variational techniques and show as first examples single-slip shear and tension/compression tests.

## 1. Introduction

The concept of crystal plasticity was first introduced by Ewing & Rosenhain in 1899 [1]. They observed that slip steps on the surface of metals under plastic deformation arise owing to slip bands inside the metal. They thus concluded that plastic deformations appear owing to shearing of certain crystal planes. This evolving microstructure (e.g. figure 1) strongly influences the macroscopic behaviour and hence it has to be transferred to appropriate simulation schemes by micromechanical modelling. Metals and alloys can be interpreted on the microstructural level as a crystal that exhibits a periodic cell structure. Therefore, plastic deformation on the microlevel can be examined within the framework of a single crystal. Detailed modelling of this complex behaviour requires much more sophisticated material models that establish the class of the so-called crystal plasticity models.

In the context of crystal plasticity, dislocations within the crystal undergo a gliding process during loading. This plastic deformation is irreversible because parts of a crystal glide along a glide plane (denoted by the normal vector **m**) in a slip direction (vector **s**). The gliding process in the crystal results from the movement of a huge number of dislocations throughout the crystal lattice. It is thus obvious that models for, e.g., perfect plasticity do not take account of the special character of crystal plasticity. Consequently, models need to be developed that resolve these complex microstructural processes with a higher degree of detail. There is already a variety of models dedicated to crystal plasticity that describe this observed elasto-plastic behaviour of a single crystal and which use different modelling approaches. For instance, we refer to [3–5]. The smallest length scale of these models is the atomistic scale. One of the atomistic approaches is molecular dynamic simulation, which discretizes crystal atoms in the sample (e.g. [3] or [6]). These models require a description of the interatomic interactions. Unfortunately, important drawbacks of these simulations are the limitations in time and space.

Another approach to modelling crystal plasticity is based on the continuum dislocation theory. Instead of modelling discrete atoms, the huge number of dislocations in the crystal is treated as a continuous quantity, which is described by the dislocation density. This dislocation network is dictated by the energy of the crystal. The continuum dislocation theory was first introduced by Kröner [7], Bilby [8] and Kondo [4]. In the case of the continuum dislocation theory for finite strain, there is still no commonly accepted definition for the dislocation density, which is the key kinematic quantity.

The last method of modelling crystal plasticity that we mention here is based on variational analysis. For this modelling concept, only two quantities have to be defined specifically for each material and process: the Helmholtz free energy and the dissipation functional. Owing to this very limited number of necessary modelling assumptions, phenomenological aspects stay minimal. In this paper, we aim to contribute to the modelling of crystal plasticity by presenting a viscous, time-dependent approach for the above-described processes, derived using the principle of the minimum of the dissipation potential. The evolving microstructure is investigated as a minimizer of a non-convex energy potential, specifically of the dissipation potential. For numerical efficiency, we restrict ourselves here to first-order laminates that are formed by the material. Of course, the order of laminates may be far higher in reality. However, application of a first-order laminate as a chosen minimizer leads to a rank-one convex approximation of the quasi-convex energy hull, which was introduced in [9].

The method of relaxing the energy of the crystal was concluded first by Ortiz & Repetto [5] and by Carstensen *et al.* [10]. The idea is that if a homogeneous material has a non-quasi-convex energy, the energy minimum is archived by the formation of microstructural patterns. Therefore, there are several models that employ this theory of relaxation to derive the incremental evolution of microstructures (see for example Oritz & Repetto [5] or Lambrecht *et al.* [11]). These models use a condensed energy functional, which may lead to some difficulties [12] because these models need explicit dissipation distances which are quite often not feasible. They also do not take account of all changes in an existing microstructure during each time step.

Kochmann and Hackl developed an incremental solution formulation based on the rank-one convex hull to simulate the evolution of the rank-one laminates in a single crystal [9]. This model captures the microstructural characteristics (see [13,14]). However, a drawback of this model is that the evolution equations for the microstructure include global minimizations in every time step. This yields very cumbersome numerical calculations and may also lead to unstable results. Hence, we expand this model in our new approach, which includes the introduction of a vanishing viscosity. This viscosity can be interpreted as a small and smooth transition zone between the laminates. By adapting the dissipation potential, we derive evolution equations without global minimizations which is a huge improvement compared with the latter model, because the calculations in our model are simplified. This reduces the computational effort and yields stable numerical results. After derivation of our modelling equations, we perform virtual experiments on a material point level for shear and tension/compression tests.

Klusemann & Kochmann published in [15] a conceptually different approach: in contrast to our ansatz of a modified dissipation, they introduced a interfacial energy contribution. However, both approaches share the characteristic of a delayed microstructural evolution. As already stated, our intention is a stabilization of the numerical behaviour. Interestingly, Klusemann and Kochmann observed that a more sophisticated numerical handling was necessary for their expanded model.

## 2. Variational framework

Our material model is based on the principle of the minimum of the dissipation potential. Details of this principle are given in references [16–22]. Several material models were designed using this variational framework, e.g. those in references [23–26]. The principle of the minimum of the dissipation potential for non-isothermal processes was presented in reference [27] and applied to shape memory alloys in reference [28].

As pointed out previously, only the Helmholtz free energy and the dissipation functional have to be specified independently. The key idea of that modelling technique is the introduction of a Lagrangian **u** is the displacement field and **z** is the tensorial generalized internal variable which carries information about the current microstructural composition of the material. Minimization of the Lagrangian yields the evolution equations for **z** without further assumptions,

The actual displacement field **u** is obtained from the principle of minimum potential,

with the total free energy of the body
*t*,**u**) represents the potential of external forces.

Application of these two minimum principles for the evolution of microstructures was first introduced by Ortiz & Repetto [5], Carstensen *et al.* [10] and Mielke [29]. For non-convex energy densities, no homogeneous deformation can be found as a minimizer, such that the previously homogeneous material forms patterns that instead reduce the energy of the crystal. The non-convex energy density is replaced by a quasi-convex hull, and the admissible microstructures are pre-accounted. These microstructures minimize the energy (see [30,31]). Because computation of the quasi-convex hull is not feasible, it is replaced here by a rank-one convex hull. The minimizing microstructure can be interpreted as first-order laminates for the rank-one convex hull [32].

## 3. Lamination

The quasi-convex hull is very complex to compute. We thus restrict ourselves to the rank-one convex energy hull for which the resulting microstructure equals a first-order laminate, as illustrated in figure 2. The phases of the laminate are described by chosen internal variables, here the volume fractions for each laminate λ_{i}, the plastic slips in each laminate *γ*_{i}, the angle of the laminate *ϕ*, and the respective hardening parameters *p*_{i}. The normal vector **b** of the laminate is expressed in terms of the angle *ϕ*.

The restriction to the formation of first-order laminates allows a kinematic simplification, because each domain of the laminate has a constant deformation gradient **F**_{i} in this case which can be defined by
**a**_{i} represent some deformation amplitude vectors.

By imposing the volume average of the deformation and ensuring the incompressibility condition of each domain,

For finite elastoplasticity, the deformation gradient **F** can be decomposed multiplicatively into its elastic and plastic contributions

The plastic deformation takes place owing to gliding processes of dislocations along an active slip system. A slip system is defined by the unit vectors **s** and **m** which denote the slip direction and the normal vector of the glide plane, respectively. For *n* slip planes, the plastic flow rule is given by
*γ* is the total amount of plastic slip within the single-slip system assumed. It thus may be calculated by simple volume averaging:

An additional flow rule (e.g. Carstensen *et al.* [10]) reads
*p*(0)=0. The hardening parameter *p* is responsible for the plastic contribution to the energy.

## 4. Material model

In this paper, we examine the evolution of microstructures in an incompressible neo-Hookean material, which can be described by the elastic free energy density
*μ* is the shear modulus and *κ* is the hardening modulus. As proposed earlier, **F**_{e} denotes the elastic deformation gradient. The parameter *α* influences the hardening behaviour of the material and is commonly set to 2, which describes linear hardening, or 4. For our calculations, we used *α*=4.

For the variational modelling, the relaxed energy (hull) of the material has to be used which was given in reference [33]. The relaxed energy was found via minimization of the weighted sum of the individual energy contributions from the respective laminates under consideration of two constraints: first, the laminates have to form compatible interfaces, expressed by
*N* denotes the number of phases in the laminate and **a**_{i} is an amplitude vector which results from the Hadamard condition. Second, each laminate is assumed to be incompressible, yielding equation (3.2). For the two constraints, Lagrange multipliers ** Λ** and

*ρ*

_{i}were introduced. Then, the relaxed energy hull reads

*Ψ*

^{rel}, we refer to [9] and recall here only the result, which is

**b**

_{i}of each phase is defined as

**C**=

**F**

^{T}⋅

**F**.

For simplicity, we assume a two-phase laminate in the first examination. Therefore, *N*=2 is set in the following. The volume fraction in phase one is consequently defined as 1−λ, and that in phase two is λ.

The intention of our model is to expand of the model in reference [33] by viscous effects. In this case, the dissipation functional considered has to be adapted in a way that the stationarity conditions of the Lagrangian *s* is the viscous parameter of a transition zone and *r* is the critical resolved shear stress. The first term is a rate-independent ansatz for the dissipation (see [10]). The second term is rate-dependent. For the two-phase laminate, the relaxed dissipation potential, based on (4.6) also has to be found. The relaxed dissipation potential owing to the change of plastic slip can be easily found by weighting the dissipation functional for each phase, given in equation (4.6), with the respective volume fraction. The summation yields the total relaxed dissipation functional owing to plastic slip as
*v* in the normal direction and is defined by the width *δ* (volume ratio) and the same viscous parameter *s* as used before. The velocity can then be expressed in terms of the rate of λ by

## 5. Numerical treatment

The evolution equations for the internal variables are found by an analysis of the stationarity conditions of *γ*_{i}, these read
*Ψ*^{rel}/∂**z** constitutes thermodynamic forces. Thus, the second law dictates that the thermodynamic fluxes always have the same sign as their respective thermodynamic forces. Employing this identity in combination with a Legendre transformation of the dissipation functional, the evolution equations can be written in terms of
*r*|*γ*_{1}−*γ*_{2}| and *r*λ_{i}, respectively). The plus sign indicates that an evolution may only take place when the respective yield function is positive. Each rate is weighted by some viscous effects. A detailed derivation of this general structure of evolution equations can be found for instance in [34].

Once a laminate exists, a rotation of the laminate is possible when energetically favourable. We check whether the amount of dissipation is smaller than the driving force, that is if
*p*_{i} has to be updated twice according to [9]. The first update is caused by the change of the volume fraction, which causes a hardening. The updated values are obtained by kinematic averaging analogously to [15] (to be precise, the same update scheme is applied for updating the accumulated plastic strains in reference [15] which differs from our approach here in which we use same update scheme for the hardening parameters)
_{n+1}=λ_{n}+Δλ. For Δλ<0, the update rule according to
*t*,

As also discussed in the recent publication of Klusemann & Kochmann [15], the update rule for the hardening parameter is questionable from a thermodynamical point of view because it underestimates the stored energy. However, the principal goal of this work is a numerical simplification of the work in reference [33]. Therefore, we do not change this part of the model, also owing to its kinematic justification.

To conclude, the equations ((5.3), (5.4), (5.7)) in combination with (5.8) or (5.9), and (5.11) and (5.12) form the final system of model equations that have to be solved, e.g. by employing a second order Runge–Kutta integration scheme for simplicity. However, a closer inspection of these equations reveals that an numerically motivated initiation of the microstructure in terms of an initial volume fraction and plastic slip is necessary. Therefore, as long as no microstructure exists, hence λ=0∨λ=1, testing values of *Ψ*^{rel}/∂λ>0 if λ=0 and sign−∂*Ψ*^{rel}/∂λ<0 if λ=1, respectively).

## 6. Numerical results

The presented material model may be applied to an arbitrary deformation. For a first analysis, we present numerical results for simple shear and tension/compression tests which are identical to the considered problems in reference [33]. Let us begin with a simple shear test, as schematically shown in figure 4. The angle *α* describes the given slip system, such that the slip direction can be obtained by *α*=135° in analogy to [33]. For a simple shear test, the given deformation gradient reads
*μ*=2 GPa (shear modulus), the hardening modulus *κ*=0.01 GPa and the critical resolved shear stress is chosen as *r*=0.001 GPa. For the viscous parameters, we choose *s*=*r***t*^{⋆} (*t*^{⋆}=0.01 s) and *δ*=*r***j* mm (*j*=0.1 mm^{3} N^{−1}).

The type of evolution equations are rate-dependent. Thus, they have to be discretized in time which is performed by introducing a time increment Δ*t*=1×10^{−4} s. The load is also discretized in terms of load increments Δ*w*=1×10^{−3} [−]. Therefore, a loading velocity *v*:=Δ*w*/Δ*t* may be defined on which the viscous contribution depends: higher-loading velocities result in larger viscous effects, slower-loading velocities result in the opposite. The following examples investigate this characteristic. The initial loading velocity is termed by *v*_{0}.

On the left-hand side of figure 5, the volume fraction of the second phase of the laminate is presented (λ). When the material is loaded, at first there is no evolving microstructure, as expected: lamination will set in only when a critical external load has been applied to the material. In our example, this is the case at approximately *w*=0.10 when the volume fraction of phase two starts to evolve. At approximately *w*=2.35, a uniform microstructure has established that only consists of the second phase of the laminate.

Simultaneously with the evolution of laminate, plastic slip evolves in both phases. The plastic slip in phase two is negative and jumps from zero to *γ*_{2}=−0.15. During further loading, its absolute value increases slowly. The plastic slip in phase one is positive and starts later, at *w*=0.35.

Of course, these internal variables strongly influence the behaviour of the stresses and the energy. The jump of *γ*_{2} may seem surprising at first glance. However, we propose a material model that also accounts for viscous effects. This, in turn, delays the evolution of plastic slip to some extent, so that we observe a very drastic evolution of plastic deformation first when the associate laminate phase λ establishes. Although the plastic slip is very huge within that phase, the total amount of plastic deformation, that is (1−λ)*γ*_{1}+λ*γ*_{2}, presented in figure 6, remains very limited (and smooth). Owing to the high external deformation that is necessary to initialize evolution of a laminate, a huge plastic deformation within the newly establish laminate is physically sound. The microstructure is also described by the angle of the laminate, presented on the left-hand side of figure 7. At the beginning of loading, the angle of the laminate is zero, which is not surprising, because the entire material is homogeneous. At the onset of lamination (at *w*=0.10), the angle starts to evolve: after a sharp jump and overshoot when the microstructures initiate to evolve, only a slight evolution of *ϕ* is observed.

The right-hand side of figure 7 shows the hardening parameters of both phases. Because the plastic slip in phase two evolves from the beginning of the lamination, there is also increasing hardening owing to the coupled evolution of hardening to the evolution of plastic slip. The hardening increases monotonously. In phase one, it is delayed and onsets when the plastic slip evolves in phase one. When phase one no longer exists, no further evolution of plastic slip nor hardening takes place in that phase. Therefore, *γ*_{1} and *p*_{1} remain constant from here on. Along with the evolving microstruture, we calculate the Cauchy stresses by
*w*>2.35, there exists no further material that could transform to phase two. Consequently, the Maxwell line is left and stresses start to increase again.

Before the microstructure initiates, the viscous zone has to be established. Therefore, minimization of the energy by the microstructure is delayed, such as the evolution of plastic slip in phase one. This is also visible in the energy plot in the right-hand side of figure 8. Until the laminate evolves, the energy increases. While the microstructure adapts to the external loading in an energy-minimizing way, the energy decreases.

We investigate the influence of the viscous effects by increasing the relaxation time, e.g. by reducing the loading velocity to *v*=0.2*v*_{0} and recalculate the shear test. All other material parameters are kept constant.

The volume fraction is presented in figure 9 at the left-hand side. Comparison with the higher-loading velocity in figure 5 reveals that now the slope of the curve is slightly increased which implies an earlier completion of the phase transformation. The evolution of plastic slip in phase two is less affected by the increased relaxation time (right-hand side of figure 9), whereas the plastic slip in phase one now onsets earlier, too. Interestingly, *γ*_{1} decreases for *w*>1.4 for the increased relaxation time. The curve of the total amount of plastic slip, presented in figure 10, also remains unaffected, because the volume fraction of phase one is too small when changes in the curve for *γ*_{1} become present.

The angle of the laminate and of the hardening are presented in figure 11. The angle follows approximately the same path as before. Comparable to the behaviour of the plastic slips, the hardening variables remain either unmodified (*p*_{2}) or show a different behaviour than before (*p*_{1}). The increase in time for the material to relax, along with the resulting ‘intensified’ microstructure, strongly influences stress and energy: the viscous peaks in the stresses and in the energy are not as high and are less distinct (figure 12). The peak in the stress for the initial pseudo-velocity has reached up to 0.65 GPa, whereas for a reduced pseudo-velocity that is only 20% of its initial value, the stress peak is only 0.25 GPa. This indicates that a viscous treatment which allows the material to relax leads to a proper solution in which the energy tends to be convex, right-hand side figure 12. The second stress drop at the end of the plateau seems to be surprising at first glance. However, same observation was made, e.g. in [36] where non-convex strain gradient plasticity was investigated.

The next presented example is a tension/compression test, such as that depicted in figure 13. We choose a slip system with an angle of *α*=70°. Furthermore, we use the same material constants as in the previous example except of *κ*=0.01 which is, again, chosen in analogy to reference [33]. The deformation gradient for this tension/compression test reads

In this test, initiation of microstructure begins almost directly (figure 14): the volume fraction increases rapidly up to a value of 30% before the slope is less pronounced. At *w*=1.75, phase transformation is completed and λ=1.

The plastic slip of phase two jumps from zero to −0.15 (left-hand side in figure 15). The plastic slip then increases up to −2.00 during formation of laminate phase two. The jump in the beginning of the plastic slip in phase two is induced, in analogy to the shear test, by a viscous effect owing to the delay of the evolution. The total amount of plastic slip (right-hand side of figure 15) decreases slowly and smoothly from zero to −2.5.

The last internal variables that describe the laminate are the angle of the laminate and the hardening parameters, presented in figure 16. The angle jumps directly at the onset of lamination to its final value. No further evolution can be observed. The hardening parameters are plotted on the right-hand side of figure 16 and show a quite interesting behaviour. They evolve in parallel with nearly the same magnitude. This specific effect results from the evolution of the plastic slip: the absolute values of *γ*_{1} and *γ*_{2} correspond to each other and only differ in their sign. Because the update of the plastic slip has the most prominent impact on the evolution of hardening, *p*_{1} and *p*_{2} are quite identical.

Figure 17 shows the resulting normal stresses *σ*_{11} and *σ*_{22}, and the energy is plotted in figure 18. The normal stress *σ*_{11} is presented on the left-hand side of figure 17. During the first loading steps, the stress increases up to a maximum value of 0.12 GPa. In these first loading steps, the volume fraction of phase two is negligible: even there is a large amount of plastic slip in phase one, the total amount of plastic slip is too small to decrease the stress. With the increase of volume fraction and the plastic slip in phase one, the stress drops drastically to 0.00 GPa and then remains constant. The *σ*_{22}-stress shows a behaviour mirrored to the abscissa. At first glance, the zero stresses during the evolution of microstructure are surprising because hardening is present. However, the hardening constitutes as an ‘inner’ hardening in each laminate (figure 16). On the other hand, plastic slip is identical in both phases, except of the sign. Therefore, a ‘global’ hardening effect is removed (see also equation 3.5).

Figure 18 shows the energy of this tension/compression test. The behaviour is similar to the energy of the shear test. Before the microstructure can evolve, the viscous zone has to be established and minimization of the energy is delayed. While the laminate evolves, it is still not adapted to the external loading, hence the energy is not minimized and increases. Thereafter, the energy is minimized by the microstructure and decreases. Thus, the energy is still non-convex.

If we change the loading velocity from *v*=*v*_{0} to *v*=0.2*v*_{0}, as in the shear test, the material has more time to relax and hence viscous effects are less intense. The previous example with higher loading velocity showed already minor viscous effects. Therefore, the evolution of microstructure is affected only to a small degree for *v*=0.2*v*_{0} (figures 19–21). However, the stress drops are further reduced (figure 22) along with the non-convexity of energy in figure 23 at *w*=0.1.

The presented results may be compared with the results in reference [33] where a quasi-static microstructure evolution by means of energy minimization was presented for the shear and the tension/compression test. There, also the results for the unrelaxed energy are investigated. The new model presented here yields results that are not completely identical to the ones in [33]. This, however, is expected: owing to the viscous behaviour of the new model, both the onset and the evolution of microstructure are delayed. For the shear test, a comparable curve both for the volume fractions (note that both λ=0 and λ=1 denote a homogeneous material) and the plastic slips is computed. However, the energy obtained by the new model is smaller than for the old one—which is physically sound owing to the viscous extension of the dissipation functional. However, there are differences for the stresses which are quite sensitive to the precise values of the internal variables.

For the tension/compression test, in the old model, the plastic slips possessed a larger difference than predicted by the new model. Thus, small hardening effects are present in the stress curve for the old model which are not observed in the new model. The argument of increased dissipated energy holds even more for this loading case.

## 7. Conclusion and outlook

This paper introduces a new approach for modelling crystal plasticity based on variational concepts. The approach is based on the addition of a vanishing viscosity to the existing model of Kochmann & Hackl [33]. Consequently, we modified the dissipation functional and derived explicit evolution equations for the internal variables that describe the underlying microstructure. These evolution equations no longer include any global minimization; however, they can be evaluated directly with standard approaches for numerical integration. This might be more beneficial in comparison with [15] where a modification of the energy yields to a more complex numerical handling. We presented first examples for numerical experiments, specifically simple shear and tension/compression tests. Both examples show that a microstructure evolves in order to minimize the energy. Introducing a vanishing viscosity yields a peak in the respective stresses, which can be reduced by decreasing the pseudo-velocity.

Owing to its stable general character, our material model may be implemented and evaluated in a finite-element framework. Artificial assumptions for the present loading state are thus obsolete and a direct comparison of the numerical results with experimental observations is feasible.

## Computational solution

The model is implemented in MATHEMATICA 8. The tolerances are considered as 10^{−8} and also an equivalent constant step-size is taken. Second-order Runge-Kutta method is used for the updating procedure.

## Competing interests

We declare we have no competing interests.

## Author contributions

K.H. conceived the mechanical model; C.G. and P.J. interpreted the computational results, and wrote the paper. C.G. implemented and performed most of the simulations. All authors gave final approval for publication.

## Funding

This work was supported by research group FOR 797 (Microplast).

- Received February 23, 2015.
- Accepted July 6, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.