2.2.3. Driving Potential Formulation

The driving potential (force) for inelastic deformation modes, slip and twinning, can be estimated through Equations (29) and (30), respectively.

$$\mathcal{G}^{a} = \left(1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^i\right) \left(\boldsymbol{\tau}^a + \boldsymbol{\Phi}^a - qG^\varepsilon \zeta^a \boldsymbol{\Psi}^a\right) \,. \tag{29}$$

$$\mathcal{P}^i = v^i \left( \mathfrak{r}^i + \Phi^i - \varrho G^\varepsilon \zeta^i \Psi^i \right) \,. \tag{30}$$

#### *2.3. Material Flow Modeling*

The material flow due to the activity of *α*-slip systems can be estimated through a power function of the shear strain rate. [54,56] given as:

$$\dot{\gamma}^a = \dot{\gamma}\_0 \left| \frac{\tau^a}{s\_r^a} \right|^{\frac{1}{m}} \text{sign}(\tau^a) \,. \tag{31}$$

In Equation (31), *γ*˙ <sup>0</sup> is the initial shear strain rate; *τ<sup>α</sup>* is the shear stress on the *α*-slip system; *s<sup>α</sup> <sup>r</sup>* is slip resistance; and *m* is the strain rate sensitivity parameter. In a similar manner, material flow due to mechanical twinning is modeled through a nonlinear function, as discussed in [25–27,57]. Accordingly, the rate of change of the twinned-martensite volume fraction of *i*-twin system is given as:

$$
\psi^i = \frac{\dot{\gamma}\_0}{\gamma^i} \left( \frac{\tau^i}{s\_r^i} \right)^{1/m},
\tag{32}
$$

where *γ*˙ <sup>0</sup> is the initial shear strain rate, *γ<sup>i</sup>* is the shear strain of *i*-twin system, and *s<sup>i</sup> <sup>r</sup>* is the twin resistance of *i*-twin system.

#### 2.3.1. Strain Hardening Rule

In the model presented by [53], a dislocation density-based hardening law is used to incorporate self, *s*˙ *α <sup>s</sup>* , and latent, *s*˙ *α <sup>l</sup>* , hardening contributions of *α*-slip systems as:

$$\dot{s}\_{r}^{a} = \dot{s}\_{s}^{a} + \dot{s}\_{l}^{a} = \sum\_{a=1}^{N\_{s}} h\_{0}^{a} \left[ 1 - \left( \frac{s\_{r}^{a} - s\_{r,0}^{a}}{s\_{r,S}^{a} - s\_{r,0}^{a}} \right) \right] |\dot{\gamma}^{a}| + \sum\_{\kappa=1}^{N\_{s}} h^{a\kappa} |\dot{\gamma}^{\kappa}|. \tag{33}$$

In Equation (33), *h<sup>α</sup>* <sup>0</sup> and *<sup>s</sup><sup>α</sup> <sup>r</sup>*,0 are the initial values of hardening rate and strength of slip system, respectively; *s<sup>α</sup> <sup>r</sup>*,*<sup>S</sup>* is the saturation value of slip strength; *<sup>h</sup>ακ* is an array of latent hardening values; and *γ*˙ *<sup>κ</sup>* is the shear strain rate of *κ*-slip system, where *κ* denotes a number of slip system except *α* (*κ* = *j*, *j* = 1, ... , *i* − 1, *i* + 1, ... , *Ns*). The hardening induced by mechanical twinning is incorporated in the model through Equation (34) as:

$$\dot{s}\_t^i = h\_{nc}^i \left(\sum\_{i=1}^{N\_t} \upsilon^i\right)^d \sum\_{\mu\_1=0}^i \gamma^i \dot{\upsilon}^{\mu\_1} + h\_{cp}^i \left(\sum\_{i=1}^{N\_t} \upsilon^i\right) \sum\_{\mu\_2=0}^i \gamma^i \dot{\upsilon}^{\mu\_2} \,. \tag{34}$$

Here, *h<sup>i</sup> nc* and *h<sup>i</sup> cp* are the initial hardening rates of non-coplanar and coplanar twin systems; *μ*<sup>1</sup> (*μ*<sup>1</sup> ∈ *i*) and *μ*<sup>2</sup> (*μ*<sup>2</sup> ∈ *i*) represent the number of non-coplanar and coplanar twin systems, respectively; *υ*˙*μ*<sup>1</sup> and *υ*˙*μ*<sup>2</sup> are the rate of change of volume fractions for non-coplanar and coplanar twin systems, respectively; and *d* is a material parameter.

#### 2.3.2. Microstrain Parameter

The microstrain, induced through crystal defects, in the case of slip and mechanical twinning can be estimated through Equations (35) and (36), respectively, as:

$$\zeta\_s = \frac{1}{\omega \mathcal{G}^\epsilon \mathcal{N}\_s} \sum\_{a=1}^{N\_s} \dot{s}\_r^a = \frac{1}{\omega \mathcal{G}^\epsilon \mathcal{N}\_s} \sum\_{a=1}^{N\_s} \sum\_{\kappa=1}^{N\_s} h^{a\kappa} |\dot{\gamma}^\kappa| \,. \tag{35}$$

$$\dot{\zeta}\_{t} = \frac{1}{\omega \mathbb{G}^{\varepsilon} N\_{t}} \sum\_{i=1}^{N\_{t}} \left[ \left\{ h\_{nc}^{i} \left( \sum\_{i=1}^{N\_{t}} \upsilon^{i} \right)^{d} \sum\_{\mu\_{1}=0}^{i} \upsilon^{\mu\_{1}} + h\_{cp}^{i} \left( \sum\_{i=1}^{N\_{t}} \upsilon^{i} \right) \sum\_{\mu\_{2}=0}^{i} \upsilon^{\mu\_{2}} \right\} \gamma^{i} \right] . \tag{36}$$

The stress-like parameters that are part of the driving potential equations for slip and twinning, Equations (29) and (30), are estimated as:

$$\Psi^{a} = \frac{\left(1 - \sum\_{i=1}^{N\_l} \upsilon^i\right)^{-1}}{\omega \,\mathrm{G}^c \mathrm{N}\_s} \sum\_{\kappa=1}^{N\_s} h^{a\kappa} \,. \tag{37}$$

$$\Psi^{i} = \frac{\left(\sum\_{i=1}^{N\_{\ell}} \dot{\upsilon}^{i}\right)^{-1}}{\omega \mathbb{G}^{c} N\_{\ell}} \sum\_{i=1}^{N\_{\ell}} \left[ h\_{nc}^{i} \left(\sum\_{i=1}^{N\_{\ell}} \upsilon^{i}\right)^{d} \sum\_{\mu\_{1}=0}^{i} \dot{\upsilon}^{\mu\_{1}} + h\_{cp}^{i} \left(\sum\_{i=1}^{N\_{\ell}} \upsilon^{i}\right) \sum\_{\mu\_{2}=0}^{i} \dot{\upsilon}^{\mu\_{2}}\right]. \tag{38}$$

#### **3. Numerical Integration Scheme of Constitutive Equations**

The development of the numerical integration scheme consists of the following:(i) identification of primary variables in constitutive formulation, (ii) discretization and numerical integration of equations in time domain, (iii) development of Newton–Raphson iterative scheme, and (iv) development of time sub-stepping algorithm to increase or decrease time step size, depending on the incremental values of primary variables.

#### *3.1. Identification of Primary Variables*

The constitutive equations of the slip- and twin-based crystal plasticity model are the first-order ordinary differential equations of the state variables mentioned in Equation (39), as follows:

