**5. Constitutive models**

### **5.1 Bulk grains**

Isotropic elasticity, anisotropic elasticity and anisotropic elasticity with crystal plasticity constitutive laws are commonly used for bulk grains. Since isotropic elasticity can not account for the effects due to different crystallographic orientations of the grains this is not covered here. Overview of the other two constitutive laws is given below.

362 grain wire 1334 grain wire

*<sup>σ</sup>El*. *ij* <sup>=</sup> *Cijkl* · *�El*. *kl* (1)

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪⎪⎪⎪⎪⎪⎨

*�*11 *�*22 *�*33 2*�*<sup>12</sup> 2*�*<sup>23</sup> 2*�*<sup>31</sup> ⎫ ⎪⎪⎪⎪⎪⎪⎬ *El*.

(2)

*i* ,

*p i*,*j* , due

⎪⎪⎪⎪⎪⎪⎭

*<sup>j</sup>* (3)

⎪⎪⎪⎪⎪⎪⎩

Solid elements 796 105 2 976 828 Cohesive elements 107 094 418 941 All elements 903 199 3 395 769

Grain-Scale Modeling Approaches for Polycrystalline Aggregates 61

Min angle < 5 [◦] 32298 (0.2886 %) 174 (0.0058 %) Max angle > 170 [◦] 45 (0.0056 %) 12 (0.0004 %) Aspect ratio > 10 2449 (0.3076 %) 152 (0.0051 %)

Worst min angle [◦] 0.09 0.47 Worst max angle [◦] 178.77 178.64 Worst aspect ratio 647.1 23.42

In general, the relation between the elastic stress tensor, *<sup>σ</sup>El*. *ij* , and the elastic strain tensor, *�El*. *kl*

Here, the *Cijkl* stands for the stiffness tensor. For materials with a cubic lattice there are only three independent values in the stiffness tensor: *Ciiii* = *C*11, *Ciijj* = *C*<sup>12</sup> and *Cijij* = *C*44. The

> *C*<sup>11</sup> *C*<sup>12</sup> *C*<sup>12</sup> 000 *C*<sup>12</sup> *C*<sup>11</sup> *C*<sup>12</sup> 000 *C*<sup>12</sup> *C*<sup>12</sup> *C*<sup>11</sup> 000 000 *C*<sup>44</sup> 0 0 0000 *C*<sup>44</sup> 0 00000 *C*<sup>44</sup>

Crystal plasticity theory (Hill & Rice, 1972; Rice, 1970) assumes that the plastic deformation in monocrystals takes place via a simple shear on a specific set of slip planes. Deformation by other mechanisms such as for example diffusion, twinning and grain boundary sliding is here not taken into account. The combination of a slip plane, denoted by its normal *m<sup>α</sup>*

to a crystallographic slip can be written using Eq. (3) (Needleman, 2000). The summation is

*γ*˙(*α*) *s* (*α*) *<sup>i</sup> <sup>m</sup>*(*α*)

performed over all active slip systems,(*α*), with *γ*˙(*α*) representing the shear rate.

*u*˙ *p <sup>i</sup>*,*<sup>j</sup>* <sup>=</sup> <sup>∑</sup>*<sup>α</sup>*

*<sup>i</sup>* , is called a slip system, (*α*). The plastic velocity gradient, *u*˙

*Number of elements*

*Values of*

for anisotropic elasticity is given by Eq. (1).

⎧ ⎪⎪⎪⎪⎪⎪⎨

*σ*11 *σ*22 *σ*33 *σ*12 *σ*23 *σ*31 ⎫ ⎪⎪⎪⎪⎪⎪⎬ *El*.

⎪⎪⎪⎪⎪⎪⎭

⎪⎪⎪⎪⎪⎪⎩

**5.1.2 Crystal plasticity**

and a slip direction, *s<sup>α</sup>*

**5.1.1 Anisotropic elasticity**

*Number of elements with*

Table 2. Wire FE models: mesh quality factors comparison.

relation between the strains and stresses is then given by Eq. (2).

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

=

Fig. 11. A FE model of the 400 *μ*m diameter stainless steel wire. Top: wire from the first experimental series, 362 grains. Bottom: wire from the second experimental series, 1334 grains.


Table 2. Wire FE models: mesh quality factors comparison.

#### **5.1.1 Anisotropic elasticity**

12 Polycrystalline Materials

Fig. 11. A FE model of the 400 *μ*m diameter stainless steel wire. Top: wire from the first experimental series, 362 grains. Bottom: wire from the second experimental series, 1334

grains.

In general, the relation between the elastic stress tensor, *<sup>σ</sup>El*. *ij* , and the elastic strain tensor, *�El*. *kl* for anisotropic elasticity is given by Eq. (1).

$$
\sigma\_{ij}^{El} = \mathbb{C}\_{ijkl} \cdot \mathbf{e}\_{kl}^{El} \tag{1}
$$

Here, the *Cijkl* stands for the stiffness tensor. For materials with a cubic lattice there are only three independent values in the stiffness tensor: *Ciiii* = *C*11, *Ciijj* = *C*<sup>12</sup> and *Cijij* = *C*44. The relation between the strains and stresses is then given by Eq. (2).

$$
\begin{Bmatrix}
\sigma\_{11} \\
\sigma\_{22} \\
\sigma\_{33} \\
\sigma\_{12} \\
\sigma\_{23} \\
\sigma\_{31} \\
\end{Bmatrix}^{El} = \begin{bmatrix}
\mathbb{C}\_{11}\mathbb{C}\_{12}\mathbb{C}\_{12} & 0 & 0 & 0 \\
\mathbb{C}\_{12}\mathbb{C}\_{11}\mathbb{C}\_{12} & 0 & 0 & 0 \\
\mathbb{C}\_{12}\mathbb{C}\_{12}\mathbb{C}\_{11} & 0 & 0 & 0 \\
0 & 0 & 0 & \mathbb{C}\_{44} & 0 & 0 \\
0 & 0 & 0 & 0 & \mathbb{C}\_{44} & 0 \\
0 & 0 & 0 & 0 & 0 & \mathbb{C}\_{44}
\end{bmatrix} \begin{Bmatrix}
\varepsilon\_{11} \\
\varepsilon\_{22} \\
\varepsilon\_{33} \\
2\varepsilon\_{12} \\
2\varepsilon\_{23} \\
2\varepsilon\_{31}
\end{Bmatrix} \tag{2}
$$

#### **5.1.2 Crystal plasticity**

