**4. Numerical approaches**

Starting from the stress equilibrium state, for a given time step and velocity gradient one has to calculate the new stress state, while at the same time considering the evolution of state variables including plastic shear amount *γ* determined by equation 29, statistically stored dislocation density *ρ*SSD determined by equation 30 and geometrically necessary dislocation density *ρ*GND determined by equation 34 for each slip system, the transformed volume fraction *f* determined by equation 36 for each transformation system and the hardening matrix **H** determined by equation 42 for the cohesive zone model of grain boundaries.

#### **4.1 Finite element method**

Based on the Abaqus platform (ABAQUS, 2009), we have developed the user material subroutines UMAT for the bulk material and UINTER for the grain boundary to solve the stress equilibrium and state variable evolution problems. In this approach, except for the plastic strain gradient used for the geometrically necessary dislocation density calculation, which is adopted from last converged time point, all of state variables are calculated by an implicit method.

#### **4.2 Discrete fast Fourier transformation method**

If the representative volume element has a very complicated microstructure and obeys a periodic boundary condition, the stress equilibrium and state variable evolution problems can be solved by the discrete fast Fourier transformation (FFT) method proposed in the literature (Lebensohn, 2001; Michel et al., 2000; 2001).

12 Will-be-set-by-IN-TECH

defined in the local coordinate [**t***I*, **t***I I*, **n**], where **n** is aligned with the normal to the interface,

According to the rate independent assumption for the loading condition, the hardening rate **s**˙

**t**˙ = **s**˙ (40)

= **Hu**˙ <sup>p</sup> (41)

(42)

**t***<sup>I</sup>* and **t***I I* in the tangent plane at the point of the interface under consideration.

**K**

infinite to consider the grain boundary glide phenomenon, see Figure 5.

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

**H**˙ = **H**˙

where **H** is the hardening matrix and its components are state variables of the cohesive zone

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

Starting from the stress equilibrium state, for a given time step and velocity gradient one has to calculate the new stress state, while at the same time considering the evolution of state variables including plastic shear amount *γ* determined by equation 29, statistically stored dislocation density *ρ*SSD determined by equation 30 and geometrically necessary dislocation density *ρ*GND determined by equation 34 for each slip system, the transformed volume fraction *f* determined by equation 36 for each transformation system and the hardening matrix **H** determined by equation 42 for the cohesive zone model of grain boundaries.

Based on the Abaqus platform (ABAQUS, 2009), we have developed the user material subroutines UMAT for the bulk material and UINTER for the grain boundary to solve the stress equilibrium and state variable evolution problems. In this approach, except for the plastic strain gradient used for the geometrically necessary dislocation density calculation, which is adopted from last converged time point, all of state variables are calculated by an

If the representative volume element has a very complicated microstructure and obeys a periodic boundary condition, the stress equilibrium and state variable evolution problems can be solved by the discrete fast Fourier transformation (FFT) method proposed in the literature

**H**, **u**˙ p

has to be fast enough to balance the load **t**˙

**4. Numerical approaches**

**4.1 Finite element method**

**4.2 Discrete fast Fourier transformation method**

(Lebensohn, 2001; Michel et al., 2000; 2001).

implicit method.

model. The evolution law of the hardening matrix

