Next Article in Journal
Initial Microstructure Effects on Hot Tensile Deformation and Fracture Mechanisms of Ti-5Al-5Mo-5V-1Cr-1Fe Alloy Using In Situ Observation
Next Article in Special Issue
Mechanical Properties of Low Carbon Alloy Steel with Consideration of Prior Fatigue and Plastic Damages
Previous Article in Journal
Characterization of ZTA Composite Ceramic/Ti6Al4V Alloy Joints Brazed by AgCu Filler Alloy Reinforced with One-Dimensional Al18B4O33 Single Crystal
Previous Article in Special Issue
Characterization of Microstructure and High Temperature Compressive Strength of Austenitic Stainless Steel (21-4N) through Powder Metallurgy Route
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Modeling and Simulations of Twinning-Induced Plasticity Using Crystal Plasticity Finite Element Method

1
Mechanical & Industrial Engineering Department, College of Engineering, Imam Mohammad Ibn Saud Islamic University, Riyadh 11432, Saudi Arabia
2
Mechanical & Industrial Engineering Department, College of Engineering, Sultan Qaboos University, Muscat 123, Oman
3
Mechanical Engineering Department, University Teklongi PETRONAS, Seri Iskandar 32610, Perak, Malaysia
*
Author to whom correspondence should be addressed.
Crystals 2022, 12(7), 930; https://doi.org/10.3390/cryst12070930
Submission received: 24 March 2022 / Revised: 25 May 2022 / Accepted: 26 May 2022 / Published: 30 June 2022
(This article belongs to the Special Issue Crystal Plasticity (Volume II))

Abstract