Crystal plasticity theory (Hill & Rice, 1972; Rice, 1970) assumes that the plastic deformation in monocrystals takes place via a simple shear on a specific set of slip planes. Deformation by other mechanisms such as for example diffusion, twinning and grain boundary sliding is here not taken into account. The combination of a slip plane, denoted by its normal *m<sup>α</sup> i* , and a slip direction, *s<sup>α</sup> <sup>i</sup>* , is called a slip system, (*α*). The plastic velocity gradient, *u*˙ *p i*,*j* , due to a crystallographic slip can be written using Eq. (3) (Needleman, 2000). The summation is performed over all active slip systems,(*α*), with *γ*˙(*α*) representing the shear rate.

$$\boldsymbol{\dot{n}}\_{i,j}^{p} = \sum\_{\alpha} \boldsymbol{\dot{\gamma}}^{(\alpha)} \boldsymbol{s}\_{i}^{(\alpha)} \boldsymbol{m}\_{j}^{(\alpha)} \tag{3}$$

**5.2 Grain boundaries with cohesive zone approach**

unloaded positions, Eq. (11).

Based on the experimental observations, e.g. (Coffman & Sethna, 2008), grain boundaries are modeled with a cohesive-zone approach using cohesive elements. The traction-separation constitutive behaviour with the damage initiation and evolution as implemented in ABAQUS are used in this work. The cohesive elements are inter-surface elements of often negligible thickness, which essentially measure the relative displacements of the surfaces of adjoining continuum elements. The strains *�* in the cohesive elements are defined using the constitutive thickness of the element *T*<sup>0</sup> (mostly different from the geometric thickness which is typically close or equal to zero) and the separations of the element nodes *δ* as compared to their initial

Grain-Scale Modeling Approaches for Polycrystalline Aggregates 63

⎧ ⎨ ⎩

tractions on the cohesive elements are then given by Eq. (12). ⎧ ⎨ ⎩

> ..... .... .... ... ... . . . . . .

...

...

..

..

..

...

..

.. .. .. .. .. .. .. .. .. .... .. .. .. .. ....................

..

...

...

..

...

...

..

..

..

..

..

.

...

...

..

..

...

.. ..

..

..

...

...

..

...

..

..

...

..

...

..

..

..

..

...

...

..

...

• •

...

...

...

..

....

.

...

..

..

....

..

....

..

...

...

..

...

.

.

.

.

.

.

.

.

.

Fig. 12. Example of traction-separation response (not to scale).

Eq. (13) for the normal direction (and both shear directions).

.

.

.

.

*δ*0

*D*(*δ*) =

then be [1 − *D*(*δ*)] *Knn* and correspondingly for the two shear directions.

.

.

.

.

.

.

.

.

.

.

⎧ ⎨ ⎩

*δ f <sup>n</sup>* (*δ*−*δ*<sup>0</sup> *n*)

*δ* � *δf <sup>n</sup>*−*δ*<sup>0</sup> *n*

.

.

.

.

.

Separation *δn*, *δs*, *δ<sup>t</sup>* [mm]

.

.

.

.

.

.............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................. ........ ...... ..... ..

Typical traction-separation response is given by Fig. 12. Damage evolution *D*(*δ*) is defined by

The actual load-carrying capability of the cohesive element in the normal direction would

.

.

.

.

.

.

*<sup>n</sup> <sup>δ</sup> <sup>f</sup>*

.

.

0 ; *δ* < *δ*<sup>0</sup>

� ; *<sup>δ</sup>* <sup>≥</sup> *<sup>δ</sup>*<sup>0</sup> *n*

.

.

• •

. ... .

.

... .

...

... .

.

.

... . ..

...

... .

. ... . .. .

.

.. . ..

... . ... . ..

... . ..

.. . ..

... . ..

. .. . ..

. ...

*n*

⎫ ⎬ ⎭

. ... . .. . ...

.. . ..

. .. . ..

..

..

. ... . .. . ..

. ...

. ... . ..

... . .. . ..

. ... . . .

... specific elastic energy ... specific fracture energy

*n*

..

. ..

...

... . . .

.. . ... . ..

..........

....................... ...... ...... ..

. ...

.

.

. ..

..................................................

. ... . .. .

.... .. . ..

...........................................................................................................................................................................................................................................................................................................................................................

*t* 0 *n*

.................................. ...... ... ... ... ... ... ... ... ... ... ... ..

Traction *tn*, *ts*, *tt*

•

..

. . . . . . ... . ...

..

... . ... . . ..

. ..

. .. ..

...

.. .

...

. ..

. .. .. . .. . ... .

*E*0

..

. ..

. .. .

. .. . .. . .. .

. ..

.

.

.

.

*n*

.

.

.

.

..

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

•

.

.

.

.

.

.

.

.

.

.

.

*δδ <sup>f</sup> <sup>δ</sup> <sup>n</sup>* <sup>0</sup>

.

*<sup>n</sup> En*

.

.

.

.

.

.

.

....................... ...... ..... .. .................................................................................................................................................................... ...................... .....................................

*δ f n*

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

. ..

............................................................................................................................................................................................................................................................................................................................... ........ ...... ...... ..

*E*0 *n En*

..

.

. . . . ... .

. .. *tn ts tt*

⎫ ⎬ ⎭ =

*�n �s �t*

> ⎡ ⎣

⎫ ⎬ <sup>⎭</sup> <sup>=</sup> <sup>1</sup> *T*0

The indices *n*, *s* and *t* denote the normal and two orthogonal shear directions of the cohesive element. The normal direction always points out of the plane of the cohesive element. The

> *Knn* 0 0 0 *Kss* 0 0 0 *Ktt*

−−−−−−−−−−−−−−−− −−−−−−−−−−−−−−

−−−−−−−−−−−− −−−−−−−−−− −−−−−−−−−−

...

...

..

...

..

...

..

..

...

..

..

...

..

..

...

..

..

.. .

..

−−−−−− −−−−−− −−−− −−

⎧ ⎨ ⎩

*δn δs δt*

> ⎤ ⎦

⎧ ⎨ ⎩

*�n �s �t*

⎫ ⎬ ⎭

⎫ ⎬ ⎭

(11)

(12)

(13)

The cumulative slip is defined as *γ* = ∑*<sup>α</sup> t* 0 *<sup>γ</sup>*˙(*α*) *dt*. From the well-known relation for small strain *�ij* = <sup>1</sup> 2 *ui*,*<sup>j</sup>* + *uj*,*<sup>i</sup>* one can obtain the plastic strain rate, Eq. (4). The constitutive relation of the elastic-plastic monocrystal is now given in terms of stress and strain rates as *<sup>σ</sup>*˙ *ij* <sup>=</sup> *Cijkl �*˙*kl* − *�*˙ *p kl* (Needleman, 2000).

