**1. Introduction**

Advanced materials have a vital place along with other key technologies in the fourth industrial revolution (4IR) [1]. Technological developments and achievements depend (and will continue to depend) on the availability of advanced materials. In addition, advanced manufacturing techniques make it possible to produce a range of products, specifically for adverse and corrosive environments and cryogenic applications. In particular, manufacturing products with geometrically complex and enhanced properties becomes possible due to the development of advanced metallic alloys. However, regardless of the availability of the class of metallic alloys, several issues and limitations still restrict advancement in product development. According to the material point of view, one of the restrictions is the simultaneous requirement of high tensile strength and ductility. These properties become extremely crucial in large-deformation applications, such as superplasticity, sheet metal forming, cold rolling, and so on, where high strength with excellent formability is required to obtain highly deformed products. In traditional metallic alloys, these properties cannot be enhanced simultaneously; rather, improvement in one can only be achieved through the other's detriment. This long-lasting issue is resolved by developing innovative and advanced high-strength steels (AHSS) [2]. One of the prime characteristics of AHSS is an excellent balance between tensile strength and ductility. This makes AHSS an optimum

**Citation:** Khan, R.; Pervez, T.; Alfozan, A.; Qamar, S.Z.; Mohsin, S. Numerical Modeling and Simulations of Twinning-Induced Plasticity Using Crystal Plasticity Finite Element Method. *Crystals* **2022**, *12*, 930. https://doi.org/10.3390/ cryst12070930

Academic Editor: Wojciech Polkowski

Received: 24 March 2022 Accepted: 26 May 2022 Published: 30 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

choice in the automotive, aerospace, and oil and gas industries [3–5]. The classification of AHSS in three generations is based on the nature of the microstructure, phases, and composition [2,4]. The prominent members of the second and third generations are transformation and twinning-induced plasticity steels, respectively.

Twinning-induced plasticity (TWIP) steels have broad range of applications due to their promising combination of work hardening, formability, and tensile strength [6–8]. These characteristics make it a dominant candidate for many advanced material applications. The special amalgamation of TWIP steels' properties is achieved through controlled microstructure, and the fraction of primary and secondary phases [9]. The primary phase in these steels is retained austenite, while the secondary may include ferrite, martensite, and sometimes bainite. The main alloying element, which plays a significant role in enhancing the properties, is manganese [10]. The weight percentage of manganese varies, but it is normally greater than 15–20%. Due to the high percentage of manganese, the stacking fault energy (SFE) of TWIP steels is relatively lower [11], in the span of 20 to 40 mJ/m2, than the other class of AHSS. The low magnitude of SFE favors activation of a secondary mode of plasticity, which is twinning [12]. The volume fraction of twinning governs the mechanical properties—more specifically, the strain hardening—of TWIP steels. Since twin regions behave as obstacles in dislocation glide, the dislocation mean free path may reduce. This eventually improves the strain hardening of TWIP steels [7,13,14]. Furthermore, in Fe–Mn–C grade, twin regions normally comprise a huge magnitude of sessile dislocations' density. The high magnitude of density results from twin nucleation and growth, which act as resilient inclusions that may hinder dislocation glide [15,16]. It can be stated that the excellent combination of strength and ductility in TWIP steels is primarily due to work hardening, which may be induced in a material due to the twinning mechanism. The secondary mode of plastic deformation, mechanical twinning, may occur in metals (in conjunction with slip), non-metals, and metallic glasses [6,7,17]. However, the fraction of twinning is highly dependent on the chemical composition of the alloy, as well as the physical conditions [18–21]. Moreover, stacking fault energy (SFE) plays an important role in defining the potential for the activation of twinning [22–24].

