TRIP Steels: A Multiscale Computational Simulation and Experimental Study of Heat Treatment and Mechanical Behavior

A multiscale investigation of the microstructure and the mechanical behavior of TRIP steels is presented. A multi-phase field model is employed to predict the microstructure of a low-alloy TRIP700 steel during a two-stage heat treatment. The resulting stability of retained austenite is examined through the Msσ temperature. The phase field results are experimentally validated and implemented into a model for the kinetics of retained austenite during strain-induced transformation. The kinetics model is calibrated by using experimental data for the evolution of the martensite volume fraction in uniaxial tension. The transformation kinetics model is used together with homogenization methods for non-linear composites to develop a constitutive model for the mechanical behavior of the TRIP steel. A methodology for the numerical integration of the constitutive equations is developed and the model is implemented in a general-purpose finite element program (ABAQUS). Necking of a bar in uniaxial tension is simulated and “forming limit diagrams” (FLDs) for sheets made of TRIP steels are calculated. The models developed provide an integrated simulation toolkit for the computer-assisted design of TRIP steels and can be used to translate mechanical property requirements into optimised microstructural characteristics and to identify the appropriate processing routes.


Introduction
TRansformation Induced Plasticity (TRIP) steels is a new generation of low alloy steels that exhibit an enhanced combination of strength and ductility. The microstructure of TRIP steels consists typically of ferrite, bainite, and retained austenite. Their remarkable mechanical properties are attributed to the TRIP phenomenon, i.e., the retained austenite transforms to martensite with plastic deformation. In the present work, a multiscale investigation of the microstructure and the mechanical properties of TRIP steels is presented. This investigation can be conducted for other steel grades applied in the automotive industry with the appropriate calibration.
The stabilization of retained austenite against martensitic transformation is a key criterion for the design of TRIP steels [1,2]. A common approach to stabilize the retained austenite, is to perform a two-stage heat treatment, consisting of intercritical annealing prior to isothermal bainitic transformation and quenching to room temperature [3,4]. During the isothermal bainitic transformation the remaining austenite is enriched with sufficient carbon and stabilized as retained austenite. Several studies were carried out in order to investigate and predict the microstructural evolution during intercritical annealing [5][6][7], bainitic isothermal transformation [8][9][10] and martensitic transformation [11].
The determination of retained austenite stability, expressed by the M σ s temperature, and the description of retained austenite strain-induced transformation kinetics are essential in predicting the mechanical performance of TRIP steels. The chemical composition and grain size of retained austenite, as well as the stress triaxiality and the strength of the austenite matrix were considered to be the primary influencing factors on M σ s [12][13][14][15]. Several models, which are reviewed by Samek et al. [16], were developed for the transformation kinetics of retained austenite. In addition, Haidemenopoulos et al. [17] took into account the effect of austenite particle size on transformation kinetics. Regarding the description of the mechanical behavior of TRIP steels, several constitutive models were developed [18][19][20]. However, these models disregard the effect of austenite particle size on transformation kinetics. Zhu et al. [21] and Shengyen [22] employed computational-based approaches for the design of TRIP steels, yet a multiscale integrated model is required for the investigation of the microstructure and the mechanical properties of the investigated steels. A brief description of each section of the present work is given below.
In Section 2, an integrated simulation model is employed to predict the microstructural characteristics of a cold rolled low-alloy TRIP700 steel subjected to a two-stage heat treatment. The heat treatment consists of intercritical annealing and an isothermal bainitic transformation at a lower temperature. A two-dimensional (2D) multi-phase field model (MPF) is developed for the simulation of the intercritical annealing process after deformation in a ferrite-pearlite initial structure. The concentration profiles, volume fractions of the phases, and microstructure obtained at the end of the intercritical annealing are implemented into the MPF model for the simulation of the isothermal bainitic treatment. The MPF model is used to determine the volume fractions and the average grain size of the phases as well as their chemical content in the resulting microstructure at the end of the heat treatment.
The experimental validation of the predictions for the heat treatment process is discussed in Section 3. The average grain size diameter of ferrite and retained austenite are determined using optical microscopy. The Saturation Magnetization (SM) technique is employed to evaluate the volume fraction of retained austenite. By applying the SS-TV-TT technique to determine the M σ s temperature, the austenite stability is measured.
The austenite stability, which influences the kinetics of transformation plasticity, is studied in Section 4 by using the model of Haidemenopoulos and Vasilakos [15] for the stability of dispersed austenite in low alloy steels. This model predicts the effects of several parameters on austenite stability. An attempt is made to distinguish the effects of individual parameters on austenite stability. The predictions of the model are consistent with the experimental results for the M σ s temperature. Also, in Section 4, the model of Haidemenopoulos et al. [17] is employed to predict the evolution of the volume fraction of martensite during the strain-induced transformation of retained austenite. The microstructural characteristics calculated from the simulation of the heat treatment are used as input to this model, which takes into consideration the chemical composition, the austenite particle size, temperature, and stress state. The calibration of the model is based on experimental data for the evolution of the martensite volume fraction in uniaxial tension tests.
In Section 5, the aforementioned transformation kinetics model is used together with non-linear homogenization methods of Ponte Castañeda and co-workers [23][24][25] to develop constitutive equations for the mechanical behavior of TRIP steels. The methodology of Papadioti et al. [26], in which TRIP steels are treated as multi-phase (composite) materials, is used to estimate the effective yield function and the average stress and strain fields in the constituent phases.
A methodology for the numerical integration of the resulting constitutive equations in the context of a displacement driven finite element formulation is presented in Appendices A and B. The model is implemented in the ABAQUS general-purpose finite element code, one-element finite element calculations for the problem of uniaxial tension are performed, and the simulation results are compared with experimental data in Section 6. In Section 7, the constitutive model is employed for the study of necking under uniaxial tension and forming limit diagrams (FLDs) are also calculated.
The scope of the present work, is to provide an integrated simulation toolkit, which will be used for the design of TRIP steels. Starting with the simulation of the microstructure evolution with respect to the heat treatment parameters, the chemical composition and the grain size of retained austenite are calculated. These results are used as an input for the description of the stability of retained austenite and the martensite volume fraction formed as a function of plastic strain in TRIP steels. After the experimental validation of the heat treatment predictions and the M σ s temperature, a constitutive model is developed for the mechanical behavior of TRIP steels with respect to the processing parameters and microstructural features. The integrated model presented in this work links composition, processing, microstructure and mechanical behavior and makes it a potential design tool for the development of optimized TRIP steels.
Throughout the text, standard notation is used, while boldface symbols correspond to tensors, the orders of which are indicated by the context. A fixed Cartesian coordinate system with base vectors e i (i = 1, 2, 3) is used for the expression of all tensor components, and Greek indices take the values 1 and 2. The summation convention is applied for repeated latin and Greek indices ( A superscript T denotes the transpose, and a superposed dot the material time derivative. Let (A, B) be second-order tensors and (C, D) fourth-order tensors; the following products are used in the The inverse C −1 of a fourth-order tensor C that has the "minor" symmetries C ijkl = C jikl = C ijlk is defined so that C : C −1 = C −1 : C = I, where I is the symmetric fourth-order identity tensor with Cartesian components I ijkl = 1 2 (δ ik δ jl + δ il δ jk ), δ ij being the Kronecker delta.

Materials and Modeling of the Heat Treatment
An integrated process chain model is employed for the prediction of the microstructural features of a cold rolled (CR) TRIP700 steel during a two-stage heat treatment, which consists of intercritical annealing prior to isothermal bainitic transformation at a lower temperature. TRIP 700 is a well-established TRIP steel, this is why it was used as a model steel in the investigation. This investigation can be carried out for other steel grades used in the automotive industry with the appropriate calibration. The chemical composition of the CR-TRIP700, provided by voestalpine Stahl GmbH Linz, is described in Table 1. After cold-rolling, the aim is to form a 50%-50% ferrite-austenite microstructure during intercritical annealing in the first stage and to stabilize the retained austenite through an isothermal bainitic transformation in the second stage.

Selection of Heat Treatment Parameters
An integrated process chain model is developed using multi-phase field modeling. The thermodynamic calculations, required to set the bounds of the process, are performed using the Thermo-Calc software [27] and the TCEF6 database for ferrous alloys. The limits for the intercritical annealing are A 1 = 693.4 • C and A 3 = 917.8 • C, and the cementite solvus temperature is A cem = 716.8 • C. During intercritical annealing, the target volume fraction of austenite is 50%. Thermodynamic calculations show that the intercritical annealing should take place above 808.3 • C. An intercritical annealing temperature at 890 • C is selected with a heating rate of 15 • C/s and an intercritical annealing time of 60 s, followed by an isothermal bainitic treatment consisting of cooling with a rate of 60 • C/s and an isothermal holding at 400 • C for 120 s.

Methodology
The multi-phase field model (MPF), developed by Steinbach et al. [28][29][30][31] and implemented in MICRESS [32], is used for the description of microstructure evolution during intercritical annealing and isothermal bainitic treatment. The adopted approach takes into account the effect of grain boundary curvature and allows the determination of the spatial distribution of phases and their grain size [33].

Stage I-Intercritical Annealing
For the simulation of the intercritical annealing process after deformation of a ferrite-pearlite initial structure, a 2D MPF model is employed. Since the heating rate is 15 • C/s, it is considered that ferrite recrystallization is completed and austenite nucleates mainly in pearlite [5,7,[34][35][36]. A very small amount of austenite nuclei is allowed to form at ferrite grain boundaries, when it is thermodynamically possible. Spheroidization of pearlite is neglected. A sufficiently large domain size is used in the calculations, so that there are no significant statistical differences in the results [37]. However the description of a fine lamellar pearlite structure is not possible in this study, since a smaller domain size is required [38,39]. For this reason, pearlite is modelled as an effective phase with uniform eutectoid composition of 0.86% wt. carbon and an equilibrium phase fraction of 22.9%. The simulations are carried out in two steps. In the first step, reheating from the A 1 temperature to the selected intercritical annealing temperature takes place for a time period sufficient for complete dissolution of pearlite. A rapid pearlite to austenite transformation occurs and austenite is assumed to inherit the carbon content of pearlite; the substitutional elements are assumed to be immobile. The growth kinetics of the phases and the redistribution of the alloying elements are described by non-partitioning local equilibrium conditions (NPLE). Austenite is allowed to nucleate at ferrite grain boundaries. The second step starts after the complete dissolution of pearlite. During the second step, ferrite transforms to austenite and the austenite formed in the first step continues to grow. The ferrite to austenite phase transformation takes place under ortho-equilibrium conditions.
The domain size is 85.5 µm × 40.5 µm. A grid of square cells with sides δx = 0.15 µm is used in the calculations. The interface thickness is set to 0.6 µm. The initial structure is constructed randomly and periodic boundary conditions are applied. The Thermo-Calc TCEF6 and MOB2 databases are used to derive diffusion data. The same diffusion coefficients are used for pearlite and ferrite. Thermodynamic data for the ferrite-pearlite and pearlite-austenite interactions are obtained using a linearized phase diagram; the TCEF6 database is used to get the corresponding data for the ferrite-austenite interaction. The local equilibrium temperature is set according to the linearized phase diagram at 744 • C, where the equilibrium carbon content is 0.86% wt. and 0.0063% wt. in pearlite and ferrite respectively. Manganese and aluminum are considered immobile during the first step of the simulations; their nominal compositions are selected as the equilibrium concentrations in the phases at 744 • C.
The interfacial energies σ ij for the phase interactions take the constant values listed in Table 2 [35,37]. The interface mobility µ ij between phases i and j is considered to be temperature-dependent: where µ o ij is the pre-exponential mobility factor, R the gas constant, and Q ij the activation energy for grain boundary migration. The pre-exponential factor is written in the form where d ij is the interatomic distance between phases i and j, v D the Debye frequency of the parent phase, and k the Boltzmann constant. The mobility parameters for the phase interactions Q ij and d ij are taken from the literature [40,41] and are listed in Table 2. In the first step of simulations, austenite nucleates within pearlite and grows until complete dissolution of pearlite. Nucleation of austenite within pearlite is assumed to be instantaneous (site-saturated). The growth of austenite is controlled by carbon diffusion in austenite and the partitioning of substitutional atoms in ferrite [42]. It is assumed that the critical surface area of a nucleus is 1 µm 2 . The number of austenite potential nucleation sites in pearlite is N n = 793 and equals the product of the pearlite phase fraction and the domain size divided by the nucleus surface area. Austenite nucleation is allowed to occur at ferrite grain boundaries both in the first and the second step, if it is thermodynamically possible. The expression that defines the heterogeneous nucleation rate of austenite in ferriteṄ S is based on classical nucleation theory under steady state conditions [35] and is of the formṄ where N n (t) is the number of potential nucleation sites as a function of time, a superposed dot denotes derivative with respect to time, ψ is a constant related to the geometry of the nucleus (shape) and the interfacial properties for the nucleus and the matrix, ∆g v is the driving force for nucleation of austenite in ferrite per unit volume and was calculated using Thermo-Calc, h and k are the Planck and Boltzmann constants accordingly, and Q D is the activation energy for nucleation, which is approximately equal to the activation energy of iron self-diffusion in ferrite and takes the value of 3.9 × 10 −19 J. The value of ψ a = 4.76 × 10 −10 J 3 /m 6 is used for nucleation on ferrite grain boundaries (Savran [40]). The value of N n is kept constant in the calculations and its value equals the product of ferrite grain boundary perimeter times the interface thickness divided by the nucleus surface area. In the 2D simulation domain, the number of austenite potential nucleation sites is N n = 2436.6 µm × 0.6 µm/1 µm 2 = 1462. A random orientation of austenite nuclei is considered with isotropic austenite and ferrite grain growth.

Stage II-Isothermal Bainitic Treatment
Bainitic ferrite is assumed to form by diffusional nucleation and displacive transformation [43]. The nucleation and growth of the entire sheaf is considered instead of the formation of sub-units. Diffusion of the alloying elements is allowed. Since excess Si and Al content in steels can suppress carbide precipitation, cementite precipitation is not considered neither in bainite nor in austenite. Plastic strains developed during the sub-unit growth are not accounted for and the plate shape morphology of bainitic ferrite is acquired artificially by using the anisotropic interfacial properties. Accordingly, the model cannot account the incomplete reaction phenomenon of bainite due to dislocation debris. Solute drag effect of substitutional atoms are ignored in the interface kinetics and Zener pinning forces due to the precipitation of carbonitrides are not included in this study.
The concentration profiles, volume fractions of the phases, and microstructure obtained at the end of the intercritical annealing are implemented into the MPF model for the simulation of the isothermal bainitic treatment. The same domain size and the interface thickness are used in the calculations. The thermodynamic data for the interaction of ferrite and austenite are derived by coupling the multi-phase field calculations with the Thermo-Calc TCEF6 and MOB2 databases. To simplify the simulations, no phase interaction between ferrite and bainite is considered. For the description of austenite-bainite phase interaction, a linearized phase diagram is calculated at the local equilibrium temperature of 744 • C. The same diffusion coefficients, interfacial energies, and interface mobilities are used both for ferrite and bainitic ferrite phases. A value of σ γ−a = 0.2 J/m 2 [44] is used for the interfacial energy between austenite and bainite; the interface mobility at 400 • C is set to µ γ−a = 9 × 10 −9 cm 4 /(J s).
Bainitic ferrite is modelled as an anisotropic phase with tetragonal symmetry. In TRIP steels there is a Kurdjumov-Sachs (KS) orientation relationship that describes the relative orientation of retained austenite and bainitic ferrite [45]. It is assumed that the misorientation angle between bainite grains and parent austenite grains follows a KS relationship with a tolerance angle of 10 • . The transformation model of Bhadeshia [46] is used for the calculation of the bainite start temperature. The required driving force per mole for nucleation in austenite ∆g v and the diffusionless driving force for bainite formation ∆G γ→a , which is regarded as the Gibbs free energy difference of ferrite and austenite, are calculated using Thermo-Calc. Bainite start temperature is calculated as B S = 616.8 • C. The nucleation rate of bainitic ferrite is calculated according to equation (3), where the activation energy for self-diffusion of carbon in austenite Q D = 142.1 kJ/mol is used [7]. For bainite nucleation at austenite grain boundaries the value ψ a = 6.33 × 10 −15 J 3 /mol 2 is used [47]. The number of bainite potential nucleation sites is represented by the total surface area of austenite grain boundaries, which is the product of the austenite grain boundary perimeter times the interface thickness divided by the nucleus surface area, and takes the value N n = 2573.3 µm × 0.6 µm/1 µm 2 = 1544. Figure 1a shows the initial structure, which consists of 99 recrystallized ferrite grains and 36 elongated pearlite grains. The initial average grain size is 5.65 µm for pearlite and 6.24 µm for ferrite, where equivalent circular mean diameters are considered. The microstructure after the complete dissolution of pearlite is depicted in Figure 1b. Symbols P, F, and A denote pearlite, ferrite and austenite grains respectively. The simulation time required for the complete dissolution of pearlite was 9 s at 828 • C. The formed austenite nuclei in pearlite are 164, whereas no austenite nuclei were observed at ferrite grain boundaries during the first step of simulations.  In the second step, the formed austenite grows at the expense of ferrite, while some austenite nuclei grow at ferrite grain boundaries until the end of intercritical annealing at 890 • C. After 60 s of intercritical annealing, a 51% volume fraction of austenite was achieved; the corresponding microstructure is shown in Figure 1c. The microstructure at the end of the isothermal bainitic treatment at 400 • C for 120 s is presented in Figure 1d. Symbols RA and B correspond to the retained austenite and the bainitic ferrite.

Simulation Results of the Heat Treatment
The evolution of the volume fraction of the phases and the average grain diameter of pearlite, ferrite, and austenite during intercritical annealing are depicted in Figure 2a. Continuous lines correspond to volume fractions and dashed lines indicate the average grain size. The vertical dotted lines indicate the end of the first step and the time (∼60 s) at which 51% volume fraction of austenite was achieved; the horizontal dotted line marks the targeted 50%-50% ferrite-austenite microstructure. It should be noted that austenite grows faster in pearlite during the first step; the rate of austenite growth in ferrite is slower in the second step. The average grain diameters of ferrite and austenite at the intercritical annealing time of t = 60 s are 5.31 µm and 3.72 µm respectively.
The average concentration of carbon, aluminum, and manganese in ferrite and austenite are plotted versus time in Figure 2b. At the beginning of the second step, there is a carbon enrichment in ferrite. However, after 13 s, when the intercritical temperature of 890 • C has been reached, the average carbon content in ferrite starts decreasing due to the decrease of the maximum soluble amount of carbon in ferrite at this temperature. This temporal carbon enrichment in ferrite is favourable for the austenite formation at the ferrite grain boundaries. Carbon depletion in austenite is still observed due to the increase of austenite phase fraction. Regarding the substitutional atoms, there is manganese depletion and aluminum enrichment in ferrite.
The evolution of the volume fractions and the average grain diameter of ferrite, retained austenite, and bainitic ferrite during the isothermal bainitic treatment are shown in Figure 2c. Continuous lines correspond to volume fractions and dashed lines indicate the average grain diameter of the phases. Bainitic ferrite starts to form after 6 s of cooling and isothermal holding. When bainite is formed, the austenite shrinks and ferrite grains grow slowly. After 120 s, the resulting microstructure in the new TRIP steel consists of 50.7% ferrite, 16.6% retained austenite, and 32.7% bainitic ferrite, and their average grain diameters are 5.39 µm, 1.62 µm, and 0.96 µm accordingly. For the rest of the paper, the term "TRIP steel A" is used to refer to this microstructure.
The average concentration of carbon, aluminum, and manganese in ferrite, retained austenite, and bainitic ferrite are depicted in Figure 2d. During bainite transformation, carbon and manganese are ejected from bainite in austenite, while aluminum diffuses into bainite and ferrite until equilibrium is reached. The carbon content in retained austenite after 120 s is 1.17% wt.

Heat Treatment
A TRIP700 steel which consists of ferrite, retained austenite, bainitic ferrite with traces of martensite, was employed for the experimental validation of the heat treatment process presented in Section 2. An austenitization process was performed above the A 3 temperature, at 950 • C, to erase its current microstructure. After the austenitization, the material is air cooled to room temperature, and a mixed microstructure consisting of ferrite, pearlite, and a small amount of bainitic ferrite was formed. It should be noted that for the simulation of the heat treatment process, an initial deformed ferrite-pearlite structure was assumed prior to the intercritical annealing process. This assumption leads to different initial microstructures for the intercritical annealing between the experimental study and the phase-field simulations. As a result, minor deviations between experimental data and simulation results, regarding the morphology and the phase volume fractions, are observed.
The heat treatment process is presented schematically in Figure 3.  The heating process consists of three steps that include austenitization, intercritical annealing, and austempering for bainitic transformation. Six specimens with dimensions 120 mm ×100 mm ×1.5 mm were preheated in a front-load furnace, at 400 • C for 60 s and then directly inserted in a second furnace for the austenitization process at 950 • C for 240 s and air cooled to room temperature. Oxidation products from the surface of the samples were carefully removed by grinding with emery paper.
Intercritical annealing and austempering treatments were performed in furnaces with molten salt baths using GS750 and AS140 salts respectively. Following the austenitization process the samples were firstly preheated for 60 s at 400 • C and annealed at 890 • C for 60 s. Subsequently the samples cooled down to 400 • C for 120 s and air cooled to room temperature.

Evaluation of the M σ s Temperature
The stability of retained austenite (RA) was evaluated by measurements of the M σ s temperature, which defines the critical boundary between the stress-assisted and strain-induced martensitic transformation in TRIP steels. M σ s was determined experimentally with the SS-TV-TT technique [48,49] using the experimental setup shown in Figure 4. Loading-unloading was carried out in an INSTRON 8801 servo-hydraulic machine at a constant crosshead velocity of 0.5 mm/min, while a Real Time Strain Sensor (RTSS) video-extensometer was used to measure the longitudinal tensile strains during testing ( Figure 4). Testing at low temperatures was performed using a plastic container with cooling bath attached at the gauge length area of the specimen. The desired temperatures, ranging from 20 • C to −30 • C, were achieved using appropriate concentrations of water, salt, ice, acetone, dry ice, and ethylene-glycol in the cooling bath.  The M σ s temperature was estimated from the stress-strain behavior during successive loading of the material up to yielding with decreasing temperature. The stress-strain curves obtained are depicted in Figure 5. They reveal an onset of stress relaxation at −19 • C, a result of the volume expansion accompanying the martensitic transformation phenomenon, i.e., the M σ s temperature is found to be −19 • C, meaning that below that temperature martensitic transformation can only be stress-assisted; above −19 • C, strain-induced martensitic transformation takes place when austenite deforms plastically at temperatures below an upper limit (M D ).

Microstructural Characterization
The microstructure after heat treatment was assessed using optical microscopy. Standard metallographic techniques of grinding and polishing were used in mounted sample. Consecutively, a stepped color tint-etching procedure was applied to reveal the microstructure. The sample was first etched with a 3% Nital solution for 3-4 s, and then with a 10% Na 2 S 2 O 5 solution for 60 s. After etching in both solutions, the specimen was thoroughly washed with water, immersed in ethanol, and dried in hot air.
The average grain diameter of retained austenite (RA) and ferrite were estimated using an image analysis software [50]. Since there is a large dispersion of phases, several micrographs were taken at different depths. This method is suitable for structures with high grain shape irregularity [50]. Approximately 1000 retained austenite and 100 ferrite grains were used for the quantification of the above microstructural characteristics.
The RA-white colour, ferrite-blue/straw brown, and bainite-dark, microstructural characteristics are displayed in Figure 6, revealing a ferritic-bainitic matrix with a fine dispersion of RA. The average size of RA particle was 0.59 µm, while the average grain size of ferrite was 6.19 µm.

Retained Austenite Measurements
The %RA volume fraction in bulk material was measured with the Saturation Magnetization (SM) technique from small extracted samples of dimensions 14 mm ×3.5 mm. The SM measurements were performed in the Steel Division Department of voestalpine in Austria. Detailed description of the performed technique is presented in [51,52].
The advantages of this method are the ability to measure the entire volume of the specimen, the good repeatability [52], and the high precision compared to other methods [53] employed for the determination of the RA volume fraction. However, small uncertainties in the measurements may occur due to the existence of alloying elements, which can affect the saturation magnetization. For this reason, certain scatter in RA volume fraction measurements should always be taken into consideration.
An initial %RA volume fraction of 15.9% was measured, which is in a good agreement with the value of 16.6% calculated from the simulation of the heat treatment described in Section 2.3.

Stability and Transformation Kinetics of Retained Austenite
A thorough investigation regarding the stability of dispersed austenite in TRIP steel A is conducted based on the work of Haidemenopoulos and Vasilakos [15]. The M s temperature is used for the characterization of austenite stability against transformation on cooling, whereas the M σ s temperature for the stability of dispersed austenite against mechanically induced transformation. During the isothermal bainitic transformation, the retained austenite is chemically enriched in carbon and manganese [54]. Primarily, the carbon content is a strong stabilizer of the retained austenite, leading to suppression of the M s temperature. However, a mechanically induced martensitic transformation can take place in dispersed austenite by strain-induced or stress-assisted nucleation [55,56]. A good combination of all the interrelated factors, which have a strong stabilizing influence on the kinetics of the mechanically induced transformation, is required so as to achieve the maximum TRIP effect.
The aforementioned interrelation makes it difficult to separate the individual effect of each factor on austenite stability. This model is used to predict the M σ s temperature and to identify and quantify the primary stabilizing factors in the TRIP steel. The model predictions (M σ s temperature) are also compared with the experimental data obtained with the SS-TV-TT technique as described in Section 3.2. Furthermore, the model developed by Haidemenopoulos et al. [17] is employed to predict the evolution of the strain-induced martensitic transformation kinetics in the TRIP steel. In the present work, the M σ s temperature is chosen as the single parameter characterizing the stability of dispersed austenite particles in the TRIP microstructure. A variation of the methodology developed by Haidemenopoulos and Vasilakos [15] is used for the calculation of M σ s . According to the Olson-Cohen theory of heterogeneous martensitic nucleation [57][58][59], the dissociation of an existing defect serves as a potential nucleation site for the formation of a martensitic embryo. The energy per unit area γ f (n) (J/m 2 ) of an embryo with a thickness of n crystal planes is where ρ is the atomic density in the close-packed fault plane (mol/m 2 ), ∆G ch (J/mol) the thermodynamic driving force for martensitic transformation, E str (J/mol) the elastic strain energy related to the distortions in the interface plane of the fault structure, W f (J/mol) the frictional work of interfacial motion accompanying the dissociation process, and γ s (J/m 2 ) the fault/matrix interfacial energy. The thickness n refers to the number of crystal planes forming the defect and its value can be interpreted as a potency measure of the nucleation site. Barrierless dissociation of an existing defect occurs when its thickness n is such that γ f ≤ 0. The condition γ f (n * ) = 0 defines the critical thickness n * of a defect: If the thickness n of the defect is greater than or equal to n * , martensitic transformation takes place at that defect. The critical nucleation thickness of a defect n * depends on temperature through ∆G ch and on the chemical enrichment of the austenite, which affects ∆G ch and W f .
To take into account the stress-dependence of n * , the mechanical driving force ∆G σ (J/mol) is introduced in the denominator of Equation (5). Then, the total driving force is ∆G = ∆G ch + ∆G σ , so that Equation (5) yields where ∆G σ (σ) is a function of the stress tensor σ, T is temperature, and X denotes dependence on chemical composition. The last equation is based on the assumption that the orientation of the operative nucleation sites is optimum for maximum interaction with the applied stress (Patel and Cohen [60]). At the other extreme of a random distribution of nucleation sites, the term ∆G σ should be replaced by ∆G σ /3 in (6) (Olson et al. [61]). In isotropic materials, ∆G σ (σ) is in general a function of the von Mises equivalent stress σ e , the "stress triaxiality" Σ, and the "Lode angle" θ: where Σ = p σ e , p = σ kk 3 is the hydrostatic stress, θ = 1 3 arcsin − 27 2 dets σ 3 e , s = σ − p δ is the stress deviator, and δ is the second-order identity tensor.
Based on the small-particle experiments of Cech and Turnbull [62] in Fe-30%Ni alloys, Cohen and Olson [63] derived the following expression for the cumulative defect potency distribution per unit austenite volume N v m −3 : where (6) was used, N 0 v m −3 is the total number of all potential nucleation sites per unit austenite volume and a is a shape factor constant. N v defines the number of nucleation sites per unit austenite volume with sufficient potency to nucleate martensite (operational sites).
A uniform distribution of austenite particles with average volume v p is considered. It is assumed that a single nucleation event transforms the particle into martensite. Therefore, the fraction f of particles to transform through sites of cumulative number density N v corresponds to the probability of finding at least one nucleation site in the particle (Cohen and Olson [63]): where (8) was used and c (a) is the available volume fraction of retained austenite for transformation. For given temperature T, composition X, and average particle size v p , Equation (9) defines f in terms of stress σ. Equation (9) can be described by the equivalent form For given triaxiality Σ and Lode angle θ, last equation defines implicitly σ e in terms of temperature T, chemical composition X, average particle size v p , and martensite volume fraction f , i.e., it yields a relationship of the form The temperature M σ s is defined as the maximum temperature at which transformation is induced by a stress below the yield stress of the parent phase (Bolling and Richman [64], Olson and Cohen [65]). For temperatures below M σ s , an elastic stress-assisted transformation takes place at the same pre-existing nucleation sites, which are responsible for the spontaneous transformation on cooling. Let σ y (T) be the temperature-dependent yield stress of austenite. The experimental data of Olson and Azrin [1] show clearly that for T ≤ M σ s , the 1% martensite curves are almost coincident with the 0.2% yield stress curves ( Figures 5 and 6 in [1]), i.e., σ = σ y and f /c (a) ≡ f 1 = 0.01 at T = M σ s . Therefore, by balancing the transformation stress in (11) with the σ y (T), the M σ s temperature can be obtained, i.e., M σ s is defined from the condition Equation (10) leads to the following solution for M σ s : Last equation defines implicitly the M σ s temperature in terms of chemical composition X, average particle size v p , and the stress state parameters (Σ, θ): A detailed calculation of M σ s for the TRIP steel under consideration is developed next. Haidemenopoulos et al. [66], based on the data of Olson and Cohen [65], suggest the following expression for ∆G σ , which depends on σ e and Σ, but is independent of θ: The individual terms of the chemical driving force ∆G ch and the frictional work of interfacial motion W f in Equation (6) were discussed in detail in the work of Haidemenopoulos and Vasilakos [15], who propose the following forms where X C , X Mn are the mole fractions of C and Mn in the retained austenite respectively, and T is the temperature in K. According to Olson and Cohen [58], the elastic strain energy is E str = 500 J/mol. The yield stress data of Samek et al. [16] for the austenite present in TRIP steels show that the austenite flow stress decreases with temperature at a rate of 1 MPa/ • C. Therefore, σ y (T) is described by where σ 0 = 530 MPa is the austenite flow stress at a temperature T 0 = 25 • C = 298.15 K and B = 1 MPa/K. Substitution of (15)- (18) in (13) leads to the following solution for M σ s : with where A is determined from (20) in J/mol, σ 0 is in MPa, B in MPa/K, T 0 in K, and M σ s is calculated from (19) in K.

Calculation of M σ s and the Effects of the Variants on Austenite Stability
Based on the small-particle experiments of Cech and Turnbull [62] in Fe-30%Ni alloys, Cohen and Olson [63] evaluated the shape factor a in (8) to be α = 0.84. The size of the particles involved in those studies was of the order of 10 µm. In the TRIP steel under consideration the austenite particles are dispersed and their size is smaller than 1 µm. Therefore, the shape factor takes the smaller value a = 0.12 (Haidemenopoulos et al. [17]). The values of the remaining parameters in Equation (20) are as follows [15,55]: γ s = 0.15 J/m 2 , ρ = 3 × 10 −5 mol/m 2 , and N 0 v = 2 × 10 17 m −3 . Under the assumption that the austenite particles are spheroidal, the austenite particle volume v p = 4 3 π R 3 was calculated using the mean austenite particle radius R = 0.8 µm, as determined in the phase field calculations of Section 2.3. The yield strength of retained austenite is equal to σ 0 = 530 MPa. The calculation of the M σ s temperature is carried out in two steps. First, the effects of various parameters on the retained austenite stability are thoroughly examined. Then, the M σ s temperature is calculated using (19) together with (20).
For the case of uniaxial tension (Σ = 1/3), M σ s is plotted as a function of the yield strength of austenite σ 0 in Figure 7. An increase in the yield strength leads to an increase in the M σ s temperature. Figure 8 shows the calculated M σ s temperature with respect to the mean austenite size in uniaxial tension. The effect of particle size refinement on the austenite stability is evident: M σ s decreases as the austenite dispersion becomes more refined. The predictions of the model are in agreement with the multitechnique investigation of Haidemenopoulos et al. [67]. The effect of stress state on M σ s is shown in Figure 9 for R = 0.8 µm. Four cases regarding the stress triaxiality are marked on the triaxiality axis: Σ equals −1/3 in uniaxial compression, 0 in pure shear, 1/3 for tension, and 1/ √ 3 = 0.5777 in plane strain tension. The destabilizing effect of stress triaxiality is evident: the higher the value of stress triaxiality Σ, the higher the M σ s . Between the two extreme cases A (uniaxial compression) and D (plane strain tension) a 30 • C rise in M σ s is observed. For the TRIP steel under consideration the M σ s temperature in uniaxial tension is predicted to be −20.3 • C. This is in reasonable agreement with the measured value of −19 • C with the SS-TV-TT technique, as reported in Section 3.2.

Transformation Kinetics of Retained Austenite
The model developed by Haidemenopoulos et al. [17], for the description of the strain-induced martensitic transformation kinetics of dispersed austenite, is calibrated using available experimental data for the evolution of martensite volume fraction f in a uniaxial tension test for the TRIP700 steel under consideration (Bellas [68], Haidemenopoulos et al. [67]).
Kuroda [69] suggested that, in strain-induced transformation, the potency distribution N v in Equation (8) should depend on both stress and strain, so that is still defined by an equation similar to (8): in which N σ0 v m −3 denotes the pre-existing nucleation sites and a σ is a constant shape factor. During plastic deformation, the austenite phase generates new more potent nucleation sites N ε0 v m −3 , which are calculated by an expression of the form (Haidemenopoulos et al. [17]) whereε (a) is the average equivalent plastic strain in austenite, N m −3 is the maximum number of sites that can be produced by plastic strain, and (k, m f ) are constants. As discussed in Haidemenopoulos et al. [17], the total number of operational sites due to the plastic strain in the where a ε is the shape factor in the strain-modified potency distribution. Then, (9) takes the form [17] f σ e , Σ, Equation (24) shows that the size of the austenite particles affects the amount of martensite produced during the transformation through the v p term.

Calibration of the Model to the Available Experimental Data
Equation (24) that describes the kinetics of strain-induced transformation was fitted non-linearly to available experimental data [67]. The maximum sites that can be formed by plastic deformation per unit volume N and the pre-existing nucleation sites N σ0 v are the parameters to be determined. The high chemical stability of the retained austenite ( Figure 2d) prevents stress assisted transformation, i.e., the high value of ∆G ch prevents the vanishing of γ f (n). In the TRIP steel under consideration, only the strain-induced transformation of retained austenite to martensite takes place. Therefore, the only parameter to be determined is N and the resulting value is 2.45 × 10 22 m −3 .
The predicted evolution of f curve in uniaxial tension together with experimental data is presented in Section 6.
The values of the mean particle radius of the retained austenite, the chemical composition, the elastic strain energy E str and the rest of the parameters and constants are the same as those used in the calculation of M σ s in Section 4.2. Since the experimental data of the volume fraction of martensite f were obtained at room temperature, the values of ∆G ch = −1828 J/mol and W f = 1710 J/mol were calculated at 25 • C using Equations (16) and (17) respectively.

Description of the Constitutive Model
In this section, a constitutive model for TRIP steels is developed. A four-phase TRIP steel that consists of a ferritic matrix with a fine dispersion of bainite and austenite is considered; due to plastic deformation, the retained austenite transforms gradually into martensite. For each phase, the following labels are used: (1) for ferrite, (2) for bainite, (3) or (a) for austenite, and (4) or (m) for martensite. The constitutive equations are developed for finite strain problems.
The strain softening which results from the strain related to the transformation process, is an important characteristic of the martensitic transformation. An additional deformation rate is introduced into the constitutive model to take into account this strain softening. The total deformation rate can be written as the sum of an elastic, a plastic, and a transformation part: The elastic properties of all phases are essentially the same and standard isotropic linear hypoelasticity of homogeneous solids is used to define the elastic behavior of the composite material. The methodology presented by Papadioti et al. [26,70] is used to describe the plastic part D p . The transformation part D TRIP has both deviatoric and volumetric parts and is in proportion to the rate of change of the volume fraction of martensite during martensitic transformation. Herein, the transformation kinetics model developed by Haidemenopoulos et al. [17] is used to describe the evolution of the volume fraction of martensite.

The Elastic Part D e of the Deformation Rate
The constitutive equation for D e is expressed as where σ is the Jaumann derivative of the stress tensor σ, M e is the elastic compliance tensor defined as where G and κ stand for the elastic shear and bulk moduli, δ and I for the second-and symmetric fourth-order identity tensors, with Cartesian components δ ij (the Kronecker delta) and I ijkl = (1/2)(δ ik δ jl + δ il δ jk ).

Yield Criterion
The TRIP steel is treated as a rate-independent composite material with a yield function of the form (Papadioti et al. [26,70] where σ is the stress tensor, s = σ − p δ the stress deviator, p = σ kk /3 the hydrostatic stress, and σ e the von Mises equivalent stress defined as c (i) (i = 1, 2, 3, 4) the volume fractions of the phases, andε (i) the average equivalent plastic strain in phase i. The effective yield stressσ 0 depends on the hardening behavior and the volume fractions of the constituent phases and is calculated from a constrained optimization problem: where N = 4 is the number of phases, and σ (i) 0 (ε (i) ) are the flow stresses of the individual phases. The methodology developed by Kaufman et al. [71] and the CONMAX software (http://www. netlib.org/opt/conmax.f) are used to solve the constrained optimization problem in (30).

The Plastic Part D p of the Deformation Rate
The plastic part of the deformation rate D p is determined from "normality" to the yield surface: whereε p is the non-negative equivalent plastic strain. Homogenization theory also provides estimates for the average value of the deformation rate D (i) in the individual phases, which are written in the form (Papadioti et al. [26,70]) whereŷ (r) are the optimal values of y (r) calculated from the optimization problem in (30).
The average equivalent plastic strain rate in the phasesε (i) is defined aṡε Then, (32a) yields˙ε The average equivalent plastic strain in the phasesε (i) , which is calculated from the time integration ofε (i) , is used to describe the hardening of the phases. In particular, the flow stress σ (i) 0 of each phase is assumed to be a function of the corresponding average equivalent plastic strain in the phase, i.e., σ

The Transformation Part D TRIP of the Deformation Rate
Finally, the transformation part D TRIP is defined as (Stringfellow et al. [18]): f = c (4) denotes the volume fraction of martensite, a superposed dot stands for the material time derivative, A 0 , A 1 and ∆ v are dimensionless parameters, and s * a is a reference stress. In (34a), the first term is deviatoric and accounts for "shape changes" caused by martensitic transformation, and the second term is volumetric (D TRIP kk =ḟ ∆ v ).

The Total Inelastic Deformation Rate D in ≡ D p + D TRIP
The inelastic deformation rate is defined as the total of the plastic and transformation parts: so that D = D e + D in .
Use of the constitutive equations for D p and D TRIP leads to

Evolution of the Volume Fraction of the Phases
According to the transformation kinetics model presented in Section 4.3, the evolution equation for the volume fraction of martensite can be expressed aṡ where c (a) = c (3) denotes the volume fraction of austenite,ε (a) =ε (3) is the average equivalent plastic strain in the austenite, and A f is determined by Haidemenopoulos et al. [17] A f ε (a) , where N ε0 v is determined in Equation (23).

The Elastoplastic Tangent Modulus
An equation relating the Jaumann derivative ∇ σ to the deformation rate D through the elastoplastic tangent modulus tensor L is derived from the elastoplastic constitutive equations. The derivation is as follows.
The elastic deformation rate D e is written as follows Substitution of the above expression for D e into the hypoelastic constitutive Equation (26) where L : N = 2 G N and L : δ = 3 κ δ were taken into consideration. Φ is an isotropic function, therefore the "consistency condition"Φ = 0 can be expressed as [72] Φ = ∂Φ ∂σ : Substitution of (42) for ∇ σ and of the evolution equations ofε (i) andċ (i) into the consistency condition leads to where N : L = 2 G N, N : N = 3/2, N : δ = 0, andε p v = mε p were taken into consideration. Last equation leads to˙ε where Substitution of the expression forε p from (46) into (42) leads to The numerical integration of the constitutive equations in the context of a finite element formulation is discussed in Appendices A and B.

Comparison of the Constitutive Model with Experimental Data
The flow stresses in the phases are written as follows where r denotes the number of the phase, σ (r) y the yield stress, andε (r) the equivalent plastic strain. The hardening behavior of the phases is determined through detailed literature search. The hardening behavior of ferrite, bainite, and martensite were acquired from experimental data presented in "Technical Steel Research" [73]. In particular, data for the annealed ferritic steel DOCOL 600 were used to determine the flow curve σ  (51) below, based on the following reasoning. In the production of TRIP steels, the final step that follows intercritical annealing is isothermal holding in the bainite transformation range. During the formation of bainitic ferrite, carbon is rejected to the retained austenite, and the carbon content of the retained austenite is raised to values above 1 wt% [74]. This provides chemical stabilization and raises the austenite yield strength considerably; values in the range of 500-550 MPa were reported [16,75]. Therefore, the stress-strain curve of retained austenite is expected to lie above that of ferrite, which exhibits a lower yield strength due to its very low carbon content. The following expression is used for the flow curve austenite at 25 • C: (51) Figure 10 illustrates the hardening curves of the individual phases as given by (50) and (51).

Experiments
The mechanical properties of the TRIP material were evaluated via uniaxial tension tests. An INSTRON 8801 servo-hydraulic machine with 100 kN load capacity was employed for the tensile tests. Mechanical properties were determined according to the ASTME8M at a constant crosshead velocity of 0.5 mm/min. Two specimens were tested in the longitudinal (L) direction using a 25-mm gage length clip-on extensometer. The geometry of the specimen used for the tensile tests is depicted in Figure 11.
The mechanical properties derived from the stress-strain curve are presented in Table 3. The yield stress and ultimate tensile strength is 530 MPa and 762 MPa, respectively. A high strain hardening rate at early stages of plastic straining is observed. The stable RA microstructure results in a gradual austenite transformation [76] and progressive hardening until necking at a uniform elongation of 25.4% [77,78].  According to the simulation results of the heat treatment presented in Section 2.3, the initial volume fractions for the individual phases are c (1) = 0.507, c (2) = 0.327, c (a) = 0.166, and c (m) = 0.0.
The constitutive model was used together with the ABAQUS general purpose finite element code to study the problem of uniaxial tension and the results are compared with the corresponding experimental data. One four-node isoparametric axisymmetric finite element (CAX4H in ABAQUS) is used to solve the uniaxial tension problem.
In the calculations, the Young's modulus is E = 220 GPa and the Poisson ratio is ν = 0.3. The relative volume change related to the martensitic transformation ∆ v , used in the constitutive equation for D TRIP (Equation (34)), is taken to be ∆ v = 0.02.
A very important part of the problem is the evolution of the volume fraction of martensite f during the uniaxial tension test. It is well known that the f -ε curve (ε = uniaxial strain) has a sigmoidal shape (Olson and Cohen [79]). The values of the parameters in the transformation kinetics model were defined in Section 4.3. Figures 12 and 13 display the calculated f -ε and σ-ε curves and the corresponding experimental data, where σ is the nominal stress. The sigmoidal shape of the strain-induced transformation is predicted quite well by the transformation kinetics model. Initially, the transformation rate d f /dε increases with strain, then reaches a rather constant rate and finally decreases at higher strains as saturation is achieved. Finally, for the temperature considered, the saturation level is lower than the value of 1, which corresponds to complete transformation. The model predictions agree well with the experimental data.

Necking of a Bar
The constitutive model developed for the TRIP steel is used together with the finite element method to study the development of a neck in an axisymmetric tension specimen with a geometric imperfection. A cylindrical specimen with aspect ratio L 0 /R 0 = 3 is considered, with 2 L 0 being the initial length and R 0 the initial radius. A cylindrical system is introduced as it is shown in Figure 14. Symmetric solutions are considered and only one half of the cylindrical specimen is studied. Figure 14 illustrates the finite element mesh used in the calculations; it comprises 1350 four-node isoparametric axisymmetric elements (CAX4H in ABAQUS) in a 15 × 90 grid. The following geometric imperfection is introduced to encourage necking: where R(z) is the perturbed radius and ξ is set equal to 0.005. Figure 14 shows the imposed geometric boundary conditions. The deformation is driven by the uniform prescribed end displacement in the z-direction on the shear-free end z = L 0 ; the lateral surface on r = R 0 is kept traction free. The initial values of the volume fractions of the constituent phases are: c The uniaxial stress-strain curves of the imperfect specimen for both a transforming and a non-transforming material are depicted in Figure 15. The points of maximum load, which coincide with the end of uniform elongation in the corresponding specimen are denoted by the arrows. For the TRIP steel, the end of uniform elongation occurs at a nominal strain of 18.4% and 750 MPa stress, while for the non-transforming steel at 16% and 705 MPa respectively. Figure 15 shows that the TRIP phenomenon, in addition to strengthening the material, increases considerably the range of uniform elongation.   Figure 16 depicts the evolution of the radius at the minimum cross section of the specimen for the transforming and non-transforming materials, and Figure 17 shows the corresponding deformed configurations. At a nominal strain of 35%, the minimum cross section in the TRIP steel contracts 39.2%, whereas in the non-transforming material it contracts 44.7%. The formation of martensite stabilizes the neck and leads to its propagation down the length of the specimen. This result is consistent with the findings of Papatriantafillou et al. [19,80].

Forming Limit Diagrams
"Forming limit diagrams" for the TRIP steel under investigation are calculated using the constitutive model presented in the previous sections. Forming limit diagrams show the maximum deformation a sheet metal can be subjected to before the material fails by flow localization in a narrow straight band. Calculations are also carried out for a non-transforming steel with the same initial values of the volume fractions of the phases.
A sheet made of TRIP steel is deformed uniformly on its plane in a way that the in-plane principal strain increments increase proportionally. The possibility of the formation of an instability in the form of a narrow straight band (Figure 18) is studied. As discussed in Section 5.7, the constitutive equations can be expressed as (Equation (48a) where L is the elastoplastic tangent modulus tensor defined in (48b). The formulation of the problem is more straightforward if "nominal" quantities are used. It can be readily shown that using the 1st Piola-Kirchhoff stress tensor t = J F −1 · σ, the rate-constitutive Equation (53) can be written as˙t A state of uniform plane stress is assumed inside and outside the band and the possibility of the formation of a neck (band) as shown in Figure 18 is investigated. Let X 1 − X 2 be the plane of the sheet and H the initial thickness of the sheet. Greek indices (α, β, γ, δ) are introduced and take values in the range (1, 2). The in-plane displacements are continuous, therefore their spatial derivatives parallel to the band remain uniform. Let ∂u α /∂X β be the displacement gradient discontinuities across the band, where [ ] denotes the difference of the field within the band and outside the band, u α is the α-component (in-plane) of the displacement field, and X β the β-component of the position vector X of a material point in the undeformed configuration. It is a well known result in Mechanics of Materials that the only discontinuities in the displacement gradient are restricted kinematically to the following form (Hadamard [81], Hill [82], Rice [83]) where N β is the β-component of the unit vector N normal to the band in the undeformed configuration, and G α the α-component of the vector G that defines the "jump" in the normal derivative of u, i.e., [∂u/∂N] ≡ [(∂u/∂X) · N] = G. The vector G takes a constant value within the neck and depends on the imposed uniform deformation outside the neck. Next, a methodology for the calculation of G is presented. Equation (54) leads to the conclusion that the in-plane components of the deformation gradient inside the band can be written as F b αβ = F αβ + G a N β , with superscript b denoting quantities inside the band, whereas quantities with no superscript refer to the uniform field outside the band. The deformation gradients in a matrix form are where λ i are the stretch ratios. The plane stress condition σ 33 = 0 insinuates t 33 = 0. The conditioṅ t 33 = 0 is solved forḞ 33 , so that the in-plane constitutive relations required for the sheet necking analysis are written in the form: Similarly, the in-plane constitutive relations within the band arė Equilibrium across the neck requires that The rate form of this equilibrium relationship iṡ Substitution of Equations (56) and (57) in (59) and use of Equation (55) for F and F b leads to the expression where Equation (60) expresses the rate of equilibrium equation across the band, i.e., the rate of (58), in terms ofĠ. In (60), the jump G in the normal derivative of u across the band is determined in terms of the imposed uniform stretching through b (b α = λ α ).
In the case of a perfect sheet (H b = H), the right hand side of (60) vanishes, since H b = H and C b = C lead to B = 0, and the deformation continues to be homogenous (Ġ = 0) until the condition for local necking bifurcation det[A] = 0 is met.
The methodology of Marciniak and Kuzyski [84] is used in the calculations. A small initial imperfection is assumed to exist in the sheet and gradual localization of the strains at the imperfection leads to necking. In particular, a straight narrow band of reduced thickness H b < H is considered (Figure 18). A state of uniform plane stress inside and outside the band is assumed and the problem is to determine the uniform state of deformation within the band that is conforming kinematically and statically with the prescribed uniform state outside the band (Tvergaard [85,86], Needleman and Tvergaard [87]). When the sheet is not perfect (H b < H), the right hand side of (60) does not eliminate (C b = C), so (60) can be solved forĠ. Since the initial sheet thickness inside and outside the band, and the imposed uniform deformation history F outside the band are known, Equations (60) are solved incrementally for ∆G =Ġ ∆t to acquire the deformation history within the band. Localization takes place when the ratio of some scalar measure of the amount of incremental straining within the band to the corresponding value outside the band becomes unbounded.
The deformation gradient outside the band F is prescribed in a way that the principal logarithmic strains ε 1 and ε 2 outside the band increase proportionally: The plane stress algorithm discussed in Appendix B is used to acquire the uniform solution outside the band. At the end of each increment, ∆G is calculated using Equation (60) and subsequently the deformation gradient within the band F b is defined. Next, the plane stress algorithm is used again to acquire the uniform solution within the band. The localization condition is met when d|G|/dλ 1 = ∞, i.e., when det[A] = 0 within the band.
To increase the precision of the calculations, instead of the rate of equilibrium Equation (60), equilibrium itself (58) is used: where the subscript n + 1 denotes values at the end of the increment. Then, the approximations T n+1 T n +Ṫ n ∆t = T n + H N ·ṫ n ∆t and T b n+1 T b n +Ṫ b n ∆t = T n + H N ·ṫ b n ∆t are used in (61). Finally, (56) and (57) are used forṫ n andṫ b n , and (55) is used for F and F b , to find which is used for the determination of ∆G instead of (60). The last term on the right hand side of (62) takes into consideration any unbalanced forces at end of the previous increment.
The initial values of the volume fractions of the constituent phases are: c The unit vector N can be expressed as N = cos Ψ e 1 + sin Ψ e 2 , with Ψ being the angle of inclination of the band relative to the X 1 axis in the undeformed configuration. For each value of ρ = dε 2 /dε 1 , calculations are conducted for values of Ψ covering the range 0 • ≤ Ψ ≤ 90 • and determine the strain level at which the localization condition det[A] = 0 is first met. The critical value Ψ cr for each ρ is defined as the one resulting in the minimum localization strain. Figure 19 depicts "forming limit curves" for proportional straining for a case without imperfection (H b /H = 1) and for two different values of the initial thickness imperfection, namely H b /H = 0.999 and H b /H = 0.99. The three solid curves denote the TRIP steel, while the dashed curves correspond to the non-transforming steel. It is observed that the TRIP phenomenon increases in general the necking localization strains. This result is consistent with the findings of Papatriantafillou et al. [19,80], who used a rate dependent constitutive model for TRIP steels (as opposed to the rate independent model used here). For the case of plane strain (ρ = 0) and no imperfection, the critical strain ε cr

Conclusions
The mechanical properties of TRIP steels are intrinsically associated with the stability of retained austenite, which results from the heat treatment design. In the current work, models for the heat treatment of TRIP steels, the austenite stability, the transformations kinetics of austenite as well as the mechanical behavior of the composite material were developed. The predictions of the models were then verified with an experimental study.
In particular, a 2D multi-phase field (MPF) model was employed for the prediction of the microstructural features of a CR-TRIP700 steel during a two-stage heat treatment, consisting of intercritical annealing, followed by an isothermal bainitic treatment. The MPF model is able to describe the temporal evolution of the volume fractions and the average grain size of the phases, as well as their average concentration in carbon, aluminum, and manganese. The phase-field results, obtained at the end of the heat treatment, were implemented in the M σ s temperature model. Both the chemical enrichment and the size refinement of the retained austenite resulted in a rather stable retained austenite dispersion.
An experimental validation of the heat treatment process was performed. The measurements concerning the volume fraction and the average grain size of retained austenite were compared to the respective phase-field results. A good agreement between the experimental data and the phase-field simulation results was observed.
The M σ s temperature model was employed to predict the stability of the retained austenite. In order to take advantage of the strain-induced transformation occurring in TRIP steels, M σ s temperature should be bellow room temperature. The value of the M σ s temperature predicted by the model is consistent with the experimental measurements of the SS-TV-TT technique. The yield strength, the mean particle size, and the stress state influence austenite stability. An increase in any of those parameters results in an increase in the M σ s temperature, leading to a decrease in the stability of the retained austenite.
The phase-field results, obtained at the end of the heat treatment, were also implemented in a model for the transformation kinetics of retained austenite. The sigmoidal shape of the strain-induced transformation was predicted quite well by this model. The transformations kinetics model and non linear homogenization techniques for non-linear composites were used to develop a constitutive model to describe the mechanical behavior of the TRIP steel under investigation. A method for the numerical integration of the constitutive model in the context of a finite element formulation was presented, and the model was introduced in a general-purpose finite element code (ABAQUS). In addition, a methodology for the numerical integration of the constitutive equations under plane stress conditions was discussed. One-element finite element calculations for the uniaxial tension problem were conducted, and the results were compared with experimental data obtained from uniaxial tension tests. The model predictions agree well with the experiments. The problem of necking under tension was analyzed thoroughly, and "forming limit diagrams" were calculated. In both cases it is evident that the TRIP phenomenon not only strengthens the material, but also increases considerably the range of uniform elongation.
The models developed provide an integrated simulation toolkit for the computer-assisted design of TRIP steels, which can be used to translate mechanical property requirements into optimised microstructural characteristics and to identify the appropriate processing routes. This methodology can be employed for other steel grades used in the automotive industry with the appropriate calibration. Funding: This work was carried out in the framework of "TRIP steels: Simulation and experimental study of heat treatment and mechanical behavior" project. This research is implemented through the Operational Program "Human Resources Development, Education and Lifelong Learning" and is co-financed by the European Union (European Social Fund) and Greek national funds.

Acknowledgments:
The authors would like to thank G.N. Haidemenopoulos and A.T. Kermanidis of the University of Thessaly for many helpful discussions and material. The authors would also like to acknowledge the help of D. Krizan from the R&D department of voestalpine Stahl GmbH Linz, Austria, for providing the TRIP700 material and carrying out the saturation magnetization measurements.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Numerical Implementation of the Constitutive Model
In a finite element environment, the solution is developed incrementally, and the constitutive equations are integrated at the element Gauss integration points. At a Gauss point, the solution F n , σ n , c (r) n ,ε (r) n at time t n as well as the deformation gradient F n+1 at time t n+1 = t n + ∆t are known, and the problem is to determine σ n+1 , c (r) n+1 ,ε (r) n+1 . Let A n and A n+1 stand for the values of A at the start t n and at the end t n+1 of the increment; the notation ∆A = A n+1 − A n is used.
The effective flow stressσ 0 (c ) of the TRIP steel is assumed to remain constant over the time increment (t n , t n+1 ) and can be calculated from the optimization problem in (30), in which the flow stresses σ (i) 0 of the phases take values (Papadioti et al. [26,70]) The optimal values ofŷ (i) n+1 , which are used to define the strain concentration factors α (i) in (32b) for the increment, are also calculated from the optimization problem (30). The calculation is implicit in general, except when β = 0 is used in (A1).
In the following, an algorithm suitable for three-dimensional, plane strain, and axisymmetric problems is presented.
Equations (A4) that determines the inelastic deformation rateĖ in , (A5) that define theε (i) , and (A9) that defines the evolution of the volume fraction of martensite f , all require numerical integration. The remaining equations are integrated exactly: whereσ e = σ n + L e : ∆E is the (known) "elastic predictor". The rest of the equations are The backward Euler method is used for the numerical integration of the inelastic deformation rate (A18), and the forward Euler scheme is used for the numerical integration of (A19) and (A20). Use of a backward Euler scheme for the numerical integration ofĖ in is imperative in order to be able to use increments of reasonable size (Aravas and Ponte Castañeda [88]). Numerical integration The non-linear Equation (A28) for ∆ε p is solved using Newton's method. After ∆ε p was found, ∆ε are calculated by using (A29)-(A34), σ e | n+1 is found from (A27), A n+1 = A 0 + (A 1 /s * a ) σ e | n+1 is calculated,σ n+1 is determined from (A24), and the integration is completed with the calculation of σ n+1 = R n+1 ·σ n+1 · R T n+1 .

Appendix B. Plane Stress Algorithm
In this appendix, a plane stress algorithm for the numerical integration of the constitutive equations is developed. In plane stress problems, the out-of-plane component of the deformation gradient is not defined kinematically and the algorithm described in Appendix A needs to be altered. The geometry analyzed is a thin plane disc of uniform thickness loaded in its plane. The mean plane of the disc coincides with the plane X 3 = 0, and the in-plane displacements are of the form The deformation gradient and the stress tensor, in isotropic materials, take the following form where (e 1 , e 2 , e 3 ) are unit vectors along the coordinate axes and the summation convention on repeated Greek indices (α, β) is used. (α, β) take values in the range (1, 2). In finite strain problems, when the in-plane displacement field is inhomogeneous, the out-of-plane displacement and the corresponding thickness variation are functions of (X 1 , X 2 ). It is assumed that, during material deformation, the resulting thickness variation is negligible, therefore the plane stress assumption is reasonable and (A36) are valid to a good approximation.
In plane stress problems, the method described in Appendix A needs to be modified since the out-of-plane component of the deformation gradient F 33 is not defined kinematically. The deformation gradient ∆F n+1 = R n+1 · U n+1 associated with the current increment is written in the form with ∆F αβ (α, β = 1, 2) being the known in-plane components, and ∆F 33 the unknown out-of-plane component. Likewise, the associated right stretch tensor U n+1 and the orthogonal tensor R n+1 from the polar decomposition of ∆F n+1 are expressed as and the logarithmic strain tensor E n+1 = ln U n+1 associated with the increment is expressed as where bared quantities ∆Ē αβ are known, and ∆E 3 is the unknown out-of-plane component of E n+1 . It is worthy of note that the only difference with the method described in Appendix A is that the values of ∆F 33 and ∆E 3 are unknown when the process of the numerical integration starts. The plane stress condition is used to calculate the value of ∆E 3 : The logarithmic strain tensor associated with the increment is expressed as where ∆Ē = ∆Ē αβ e a e β is the known in-plane part of E, and a is the deviatoric part of a: a = a − 1 3 a kk δ = − 1 3 (e 1 e 1 + e 2 e 2 − 2 e 3 e 3 ) .
Equations (A11)-(A20) are now expressed as follows: where bared quantities are known, andσ e = σ n + L e : ∆Ē is the "elastic predictor" corresponding to the known part of ∆E.