$$\dot{\varepsilon}\_{ij}^p = \sum\_a \frac{1}{2} \dot{\gamma}^{(a)} \left( s\_i^{(a)} m\_j^{(a)} + s\_j^{(a)} m\_i^{(a)} \right) \tag{4}$$

It is assumed that the shear rate *γ*˙(*α*) depends on the stress only via the Schmid resolved shear stress, *τ*(*α*), Eq. (5) and Eq. (6). This is a reasonable approximation at room temperature and for ordinary strain rates and pressures (Needleman, 2000). Yielding is then assumed to take place when the Schmid resolved shear stress exceeds the critical shear stress *τ*0.

$$
\gamma^{(a)} = \dot{a}^{(a)} \left( \frac{\tau^{(a)}}{\mathcal{S}^{(a)}} \right) \left| \frac{\tau^{(a)}}{\mathcal{S}^{(a)}} \right|^{n-1} \tag{5}
$$

$$
\pi^{(a)} = s\_i^{(a)} \sigma\_{i\dot{j}} m\_{\dot{j}}^{(a)} \tag{6}
$$

*a*˙ (*α*) represents the reference strain rate, *n* the strain-rate-sensitivity parameter and *g*(*α*) the current strain-hardened state of the crystal. In the limit, as *n* approaches infinity, this power law approaches that of a rate-independent material. The current strain-hardened state *g*(*α*) can be derived from Eq. (7), where *hαβ* are the slip-hardening moduli defined by the adopted hardening law.

$$\boldsymbol{\dot{\xi}}^{(\alpha)} = \sum\_{\beta} \boldsymbol{h}\_{\alpha \beta} \quad \boldsymbol{\dot{\gamma}}^{(\beta)} \tag{7}$$

$$h\_{aa} = \left\{ (h\_0 - h\_s) \operatorname{sech}^2 \left[ \frac{(h\_0 - h\_s) \,\gamma^{(a)}}{\tau\_s - \tau\_0} \right] + h\_s \right\} G \left( \gamma^{(\beta)}; \beta \neq a \right) \tag{8}$$

$$
\mathfrak{h}\_{\mathfrak{a}\mathfrak{F}} = \mathfrak{q}\mathfrak{h}\_{\mathfrak{a}\mathfrak{a}\mathfrak{v}} \qquad (\mathfrak{a} \neq \mathfrak{B}) , \tag{9}
$$

$$G\left(\gamma^{(\beta)}; \beta \neq a\right) = 1 + \sum\_{\beta \neq a} f\_{a\beta} \tanh\left(\frac{\gamma^{(\beta)}}{\gamma\_0}\right) \tag{10}$$

In this work the Bassani (Bassani & Wu, 1991) hardening law is used with the hardening moduli defined with Eqs. (8, 9, 10). Here *h*<sup>0</sup> stands for the initial hardening modulus, *τ*<sup>0</sup> the yield stress (equal to the initial value of the current strength *g*(*α*)(0)) and *τ<sup>s</sup>* a reference stress where large plastic flow initiates (Huang, 1991). *hs* is hardening modulus during easy glide within stage I hardening and *q* is a hardening factor. The function *G* is associated with cross hardening where *γ*<sup>0</sup> is the amount of slip after which the interaction between slip systems reaches the peak strength, and each component *fαβ* represents the magnitude of the strength of a particular slip interaction. This model was implemented as a user-subroutine into the finite element code ABAQUS. Further details on the applicable theory and implementation can be found in (Huang, 1991).

#### **5.2 Grain boundaries with cohesive zone approach**

14 Polycrystalline Materials

relation of the elastic-plastic monocrystal is now given in terms of stress and strain rates as

It is assumed that the shear rate *γ*˙(*α*) depends on the stress only via the Schmid resolved shear stress, *τ*(*α*), Eq. (5) and Eq. (6). This is a reasonable approximation at room temperature and for ordinary strain rates and pressures (Needleman, 2000). Yielding is then assumed to take

> 

*τ*(*α*) *g*(*α*)   *n*−1

*dt*. From the well-known relation for small

*<sup>j</sup>* (6)

*<sup>h</sup>αβ <sup>γ</sup>*˙(*β*) (7)

*hαβ* = *qhαα*, (*α* �= *β*), (9)

 *γ*(*β*) *γ*0

; *β* �= *α*  (4)

(5)

(8)

(10)

one can obtain the plastic strain rate, Eq. (4). The constitutive

 *t* 0 *<sup>γ</sup>*˙(*α*) 

The cumulative slip is defined as *γ* = ∑*<sup>α</sup>*

*ui*,*<sup>j</sup>* + *uj*,*<sup>i</sup>*

(Needleman, 2000).

1 2 *γ*˙(*α*) *s* (*α*) *<sup>i</sup> <sup>m</sup>*(*α*) *<sup>j</sup>* + *s* (*α*) *<sup>j</sup> <sup>m</sup>*(*α*) *i* 

place when the Schmid resolved shear stress exceeds the critical shear stress *τ*0.

(*α*) *τ*(*α*) *g*(*α*)

*τ*(*α*) = *s*

(*α*) *<sup>i</sup> <sup>σ</sup>ijm*(*α*)

(*α*) represents the reference strain rate, *n* the strain-rate-sensitivity parameter and *g*(*α*) the current strain-hardened state of the crystal. In the limit, as *n* approaches infinity, this power law approaches that of a rate-independent material. The current strain-hardened state *g*(*α*) can be derived from Eq. (7), where *hαβ* are the slip-hardening moduli defined by the adopted

> (*h*<sup>0</sup> <sup>−</sup> *hs*) *<sup>γ</sup>*(*α*) *τ<sup>s</sup>* − *τ*<sup>0</sup>

= 1 + ∑ *β*�=*α*

In this work the Bassani (Bassani & Wu, 1991) hardening law is used with the hardening moduli defined with Eqs. (8, 9, 10). Here *h*<sup>0</sup> stands for the initial hardening modulus, *τ*<sup>0</sup> the yield stress (equal to the initial value of the current strength *g*(*α*)(0)) and *τ<sup>s</sup>* a reference stress where large plastic flow initiates (Huang, 1991). *hs* is hardening modulus during easy glide within stage I hardening and *q* is a hardening factor. The function *G* is associated with cross hardening where *γ*<sup>0</sup> is the amount of slip after which the interaction between slip systems reaches the peak strength, and each component *fαβ* represents the magnitude of the strength of a particular slip interaction. This model was implemented as a user-subroutine into the finite element code ABAQUS. Further details on the applicable theory and implementation

 + *hs G γ*(*β*)