One of researchers' foremost challenges was to computationally couple two modes of plasticity, slip and twinning, in predicting thedeformation behavior of metals [25–27]. Among these obstacles, one was to computationally account for the huge number of twin orientations that may occur during the course of deformation [25,26]. Multiple solutions have been proposed to overcome this long-lasting issue. These mainly include Taylor's least work hypothesis [28] to minimize the orientation factor [29]; utilizing the statistical technique in determining the re-orientation of whole grain to account for the total number of grain orientations in computation [26]; and employing the condition of weighted grain orientations to avoid the generation of new grains [25]. Kalidindi et al. [27] presented another possible solution to incorporate mechanical twinning in crystal plasticity theory. It is based on the concept of multiplicative decomposition of the total deformation gradient into elastic and plastic, as initially proposed by Asaro and Rice [30]. Furthermore, plastic deformation gradient and strain hardening effects are defined according to slip and mechanical twinning. In addition, the rate of change of twin volume fraction depends on the resolved shear stress and twin resistance of potentially active twin systems, as in the definition of slip system hardening by Asaro and Rice [30]. In subsequent work, several attempts have been made to predict strain-hardening effects of mechanical twinning [31–33]. In these models, twin-related strain-hardening effects of *α*-titanium alloys are incorporated in crystal plasticity theory. The models have utilized the crystal plasticity finite element method to predict metallic alloys' deformation behavior, which undergoes slip and twinning.

The crystal plasticity finite element method (CPFEM) is frequently utilized to model mechanical twinning in shape-memory alloys [34–36], advanced high-strength steels [37–39], magnesium alloys [40], and high-Mn austenitic steels [19,41,42]. It has also found huge application in modeling and simulating the deformation and damage behaviors of transformationinduced plasticity (TRIP), twinning-induced plasticity (TWIP), and multiphase steels [43–45]. More specifically, Gui et al. [43] have presented a multi-mechanism and microstructurebased crystal plasticity model for estimating the shear deformation behavior of TRIP steel under cyclic load. In this, the simulation results show that the samples subjected to high strain magnitude exhibit stronger cyclic shear hardening and that activation of martensitic transformation promotes cyclic hardening. Furthermore, Qayyum et al. [44] have presented a novel physics-based crystal plasticity model with ductile damage criterion for predicting the damage behavior of austenite-based TRIP steel (X8CrMnNi16-6-6) with 10% zirconia particles. Among other conclusions, it was found that the energy absorbed by the zirconia particles in the course of deformation is comparatively high, and there is substantial stiffness degradation in the bulk material. These factors significantly influenced the composite material behavior. In a continuation of this work, Qayyum et al. [45] have utilized a similar physics-based crystal plasticity numerical modeling technique to create a semi-automatic virtual laboratory to analyze and create customized functional multiphase materials. The CPFEM has also successfully implemented in modeling the behavior of TWIP steels. More specifically, it has been implemented to estimate: elastic properties of single-crystal TWIP steel through nano-indentation [46]; deformation behavior, texture evolution, and earing mechanism in TWIP steels [47,48]; and combined effects of slip, mechanical twinning, and martensitic transformation on the overall behavior of high-Mn steels [49–51]. In particular, Madivala et al. [52] have investigated the strain-hardening and fracture behavior of high-manganese austenitic TWIP steel at temperatures ranging from 123 K to 773 K. It was observed that twinning becomes the dominant deformation mechanism at 298 K, and the twin fraction increases with temperature until a transition temperature of about 473 K. Beyond this, dislocation glide alone becomes dominant, instead of twinning and dislocation glide. Recently, Khan et al. [53] presented a micromechanical model of twinning-induced plasticity using crystal plasticity and thermodynamic frameworks. The deformation gradients resulting from crystallographic slip and mechanical twinning are modeled through the kinematic decomposition of the total deformation gradient. The constitutive formulation of dissipated energy and Helmholtz free energy and the driving potentials for inelastic deformation modes are represented through a thermodynamic framework. The deformation gradients resulting from crystallographic slip and mechanical twinning are modeled through the kinematic decomposition of the total deformation gradient. Finally, a numerical integration scheme is used to incorporate the constitutive formulation in commercial finite element code ABAQUS through the user-defined material subroutine. It was observed that when the material point in the single crystal is subjected to tension, twin deformation plays a dominant role, while the reverse is observed in compression. Furthermore, in both tension and compression, the variation in the volume fraction of twinned martensite is found in all crystallographic directions (i.e., [100], [110], and [111]), but with different magnitudes.

