Strain simulation of steel during a heating-cooling cycle including solid-solid phase change

Designed for the simulation of transformation induced plasticity in steels occurring on cooling, the model ﬁrst proposed by Leblond in the middle of 80’s works well for constant uni-axial stress test. But further experimental investigation have shown that for non-proportional multi-axial test, the quality of the model is worse. The purpose of this paper is to improve this model through the removal of the hypothesis that plastic ﬂow for the Greenwood–Johnson mechanism occurs only in the ‘weaker’ phase. The improvement, leading to better agreement, enables to use the model on cooling, and shows that for TRIP plastic ﬂowing in harder phase is almost important, even if it is lower than in the ‘weak’ phase


Introduction
The model presented below is applied to the simulation of transformation induced plasticity (TRIP) for low carbon steels.Material behavior taken into account are elastic with either ideal plasticity or linear kinematic hardening.The focus of this model is to relate stress, temperature and phase proportion to strain during the transformation.
For steels undergoing a solid-solid phase transformation, TRIP appears.This phenomenon is a plastic flow occurring when an external load is applied during the transformation.This occurs even if this load is small regarding yield stress of the 'weaker' phase.Greenwood and Johnson (1965) and Magee (1970) give two complementary explanations to this phenomenon: -orientation of the local plastic flow due to phase volume incompatibility by external loading (Greenwood and Johnson mechanism); -preferred orientation of the martensite plates arising from the external loading (Magee mechanism).
From a metallurgical point of view, Greenwood and Johnson mechanism is related to diffusional transformation and Magee effect is related to displacive transformation occurring for martensitic transformation involving Bain strain.
Of course these mechanisms are not exclusive, but experimental works show preponderant relation between them (Bhadeshia, 1995).
Modeling of Greenwood and Johnson transformation plasticity phenomenon has mainly been developed by Leblond et al. (1986aLeblond et al. ( , 1986bLeblond et al. ( , 1989bLeblond et al. ( , 1989a) ) and by Fischer (1990).According to recent experimental work of Coret et al. (2003) and to analytical work of Taleb and Sidoroff (2002), Leblond model give accurate results.However this model is not designed to describe the behavior of the 'harder' phase.Therefore, the goal of the following work is to up grade the so called Leblond model by taking into account the possibility of plastic flowing in the classical named 'harder' phase.Indeed, in agreement with Greenwood and Johnson (1965) work on TRIP, Leblond et al. (1986b) has removed the possibility of plastic flow of the 'harder' phase.
Since model used is similar to Leblond's one, we will first present how to relate macroscopic (structure) scale to microscopic scale (grain scale) in the sense of Leblond.Then we will present the micromodel which enable to take plasticity into account in the harder phase.Application to ideal elasto-plastic case and how to mix it up with the precedent approach is done.Then, an extension to the linear kinematic hardening case is also introduced and at last, a comparison with experiments is done.

A summary of the Leblond model
Since this paper is strongly linked to Leblond model, this part presents the main feature of the model (Leblond et al. 1986a(Leblond et al. , 1986b(Leblond et al. , 1989b(Leblond et al. , 1989a)).We will take for instance a transformation occurring on cooling from a weak phase γ to a hard phase α, the notation λ will be reserved for denoting either the α or γ phase.In the whole paper z will denote the α phase proportion.In his model, Leblond et al. (1986a) uses a macroscopic scale, capital notations in this paper, and a microscopic scale, lowercase notations.The macroscopic one enables to separate the contribution of different phenomenon that are: -classical plasticity related to mechanical load variation; -classical plasticity related to thermal load variation; -transformation plasticity related to phase proportion variation implying Greenwood and Johnson and Magee effects.
Let us now introduce hypothesis: Hypothesis 1.The microscopic elastic compliance tensor m may be equated to the macroscopic overall elastic compliance tensor M.
Under Hypothesis 1, one obtains macroscopic plastic strain rate Ėp decomposition given by Eq. (2.1), into contribution of classical plasticity strain rate due to external load Σ (stress tensor) variation, Ėcp Σ , classical plasticity strain rate due to temperature variation Ṫ , Ėcp T , and transformation plasticity due to phase variation, Ėtp .Next, from a thermodynamical approach, we could get homogenization relations between macroscopic strain rate and microscopic plastic strain p variation, given by Eq. (2.1).
The last term of Ėtp , which is an average on the front, is linked to Magee mechanism and the two first, average on α, V α , and γ , V γ , phase volume, explain the Greenwood and Johnson effect.Since the model is only interesting in diffusional transformation, the Magee mechanism will be neglected for the model.Next step is the removing of 'weak' phase terms under Hypothesis 2: Hypothesis 2. For small or moderately high applied stresses, the γ phase is entirely plastic, but the α phase remains elastic, or its plastic strain rate remains always much smaller than the γ phase's one.
This hypothesis first appears in Greenwood and Johnson (1965).It is justified by the fact that yield stress of γ phase is much smaller than that of α phase.But, as will be discussed in Section 3, this justification could not be sufficient to remove this term.Finally, these hypothesis lead to Eqs. (2.3): (2.3) From this point, Leblond begins a homogenization step to replace the preceding set of equations by a new one which only depends on the microscopic equivalent strain in order to get useful equations for the micro-mechanic model.
Hypothesis 3.Both phases obey the von Mises criterion and the Prandtl-Reuss flow rule.
Hypothesis 4. Interactions between δ eq γ /δΣ ij and s γ which is the microscopic stress deviator in γ phase, can be neglected.
Hypothesis 5.For small applied stresses, the macroscopic stress deviators S γ and S α in phase γ and α are almost equal and identical to the total macroscopic stress deviator S, i.e., the average value of the microscopic deviator over the whole representative volume V : Hypothesis 6. Σ eq = ( 3 2 S ij S ij ) 1/2 denoting the macroscopic von Mises equivalent stress, Ėcp Σ is non-zero only if Σ eq varies, i.e.: V γ S Σeq . (2.4) The last set of equation enables to relate macroscopic strain rate and microscopic strain variation.It remains to evaluate microscopic strain variation.It could be done by submitting a representative local zone to a variation of z, T , and Σ .
Since z and T are scalar quantities, the micro-model used to estimate Ėtp and Ėcp T is based on a spherical geometry (Leblond et al., 1989b).On the other hand, since Σ eq is related to a tensorial quantity, a shearing micro-model and a traction micro-model are introduced (Leblond et al., 1986b).Micro-model for the estimation of δ eq γ δz V γ and δ eq δT V γ .Micro-model (cf.Fig. 1) to estimate is made of two inclusive spheres of each materials which radius is evolving proportionally to phase proportion.Mechanical problem to solve on these micro-model is defined through hypothesis and boundary conditions: Hypothesis 7. In spherical micro-models, hypothesis used are listed below: -external stress is neglected in comparison to the stress induced by transformation induced plasticity; -both phases are expected to be ideal elasto-plastic materials.
Boundary conditions 1. Boundary conditions used to obtain eq in the micro model are: -behavior at the center is smooth.Hence, no plastic flow occurs in the α phase without significant external load; -continuity of the radial displacement on the boundary α − γ p and γ p − γ el ; -continuity of the radial stress on the boundary α − γ p and γ p − γ el ; -continuity of the circumferential stress on boundary α − γ p ; -null radial stress on the external boundary of the γ phase.
Solving the mechanical problem using equilibrium, the variation of local equivalent strain δ eq in the plastified zone writes: where x denote the radial coordinate and, th γ α = th α − th γ is the difference of thermal strain between the two phases.As the micro-model is assumed to be representative of a local zone, it allows access to homogenized quantities through integration.Taleb's work (Taleb and Sidoroff, 2002) on the plastification radius for the TRIP term, enables to reduce the integration on the only plastified zone, so that one gets Eqs.(2.6) and (2.7) for TRIP strain rates.
with α λ the expansion coefficient of the λ phase, and K stands for the elastic bulk moduli.

Micro-model for the estimation of
Leblond proposes from a systematic numerical analysis (Leblond et al., 1986b) the following model for the classical plasticity strain rate at constant stress: Of course the g function is an approximation, but it appears that this term accounts only for a small part of the plastic strain rate and does not need an accurate estimation.

Symmetrical extension of Leblond model in ideal elasto-plastic case
We now propose a symmetrical extension of the Leblond model taking into account the possibility of plastic flow in the harder phase.The main idea here is simply to consider that not only the α phase could be totally include into the γ phase, but that the reverse situation can also occur.When this situation happens, due to the isotropy of the local problem, the γ phase could not sustain plastic flow.Hence, in this case, plastic flow is reported to harder phase.We shall discard Hypothesis 2, and even if the γ phase elastic yield limit is small with respect to that of α phase we shall take into account a possible plastic flow in the α phase.

Treatment of the 'harder' phase associated macro term
First homogenization step from the Leblond is kept, we only do not remove terms associated to the α phase.the set of equations obtained without mathematical difficulties is presented below: (3.1) We hence obtain a set of homogenized equation able to take into account plasticity in the 'weaker' phase, so that it only remains to determine the relevant micro-model for that part.Since we are keeping the two terms of plastic strain rate in each phase, we have to introduce two micro-model called γ and α configuration for calculating respectively plastic flow in γ phase and in α phase.In the same manner, we denote by v κ λ the volume of the λ phase in the κ configuration, where κ is either α or γ .

Inversion of the micro-mechanical model for the TRIP term
The micro-model used for the estimation of the equivalent strain in the α phase is displayed on Fig. 1(left).It's main difference with Leblond one (Leblond et al., 1989b) is that the γ phase is now included within the α phase.For this micro-model, even if yield stress of the γ phase is lower than that of the α phase, plastic flow can occur only in the α due to the necessary isotropic volume increase of γ phase.
Hypothesis 7 and Boundary conditions 2, are used to define the mechanical problem to solve on the configuration, which remain in the same spirit than the preceding case.
Boundary conditions 2. The boundary conditions used to obtain the equivalent strain of the inverted micro-model are: -behavior at the center is regular.Hence, no plastic flow occurs in the γ phase without significant external load; -continuity of the radial displacement on the boundaries γ − α p and α p − α el ; -continuity of the radial stress on the boundaries γ − α p and α p − α el ; -continuity of the circumferential stress on the boundary γ − α p ; -null radial stress on the external boundary of the α phase.
For this micro-model, Leblond et al. (1989b) andTaleb's (Taleb andSidoroff, 2002) work remain applicable.We can use it to get first a plastification radius, next the equivalent plastic strain rate: Let us underline that the inversion of phase relative position dot not change the growth of plastic zones which always start to develop from the inner towards the outer radius, due to the constant sign of the difference in compactness γ α .

Inversion of the micro-mechanical model for the temperature classical plasticity term
The situation is opposite for this part of the strain rate.We have now to introduce a new hypothesis, because the plasticity now starts from the outer towards the inner radius.We ensure that the two phase remain a continuous medium through assuming that the two phases are perfectly bonded.Under this hypothesis, and working like in the standard case, we obtain: (3.3)

Inversion for the mechanic load classical plasticity term
For Ėcp Σ , a different point of view is needed.This term is not only related to isotropic volume variation.It is then not obvious that any plastic flow will occur in the harder phase.Plastic flow may occur in a totally included phase if the mechanical load is not isotropic.Considering this fact it was decided to keep an elastic behavior of the α phase in case of variation of external mechanic load.

Consequence for macroscopic terms
All preceding results can now be combined to obtain the formulation of the macroscopic terms.We consider a representative macroscopic zone that contains the two phases α and γ .These phases are represented through a set of N − n configuration α, and n configurations γ .The problem is to get a model which respects the global phase proportion independently of the number n and N .For that, we will construct configurations that always respect the global proportions, so that, with V λ the whole volume of λ phase: Since we consider an infinite number of configurations, the change from the γ dominant to a α dominant configuration case can be considered in a continuous way.Introducing V λ λ to denote the whole volume of the λ in the λ configuration, we obtain: Finally, using the preceding mixture rule, we obtain the macroscopic strain variations:

Using the inversion for simulating heating transformation
Theoretically, the Leblond model was not design to simulate a transformation on heating.Indeed, the main part of the plastic flow is considered to happen at the beginning of the transformation.So, for modeling TRIP phenomenon on heating, it is necessary to consider plastic flow in the α phase, exactly what has been done here.Hence, assuming that interface between α and γ phase is perfectly bonded, we obtain evolution equations:

β(z) law for the change of configuration
The symmetrical extension of the Leblond model introduces a new undetermined function, β(z).This law is directly related to the topology of the phase transformation.We could then imagine three main steps, first, at beginning child phase is completely include into parent phase, so β is valued to one.Second, at last, the parent phase is completely include into child phase, so β vanish.Between these two steps, and behind the lack of experimental observation we take a linear law function of the phase proportion, as illustrated on Fig. 2. So the determination of the β(z) law is replaced by the determination of two proportions limits, z 0 and z 1 .These limits represent phase proportion such that one phase is totally included in the other one.In order to give an estimation of z 0 and z 1 , we use simplified polyedric geometry.The growth of the other phase is simulated by sphere on the corners.We consider that one phase is included into the other when the radius of these sphere is such that they are in contact, as shown on Fig. 3.This approach is simplified, but give an approximation for z 0 and z 1 .Results have been gathered in Table 1.For the simulation, we take z 0 = 74.05%,like in the tetraedra case, and assuming that there is no reason that behavior acts differently in the other direction, we take z 1 = 1 − z 0 = 25.95%.

Coupling with linear kinematic strain hardening phenomena
Coupling with kinematic hardening is relevant for such a model since experimentation conducted by Coret et al. (2003) under load reversals shows a typically hardening behavior.The main difficulty here is the estimation of the variation of the hardening parameters.Leblond's work (Leblond et al., 1989a) can be used to predict variation of strain rate.

Evolution equations for transformation and classical plastic strain
The determination of the plastic strain rate remains unchanged.We denote by A κ λ the center of the macroscopic yield surface of the λ phase in the κ configuration.Reasoning in the same manner as Leblond et al. (1989a), we obtain constitutive equations: (5.2) (5.4)

Evolution equations for the hardening parameters
The last step is the estimation of the homogenized hardening parameter evolution for linear kinematic hardening during phase transformation.Each configuration λ undergoes a different history.Hence a fictitious front appears between the α and γ configurations.Expression of Ḃλ , the homogenized hardening parameter in the phase λ, will be obtained as a function of ḃλ , the microscopic hardening parameter in the phase λ, by time differentiation.Hence, taking into account that the front between the two phases is moving: where F denote the transformation front and U n the normal velocity.Hence, we could get for each phase: In upper equations, F λ represents the transformation front between the two phases in the λ configuration.As said before, F * denotes a fictitious front between the γ configuration and the α configuration.On the other hand, we have the expressions: .
Since determination of other terms and of microscopic parameters is not of first interest, they are reported in Appendix A. We summarize here main results that are:

Heating case
Always keeping on mind that interface between α and γ phase is assumed to be perfectly bonded, and working in the same manner than above, we successively get constitutive equations: (5.7) and then evolution equations: (5.9)

Experimental comparison
The law we obtain will now be used to compare simulation and experimental results due to Coret et al. (2003).We shall present one cooling bi-axial, for a bainitic transformation, and one heating uni-axial example.Material and test parameters could be find in Appendix B.
For using the symmetric extension, and also Leblond's model for simulation, we should have access to the evolution of phase proportion during test.In that goal, phase proportion has been obtained from a free dilatometry test.Imposed load coming from test is then shifted so that beginning of the load imposition corresponds to the beginning of the phase change given by free dilatometry.

Bi-axial stress and cooling
This section compares results with a bi-axial stress case.During cooling the test bar is submitted to a constant tensile stress and a variable sinusoidal shear stress.Strain versus temperature is shown on Fig. 4. It can be observed in this figure that, for the longitudinal deformation the new extension give as good results as those obtained with the original Leblond's one but as far as the shear strain is concerned, the new extension give a better agreement with experiments.
The proposed β law gives a way to comment obtained results.We would focus on the shearing test, since it's the one which is the more different for the 2 models.The main difference we found is not on the beginning of the transformation, which seems quite reasonable, since the old and the new model response are identical at this time of the transformation.
As the transformation is going on, below 420 • C, the first difference appears on the maximal shear which agree better with experimental data.At that time the behavior of the two phase is take into account, which is certainly responsible for such a behavior.
After, below 350 • C, TRIP is going faster with the new model.This behavior is clearly due to plastification appearing in the harder phase and enable a better agreement with the final deformation.
For the new model, it could be found that characteristic points and tendencies are shifted in temperature.This could be due to some difference of evolution of phase proportion between the stressed test versus the free test giving the phase proportion for the simulation.

Heating case
In this case, the sample is submitted to a constant traction stress during transformation while heating.Comparison between experiment and symmetrical extension model are shown on Fig. 5.It can be seen that even if the value during the test do not agree perfectly with the experiment, the final response, with the ideal elasto-plastic behavior, is near the experimental result.

Conclusion
As it has been shown in this paper that the Greenwood-Johnson hypothesis about the non-plastification of the weaker phase, in a Leblond model point of view, can be removed.The two first part show how to do it and its implication on macroscopic model.Then, using a mixture rule, all preceding equations can be combined to introduce a model with linear kinematic hardening.This improved model leads to a better understanding of the TRIP phenomenon since it shows that the plastic behavior of the two phases must be take into account Preceding figures show that the proposed improvement, using the symmetrical extension, do not improve drastically the results.But, for the cases for which this extension is not very efficient, the results remain close to those obtained by classical Leblond model.When the results are different from the Leblond's one, the predictions by the symmetrical extension model give definitely more accurate results.Another interesting feature is that the symmetrical extension can be used for heating simulations.This feature is of interest for simulation of multi-pass welding and so on.Coupling with viscosity following Vincent (2002) has been tried but has not led to significant improvements.

Acknowledgement
Many thanks to Professor R. de Borst for his help to finalize the paper.

Appendix A. Determination of hardening parameters
Details of the time differentiation for each case are given by: To conclude, it is now necessary to estimate hardening parameter for each phase and each configuration at a microscopic level.
Hardening parameter for the γ configuration.Directly from Leblond's (Leblond et al., 1989b) work, we can write: -at a point which is not undergoing the transformation, said didn't belong to transformation front F γ , evolution law is take as usual ḃλ = ˙ Hardening parameter for the α configuration.
-at a point which is not undergoing the transformation, said didn't belong to transformation front F α , evolution law is take as usual ḃλ = ˙ Hardening parameter when changing of configuration (from γ to α).
-the γ phase in the γ configuration is not changed by the inversion of configuration b

Fig. 1 .
Fig. 1.Configuration γ and α: schema of micro-models for the estimation of plastic flow in γ and α phases.
the influence of the disappearing austenite on the position and velocity of transformation front b γ F γ = b γ V γ γ ; -memory effect is represented by a memory coefficient θ which simulate the fact that a part of hardening of disappearing austenite is given to new phase b α F γ = θ b γ F γ = θB γ γ .
the influence of the disappearing austenite on the position and velocity of transformation front b γ F α = b γ V α γ ; -memory effect is represented by a memory coefficient θ which simulate the fact that a part of hardening of disappearing austenite is given to new phase b α F α = θ b γ F α = θB α γ .
γ γ F * = b γ V γ γ ; -the α phase in the γ configuration is not changed by the inversion b γ α F * = b α V γ α ; -the γ phase in the α configuration get exactly hardening it had in the precedent configuration b α γ F * = b γ V γ γ ; -the α phase in the α configuration get exactly hardening it had in the precedent configuration b α α F * = b α V γ α .

Fig
Fig. B.1.γ phase yield strength (MPa) versus temperature from Vincent and Ahrens experimental data.

Fig
Fig. B.2. α phase yield strength (MPa) versus temperature from Vincent experimental data.

Table 1
Evolution of z 0 for three-dimensional geometry