**Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals**

Anxin Ma and Alexander Hartmaier

*Interdisciplinary Centre for Advanced Materials Simulation, Ruhr-University Bochum Germany*

### **1. Introduction**

Plastic deformation of polycrystalline materials includes dislocation slip, twinning, grain boundary sliding and eigenstrain produced by phase transformations and diffusion. These mechanisms are often alternative and competing in different loading conditions described by stress level, strain rate and temperature. Modelling of plasticity in polycrystalline materials has a clear multiscale character, such that plastic deformation has been widely studied on the macro-scale by the finite element methods, on the meso-scale by representative volume element approaches, on the micro-scale by dislocation dynamics methods and on the atomic scale by molecular dynamics simulations. Advancement and further improvement of the reliability of macro-scale constitutive models is expected to originate from developments at microstructural or even smaller length scales by transfering the observed mechanisms to the macro-scale in a suited manner. Currently efficient modelling tools have been developed for different length scales and there still exists a challenge in passing relevant information between models on different scales. This chapter aims at overviewing the current stage of modelling tools at different length scales, discussing the possible approaches to bridge different length scales, and reporting successful multiscale modelling applications.

Fig. 1. Multiphase polycrystalline RVE (right) with 90% matrix and 10% precipitate. The grain size has a normal distribution (middle) and the [111] polfigure (left) shows a random texture.

Because **v**� is a arbitrary vector in the [**x**�

individual grain with coordinate

**3.1 Kinematics**

deformation gradient tensor **F**

twinned material point with coordinate

(2*<sup>k</sup>* <sup>−</sup> <sup>1</sup>) <sup>−</sup> <sup>1</sup>

2 · *<sup>d</sup>* <sup>&</sup>lt; *z*� <sup>1</sup> − *z*� 0 ≤ <sup>2</sup>*<sup>k</sup>* <sup>−</sup> <sup>1</sup> 2 

deformation is accommodated by dislocation slip.

**3. Constitutive models based on continuum mechanics**

**<sup>F</sup>** <sup>=</sup> *<sup>∂</sup>***<sup>x</sup>**

can be derived for the elastic and the plastic deformation gradients as

*<sup>∂</sup>***<sup>X</sup>** <sup>=</sup> *<sup>∂</sup>***<sup>x</sup>** *∂***x** *∂***x**

Otherwise orientation **Q**<sup>I</sup> will be assigned to this material point.

and we find the relationship between **Q**<sup>I</sup> and **Q**II as

, **y**� , **z**�

**Q**II = **R**<sup>T</sup>

**x**� <sup>1</sup>, **y**� <sup>1</sup>, **z**� 1 

the grain center along the habit plane normal direction and the lamella thickness *d* satisfy

Fig. 2. The multiplicative decomposition of the deformation gradient where the plastic

For large strains the elastic and plastic deformation can be separated consistently. We follow the well-known multiplicative decomposition proposed by Lee (1969) (see Figure 2) of the

where **F**e is the elastic part comprising the stretch and rotation of the lattice, and **F**p corresponds to the plastic deformation. The lattice rotation **R**e and stretch **U**e are included in the mapping **F**e. They can be calculated by the polar decomposition **F**<sup>e</sup> = **R**e**U**e, i.e., the texture evolution is included in this part of the deformation. Furthermore, two rate equations

<sup>L</sup> = **R**<sup>M</sup>

Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals 5

L**R**-1

The crystal orientation **Q**<sup>I</sup> has been assigned to the material point at the centre of the

**R**L**Q**II**R**<sup>T</sup> L 

**R**L**Q**I**R**<sup>T</sup>

**x**� <sup>0</sup>, **y**� <sup>0</sup>, **z**� 0  ] coordinate system equation 3 will reduce to

<sup>M</sup>**R**L**Q**I. (5)

. The orientation **Q**II will be assigned to the

when the distance between this point and

*<sup>∂</sup>***<sup>X</sup>** <sup>=</sup> **<sup>F</sup>**e**F**<sup>p</sup> (7)

**<sup>F</sup>**˙ <sup>e</sup> <sup>=</sup> **LF**<sup>e</sup> <sup>−</sup> **<sup>F</sup>**e**L**<sup>p</sup> (8)

**F**˙ <sup>p</sup> = **L**p**F**<sup>p</sup> (9)

· *d* with *k* = 1, 2, 3, ... (6)

(4)

#### **2. Generating realistic material microstructures**

The current advanced high strength steels (AHSS) such as dual phase steels, transformation induced plasticity (TRIP) steels, twin induced plasticity (TWIP) steels and martensitic steels are all multiphase polycrystalline materials. In order to model the the macroscopic mechanical properties such as yield stress, work hardening rate and elongation to fracture, one has to build a representative volume element (RVE) for each macroscale material point and investigate the local deformation of each material point within the RVE, and then make a volume average. In this micro-macro-transition procedure, in order to reduce the computational costs the statistically similar representative volume elements (SSRVEs) have been developed to replace real microstructures from metallurgical images by Schröder et al. (2010).

Considering the real microstructure of multiphase materials, during the representative volume element generation one should consider grain shape distribution, crystalline orientation distribution, grain boundary misorientation angle distribution and volume fraction of different phases. Figure 1 is an example of the RVE we have generated for TRIP steels where the Voronoi tessellation algorithm has been used.

Recent studies (Lu et al., 2009; 2004) show bulk specimens comprising nanometer sized grains with embedded lamella with coherent, thermal and mechanical stable twin boundaries exhibiting high strength and considerable ductility at the same time. These materials have higher loading rate sensitivity, better tolerance to fatigue crack initiation, and greater resistance to deformation. Under this condition, the RVE with nanometer sized twin lamella inside nanometer sized twin lamella inside nanometer sized grains will help us to understand existing material behavior and design new materials.

Assume two orientations **Q**<sup>I</sup> and **Q**II have the twin relationship in (1, 1, 1) habit plane along [1, 1, −2] twinning direction. For any vector **V**, these two tensors will map as **v**<sup>I</sup> = **Q**I**V** and **v**II = **Q**II**V**. The twin relationship between **v**<sup>I</sup> and **v**II is easier to see in the local twin coordinate system with **x**� //[1, 1, −2] and **z**� //[1, 1, 1] rather than in the global coordinate system [**x**, **y**, **z**]. We define a orthogonal tensor **R**<sup>L</sup> for the mapping from global coordinate system to the local coordinate system

$$R\_{\rm Li/j} = \begin{bmatrix} \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{3}}\\ \frac{1}{\sqrt{6}} & \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{3}}\\ \frac{-2}{\sqrt{6}} & 0 & \frac{1}{\sqrt{3}} \end{bmatrix} \tag{1}$$

and another tensor **R**<sup>M</sup> for the mirror symmetry operation

$$R\_{\mathbf{M}ij} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1 \end{bmatrix} \tag{2}$$

and get the twin relationship in the local twin coordinate system

$$\left(\mathbf{R}\_{\mathrm{L}}\mathbf{Q}\_{\mathrm{I}}\mathbf{R}\_{\mathrm{L}}^{\mathrm{T}}\right)\mathbf{v}' = \mathbf{R}\_{\mathrm{M}}\left(\mathbf{R}\_{\mathrm{L}}\mathbf{Q}\_{\mathrm{II}}\mathbf{R}\_{\mathrm{L}}^{\mathrm{T}}\right)\mathbf{v}'.\tag{3}$$

Because **v**� is a arbitrary vector in the [**x**� , **y**� , **z**� ] coordinate system equation 3 will reduce to

$$\mathbf{R}\_{\rm L} \mathbf{Q}\_{\rm I} \mathbf{R}\_{\rm L}^{\rm T} = \mathbf{R}\_{\rm M} \left( \mathbf{R}\_{\rm L} \mathbf{Q}\_{\rm II} \mathbf{R}\_{\rm L}^{\rm T} \right) \tag{4}$$

and we find the relationship between **Q**<sup>I</sup> and **Q**II as

$$\mathbf{Q}\_{\rm II} = \mathbf{R}\_{\rm L}^{\rm T} \mathbf{R}\_{\rm M}^{-1} \mathbf{R}\_{\rm L} \mathbf{Q}\_{\rm I}. \tag{5}$$

The crystal orientation **Q**<sup>I</sup> has been assigned to the material point at the centre of the individual grain with coordinate **x**� <sup>0</sup>, **y**� <sup>0</sup>, **z**� 0 . The orientation **Q**II will be assigned to the twinned material point with coordinate **x**� <sup>1</sup>, **y**� <sup>1</sup>, **z**� 1 when the distance between this point and the grain center along the habit plane normal direction and the lamella thickness *d* satisfy