*fαβ* tanh

*γ*˙(*α*) = *a*˙

*g*˙ (*α*) <sup>=</sup> ∑ *β*

; *β* �= *α* 

*�*˙ *p ij* <sup>=</sup> <sup>∑</sup>*<sup>α</sup>*

strain *�ij* = <sup>1</sup>

*σ*˙ *ij* = *Cijkl*

*a*˙

hardening law.

*hαα* =

can be found in (Huang, 1991).

*G γ*(*β*)

(*h*<sup>0</sup> <sup>−</sup> *hs*) sech2

2 

 *�*˙*kl* − *�*˙ *p kl*  Based on the experimental observations, e.g. (Coffman & Sethna, 2008), grain boundaries are modeled with a cohesive-zone approach using cohesive elements. The traction-separation constitutive behaviour with the damage initiation and evolution as implemented in ABAQUS are used in this work. The cohesive elements are inter-surface elements of often negligible thickness, which essentially measure the relative displacements of the surfaces of adjoining continuum elements. The strains *�* in the cohesive elements are defined using the constitutive thickness of the element *T*<sup>0</sup> (mostly different from the geometric thickness which is typically close or equal to zero) and the separations of the element nodes *δ* as compared to their initial unloaded positions, Eq. (11).

$$\begin{Bmatrix} \begin{cmatrix} \varepsilon\_{ll} \\ \varepsilon\_{s} \\ \varepsilon\_{l} \end{cmatrix} \end{Bmatrix} = \frac{1}{T\_{0}} \begin{Bmatrix} \delta\_{n} \\ \delta\_{s} \\ \delta\_{t} \end{Bmatrix} \tag{11}$$

The indices *n*, *s* and *t* denote the normal and two orthogonal shear directions of the cohesive element. The normal direction always points out of the plane of the cohesive element. The tractions on the cohesive elements are then given by Eq. (12).