The current work is an extension of the micromechanical model presented by Khan et al. [53], where a novel numerical integration scheme is developed. In this, firstly, the constitutive equations of the micromechanical model of twinning-induced plasticity are presented. Secondly, a numerical integration scheme is utilized to update the constitutive equations using a fully implicit time integration procedure based on the backward Euler method. Thirdly, a three-level iterative scheme is developed to solve the coupled nonlinear system of equations through the Newton–Raphson method. In this, a L<sup>2</sup> (two norm) convergence criterion is used to estimate the response of incremental values of state variables. A time sub-stepping algorithm is incorporated with the numerical integration scheme to improve the convergence and reduce computational time. The developed numerical scheme is then implemented as a user-defined material subroutine in the finite element software ABAQUS 6.14. The model is then validated through published experimental observations of TWIP steels. Finally, the deformation behavior of TWIP steel under different loading conditions is estimated through finite element simulations.

#### **2. Constitutive Modeling**

In this part, constitutive equations of the model developed by [53] are presented, where twinning-induced plasticity is incorporated with slip-based crystal plasticity theory. The constitutive equations include the formulations of: (i) multiplicative decomposition of the total deformation gradient into elastic and plastic parts, (ii) driving potentials for slip and twining, (iii) dissipated energy as a result of plastic flow, (iv) recovered energy due to elastic deformation, (v) plastic flow rule due to slip and twinning, and (vi) hardening law.

In the subsequent sections, standard symbols are employed for designating tensors and their operations. The tensor and vector quantities are, respectively, expressed through the capital and small bold letters. Fourth-, second-, and first-order tensors are symbolized as C, **A**, and **a**, respectively. The notations **a** ⊗ **b** and **Aa** represent the dyadic product of vectors and the contraction of second-order tensor with vector, respectively. The mathematical operations of second-order tensors are illustrated as: (i) inner product: **AB**, (ii) dyadic product: **A** ⊗ **B**, and (iii) scalar product: **A**:**B**. The fourth- and second-order tensors' contraction is expressed as C:**A**. Any non-standardized notation will be defined explicitly.

#### *2.1. Kinematic-Based Modeling*

In all successive sections of mathematical modeling, only the final equations of each constitutive formulation are presented. For detailed derivations, the readers are advised to refer to [53,54].

#### Multiplicative Decomposition of Deformation Gradient

The decomposition of the total deformation gradient into elastic and plastic parts, as discussed by [55], can be represented as:

$$\mathbf{F} = \mathbf{F}^c \mathbf{F}^p,\tag{1}$$

where **F** is the total deformation gradient, while **F***<sup>e</sup>* and **F***<sup>p</sup>* are, respectively, elastic and plastic deformation gradients. The elastic part is further categorized in symmetric left stretch **V***<sup>e</sup>* and orthogonal rotation **R***<sup>e</sup>* tensors, **F***<sup>e</sup>* = **V***e***R***<sup>e</sup>* . The plastic deformation gradient incorporates crystallographic slip and mechanical twinning. The rotation and plastic deformation gradient tensors can be combined into a plastic rigid rotation tensor, **F**∗, as discussed in [54]. The overall decomposition of deformation gradients can be represented by Equation (2) as:

$$\mathbf{F} = \mathbf{V}^{\varepsilon} \mathbf{F}^{\*} \,, \quad \text{where} \quad \mathbf{F}^{\*} = \mathbf{R}^{\varepsilon} \mathbf{F}^{p} = \mathbf{R}^{\varepsilon} \mathbf{F}\_{t}^{p} \mathbf{F}\_{s}^{p} \,. \tag{2}$$