$$\left\{\mathbf{T}^{\varepsilon},\ s\_{r}^{a},\ s\_{t'}^{i},\ \mathbf{R}^{\varepsilon},\ v^{i}\right\},\tag{39}$$

where **T***<sup>e</sup>* represents second Piola–Kirchhoff stress tensor, {*s<sup>α</sup> <sup>r</sup>* | *α* = 1, ... , *Nsl*} denotes slip resistance of *α*-slip system, {*s<sup>i</sup> <sup>t</sup>* | *i* = 1, ... , *Ntw*} shows twin resistance of *i*-twin system, **R***<sup>e</sup>* is the rigid body rotation tensor, and {*υ<sup>i</sup>* | *i* = 1, ... , *Ntw*} is twinned martensite volume fraction.

#### *3.2. Discretization of Constitutive Equations*

A time integration scheme is executed in sample coordinate axes through discretizing the deformation history in the time domain and subsequently numerically integrating constitutive equations for each time step. In order to define a deformation time history, the configurations of a material point are considered at time *tn* and *tn*+1. In this, a relation of *tn*+<sup>1</sup> = *tn* + Δ*t* is used, where *tn* and *tn*+<sup>1</sup> represent time at the start and end of the time step, respectively. Afterwards, the magnitudes of all variables are evaluated at *tn* and *tn*+1, and denoted with the subscripts *n* and *n* + 1, respectively. The numerical time integration scheme is based on the following assumptions:


The output of the numerical integration scheme provides updated values of variables as: **T***<sup>e</sup> <sup>n</sup>*+1, *<sup>s</sup><sup>α</sup> <sup>r</sup>*,*n*+1, *<sup>s</sup><sup>i</sup> <sup>t</sup>*,*n*+1, **<sup>R</sup>***<sup>e</sup> <sup>n</sup>*+1, and *<sup>υ</sup><sup>i</sup> <sup>n</sup>*+<sup>1</sup> at time *tn*+1. The constitutive equations are discretized through fully implicit time integration procedure using backward Euler scheme. According to the kinematic formulation, an incremental form of Green strain tensor in terms of the symmetric part of the velocity gradient, Equation (25) in [53], can be written as:

$$
\Delta \hat{\mathbf{E}}^c = \Delta t (\hat{\mathbf{D}} - \hat{\mathbf{D}}^\*) \,. \tag{40}
$$

Integration of Equation (40) using backward Euler scheme provides an updated value of Green strain tensor, as follows:

$$
\hat{\mathbf{E}}\_{n+1}^{\mathbf{c}} = \hat{\mathbf{E}}\_{n}^{\mathbf{c}} + \Delta \hat{\mathbf{E}}\_{n}^{\mathbf{c}} = \hat{\mathbf{E}}\_{n}^{\mathbf{c}} + \Delta t (\hat{\mathbf{D}}\_{n+1} - \hat{\mathbf{D}}\_{n+1}^{\mathbf{s}}) \,. \tag{41}
$$

where **<sup>D</sup>** *<sup>n</sup>*+<sup>1</sup> and **<sup>D</sup>** <sup>∗</sup> *<sup>n</sup>*+<sup>1</sup> are the symmetric parts of velocity gradients **<sup>L</sup>***n*+<sup>1</sup> and **<sup>L</sup>**<sup>∗</sup> *n*+1, respectively. The incremental values of these parameters can be illustrated as:

$$
\hat{\mathbf{D}}\_{n+1} = \frac{1}{2} [\check{\mathbf{L}}\_{n+1} + \hat{\mathbf{L}}\_{n+1}^{\mathrm{T}}] = \mathrm{sym}(\check{\mathbf{L}}\_{n+1})\,. \tag{42}
$$

$$
\hat{\mathbf{D}}\_{n+1}^{\*} = \text{sym}(\hat{\mathbf{C}}\_{n+1}^{\varepsilon}\hat{\mathbf{O}}\_{n+1}^{\varepsilon}) + \mathbf{R}\_{n+1}^{\varepsilon}\mathbf{D}\_{n+1}^{p}(\mathbf{R}\_{n+1}^{\varepsilon})^T,\tag{43}
$$

where an updated value of right Cauchy–Green deformation tensor in third intermediate configuration can be defined as:

$$
\hat{\mathbf{C}}\_{n+1}^{\varepsilon} = \mathbf{V}\_{n+1}^{\varepsilon} (\mathbf{V}\_{n+1}^{\varepsilon})^{\mathrm{T}} \,, \tag{44}
$$

Furthermore, an incremental magnitude of spin tensor **<sup>Θ</sup>** *<sup>e</sup> <sup>n</sup>*+<sup>1</sup> is represented as:

$$
\widetilde{\boldsymbol{\Theta}}\_{n+1}^{\varepsilon} = \dot{\mathbf{R}}\_{n+1}^{\varepsilon} (\mathbf{R}\_{n+1}^{\varepsilon})^T = \boldsymbol{\Delta t} \mathbf{R}\_{n+1}^{\varepsilon} (\mathbf{R}\_{n+1}^{\varepsilon})^T,\tag{45}
$$

In addition, **D**˘ *<sup>p</sup> <sup>n</sup>*+<sup>1</sup> is an updated symmetric part of (**C**˘ *<sup>e</sup> n*+1**L**˘ *<sup>p</sup> <sup>n</sup>*+1). Here, **<sup>C</sup>**˘ *<sup>e</sup> <sup>n</sup>*+<sup>1</sup> is an incremental form of right Cauchy–Green deformation tensor in second intermediate configuration, which can be expressed as:

$$\mathbf{\widetilde{C}}\_{n+1}^{\varepsilon} = \left(\mathbf{F}\_{n+1}^{\varepsilon}\right)^{T} \left(\mathbf{F}\_{n+1}^{\varepsilon}\right) \,. \tag{46}$$

The incremental value of the plastic velocity gradient in second intermediate configuration can be defined as:

$$\mathbf{L}\_{n+1}^{p} = \left(1 - \sum\_{i=1}^{N\_{\text{tr}}} \boldsymbol{\upsilon}\_{n+1}^{i}\right) \sum\_{a=1}^{N\_{sl}} \Delta t \boldsymbol{\dot{\gamma}\_{n+1}^{a}} \mathbf{S}\_{n+1}^{a} \tag{47}$$

$$+ \mathbf{F}\_{s,n+1}^{p} \left\{\sum\_{i=1}^{N\_{\text{tr}}} \boldsymbol{\upsilon}\_{n+1}^{i} \Delta t \boldsymbol{\dot{\gamma}\_{n+1}^{i}} \mathbf{S}\_{n+1}^{i}\right\} \left(\mathbf{F}\_{s,n+1}^{p}\right)^{-1} \ .$$

where *υ<sup>i</sup> <sup>n</sup>*+<sup>1</sup> represents an incremental change in the volume fraction of *i th* twin system, *γ*˙ *α <sup>n</sup>*+<sup>1</sup> is an updated value of the shear strain rate of *<sup>α</sup>*-slip system, *<sup>γ</sup>*˙ *<sup>i</sup> <sup>n</sup>*+<sup>1</sup> is an updated value of shear strain induced by *i*-twin system, **S**¯ *<sup>α</sup> <sup>n</sup>*+<sup>1</sup> and **<sup>S</sup>**ˇ*<sup>i</sup> <sup>n</sup>*+<sup>1</sup> are the updated magnitudes of Schmid and twin orientation tensors in second intermediate configuration, and **F***<sup>p</sup> <sup>s</sup>*,*n*+<sup>1</sup> is an incremental change in the value of the plastic deformation gradient. Substitution of **D**˘ *<sup>p</sup> n*+1 in Equation (43) gives the following:

$$\begin{split} \tilde{\mathbf{D}}\_{n+1}^{\epsilon} &= \text{sym}(\tilde{\mathbf{C}}\_{n+1}^{\epsilon} \tilde{\mathbf{G}}\_{n+1}^{\epsilon}) \\ &+ (\mathbf{R}\_{n+1}^{\epsilon}) \text{sym} \left[ \check{\mathbf{C}}\_{n+1}^{\epsilon} \left\{ \left( 1 - \sum\_{i=1}^{N\_{ltr}} \upsilon\_{n+1}^{i} \right) \sum\_{a=1}^{N\_{l}} \Delta t \boldsymbol{\gamma}\_{n+1}^{a} \bar{\mathbf{S}}^{a} \right. \\ &+ \mathbf{F}\_{s,n+1}^{p} \left\{ \sum\_{i=1}^{N\_{ltr}} \upsilon\_{n+1}^{i} \Delta t \boldsymbol{\gamma}\_{n+1}^{i} \check{\mathbf{S}}^{i} \right\} \left( \mathbf{F}\_{s,n+1}^{p} \right)^{-1} \right] \left( \mathbf{R}\_{n+1}^{\epsilon} \right)^{T} . \end{split} \tag{48}$$

Considering the effects of rigid body rotation tensor **R***<sup>e</sup> <sup>n</sup>*+<sup>1</sup> on Schmid and twin orientation tensors, Equation (48) can also be written as:

$$\begin{split} \tilde{\mathbf{D}}\_{n+1}^{\star} &= \text{sym} \left( \tilde{\mathbf{C}}\_{n+1}^{\epsilon} \tilde{\mathbf{B}}\_{n+1}^{\epsilon} \right) \\ &+ \text{sym} \left[ \mathbf{C}\_{n+1}^{\epsilon} \left\{ \left( 1 - \sum\_{i=1}^{N\_{\text{ltr}}} \boldsymbol{\nu}\_{n+1}^{i} \right) \sum\_{a=1}^{N\_{\text{ltr}}} \Delta t \boldsymbol{\gamma}\_{n+1}^{a} \left( \mathbf{R}\_{n+1}^{\epsilon} \right) \left( \mathbf{m}^{a} \otimes \mathbf{\dot{n}}^{a} \right) \left( \mathbf{R}\_{n+1}^{\epsilon} \right)^{T} \right. \\ &+ \mathbf{F}\_{s,n+1}^{p} \left\{ \sum\_{i=1}^{N\_{\text{ltr}}} \boldsymbol{\nu}\_{n+1}^{i} \Delta t \boldsymbol{\gamma}\_{n+1}^{i} \left( \mathbf{R}\_{n+1}^{\epsilon} \right) \left( \mathbf{m}^{i} \otimes \mathbf{\dot{n}}^{i} \right) \left( \mathbf{R}\_{n+1}^{\epsilon} \right)^{T} \right\} \left( \mathbf{F}\_{s,n+1}^{p} \right)^{-1} \right\} \,. \end{split} \tag{49}$$

The forward mapping approach for Schmid and twin orientation tensors can be implemented into Equation (49) as:

$$\begin{split} \left\{ \left( \mathbf{R}\_{n+1}^{\varepsilon} . (\mathbf{m}^{a} \otimes \mathbf{n}^{a}) \right) . (\mathbf{R}\_{n+1}^{\varepsilon})^{T} = \left\{ \mathbf{R}\_{n+1}^{\varepsilon} . (\mathbf{Q}\_{0} . \mathbf{m}\_{0}^{a} \otimes \mathbf{Q}\_{0} . \mathbf{n}\_{0}^{a}) \right\} . (\mathbf{R}\_{n+1}^{\varepsilon})^{T} \\ = \left\{ \mathbf{R}\_{n+1}^{\varepsilon} . (\mathbf{Q}\_{0} . \mathbf{m}\_{0}^{a}) \right\} \otimes \left\{ \mathbf{R}\_{n+1}^{\varepsilon} . (\mathbf{Q}\_{0} . \mathbf{n}\_{0}^{a}) \right\} . \end{split} \tag{50}$$

where **Q**<sup>0</sup> is the initial rotation matrix that is used to transform crystal coordinates to sample coordinate systems through Euler angles, and **m***<sup>α</sup>* <sup>0</sup> and **<sup>n</sup>***<sup>α</sup>* <sup>0</sup> are the initial vectors representing *α*-slip system in reference configuration. It was stated previously that the rotation matrix (Euler angles) can be updated through a rigid body rotation tensor as follows: **Q***n*+<sup>1</sup> = **R***e <sup>n</sup>*+1.**Q**0. The rotation matrix can also be used to transform Schmid orientation vectors from reference to first intermediate configuration as follows: **m**¯ *<sup>α</sup>* = **Q**0.**m***<sup>α</sup>* <sup>0</sup> and **<sup>n</sup>**¯ *<sup>α</sup>* = **<sup>Q</sup>**0.**n***<sup>α</sup>* 0. In this condition, Equation (50) can be rewritten as:

$$\left\{ \left( \mathbf{R}\_{n+1}^{\varepsilon} \right) . (\bar{\mathbf{m}}^a \otimes \bar{\mathbf{n}}^a) \right\} \left( \mathbf{R}\_{n+1}^{\varepsilon} \right)^T = \left( \mathbf{Q}\_{n+1} . \mathbf{m}\_0^a \right) \otimes \left( \mathbf{Q}\_{n+1} . \mathbf{n}\_0^a \right) = \widetilde{\mathbf{m}}\_{n+1}^a \otimes \widetilde{\mathbf{n}}\_{n+1}^a \,, \tag{51}$$

In Equation (51), **<sup>m</sup>** *<sup>α</sup> <sup>n</sup>*+<sup>1</sup> and **<sup>n</sup>***<sup>α</sup> <sup>n</sup>*+<sup>1</sup> are the slip incremental values of direction and area normal vectors of *α*-slip system in third intermediate configuration. It is clear from Equation (51) that the updated rotation matrix **Q***n*+<sup>1</sup> is used to transform Schmid orientation vectors from reference to third intermediate configuration. Similarly, twin orientation vectors are given as:

$$\{\mathbf{R}\_{n+1}^{\varepsilon}.(\mathbf{n}^{\hat{i}} \odot \hat{\mathbf{n}}^{\hat{i}})\}.(\mathbf{R}\_{n+1}^{\varepsilon})^T = \tilde{\mathbf{m}}\_{n+1}^{\hat{i}} \otimes \tilde{\mathbf{n}}\_{n+1}^{\hat{i}}.\tag{52}$$

Substitution of Equations (51) and (52) in (49) provides

$$\begin{split} \widetilde{\mathbf{D}}\_{n+1}^{\*} &= \text{sym} \left( \widetilde{\mathbf{C}}\_{n+1}^{c} \widetilde{\mathbf{O}}\_{n+1}^{c} \right) \\ &+ \text{sym} \left[ \widetilde{\mathbf{C}}\_{n+1}^{c} \cdot \left\{ \left( 1 - \sum\_{i=1}^{N\_{lw}} \upsilon\_{n+1}^{i} \right) \sum\_{a=1}^{N\_{al}} \Delta t \dot{\gamma}\_{n+1}^{a} \widetilde{\mathbf{S}}\_{n+1}^{a} \right. \\ &+ \mathbf{F}\_{s,n+1}^{p} \left\{ \sum\_{i=1}^{N\_{lw}} \upsilon\_{n+1}^{i} \Delta t \dot{\gamma}\_{n+1}^{i} \widetilde{\mathbf{S}}\_{n+1}^{i} \right\} (\mathbf{F}\_{s,n+1}^{p})^{-1} \right] . \end{split} \tag{53}$$

Here, **<sup>S</sup>***<sup>α</sup> <sup>n</sup>*+<sup>1</sup> <sup>=</sup> **<sup>m</sup>** *<sup>α</sup> <sup>n</sup>*+<sup>1</sup> <sup>⊗</sup> **<sup>n</sup>***<sup>α</sup> <sup>n</sup>*+<sup>1</sup> and **<sup>S</sup>***<sup>i</sup> <sup>n</sup>*+<sup>1</sup> <sup>=</sup> **<sup>m</sup>** *<sup>i</sup> <sup>n</sup>*+<sup>1</sup> <sup>⊗</sup> **<sup>n</sup>***<sup>i</sup> <sup>n</sup>*+<sup>1</sup> are the updated Schmid and twin orientation tensors in third intermediate configuration. The updated values of shear strain for *α*-slip system and volume fraction of twinned martensite at time *tn*+<sup>1</sup> can be defined by Equations (54) and (55), respectively, as:

$$
\gamma\_{n+1}^{\mathfrak{a}} = \gamma\_n^{\mathfrak{a}} + \Delta \gamma^{\mathfrak{a}} = \gamma\_n^{\mathfrak{a}} + \gamma\_0 \left| \frac{\mathfrak{T}\_{n+1}^{\mathfrak{a}}}{s\_{r,n+1}^{\mathfrak{a}}} \right|^{\frac{1}{m}} \sin(\mathfrak{r}\_{n+1}^{\mathfrak{a}}) \,. \tag{54}
$$

$$\boldsymbol{\upsilon}\_{n+1}^{i} = \boldsymbol{\upsilon}\_{n}^{i} + \Delta \boldsymbol{\upsilon}^{i} = \boldsymbol{\upsilon}\_{n}^{i} + \frac{\gamma\_{0}}{\gamma\_{n+1}^{i}} \left| \frac{\boldsymbol{\tau}\_{n+1}^{i}}{\boldsymbol{s}\_{t,n+1}^{i}} \right|^{\frac{1}{m}}.\tag{55}$$

As discussed previously, the rigid body rotation tensor **R***<sup>e</sup> <sup>n</sup>*+<sup>1</sup> updates the crystal orientation (Euler angles) matrix **Q***n*+1, which can be used to transform the fourth order elasticity tensors for slip <sup>C</sup>*<sup>s</sup> <sup>n</sup>*+<sup>1</sup> and twinned <sup>C</sup>*<sup>t</sup> <sup>n</sup>*+<sup>1</sup> regions to the third intermediate configuration as follows:

$$
\check{\mathbf{C}}\_{n+1}^s = \left(\mathbf{Q}\_{n+1} \otimes \mathbf{Q}\_{n+1}\right) : \check{\mathbb{C}}\_0^s : \left(\mathbf{Q}\_{n+1} \otimes \mathbf{Q}\_{n+1}\right)^T. \tag{56}
$$

$$
\widetilde{\mathbb{C}}\_{n+1}^{t} = \left(\mathbf{Q}\_{n+1} \otimes \mathbf{Q}\_{n+1}\right) : \widetilde{\mathbb{C}}\_{0}^{t} : \left(\mathbf{Q}\_{n+1} \otimes \mathbf{Q}\_{n+1}\right)^{T} \,. \tag{57}
$$

Any one of the elasticity tensors can be obtained in terms of another by using coordinate transformation rule defined as:

$$
\widetilde{\mathbb{C}}\_{ijkl}^t = \widetilde{\mathbb{C}}\_{pqrs}^s \mathbf{P}\_{ip} \mathbf{P}\_{jq} \mathbf{P}\_{kr} \mathbf{P}\_{ls} = (\mathbf{P} \otimes \mathbf{P}) : \widetilde{\mathbb{C}}^s : (\mathbf{P} \otimes \mathbf{P})^T. \tag{58}
$$

Here, [**P**] is the transformation matrix that relates lattice orientations in twinned and untwinned (slipped) regions. The components of transformation matrix [**P**] can be defined through Equation (59), given by [26], as:

$$\mathbf{P}\_{\text{jk}} = 2n\_{\text{j}}n\_{k} - \delta\_{\text{jk}} \quad , \quad j,k = 1,2,3 \tag{59}$$

where **n** is the area normal vector of the twin plane, and *δ* is the Kroneker delta. An equivalent elasticity tensor can be calculated as:

$$
\breve{\mathbb{C}}\_{n+1}^{\mathfrak{c}} = \left(1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}\_{n+1}^{i}\right) \breve{\mathbb{C}}\_{n+1}^{\mathfrak{s}} + \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}\_{n+1}^{i} \breve{\mathbb{C}}\_{n+1}^{\mathfrak{t}} \,. \tag{60}
$$

Furthermore, an updated form of equivalent second Piola–Kirchhoff stress tensor can be estimated as:

$$\mathbf{T}\_{n+1}^{\mathbf{c}} = \hat{\mathbb{C}}\_{n+1}^{\mathbf{c}} : \hat{\mathbb{E}}\_{n+1}^{\mathbf{c}} \ . \tag{61}$$

An evolution equation for rigid body rotation tensor **R***<sup>e</sup>* is integrated through the exponential map discussed by [58] as follows:

$$\mathbf{R}\_{n+1}^{\varepsilon} = \exp(\Delta t \tilde{\Theta}\_{n+1}^{\varepsilon}) . \mathbf{R}\_{n}^{\varepsilon} \; . \tag{62}$$

where an updated value of the spin of lattice **<sup>Θ</sup>** *<sup>e</sup> <sup>n</sup>*+<sup>1</sup> can be estimated through the skewsymmetric component -**W**<sup>∗</sup> *<sup>n</sup>*+<sup>1</sup> of the velocity gradient **<sup>L</sup>**<sup>∗</sup> *<sup>n</sup>*+<sup>1</sup> as:

$$
\bar{\mathbf{W}}\_{n+1}^{\*} = \text{skew}(\hat{\mathbf{C}}\_{n+1}^{\epsilon} \hat{\mathbf{O}}\_{n+1}^{\epsilon}) + \sum\_{a=1}^{N\_s} \Delta t \dot{\gamma}\_{n+1}^{a} \text{skew}(\hat{\mathbf{C}}\_{n+1}^{\epsilon} \hat{\mathbf{S}}\_{n+1}^{a}) \,. \tag{63}
$$

By using backward mapping, the skew-symmetric component of velocity gradient **L***n*+<sup>1</sup> can be estimated as:

$$\widehat{\boldsymbol{\mathbf{W}}}\_{n+1} = \left\{ (\mathbf{V}\_{n+1}^{\varepsilon})^T \mathbf{W}\_{n+1} \right\} \mathbf{V}\_{n+1}^{\varepsilon} = (\mathbf{V}^{\varepsilon})^T \text{skew} \{ (\mathbf{V}\_{n+1}^{\varepsilon})^T \{ \boldsymbol{\Delta t} \mathbf{V}\_{n+1}^{\varepsilon} \} \right\} \mathbf{V}^{\varepsilon} + \overleftarrow{\mathbf{W}}\_{n+1}^{\*} \,. \tag{64}$$

Here, **<sup>W</sup>***n*+<sup>1</sup> is the updated skew-symmetric component of **<sup>L</sup>***n*+1. Substitution of -**W**<sup>∗</sup> *n*+1 from Equation (63) in (64) provides

$$\begin{split} \text{skew}(\widetilde{\mathbf{C}}\_{n+1}^{\varepsilon}\widetilde{\mathbf{O}}\_{n+1}^{\varepsilon}) &= \widetilde{\mathbf{W}}\_{n+1} - \left(\mathbf{V}^{\varepsilon}\right)^{T} \text{skew}\{\left(\mathbf{V}\_{n+1}^{\varepsilon}\right)^{T}\Delta t \mathbf{V}\_{n+1}^{\varepsilon}\} \mathbf{V}^{\varepsilon} \\ &- \sum\_{a=1}^{N\_{s}} \Delta t \dot{\gamma}\_{n+1}^{a} \text{skew}(\widetilde{\mathbf{C}}\_{n+1}^{\varepsilon}\widetilde{\mathbf{S}}\_{n+1}^{a}) \ . \end{split} \tag{65}$$

For small elastic strain problems, the value of the right Cauchy–Green deformation tensor **<sup>C</sup>***<sup>e</sup> <sup>n</sup>*+<sup>1</sup> is typically small. In this case, Equation (65) can be written as:

$$\text{skew}(\check{\Theta}\_{n+1}^{\varepsilon}) = \overline{\mathbf{W}}\_{n+1} - (\mathbf{V}^{\varepsilon})^T \text{skew}\{ (\mathbf{V}\_{n+1}^{\varepsilon})^T \Delta t \mathbf{V}\_{n+1}^{\varepsilon} \} \mathbf{V}^{\varepsilon} - \sum\_{a=1}^{N\_b} \Delta t \dot{\gamma}\_{n+1}^a \text{skew}(\check{\mathbf{S}}\_{n+1}^a) \tag{66}$$

An incremental value of the slip resistance of *α*-slip system can be evaluated using a backward Euler scheme as follows:

$$s\_{r,n+1}^a = s\_{r,n}^a + \Delta s\_r^a = s\_{r,n}^a + \sum\_{a=1}^{N\_s} h\_0^a \left[ 1 - \left( \frac{s\_{r,n+1}^a - s\_{r,0}^a}{s\_{r,S\_{n+1}}^a - s\_{r,0}^a} \right) \right] \left| (\Delta t) \gamma\_{n+1}^a \right| \,. \tag{67}$$

In the present model, slip resistances for all *α*-slip systems are considered to be similar. Therefore, Equation (67) can be modified as follows:

$$s\_{r,n+1} = s\_{r,n} + h\_0 \left[ 1 - \left( \frac{s\_{r,n+1} - s\_{r,0}}{s\_{r,S\_{n+1}} - s\_{r,0}} \right) \right] \sum\_{a=1}^{N\_s} \left| (\Delta t) \dot{\gamma}\_{n+1}^a \right| \,\tag{68}$$

In Equation (68), *sr*,*Sn*+<sup>1</sup> is an incremental saturation value of slip resistance, which can be calculated as:

$$s\_{r, \mathcal{S}\_{n+1}} = s\_{r, \mathcal{S}\_n} + \Delta s\_{r, \mathcal{S}} = s\_{r, \mathcal{S}\_n} + s\_{r, \mathcal{S}\_0} \left[ \frac{\sum\_{\mathcal{S}}^{\mathcal{N}\_\mathbf{s}} |\gamma\_{n+1}^{\mathcal{S}}|}{\gamma\_{\mathcal{S}\_0}} \right]^{\frac{k\theta}{(\mathcal{G}\_{\mathbf{n}+1} + \mathbf{p}^{\mathcal{S}\_Z})}} + s\_{r, p} \left( \sum\_{\lambda=0}^i \upsilon\_{n+1}^{\lambda} \right)^{a\_1}, \tag{69}$$

where *sr*,*S*<sup>0</sup> is the initial value of saturation slip resistance, Δ*γS*<sup>0</sup> is an incremental initial value of the slip system shear strain at the initial value of saturation slip resistance, *a* is the material parameter, *sr*,*<sup>p</sup>* is the material slip-hardening parameter, *λ* (*λ* ∈ *i*) represents the number of non-coplanar twin systems to slip system, and *a*<sup>1</sup> is a material parameter. In addition, the twin resistance can be estimated using Equation (34) as:

$$\begin{split} s\_{t,n+1}^{\dot{i}} = s\_{t,n}^{\dot{i}} + \Delta s\_t^{\dot{i}} = s\_{t,n}^{\dot{i}} + \left[ h\_{nc}^{\dot{i}} \left( \sum\_{i=1}^{N\_l} \upsilon\_{n+1}^{\dot{i}} \right)^d \sum\_{\mu\_1=0}^{\dot{i}} \gamma\_{n+1}^{\dot{i}} \Delta t \dot{\upsilon}\_{n+1}^{\mu\_1} \\ + h\_{cp}^{\dot{i}} \left( \sum\_{i=1}^{N\_l} \upsilon\_{n+1}^{\dot{i}} \right) \sum\_{\mu\_2=0}^{\dot{i}} \gamma\_{n+1}^{\dot{i}} \Delta t \dot{\upsilon}\_{n+1}^{\mu\_2} \right] \end{split} \tag{70}$$

The resistance of all twin systems are assumed to be identical. Therefore, Equation (70) can be expressed as:

$$\begin{split} s\_{t, \mathbf{u} + 1} = s\_{t, \mathbf{u}} + \left[ h\_{\mathbf{nc}} \left( \sum\_{i=1}^{N\_t} \boldsymbol{\upsilon}\_{n+1}^i \right)^d \sum\_{\mu\_1 = 0}^i \boldsymbol{\upsilon}\_{n+1}^i \Delta t \boldsymbol{\upsilon}\_{n+1}^{\mu\_1} \\ + h\_{\mathbf{cp}} \left( \sum\_{i=1}^{N\_t} \boldsymbol{\upsilon}\_{n+1}^i \right) \sum\_{\mu\_2 = 0}^i \boldsymbol{\upsilon}\_{n+1}^i \Delta t \boldsymbol{\upsilon}\_{n+1}^{\mu\_2} \right]. \end{split} \tag{71}$$

The updated driving potentials for slip and twin mode of deformations can be estimated using Equations (29) and (30), respectively as:

$$\mathbf{G}\_{\mathfrak{s},n+1}^{\mathfrak{a}} = \mathbf{G}\_{\mathfrak{s},\mathfrak{a}}^{\mathfrak{a}} + \boldsymbol{\Delta}\mathbf{G}\_{\mathfrak{s}}^{\mathfrak{a}} = \mathbf{G}\_{\mathfrak{s},\mathfrak{a}}^{\mathfrak{a}} + \left(1 - \sum\_{i=1}^{N\_{\mathfrak{l}}} \boldsymbol{\upsilon}\_{n+1}^{i}\right) \left(\boldsymbol{\pi}\_{n+1}^{\mathfrak{a}} + \boldsymbol{\Phi}^{\mathfrak{a}} - q\boldsymbol{G}\_{n+1}^{\mathfrak{c}}\boldsymbol{\zeta}\_{n+1}\boldsymbol{\Psi}\_{n+1}^{\mathfrak{a}}\right) . \tag{72}$$

$$\mathbf{G}\_{t,n+1}^{i} = \mathbf{G}\_{t,n}^{i} + \Delta \mathbf{G}\_{t}^{i} = \mathbf{G}\_{t,n}^{i} + \boldsymbol{\nu}\_{n+1}^{i} \left(\boldsymbol{\tau}\_{n+1}^{i} + \Phi^{i} - \boldsymbol{\varphi} \mathbf{G}\_{n+1}^{\boldsymbol{\varepsilon}} \boldsymbol{\zeta}\_{n+1} \boldsymbol{\Psi}\_{n+1}^{i}\right). \tag{73}$$

The parameters Φ*<sup>α</sup>* and Φ*<sup>i</sup>* are the thermal analogues for *α*-slip and *i*-twin systems, respectively. For an isothermal process, these factors will remain constant during the deformation process. *ϕ* is a dimensionless dislocation parameter, and it is assumed to be constant throughout the deformation. An equivalent modulus of rigidity *G<sup>e</sup> <sup>n</sup>*+<sup>1</sup> can be estimated as:

$$\mathbf{G}\_{n+1}^{\varepsilon} = \left(1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}\_{n+1}^{i}\right) \mathbf{G}\_{\mathbf{s}} + \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}\_{n+1}^{i} \mathbf{G}\_{\mathbf{t}} \,. \tag{74}$$

The incremental forms of crystal defect microstrain parameters for the slip and twin modes of deformations can be calculated using Equations (35) and (36), respectively, as:

$$\mathcal{L}\_{s,n+1} = \mathcal{L}\_{s,n} + \Delta \mathcal{L}\_s = \mathcal{L}\_{s,n} + \frac{1}{\omega \sigma \mathcal{G}\_{n+1}^{\varepsilon} N\_s} \sum\_{a=1}^{N\_s} \sum\_{\beta=1}^{N\_s} h\_{n+1}^{a\beta} |\gamma\_{n+1}^{\beta}| \,. \tag{75}$$

$$\begin{split} \zeta\_{t,n+1} &= \zeta\_{t,n} + \Delta \zeta\_t = \zeta\_{t,n} + \frac{1}{\omega \circ G\_{n+1}^{\epsilon} \, N\_l} \Bigg[ \Bigg( h\_{n\epsilon} \left( \sum\_{i=1}^{N\_l} \upsilon\_{n+1}^i \right)^d \sum\_{\mu\_1=0}^i (\Delta t) \upsilon\_{n+1}^{\mu\_1} \\ &+ h\_{\epsilon j} \left( \sum\_{i=1}^{N\_l} \upsilon\_{n+1}^i \right) \sum\_{\mu\_2=0}^i (\Delta t) \upsilon\_{n+1}^{\mu\_2} \Bigg) \sum\_{i=1}^{N\_l} \gamma\_{n+1}^i \Bigg] . \end{split} \tag{76}$$

The updated stress-like functions Ψ*<sup>α</sup> <sup>n</sup>*+<sup>1</sup> and <sup>Ψ</sup>*<sup>i</sup> <sup>n</sup>*+<sup>1</sup> for slip and twin modes of deformation can be calculated using Equations (37) and (38), respectively, as:

$$\Psi\_{n+1}^{a} = \frac{\left(1 - \sum\_{i=1}^{N\_l} \nu\_{n+1}^i\right)^{-1}}{\omega \,\mathrm{G}\_{n+1}^{\varepsilon} \,\mathrm{N}\_s} \sum\_{\beta=1}^{N\_b} h\_{n+1}^{a\beta} \,. \tag{77}$$

$$\begin{split} \Psi\_{n+1}^{i} = \frac{\left(\sum\_{i=1}^{N\_{t}} \boldsymbol{\upsilon}\_{n+1}^{i}\right)^{-1}}{\omega \mathbf{G}\_{n+1}^{c} N\_{t}} \left[ h\_{\text{nc}} \left( \sum\_{i=1}^{N\_{t}} \boldsymbol{\upsilon}\_{n+1}^{i} \right)^{d} \sum\_{\mu\_{1}=0}^{i} \boldsymbol{\upsilon}\_{n+1}^{\mu\_{1}} \\ &+ h\_{\text{cp}} \left( \sum\_{i=1}^{N\_{t}} \boldsymbol{\upsilon}\_{n+1}^{i} \right) \sum\_{\mu\_{2}=0}^{i} \boldsymbol{\upsilon}\_{n+1}^{\mu\_{2}} \right]. \end{split} \tag{78}$$

Furthermore, the updated Cauchy stress tensor **T***n*+<sup>1</sup> can be calculated using second Piola–Kirchhoff stress tensor, **T***<sup>e</sup> <sup>n</sup>*+1, as:

$$\mathbf{T}\_{n+1} = \mathbf{F}\_{n+1}^{\mathbf{c}} \left\{ \det \mathbf{F}\_{n+1}^{\mathbf{c}} . \mathbf{T}\_{n+1}^{\mathbf{c}} \right\} \left( \mathbf{F}\_{n+1}^{\mathbf{c}} \right)^{T} . \tag{79}$$

#### *3.3. Newton–Raphson Iterative Scheme*

In the preceding section, a set of couple nonlinear algebraic Equations (61), (62), (68), (71) and (55) for the primary variables **T***<sup>e</sup> <sup>n</sup>*+1, **<sup>R</sup>***<sup>e</sup> <sup>n</sup>*+1, *sr*,*n*<sup>+</sup>1, *st*,*n*<sup>+</sup>1, and *<sup>υ</sup><sup>i</sup> <sup>n</sup>*+1, respectively, were developed. In this, the updated form of these primary variables are calculated. A set of primary variables can be represented by a vector {*p<sup>v</sup> <sup>i</sup>* |*i* = 1, 2, . . . , 5}, where

$$p\_1^\upsilon = \Delta \mathbf{T}^\varepsilon \,, \ p\_2^\upsilon = \Delta \mathbf{s}\_r \,, \ p\_3^\upsilon = \Delta v^i \,, \ p\_4^\upsilon = \Delta \mathbf{R}^\varepsilon \,, \ p\_5^\upsilon = \Delta \mathbf{s}\_t \, \,. \tag{80}$$

The primary variables are the main constituents (directly or indirectly) of the constitutive model. An elastic-plastic response of a system mainly governs by these variables. In order to obtain an overall response of a material, a set of five equations in terms of primary variables are constructed in a residual format. A residue of the system of equations can also be represented in a vector form as {R*<sup>i</sup>* |*i* = 1, 2, ... , 5}. The components of R*<sup>i</sup>* are given through Equations (81)–(85) as:

$$\begin{split} \mathcal{R}\_1 &= \hat{\mathcal{R}}\_1 (\mathbf{T}\_{n+1}^{\varepsilon}, s\_{r,n+1}^{\mathfrak{a}}, \boldsymbol{\upsilon}\_{n+1}^{\mathfrak{l}}, \mathbf{R}\_{n+1}^{\varepsilon}, \mathbf{s}\_{t,n+1}^{\mathfrak{l}}) \\ &= (\breve{\mathbf{C}}\_{n+1}^{\varepsilon})^{-1} : \mathbf{T}\_{n+1}^{\varepsilon} - \breve{\mathbf{E}}\_{n+1}^{\varepsilon} \\ &= (\breve{\mathbf{C}}\_{n+1}^{\varepsilon})^{-1} : \mathbf{T}\_{n+1}^{\varepsilon} - \breve{\mathbf{E}}\_{n}^{\varepsilon} - \Delta t (\breve{\mathbf{D}}\_{n+1} - \breve{\mathbf{D}}\_{n+1}^{\ast}) = \boldsymbol{0} \ . \end{split} \tag{81}$$

$$\begin{split} \mathcal{R}\_2 &= \mathring{\mathcal{R}}\_2(\mathbf{T}\_{n+1}^{\epsilon}, s\_{r,n+1}^a, \boldsymbol{v}\_{n+1}^i, \mathbf{R}\_{n+1}^{\epsilon}, s\_{t,n+1}^i) \\ &= s\_{r,n+1}^a - s\_{r,n}^a - \sum\_{a=1}^{N\_a} h\_0^a \left[ 1 - \left( \frac{s\_{r,n+1}^a - s\_{r,0}^a}{s\_{r,S\_{n+1}}^a - s\_{r,0}^a} \right) \right] |(\boldsymbol{\Delta t}) \dot{\boldsymbol{\gamma}}\_{n+1}^a| = 0 \ . \end{split} \tag{82}$$

$$\begin{split} \mathcal{R}\_3 &= \mathring{\mathcal{R}}\_3(\mathbf{T}\_{n+1}^{\epsilon}, s\_{r,n+1}^{\bar{a}}, \boldsymbol{\upsilon}\_{n+1}^{\bar{i}}, \mathbf{R}\_{n+1}^{\epsilon}, s\_{t,n+1}^{\bar{i}}) \\ &= \boldsymbol{\upsilon}\_{n+1}^{\bar{i}} - \boldsymbol{\upsilon}\_{n}^{\bar{i}} - \frac{\gamma\_0}{\gamma\_{n+1}^{\bar{i}}} \left| \frac{\boldsymbol{\tau}\_{n+1}^{\bar{i}}}{s\_{t,n+1}^{\bar{i}}} \right|^{1/m} = 0 \ . \end{split} \tag{83}$$

$$\begin{split} \mathcal{R}\_4 &= \mathcal{R}\_4(\mathbf{T}\_{n+1'}^\varepsilon s\_{r,n+1'}^\varepsilon \upsilon\_{n+1'}^\varepsilon \mathbf{R}\_{n+1'}^\varepsilon s\_{t,n+1}^\varepsilon) \\ &= \mathbf{R}\_{n+1}^\varepsilon - \exp\left\{\Delta t \check{\Theta}\_{n+1}^\varepsilon\right\} . \mathbf{R}\_n^\varepsilon = 0 \ . \end{split} \tag{84}$$

$$\begin{split} \mathcal{R}\_{5} &= \mathring{\mathcal{R}}\_{5} (\mathbf{T}\_{n+1}^{\varepsilon}, s\_{r,n+1}^{i}, \boldsymbol{\upsilon}\_{n+1}^{i}, \mathbf{R}\_{n+1}^{\varepsilon}, s\_{t,n+1}^{i}) \\ &= s\_{t,n+1}^{i} - s\_{t,n}^{i} - \left[ h\_{\mathrm{nc}}^{i} \left( \sum\_{i=1}^{N\_{t}} \boldsymbol{\upsilon}\_{n+1}^{i} \right)^{d} \sum\_{\mu\_{1}=0}^{i} \boldsymbol{\gamma}\_{n+1}^{i} (\Delta t) \boldsymbol{\dot{\upsilon}}\_{n+1}^{\mu\_{1}} \\ &\quad + h\_{\mathrm{cp}}^{i} \left( \sum\_{i=1}^{N\_{t}} \boldsymbol{\upsilon}\_{n+1}^{i} \right) \sum\_{\mu\_{2}=0}^{i} \boldsymbol{\gamma}\_{n+1}^{i} (\Delta t) \boldsymbol{\dot{\upsilon}}\_{n+1}^{\mu\_{2}} \right] = 0 \end{split} \tag{85}$$

In the next step, Equations (81) and (82) are solved using a full Newton–Raphson (N-R) method, since these two are implicit in nature; however, the rest are explicit. In the current work, a two-level iterative scheme, similar as that presented by [56], is used to obtain the values of primary variables. In the first level of iteration, the N-R method is used to solve Equation (81) for **T***<sup>e</sup> <sup>n</sup>*+<sup>1</sup> by assuming the best possible values of the other primary variables (**R***<sup>e</sup> n*+1,*s<sup>α</sup> r*,*n*+1,*s<sup>i</sup> <sup>t</sup>*,*n*+1, *<sup>υ</sup><sup>i</sup> <sup>n</sup>*+1). Once the updated value of the second Piola–Kirchhoff stress tensor is obtained, the second level of the iterative procedure is performed, which includes an N-R solution of the slip resistance *s<sup>α</sup> <sup>r</sup>*,*n*+<sup>1</sup> from Equation (82). This considers an updated value of **T***<sup>e</sup> <sup>n</sup>*+<sup>1</sup> and the estimated values of **<sup>R</sup>***<sup>e</sup> <sup>n</sup>*+1, *<sup>s</sup><sup>i</sup> <sup>t</sup>*,*n*+1, and *<sup>υ</sup><sup>i</sup> <sup>n</sup>*+1. Finally, updated values of the twinned martensite volume fraction *υ<sup>i</sup> <sup>n</sup>*+1, rigid body rotation tensor **<sup>R</sup>***<sup>e</sup> <sup>n</sup>*+1, and twin resistance *s<sup>i</sup> <sup>t</sup>*,*n*+<sup>1</sup> are calculated from Equations (83)–(85), respectively.

#### Convergence Criterion

The convergence criterion is required to terminate the iterative loop once the solution is assumed to be sufficiently accurate. The convergence criterion for the presented twolevel iterative procedure is based on the variation of L2-norm for **T***<sup>e</sup> <sup>n</sup>*+1, and *<sup>s</sup><sup>α</sup> <sup>r</sup>*,*n*+1. If the L2-norm of the residuals is less than an imposed tolerance, then the incremental solution of the updated primary variables is converged (fully elastic). Otherwise, the trial value must be updated iteratively (based on the calculated value) until the residual satisfies the convergence criterion, given as:

$$\|\|\mathcal{R}\_{\text{trial}}\|\|\_{2} < \text{Tol.}\tag{86}$$

In the current Newton–Raphson iterative scheme, iterations are carried out unless and until the variation in the L2-norm of the residuals for **T***<sup>e</sup> <sup>n</sup>*+<sup>1</sup> and *<sup>s</sup><sup>α</sup> <sup>r</sup>*,*n*+<sup>1</sup> satisfy the conditions in Equations (87) and (88), respectively, as follows:

$$\|\Delta \mathbf{T}^c\|\_2 < (10^{-4}) \|\mathbf{T}^c\|\_{\text{trial}} \tag{87}$$

$$\|\Delta s\_r^{\alpha}\|\_2 < (10^{-4})|s\_r^{\alpha}|\_{\text{trial}}\tag{88}$$

In the current work, a time sub-stepping algorithm (TSSA) is used in the numerical integration scheme to control the time step size. The advantage of using time sub-stepping in an iterative procedure is to improve the convergence of a solution by reducing the time increment once needed. The TSSA must also be capable of increasing the time step size at a material point where convergence can easily be achieved to reduce computational time. The TSSA can reduce and increase the time step size based on the incremental variation. In this algorithm, if any of the convergence conditions, Equations (87) and (88), are not

satisfied, then the N-R iterative scheme will call TSSA, which controls the time step size according to the incremental variation in the values of primary variables.

#### *3.4. Time Sub-Stepping Algorithm*

The proposed TSSA is based on the ratio of the maximum and desired incremental values of the primary variable. The ratio is represented by the parameter K as:

$$
\mathcal{K} = \frac{\Delta \mathcal{A}\_{\text{max}}}{\Delta \mathcal{A}\_{\text{x}}} \,, \tag{89}
$$

where ΔAmax is the maximum value of ΔA over all the crystals, all the integration points are in finite element mesh, and ΔA*<sup>x</sup>* is a desired incremental value of ΔA in the numerical algorithm. The value of ΔA is normally regulated by the computational performance and accuracy of the numerical integration procedure. In a fully implicit numerical scheme, the computational performance does not usually becomes a challenge, but the accuracy does. Therefore, the value of ΔA is mainly controlled by the accuracy of the numerical solution. In accordance with the two-level iterative scheme, A is defined as a set of three primary variables A = {A<sup>1</sup> A2}, where

$$\mathcal{A}\_1 = \mathbf{T}^{\epsilon} \; , \; \mathcal{A}\_2 = s\_r^{\alpha} \; . \tag{90}$$

In the proposed time sub-stepping algorithm, three ranges of values are defined for the parameter K: (i) If the parameter K exceeds 1.25, the estimated value of A is rejected, and the new time increment (25% smaller than the previous) will be defined. This condition makes sure that the difference between the maximum and desired value of A will not reach beyond 25%, which could produce an inaccurate solution. (ii) If K lies in the range of 0.8 to 1.25, then the estimated solution is accepted, defining the new time step as Δ*tn*+<sup>1</sup> = (Δ*t*)*n*/K. In this condition, the new time step size Δ*tn*+<sup>1</sup> is more or less identical to the previous Δ*tn*. (iii) If K is less than 0.8, then the estimated solution is assumed to be converged, and a new time step size is defined that is 25% larger than the previous. This condition ensures that the solution is well converged, so the time increment could be increased to reduce computational time. A summary of the time sub-stepping algorithm is given in Table 1. The numerical integration scheme for elastic-plastic deformation of a crystal based on crystal plasticity formulations is summarized in Figure 2.

**Table 1.** Time sub-stepping algorithm for Newton–Raphson Iterative Scheme.


**Figure 2.** Numerical integration algorithm for crystal plasticity model.

#### *3.5. Summary of Numerical Integration Algorithm*

The presented numerical integration scheme for the crystal plasticity model is developed in a generalized framework. The integration scheme can also be used for a slip-based crystal plasticity model. In this case, the number of primary variables will reduce from five to three, that is, {**T***<sup>e</sup>* , *s<sup>α</sup> <sup>r</sup>* , *υ<sup>i</sup>* , *s<sup>i</sup> <sup>t</sup>*, **R***<sup>e</sup>* , } to {*T<sup>e</sup>* , *s<sup>α</sup> <sup>r</sup>* , **R***e*}, by omitting twin resistance *s<sup>i</sup> <sup>t</sup>* and twinned martensite volume fraction *υ<sup>i</sup>* . However, the two-level Newton–Raphson iterative scheme remains two-level.

#### **4. Finite Element Modeling**

The numerical integration scheme of the twin-based crystal plasticity model is validated and further used to predict the deformation behavior of metals through finite element simulations. For this, numerical simulations are performed for single-crystal and polycrystal FCC-austenite subjected to biaxial and combined tension-shear loading using finite element software ABAQUS. The material model's constitutive equations are incorporated in finite element simulations through a user-defined material subroutine (UMAT). A material point in a single crystal of austenite is modeled through the eight-node brick element of unit side length with reduced integration (C3D8R). For polycrystalline simulations, each finite element represents 500 grains of random crystallographic texture. The random texture of grains, expressed in Euler angles, is developed using Kocks convention [59]. In addition, a weighted average procedure is utilized to estimate the cumulative response of polycrystalline austenite material. Furthermore, published experimental results of TWIP steel are referred to for validating the developed numerical model. Consequently, the deformation behavior of austenite-based TWIP steel, subjected to different loading conditions, is estimated and analyzed through finite element simulations.

## *4.1. Geometry and Boundary Conditions*

In finite element simulations, two modes of loadings are considered: (i) uniaxial tension, and (ii) uniaxial compression. In uniaxial tension, a displacement of +0.15 mm is applied on a cube face, as shown in Figure 3a. The planar symmetric boundary condition is employed on three faces, while the two remaining surfaces are traction-free. Similar boundary conditions are adopted in uniaxial compression, except for a negative displacement of 0.15 mm on an element surface, as illustrated in Figure 3b. All loading conditions follow an analogous displacement rate, 1000 increments in a logical time bound of 0 to 1, of 1.5 × 10−<sup>4</sup> mm/time. The effect of texture on the deformation pattern is incorporated through simulations of three crystallographic orientations, as illustrated in Figure 4. These are represented in two domains, that is, crystal and specimen coordinate systems, which are, respectively represented by the (*e<sup>c</sup>* 1,*e<sup>c</sup>* 2,*e<sup>c</sup>* <sup>3</sup>) and (*e<sup>s</sup>* 1,*e<sup>s</sup>* 2,*e<sup>s</sup>* <sup>3</sup>) axes. The crystal direction [100], corresponding to Euler angles (*φ*<sup>1</sup> <sup>1</sup>, *<sup>φ</sup>*<sup>1</sup> <sup>2</sup>, *<sup>φ</sup>*<sup>1</sup> <sup>3</sup>), is equal to (−90o, 0, 90o) (see Figure 4a). Likewise, [110] and [111], corresponding to (*φ*<sup>2</sup> <sup>1</sup>, *<sup>φ</sup>*<sup>2</sup> <sup>2</sup>, *<sup>φ</sup>*<sup>2</sup> <sup>3</sup>) and (*φ*<sup>3</sup> <sup>1</sup>, *<sup>φ</sup>*<sup>3</sup> <sup>2</sup>, *<sup>φ</sup>*<sup>3</sup> <sup>3</sup>), are equivalent to (−45o, 0, 90o) and (−45o, 35.26o, 54.74o), respectively, as represented in Figure 4b,c.

**Figure 3.** Finite element models subjected to (**a**) tension and (**b**) compression.

**Figure 4.** Loaded directions (LDs) in single crystal with respect to specimen orientation (**a**) [100]-LD , (**b**) [110]-LD, and (**c**) [111]-LD.

#### *4.2. Modeling Parameters*

In the present work, simulation parameters (material, hardening, and thermal) are considered as in [51,53,57], as summarized in Table 2.

**Table 2.** General modeling parameters of TWIP steels.


#### **5. Results and Discussion**

In the current section, reported experimental investigations are utilized to validate the developed numerical integration scheme. Afterward, further finite element simulations of single-crystal and polycrystal austenite-based TWIP steel, subjected to uniaxial tension and compression, are executed and investigated. The prime purpose of these simulations is to test and verify the developed numerical scheme under different material deformation behaviors.

#### *5.1. Model Validation*

The developed numerical integration scheme is validated through the published experimental results of twinning-induced plasticity (TWIP) steels. In this, the uniaxial tensile test results, as reported in [60], of three TWIP steels with different chemical compositions are utilized. These steels are: TWIP 22% Mn-0.6% C, TWIP 30% Mn-0.5% C, and TWIP

29% Mn-0.8% C. In finite element simulations of uniaxial tension, displacement of 0.45 to 0.6 mm is applied on a surface of a cubic element with outward normal parallel to *e<sup>s</sup>* 1 axis, as shown in Figure 4. In uniaxial compression, displacement of 35.0 × 10−<sup>2</sup> mm is applied along *e<sup>s</sup>* <sup>1</sup> vector on a surface with *<sup>e</sup><sup>s</sup>* <sup>1</sup> as area normal vector. It is assumed that the cubic element of a material point is comprised of 500 grains with random texture. The orientation distribution function of grains is expressed in terms of Euler angles through the Kocks convention. The general modeling parameters of TWIP steels are summarized in Table 2, while parameters specific to TWIP steels 1, 2, and 3 are listed in Table 3.

The simulation results of uniaxial tension are in good agreement with the experimental observations of TWIP steels, as evident from Figure 5. However, it is noted that for TWIP 22% Mn-0.6% C and TWIP 29% Mn-0.8% C, the experimental results are under-predicted at higher strain (>0.35). For instance, the maximum absolute percent relative error ((experimental stress − simulation value)/experimental stress) in equivalent stress for TWIP 22% Mn-0.6% C at 0.43 equivalent strain is 11.63%. In TWIP 29% Mn-0.8% C, it is 7.44% at 0.45 strain. The deformation behavior of TWIP steels is solely dependent on the complex interactions of slip and twin planes, especially at higher strains. These interactions are the function of slip and twin planes' resistances and orientations. At higher strains, many slip and twin planes may activate and interact, which may cause a higher magnitude of strain hardening, as shown in the experimental results of Figure 5. This complicated phenomenon of higher strain hardening may not be fully encapsulated in numerical simulations. Nevertheless, the numerical model predicts experimental results of TWIP steels well. Therefore, it could be further used to simulate the deformation behavior of these steels, using similar modeling parameters, subjected to multiple types of loading conditions.


**Table 3.** Specific modeling parameters of TWIP 22% Mn-0.6% C, TWIP 30% Mn-0.5% C, and TWIP 29% Mn-0.8% C steels.

**Figure 5.** Comparison of experimental and simulation results of TWIP steels subjected to uniaxial tension: (**a**) TWIP 22% Mn-0.6% C, (**b**) TWIP 30% Mn-0.5% C, (**c**) TWIP 29% Mn-0.8% C.

#### *5.2. Finite Element Simulations*

After establishing the modeling parameters, further finite element simulations are performed to evaluate the deformation pattern of single-crystal and polycrystal TWIP steels subjected to uniaxial tension and compression. The other reasons for the simulations are to observe the effects of loading directions on the (i) activity of slip and twin-systems, (ii) magnitude of slip and twin shear strain, and (iii) volume fraction of the twinned region. In the subsequent sections, the variation of these parameters for three crystallographic directions ([100], [110], and [111]) in a single crystal of TWIP steels are analyzed. For polycrystal, the material point is represented by a set of five hundred randomly oriented grains. A detailed discussion about the choice of the optimum number of grains is presented in [53].

#### 5.2.1. Slip and Twin Planes' Activity

In twinning- and transformation-induced plasticity steels, it would be worthwhile to investigate the contribution of slip and twinning in overall plastic strain. Moreover, their effects on deformation and hardening behaviors are equally essential. Therefore, the activity of slip and twin systems during the course of single crystals' deformation is evaluated through the current model in three crystallographic orientations, as illustrated in Figures 6 and 7. It is noted that slip and twin systems of austenite single crystal, as mentioned in Appendix A, are designated through numbers, as shown in these figures. For example, slip system 1 is the combination of slip plane normal, *n<sup>α</sup>* and slip direction, and *mα <sup>k</sup>* vectors, where *<sup>α</sup>* and *<sup>k</sup>* are equal to 1, as shown in Table A1. In addition, vectors *<sup>n</sup>*<sup>1</sup> and *m*1 <sup>2</sup> form slip system 2. A similar convention is used for twin systems; see Table A2. The activity of slip and twin systems is assessed based on the ratios of resolved shear stress and slip or twin resistance. The slip or twin plane becomes active once these ratios are greater than 1.

**Figure 6.** Activity of slip planes in single crystal of TWIP 1, TWIP 2, and TWIP 3 steels under tension and compression in three crystallographic directions.

**Figure 7.** Activity of twin planes in single crystal of TWIP 1, TWIP 2, and TWIP 3 steels under tension and compression in three crystallographic directions.

It is evident from Figures 6 and 7 that all TWIP steels show similar activity of slip and twin systems, regardless of their varying compositions. Furthermore, a similar number of slip and twin systems are active while the crystal is subjected to tension and compression in the [100] and [110] directions. However, this is not the case in [111], where a higher number of systems show activity in tension (slip systems: 1, 3–5, 8, and 9; twin systems: 2, 6, and 7) than compression (slip systems: 4, 8, 9, and 12; twin systems: 6 and 7). It is also found that some slip and twin systems become activated at lower strain but deactivated at higher, or vice versa. This is probably due to the interaction among slip–slip, slip–twin, or twin–twin systems, and/or reorientation of the slip or twin planes during the course of deformation. These interactions are one of the main causes of material softening or hardening at the micro-level. Another important finding from Figures 6 and 7 is related to the contribution of the slip and twin modes in overall plasticity. As a whole, the contribution of slip is higher for all crystallographic directions. However, the highest level of twin activity is observed in the [100] direction.

#### 5.2.2. Deformation and Hardening Behavior

The stress-strain response of single-crystal and polycrystal TWIP steels under tension and compression are represented, respectively, in Figures 8 and 9. It is noted that from now on, the conventions of TWIP 1, TWIP 2, and TWIP 3 are used for TWIP 22% Mn-0.6% C, TWIP 30% Mn-0.5% C, and TWIP 29% Mn-0.8% C, respectively. A prominent variation in deformation pattern, both in tension and compression, is seen in crystallographic orientations and polycrystal.

**Figure 8.** Deformation behavior of single-crystal and polycrystal TWIP steels subjected to tension.

It is noted in Figures 8 and 9 that orientations [100] and [110] show similar behavior in tension and compression, except for tension in TWIP 1. In this, a sample loaded in [110] initially shows hardening, and then softening after equivalent strain 0.4. The softening behavior may be induced due to the activation of mechanical twinning; however, this does not comply with twin systems' activity, as shown in Figure 6b, where a similar number of twin systems are active at strain 0.4. The other possible reasons may include: (i) reorientation of slip and/or twin systems that may enhance the overall shear strain rates of both modes, or (ii) variation in crystal defect energy as a result of dislocations' interaction. In all TWIP steels, crystals subjected to tensile load in the [111] direction present the largest magnitude of stress; however, in compression, a similar pattern is observed until the equivalent strain 0.4. After this, the [111] direction shows a lower

magnitude of equivalent stress. Furthermore, it is also noted that all polycrystal TWIP steels represent a higher magnitude of stress in tension than compression. This may indicate a greater dominancy of slip and twin systems' reorientation (primary hardening) and interactions (latent hardening) in tension than compression. A quantitative comparison of stress magnitudes in single-crystal and polycrystal TWIP steels is presented in Table 4.

**Figure 9.** Deformation behavior of single-crystal and polycrystal TWIP steels subjected to compression.

**Table 4.** Comparison of equivalent stress at 0.4 equivalent strain of single-crystal and polycrystal TWIP steels under tension and compression.


#### 5.2.3. Slip and Twin Shear Strain

The magnitude of shear strain is another important parameter in crystal plasticity, especially if modes other than slip are also favorable. In addition, the slip and twin contribution in overall plasticity can only be quantitatively represented through shear strain magnitude, not by activity of slip and twin planes, as it is a qualitative measurement. In this regard, the magnitudes of shear strain in slip and twin modes are analyzed. For almost all TWIP steels, a linear variation of shear strain, under tension and compression, is observed in slip mode; however, nonlinearity is dominant in twin, as illustrated in Figures 10 and 11. Moreover, the magnitude of the shear strain in slip is far greater than twin in all steels' crystallographic directions, which refers to the dominancy of slip over twin mode in plasticity. In slip, all crystallographic directions represent, with few exceptions, a similar magnitude of shear strain under the same kind of loading; however, a significant

variation is observed in twin mode. It is evident from Figure 11 that shear strain becomes nearly constant for all crystallographic directions of TWIP 1 and 2 after specific equivalent strain. On the contrary, shear strain's magnitude continuously increases in TWIP 3. This peculiar behavior could be an indication of latent hardening and planes' interaction effects. A quantifiable observation of shear strain in slip and twin modes is exhibited in Table 5.

**Figure 10.** Magnitude of commulative slip shear strain of a single crystal of TWIP steels subjected to tension and compression.

**Figure 11.** Magnitude of commulative twin shear strain of a single crystal of TWIP steels subjected to tension and compression.


**Table 5.** Shear strain's magnitude at 0.4 equivalent strain in slip and twin modes of TWIP 1, 2, and 3 steels under tension and compression.

An obvious observation from Table 5 is the highest magnitude of slip shear strain of TWIP 3 under tension and compression in all crystallographic directions. On the other hand, twin shear strain shows the uppermost values in TWIP 2 for all directions under both types of loading.