:
In the current work, a fully implicit numerical integration scheme is developed for modeling twinning-induced plasticity using a crystal plasticity framework. Firstly, the constitutive formulation of a twin-based micromechanical model is presented to estimate the deformation behavior of steels with low stacking fault energy. Secondly, a numerical integration scheme is developed for discretizing constitutive equations through a fully implicit time integration scheme using the backward Euler method. A time sub-stepping algorithm and the two-norm convergence criterion are used to regulate time step size and stopping criterion. Afterward, a numerical scheme is implemented in finite element software ABAQUS as a user-defined material subroutine. Finally, finite element simulations are executed for observing the validity, performance, and limitations of the numerical scheme. It is observed that the simulation results are in good agreement with the experimental observations with a maximum error of 16% in the case of equivalent stress and strain. It is also found that the developed model is able to estimate well the deformation behavior, magnitude of slip and twin shear strains, and twin volume fraction of three different TWIP steels where the material point is subjected to tension and compression.

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 choice in the automotive, aerospace, and oil and gas industries [3,4,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,7,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/m 2 , 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,19,20,21]. Moreover, stacking fault energy (SFE) plays an important role in defining the potential for the activation of twinning [22,23,24].
One of researchers’ foremost challenges was to computationally couple two modes of plasticity, slip and twinning, in predicting thedeformation behavior of metals [25,26,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,32,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,35,36], advanced high-strength steels [37,38,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 transformation-induced plasticity (TRIP), twinning-induced plasticity (TWIP), and multiphase steels [43,44,45]. More specifically, Gui et al. [43] have presented a multi-mechanism and microstructure-based 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,50,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 2 (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 A a 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: A B , (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:
F = F e F p ,
where F is the total deformation gradient, while F e and F p are, respectively, elastic and plastic deformation gradients. The elastic part is further categorized in symmetric left stretch V e and orthogonal rotation R e tensors, F e = V e R e . 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:
F = V e F * , where F * = R e F p = R e F t p F s p .
In Equation (2), the slip and twinning parts of F p are, respectively, represented as F s p and F t p ; V e is the symmetric left stretch tensor; and R e 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 Ω 0 and Ω f , respectively. In this, total deformation gradient is disintegrated into three intermediate configurations Ω ¯ I , Ω ˘ II , and Ω ˜ 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 e . 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 ) 1 without rotation, from the current (deformed) state to a stress-free configuration [54].
Referring to the kinematic decomposition and definition of the velocity gradient, the velocity gradient in the current (deformed) configuration can be presented as:
L = F ˙ F 1 = V e ˙ ( V e ) 1 + V e L ˜ * ( V e ) 1 .
Here, L ˜ * is the velocity gradient in third intermediate configuration. By considering Equations (1) and (2), it can be written as:
L ˜ * = F * ˙ ( F * ) 1 = R e ˙ ( R e ) 1 + R e L ˘ p ( R e ) 1 ,
where L ˘ p , 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 p = F t p F s p , it can be written as:
L ˘ p = F ˙ p ( F p ) 1 = L ¯ s p + F s p L ˘ t p ( F s p ) 1 .
In Equation (5), L ¯ s p and L ˘ t p represent the plastic velocity gradients of slip and twinning, respectively. By using the definition of the velocity gradient, Equation (5) can be transformed as:
L ¯ p = α = 1 N s γ ˙ α S ¯ α + F s p { i = 1 N t γ ˙ i S ˘ i } ( F s p ) 1 ,
where α represents the slip system’s number as ( α = 1 , , N s ) ; N s is the total number of slip systems; γ ˙ α is the plastic shear strain rate of α -slip system; S ¯ α is the Schmid orientation tensor in first intermediate configuration (represented by slip direction, m ¯ α , and slip plane area normal vectors, n ¯ α , as S ¯ α = m ¯ α n ¯ α , where i denotes the twin system ( i = 1 , , N t ) , and N t represents the entire number of twin systems); γ ˙ i is the plastic shear strain rate of i-twin system; and S ˘ i is the twin orientation tensor in second intermediate state (expressed through twin direction, m ˘ i , and twin plane area normal, n ˘ i , vectors as S ˘ i = m ˘ i n ˘ i ). 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:
L ˘ p = ( 1 i = 1 N t υ i ) α = 1 N s γ ˙ α S ¯ α + F s p { i = 1 N t υ i γ ˙ i S ˘ i } ( F s p ) 1 .
Here, υ i is the volume fraction of i-twin system in second configuration. Similarly, the plastic velocity gradient in third intermediate configuration can be expressed as:
L ˜ * = Θ ˜ e + R e [ ( 1 i = 1 N t υ i ) α = 1 N s γ ˙ α S ¯ α + F s p { i = 1 N t υ i γ ˙ i S ˘ i } ( F s p ) 1 ] ( R e ) T ,
where Θ ˜ e ( = R e ˙ ( R e ) T ) is the spin tensor. The Schmid tensor— S ¯ α from first and twin orientation tensor, S ˘ i from second states—can be transformed to the third intermediate configuration through forward conversion procedure using rigid rotation tensor, to yield the following:
S ˜ α = R e S ¯ α ( R e ) T ,
S ˜ i = R e S ˘ i ( R e ) T .
The final form of velocity gradient in the third intermediate configuration can be expressed as:
L ˜ * = Θ ˜ e + ( 1 i = 1 N t υ i ) α = 1 N s γ ˙ α S ˜ α + F s p { i = 1 N t υ i γ ˙ i S ˜ i } ( F s p ) 1 .
The symmetric part of the velocity gradient, as mentioned in Equation (3), is given as:
D = D ˜ e + 1 2 [ ( V e ) T C ˜ e L ˜ * ( V e ) 1 + ( V e ) T ( C ˜ e L ˜ * ) T ( V e ) 1 ] ,
where D ˜ e is the symmetric component of V ˙ e ( V e ) 1 , and C ˜ e = ( V e ) ( V e ) T = ( V e ) 2 is the right Cauchy–Green deformation tensor in the third intermediate configuration. Equation (12) can also be represented as:
D = D ˜ e + ( V e ) T D ˜ * ( V e ) 1 .
The symmetric part of C ˜ e L ˜ * is given as:
D ˜ * = sym ( C ˜ e Θ ˜ e ) + R e D ¯ p ( R e ) T ,
where D ¯ p is the symmetric component of C ¯ e L ¯ p . Similarly, the skew-symmetric component of L can be represented as:
W = 1 2 [ L L T ] = W ˜ e + ( V e ) T W ˜ * ( V e ) 1 ,
where W ˜ e and W ˜ * are the skew-symmetric components of V ˙ e ( V e ) 1 and C ˜ e L ˜ * , respectively. The component W ˜ * can be derived as:
W ˜ * = skew ( C ˜ e Θ ˜ e ) + R e W ¯ p ( R e ) T .
Here, W ¯ p is the anti-symmetric part of C ¯ e L ¯ p . 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:
L ˜ = ( V e ) 1 L V e = ( V e ) 1 V ˙ e + L ˜ * .
After substituting the value of L ˜ * from Equation (11), Equation (17) becomes
L ˜ = ( V e ) 1 V ˙ e + Θ ˜ e + [ ( 1 i = 1 N t υ i ) α = 1 N s γ ˙ α S ˜ α + F s p { i = 1 N t υ i γ ˙ i S ˜ i } ( F s p ) 1 ] .
The finite Green strain tensor in reference configuration (third intermediate state with reference to current configuration) can be represented as:
E ˜ e = 1 2 ( C ˜ e I ) , E ˜ ˙ e = 1 2 C ˜ ˙ e .
Here, C ˜ e is the right Cauchy–Green strain tensor, and I is a second-order identity tensor. The 2nd Piola–Kirchhoff (PK2) stress, T e , 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:
T s e = C ˜ s : E ˜ e ,
T t e = C ˜ t : E ˜ e ,
where T s e and T t e are the PK2 stress tensors, and C ˜ s and C ˜ t 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:
T e = C ˜ e : E ˜ e .
Here, T e is an equivalent PK2 stress, and C ˜ e is the equivalent elasticity tensor. An equivalent form of elasticity tensor can be expressed through volume averaging technique as:
C ˜ e = ( 1 i = 1 N t υ i ) C ˜ s + i = 1 N t υ i C ˜ t .
Similarly, Cauchy stress can also be approximated on the basis of volume averaging technique as:
T = ( 1 i = 1 N t υ i ) T s + i = 1 N t υ i T t .

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:
E ˙ d = σ : F ˙ + ρ 0 θ Π ˙ m ρ 0 ϵ ˙ q ,
where σ is PK1 (1 st Piola–Kirchhoff) stress, F ˙ is the rate of change of F , ρ 0 is the density, θ is an absolute temperature, Π ˙ m 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:
E ˙ d = ( 1 i = 1 N t υ i ) α = 1 N s ( τ α + Φ α ρ 0 ϵ ζ Ψ α ) γ ˙ α + i = 1 N t υ ˙ i ( τ i + Φ i ρ 0 ϵ ζ Ψ i ) γ i ( θ ) Π q .
where τ α and τ i are the resolved shear stresses on α -slip and i-twin systems, respectively; Φ α and Φ i are the thermal equivalents of resolved shear stresses τ α and τ i , respectively; ϵ is the internal energy density; ζ is the crystal defect microstrain parameter; Ψ α and Ψ i 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 Π q 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 ( F e , θ , ζ , υ ) = E h m ( F e ) + E h t ( θ ) + E h d ( ζ ) + E h s ( υ ) ,
where E h m , E h t , E h d , and E h s 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 e , absolute temperature θ , crystal defect microstrain parameter ζ , and twin martensite volume fraction υ . In consideration of this, HFE can be derived as:
E h = 1 ρ 0 { ( F e F p ) : ( C ˜ e : E ˜ e E ˜ e ) } ( F p ) T : [ F e + ( F e ) T ] 1 + θ [ h e ln ( θ θ r ) + ( h e q Π m , 0 e ) ] + 1 2 ρ 0 φ G e ζ 2 + χ ρ 0 l 0 i = 1 N t υ i ( 1 i = 1 N t υ i ) .
In Equation (28), h e is an equivalent specific heat; θ r is the reference temperature; Π m , 0 e is the initial entropy density; φ is a dimensionless dislocation interaction parameter, which incorporates the effects of dislocations’ mobility and their interactions in plasticity; G e is the equivalent modulus of rigidity; l 0 is the length scale parameter; and χ is the interfacial energy per unit area.

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.
G α = ( 1 i = 1 N t υ i ) ( τ α + Φ α φ G e ζ Ψ α ) .
P i = υ i ( τ i + Φ i φ G e ζ Ψ i ) .

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:
γ ˙ α = γ ˙ 0 | τ α s r α | 1 m sign ( τ α ) .
In Equation (31), γ ˙ 0 is the initial shear strain rate; τ α is the shear stress on the α -slip system; s r α 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,26,27,57]. Accordingly, the rate of change of the twinned-martensite volume fraction of i-twin system is given as:
υ ˙ i = γ ˙ 0 γ i ( τ i s r i ) 1 / m ,
where γ ˙ 0 is the initial shear strain rate, γ i is the shear strain of i-twin system, and s r i 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 ˙ s α , and latent, s ˙ l α , hardening contributions of α -slip systems as:
s ˙ r α = s ˙ s α + s ˙ l α = α = 1 N s h 0 α [ 1 ( s r α s r , 0 α s r , S α s r , 0 α ) ] | γ ˙ α | + κ = 1 N s h α κ | γ ˙ κ | .
In Equation (33), h 0 α and s r , 0 α are the initial values of hardening rate and strength of slip system, respectively; s r , S α is the saturation value of slip strength; h α κ is an array of latent hardening values; and γ ˙ κ is the shear strain rate of κ -slip system, where κ denotes a number of slip system except α ( κ = j , j = 1 , , i 1 , i + 1 , , N s ) . The hardening induced by mechanical twinning is incorporated in the model through Equation (34) as:
s ˙ t i = h n c i ( i = 1 N t υ i ) d μ 1 = 0 i γ i υ ˙ μ 1 + h c p i ( i = 1 N t υ i ) μ 2 = 0 i γ i υ ˙ μ 2 .
Here, h n c i and h c p i are the initial hardening rates of non-coplanar and coplanar twin systems; μ 1 ( μ 1 i ) and μ 2 ( μ 2 i ) represent the number of non-coplanar and coplanar twin systems, respectively; υ ˙ μ 1 and υ ˙ μ 2 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:
ζ ˙ s = 1 ω G e N s α = 1 N s s ˙ r α = 1 ω G e N s α = 1 N s κ = 1 N s h α κ | γ ˙ κ | .
ζ ˙ t = 1 ω G e N t i = 1 N t [ { h n c i ( i = 1 N t υ i ) d μ 1 = 0 i υ ˙ μ 1 + h c p i ( i = 1 N t υ i ) μ 2 = 0 i υ ˙ μ 2 } γ i ] .
The stress-like parameters that are part of the driving potential equations for slip and twinning, Equations (29) and (30), are estimated as:
Ψ α = ( 1 i = 1 N t υ i ) 1 ω G e N s κ = 1 N s h α κ .
Ψ i = ( i = 1 N t υ ˙ i ) 1 ω G e N t i = 1 N t [ h n c i ( i = 1 N t υ i ) d μ 1 = 0 i υ ˙ μ 1 + h c p i ( i = 1 N t υ i ) μ 2 = 0 i υ ˙ μ 2 ] .

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:
{ T e , s r α , s t i , R e , υ i } ,
where T e represents second Piola–Kirchhoff stress tensor, { s r α | α = 1 , , N s l } denotes slip resistance of α -slip system, { s t i | i = 1 , , N t w } shows twin resistance of i-twin system, R e is the rigid body rotation tensor, and { υ i | i = 1 , , N t w } 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 t n and t n + 1 . In this, a relation of t n + 1 = t n + Δ t is used, where t n and t n + 1 represent time at the start and end of the time step, respectively. Afterwards, the magnitudes of all variables are evaluated at t n and t n + 1 , and denoted with the subscripts n and n + 1 , respectively. The numerical time integration scheme is based on the following assumptions:
  • The total deformation gradient F n + 1 and velocity gradient L n + 1 are given.
  • The values of variables T n e , s r , n α , s t , n i , R n e , and υ n i at time t n are known.
  • The initial values of time-independent slip ( m 0 α , n 0 α ) and twin ( m 0 i , n 0 i ) systems’ vectors, elasticity tensors C 0 s l and C 0 t w , initial crystallographic orientation of crystal Q 0 , initial kinetic flow rule (m and γ ˙ 0 ), and hardening parameters ( h 0 α , s r , 0 α , s r , S 0 α , γ ˙ S 0 ) are used as an input.
The output of the numerical integration scheme provides updated values of variables as: T n + 1 e , s r , n + 1 α , s t , n + 1 i , R n + 1 e , and υ n + 1 i at time t n + 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:
Δ E ˜ e = Δ t ( D ˜ D ˜ * ) .
Integration of Equation (40) using backward Euler scheme provides an updated value of Green strain tensor, as follows:
E ˜ n + 1 e = E ˜ n e + Δ E ˜ n e = E ˜ n e + Δ t ( D ˜ n + 1 D ˜ n + 1 * ) .
where D ˜ n + 1 and D ˜ n + 1 * are the symmetric parts of velocity gradients L ˜ n + 1 and L ˜ n + 1 * , respectively. The incremental values of these parameters can be illustrated as:
D ˜ n + 1 = 1 2 [ L ˜ n + 1 + L ˜ n + 1 T ] = sym ( L ˜ n + 1 ) .
D ˜ n + 1 * = sym ( C ˜ n + 1 e Θ ˜ n + 1 e ) + R n + 1 e D ˘ n + 1 p ( R n + 1 e ) T ,
where an updated value of right Cauchy–Green deformation tensor in third intermediate configuration can be defined as:
C ˜ n + 1 e = V n + 1 e ( V n + 1 e ) T ,
Furthermore, an incremental magnitude of spin tensor Θ ˜ n + 1 e is represented as:
Θ ˜ n + 1 e = R ˙ n + 1 e ( R n + 1 e ) T = Δ t R n + 1 e ( R n + 1 e ) T ,
In addition, D ˘ n + 1 p is an updated symmetric part of ( C ˘ n + 1 e L ˘ n + 1 p ) . Here, C ˘ n + 1 e is an incremental form of right Cauchy–Green deformation tensor in second intermediate configuration, which can be expressed as:
C ˘ n + 1 e = ( F n + 1 e ) T ( F n + 1 e ) .
The incremental value of the plastic velocity gradient in second intermediate configuration can be defined as:
L ˘ n + 1 p = ( 1 i = 1 N t w υ n + 1 i ) α = 1 N s l Δ t γ ˙ n + 1 α S ¯ n + 1 α + F s , n + 1 p { i = 1 N t w υ n + 1 i Δ t γ ˙ n + 1 i S ˘ n + 1 i } ( F s , n + 1 p ) 1 .
where υ n + 1 i represents an incremental change in the volume fraction of i t h twin system, γ ˙ n + 1 α is an updated value of the shear strain rate of α -slip system, γ ˙ n + 1 i is an updated value of shear strain induced by i-twin system, S ¯ n + 1 α and S ˇ n + 1 i are the updated magnitudes of Schmid and twin orientation tensors in second intermediate configuration, and F s , n + 1 p is an incremental change in the value of the plastic deformation gradient. Substitution of D ˘ n + 1 p in Equation (43) gives the following:
D ˜ n + 1 * = sym ( C ˜ n + 1 e Θ ˜ n + 1 e ) + ( R n + 1 e ) sym [ C ˘ n + 1 e { ( 1 i = 1 N t w υ n + 1 i ) α = 1 N s l Δ t γ ˙ n + 1 α S ¯ α + F s , n + 1 p { i = 1 N t w υ n + 1 i Δ t γ ˙ n + 1 i S ˘ i } ( F s , n + 1 p ) 1 } ] ( R n + 1 e ) T .
Considering the effects of rigid body rotation tensor R n + 1 e on Schmid and twin orientation tensors, Equation (48) can also be written as:
D ˜ n + 1 * = sym ( C ˜ n + 1 e Θ ˜ n + 1 e ) + sym [ C ¯ n + 1 e { ( 1 i = 1 N t w υ n + 1 i ) α = 1 N s l Δ t γ ˙ n + 1 α ( R n + 1 e ) ( m ¯ α n ¯ α ) ( R n + 1 e ) T + F s , n + 1 p { i = 1 N t w υ n + 1 i Δ t γ ˙ n + 1 i ( R n + 1 e ) ( m ˘ i n ˘ i ) ( R n + 1 e ) T } ( F s , n + 1 p ) 1 } ] .
The forward mapping approach for Schmid and twin orientation tensors can be implemented into Equation (49) as:
{ R n + 1 e . ( m ¯ α n ¯ α ) } . ( R n + 1 e ) T = { R n + 1 e . ( Q 0 . m 0 α Q 0 . n 0 α ) } . ( R n + 1 e ) T = { R n + 1 e . ( Q 0 . m 0 α ) } { R n + 1 e . ( Q 0 . n 0 α ) } .
where Q 0 is the initial rotation matrix that is used to transform crystal coordinates to sample coordinate systems through Euler angles, and m 0 α and n 0 α 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 + 1 = R n + 1 e . Q 0 . The rotation matrix can also be used to transform Schmid orientation vectors from reference to first intermediate configuration as follows: m ¯ α = Q 0 . m 0 α and n ¯ α = Q 0 . n 0 α . In this condition, Equation (50) can be rewritten as:
{ ( R n + 1 e ) . ( m ¯ α n ¯ α ) } ( R n + 1 e ) T = ( Q n + 1 . m 0 α ) ( Q n + 1 . n 0 α ) = m ˜ n + 1 α n ˜ n + 1 α .
In Equation (51), m ˜ n + 1 α and n ˜ n + 1 α 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 + 1 is used to transform Schmid orientation vectors from reference to third intermediate configuration. Similarly, twin orientation vectors are given as:
{ R n + 1 e . ( m ˘ i n ˘ i ) } . ( R n + 1 e ) T = m ˜ n + 1 i n ˜ n + 1 i .
Substitution of Equations (51) and (52) in (49) provides
D ˜ n + 1 * = sym ( C ˜ n + 1 e . Θ ˜ n + 1 e ) + sym [ C ˘ n + 1 e . { ( 1 i = 1 N t w υ n + 1 i ) α = 1 N s l Δ t γ ˙ n + 1 α S ˜ n + 1 α + F s , n + 1 p . { i = 1 N t w υ n + 1 i Δ t γ ˙ n + 1 i S ˜ n + 1 i } ( F s , n + 1 p ) 1 } ] .
Here, S ˜ n + 1 α = m ˜ n + 1 α n ˜ n + 1 α and S ˜ n + 1 i = m ˜ n + 1 i n ˜ n + 1 i 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 t n + 1 can be defined by Equations (54) and (55), respectively, as:
γ n + 1 α = γ n α + Δ γ α = γ n α + γ 0 | τ n + 1 α s r , n + 1 α | 1 m sin ( τ n + 1 α ) .
υ n + 1 i = υ n i + Δ υ i = υ n i + γ 0 γ n + 1 i | τ n + 1 i s t , n + 1 i | 1 m .
As discussed previously, the rigid body rotation tensor R n + 1 e updates the crystal orientation (Euler angles) matrix Q n + 1 , which can be used to transform the fourth order elasticity tensors for slip C ˜ n + 1 s and twinned C ˜ n + 1 t regions to the third intermediate configuration as follows:
C ˜ n + 1 s = ( Q n + 1 Q n + 1 ) : C ˜ 0 s : ( Q n + 1 Q n + 1 ) T .
C ˜ n + 1 t = ( Q n + 1 Q n + 1 ) : C ˜ 0 t : ( Q n + 1 Q n + 1 ) T .
Any one of the elasticity tensors can be obtained in terms of another by using coordinate transformation rule defined as:
C ˜ i j k l t = C ˜ p q r s s P i p P j q P k r P l s = ( P P ) : C ˜ s : ( P P ) T .
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:
P j k = 2 n j n k δ j k , j , k = 1 , 2 , 3
where n is the area normal vector of the twin plane, and δ is the Kroneker delta. An equivalent elasticity tensor can be calculated as:
C ˜ n + 1 e = ( 1 i = 1 N t υ n + 1 i ) C ˜ n + 1 s + i = 1 N t υ n + 1 i C ˜ n + 1 t .
Furthermore, an updated form of equivalent second Piola–Kirchhoff stress tensor can be estimated as:
T n + 1 e = C ˜ n + 1 e : E ˜ n + 1 e .
An evolution equation for rigid body rotation tensor R e is integrated through the exponential map discussed by [58] as follows:
R n + 1 e = exp ( Δ t Θ ˜ n + 1 e ) . R n e ,
where an updated value of the spin of lattice Θ ˜ n + 1 e can be estimated through the skew-symmetric component W ˜ n + 1 * of the velocity gradient L ˜ n + 1 * as:
W ˜ n + 1 * = skew ( C ˜ n + 1 e Θ ˜ n + 1 e ) + α = 1 N s Δ t γ ˙ n + 1 α skew ( C ˜ n + 1 e S ˜ n + 1 α ) .
By using backward mapping, the skew-symmetric component of velocity gradient L ˜ n + 1 can be estimated as:
W ˜ n + 1 = { ( V n + 1 e ) T W n + 1 } V n + 1 e = ( V e ) T skew { ( V n + 1 e ) T { Δ t V n + 1 e } } V e + W ˜ n + 1 * .
Here, W n + 1 is the updated skew-symmetric component of L n + 1 . Substitution of W ˜ n + 1 * from Equation (63) in (64) provides
skew ( C ˜ n + 1 e Θ ˜ n + 1 e ) = W ˜ n + 1 ( V e ) T skew { ( V n + 1 e ) T Δ t V n + 1 e } V e α = 1 N s Δ t γ ˙ n + 1 α skew ( C ˜ n + 1 e S ˜ n + 1 α ) .
For small elastic strain problems, the value of the right Cauchy–Green deformation tensor C ˜ n + 1 e is typically small. In this case, Equation (65) can be written as:
skew ( Θ ˜ n + 1 e ) = W ˜ n + 1 ( V e ) T skew { ( V n + 1 e ) T Δ t V n + 1 e } V e α = 1 N s Δ t γ ˙ n + 1 α skew ( S ˜ n + 1 α ) .
An incremental value of the slip resistance of α -slip system can be evaluated using a backward Euler scheme as follows:
s r , n + 1 α = s r , n α + Δ s r α = s r , n α + α = 1 N s h 0 α [ 1 ( s r , n + 1 α s r , 0 α s r , S n + 1 α s r , 0 α ) ] | ( Δ t ) γ ˙ n + 1 α | .
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 [ 1 ( s r , n + 1 s r , 0 s r , S n + 1 s r , 0 ) ] α = 1 N s | ( Δ t ) γ ˙ n + 1 α | ,
In Equation (68), s r , S n + 1 is an incremental saturation value of slip resistance, which can be calculated as:
s r , S n + 1 = s r , S n + Δ s r , S = s r , S n + s r , S 0 [ β N s | γ n + 1 β | γ S 0 ] k θ ( G e ) n + 1 b 3 Z + s r , p ( λ = 0 i υ n + 1 λ ) a 1 ,
where s r , S 0 is the initial value of saturation slip resistance, Δ γ S 0 is an incremental initial value of the slip system shear strain at the initial value of saturation slip resistance, a is the material parameter, s r , p is the material slip-hardening parameter, λ ( λ i ) represents the number of non-coplanar twin systems to slip system, and a 1 is a material parameter. In addition, the twin resistance can be estimated using Equation (34) as:
s t , n + 1 i = s t , n i + Δ s t i = s t , n i + [ h n c i ( i = 1 N t υ n + 1 i ) d μ 1 = 0 i γ n + 1 i Δ t υ ˙ n + 1 μ 1 + h c p i ( i = 1 N t υ n + 1 i ) μ 2 = 0 i γ n + 1 i Δ t υ ˙ n + 1 μ 2 ] ,
The resistance of all twin systems are assumed to be identical. Therefore, Equation (70) can be expressed as:
s t , n + 1 = s t , n + [ h n c ( i = 1 N t υ n + 1 i ) d μ 1 = 0 i γ n + 1 i Δ t υ ˙ n + 1 μ 1 + h c p ( i = 1 N t υ n + 1 i ) μ 2 = 0 i γ n + 1 i Δ t υ ˙ n + 1 μ 2 ] .
The updated driving potentials for slip and twin mode of deformations can be estimated using Equations (29) and (30), respectively as:
G s , n + 1 α = G s , n α + Δ G s α = G s , n α + ( 1 i = 1 N t υ n + 1 i ) ( τ n + 1 α + Φ α φ G n + 1 e ζ n + 1 Ψ n + 1 α ) .
G t , n + 1 i = G t , n i + Δ G t i = G t , n i + υ n + 1 i ( τ n + 1 i + Φ i φ G n + 1 e ζ n + 1 Ψ n + 1 i ) .
The parameters Φ α and Φ i 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 n + 1 e can be estimated as:
G n + 1 e = ( 1 i = 1 N t υ n + 1 i ) G s + i = 1 N t υ n + 1 i G t .
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:
ζ s , n + 1 = ζ s , n + Δ ζ s = ζ s , n + 1 ω G n + 1 e N s α = 1 N s β = 1 N s h n + 1 α β | γ n + 1 β | .
ζ t , n + 1 = ζ t , n + Δ ζ t = ζ t , n + 1 ω G n + 1 e N t [ { h n c ( i = 1 N t υ n + 1 i ) d μ 1 = 0 i ( Δ t ) υ ˙ n + 1 μ 1 + h c p ( i = 1 N t υ n + 1 i ) μ 2 = 0 i ( Δ t ) υ ˙ n + 1 μ 2 } i = 1 N t γ n + 1 i ] .
The updated stress-like functions Ψ n + 1 α and Ψ n + 1 i for slip and twin modes of deformation can be calculated using Equations (37) and (38), respectively, as:
Ψ n + 1 α = ( 1 i = 1 N t υ n + 1 i ) 1 ω G n + 1 e N s β = 1 N s h n + 1 α β .
Ψ n + 1 i = ( i = 1 N t υ n + 1 i ) 1 ω G n + 1 e N t [ h n c ( i = 1 N t υ n + 1 i ) d μ 1 = 0 i υ n + 1 μ 1 + h c p ( i = 1 N t υ n + 1 i ) μ 2 = 0 i υ n + 1 μ 2 ] .
Furthermore, the updated Cauchy stress tensor T n + 1 can be calculated using second Piola–Kirchhoff stress tensor, T n + 1 e , as:
T n + 1 = F n + 1 e { det F n + 1 e . T n + 1 e } ( F n + 1 e ) T .

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 n + 1 e , R n + 1 e , s r , n + 1 , s t , n + 1 , and υ n + 1 i , 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 i v | i = 1 , 2 , , 5 } , where
p 1 v = Δ T e , p 2 v = Δ s r , p 3 v = Δ υ i , p 4 v = Δ R e , p 5 v = Δ s t .
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 i | i = 1 , 2 , , 5 } . The components of R i are given through Equations (81)–(85) as:
R 1 = R ^ 1 ( T n + 1 e , s r , n + 1 α , υ n + 1 i , R n + 1 e , s t , n + 1 i ) = ( C ˜ n + 1 e ) 1 : T n + 1 e E ˜ n + 1 e = ( C ˜ n + 1 e ) 1 : T n + 1 e E ˜ n e Δ t ( D ˜ n + 1 D ˜ n + 1 * ) = 0 .
R 2 = R ^ 2 ( T n + 1 e , s r , n + 1 α , υ n + 1 i , R n + 1 e , s t , n + 1 i ) = s r , n + 1 α s r , n α α = 1 N s h 0 α [ 1 ( s r , n + 1 α s r , 0 α s r , S n + 1 α s r , 0 α ) ] | ( Δ t ) γ ˙ n + 1 α | = 0 .
R 3 = R ^ 3 ( T n + 1 e , s r , n + 1 α , υ n + 1 i , R n + 1 e , s t , n + 1 i ) = υ n + 1 i υ n i γ 0 γ n + 1 i | τ n + 1 i s t , n + 1 i | 1 / m = 0 .
R 4 = R ^ 4 ( T n + 1 e , s r , n + 1 α , υ n + 1 i , R n + 1 e , s t , n + 1 i ) = R n + 1 e exp { Δ t Θ ˜ n + 1 e } . R n e = 0 .
R 5 = R ^ 5 ( T n + 1 e , s r , n + 1 α , υ n + 1 i , R n + 1 e , s t , n + 1 i ) = s t , n + 1 i s t , n i [ h n c i ( i = 1 N t υ n + 1 i ) d μ 1 = 0 i γ n + 1 i ( Δ t ) υ ˙ n + 1 μ 1 + h c p i ( i = 1 N t υ n + 1 i ) μ 2 = 0 i γ n + 1 i ( Δ t ) υ ˙ n + 1 μ 2 ] = 0 .
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 n + 1 e by assuming the best possible values of the other primary variables ( R n + 1 e , s r , n + 1 α , s t , n + 1 i , υ n + 1 i ) . 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 r , n + 1 α from Equation (82). This considers an updated value of T n + 1 e and the estimated values of R n + 1 e , s t , n + 1 i , and υ n + 1 i . Finally, updated values of the twinned martensite volume fraction υ n + 1 i , rigid body rotation tensor R n + 1 e , and twin resistance s t , n + 1 i 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 two-level iterative procedure is based on the variation of L 2 -norm for T n + 1 e , and s r , n + 1 α . If the L 2 -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:
R trial 2 < Tol .
In the current Newton–Raphson iterative scheme, iterations are carried out unless and until the variation in the L 2 -norm of the residuals for T n + 1 e and s r , n + 1 α satisfy the conditions in Equations (87) and (88), respectively, as follows:
Δ T e 2 < ( 10 4 ) | T e | trial
Δ s r α 2 < ( 10 4 ) | s r α | trial
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:
K = Δ A max Δ A x ,
where Δ A max is the maximum value of Δ A over all the crystals, all the integration points are in finite element mesh, and Δ A x 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 1 A 2 } , where
A 1 = T e , A 2 = s r α .
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 Δ t n + 1 = ( Δ t ) n / K . In this condition, the new time step size Δ t n + 1 is more or less identical to the previous Δ t n . (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.

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 e , s r α , υ i , s t i , R e , } to { T e , s r α , R e } , by omitting twin resistance s t i and twinned martensite volume fraction υ i . 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 4 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 1 c , e 2 c , e 3 c ) and ( e 1 s , e 2 s , e 3 s ) axes. The crystal direction [100], corresponding to Euler angles ( ϕ 1 1 , ϕ 2 1 , ϕ 3 1 ) , is equal to ( 90 o , 0 , 90 o ) (see Figure 4a). Likewise, [110] and [111], corresponding to ( ϕ 1 2 , ϕ 2 2 , ϕ 3 2 ) and ( ϕ 1 3 , ϕ 2 3 , ϕ 3 3 ) , are equivalent to ( 45 o , 0 , 90 o ) and ( 45 o , 35 . 26 o , 54 . 74 o ) , respectively, as represented in Figure 4b,c.

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.

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 1 s axis, as shown in Figure 4. In uniaxial compression, displacement of 35.0 × 10 2 mm is applied along e 1 s vector on a surface with e 1 s 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.

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 Figure 6 and Figure 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 α and slip direction, and m k α vectors, where α and k are equal to 1, as shown in Table A1. In addition, vectors n 1 and m 2 1 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.
It is evident from Figure 6 and Figure 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 Figure 6 and Figure 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 Figure 8 and Figure 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.
It is noted in Figure 8 and Figure 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.

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 Figure 10 and Figure 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.
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.