In Equation (2), the slip and twinning parts of **F***<sup>p</sup>* are, respectively, represented as **F***<sup>p</sup> s* and **F***<sup>p</sup> <sup>t</sup>* ; **<sup>V</sup>***<sup>e</sup>* is the symmetric left stretch tensor; and **<sup>R</sup>***<sup>e</sup>* is an orthogonal rotation tensor. The elastic-plastic behavior of a material point involves slip and twin deformation modes, which can be elaborated through kinematic decomposition, as illustrated in Figure 1.

The undeformed and deformed configurations of a material point are represented by **Ω<sup>0</sup>** and **Ωf**, respectively. In this, total deformation gradient is disintegrated into three intermediate configurations **<sup>Ω</sup>I**, **<sup>Ω</sup>˘ II**, and **<sup>Ω</sup> III**. The first and second intermediate configurations represent plastic deformation due to slip and mechanical twinning, respectively, while the third shows rigid body rotation. The deformed state, **Ωf**, is projected from the third intermediate configuration through the stretch tensor **V***<sup>e</sup>* . In the present model, the relaxed (third) intermediate configuration **Ω III** is adopted for representing the constitutive equations. Theoretically, this configuration is obtained by elastically unloading, using (**V***e*)−<sup>1</sup> without rotation, from the current (deformed) state to a stress-free configuration [54].

**Figure 1.** Kinematic decomposition of a material point into three intermediate configurations: Ω¯ *<sup>I</sup>*, Ω˘ *I I*, and Ω˜ *III*.

Referring to the kinematic decomposition and definition of the velocity gradient, the velocity gradient in the current (deformed) configuration can be presented as:

$$\mathbf{L} = \dot{\mathbf{F}} \mathbf{F}^{-1} \quad = \dot{\mathbf{V}}^{\varepsilon} (\mathbf{V}^{\varepsilon})^{-1} + \mathbf{V}^{\varepsilon} \widetilde{\mathbf{L}}^{\*} (\mathbf{V}^{\varepsilon})^{-1}. \tag{3}$$

Here, **<sup>L</sup>**<sup>∗</sup> is the velocity gradient in third intermediate configuration. By considering Equations (1) and (2), it can be written as:

$$
\widetilde{\mathbf{L}}^{\*} = \dot{\mathbf{F}}^{\*} (\mathbf{F}^{\*})^{-1} \quad = \dot{\mathbf{R}}^{c} (\mathbf{R}^{c})^{-1} + \mathbf{R}^{c} \check{\mathbf{L}}^{p} (\mathbf{R}^{c})^{-1}, \tag{4}
$$

where **L**˘ *<sup>p</sup>*, the plastic velocity gradient in second intermediate configuration **Ω˘ II**, incorporates plastic flow due to crystallographic slip and mechanical twinning in the constitutive model. Since the plastic deformation gradient can be divided into slip and twin contributions, **F***<sup>p</sup>* = **F***<sup>p</sup> t* **F***p <sup>s</sup>* , it can be written as:

$$\mathbb{L}^p = \mathbb{F}^p(\mathbb{F}^p)^{-1} = \mathbb{L}\_s^p + \mathbb{F}\_s^p \mathbb{L}\_t^p (\mathbb{F}\_s^p)^{-1} \,. \tag{5}$$

In Equation (5), **L**¯ *<sup>p</sup> <sup>s</sup>* and **<sup>L</sup>**˘ *<sup>p</sup> <sup>t</sup>* represent the plastic velocity gradients of slip and twinning, respectively. By using the definition of the velocity gradient, Equation (5) can be transformed as:

$$\mathbf{L}^p = \sum\_{a=1}^{N\_s} \dot{\gamma}^a \mathbf{S}^a + \mathbf{F}\_s^p \left\{ \sum\_{i=1}^{N\_l} \gamma^i \mathbf{S}^i \right\} (\mathbf{F}\_s^p)^{-1} \,, \tag{6}$$