or

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 the inclusion at location **x**� has stress increment Δσ*<sup>I</sup>* , strain increment Δ�*<sup>I</sup>* and stiffness **C***<sup>I</sup>* = *∂*Δσ*I*/*∂*Δ�*<sup>I</sup>* , and the matrix material point at the same location has stress increment Δσ*M*, strain increment Δ�*<sup>M</sup>* and stiffness **C**¯ .

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

$$
\Delta \epsilon^I = \Delta \epsilon^M = \overline{\Delta \epsilon} \tag{43}
$$

where Δ� is the given fixed strain increment. Because each material point has the same volume and shape, the matrix stiffness can be determined simply as

$$
\tilde{\mathbf{C}} = \sum\_{\mathbf{x'} \in \mathbb{R}VE} \frac{1}{N} \mathbf{C}^{I}. \tag{44}
$$

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.

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

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 and discrete Fourier transformation. With the help of the delta function

$$\delta(\mathbf{x} - \mathbf{x}') = \begin{cases} 1 & \mathbf{x} = \mathbf{x}' \\\\ 0 & \mathbf{x} \neq \mathbf{x}' \end{cases} \tag{45}$$

a unit force *δmi*(**x** − **x**� ) in *m* plane along *i* direction applying at **x**� will cause a displacement field *Gkm*(**x** − **x**� ) at **x** satisfying the stress equilibrium

$$
\vec{\mathbb{C}}\_{ijkl}\mathbf{G}\_{km\,lj} + \delta\_{mi} = 0. \tag{46}
$$

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

$$-\vec{\mathsf{C}}\_{ijkl}\hat{\mathsf{G}}\_{km}\mathfrak{xi}\_{l}\mathfrak{f}\_{\mathfrak{f}} + \delta\_{mi} = 0\tag{47}$$

where *<sup>ξ</sup>* represents the frequency. Through defining a second order tensor *Aik* = **<sup>C</sup>**¯ *ijklξlξ<sup>j</sup>* we find the displacement *G*ˆ *km* in the frequency space

$$A\_{ik}\hat{\mathbf{G}}\_{km} = \delta\_{im} \quad or \quad \hat{\mathbf{G}} = \mathbf{A}^{-1}.\tag{48}$$

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 displacement field caused by the polarised stress in the frequency space.

Because each inclusion only has interaction with the matrix, for a fixed matrix property,

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

material point we can calculate the local strain increment which can keep the total strain

and with help of equation 57, the stress increment Δσ of inclusion at **x**� can be easily

In Lebensohn (2001) the deformation compatibility problem and the energy minimisation problems are joined together through adopting the Lagrange multiplier method. In this work, after the stress and strain calculation for inclusions we simply set the matrix material point

Δ�*<sup>M</sup>* = Δ�*<sup>I</sup>*

The final stress and strain solution for inclusions satisfying stress equilibrium and strain

Fig. 6. Local stress and strain patterns of RVEs having grain boundaries with with different

**<sup>C</sup>**¯ <sup>Δ</sup>�*<sup>M</sup>* <sup>−</sup> <sup>σ</sup>*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>σ*<sup>I</sup>*

−<sup>1</sup>

. Under this condition, for each

. (58)

<sup>≤</sup> *Crsd* (59)

(57)

material points at **x** and at **x**� are independent if **x** �= **x**�

Δ�*<sup>I</sup>* = **C***<sup>I</sup>* + **C**¯

> ∑ **x**�∈*RVE*

where *k* and *k* + 1 are iteration numbers and *Crsd* the critical residual.

1 *N* Δ�*<sup>I</sup> <sup>k</sup>* <sup>−</sup> <sup>Δ</sup>�*<sup>I</sup> k*+1 

energy increment minimisation

recalculated by the constitutive law.

compatibility will not be achieved until

**4.2.3 Deformation compatibility**

deformation increment as

**5. Instructive examples**

properties.

Based on the solution of equation 46, the displacement Δ**u***<sup>M</sup>* with respect to the polarised stress field <sup>Δ</sup>σ*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>σ*<sup>M</sup>* can be calculated with the help of Gauss's Theorem and the periodic arrangement of the RVE

$$
\Delta u\_k^M = \sum\_{\mathbf{x'} \in RVE} \mathcal{G}\_{km,n} \left( \Delta \sigma\_{mn}^I - \Delta \sigma\_{mn}^M \right) \Omega \tag{49}
$$

where Ω is the volume of the regular element inside the RVE. Furthermore we can calculate the gradient of the displacement

$$
\Delta u\_{k\rho}^{M} = \sum\_{\mathbf{x}' \in RVE} \mathbb{G}\_{km, \text{no}} \left( \Delta \sigma\_{mn}^{I} - \Delta \sigma\_{mn}^{M} \right) \Omega \tag{50}
$$

and transfer equation 50 into the frequency space, with the help of the convolution theorem:

$$
\Delta\boldsymbol{\hat{u}}\_{k,o}^{M} = -\sum\_{\mathbf{x'} \in \mathcal{R}VE} \hat{\mathbf{G}}\_{km} \mathbb{E}\_{n} \mathbb{E}\_{o} \left( \Delta\boldsymbol{\hat{\sigma}}\_{mn}^{I} - \Delta\boldsymbol{\hat{\sigma}}\_{mn}^{M} \right) \boldsymbol{\Omega}.\tag{51}
$$

Through taking equation 48 into equation 51 we can calculate the displacement gradient in the frequency space

$$\begin{split} \Delta \hat{\mu}\_{k,o}^{M} &= -\sum\_{\mathbf{x'} \in RVE} A\_{km}^{-1} \mathbb{S}\_{n} \mathbb{S}\_{o} \left( \Delta \hat{\sigma}\_{mn}^{I} - \Delta \hat{\sigma}\_{mn}^{M} \right) \Omega \\ &= -\sum\_{\mathbf{x'} \in RVE} \hat{\mathbf{S}}\_{komn} \left( \Delta \hat{\sigma}\_{mn}^{I} - \Delta \hat{\sigma}\_{mn}^{M} \right) \Omega. \end{split} \tag{52}$$

where **S**ˆ *komn* is the compliance tensor in the frequency space for the material point at the real space **x**� . Finally, after one transfers equation 52 from frequency space back to the real physical space one gets the stress and strain increments for the matrix material points

$$
\Delta \epsilon\_{i\dot{j}}^M = \frac{1}{2} (\Delta u\_{i,\dot{j}}^M + \Delta u\_{j,\dot{i}}^M) \tag{53}
$$

$$
\Delta \sigma\_{ij}^{M} = \vec{\mathbb{C}}\_{ijkl}^{-1} \Delta \epsilon\_{kl}^{M} \,. \tag{54}
$$

#### **4.2.2 Stress and strain increment of inclusions**

For each material point at **<sup>x</sup>** when there is a strain misfit between matrix and inclusion <sup>Δ</sup>�*<sup>I</sup>* �<sup>=</sup> Δ�*<sup>M</sup>* there will be a internal misfit stress **C**¯ <sup>Δ</sup>�*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>�*<sup>M</sup>* . The total strain energy increment amounts to

$$
\Delta E = \sum\_{\mathbf{x'} \in \mathbb{R}VE} \left[ \left( \sigma^I + \Delta \sigma^I \right) \cdot \Delta \epsilon^I + \bar{\mathbb{C}} \left( \Delta \epsilon^I - \Delta \epsilon^M \right) \cdot \left( \Delta \epsilon^I - \Delta \epsilon^M \right) \Omega \right]. \tag{55}
$$

Now the procedure of finding stress equilibrium can be replaced with the procedure of achieving the total strain energy increment minimisation

$$\nabla \cdot \sigma^I = \mathbf{0} \iff \frac{\partial \Delta E}{\partial \Delta \epsilon^I} = 0. \tag{56}$$

Because each inclusion only has interaction with the matrix, for a fixed matrix property, material points at **x** and at **x**� are independent if **x** �= **x**� . Under this condition, for each material point we can calculate the local strain increment which can keep the total strain energy increment minimisation

$$
\Delta \boldsymbol{\epsilon}^{I} = \left(\boldsymbol{\mathsf{C}}^{I} + \boldsymbol{\mathsf{\bar{C}}}\right)^{-1} \left(\boldsymbol{\mathsf{\bar{C}}} \Delta \boldsymbol{\epsilon}^{M} - \boldsymbol{\sigma}^{I} - \Delta \boldsymbol{\sigma}^{I}\right) \tag{57}
$$

and with help of equation 57, the stress increment Δσ of inclusion at **x**� can be easily recalculated by the constitutive law.

#### **4.2.3 Deformation compatibility**

14 Will-be-set-by-IN-TECH

Based on the solution of equation 46, the displacement Δ**u***<sup>M</sup>* with respect to the polarised stress field <sup>Δ</sup>σ*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>σ*<sup>M</sup>* can be calculated with the help of Gauss's Theorem and the periodic

> *Gkm*,*<sup>n</sup>* Δ*σ<sup>I</sup>*

*Gkm*,*no* Δ*σ<sup>I</sup>*

*G*ˆ *kmξnξ<sup>o</sup>*

and transfer equation 50 into the frequency space, with the help of the convolution theorem:

Through taking equation 48 into equation 51 we can calculate the displacement gradient in

*A*−<sup>1</sup> *km ξnξ<sup>o</sup>*

**S**ˆ *komn* Δ*σ*ˆ *<sup>I</sup>*

 Δ*σ*ˆ *<sup>I</sup>*

> Δ*σ*ˆ *<sup>I</sup>*

*komn* is the compliance tensor in the frequency space for the material point at the real

. Finally, after one transfers equation 52 from frequency space back to the real physical

*<sup>i</sup>*,*<sup>j</sup>* <sup>+</sup> <sup>Δ</sup>*u<sup>M</sup>*

*ijkl*Δ*�<sup>M</sup>*

<sup>Δ</sup>�*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>�*<sup>M</sup>*

<sup>Δ</sup>�*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>�*<sup>M</sup>*

*∂*Δ*E*

· 

where Ω is the volume of the regular element inside the RVE. Furthermore we can calculate

*mn* <sup>−</sup> <sup>Δ</sup>*σ<sup>M</sup> mn* 

*mn* <sup>−</sup> <sup>Δ</sup>*σ<sup>M</sup> mn* 

> *mn* <sup>−</sup> <sup>Δ</sup>*σ*<sup>ˆ</sup> *<sup>M</sup> mn*

*mn* <sup>−</sup> <sup>Δ</sup>*σ*<sup>ˆ</sup> *<sup>M</sup> mn* Ω

*mn* <sup>−</sup> <sup>Δ</sup>*σ*<sup>ˆ</sup> *<sup>M</sup> mn*  Ω (49)

Ω (50)

Ω. (51)

Ω. (52)

*<sup>j</sup>*,*i*) (53)

. The total strain energy increment

Ω 

. (55)

<sup>Δ</sup>�*<sup>I</sup>* <sup>−</sup> <sup>Δ</sup>�*<sup>M</sup>*

*<sup>∂</sup>*Δ�*<sup>I</sup>* <sup>=</sup> 0. (56)

*kl* . (54)

arrangement of the RVE

the frequency space

where **S**ˆ

space **x**�

amounts to

the gradient of the displacement

Δ*u<sup>M</sup>*

Δ*u<sup>M</sup>*

Δ*u*ˆ*<sup>M</sup>*

Δ*u*ˆ*<sup>M</sup>*

**4.2.2 Stress and strain increment of inclusions**

Δ�*<sup>M</sup>* there will be a internal misfit stress **C**¯

achieving the total strain energy increment minimisation

σ*<sup>I</sup>* + Δσ*<sup>I</sup>*

**x**�∈*RVE*

Δ*E* = ∑

*<sup>k</sup>* = ∑ **x**�∈*RVE*

*<sup>k</sup>*,*<sup>o</sup>* = ∑ **x**�∈*RVE*

*<sup>k</sup>*,*<sup>o</sup>* = − ∑

*<sup>k</sup>*,*<sup>o</sup>* = − ∑

= − ∑ **x**�∈*RVE*

**x**�∈*RVE*

**x**�∈*RVE*

space one gets the stress and strain increments for the matrix material points

Δ*σ<sup>M</sup>*

*ij* <sup>=</sup> **<sup>C</sup>**¯ <sup>−</sup><sup>1</sup>

For each material point at **<sup>x</sup>** when there is a strain misfit between matrix and inclusion <sup>Δ</sup>�*<sup>I</sup>* �<sup>=</sup>

Now the procedure of finding stress equilibrium can be replaced with the procedure of

· <sup>Δ</sup>�*<sup>I</sup>* <sup>+</sup> **<sup>C</sup>**¯

∇ · <sup>σ</sup>*<sup>I</sup>* <sup>=</sup> **<sup>0</sup>** ⇐⇒

Δ*�<sup>M</sup> ij* <sup>=</sup> <sup>1</sup> 2 (Δ*u<sup>M</sup>* In Lebensohn (2001) the deformation compatibility problem and the energy minimisation problems are joined together through adopting the Lagrange multiplier method. In this work, after the stress and strain calculation for inclusions we simply set the matrix material point deformation increment as

$$
\Delta \boldsymbol{\epsilon}^{M} = \Delta \boldsymbol{\epsilon}^{I}.\tag{58}
$$

The final stress and strain solution for inclusions satisfying stress equilibrium and strain compatibility will not be achieved until

$$\sum\_{\mathbf{x'} \in RVE} \frac{1}{N} \left\| \Delta \mathbf{e}\_k^I - \Delta \mathbf{e}\_{k+1}^I \right\| \le \mathcal{C}\_{rsd} \tag{59}$$

where *k* and *k* + 1 are iteration numbers and *Crsd* the critical residual.