5.2.4. Twin Volume Fraction

The magnitude of shear strain alone does not present a complete quantitative estimation of twinning. A complete estimation requires the magnitude of the volume of twin in the overall plasticity. In view of this requirement, the twin volume fraction is estimated during the course of deformation, as represented through Figure 12. In conjunction with twin systems’ activity, the twin volume fraction represents the ratio of active to total number of twin planes. A significant variation of twin volume is observed under tension and compression in all steels and crystallographic directions. As in the case of shear strain, the twin volume fraction becomes nearly steady for all crystallographic directions of TWIP 1 and 2, but not for TWIP 3, after a certain equivalent strain. The magnitude of twin volume fraction for TWIP steels is presented in Table 6. As is evident, the largest magnitude of twin volume is observed in the [100] crystallographic direction under tension and compression for all three TWIP steels. Overall, the twin volume fraction is not exceeded by 0.333 (33.33 %) in all directions and steels. This shows that the contribution of twinning in the plasticity of TWIP steels 1, 2, and 3 is limited to 33.33 % for a single crystal.

6. Conclusions

A numerical scheme is developed and implemented for modeling the elastic-plastic deformation behavior of twinning-induced plasticity steel, in which mechanical twinning contributes significantly to plasticity along with a crystallographic slip. Initially, a constitutive formulation of equations, reported in earlier work, is briefly discussed. Afterward, a numerical integration procedure is established by identifying primary variables, discretizing constitutive equations in the time domain, developing an iterative scheme, and introducing a time sub-stepping algorithm. A numerical integration scheme is then incorporated in finite element software ABAQUS through a user-defined material subroutine. Finite element models of single-crystal and polycrystal material points are developed and simulation results are compared with the published experimental observations. They are in close accord with the maximum error of 16.15% in TWIP steels for equivalent stress. However, the error becomes higher beyond 0.3 strain. This puts limitations on the current model to be implemented, with more accuracy, at high strain deformation. This must be further explored in future work. In addition, further simulations are executed to quantify slip and twin systems’ activity, deformation behavior, shear strain pattern, and twin volume fraction of three TWIP steels subjected to uniaxial tension and compression. The slip and twin systems’ activity shows that the slip contribution is higher for all crystallographic directions in all steels; however, the highest level of twin activity (number of active twin systems) is observed in the [100] direction. In addition, a higher magnitude of stress at 0.45 equivalent strain is observed in tension (1150 to 1300 MPa) than compression (780 to 850 MPa) in all polycrystal TWIP steels. It may indicate a more prominent role of slip and twin systems’ reorientation and interactions in tension. Moreover, the twin shear strain becomes nearly constant for all crystallographic directions of TWIP 1 and 2 after 0.2 equivalent strain, but not for TWIP 3. It is also found that the fraction of twin volume is not surpassed by 33.33% for all directions and TWIP steels. The development and implementation of a fully implicit numerical integration scheme in modeling twinning-induced plasticity provide an effective platform to estimate the deformation behavior of TWIP steels. The current work can be enhanced to incorporate martensitic phase transformation and damage criterion in a coupled slip, twinning, and transformation-induced plasticity model.