$$\begin{Bmatrix} t\_{\boldsymbol{\eta}} \\ t\_{\boldsymbol{s}} \\ t\_{\boldsymbol{t}} \end{Bmatrix} = \begin{bmatrix} K\_{\boldsymbol{m}\boldsymbol{n}} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & K\_{\boldsymbol{s}\boldsymbol{s}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & K\_{\boldsymbol{t}\boldsymbol{t}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\epsilon}\_{\boldsymbol{n}} \\ \boldsymbol{\epsilon}\_{\boldsymbol{s}} \\ \boldsymbol{\epsilon}\_{\boldsymbol{t}} \end{Bmatrix} \tag{12}$$

Fig. 12. Example of traction-separation response (not to scale).

Typical traction-separation response is given by Fig. 12. Damage evolution *D*(*δ*) is defined by Eq. (13) for the normal direction (and both shear directions).

$$D(\delta) = \begin{cases} 0 & \text{; } \delta < \delta\_n^0 \\ \frac{\delta\_n^{\ell} \left( \delta - \delta\_n^0 \right)}{\delta \left( \delta\_n^{\ell} - \delta\_n^0 \right)} \text{; } \delta \ge \delta\_n^0 \end{cases} \tag{13}$$

The actual load-carrying capability of the cohesive element in the normal direction would then be [1 − *D*(*δ*)] *Knn* and correspondingly for the two shear directions.

The resulting stress tensor is given by Eq. (19).

grain boundaries.

(Avg: 75%) S, S33

> +1.340e+02 +1.422e+02 +1.503e+02 +1.585e+02 +1.667e+02 +1.748e+02 +1.830e+02 +1.912e+02 +1.994e+02 +2.075e+02 +2.157e+02 +2.239e+02 +2.320e+02

(Avg: 75%) S, S33

> +1.339e+02 +1.421e+02 +1.503e+02 +1.585e+02 +1.667e+02 +1.749e+02 +1.831e+02 +1.913e+02 +1.995e+02 +2.077e+02 +2.159e+02 +2.241e+02 +2.323e+02

the triangular prism cohesive elements.

*σij* =

⎡ ⎣

*nGrain*<sup>1</sup>*Grain*<sup>2</sup> =

*nGrain*<sup>2</sup>*Grain*<sup>3</sup> =

*nGrain*<sup>1</sup>*Grain*<sup>3</sup> =

200 0 0 0 100 0 0 0 90

Grain-Scale Modeling Approaches for Polycrystalline Aggregates 65

The boundaries between the grain are defined with the vectors, normal to the planes of the

⎡ ⎣ 1 1 0

⎡ ⎣ -1 1 0

⎡ ⎣ 1 0 0

Since we know the stress tensor and the normals for these three planes, we can compute the stresses in the normal direction for each of them. For the plane between the Grain1 and Grain2

> (Avg: 75%) S, S33

> > +1.500e+02 +1.542e+02 +1.583e+02 +1.625e+02 +1.667e+02 +1.708e+02 +1.750e+02 +1.792e+02 +1.833e+02 +1.875e+02 +1.917e+02 +1.958e+02 +2.000e+02

(Avg: 75%) S, S33

> +1.500e+02 +1.542e+02 +1.583e+02 +1.625e+02 +1.667e+02 +1.708e+02 +1.750e+02 +1.792e+02 +1.833e+02 +1.875e+02 +1.917e+02 +1.958e+02 +2.000e+02

Fig. 14. Normal stresses in the cohesive elements. Triangular prisms (left) and rectangular prisms (right) cohesive elements. Variation in stress at the Y triple points can be observed for

⎤ <sup>⎦</sup> · <sup>1</sup>

⎤ <sup>⎦</sup> · <sup>1</sup>

⎤ <sup>⎦</sup> · <sup>1</sup>

⎤

⎦ MPa (19)

<sup>√</sup><sup>2</sup> (20)

<sup>√</sup><sup>2</sup> (21)

<sup>√</sup><sup>2</sup> (22)

#### **5.2.1 Cohesive elements issues**

In this section we explore the response of the cohesive elements in their normal direction. Let us have a cuboid, divided into three grains as depicted in the Fig. 13 by the three colors. Let us put 100 MPa of tensile stress on the top and the bottom surface and 200 MPa of tensile stress on the left and right surface. Let us constrain the front and the back surface in the Z direction, resulting in *�*33=0. Furthermore, let us assume that we are dealing with isotropic elastic material with Young modulus *E*=200 000 MPa and Poisson ratio of *ν*=0.3. For isotropic elastic material the Eqs. (14,15,16) are valid.

Fig. 13. A simple Y model.

**Z**

$$
\sigma\_{11} = \frac{E}{1+\nu} \epsilon\_{11} + \frac{E\nu}{(1+\nu)(1-2\nu)} (\epsilon\_{11} + \epsilon\_{22} + \epsilon\_{33}) \tag{14}
$$

$$
\sigma\_{22} = \frac{E}{1+\nu} \epsilon\_{22} + \frac{E\nu}{(1+\nu)(1-2\nu)} (\epsilon\_{11} + \epsilon\_{22} + \epsilon\_{33}) \tag{15}
$$

$$
\omega\_{33} = \frac{E}{1+\nu} \epsilon\_{33} + \frac{E\nu}{(1+\nu)(1-2\nu)} (\epsilon\_{11} + \epsilon\_{22} + \epsilon\_{33}) \tag{16}
$$

From Eqs. (14,15), expressions Eq. (17,18) can be derived as *�*33=0 due to the boundary conditions. Using *<sup>σ</sup>*11=200 MPa and *<sup>σ</sup>*22=100 MPa we obtain values *�*<sup>22</sup> <sup>=</sup> 6.5 · <sup>10</sup>−<sup>5</sup> and *�*<sup>11</sup> <sup>=</sup> 7.15 · <sup>10</sup>−4. Using these two values in Eq. (16) we obtain *<sup>σ</sup>*33=90 MPa.

$$\varepsilon\_{22} = \frac{\sigma\_{11}(1+\nu)(1-2\nu) - \sigma\_{22}\frac{(1+\nu)(1-2\nu)(1-\nu)}{\nu}}{E\nu - E\frac{(1-\nu)^2}{\nu}}\tag{17}$$

$$
\varepsilon\_{11} = \frac{\sigma\_{22}(1+\nu)(1-2\nu) - \varepsilon\_{22}E(1-\nu)}{E\nu} \tag{18}
$$

The resulting stress tensor is given by Eq. (19).

16 Polycrystalline Materials

In this section we explore the response of the cohesive elements in their normal direction. Let us have a cuboid, divided into three grains as depicted in the Fig. 13 by the three colors. Let us put 100 MPa of tensile stress on the top and the bottom surface and 200 MPa of tensile stress on the left and right surface. Let us constrain the front and the back surface in the Z direction, resulting in *�*33=0. Furthermore, let us assume that we are dealing with isotropic elastic material with Young modulus *E*=200 000 MPa and Poisson ratio of *ν*=0.3. For isotropic

**100 MPa**

**Grain 3**

**Grain 2**

**200 MPa**

(*�*<sup>11</sup> + *�*<sup>22</sup> + *�*33) (14)

(*�*<sup>11</sup> + *�*<sup>22</sup> + *�*33) (15)

(*�*<sup>11</sup> + *�*<sup>22</sup> + *�*33) (16)

*<sup>E</sup><sup>ν</sup>* (18)

(17)

**5.2.1 Cohesive elements issues**

**X**

**Y**

Fig. 13. A simple Y model.

**Z**

elastic material the Eqs. (14,15,16) are valid.

**Grain 1**

*<sup>σ</sup>*<sup>11</sup> <sup>=</sup> *<sup>E</sup>* 1 + *ν*

*<sup>σ</sup>*<sup>22</sup> <sup>=</sup> *<sup>E</sup>* 1 + *ν*

*<sup>σ</sup>*<sup>33</sup> <sup>=</sup> *<sup>E</sup>* 1 + *ν* *�*<sup>11</sup> +

*�*<sup>22</sup> +

*�*<sup>33</sup> +

*�*<sup>11</sup> <sup>=</sup> 7.15 · <sup>10</sup>−4. Using these two values in Eq. (16) we obtain *<sup>σ</sup>*33=90 MPa.

*Eν* (1 + *ν*)(1 − 2*ν*)

*Eν* (1 + *ν*)(1 − 2*ν*)

*Eν* (1 + *ν*)(1 − 2*ν*)

From Eqs. (14,15), expressions Eq. (17,18) can be derived as *�*33=0 due to the boundary conditions. Using *<sup>σ</sup>*11=200 MPa and *<sup>σ</sup>*22=100 MPa we obtain values *�*<sup>22</sup> <sup>=</sup> 6.5 · <sup>10</sup>−<sup>5</sup> and

*�*<sup>22</sup> <sup>=</sup> *<sup>σ</sup>*11(<sup>1</sup> <sup>+</sup> *<sup>ν</sup>*)(<sup>1</sup> <sup>−</sup> <sup>2</sup>*ν*) <sup>−</sup> *<sup>σ</sup>*<sup>22</sup> (1+*ν*)(1−2*ν*)(1−*ν*)

*�*<sup>11</sup> <sup>=</sup> *<sup>σ</sup>*22(<sup>1</sup> <sup>+</sup> *<sup>ν</sup>*)(<sup>1</sup> <sup>−</sup> <sup>2</sup>*ν*) <sup>−</sup> *�*22*E*(<sup>1</sup> <sup>−</sup> *<sup>ν</sup>*)

*<sup>E</sup><sup>ν</sup>* <sup>−</sup> *<sup>E</sup>* (1−*ν*)<sup>2</sup> *ν*

*ν*

$$
\sigma\_{ij} = \begin{bmatrix} 200 & 0 & 0 \\ 0 & 100 & 0 \\ 0 & 0 & 90 \end{bmatrix} \text{ MPa} \tag{19}
$$

The boundaries between the grain are defined with the vectors, normal to the planes of the grain boundaries.

$$n\_{Grain1Grain2} = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} \cdot \frac{1}{\sqrt{2}} \tag{20}$$

$$m\_{Grain2Grain3} = \begin{bmatrix} \cdot 1 \\ 1 \\ 0 \end{bmatrix} \cdot \frac{1}{\sqrt{2}} \tag{21}$$

$$m\_{Grain1Grain3} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} \cdot \frac{1}{\sqrt{2}} \tag{22}$$

Since we know the stress tensor and the normals for these three planes, we can compute the stresses in the normal direction for each of them. For the plane between the Grain1 and Grain2

Fig. 14. Normal stresses in the cohesive elements. Triangular prisms (left) and rectangular prisms (right) cohesive elements. Variation in stress at the Y triple points can be observed for the triangular prism cohesive elements.

the stress vector on the plane, *p*, and the normal stress, *σn*, are given by Eq. (23) and Eq. (24).

$$\begin{aligned} \text{Grain1Grain2}: \quad p &= \underbrace{\begin{bmatrix} 200 & 0 & 0 \\ 0 & 100 & 0 \\ 0 & 0 & 90 \end{bmatrix}}\_{\sigma\_{ij}} \cdot \underbrace{\begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix}}\_{n\_{Gain1Gain2}} = \begin{bmatrix} 200 \\ 100 \\ 0 \end{bmatrix} \cdot \frac{1}{\sqrt{2}} \text{MPa} \tag{23} \end{aligned} \tag{23}$$

$$\text{Grain1Grain2}: \quad \sigma\_{\text{n}} = p\_{\text{i}} \cdot n\_{\text{i}} = 150 \,\text{MPa} \tag{24}$$

Similarly, we obtain the following normal stress values for the other two planes:

$$\text{Grain2Grain3}: \quad \sigma\_{\text{ll}} = 150 \,\text{MPa} \tag{25}$$

$$
\sigma\_{\text{Grain1Grain3}} : \quad \sigma\_{\text{ll}} = 200 \,\text{MPa} \tag{26}
$$

<sup>−</sup><sup>10</sup> <sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> <sup>100</sup> <sup>−</sup><sup>10</sup>

[MPa]

Theoretical normal stress σ<sup>n</sup>

Grain-Scale Modeling Approaches for Polycrystalline Aggregates 67

Fig. 15. Theoretical and computed normal stresses at integration points of the cohesive

Fig. 16. Cohesive elements with more than 50 % difference between the theoretical and computed normal stresses (in red). 3D Voronoi, 100 grains. Left: element size=0.025. Right:

boundaries, depending upon the crystallographic orientation of the neighboring grains. In this work a simplification is used where resistant grain boundaries are defined as coincidence site lattice (Grimmer et al., 1974) (Σ3 through Σ29) grain boundaries and low angle grain boundaries with misorientation angle between the neighboring grains below 15 o, following (Marrow et al., 2006). Σ value is computed as the ratio enclosed by a unit cell of the

element size=0.0125.

elements. 3D Voronoi, 100 grains, element size=0.025.

Cohesive normal stress σn [MPa]

Fig. 14 displays the normal stresses in the cohesive elements as calculated from the ABAQUS finite element models. On the left hand side triangular prism cohesive elements are used, whereas on the right hand side rectangular prism cohesive elements are used. One can see that in the case of triangular prisms there is a variation in the normal stress for a given plane, in particular at the Y triple points. For the rectangular prism cohesive elements no such variations are observed and the values from the FE model match exactly with the theoretically computed values. Rectangular prism cohesive elements should therefore be preferentially used.

Unfortunately, meshing complicated shapes with rectangular prism cohesive elements also requires that the grains need to be meshed with rectangular prisms-hexahedral elements. With shapes as complicated as seen in Fig. 5 this is extremely difficult and one therefore uses triangular prisms-tetrahedral elements. This then requires the use of triangular prism cohesive elements if one is to obtain a conformal mesh between the structural and cohesive elements. One therefore automatically introduces a degree of discrepancy.

Fig. 15 displays the computed versus the theoretical normal stresses for the cohesive elements for the 3D Voronoi tessellation with 100 grains with external load of 70 MPa. Isotropic elastic constitutive law is used for the grains with the wire loaded in tension. The solid line represents the ideal response. One can observe a significant scatter from the ideal response. Fig. 16 displays the cohesive elements with red color indicating elements with more than 50 % difference between the theoretical and computed normal stresses. One can see that a significant scatter exists. Similar scatter has been reported on the grain structure simulated using 3D Voronoi tessellations (Kamaya & Itakura, 2009). Vast majority of the problematic cohesive elements are located on the triple lines between the grains. The scatter of the normal stresses of the cohesive elements not lying on the triple lines is significantly smaller, with most values within the ±20 % deviation. Increasing mesh density helps to alleviate the issue to some degree by reducing the area of the problematic elements. Other factors such as the stiffness of the cohesive element and its thickness have negligible effect on the scatter (Simonovski & Cizelj, 2011c).

#### **6. Examples**

In this section we look at the some results for the wire model with initialization and propagation of intergranular stress corrosion cracks. The grain boundaries are modeled using the above described cohesive zone approach and classified into resistant and susceptible grain 18 Polycrystalline Materials

the stress vector on the plane, *p*, and the normal stress, *σn*, are given by Eq. (23) and Eq. (24).

Fig. 14 displays the normal stresses in the cohesive elements as calculated from the ABAQUS finite element models. On the left hand side triangular prism cohesive elements are used, whereas on the right hand side rectangular prism cohesive elements are used. One can see that in the case of triangular prisms there is a variation in the normal stress for a given plane, in particular at the Y triple points. For the rectangular prism cohesive elements no such variations are observed and the values from the FE model match exactly with the theoretically computed values. Rectangular prism cohesive elements should therefore be preferentially

Unfortunately, meshing complicated shapes with rectangular prism cohesive elements also requires that the grains need to be meshed with rectangular prisms-hexahedral elements. With shapes as complicated as seen in Fig. 5 this is extremely difficult and one therefore uses triangular prisms-tetrahedral elements. This then requires the use of triangular prism cohesive elements if one is to obtain a conformal mesh between the structural and cohesive

Fig. 15 displays the computed versus the theoretical normal stresses for the cohesive elements for the 3D Voronoi tessellation with 100 grains with external load of 70 MPa. Isotropic elastic constitutive law is used for the grains with the wire loaded in tension. The solid line represents the ideal response. One can observe a significant scatter from the ideal response. Fig. 16 displays the cohesive elements with red color indicating elements with more than 50 % difference between the theoretical and computed normal stresses. One can see that a significant scatter exists. Similar scatter has been reported on the grain structure simulated using 3D Voronoi tessellations (Kamaya & Itakura, 2009). Vast majority of the problematic cohesive elements are located on the triple lines between the grains. The scatter of the normal stresses of the cohesive elements not lying on the triple lines is significantly smaller, with most values within the ±20 % deviation. Increasing mesh density helps to alleviate the issue to some degree by reducing the area of the problematic elements. Other factors such as the stiffness of the cohesive element and its thickness have negligible effect on the scatter

In this section we look at the some results for the wire model with initialization and propagation of intergranular stress corrosion cracks. The grain boundaries are modeled using the above described cohesive zone approach and classified into resistant and susceptible grain

elements. One therefore automatically introduces a degree of discrepancy.

⎤ ⎦ · ⎡ ⎣ 1 1 0 ⎤ <sup>⎦</sup> · <sup>1</sup> √2

� �� � *nGrain*<sup>1</sup>*Grain*<sup>2</sup>

= ⎡ ⎣ 200 100 0

*Grain*1*Grain*2 : *σ<sup>n</sup>* = *pi* · *ni* = 150 MPa (24)

*Grain*2*Grain*3 : *σ<sup>n</sup>* = 150 MPa (25)

*Grain*1*Grain*3 : *σ<sup>n</sup>* = 200 MPa (26)

⎤ <sup>⎦</sup> · <sup>1</sup> √2

MPa (23)

200 0 0 0 100 0 0 0 90

� �� � *<sup>σ</sup>ij*

Similarly, we obtain the following normal stress values for the other two planes:

⎡ ⎣

*Grain*1*Grain*2 : *p* =

used.

(Simonovski & Cizelj, 2011c).

**6. Examples**

Fig. 15. Theoretical and computed normal stresses at integration points of the cohesive elements. 3D Voronoi, 100 grains, element size=0.025.

Fig. 16. Cohesive elements with more than 50 % difference between the theoretical and computed normal stresses (in red). 3D Voronoi, 100 grains. Left: element size=0.025. Right: element size=0.0125.

boundaries, depending upon the crystallographic orientation of the neighboring grains. In this work a simplification is used where resistant grain boundaries are defined as coincidence site lattice (Grimmer et al., 1974) (Σ3 through Σ29) grain boundaries and low angle grain boundaries with misorientation angle between the neighboring grains below 15 o, following (Marrow et al., 2006). Σ value is computed as the ratio enclosed by a unit cell of the

(Avg: 75%) Mises [MPa] (Avg: 75%) SDEG

> 0.00 0.08 0.17 0.25 0.33 0.41 0.50 0.58 0.66 0.74 0.82 0.91 0.99

Grains: AE+CP Final time increment

 5.00 26.45 47.89 69.34 90.79 112.23 133.68 155.13 176.58 198.02 219.47 240.92 262.36

Grains: AE+CP Final time increment

Grains: IE Early time increment Grains: IE Early time increment

Grain-Scale Modeling Approaches for Polycrystalline Aggregates 69

Grains: AE Early time increment Grains: AE Early time increment

Large compressive stress

Fig. 17. Mises stress (left) and damage (right) development. Isotropic elastic (top), anisotropic elastic (middle) and anisotropic elastic + crystal plasticity model (bottom).

coincidence sites and the standard unit cell (Bollman, 1982) from the Rodrigues vectors of the two neighboring grains. Brandon (Brandon, 1996) criterion for a proximity to a coincidence site lattice structure with proportionality constant of 10 <sup>o</sup> is used. All other grain boundaries are defined as susceptible grain boundaries.

First, we constrain the nodes on the back surface in all three directions. Next, we apply tensile stress of 60 MPa to the wire's front surface. This load is equal to the one in the experiment (King et al., 2008a;b). After the application of the load we assume that the wire is exposed to an acidified solution of potassium tetrathionate (*K*2*S*4*O*6) that penetrates into the wire and degrades the susceptible grain boundaries. This is done to mimic the experiment and is accomplished by employing a user-defined field variable of a cylindrical shape with it's axis aligned with the wire's axis.

The radius of the user-defined field variable cylinder is progressively decreased. Once a cohesive element of a susceptible grain boundary lies outside the user-defined field variable cylinder's radius, its *δ*<sup>0</sup> *<sup>n</sup>* and *δδ <sup>f</sup> <sup>n</sup>* values are decreased, resulting in practically instantaneous full degradation.

Damage of a cohesive element of a susceptible grain boundary that is inside the user-defined field variable cylinder's radius is caused only by mechanical overload, as depicted in Fig. 12. The described approach is purely mechanical. To improve convergence, the decrease of *δ*<sup>0</sup> *n* and *δδ <sup>f</sup> <sup>n</sup>* values is not done abruptly. Additionally, a viscous damping of 0.01 is applied to the damage function to improve the convergence during the cohesive elements damage evolution. Susceptible grain boundaries at initial and end 10 % of the wire's length are not allowed to be affected by the user-variable to reduce the edge effect and improve the numerical stability.

The rate at which the radius of the user-defined field variable decreases (i.e. degradation rate) is linked to the stability of the computation. If a large value is selected, then within one computational time increment the separation of the opposite faces in a cohesive element can reach the critical value of *<sup>δ</sup> <sup>f</sup> <sup>n</sup>* at which the element completely degrades. If this degradation process occurs within one computational time increment, the convergence is degraded, even more so when this occurs simultaneously in several cohesive elements. Small enough degradation rates therefore need to be used, so that the process of degradation of cohesive elements is captured within the computational time increments. Alternatively, one can select very small computational time increments or increase the step time. Similarly, the *δ*<sup>0</sup> *<sup>n</sup>* and *δδ <sup>f</sup> n* values should not be very small since at already small load increment the resulting separation of the cohesive element faces could be high enough to instantaneously completely damage

Fig. 17 shows the Mises stresses in bulk grains (left part) and damage evolution on the grain boundaries (right part) for the three constitutive models for grains. 30K case with 10.0 *μ*m element size is presented. The first two rows (with the exception of the legend) display the results at an early time increment where the damage due to the corrosion is very limited. For the Mises stress significant differences can be observed between the IE (isotropic elasticity) and AE (anisotropic elasticity) constitutive laws. This is to be expected since the crystallographic orientations are disregarded in IE approach. Due to the low applied load (less than 1/3rd of the yield stress) almost no difference can be observed between the Mises stresses for the AE and AE+CP (anisotropic elasticity+crystal plasticity) constitutive models. This was true even at the end of the simulation where the damage of the grain

the element, again causing convergence issues.

20 Polycrystalline Materials

coincidence sites and the standard unit cell (Bollman, 1982) from the Rodrigues vectors of the two neighboring grains. Brandon (Brandon, 1996) criterion for a proximity to a coincidence site lattice structure with proportionality constant of 10 <sup>o</sup> is used. All other grain boundaries

First, we constrain the nodes on the back surface in all three directions. Next, we apply tensile stress of 60 MPa to the wire's front surface. This load is equal to the one in the experiment (King et al., 2008a;b). After the application of the load we assume that the wire is exposed to an acidified solution of potassium tetrathionate (*K*2*S*4*O*6) that penetrates into the wire and degrades the susceptible grain boundaries. This is done to mimic the experiment and is accomplished by employing a user-defined field variable of a cylindrical shape with it's

The radius of the user-defined field variable cylinder is progressively decreased. Once a cohesive element of a susceptible grain boundary lies outside the user-defined field variable

Damage of a cohesive element of a susceptible grain boundary that is inside the user-defined field variable cylinder's radius is caused only by mechanical overload, as depicted in Fig. 12. The described approach is purely mechanical. To improve convergence, the decrease of *δ*<sup>0</sup>

process occurs within one computational time increment, the convergence is degraded, even more so when this occurs simultaneously in several cohesive elements. Small enough degradation rates therefore need to be used, so that the process of degradation of cohesive elements is captured within the computational time increments. Alternatively, one can select

values should not be very small since at already small load increment the resulting separation of the cohesive element faces could be high enough to instantaneously completely damage

Fig. 17 shows the Mises stresses in bulk grains (left part) and damage evolution on the grain boundaries (right part) for the three constitutive models for grains. 30K case with 10.0 *μ*m element size is presented. The first two rows (with the exception of the legend) display the results at an early time increment where the damage due to the corrosion is very limited. For the Mises stress significant differences can be observed between the IE (isotropic elasticity) and AE (anisotropic elasticity) constitutive laws. This is to be expected since the crystallographic orientations are disregarded in IE approach. Due to the low applied load (less than 1/3rd of the yield stress) almost no difference can be observed between the Mises stresses for the AE and AE+CP (anisotropic elasticity+crystal plasticity) constitutive models. This was true even at the end of the simulation where the damage of the grain

very small computational time increments or increase the step time. Similarly, the *δ*<sup>0</sup>

*<sup>n</sup>* values is not done abruptly. Additionally, a viscous damping of 0.01 is applied to the damage function to improve the convergence during the cohesive elements damage evolution. Susceptible grain boundaries at initial and end 10 % of the wire's length are not allowed to be affected by the user-variable to reduce the edge effect and improve the numerical stability. The rate at which the radius of the user-defined field variable decreases (i.e. degradation rate) is linked to the stability of the computation. If a large value is selected, then within one computational time increment the separation of the opposite faces in a cohesive element can

*<sup>n</sup>* values are decreased, resulting in practically instantaneous

*<sup>n</sup>* at which the element completely degrades. If this degradation

*n*

*<sup>n</sup>* and *δδ <sup>f</sup> n*

are defined as susceptible grain boundaries.

*<sup>n</sup>* and *δδ <sup>f</sup>*

axis aligned with the wire's axis.

cylinder's radius, its *δ*<sup>0</sup>

reach the critical value of *<sup>δ</sup> <sup>f</sup>*

the element, again causing convergence issues.

full degradation.

and *δδ <sup>f</sup>*

Fig. 17. Mises stress (left) and damage (right) development. Isotropic elastic (top), anisotropic elastic (middle) and anisotropic elastic + crystal plasticity model (bottom).

elasticity and anisotropic elasticity+crystal plasticity approaches, with the latter resulting in

Grain-Scale Modeling Approaches for Polycrystalline Aggregates 71

In all cases almost uniform degradation of the grain boundaries is observed. This is attributed to a) a missing link between the grain boundary load and rate at which the corrosion penetrates the grain boundaries, b) uniform grain boundary stiffnesses whereas stiffness distribution based upon the properties of adjacent grains is expected (Coffman & Sethna, 2008) and c) the selected corrosive environment penetration approach does not account for the topology of the grain boundary network. These issues are the subject of further work.

The numerical stability of the simulation including damage is reasonable, with slightly better convergence for the anisotropic elasticity+crystal plasticity approach. The degradation of a cohesive element is linked to the stability of a simulation. If this degradation process occurs within one computational time increment, the convergence is degraded, even more so when this occurs simultaneously in several cohesive elements. Small computational time increments

*<sup>n</sup>* and *δδ <sup>f</sup>*

already small load increment the resulting separation of the cohesive element faces could be high enough to instantaneously completely damage the element, again causing convergence

The authors would like to thank Prof. James Marrow, (formerly of the University of Manchester, now University of Oxford) and Dr. Andrew King from the European Synchrotron Radiation Facility for kindly providing the results of DCT measurements and also useful discussions. The authors also gratefully acknowledge the financial support from the Slovenian

Aurenhammer, F. (1991). Voronoi diagrams-a survey of a fundamental geometric data

Bassani, J. L. & Wu, T.-Y. (1991). Latent hardening in single crystals. II. Analytical

Bennett, V. P. & McDowell, D. L. (2003). Crack tip displacements of microstructurally small

Bollman, W. (1982). *Crystal Lattices, Interfaces, Matrices*, Polycrystal Book Service, ISBN

Brandon, D. (1996). The structure of high-angle grain boundaries, *Acta Metallurgica*

Cailletaud, G., Forest, S., Jeulin, D., Feyel, F., Galliet, I., Mounoury, V. & Quilici, S.

Coffman, V. & Sethna, J. (2008). Grain boundary energies and cohesive strength as a function

characterization and predictions, *Proceedings: Mathematical and Physical Sciences*

surface cracks in single phase ductile polycrystals, *Engineering Fracture Mechanics*

(2003). Some elements of microstructural mechanics, *Computational Materials Science*

of geometry (http://link.aps.org/doi/10.1103/physrevb.77.144111), *Physical review*

research agency through research programmes P2-0026 and J2-9168.

structure, *ACM Computing Surveys* 23(3): 345–405.

*<sup>n</sup>* values should not be very small since at

more than twice as long computation times.

should therefore be used. Similarly, the *δ*<sup>0</sup>

issues.

**8. Acknowledgments**

**9. References**

435(1893): 21–41.

70(2): 185–207.

2-88105-000-X.

27(3): 351–374.

*B* 77(14): 1–11.

14(11): 1479–1484.

boundaries was at its highest, resulting in stress redistribution from the failed areas of grain boundaries to the neighboring areas of grains. The last row therefore displays only the Mises stress for the AE+CP model. Redestribution of the Mises stress can be observed due to the degradation of the susceptible, tensile-loaded grain boundaries which decreases the amount of stress transferred from one side of the boundary to the other. Stress is redistributed in the neighboring areas of grain boundaries and grains, increasing the compressive loading. This can be seen as dark patches in the last row.

The presented approach is, however, not without its deficiencies. First, all tensile-loaded grain boundaries degrade at the same rate. In reality, higher degradation rates might be expected for the susceptible grain boundaries with higher tensile load. Also, the initial grain boundary stiffness is taken to be uniform whereas stiffness distribution based upon the properties of adjacent grains is expected (Coffman & Sethna, 2008). Lastly, the selected corrosive environment penetration approach does not account for the topology of the grain boundary network. These issues will therefore need to be addressed in future work. Slightly better convergence of the AE+CP model was observed. However, the computational times for the AE+CP model were more than twice those for the AE model, see Table 3.


Table 3. Model performance data. 10.0 *μ*m element size.