where *α* represents the slip system's number as (*α* = 1, ... , *Ns*); *Ns* is the total number of slip systems; *γ*˙ *<sup>α</sup>* is the plastic shear strain rate of *α*-slip system; **S**¯ *<sup>α</sup>* is the Schmid orientation tensor in first intermediate configuration (represented by slip direction, **m**¯ *<sup>α</sup>*, and slip plane area normal vectors, **n**¯ *<sup>α</sup>*, as **S**¯ *<sup>α</sup>* = **m**¯ *<sup>α</sup>* ⊗ **n**¯ *<sup>α</sup>*, where *i* denotes the twin system (*i* = 1, ... , *Nt*), and *Nt* represents the entire number of twin systems); *γ*˙ *<sup>i</sup>* is the plastic shear strain rate of *i*-twin system; and **S**˘*<sup>i</sup>* is the twin orientation tensor in second intermediate state (expressed through twin direction, **m**˘ *<sup>i</sup>* , and twin plane area normal, **n**˘*<sup>i</sup>* , vectors as **S**˘*<sup>i</sup>* = **m**˘ *<sup>i</sup>* ⊗ **n**˘*<sup>i</sup>* ). In the presented model, mechanical twinning is assumed to occur within the plastically deformed region by crystallographic slip. This leads to the inclusion of the volume fraction of each region (slip and twin) in the velocity gradient equation as folllows:

$$\check{\mathbf{L}}^p = \left(1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^i\right) \sum\_{a=1}^{N\_s} \boldsymbol{\gamma}^a \bar{\mathbf{S}}^a + \mathbf{F}\_s^p \left\{ \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^i \boldsymbol{\dot{\gamma}}^i \bar{\mathbf{S}}^i \right\} \left(\mathbf{F}\_s^p\right)^{-1} \,. \tag{7}$$

Here, *υ<sup>i</sup>* is the volume fraction of *i*-twin system in second configuration. Similarly, the plastic velocity gradient in third intermediate configuration can be expressed as:

$$\mathbf{\tilde{L}}^{\*} = \mathbf{\tilde{O}}^{\varepsilon} + \mathbf{R}^{\varepsilon} \left[ \left( 1 - \sum\_{i=1}^{N\_{\ell}} \boldsymbol{\upsilon}^{i} \right) \sum\_{a=1}^{N\_{\ell}} \boldsymbol{\gamma}^{a} \mathbf{\tilde{S}}^{a} + \mathbf{F}\_{s}^{p} \left\{ \sum\_{i=1}^{N\_{\ell}} \boldsymbol{\upsilon}^{i} \boldsymbol{\dot{\gamma}}^{i} \mathbf{\tilde{S}}^{i} \right\} \left( \mathbf{F}\_{s}^{p} \right)^{-1} \right] \left( \mathbf{R}^{\varepsilon} \right)^{T},\tag{8}$$

where **<sup>Θ</sup>** *<sup>e</sup>*(= **<sup>R</sup>**˙*e*(**R***e*)*T*) is the spin tensor. The Schmid tensor— **<sup>S</sup>**¯ *<sup>α</sup>* from first and twin orientation tensor, **S**˘*<sup>i</sup>* from second states—can be transformed to the third intermediate configuration through forward conversion procedure using rigid rotation tensor, to yield the following:

$$
\breve{\mathbf{S}}^{\alpha} = \mathbf{R}^{\varepsilon} \mathbf{S}^{\alpha} (\mathbf{R}^{\varepsilon})^{T} \,, \tag{9}
$$

$$
\tilde{\mathbf{S}}^i = \mathbf{R}^\varepsilon \tilde{\mathbf{S}}^i (\mathbf{R}^\varepsilon)^T. \tag{10}
$$

The final form of velocity gradient in the third intermediate configuration can be expressed as:

$$\check{\mathbf{L}}^{\*} = \check{\mathbf{G}}^{\varepsilon} + \left(1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^{i}\right) \sum\_{a=1}^{N\_b} \boldsymbol{\gamma}^{a} \tilde{\mathbf{S}}^{a} + \mathbf{F}\_s^p \left\{ \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^{i} \boldsymbol{\gamma}^{i} \tilde{\mathbf{S}}^{i} \right\} (\mathbf{F}\_s^p)^{-1} \,. \tag{11}$$

The symmetric part of the velocity gradient, as mentioned in Equation (3), is given as:

$$\mathbf{D} = \breve{\mathbf{D}}^{\varepsilon} + \frac{1}{2} \left[ (\mathbf{V}^{\varepsilon})^{-T} \breve{\mathbf{C}}^{\varepsilon} \breve{\mathbf{L}}^{\ast} (\mathbf{V}^{\varepsilon})^{-1} + (\mathbf{V}^{\varepsilon})^{-T} (\breve{\mathbf{C}}^{\varepsilon} \breve{\mathbf{L}}^{\ast})^{T} (\mathbf{V}^{\varepsilon})^{-1} \right], \tag{12}$$

where **<sup>D</sup>** *<sup>e</sup>* is the symmetric component of **<sup>V</sup>**˙ *<sup>e</sup>*(**V***e*)<sup>−</sup>1, and **<sup>C</sup>***<sup>e</sup>* = (**V***e*)(**V***e*)*<sup>T</sup>* = (**V***e*)<sup>2</sup> is the right Cauchy–Green deformation tensor in the third intermediate configuration. Equation (12) can also be represented as:

$$\mathbf{D} = \widetilde{\mathbf{D}}^{\varepsilon} + (\mathbf{V}^{\varepsilon})^{-T} \widetilde{\mathbf{D}}^{\*} (\mathbf{V}^{\varepsilon})^{-1} \,. \tag{13}$$

The symmetric part of **<sup>C</sup>***e***L**<sup>∗</sup> is given as:

$$\tilde{\mathbf{D}}^{\*} = \text{sym}(\tilde{\mathbf{C}}^{\varepsilon}\tilde{\mathbf{O}}^{\varepsilon}) + \mathbf{R}^{\varepsilon}\mathbf{D}^{p}(\mathbf{R}^{\varepsilon})^{T},\tag{14}$$

where **D**¯ *<sup>p</sup>* is the symmetric component of **C**¯ *<sup>e</sup>***L**¯ *<sup>p</sup>*. Similarly, the skew-symmetric component of **L** can be represented as:

$$\mathbf{W} = \frac{1}{2} \left[ \mathbf{L} - \mathbf{L}^T \right] = \overline{\mathbf{W}}^c + (\mathbf{V}^c)^{-T} \overline{\mathbf{W}}^\* (\mathbf{V}^c)^{-1} \,, \tag{15}$$

where -**W***<sup>e</sup>* and -**W**<sup>∗</sup> are the skew-symmetric components of **<sup>V</sup>**˙ *<sup>e</sup>*(**V***e*)−<sup>1</sup> and **<sup>C</sup>***e***L**∗, respectively. The component -**W**<sup>∗</sup> can be derived as:

$$
\hat{\mathbf{W}}^{\epsilon} = \text{skew}(\hat{\mathbf{C}}^{\epsilon}\hat{\mathbf{G}}^{\epsilon}) + \mathbf{R}^{\epsilon}\hat{\mathbf{W}}^{p}(\mathbf{R}^{\epsilon})^{T}.\tag{16}
$$

Here, **W**¯ *<sup>p</sup>* is the anti-symmetric part of **C**¯ *<sup>e</sup>***L**¯ *<sup>p</sup>*. The finite Green strain tensor in reference configuration (third intermediate state with reference to current configuration) can be represented in terms of the symmetric part of the velocity gradient through Equation (19).

A procedure of backward mapping, from current to third intermediate configurations, is adopted for developing an equation of the velocity gradient, **L**, that incorporates elastic stretch along with plasticity and rigid body rotation. This can be expressed through Equation (17) as:

$$
\widetilde{\mathbf{L}} = (\mathbf{V}^{\varepsilon})^{-1} \mathbf{L} \mathbf{V}^{\varepsilon} = (\mathbf{V}^{\varepsilon})^{-1} \mathbf{\tilde{V}}^{\varepsilon} + \widetilde{\mathbf{L}}^{\*} \,. \tag{17}
$$

After substituting the value of **<sup>L</sup>**<sup>∗</sup> from Equation (11), Equation (17) becomes

$$\widetilde{\mathbf{L}} = (\mathbf{V}^{\varepsilon})^{-1}\dot{\mathbf{V}}^{\varepsilon} + \widetilde{\mathbf{\Theta}}^{\varepsilon} + \left[ \left( 1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^{i} \right) \sum\_{a=1}^{N\_s} \dot{\boldsymbol{\gamma}}^{a} \widetilde{\mathbf{S}}^{a} + \mathbf{F}\_s^{p} \left\{ \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^{i} \dot{\boldsymbol{\gamma}}^{i} \widetilde{\mathbf{S}}^{i} \right\} (\mathbf{F}\_s^{p})^{-1} \right]. \tag{18}$$

The finite Green strain tensor in reference configuration (third intermediate state with reference to current configuration) can be represented as:

$$
\check{\mathbf{E}}^{\varepsilon} = \frac{1}{2} (\check{\mathbf{C}}^{\varepsilon} - \mathbf{I}) \; , \qquad \dot{\mathbf{E}}^{\varepsilon} = \frac{1}{2} \mathring{\mathbf{C}}^{\varepsilon} . \tag{19}
$$

Here, **<sup>C</sup>***<sup>e</sup>* is the right Cauchy–Green strain tensor, and **<sup>I</sup>** is a second-order identity tensor. The 2nd Piola–Kirchhoff (PK2) stress, **T***<sup>e</sup>* , in reference configuration (third intermediate configuration) can be represented in terms of finite Green strain tensor. The stresses in slipped and twinned regions can be written in the form of constitutive formulation as:

$$\mathbf{T}\_s^c = \widetilde{\mathbb{C}}^s : \mathbb{B}^c \,, \tag{20}$$

$$\mathbf{T}\_{l}^{\epsilon} = \breve{\mathbb{C}}^{t} : \breve{\mathbb{E}}^{\epsilon} \,, \tag{21}$$

where **T***<sup>e</sup> <sup>s</sup>* and **T***<sup>e</sup> <sup>t</sup>* are the PK2 stress tensors, and <sup>C</sup>*<sup>s</sup>* and <sup>C</sup>*<sup>t</sup>* are the fourth-order elasticity tensors for slipped and twinned regions, respectively. In view of incorporating the effects of slip and twinning in the PK2 stress tensor, an equivalent form can be defined as:

$$\mathbf{T}^{\mathbf{c}} = \mathbb{\hat{C}}^{\mathbf{c}} : \mathbb{\hat{E}}^{\mathbf{c}} \text{ . \tag{22}$$

Here, **<sup>T</sup>***<sup>e</sup>* is an equivalent PK2 stress, and <sup>C</sup>*<sup>e</sup>* is the equivalent elasticity tensor. An equivalent form of elasticity tensor can be expressed through volume averaging technique as:

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

Similarly, Cauchy stress can also be approximated on the basis of volume averaging technique as:

$$\mathbf{T} = \left(1 - \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^i\right) \mathbf{T}\_s + \sum\_{i=1}^{N\_l} \boldsymbol{\upsilon}^i \mathbf{T}\_t. \tag{24}$$

#### *2.2. Kinetic-Based Modeling*

2.2.1. Dissipated Energy Formulation

The dissipated energy density (energy per unit reference volume) in the form of rate of change of entropy is estimated as:

$$\hat{E}\_d = \sigma : \mathbf{\hat{F}} + \rho\_0 \theta \mathbf{\hat{I}} \mathbf{\hat{I}}\_m - \rho\_0 \dot{\epsilon} - \nabla \vec{\mathbf{q}} \tag{25}$$

where *σ* is PK1 (1st Piola–Kirchhoff) stress, **F**˙ is the rate of change of **F**, *ρ*<sup>0</sup> is the density, *θ* is an absolute temperature, Π˙ *<sup>m</sup>* is the rate of change of entropy as a consequence of an external thermomechanical load, ˙ is the rate of change of internal energy per unit mass, and **q** is the heat flux due to temperature variation. The final form of the dissipated energy density rate is given as:

$$\begin{split} E\_d &= \left( 1 - \sum\_{i=1}^{N\_l} \upsilon^i \right) \sum\_{a=1}^{N\_s} \left( \tau^a + \Phi^a - \rho\_0 \frac{\partial \varepsilon}{\partial \check{\zeta}} \Psi^a \right) \dot{\gamma}^a \\ &+ \sum\_{i=1}^{N\_l} \dot{\upsilon}^i \left( \tau^i + \Phi^i - \rho\_0 \frac{\partial \varepsilon}{\partial \check{\zeta}} \Psi^i \right) \dot{\gamma}^i - (\theta) \nabla \Pi\_q \ . \end{split} \tag{26}$$

where *τ<sup>α</sup>* and *τ<sup>i</sup>* are the resolved shear stresses on *α*-slip and *i*-twin systems, respectively; Φ*<sup>α</sup>* and Φ*<sup>i</sup>* are the thermal equivalents of resolved shear stresses *τ<sup>α</sup>* and *τ<sup>i</sup>* , respectively;  is the internal energy density; *ζ* is the crystal defect microstrain parameter; Ψ*<sup>α</sup>* and Ψ*<sup>i</sup>* are stress-like terms and functions of slip resistance of *α*-slip system and twin resistance of *i*-twin system, respectively; ∇ is the del operator; and Π*<sup>q</sup>* is the entropy flux (*q* = Π*qθ*).

#### 2.2.2. Helmholtz Free Energy Formulation

In the micromechanical model presented by [53], Helmholtz free energy (HFE) is evaluated through an additive decomposition of four forms of energies as:

$$E\_h(\mathbf{F}^\varepsilon, \theta, \zeta, \upsilon) = E\_{hm}(\mathbf{F}^\varepsilon) + E\_{ht}(\theta) + E\_{hd}(\zeta) + E\_{hs}(\upsilon) \,, \tag{27}$$

where *Ehm*, *Eht*, *Ehd*, and *Ehs* represent mechanical, thermal, crystal defect, and surface energy components of HFE. In this formulation, HFE depends on four state variables, which are elastic deformation gradient **F***<sup>e</sup>* , absolute temperature *θ*, crystal defect microstrain parameter *ζ*, and twin martensite volume fraction *υ*. In consideration of this, HFE can be derived as:

$$\begin{split} E\_h &= \frac{1}{\rho\_0} \left\{ (\mathbf{F}^c \mathbf{F}^p) : (\mathbf{\tilde{C}}^c : \mathbf{\tilde{E}}^c \mathbf{\tilde{E}}^c) \right\} (\mathbf{F}^p)^T : \left[ \mathbf{F}^c + (\mathbf{F}^c)^T \right]^{-1} \\ &+ \theta \left[ -h\_\epsilon \ln \left( \frac{\theta}{\theta\_r} \right) + \left( h\_{c\eta} - \Pi\_{m,0}^\epsilon \right) \right] + \frac{1}{2\rho\_0} \rho G^\epsilon \zeta^2 \\ &+ \frac{\chi}{\rho\_0 l\_0} \sum\_{i=1}^{N\_l} \upsilon^i \left( 1 - \sum\_{i=1}^{N\_l} \upsilon^i \right) . \end{split} \tag{28}$$

In Equation (28), *he* is an equivalent specific heat; *θ<sup>r</sup>* is the reference temperature; Π*<sup>e</sup> m*,0 is the initial entropy density; *ϕ* is a dimensionless dislocation interaction parameter, which incorporates the effects of dislocations' mobility and their interactions in plasticity; *G<sup>e</sup>* is the equivalent modulus of rigidity; *l*<sup>0</sup> is the length scale parameter; and *χ* is the interfacial energy per unit area.