Author Contributions

R.K.: conceptualization, methodology, investigation, and writing—original draft preparation; T.P.: supervision, and writing—review and editing; A.A.: validation and visualization; S.Z.Q.: supervision, formal analysis; S.M.: software and visualization. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors acknowledge the support of Imam Mohammad Ibn Saud Islamic University for this research.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Slip and Twin Systems of Austenite Crystal

Face-centered cubic (FCC) austenite crystal’s slip and twin systems are represented in terms of Miller indices and unit vectors in Table A1 and Table A2, respectively.
Table A1. Slip systems of FCC crystal.
Table A1. Slip systems of FCC crystal.
Slip Plane Normal n α Slip Direction m k α
α = 1–4, k = 1 , 2 , 3
MillerUnit Vector MillerUnit Vector
n 1 (111)( 1 3 , 1 3 , 1 3 ) m 1 1 [ 01 1 ¯ ] (0, 1 2 , 1 2 )
m 2 1 [ 10 1 ¯ ] ( 1 2 , 0, 1 2 )
m 3 1 [ 1 1 ¯ 0 ] ( 1 2 , 1 2 , 0)
n 2 ( 1 ¯ 11 )( 1 3 , 1 3 , 1 3 ) m 1 2 [ 01 1 ¯ ] (0, 1 2 , 1 2 )
m 2 2 [ 101 ] ( 1 2 , 0, 1 2 )
m 3 2 [ 110 ] ( 1 2 , 1 2 , 0)
n 3 ( 1 ¯ 1 ¯ 1 )( 1 3 , 1 3 , 1 3 ) m 1 3 [ 011 ] (0, 1 2 , 1 2 )
m 2 3 [ 101 ] ( 1 2 , 0, 1 2 )
m 3 3 [ 1 1 ¯ 0 ] ( 1 2 , 1 2 , 0)
n 4 ( 1 1 ¯ 1 )( 1 3 , 1 3 , 1 3 ) m 1 4 [ 011 ] (0, 1 2 , 1 2 )
m 2 4 [ 10 1 ¯ ] ( 1 2 , 0, 1 2 )
m 3 4 [ 110 ] ( 1 2 , 1 2 , 0)
Table A2. Twin systems of FCC crystal.
Table A2. Twin systems of FCC crystal.
Twin Plane Normal n i Twin Direction m k i
i = 1–4, k = 1 , 2 , 3
MillerUnit Vector MillerUnit Vector
n 1 (111)( 1 3 , 1 3 , 1 3 ) m 1 1 [ 2 ¯ 11 ] ( 2 6 , 1 6 , 1 6 )
m 2 1 [ 1 2 ¯ 1 ] ( 1 6 , 2 6 , 1 6 )
m 3 1 [ 11 2 ¯ ] ( 1 6 , 1 6 , 2 6 )
n 2 ( 1 ¯ 11 )( 1 3 , 1 3 , 1 3 ) m 1 2 [ 211 ] ( 2 6 , 1 6 , 1 6 )
m 2 2 [ 1 ¯ 2 ¯ 1 ] ( 1 6 , 2 6 , 1 6 )
m 3 2 [ 1 ¯ 1 2 ¯ ] ( 1 6 , 1 6 , 2 6 )
n 3 ( 1 ¯ 1 ¯ 1 )( 1 3 , 1 3 , 1 3 ) m 1 3 [ 2 1 ¯ 1 ] ( 2 6 , 1 6 , 1 6 )
m 2 3 [ 1 ¯ 21 ] ( 1 6 , 2 6 , 1 6 )
m 3 3 [ 1 ¯ 1 ¯ 2 ¯ ] ( 1 6 , 1 6 , 2 6 )
n 4 ( 1 1 ¯ 1 )( 1 3 , 1 3 , 1 3 ) m 1 4 [ 2 ¯ 1 ¯ 1 ] ( 2 6 , 1 6 , 1 6 )
m 2 4 [ 121 ] ( 1 6 , 2 6 , 1 6 )
m 3 4 [ 1 1 ¯ 2 ¯ ] ( 1 6 , 1 6 , 2 6 )