$$\left[\left(2k-1\right)-\frac{1}{2}\right] \cdot d < \left|z\_1' - z\_0'\right| \le \left[2k - \frac{1}{2}\right] \cdot d \text{ with } k = 1, 2, 3, \dots \tag{6}$$

Otherwise orientation **Q**<sup>I</sup> will be assigned to this material point.

Fig. 2. The multiplicative decomposition of the deformation gradient where the plastic deformation is accommodated by dislocation slip.

#### **3. Constitutive models based on continuum mechanics**

#### **3.1 Kinematics**

2 Will-be-set-by-IN-TECH

The current advanced high strength steels (AHSS) such as dual phase steels, transformation induced plasticity (TRIP) steels, twin induced plasticity (TWIP) steels and martensitic steels are all multiphase polycrystalline materials. In order to model the the macroscopic mechanical properties such as yield stress, work hardening rate and elongation to fracture, one has to build a representative volume element (RVE) for each macroscale material point and investigate the local deformation of each material point within the RVE, and then make a volume average. In this micro-macro-transition procedure, in order to reduce the computational costs the statistically similar representative volume elements (SSRVEs) have been developed to replace real microstructures from metallurgical images by Schröder et al.

Considering the real microstructure of multiphase materials, during the representative volume element generation one should consider grain shape distribution, crystalline orientation distribution, grain boundary misorientation angle distribution and volume fraction of different phases. Figure 1 is an example of the RVE we have generated for TRIP

Recent studies (Lu et al., 2009; 2004) show bulk specimens comprising nanometer sized grains with embedded lamella with coherent, thermal and mechanical stable twin boundaries exhibiting high strength and considerable ductility at the same time. These materials have higher loading rate sensitivity, better tolerance to fatigue crack initiation, and greater resistance to deformation. Under this condition, the RVE with nanometer sized twin lamella inside nanometer sized twin lamella inside nanometer sized grains will help us to understand

Assume two orientations **Q**<sup>I</sup> and **Q**II have the twin relationship in (1, 1, 1) habit plane along [1, 1, −2] twinning direction. For any vector **V**, these two tensors will map as **v**<sup>I</sup> = **Q**I**V** and **v**II = **Q**II**V**. The twin relationship between **v**<sup>I</sup> and **v**II is easier to see in the local twin

system [**x**, **y**, **z**]. We define a orthogonal tensor **R**<sup>L</sup> for the mapping from global coordinate

√ 1 <sup>6</sup> <sup>√</sup> 1 <sup>2</sup> <sup>√</sup> 1 3 ⎤ ⎥ ⎥ ⎦

⎤

**R**L**Q**II**R**<sup>T</sup> L � **v**�

√ 1 <sup>6</sup> √−<sup>1</sup> <sup>2</sup> <sup>√</sup> 1 3

√−2 <sup>6</sup> <sup>0</sup> <sup>√</sup> 1 3

⎡ ⎣

**v**� = **R**<sup>M</sup>

�

⎡ ⎢ ⎢ ⎣ //[1, 1, 1] rather than in the global coordinate

⎦ (2)

. (3)

(1)

//[1, 1, −2] and **z**�

*R*L*ij* =

*R*M*ij* =

**2. Generating realistic material microstructures**

steels where the Voronoi tessellation algorithm has been used.

existing material behavior and design new materials.

and another tensor **R**<sup>M</sup> for the mirror symmetry operation

and get the twin relationship in the local twin coordinate system

**R**L**Q**I**R**<sup>T</sup> L �

�

coordinate system with **x**�

system to the local coordinate system

(2010).

For large strains the elastic and plastic deformation can be separated consistently. We follow the well-known multiplicative decomposition proposed by Lee (1969) (see Figure 2) of the deformation gradient tensor **F**

$$\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \frac{\partial \mathbf{x}}{\partial \widetilde{\mathbf{x}}} \frac{\partial \widetilde{\mathbf{x}}}{\partial \mathbf{X}} = \mathbf{F}\_{\mathbf{e}} \mathbf{F}\_{\mathbf{P}} \tag{7}$$

where **F**e is the elastic part comprising the stretch and rotation of the lattice, and **F**p corresponds to the plastic deformation. The lattice rotation **R**e and stretch **U**e are included in the mapping **F**e. They can be calculated by the polar decomposition **F**<sup>e</sup> = **R**e**U**e, i.e., the texture evolution is included in this part of the deformation. Furthermore, two rate equations can be derived for the elastic and the plastic deformation gradients as

$$
\dot{\mathbf{F}}\_{\mathbf{e}} = \mathbf{L}\mathbf{F}\_{\mathbf{e}} - \mathbf{F}\_{\mathbf{e}}\mathbf{L}\_{\mathbf{p}} \tag{8}
$$

$$\mathbf{F\_P = L\_P F\_P} \tag{9}$$

the Cauchy stress σ amount to

By choosing **F**e0 to satisfy

⎡ ⎣

*F*e0 =

�**l**

the volume change is always purely elastic.

cos(*ϕ*2) − sin(*ϕ*2) 0 sin(*ϕ*2) cos(*ϕ*2) 0 0 01

**3.3 The dislocation slip based plastic deformation**

conventional metal forming processes.

is useful in the later GND calculations.

⎤ ⎦ ⎡ ⎣

*<sup>α</sup>* <sup>=</sup> **<sup>d</sup>**�*<sup>α</sup>* <sup>×</sup> **<sup>n</sup>**�*<sup>α</sup>* we can define one local coordinate system for slip system *<sup>α</sup>* as [

**S** = **F**−<sup>1</sup>

<sup>σ</sup> <sup>=</sup> <sup>1</sup> *J* **<sup>F</sup>**<sup>e</sup> **S F** � *<sup>T</sup>*

where *J* = *det*(**F**) = *det*(**F**e), which means the isochoric plastic deformation is assumed, i.e.,

Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals 7

In a polycrystal the different grains have different initial orientations. Therefore, in the global coordinate system, different stiffness tensors, slip plane normals and slip directions should be specified for every crystal. In order to use only one set of data, a virtual deformation step is introduced before the calculation in the following form: **F**p0 is set as the initial value for **F**p.

the starting value for **F** amounts to **I** as desired. If one adopts the Bunge Euler angle (*ϕ*1, Φ, *ϕ*2) to define the crystal orientation, the matrix of the elastic deformation gradient amounts to

The plastic deformation mechanism discussed here is the slip mechanism where dislocations slip in certain crystallographic planes along certain crystallographic directions to accommodate shape changes of the crystal. This is the most common mechanism in

The concept for describing displacement fields around dislocations in crystals was developed mathematically by Volterra and used for calculating elastic deformation fields by Orowan and Taylor in 1934 (Hirth & Lothe, 1992). Dislocations are one dimensional lattice defects which can not begin or end inside a crystal, but must intersect a free surface, form a closed loop or make junctions with other dislocations. Due to energetic reasons there is a strong tendency for dislocations to assume a minimum Burgers vector, and to slide in the planes with maximum interplanar separation and along the most densely packed directions. Under the applied stress, the lattice deforms elastically, until stretched bonds near the dislocation core break and new bonds form. The dislocation moves step by step by one Burgers vector. It is the dislocation slip mechanism that can explain why the actually observed strength of most crystalline materials is between one to four orders of magnitude smaller than the intrinsic or theoretical strength required for breaking the atom bonds without the presence of dislocations. In order to set the connection between the continuum variables and the process of dislocation slip, we have to determine the shear amount of individual slip systems. The slip systems are mathematically described by the Schmid tensor �**M***<sup>α</sup>* <sup>=</sup> **<sup>d</sup>**�*<sup>α</sup>* <sup>⊗</sup> **<sup>n</sup>**�*<sup>α</sup>* where **<sup>d</sup>**�*<sup>α</sup>* <sup>=</sup> **<sup>b</sup>**/*<sup>b</sup>* expresses the slip direction, which is parallel to the Burgers vector **<sup>b</sup>**, but normalised, and **<sup>n</sup>**�*α*, the slip plane normal with respect to the undistorted configuration. Through calculating the line vector

10 0 0 cos(Φ) − sin(Φ) 0 sin(Φ) cos(Φ)

<sup>p</sup> **S F** � <sup>−</sup>*<sup>T</sup>*

<sup>p</sup> (15)

<sup>e</sup> (16)

cos(*ϕ*1) − sin(*ϕ*1) 0 sin(*ϕ*1) cos(*ϕ*1) 0 0 01

�**l**

*<sup>α</sup>*, **<sup>d</sup>**�*α*, **<sup>n</sup>**�*α*], which

⎤ ⎦ .

**F**e0 **F**p0 = **I**, **F**e0, **F**p0 ∈ *Orth*. (17)

⎤ ⎦ ⎡ ⎣

where **L** = **FF**˙ <sup>−</sup><sup>1</sup> and **L**<sup>p</sup> = **F**˙ <sup>p</sup>**F**−<sup>1</sup> <sup>p</sup> are the total and the plastic velocity gradients defined in the current and the unloaded configuration respectively. Because the stress produced by the elastic deformation can supply driving forces for dislocation slip, twinning formation and phase transformation which can accommodate the plastic deformation, **F**e and **F**p are not independent. If the total deformation process is known, the elastic and plastic deformation evolutions can be determined through solving equations 7,8 and 9.

When the eigen-deformation **F**t of phase transformation and the plastic deformation **F**p of dislocation slip coexist, the multiply decomposition (see Figure 3) should be reformulated as the following

$$\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \widetilde{\mathbf{x}}'} \frac{\partial \widetilde{\mathbf{x}}'}{\partial \widetilde{\mathbf{x}}} \frac{\partial \widetilde{\mathbf{x}}}{\partial \mathbf{X}} = \mathbf{F}\_{\mathbf{e}} \mathbf{F}\_{\mathbf{t}} \mathbf{F}\_{\mathbf{P}}.\tag{10}$$

The evolution of **F**<sup>t</sup> is controlled by the transformed volume fraction *f<sup>α</sup>* because the eigen-deformation **<sup>H</sup>** *<sup>α</sup>* <sup>t</sup> of each transformation system with unit volume fraction is a constant tensor

$$\mathbf{F}\_{\mathbf{t}} = \mathbf{I} + \sum\_{\alpha=1}^{N\_{\mathbf{T}}} f^{\alpha} \widetilde{\mathbf{H}}\_{\mathbf{t}}^{\alpha} \tag{11}$$

$$\mathbf{F}\_{\mathbf{l}} = \sum\_{\alpha=1}^{N\_T} f^{\alpha} \tilde{\mathbf{H}}\_{\mathbf{l}}^{\alpha}. \tag{12}$$

where *NT* is the total number of transformation system.

Fig. 3. The multiplicative decomposition of the deformation gradient when dislocation slip and phase transformation coexist.

#### **3.2 The elastic deformation**

For the dislocation slip case, the plastic deformation **F**p will not change the lattice orientation, i.e., we can use a constant stiffness tensor **K** <sup>0</sup> for the stress calculations and define the elastic law in the unloaded configuration. By defining the second Piola-Kirchhoff stress tensor **S** in the unloaded configuration and its work conjugated elastic Green strain tensor **E**˜, the elastic law is derived as

$$
\bar{\mathbf{S}} = \bar{\mathbb{K}}\_0 \bar{\mathbf{E}} \tag{13}
$$

with

$$\tilde{\mathbf{E}} = \frac{1}{2} \left( \mathbf{F}\_{\mathbf{e}}^{\mathrm{T}} \mathbf{F}\_{\mathbf{e}} - \mathbf{I} \right) \tag{14}$$

where **I** is the second order unity tensor. When the variational principle of the FEM is formulated in the reference or current configuration, the second Piola-Kirchhoff stress **S** or the Cauchy stress σ amount to

4 Will-be-set-by-IN-TECH

the current and the unloaded configuration respectively. Because the stress produced by the elastic deformation can supply driving forces for dislocation slip, twinning formation and phase transformation which can accommodate the plastic deformation, **F**e and **F**p are not independent. If the total deformation process is known, the elastic and plastic deformation

When the eigen-deformation **F**t of phase transformation and the plastic deformation **F**p of dislocation slip coexist, the multiply decomposition (see Figure 3) should be reformulated as

The evolution of **F**<sup>t</sup> is controlled by the transformed volume fraction *f<sup>α</sup>* because the

*NT* ∑ *α*=1

*NT* ∑ *α*=1 ˙ *<sup>f</sup> <sup>α</sup>***<sup>H</sup>** *<sup>α</sup>*

Fig. 3. The multiplicative decomposition of the deformation gradient when dislocation slip

For the dislocation slip case, the plastic deformation **F**p will not change the lattice orientation, i.e., we can use a constant stiffness tensor **K** <sup>0</sup> for the stress calculations and define the elastic law in the unloaded configuration. By defining the second Piola-Kirchhoff stress tensor **S** in the unloaded configuration and its work conjugated elastic Green strain tensor **E**˜, the elastic

where **I** is the second order unity tensor. When the variational principle of the FEM is formulated in the reference or current configuration, the second Piola-Kirchhoff stress **S** or

**<sup>E</sup>** <sup>=</sup> <sup>1</sup> 2 **F**T <sup>e</sup> **F**<sup>e</sup> − **I** 

*<sup>f</sup> <sup>α</sup>***<sup>H</sup>** *<sup>α</sup>*

<sup>t</sup> of each transformation system with unit volume fraction is a constant

*∂***x**� *∂***x** *∂***x**

**F**<sup>t</sup> = **I** +

**F**˙ <sup>t</sup> =

evolutions can be determined through solving equations 7,8 and 9.

where *NT* is the total number of transformation system.

**<sup>F</sup>** <sup>=</sup> *<sup>∂</sup>***<sup>x</sup>** *∂***x**�

<sup>p</sup> are the total and the plastic velocity gradients defined in

*<sup>∂</sup>***<sup>X</sup>** <sup>=</sup> **<sup>F</sup>**e**F**t**F**p. (10)

<sup>t</sup> (11)

<sup>t</sup> . (12)

**S** = **K** <sup>0</sup>**E** (13)

(14)

where **L** = **FF**˙ <sup>−</sup><sup>1</sup> and **L**<sup>p</sup> = **F**˙ <sup>p</sup>**F**−<sup>1</sup>

the following

tensor

eigen-deformation **<sup>H</sup>** *<sup>α</sup>*

and phase transformation coexist.

**3.2 The elastic deformation**

law is derived as

with

$$\mathbf{S} = \mathbf{F}\_{\mathbf{P}}^{-1} \, \widetilde{\mathbf{S}} \, \mathbf{F}\_{\mathbf{P}}^{-T} \tag{15}$$

$$
\sigma = \frac{1}{J} \mathbf{F}\_{\mathbf{e}} \,\, \widetilde{\mathbf{S}} \,\, \mathbf{F}\_{\mathbf{e}}^{T} \tag{16}
$$

where *J* = *det*(**F**) = *det*(**F**e), which means the isochoric plastic deformation is assumed, i.e., the volume change is always purely elastic.

In a polycrystal the different grains have different initial orientations. Therefore, in the global coordinate system, different stiffness tensors, slip plane normals and slip directions should be specified for every crystal. In order to use only one set of data, a virtual deformation step is introduced before the calculation in the following form: **F**p0 is set as the initial value for **F**p. By choosing **F**e0 to satisfy

$$\mathbf{F}\_{\rm e0} \mathbf{F}\_{\rm p0} = \mathbf{I}, \quad \mathbf{F}\_{\rm e0}, \mathbf{F}\_{\rm p0} \in \operatorname{Orth}.\tag{17}$$

the starting value for **F** amounts to **I** as desired. If one adopts the Bunge Euler angle (*ϕ*1, Φ, *ϕ*2) to define the crystal orientation, the matrix of the elastic deformation gradient amounts to

$$F\_{\mathbf{e}0} = \begin{bmatrix} \cos(\varphi\_2) - \sin(\varphi\_2) \ 0 \\ \sin(\varphi\_2) & \cos(\varphi\_2) \ 0 \\ 0 & 0 \ 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 \cos(\Phi) - \sin(\Phi) & \\ 0 \sin(\Phi) & \cos(\Phi) \end{bmatrix} \begin{bmatrix} \cos(\varphi\_1) - \sin(\varphi\_1) \ 0 \\ \sin(\varphi\_1) & \cos(\varphi\_1) \ 0 \\ 0 & 0 \ 1 \end{bmatrix} \dots$$

#### **3.3 The dislocation slip based plastic deformation**

The plastic deformation mechanism discussed here is the slip mechanism where dislocations slip in certain crystallographic planes along certain crystallographic directions to accommodate shape changes of the crystal. This is the most common mechanism in conventional metal forming processes.

The concept for describing displacement fields around dislocations in crystals was developed mathematically by Volterra and used for calculating elastic deformation fields by Orowan and Taylor in 1934 (Hirth & Lothe, 1992). Dislocations are one dimensional lattice defects which can not begin or end inside a crystal, but must intersect a free surface, form a closed loop or make junctions with other dislocations. Due to energetic reasons there is a strong tendency for dislocations to assume a minimum Burgers vector, and to slide in the planes with maximum interplanar separation and along the most densely packed directions. Under the applied stress, the lattice deforms elastically, until stretched bonds near the dislocation core break and new bonds form. The dislocation moves step by step by one Burgers vector. It is the dislocation slip mechanism that can explain why the actually observed strength of most crystalline materials is between one to four orders of magnitude smaller than the intrinsic or theoretical strength required for breaking the atom bonds without the presence of dislocations.

In order to set the connection between the continuum variables and the process of dislocation slip, we have to determine the shear amount of individual slip systems. The slip systems are mathematically described by the Schmid tensor �**M***<sup>α</sup>* <sup>=</sup> **<sup>d</sup>**�*<sup>α</sup>* <sup>⊗</sup> **<sup>n</sup>**�*<sup>α</sup>* where **<sup>d</sup>**�*<sup>α</sup>* <sup>=</sup> **<sup>b</sup>**/*<sup>b</sup>* expresses the slip direction, which is parallel to the Burgers vector **<sup>b</sup>**, but normalised, and **<sup>n</sup>**�*α*, the slip plane normal with respect to the undistorted configuration. Through calculating the line vector �**l** *<sup>α</sup>* <sup>=</sup> **<sup>d</sup>**�*<sup>α</sup>* <sup>×</sup> **<sup>n</sup>**�*<sup>α</sup>* we can define one local coordinate system for slip system *<sup>α</sup>* as [ �**l** *<sup>α</sup>*, **<sup>d</sup>**�*α*, **<sup>n</sup>**�*α*], which is useful in the later GND calculations.

pile-up size, *L*, and of the temperature, *θ*; i.e.,

case of infinitesimally small elastic stretches

(1992)

and

cutting event

and

*v* = *v* 

plane. Both *ρ*SSD and *ρ*GND are contributing to *ρ*<sup>F</sup> and *ρ*<sup>P</sup>

*ρα* <sup>F</sup> =

*ρα* <sup>P</sup> =

complicated than the schematic pictures we use here.

*NS* ∑ *β*=1

*NS* ∑ *β*=1 *<sup>χ</sup>αβ ρ β* SSD + *ρ β* GND 

*<sup>χ</sup>αβ ρ β* SSD + *ρ β* GND sin(**<sup>n</sup>***<sup>α</sup>*,**<sup>t</sup>**

FCC crystal, and use a single set of interaction strengths for both SSDs and GNDs.

where we introduce the interaction strength, *χαβ*, between different slip systems, which includes the self interaction strength, coplanar interaction strength, cross slip strength, glissile junction strength, Hirth lock strength, and Lomer-Cottrell lock strength. One can go further to see the definition of these interactions in literature (Devincre et al., 2008; Madec et al., 2008). In this formulation we only consider edge dislocations owing to their limited mobility for the

With the help of the forest dislocation density *ρ*F, we can determine the average jump distance of the mobile dislocation and the activation volume for the thermal activated forest dislocation

> *<sup>λ</sup>* <sup>=</sup> *<sup>c</sup>*<sup>1</sup> <sup>√</sup>*ρ<sup>F</sup>*

where *c*<sup>1</sup> and *c*<sup>2</sup> are constants to reflect the real dislocation line configuration which is more

With the help of the parallel dislocation density *<sup>ρ</sup>*<sup>P</sup> and the gradient of GND density *∂ρ*GND

we can determine the average athermal passing stress *τ*<sup>p</sup> and back stress *τ*<sup>b</sup> as following

*<sup>τ</sup>*, *<sup>ρ</sup>*SSD, *<sup>ρ</sup>*GND, *∂ρ*GND

Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals 9

The resolved shear stress, *τ*, is the projection of the stress measure onto the slip system. In the

*<sup>τ</sup>*, within the intermediate configuration **<sup>x</sup>** can be approximated by following Kalidindi et al.

In order to accommodate a part of the external plastic deformation, the mobile dislocations, *ρ*M, must overcome the stress field of the parallel dislocations, *ρ*P, which cause the passing stress. They must also cut the forest dislocations, *ρ*F, with the aid of thermal activation. We define the parallel dislocation density and the forest dislocation density as: *ρ*<sup>P</sup> for all dislocations parallel to the slip plane, and *ρ*<sup>F</sup> for the dislocations perpendicular to the slip

 **C** <sup>=</sup> **F**T <sup>e</sup> **F**<sup>e</sup> 

*<sup>∂</sup><sup>X</sup>* , *<sup>L</sup>*, *<sup>θ</sup>*

*<sup>τ</sup><sup>α</sup>* <sup>=</sup> **<sup>S</sup> <sup>C</sup>** · **M***<sup>α</sup>* <sup>∼</sup><sup>=</sup> **<sup>S</sup>**˜ · **M***<sup>α</sup>* (22)

cos(**<sup>n</sup>***<sup>α</sup>*,**<sup>t</sup>**

*β*) 

*β*) 

*V* = *c*2*b*2*λ* (26)

*<sup>τ</sup>*<sup>p</sup> <sup>=</sup> *<sup>c</sup>*3*Gb*√*ρ<sup>P</sup>* (27)

<< 1, the resolved shear stress,

(23)

(24)

(25)

*<sup>∂</sup><sup>X</sup>* ,

(21)

For the FCC crystal structure, the close-packed planes {111} and close-packed directions �110� form 12 slip systems. For the BCC crystal structure, the pencil glide phenomenon is observed, which resembles slip in a fixed direction on apparently random planes. In literature, experimental studies have shown that for BCC crystals the slip direction is along �111�, and three groups of slip planes exist, including {110},{112} and the less common {123} planes. Totally there are 48 slip systems for BCC crystals. Therefore, for FCC and BCC crystals it is possible to supply five independent slip systems to accommodate any arbitrary external plastic deformation, and in the middle temperature range the slip is the main mechanism for plastic deformations. For the HCP crystal, the slip system number is dependent on the axis ratio of the HCP unit cell. When this ratio assumes values in a certain range as discussed by Gottstein (2004), only two independent slip systems exist, and it is impossible to accommodate a arbitrary deformation by slip steps. As a result, mechanical twins appear during plastic deformation.

Among all of the dislocations in one slip system, only the mobile dislocations can produce plastic deformation, and their speed can be determined by the forces acting on them. In general, the driving force is related to external loads, short range isotropic resistance of dislocation interactions and long range back stress of dislocation pile-ups and lattice frictions. The widely-adopted constitutive assumption for crystal plasticity reads

$$\mathbf{L}\_{\mathbf{P}} = \sum\_{a=1}^{N\_{\mathbf{S}}} \boldsymbol{\gamma}^{a} \widetilde{\mathbf{M}}^{a} \,, \qquad \dot{\mathbf{F}}\_{\mathbf{P}} = \sum\_{a=1}^{N\_{\mathbf{S}}} \boldsymbol{\gamma}^{a} \widetilde{\mathbf{M}}^{a} \mathbf{F}\_{\mathbf{P}} \tag{18}$$

where *<sup>γ</sup>*˙ *<sup>α</sup>* is the slip rate on slip system *<sup>α</sup>* within the intermediate configuration **<sup>x</sup>**, and *NS*, the total number of slip systems.

When a part of material is transferred into another lattice with volume fraction *f* = ∑*NT <sup>α</sup>*=<sup>1</sup> *<sup>f</sup> <sup>α</sup>* which can only deform elastically, equation 18 has been modified as

$$\mathbf{L}\_{\mathbf{P}} = (1 - f) \sum\_{a=1}^{N\_{\mathbf{S}}} \dot{\gamma}^a \widetilde{\mathbf{M}}^a \,, \qquad \dot{\mathbf{F}}\_{\mathbf{P}} = (1 - f) \sum\_{a=1}^{N\_{\mathbf{S}}} \dot{\gamma}^a \widetilde{\mathbf{M}}^a \mathbf{F}\_{\mathbf{P}}.\tag{19}$$

#### **3.3.1 The Orowan equation**

Commonly used expressions for the relation of the shear rate, *γ*˙ , and the resolved shear stress, *τ*, include a phenomenological viscoplastic law in the form of a power law by Peirce et al. (1982), and more physically-based ones such as those of Kocks et al. (1975) and Nemat-Nasser et al. (1998), which can take account of rate and temperature dependencies. In this paper we use the Orowan equation to calculate the plastic shear rate *γ*˙ of each slip system as a function of the mobile dislocation density, *ρ*m, on that slip system

$$
\dot{\gamma} = \rho\_{\rm m} \,\, b \,\, v \tag{20}
$$

where the average velocity of the mobile dislocations, *v*, is a function of the resolved shear stress, *τ*, of the dislocation densities, *ρ*M, *ρ*SSD and *ρ*GND and its gradient, the average GND pile-up size, *L*, and of the temperature, *θ*; i.e.,

$$v = v\left(\tau, \rho\_{\text{SSD}\prime} \rho\_{\text{GND}\prime} \frac{\partial \rho\_{\text{GND}}}{\partial X}, L, \theta\right) \tag{21}$$

The resolved shear stress, *τ*, is the projection of the stress measure onto the slip system. In the case of infinitesimally small elastic stretches **C** <sup>=</sup> **F**T <sup>e</sup> **F**<sup>e</sup> << 1, the resolved shear stress, *<sup>τ</sup>*, within the intermediate configuration **<sup>x</sup>** can be approximated by following Kalidindi et al. (1992)

$$
\tau^{\alpha} = \mathbf{S} \mathbf{C} \cdot \dot{\mathbf{M}}^{\alpha} \cong \bar{\mathbf{S}} \cdot \bar{\mathbf{M}}^{\alpha} \tag{22}
$$

In order to accommodate a part of the external plastic deformation, the mobile dislocations, *ρ*M, must overcome the stress field of the parallel dislocations, *ρ*P, which cause the passing stress. They must also cut the forest dislocations, *ρ*F, with the aid of thermal activation. We define the parallel dislocation density and the forest dislocation density as: *ρ*<sup>P</sup> for all dislocations parallel to the slip plane, and *ρ*<sup>F</sup> for the dislocations perpendicular to the slip plane. Both *ρ*SSD and *ρ*GND are contributing to *ρ*<sup>F</sup> and *ρ*<sup>P</sup>

$$\rho\_{\mathbf{F}}^{a} = \sum\_{\beta=1}^{N\_{\mathbf{S}}} \chi^{a\beta} \left( \rho\_{\mathbf{SSD}}^{\beta} + \rho\_{\mathbf{GND}}^{\beta} \right) \left| \cos(\tilde{\mathbf{n}}^{a}, \tilde{\mathbf{f}}^{\beta}) \right| \tag{23}$$

and

6 Will-be-set-by-IN-TECH

For the FCC crystal structure, the close-packed planes {111} and close-packed directions �110� form 12 slip systems. For the BCC crystal structure, the pencil glide phenomenon is observed, which resembles slip in a fixed direction on apparently random planes. In literature, experimental studies have shown that for BCC crystals the slip direction is along �111�, and three groups of slip planes exist, including {110},{112} and the less common {123} planes. Totally there are 48 slip systems for BCC crystals. Therefore, for FCC and BCC crystals it is possible to supply five independent slip systems to accommodate any arbitrary external plastic deformation, and in the middle temperature range the slip is the main mechanism for plastic deformations. For the HCP crystal, the slip system number is dependent on the axis ratio of the HCP unit cell. When this ratio assumes values in a certain range as discussed by Gottstein (2004), only two independent slip systems exist, and it is impossible to accommodate a arbitrary deformation by slip steps. As a result, mechanical twins appear during plastic

Among all of the dislocations in one slip system, only the mobile dislocations can produce plastic deformation, and their speed can be determined by the forces acting on them. In general, the driving force is related to external loads, short range isotropic resistance of dislocation interactions and long range back stress of dislocation pile-ups and lattice frictions.

*<sup>γ</sup>*˙ *<sup>α</sup>***<sup>M</sup>***<sup>α</sup>* , **<sup>F</sup>**˙ <sup>p</sup> <sup>=</sup>

When a part of material is transferred into another lattice with volume fraction *f* = ∑*NT*

where *<sup>γ</sup>*˙ *<sup>α</sup>* is the slip rate on slip system *<sup>α</sup>* within the intermediate configuration **<sup>x</sup>**, and *NS*, the

*<sup>γ</sup>*˙ *<sup>α</sup>***<sup>M</sup>***<sup>α</sup>* , **<sup>F</sup>**˙ <sup>p</sup> = (<sup>1</sup> <sup>−</sup> *<sup>f</sup>*)

Commonly used expressions for the relation of the shear rate, *γ*˙ , and the resolved shear stress, *τ*, include a phenomenological viscoplastic law in the form of a power law by Peirce et al. (1982), and more physically-based ones such as those of Kocks et al. (1975) and Nemat-Nasser et al. (1998), which can take account of rate and temperature dependencies. In this paper we use the Orowan equation to calculate the plastic shear rate *γ*˙ of each slip system as a function

where the average velocity of the mobile dislocations, *v*, is a function of the resolved shear stress, *τ*, of the dislocation densities, *ρ*M, *ρ*SSD and *ρ*GND and its gradient, the average GND

*NS* ∑ *α*=1

> *NS* ∑ *α*=1

*γ*˙ = *ρ*<sup>m</sup> *b v* (20)

*<sup>γ</sup>*˙ *<sup>α</sup>***<sup>M</sup>***α***F**<sup>p</sup> (18)

*<sup>γ</sup>*˙ *<sup>α</sup>***<sup>M</sup>***α***F**p. (19)

*<sup>α</sup>*=<sup>1</sup> *<sup>f</sup> <sup>α</sup>*

The widely-adopted constitutive assumption for crystal plasticity reads

*NS* ∑ *α*=1

which can only deform elastically, equation 18 has been modified as

*NS* ∑ *α*=1

**L**<sup>p</sup> =

**L**<sup>p</sup> = (1 − *f*)

of the mobile dislocation density, *ρ*m, on that slip system

deformation.

total number of slip systems.

**3.3.1 The Orowan equation**

$$\rho^a\_\mathbf{P} = \sum\_{\beta=1}^{N\_\mathbf{S}} \chi^{a\beta} \left( \rho^\beta\_{\mathbf{SSD}} + \rho^\beta\_{\mathbf{GND}} \right) \left| \sin(\tilde{\mathbf{n}}^a, \tilde{\mathbf{t}}^\beta) \right| \tag{24}$$

where we introduce the interaction strength, *χαβ*, between different slip systems, which includes the self interaction strength, coplanar interaction strength, cross slip strength, glissile junction strength, Hirth lock strength, and Lomer-Cottrell lock strength. One can go further to see the definition of these interactions in literature (Devincre et al., 2008; Madec et al., 2008). In this formulation we only consider edge dislocations owing to their limited mobility for the FCC crystal, and use a single set of interaction strengths for both SSDs and GNDs.

With the help of the forest dislocation density *ρ*F, we can determine the average jump distance of the mobile dislocation and the activation volume for the thermal activated forest dislocation cutting event

$$
\lambda = \frac{c\_1}{\sqrt{\rho\_F}} \tag{25}
$$

and

$$V = c\_2 b^2 \lambda \tag{26}$$

where *c*<sup>1</sup> and *c*<sup>2</sup> are constants to reflect the real dislocation line configuration which is more complicated than the schematic pictures we use here.

With the help of the parallel dislocation density *<sup>ρ</sup>*<sup>P</sup> and the gradient of GND density *∂ρ*GND *<sup>∂</sup><sup>X</sup>* , we can determine the average athermal passing stress *τ*<sup>p</sup> and back stress *τ*<sup>b</sup> as following

$$
\pi\_{\rm P} = c\_3 Gb\sqrt{\rho\_P} \tag{27}
$$

calculate the net Burgers vector for an area

calculation

and

**<sup>Λ</sup>** <sup>=</sup> **<sup>b</sup>**¯ <sup>⊗</sup> ¯**<sup>l</sup>** <sup>=</sup> <sup>−</sup>

**b***α* GNDe =

**b***α* GNDs =

Furthermore we also can calculate the GND density as the following

GND <sup>=</sup> (�**b***<sup>α</sup>*

Fig. 4. A transformation system of the austenite-martensite phase transformation.

The transformation-induced plasticity (TRIP) assisted steels are mixtures of allotriomorphic ferrite, bainite and retained austenite. Experimental and modelling publications have highlighted that the transformation of retained austenite to martensite under the influence of a applied stress or strain can improve material ductility and strength efficiently, as shown

According to the geometrically nonlinear theory of martensitic transformations (Bhattacharya, 1993; Hane & Shield, 1998; 1999) there are 24 transformation systems and they are constructed by two body-centered tetragonal (BCT) variants with relative rotations and volume fractions, in order to produce habit planes between austenite and martensite arrays and pairwise arranging twin related variant lamellas. Each transformation system corresponds to one

*ρα*

**3.4 Eigenstrain of phase transformations**

by Bhadeshia (2002).

where ∇<sup>X</sup> = *∂*/*∂***X**, is defined as the derivative with respect to the reference coordinates and **b**¯ and ¯**l** are, respectively, the net Burgers vector and net line vector after an volume average operation. Using equation (31) the resulting Burgers vector for a circuit with an arbitrary orientation can be calculated. In general this tensor is non-symmetric and it can be mapped to nine independent slip systems in a unique fashion. For the FCC crystal structure with its 12 slip systems, only 5 systems are independent according to the von Mises-Taylor constraint. This implies that it is impossible to calculate the exact amount of GNDs for every slip system in a unique way. Nevertheless, we can project **Λ** to each of the slip systems to determine the Burgers vector of the edge and screw type GNDs for the pass stress and backing stress

Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals 11

**<sup>d</sup>***<sup>α</sup>* · **<sup>Λ</sup>** ·**<sup>l</sup>** *α* 

**<sup>d</sup>***<sup>α</sup>* · **<sup>Λ</sup>** · **<sup>d</sup>***<sup>α</sup>*

GNDe� <sup>+</sup> <sup>|</sup>**b***<sup>α</sup>*

∇<sup>X</sup> × **F**<sup>p</sup>

<sup>T</sup> (31)

**<sup>d</sup>***<sup>α</sup>* (32)

**<sup>d</sup>***<sup>α</sup>*. (33)

GNDs|) /*b*. (34)

and

$$
\pi\_{\mathsf{b}} = GbL^2 \frac{\partial \rho\_{\mathsf{GND}}}{\partial X} \tag{28}
$$

where *c*<sup>3</sup> is the constant for the Taylor hardening mechanism. For reasons of simplicity, in equation 28 the back stress of one slip system only comes from the GND pile-up of this slip system. This equation can be easily extended to consider the back stress from all of the slip system at the same time.

Compared with flow rules in the literature which contain a constant reference shear rate and a constant rate sensitivity exponent, here a flow rule is derived based on the dislocation slip mechanism

$$\dot{\gamma} = \begin{cases} \rho\_{\rm m} b \lambda \nu\_{0} \exp\left(-\frac{Q\_{\rm slip}}{\rm kg\theta}\right) \exp\left(\frac{|\tau + \tau\_{\rm b}| - \tau\_{\rm p}}{\rm kg\theta} V\right) \text{sign}\left(\tau + \tau\_{\rm b}\right) & |\tau + \tau\_{\rm b}| > \tau\_{\rm p} \\\\ 0 & |\tau + \tau\_{\rm b}| \le \tau\_{\rm p} \end{cases} \tag{29}$$

where kB is the Boltzmann constant, *ν*<sup>0</sup> the attempt frequency and *Qslip* the effective activation energy.

Inside the flow rule given by equation 29 determination of the mobile dislocation density is a hard task. In some research work, the mobile dislocation density was found to be a small fraction of total dislocation density and is even treated as a constant. The more sophisticated model to deal with this dislocation density based on energy minimization can be found in Ma & Roters (2004). For reasons of simplicity here the mobile dislocation density is treated as a constant number.

#### **3.3.2 Evolution of the dislocation densities**

There are four processes contributing to the evolution of the SSD density as discussed by Ma (2006). The lock forming mechanism between mobile dislocations and forest dislocations, the dipole forming mechanism between mobile dislocations with parallel line vectors, and anti-parallel Burgers vector determine the multiplication terms, while the annihilation term includes annihilation between one mobile dislocation with another immobile one and annihilation between two immobile dislocations. The often used Kocks-Mecking model, as discussed in Roters (1999), only adopts the locks formation and mobile-immobile annihilation mechanisms for the SSD evolution

$$\phi\_{\rm SSD} = \left(c\_4 \sqrt{\rho\_{\rm F}} - c\_5 \,\rho\_{\rm SSD}\right) \dot{\gamma} \tag{30}$$

Where *c*<sup>4</sup> and *c*<sup>5</sup> are constants used to adjust the locks and annihilation radius.

When plastic deformation gradients are present in a volume portion, GNDs must be introduced to preserve the continuity of the crystal lattice. A relation between a possible GND measure and the plastic deformation gradient has been proposed by Nye (1953). This approach has been later extended to a more physically motivated continuum approach to generally account for strain gradient effects by Dai & Parks (1997). Following these pioneering approaches, we use as a dislocation density tensor, **Λ**, for a selected volume portion to calculate the net Burgers vector for an area

$$\mathbf{A} = \bar{\mathbf{b}} \otimes \bar{\mathbf{l}} = -\left(\nabla\_{\mathbf{X}} \times \mathbf{F}\_{\mathbf{p}}\right)^{\mathrm{T}} \tag{31}$$

where ∇<sup>X</sup> = *∂*/*∂***X**, is defined as the derivative with respect to the reference coordinates and **b**¯ and ¯**l** are, respectively, the net Burgers vector and net line vector after an volume average operation. Using equation (31) the resulting Burgers vector for a circuit with an arbitrary orientation can be calculated. In general this tensor is non-symmetric and it can be mapped to nine independent slip systems in a unique fashion. For the FCC crystal structure with its 12 slip systems, only 5 systems are independent according to the von Mises-Taylor constraint. This implies that it is impossible to calculate the exact amount of GNDs for every slip system in a unique way. Nevertheless, we can project **Λ** to each of the slip systems to determine the Burgers vector of the edge and screw type GNDs for the pass stress and backing stress calculation

$$\mathbf{b}\_{\rm GNNDe}^{\boldsymbol{\alpha}} = \left(\boldsymbol{\tilde{\mathbf{d}}}^{\boldsymbol{\alpha}} \cdot \boldsymbol{\Lambda} \cdot \boldsymbol{\tilde{\mathbf{l}}}^{\boldsymbol{\alpha}}\right) \boldsymbol{\tilde{\mathbf{d}}}^{\boldsymbol{\alpha}} \tag{32}$$

and

8 Will-be-set-by-IN-TECH

*<sup>τ</sup>*<sup>b</sup> <sup>=</sup> *GbL*<sup>2</sup> *∂ρ*GND

where *c*<sup>3</sup> is the constant for the Taylor hardening mechanism. For reasons of simplicity, in equation 28 the back stress of one slip system only comes from the GND pile-up of this slip system. This equation can be easily extended to consider the back stress from all of the slip

Compared with flow rules in the literature which contain a constant reference shear rate and a constant rate sensitivity exponent, here a flow rule is derived based on the dislocation slip

> � <sup>|</sup>*τ*+*τ*b|−*τ*<sup>p</sup> kB*<sup>θ</sup> V*

where kB is the Boltzmann constant, *ν*<sup>0</sup> the attempt frequency and *Qslip* the effective activation

Inside the flow rule given by equation 29 determination of the mobile dislocation density is a hard task. In some research work, the mobile dislocation density was found to be a small fraction of total dislocation density and is even treated as a constant. The more sophisticated model to deal with this dislocation density based on energy minimization can be found in Ma & Roters (2004). For reasons of simplicity here the mobile dislocation density is treated as a

There are four processes contributing to the evolution of the SSD density as discussed by Ma (2006). The lock forming mechanism between mobile dislocations and forest dislocations, the dipole forming mechanism between mobile dislocations with parallel line vectors, and anti-parallel Burgers vector determine the multiplication terms, while the annihilation term includes annihilation between one mobile dislocation with another immobile one and annihilation between two immobile dislocations. The often used Kocks-Mecking model, as discussed in Roters (1999), only adopts the locks formation and mobile-immobile annihilation

When plastic deformation gradients are present in a volume portion, GNDs must be introduced to preserve the continuity of the crystal lattice. A relation between a possible GND measure and the plastic deformation gradient has been proposed by Nye (1953). This approach has been later extended to a more physically motivated continuum approach to generally account for strain gradient effects by Dai & Parks (1997). Following these pioneering approaches, we use as a dislocation density tensor, **Λ**, for a selected volume portion to

*ρ*˙SSD = (*c*<sup>4</sup>

Where *c*<sup>4</sup> and *c*<sup>5</sup> are constants used to adjust the locks and annihilation radius.

�

*<sup>∂</sup><sup>X</sup>* (28)

sign (*τ* + *τ*b) |*τ* + *τ*b| > *τ*<sup>p</sup>

<sup>√</sup>*ρ*<sup>F</sup> <sup>−</sup> *<sup>c</sup>*<sup>5</sup> *<sup>ρ</sup>*SSD) *<sup>γ</sup>*˙ (30)

(29)

0 |*τ* + *τ*b| ≤ *τ*<sup>p</sup>

and

system at the same time.

⎧ ⎪⎪⎨

*ρ*m*bλν*<sup>0</sup> exp

**3.3.2 Evolution of the dislocation densities**

mechanisms for the SSD evolution

� <sup>−</sup> *<sup>Q</sup>*slip kB*θ* � exp

⎪⎪⎩

*γ*˙ =

constant number.

mechanism

energy.

$$\mathbf{b}^{a}\_{\text{GNDs}} = \left(\widetilde{\mathbf{d}}^{a} \cdot \mathbf{A} \cdot \widetilde{\mathbf{d}}^{a}\right) \widetilde{\mathbf{d}}^{a}.\tag{33}$$

Furthermore we also can calculate the GND density as the following

$$
\rho\_{\rm GND}^{\mathfrak{a}} = \left( ||\mathbf{b}\_{\rm GNDe}^{\mathfrak{a}}|| + |\mathbf{b}\_{\rm GNDs}^{\mathfrak{a}}| \right) / b. \tag{34}
$$

Fig. 4. A transformation system of the austenite-martensite phase transformation.

#### **3.4 Eigenstrain of phase transformations**

The transformation-induced plasticity (TRIP) assisted steels are mixtures of allotriomorphic ferrite, bainite and retained austenite. Experimental and modelling publications have highlighted that the transformation of retained austenite to martensite under the influence of a applied stress or strain can improve material ductility and strength efficiently, as shown by Bhadeshia (2002).

According to the geometrically nonlinear theory of martensitic transformations (Bhattacharya, 1993; Hane & Shield, 1998; 1999) there are 24 transformation systems and they are constructed by two body-centered tetragonal (BCT) variants with relative rotations and volume fractions, in order to produce habit planes between austenite and martensite arrays and pairwise arranging twin related variant lamellas. Each transformation system corresponds to one constant shape strain vector, **<sup>v</sup>**�*<sup>α</sup> <sup>s</sup>* , and one constant habit plane normal vector, �**v***<sup>α</sup> <sup>n</sup>*, see Figure 4. Following the classical Kurdjumov-Kaminsky relations, these two vectors are influenced by the carbon concentration through the lattice parameter magnitude variation (Hane & Shield, 1998; Wechsler et al., 1953). The eigenstrain of the transformation system *α* amounts to

$$
\tilde{\mathbf{H}}\_{\mathbf{t}}^{\alpha} = \tilde{\mathbf{v}}\_{s}^{\alpha} \otimes \tilde{\mathbf{v}}\_{n}^{\alpha} \tag{35}
$$

Fig. 5. Schematic drawings of cohesive behaviour of grain boundaries along normal (left)

Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals 13

In experimental works and atomistic simulations with respect to deformation of nanocrystalline materials, dislocation glide inside the grains and grain boundary sliding have both been reported. It is obvious that grain-boundary sliding and separation mechanisms begin to play important roles in the overall inelastic response of a polycrystalline material when the grain-size decreases and dislocation activity within the grain interior becomes more difficult. In recent work, the atomic bonds across grain boundaries have been characterized with ab initio calculations within the framework of the density functional theory (Janisch et al., 2010). In this work not only the energetics of grain boundaries have been characterized, but also the mechanical response of a grain boundary to applied loads is studied. Such information can be used to parameterize cohesive zone models based on ab initio calculations. The cohesive zone model is useful for RVE models of polycrystals, in situations when grain boundary deformation needs to be taken into account explicitely, e.g. when grain boundary sliding or damage initiation at grain boundaries or phase bounaries has to be considered. By adjusting the cohesive zone parameters for grain boundary sliding and opening the competing mechanisms of bulk material deformation and grain boundary accommodated deformation can be studied. Furthermore, it is also possible to investigate damage nucleation

We follow Wei & Anand (2004) to generate a rate independent cohesive zone (CZ) modelling approach for the reasons of simplicity. The velocity jump across a cohesive surface has been

The elastic relative velocities are connected with its power-conjugate traction rate by the

For some special grain boundaries there may exist glide anisotropy inside the grain boundary plane, although our knowledge about this topic is far away from formulating this anisotropy for general grain boundaries. So, we have to assume isotropic plastic deformation property inside the grain boundary, and the trace vector, displacement vector and resistance vector are

**u**˙ − **u**˙ <sup>p</sup>

**t**˙ = **Ku**˙ <sup>e</sup> = **K**

**u**˙ = **u**˙ <sup>e</sup> + **u**˙ p. (38)

(39)

additively decomposed into a elastic and a plastic part as follow

and tangential (right) directions.

at GB triple junctions.

interface elastic stiffness tensor

**3.5 Cohesive zone model for grain and phase boundaries**

As an example in the literature (Kouznetsova & Geers, 2008; Tjahjanto, 2008), the shape strain vector and the habit plane normal vector with respect to specific carbon concentrations have been determined and listed. Because the shape strain vector and the habit plane normal vector are explicit functions of austensite and martensite lattice parameters, the components of the two vectors are irrational. **<sup>H</sup>**� *<sup>α</sup>* <sup>t</sup> is similar to the Schmid tensor �**M***<sup>α</sup>* of dislocation slip except that there is volume change *det*(**H**� *<sup>α</sup>* <sup>t</sup> ) > 1.

In literature, the austenite-martensite phase transformation has been formulated as stress and strain driving mechanisms among different temperature regions. The stress controlled transformation often occurs at lower temperatures where the chemical driving force is so high that a external load below austenite yield stress can help the already existing martensite nuclei to grow. At strain controlled transformation regions at higher temperatures, the chemical driving force is so low that additional loads which are higher than the yield stress are needed in order for existing nuclei to continue grow. Due to the fact that plastic deformation in austenite is easier than spontaneous martensite formation, the transformation has to be continued by new nuclei formation at the shear band intersection region according to Olson & Cohen (1972).

Following the Olson-Cohen model, the transformation kinetics are formulated in meso-scale on the micro-band level. At first the shear band density is estimated, then the intersection frequency of shear bands is calculated and lastly the nucleation producing probability is evaluated. This governing equation for martensite volume fraction reads

$$\dot{f}^{a} = \begin{cases} c\_{6} (1 - \sum\_{\tilde{\mathcal{P}}=1}^{N\_{\rm T}} f\_{\tilde{\mathcal{P}}}) \left| \tau\_{\mathbf{a}} \right|^{c\_{7}} \sum\_{i,j=1}^{N\_{\rm S}} \left( 1 - \left| \tilde{\mathbf{1}}\_{i} \cdot \tilde{\mathbf{1}}\_{j} \right| \right) \left( 1 - \left| \tilde{\mathbf{n}}\_{i} \cdot \tilde{\mathbf{n}}\_{j} \right| \right) \sqrt{\left| \dot{\gamma}\_{i} \dot{\gamma}\_{j} \right|} & \tau\_{\mathbf{a}} > 0 \\\\ 0 & \tau\_{\mathbf{a}} \le 0 \end{cases} \tag{36}$$

where *N*<sup>T</sup> and *N*<sup>S</sup> are the total number of transformation systems and slip systems with, respectively, *c*<sup>6</sup> and *c*<sup>7</sup> as two fitting parameters to control the transformation kinetics; and

$$
\tau^{\alpha} \cong \tilde{\mathbf{S}} \cdot \tilde{\mathbf{H}}\_{\text{t}}^{\alpha} \tag{37}
$$

the resolved stress in the transformation system *α* with includes the shear part and the tensile or compression part at the same time. Indeed, in equation 36 the phase transformation is controlled by the external load potential minimisation.

Because the transformed martensitic phase includes twinned wedge microstructures, the Frank-Read dislocation source may suffer higher resistance compared with the original austenitic phase. The dislocation slip based plasticity of martensite has been neglected here.

Fig. 5. Schematic drawings of cohesive behaviour of grain boundaries along normal (left) and tangential (right) directions.

#### **3.5 Cohesive zone model for grain and phase boundaries**

10 Will-be-set-by-IN-TECH

Following the classical Kurdjumov-Kaminsky relations, these two vectors are influenced by the carbon concentration through the lattice parameter magnitude variation (Hane & Shield, 1998; Wechsler et al., 1953). The eigenstrain of the transformation system *α* amounts to

As an example in the literature (Kouznetsova & Geers, 2008; Tjahjanto, 2008), the shape strain vector and the habit plane normal vector with respect to specific carbon concentrations have been determined and listed. Because the shape strain vector and the habit plane normal vector are explicit functions of austensite and martensite lattice parameters, the components of the

In literature, the austenite-martensite phase transformation has been formulated as stress and strain driving mechanisms among different temperature regions. The stress controlled transformation often occurs at lower temperatures where the chemical driving force is so high that a external load below austenite yield stress can help the already existing martensite nuclei to grow. At strain controlled transformation regions at higher temperatures, the chemical driving force is so low that additional loads which are higher than the yield stress are needed in order for existing nuclei to continue grow. Due to the fact that plastic deformation in austenite is easier than spontaneous martensite formation, the transformation has to be continued by new nuclei formation at the shear band intersection region according to Olson

Following the Olson-Cohen model, the transformation kinetics are formulated in meso-scale on the micro-band level. At first the shear band density is estimated, then the intersection frequency of shear bands is calculated and lastly the nucleation producing probability is

where *N*<sup>T</sup> and *N*<sup>S</sup> are the total number of transformation systems and slip systems with, respectively, *c*<sup>6</sup> and *c*<sup>7</sup> as two fitting parameters to control the transformation kinetics; and

*<sup>τ</sup><sup>α</sup>* <sup>∼</sup><sup>=</sup> **<sup>S</sup>**˜ · **<sup>H</sup>**� *<sup>α</sup>*

the resolved stress in the transformation system *α* with includes the shear part and the tensile or compression part at the same time. Indeed, in equation 36 the phase transformation is

Because the transformed martensitic phase includes twinned wedge microstructures, the Frank-Read dislocation source may suffer higher resistance compared with the original austenitic phase. The dislocation slip based plasticity of martensite has been neglected here.

evaluated. This governing equation for martensite volume fraction reads

*<sup>c</sup>*<sup>7</sup> ∑*N*<sup>S</sup> *i*,*j*=1 � 1 − � � � �**l***<sup>i</sup>* ·�**l***<sup>j</sup>* � � � � � 1 − � � � **<sup>n</sup>**�*<sup>i</sup>* · **<sup>n</sup>**�*<sup>j</sup>* � � � � �� � �*γ*˙ *<sup>i</sup>γ*˙ *<sup>j</sup>* � � �

*<sup>s</sup>* <sup>⊗</sup> **<sup>v</sup>**�*<sup>α</sup>*

**H**� *α* <sup>t</sup> <sup>=</sup> **<sup>v</sup>**�*<sup>α</sup>*

<sup>t</sup> ) > 1.

*<sup>s</sup>* , and one constant habit plane normal vector, �**v***<sup>α</sup>*

<sup>t</sup> is similar to the Schmid tensor �**M***<sup>α</sup>* of dislocation slip except that

0 *τα* ≤ 0

<sup>t</sup> (37)

*<sup>n</sup>*, see Figure 4.

*τα* > 0

(36)

*<sup>n</sup>* (35)

constant shape strain vector, **<sup>v</sup>**�*<sup>α</sup>*

two vectors are irrational. **<sup>H</sup>**� *<sup>α</sup>*

& Cohen (1972).

⎧ ⎪⎪⎨

*<sup>c</sup>*6(<sup>1</sup> <sup>−</sup> <sup>∑</sup>*N*<sup>T</sup>

*<sup>β</sup>*=<sup>1</sup> *fβ*)|*τα*|

controlled by the external load potential minimisation.

⎪⎪⎩

˙ *f <sup>α</sup>* =

there is volume change *det*(**H**� *<sup>α</sup>*

In experimental works and atomistic simulations with respect to deformation of nanocrystalline materials, dislocation glide inside the grains and grain boundary sliding have both been reported. It is obvious that grain-boundary sliding and separation mechanisms begin to play important roles in the overall inelastic response of a polycrystalline material when the grain-size decreases and dislocation activity within the grain interior becomes more difficult. In recent work, the atomic bonds across grain boundaries have been characterized with ab initio calculations within the framework of the density functional theory (Janisch et al., 2010). In this work not only the energetics of grain boundaries have been characterized, but also the mechanical response of a grain boundary to applied loads is studied. Such information can be used to parameterize cohesive zone models based on ab initio calculations.

The cohesive zone model is useful for RVE models of polycrystals, in situations when grain boundary deformation needs to be taken into account explicitely, e.g. when grain boundary sliding or damage initiation at grain boundaries or phase bounaries has to be considered. By adjusting the cohesive zone parameters for grain boundary sliding and opening the competing mechanisms of bulk material deformation and grain boundary accommodated deformation can be studied. Furthermore, it is also possible to investigate damage nucleation at GB triple junctions.

We follow Wei & Anand (2004) to generate a rate independent cohesive zone (CZ) modelling approach for the reasons of simplicity. The velocity jump across a cohesive surface has been additively decomposed into a elastic and a plastic part as follow

$$
\dot{\mathbf{u}} = \dot{\mathbf{u}}\_{\mathbf{e}} + \dot{\mathbf{u}}\_{\mathbf{P}}.\tag{38}
$$

The elastic relative velocities are connected with its power-conjugate traction rate by the interface elastic stiffness tensor

$$
\dot{\mathbf{t}} = \mathbf{K}\dot{\mathbf{u}}\_{\mathbf{e}} = \mathbf{K} \left( \dot{\mathbf{u}} - \dot{\mathbf{u}}\_{\mathbf{P}} \right) \tag{39}
$$

For some special grain boundaries there may exist glide anisotropy inside the grain boundary plane, although our knowledge about this topic is far away from formulating this anisotropy for general grain boundaries. So, we have to assume isotropic plastic deformation property inside the grain boundary, and the trace vector, displacement vector and resistance vector are

According to this approach the material points of the real RVE are approximated to inclusions inside one homogeneous matrix, the property of which can be determined as volume average of these inclusions. After the regular discretisation of the RVE, in the current configuration

Scale Bridging Modeling of Plastic Deformation and Damage Initiation in Polycrystals 15

Assuming we deal with a deformation control process, at the beginning of the iteration loop

where Δ� is the given fixed strain increment. Because each material point has the same volume

where *N* is the total number of inclusions. At this stage, for the inclusions the strain field satisfies deformation compatibility while the stress field does't satisfy the stress equilibrium.

The polarized stress increment field <sup>Δ</sup>σ*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>σ*<sup>M</sup>* can cause a strain increment field in the matrix. Because the matrix material is homogeneous and suffering a periodic boundary condition, this strain increment can be calculated efficiently with help of Green's function

> ⎧ ⎪⎨

1 **x** = **x**�

0 **x** �= **x**�

) in *m* plane along *i* direction applying at **x**� will cause a displacement

**<sup>C</sup>**¯ *ijklGkm*,*lj* + *<sup>δ</sup>mi* = 0. (46)

*km* = *<sup>δ</sup>im or* **<sup>G</sup>**<sup>ˆ</sup> = **<sup>A</sup>**<sup>−</sup>1. (48)

*kmξlξ<sup>j</sup>* + *δmi* = 0 (47)

⎪⎩

where *<sup>ξ</sup>* represents the frequency. Through defining a second order tensor *Aik* = **<sup>C</sup>**¯ *ijklξlξ<sup>j</sup>* we

When one transfers back **G**ˆ from the frequency space to the real physical space one can get the solution **G** of equation 46. However, this is not necessary because we will solve the

1 *<sup>N</sup>* **<sup>C</sup>***<sup>I</sup>*

**<sup>C</sup>**¯ <sup>=</sup> ∑ **x**�∈*RVE*

, and the matrix material point at the same location has stress increment Δσ*M*,

, strain increment Δ�*<sup>I</sup>* and stiffness **C***<sup>I</sup>* =

. (44)

(45)

Δ�*<sup>I</sup>* = Δ�*<sup>M</sup>* = Δ� (43)

the inclusion at location **x**� has stress increment Δσ*<sup>I</sup>*

and shape, the matrix stiffness can be determined simply as

**4.2.1 Stress and strain increment of matrix material points**

and discrete Fourier transformation. With the help of the delta function

*δ*(**x** − **x**�

In order to solve equation 46 we have to transfer it into the frequency space

<sup>−</sup> **<sup>C</sup>**¯ *ijklG*<sup>ˆ</sup>

*km* in the frequency space

displacement field caused by the polarised stress in the frequency space.

*AikG*ˆ

) at **x** satisfying the stress equilibrium

) =

strain increment Δ�*<sup>M</sup>* and stiffness **C**¯ .

*∂*Δσ*I*/*∂*Δ�*<sup>I</sup>*

a unit force *δmi*(**x** − **x**�

find the displacement *G*ˆ

field *Gkm*(**x** − **x**�

defined in the local coordinate [**t***I*, **t***I I*, **n**], where **n** is aligned with the normal to the interface, **t***<sup>I</sup>* and **t***I I* in the tangent plane at the point of the interface under consideration.

According to the rate independent assumption for the loading condition, the hardening rate **s**˙ has to be fast enough to balance the load **t**˙

$$
\dot{\mathbf{t}} = \dot{\mathbf{s}} \tag{40}
$$

or

$$\mathbf{K} \left( \dot{\mathbf{u}} - \dot{\mathbf{u}}\_{\rm P} \right) = \mathbf{H} \dot{\mathbf{u}}\_{\rm P} \tag{41}$$

where **H** is the hardening matrix and its components are state variables of the cohesive zone model. The evolution law of the hardening matrix

$$
\dot{\mathbf{H}} = \dot{\mathbf{H}} \left( \mathbf{H}, \dot{\mathbf{u}}\_{\mathbb{P}} \right) \tag{42}
$$

can be obtained by experimental date fitting or can directly come from the molecular dynamic and ab-initial calculations. Here we adopted the phenomenological hardening rule proposed by Wei & Anand (2004) for the cohesive zone model along the normal direction. For the cohesive zone model along the tangential direction, the failure displacement has been set to infinite to consider the grain boundary glide phenomenon, see Figure 5.