References

  1. Xu, M.; David, J.; Kim, S. The Fourth Industrial Revolution: Opportunities and Challenges. Int. J. Financ. Res. 2018, 9, 90. [Google Scholar] [CrossRef] [Green Version]
  2. Demeri, M.Y. Advanced High-Strength Steels: Science, Technology, and Applications; ASM International: Almere, The Netherlands, 2013. [Google Scholar]
  3. Kuziak, R.; Kawalla, R.; Waengler, S. Advanced High Strength Steels for Automotive Industry. Arch. Civ. Mech. Eng. 2008, 8, 103–117. [Google Scholar] [CrossRef]
  4. Tamarelli, C.M. The Evolving Use of Advanced High Strength Steels for Automotive Applications; Technology Report; University of Michigan, USA, Materials Science and Engineering Department: Ann Arbor, MI, USA, 2011. [Google Scholar]
  5. Al-Abri, O.S.; Pervez, T.; Qamar, S.Z.; Khan, R. On the performance analysis of AHSS with an application to SET technology—FEM simulations and experimental measurements. Thin-Walled Struct. 2016, 101, 58–74. [Google Scholar] [CrossRef]
  6. Cooman, B.C.D.; Estrin, Y.; Kim, S.K. Twinning-induced plasticity (TWIP) steels. Acta Mater. 2018, 142, 283–362. [Google Scholar] [CrossRef]
  7. Fan, J.; Qiao, J.W.; Wang, Z.H.; Rao, W.; Kang, G.Z. Twinning-induced plasticity (TWIP) and work hardening in Ti-based metallic glass matrix composites. Sci. Rep. 2017, 7, 1877. [Google Scholar] [CrossRef] [Green Version]
  8. Grässel, O.; Kröger, L.; Frommeyer, G.; Meyer, L.W. High strength Fe-Mn-Al-Si TRIP/TWIP steels development: Properties & Application. Int. J. Plast. 2000, 16, 1391–1409. [Google Scholar]
  9. Yang, G.; Kim, J.-K. An Overview of High Yield Strength Twinning-Induced Plasticity Steels. Metals 2021, 11, 124. [Google Scholar] [CrossRef]
  10. Li, D.; Qian, L.; Wei, C.; Liu, S.; Zhang, F.; Meng, J. The role of Mn on twinning behavior and tensile properties of coarse- and fine-grained Fe–Mn–C twinning-induced plasticity steels. Mater. Sci. Eng. A 2020, 789, 139586. [Google Scholar] [CrossRef]
  11. Pierce, D.; Jime´nez, J.; Bentley, J.; Raabe, D.; Oskay, C.; Wittig, J. The influence of manganese content on the stacking fault and austenite-ϵ martensite interfacial energies in Fe-Mn-(Al-Si steels investigated by experiment and theory. Acta Mater. 2014, 68, 238–253. [Google Scholar] [CrossRef]
  12. Lee, Y.K.; Lee, S.J.; Han, J. Critical assessment 19: Stacking fault energies of austenitic steels. Mater. Sci. Technol. 2016, 32, 1–8. [Google Scholar] [CrossRef]
  13. Bouaziz, O.; Allain, S.; Scott, C.P.; Cugy, P.; Barbier, D. High manganese austenitic twinning induced plasticity steels: A review of the microstructure properties relationships. Curr. Opin. Solid State Mater. Sci. 2011, 15, 141–168. [Google Scholar] [CrossRef]
  14. Zhi, H.; Zhang, C.; Antonov, S.; Yu, H.; Guo, T.; Su, Y. Investigations of dislocation-type evolution and strain hardening during mechanical twinning in Fe-22Mn-0.6C twinning-induced plasticity steel. Acta Mater. 2020, 195, 371–382. [Google Scholar] [CrossRef]
  15. Idrissi, H.; Renard, K.; Ryel, L.; Schryvers, D.; Jacques, P.J. On the mechanism of twin formation in Fe–Mn–C TWIP steels. Acta Mater. 2010, 58, 2464–2476. [Google Scholar] [CrossRef]
  16. Wang, S.; Liu, Z.; Wang, G.; Liu, J.; Liang, G.; Li, Q. Effects of Twin-Dislocation and Twin-Twin Interactions on the Strain Hardening Behavior of TWIP Steels. J. Iron Steel Res. Int. 2010, 17, 70–74. [Google Scholar] [CrossRef]
  17. Pauly, S.; Gorantla, S.; Wang, G.; Kühn, U.; Eckert, J. Transformation-mediated ductility in CuZr-based bulk metallic glasses. Nat. Mater. 2010, 9, 473–477. [Google Scholar] [CrossRef]
  18. Bouaziz, O.; Zurob, H.; Chehab, B.; Embury, J.D.; Allain, S.; Huang, M. Effect of chemical composition on work hardening of Fe—Mn—C TWIP steels. Mater. Sci. Technol. 2011, 27, 707–709. [Google Scholar] [CrossRef]
  19. Peng, X.; Zhu, D.; Hu, Z.; Yi, W.; Liu, H.; Wang, M. Stacking fault energy and tensile deformation behavior of high-carbon twinning-induced plasticity steels: Effect of Cu addition. Mater. Des. 2013, 45, 518–523. [Google Scholar] [CrossRef]
  20. Martin, S.; Wolf, S.; Martin, U.; Krüger, L.; Rafaja, D. Deformation Mechanisms in Austenitic TRIP/TWIP Steel as a Function of Temperature. Metall. Mater. Trans. 2016, 47, 49–58. [Google Scholar] [CrossRef] [Green Version]
  21. Fei, L.; Weigang, Z.; Wenjiao, D. Stress-strain Response for Twinning-induced Plasticity Steel with Temperature. Procedia Eng. 2014, 81, 1330–1335. [Google Scholar]
  22. Allain, S.; Chateau, J.P.; Bouaziz, O.; Migot, S.; Guelton, N. Correlations between the calculated stacking fault energy and the plasticity mechanisms in Fe-Mn-C alloys. Mater. Sci. Eng. 2004, 387–389, 158–162. [Google Scholar] [CrossRef]
  23. Curtze, S.; Kuokkala, V.T.; Oikari, A.; Talonen, J.; Hänninen, H. Thermodynamic modeling of the stacking fault energy of austenitic steels. Acta Mater. 2011, 59, 1068–1076. [Google Scholar] [CrossRef]
  24. Idrissi, H.; Ryelandt, L.; Veron, M.; Schryvers, D.; Jacques, P. Is there a relationship between the stacking fault character and the activated mode of plasticity of fe-mn-based austenitic steels? Scr. Mater. 2009, 60, 941–944. [Google Scholar] [CrossRef]
  25. Tomé, C.N.; Lebensohn, R.A.; Kocks, U.F. A model for texture development dominated by deformation twinning: Application to zirconium alloys. Acta Metall. Mater. 1991, 39, 2667–2680. [Google Scholar] [CrossRef] [Green Version]
  26. Houtte, P.V. Simulation of the rolling and shear texture of brass by the Taylor theory adapted for mechanical twinning. Acta Metall. 1978, 26, 591–604. [Google Scholar] [CrossRef]
  27. Kalidindi, S.R. Incorporation of deformation twinning in crystal plasticity models. J. Mech. Phys. Solids 1998, 46, 267–290. [Google Scholar] [CrossRef]
  28. Taylor, G.I. Plastic strain in metals. J. Inst. Met. 1938, 62, 307–324. [Google Scholar]
  29. Chin, G.; Hosford, W.; Mendorf, D. Accommodation of constrained deformation in F. C. C. metals by slip and twinning. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 1969, 309, 433–456. [Google Scholar]
  30. Asaro, R.J.; Rice, J.R. Strain localization in ductile single crystals. J. Mech. Phys. Solids 1977, 25, 309–338. [Google Scholar] [CrossRef] [Green Version]
  31. Salem, A.A.; Kalidindi, S.R.; Doherty, R.D. Strain hardening of titanium: Role of deformation twinning. Acta Mater. 2003, 51, 4225–4237. [Google Scholar] [CrossRef]
  32. Salem, A.A.; Kalidindi, S.R.; Semiatin, S.L. Strain hardening due to deformation twinning in α-titanium: Constitutive relations and crystal-plasticity modeling. Acta Mater. 2005, 53, 3495–3502. [Google Scholar] [CrossRef]
  33. Salem, A.A.; Kalidindi, S.R.; Doherty, R.D.; Semiatin, S.L. Strain hardening due to deformation twinning in α-titanium: Mechanisms. Metall. Mater. Trans. A 2006, 37, 259–268. [Google Scholar] [CrossRef]
  34. Bhattacharya, K. Wedge-like microstructure in martensites. Acta Metall. Et Mater. 1991, 39, 2431–2444. [Google Scholar] [CrossRef]
  35. Wang, J.; Sehitoglu, H. Modelling of martensite slip and twinning in NiTiHf shape memory alloys. Philos. Mag. 2014, 94, 2297–2317. [Google Scholar] [CrossRef]
  36. Ostadrahimi, A.; Taheri-Behrooz, F. Analytical solution for twinning deformation effect of pre-strained shape memory effect beam-columns. J. Intell. Mater. Syst. Struct. 2019, 30, 2147–2165. [Google Scholar] [CrossRef]
  37. Allain, S.; Chateau, J.P.; Bouaziz, O. A physical model of the twinning-induced plasticity effect in a high manganese austenitic steel. Mater. Sci. Eng. 2004, 387–389, 143–147. [Google Scholar] [CrossRef]
  38. Bouaziz, O. Strain-hardening of twinning-induced plasticity steels. Scr. Mater. 2012, 66, 982–985. [Google Scholar] [CrossRef]
  39. Lee, M.G.; Kim, S.J.; Han, H.N. Crystal plasticity finite element modeling of mechanically induced martensitic transformation (MIMT) in metastable austenite. Int. J. Plast. 2010, 26, 688–710. [Google Scholar] [CrossRef]
  40. Khosravani, A.; Scott, J.; Miles, M.P.; Fullwood, D.; Adams, B.L.; Mishra, R.K. Twinning in magnesium alloy AZ31B under different strain paths at moderately elevated temperatures. Int. J. Plast. 2013, 45, 160–173. [Google Scholar] [CrossRef]
  41. Shiekhelsouk, M.N.; Favier, V.; Inal, K.; Cherkaoui, M. Modelling the behaviour of polycrystalline austenitic steel with twinning-induced plasticity effect. Int. J. Plast. 2009, 25, 105–133. [Google Scholar] [CrossRef]
  42. Khan, R.; Zahedi, F.I.; Siqqiui, A.K. Numerical modeling of twinning induced plasticity in austenite based advanced high strength steels. Procedia Manuf. 2016, 5, 772–786. [Google Scholar] [CrossRef] [Green Version]
  43. Yang, G.; Dayong, A.; Fengbo, H.; Liu, X.; Guozheng, K.; Xu, Z. Multiple-mechanism and microstructure-based crystal plasticity modeling for cyclic shear deformation of TRIP steel. Int. J. Mech. Sci. 2022, 222, 107269. [Google Scholar]
  44. Qayyum, F.; Guk, S.; Prahl, U. Studying the Damage Evolution and the Micro-Mechanical Response of X8CrMnNi16-6-6 TRIP Steel Matrix and 10% Zirconia Particle Composite Using a Calibrated Physics and Crystal-Plasticity-Based Numerical Simulation Model. Crystals 2021, 11, 759. [Google Scholar] [CrossRef]
  45. Qayyum, F.; Guk, S.; Kawalla, R.; Prahl, U. On Attempting to Create a Virtual Laboratory for Application-Oriented Microstructural Optimization of Multi-Phase Materials. Appl. Sci. 2021, 11, 1506. [Google Scholar] [CrossRef]
  46. Ching-Tun, P.; Mao, L.; Cheng, L.; Huijun, L. The determination of self hardening parameters of twinning induced plasticity steel via crystal plasticity modeling. J. Comput. Theor. Nanosci. 2015, 12, 2523–2530. [Google Scholar]
  47. Li, Y.; Zhu, L.; Liu, Y.; Wei, Y.; Wu, Y.; Tang, D.; Mi, Z. On the strain hardening and texture evolution in high manganese steels: Experiments and numerical investigation. J. Mech. Phys. Solids 2013, 61, 2588–2604. [Google Scholar] [CrossRef] [Green Version]
  48. Sun, C.; Wang, B.; Politis, D.J.; Wang, L.; Cai, Y.; Guo, X.; Guo, N. Prediction of earing in TWIP steel sheets based on coupled twinning crystal plasticity model. Int. J. Adv. Manuf. Technol. 2017, 89, 3037–3047. [Google Scholar] [CrossRef]
  49. Wong, S.L.; Madivala, M.; Prahl, U.; Roters, F.; Raabe, D. A crystal plasticity model for twinning and transformation-induced plasticity. Acta Mater. 2016, 118, 140–151. [Google Scholar] [CrossRef]
  50. Sun, C.; Guo, N.; Fu, M.; Wang, S. Modeling of slip, twinning and transformation induced plastic deformation for TWIP steel based on crystal plasticity. Int. J. Plast. 2016, 76, 186–212. [Google Scholar] [CrossRef]
  51. Khan, R.; Pervez, T.; Qamar, S.Z. Modeling and simulations of transformation and twinning induced plasticity in advanced high strength austenitic steels. Mech. Mater. 2016, 95, 83–101. [Google Scholar] [CrossRef]
  52. Manjunatha, M.; Alexander, S.; Su, L.W.; Franz, R.; Ulrich, P.; Wolfgang, B. Temperature dependent strain hardening and fracture behavior of TWIP steel. Int. J. Plast. 2018, 104, 80–103. [Google Scholar]
  53. Khan, R.; Alfozan, A. Modeling of twinning-induced plasticity using crystal plasticity and thermodynamic framework. Acta Mech. 2019, 230, 2687–2715. [Google Scholar] [CrossRef]
  54. Marin, E.B. On the Formulation of a Crystal Plasticity Model; Tech. Rep.; Sandia National Laboratories: Livermore, NM, USA, 2006.
  55. Lee, E.H. Elastic-plastic deformation at finite strains. J. Appl. Mech. 1969, 36, 1–6. [Google Scholar] [CrossRef]
  56. Kalidindi, S.R.; Bronkhorst, C.A.; Anand, L. Crystallographic texture evolution in bulk deformation processing of FCC metals. J. Mech. Phys. Solids 1992, 40, 537–569. [Google Scholar] [CrossRef]
  57. Turteltaub, S.; Suiker, A.S.J. A multiscale thermomechanical model for cubic to tetragonal martensitic phase transformations. Int. J. Solids Struct. 2006, 43, 4509–4545. [Google Scholar] [CrossRef] [Green Version]
  58. Simo, J.C.; Hughes, T.J.R. Computational Inelasticity; Springer: Berlin/Heidelberg, Germany, 1998; Volume 7. [Google Scholar]
  59. Kocks, U.F.; Tome, C.N.; Wenk, H.R. Texture and Anisotropy: Preferred Orientations in Polycrystals and Their Effect on Materials Properties; Cambridge University Press: Cambridge, CA, USA, 1998. [Google Scholar]
  60. Schmitt, J.H.; Iung, T. New developments of advanced high-strength steels for automotive applications. Comptes Rendus Phys. 2018, 19, 641–656. [Google Scholar] [CrossRef]
Figure 1. Kinematic decomposition of a material point into three intermediate configurations: Ω ¯ I , Ω ˘ I I , and Ω ˜ I I I .
Figure 1. Kinematic decomposition of a material point into three intermediate configurations: Ω ¯ I , Ω ˘ I I , and Ω ˜ I I I .
Crystals 12 00930 g001
Figure 2. Numerical integration algorithm for crystal plasticity model.
Figure 2. Numerical integration algorithm for crystal plasticity model.
Crystals 12 00930 g002
Figure 3. Finite element models subjected to (a) tension and (b) compression.
Figure 3. Finite element models subjected to (a) tension and (b) compression.
Crystals 12 00930 g003
Figure 4. Loaded directions (LDs) in single crystal with respect to specimen orientation (a) [100]-LD, (b) [110]-LD, and (c) [111]-LD.
Figure 4. Loaded directions (LDs) in single crystal with respect to specimen orientation (a) [100]-LD, (b) [110]-LD, and (c) [111]-LD.
Crystals 12 00930 g004
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.
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.
Crystals 12 00930 g005
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 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.
Crystals 12 00930 g006
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.
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.
Crystals 12 00930 g007
Figure 8. Deformation behavior of single-crystal and polycrystal TWIP steels subjected to tension.
Figure 8. Deformation behavior of single-crystal and polycrystal TWIP steels subjected to tension.
Crystals 12 00930 g008
Figure 9. Deformation behavior of single-crystal and polycrystal TWIP steels subjected to compression.
Figure 9. Deformation behavior of single-crystal and polycrystal TWIP steels subjected to compression.
Crystals 12 00930 g009
Figure 10. Magnitude of commulative slip shear strain of a single crystal of TWIP steels subjected to tension and compression.
Figure 10. Magnitude of commulative slip shear strain of a single crystal of TWIP steels subjected to tension and compression.
Crystals 12 00930 g010
Figure 11. Magnitude of commulative twin 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.
Crystals 12 00930 g011
Figure 12. Twin volume fraction of a single crystal of TWIP steels subjected to tension and compression.
Figure 12. Twin volume fraction of a single crystal of TWIP steels subjected to tension and compression.
Crystals 12 00930 g012
Table 1. Time sub-stepping algorithm for Newton–Raphson Iterative Scheme.
Table 1. Time sub-stepping algorithm for Newton–Raphson Iterative Scheme.
1. Calculate the values of K for each component of vector A
                                     K = Δ A max Δ A x
2. IF   K > 1.25
  THEN: Solution is rejected and define new time increment as:
                                     Δ t n + 1 = 0.75 ( Δ t ) n
    GOTO N-R iterative algorithm
  ELSE GOTO step 3
3. IF 0.8 K 1.25
  THEN: Solution is accepted and define new time increment as:
                                     Δ t n + 1 = ( Δ t ) n / K
    GOTO N-R iterative algorithm
  ELSE GOTO step 4
4. IF   K < 0.8
     THEN: Solution is accepted and define new time increment as:
                                     Δ t n + 1 = 1.25 ( Δ t ) n
        GOTO N-R iterative algorithm
     ELSE END
Table 2. General modeling parameters of TWIP steels.
Table 2. General modeling parameters of TWIP steels.
TypeParameter
Material
Moduli of rigidity G s = 111.0 , G t = 98.4 (GPa)
Bulk modulus K a = 206.5 (GPa)
Flow rule
Initial shear strain rate γ ˙ 0 = 0.001 s 1
Hardening rule
Initial hardening rate h 0 α = 200 (MPa)
Initial saturation value of slip resistance s r , S 0 α = 120 (MPa)
Saturation slip resistance exponent a = 0.005
Shear strain rate at saturation slip resistance γ ˙ S 0 = 5 × 10 10 s 1
Latent hardening coefficient q α κ = 1.4
Slip hardening parameter s r , p = 350 (MPa)
Boltzman’s constant, Equation (80) [53]k = 1.38 × 10 23 (J/K)
Absolute temperature, Equation (80) [53] θ = 298 (K)
Product of Burger’s vector and material parameter, Equation (80) [53] b 3 B = 0.005
Material parameter, Equation (80) [53]x = 0.5
Initial hardening rate of non-coplanar twin systems h n c i = 800 (MPa)
Initial hardening rate of coplanar twin systems h c p i = 8000 (MPa)
Defect energy
Crystal defect energy parameters φ = 10 , ω = 5
Initial crystal defect energy ζ 0 = 4.5 × 10 4 s 1
Thermal energy
Thermal analogous of resolved shear stresses Φ α = 12 , Φ i = 12 (MPa)
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.
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.
TWIP 22% Mn-0.6% C
TypeParameter
Material
Elasticity tensor components E 11 a = 286.80 , E 12 a = 166.40 , E 44 a = 145.10 (GPa)
Flow rule
Rate sensitivity parameter m = 0.03
Hardening rule
Initial slip resistance s r , 0 α = 70 (MPa)
Initial twin resistance s r , 0 i = 80 (MPa)
TWIP 30% Mn-0.5% C
TypeParameter
Material
Elasticity tensor components E 11 a = 286.80 , E 12 a = 166.40 , E 44 a = 145.10 (GPa)
Flow rule
Rate sensitivity parameter m = 0.02
Hardening rule
Initial slip resistance s r , 0 α = 80 (MPa)
Initial twin resistance s r , 0 i = 80 (MPa)
TWIP 29% Mn-0.8% C
TypeParameter
Material
Elasticity tensor components E 11 a = 286.80 , E 12 a = 166.40 , E 44 a = 145.10 (GPa)
Flow rule
Rate sensitivity parameter m = 0.02
Hardening rule
Initial slip resistance s r , 0 α = 90 (MPa)
Initial twin resistance s r , 0 i = 120 (MPa)
Table 4. Comparison of equivalent stress at 0.4 equivalent strain of single-crystal and polycrystal TWIP steels under tension and compression.
Table 4. Comparison of equivalent stress at 0.4 equivalent strain of single-crystal and polycrystal TWIP steels under tension and compression.
Equivalent Stress (MPa)
TensionCompression
Poly[100][110][111]Poly[100][110][111]
TWIP 11094549626839688542542540
TWIP 21138510512820678560555552
TWIP 31230575577873712570568565
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.
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.
Twin Volume Fraction
TensionCompression
[100][110][111][100][110][111]
TWIP 10.2190.1600.2460.2110.1920.191
SlipTWIP 20.1320.1070.2130.1420.1120.162
TWIP 30.2350.2110.2450.2380.2080.196
TWIP 10.0150.0160.0080.0170.0160.007
TwinTWIP 20.0910.0870.0330.0910.0860.024
TWIP 32.4 × 10 5 1.4 × 10 5 9.8 × 10 5 2.4 × 10 5 1.3 × 10 5 5.1 × 10 6
Table 6. Highest values of twin volume fraction for TWIP 1, 2, and 3 steels under tension and compression.
Table 6. Highest values of twin volume fraction for TWIP 1, 2, and 3 steels under tension and compression.
Twin Volume Fraction
TensionCompression
[100][110][111][100][110][111]
TWIP 10.3330.1660.2490.3330.1660.166
TWIP 20.2490.1660.1660.3330.1660.166
TWIP 30.2490.1660.1660.2000.1100.080
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

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

AMA Style

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

Chicago/Turabian Style

Khan, Rashid, Tasneem Pervez, Adel Alfozan, Sayyad Zahid Qamar, and Sumiya Mohsin. 2022. "Numerical Modeling and Simulations of Twinning-Induced Plasticity Using Crystal Plasticity Finite Element Method" Crystals 12, no. 7: 930. https://doi.org/10.3390/cryst12070930

APA Style

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

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop