An implicit integration scheme with consistent tangent modulus for Leblond’s model of transformation-induced plasticity in steels

Thermomechanical treatments involving solid-state phase transformations play an important role for the manufacturing of functional and reliable components in many engineering applications. Accordingly, numerical investigation and optimization of such processes require considering thermoelastoplasticity under the influence of ongoing transformations and in particular the impact of transformation-induced plasticity (TRIP). While a number of elaborate plasticity models have been proposed for the description of TRIP, none of them seem to have received much prevalence in applications due to their complexity or hard to determine model parameters. Instead, the overwhelming majority of applied research either relies on simplistic formulations dating back to early phenomenological approaches or neglects TRIP altogether. In this work, we therefore provide an accessible, straightforward and easy-to-implement solution scheme for the TRIP model proposed by Leblond et al. which, despite being widely recognized, is hardly ever employed in full form. Specifically, we employ implicit backward-Euler integration and an elastic–plastic operator split approach to update the stresses in order to obtain a simple and concise algorithm for which we then derive the corresponding consistent tangent modulus. Furthermore, the work contains an application of the solution scheme to a symmetrically cooled plate and an in-depth discussion of the influence of TRIP by means of this tractable numerical example. Specifically, we highlight the discrepancies arising in transient and residual stresses and strains compared to the conventional J2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2$$\end{document}-plasticity approach where the phase transformation is accounted for merely by adapting the yield strength of the compound.

may also arise in cases where the phase transformation is induced thermally and instead of external mechanical loading, internal stresses, which result from transient inhomogeneous temperature distributions and transformation progress, drive plasticity. This applies e.g. to certain heat treatment and welding processes and will be the focus of this work. With regard to the theoretical understanding and description of TRIP, Greenwood and Johnson [10] and Magee and Paxton [23] were among the first to publish pioneering mechanical interpretations and models to explain their experimental investigations and relate them to the microstructural mechanisms attributed to TRIP. However, initial efforts to model TRIP were restricted to uniaxial loading [1,10,23,26], while the advent of computational plasticity and its application in numerical simulations of industrial applications soon called for general, three-dimensional formulations. Franitza [8] was the first to consider TRIP in such an application by simulating residual stresses resulting from heat treatment using an ad hoc generalization of Greenwood and Johnson's model to three dimensions. Similarly, up until today the majority of the rather few authors who incorporate TRIP in simulations of thermo-mechanical-metallurgical manufacturing processes do so by imposing an additional inelastic transformation plastic strain rate of the forṁ tp = 3 2 κφ (z)żσ (1) while applying J 2 -plasticity [6,21]. Here, z and σ denote the local volume fraction of the product phase and the stress deviator, respectively, while φ (z) is some nonnegative saturation function for TRIP satisfying φ (0) = 0 and φ (1) = 1, see [6]. Besides lacking theoretical justification, this phenomenological approach suffers the downside of requiring a careful selection of the parameter κ either based on experiments, see e.g. [34] and [28], or empirical formulae with limited validity across materials and processing conditions. Furthermore, as Leblond et al. [15] show, plastic strains in the weaker phase driven by variations of temperature and applied stress may contribute to the material's inelastic response additionally. In a series of very well received publications at the end of the last century, Leblond et al. [14,16,17] thus proposed a constitutive model free of additional parameters to describe TRIP, which accounts for the underlying mechanisms by first adopting J 2thermoelastoplasticity on the microstructural scale and then deriving a macroscopic constitutive formulation by analytic homogenization. While refined models-see e.g. [3,7,25,29] among others-as well as generalizations of Leblond's reasoning to multiple product phases [18,39] have been published since then, only Leblond's model received noteworthy attention from researchers other than the respective original authors. This is most likely due to its comparatively concise form, the experimental scrutiny it received [34], the absence of difficult to obtain model parameters and an undisclosed implementation being available for use in a domain specific commercial finite element code designed for welding applications, see e.g. [43]. Nevertheless, apart from few applications of recent generalized formulations of the model by their respective originators, e.g. [19,41], Leblond's model is almost always not considered in full but is instead reduced to a special case of Eq. (1) to lessen the implementation effort [24,30,33]. Furthermore, as hinted above, TRIP is sometimes disregarded completely in heat treatment simulations and other applications of potential relevance (see e.g. [5,11,27] for such studies on quenching, welding and additive manufacturing) and instead, a modified yield stress depending on the local phase proportions is often used in conjunction with conventional J 2 -plasticity, as already suggested by Mitter [26] for ease of implementation. Hence, the purpose and structure of this work is two-fold: Firstly, following a concise review of the mechanisms of TRIP and Leblond's model in Sect. 2, we provide a fully implicit, straightforward and easy-to-implement scheme with consistent tangent modulus for Leblond's model in a notation suitable for immediate implementation in Sect. 3. To the best of our knowledge, the only algorithm published in detail for Leblond's original model is the one by Kim et al. [12] based on solving for the equivalent plastic strain update in consideration of isotropic hardening and thereby updating the stresses. However, since we consider ideal plastic phases and small strains, the averaged equivalent plastic strains of the phases play no role in the constitutive model and a considerably simpler and straightforward algorithm in terms of both derivation and implementation will be presented for this special case. This assessment also applies to the algorithm proposed by Lee et al. [18] for a generalization to multiple product phases. This algorithm may be reduced to our case of a single product phase. However, the reduction of their algorithm to ideal plastic phases is not obvious since they alternate between a Newton-Raphson iteration for the phases' equivalent plastic strain increments and a stress update step reliant on these increments. Additionally, computational efficiency may be impaired by multiple Newton-Raphson passes being required in their alternating scheme and especially the lack of an analytic consistent tangent modulus. The latter forces Lee et al. [18] to apply a finite-difference scheme that requires reevaluation of the stress response for every strain component perturbation.
Secondly, we apply the solution scheme to a tractable thermo-mechanical-metallurgical problem in Sect. 4 to discuss the differences to the previously mentioned simplified approach based on classical J 2 -plasticity and (a) Greenwood and Johnson's mechanism (b) Magee and Paxton's mechanism Fig. 1 Simplified representation of the microstructural mechanisms underlying transformation-induced plasticity in steels with austenite, ferrite/pearlite and martensite denoted by γ , α and α , respectively. Adapted from [9] a modified yield stress to highlight the potential impact of TRIP on residual stresses and strains and end with a conclusion in Sect. 5.

Mechanisms of transformation-induced plasticity in steels
According to Mitter [26], TRIP can be observed as a significantly increased plastic strain during phase transformations under applied loads for which the corresponding equivalent stress is small compared to the softer phase's yield strength. This macroscopic phenomenon is linked to the eigenstrain accompanying a phase transformation γ → α which in turn relates to changes of the crystal structure associated with the transformation. For example, the decomposition of closed-packed face centered cubic austenite into body centered cubic ferrite causes a local increase of the specific volume at the transformation site and the austenite to martensite decomposition is furthermore accompanied by a deviatoric eigenstrain. Two mechanisms driven by both the transformation eigenstrain and the applied stress are attributed to TRIP in steels [6,15], see Fig. 1: -Greenwood and Johnson's mechanism [10]: The dilatation at the transformation site must be accommodated by strains (illustrated red in Fig. 1a) in the surrounding matrix for reasons of compatibility, leading to internal stresses in the latter. These promote plasticity in the softer austenitic matrix if additionally an applied stress (resulting e.g. from external loading or inhomogeneous thermal strain of the component) is imposed. -Magee and Paxton's mechanism [23]: The applied stress may promote martensite to form with a preferred crystallographic orientation such that the associated deviatoric eigenstrains average to a nonzero strain on the component scale.
Since both mechanisms correspond to the microstructural scale, any constitutive model of transformation plasticity based on them must incorporate microstructural information of some kind into its formulation. As outlined in the following section, this holds in particular for Leblond's model which is based on Greenwood and Johnson's mechanism.

Constitutive modeling of transformation-induced plasticity in steels-Leblond's model
In engineering applications, the computational solution of problems involving metal plasticity commonly requires mathematical models and corresponding numerical solution schemes that are efficient enough to solve boundary value problems on geometrically complex components which are large in comparison with the microstructural features of the material. Hence, spatially resolving the microstructure in the numerical model often is computationally prohibitively expensive and the impact of the mechanisms described in the preceding section on the macroscopic inelastic response must be modeled by constitutive equations formulated exclusively on the component level continuum scale.
To this end, the plasticity model considered in this work and originally proposed by Leblond and co-workers [15][16][17] relies on an approximate analytic homogenization scheme to bridge the gap between the microstructural and the component scales, see Fig. 2. Postulating ideal small-strain thermoelastoplasticity governed by classical . 2 Schematic representation of the two-scale approach underlying Leblond's constitutive modeling for TRIP during a phase transformation γ → α and overview of the strain components on both the microstructural and the component scales. The strains on the component scale are given as volume averages of the strains on the microstructural scale over an imaginary representative volume V both sufficiently large to resolve local microstructural features and sufficiently small to neglect gradients of componentlevel quantities (e.g. the temperature gradient) in the derivation [15]. On the component scale, z = V α V −1 denotes the local volume fraction of the product phase α that occupies the subvolume V α ⊂ V J 2 -theory within the individual phases' domains, the total strain ε t on the microstructural scale is given by ε t = ε e + ε th + ε p . Here, ε e , ε th and ε p denote the elastic, the thermal and the plastic strain, respectively. Additionally, wherever transformation progresses at the phase boundary, the hydrostatic part Δε th γ →α I of the local eigenstrain jump Δε γ →α associated with the γ → α phase transformation is included in ε th while the corresponding deviator Δε p γ →α is included in ε p [15]. The hydrostatic part Δε th γ →α I corresponds to the difference of specific volume between the phases while the deviator Δε p γ →α comprises the shape-altering part of the transformation strain, e.g. shear strains related to martensitic transformation 1 . The corresponding strains e , thm and p on the component scale (see Fig. 2 for their definition) are related to the strains on the microstructural scale by means of volume averaging on the component scale. Based on the well-known constitutive equations of J 2 -thermoelastoplasticity applied at the microstructural scale and an analytic, approximate evaluation of the resulting averages, the governing equations for Leblond's model may moreover be formulated entirely on the component scale with the volume fraction z = V α V −1 ∈ [0, 1] of the product phase remaining as the only microstructural variable in the description. The analytic treatment is enabled by the assumption of plasticity being confined to the weaker γ -phase and the consideration of geometrically simplified surrogate models for the microstructure in parts of the derivation, among other more technical assumptions, see [16,17]. Thus, averaging of ε th -which is hydrostatic provided the thermal expansion of both phases is isotropicyields the "thermo-metallurgical" strain where th γ (T ) and th α (T ) denote the thermal strains obtained for the respective phases by means of a quenching dilatometry experiment [15]. Furthermore, as a result of the averaging, the inelastic part p of the total strain decomposes into three parts, where cp σ and cp T describe plastic strains associated to the variation of stress σ and temperature T , respectively 2 , and tp refers to transformation plasticity driven by the variation of the phase proportion z [15].
Since the total strain t is usually either by itself a primary variable of the solution scheme employed for the boundary value problem or directly linked to another primary variable (e.g. to the displacements in a displacement-based formulation) and the thermo-metallurgical strain thm is known a priori for given temperature T and phase proportion z, only p remains to be determined in typical engineering problems where the stress σ is to be computed based on the elastic strains e = t − thm − p . To this end, the following set of governing equations resulting from the aforementioned averaging of the common constitutive equations of J 2 -plasticity assumed for the γ -phase were derived by Leblond et al. [16,17]: 1 Leblond et al. [15] assume Δε p γ →α to be randomly orientated over the portion of the phase boundary in the representative volume V (seen in Fig. 2) where transformation progresses. Hence, Magee and Paxton's mechanism is disregarded and the deviatoric transformation eigenstrain Δε p γ →α cancels out during averaging so that only the hydrostatic part Δε th γ →α I remains in the resulting constitutive equations, Eqs. (5)-(9). 2 The superscript "cp" is used according to the original work of Leblond et al. [15][16][17].
If σ eq < σ Y ("TRIP-branch"): Else ("J 2 -branch"): Here, H (z) denotes the Heaviside function, z c the minimum volume fraction of the product phase 3 for the onset of the TRIP strain rate˙ tp , and g(z) the piece-wise linear function defined by interpolation of the values in Table 1. Total derivatives w.r.t. time are indicated by over-setting a dot. Furthermore, the hydrostatic part of the transformation strain Δ th γ →α and the function h(x) are given, respectively, by [17,20] and The von Mises equivalent stress σ eq as well as the yield strength σ Y , which is assumed as threshold for plastic flow of the whole phase mixture, are given as σ eq (σ ) = 3 2 σ : σ (12) and respectively [17], where we write σ for the deviator of the stress tensor. The function f (z) is again the piecewise linear interpolant of the values in Table 1. Finally, the material parameters involved are Young's modulus E and Poisson's ratio ν (assumed to be equal for both phases), the coefficients of thermal expansion α γ and α α and the yield strengths σ y γ and σ y α of the respective phases. Concluding the presentation of the model, we highlight that the subset of governing equations (5)-(9) in effect for the plastic strain rate depends on whether the current equivalent stress σ eq is below the mixed yield strength σ Y (z). The reasoning behind this is that only then the main assumption underlying the equations of the "TRIP-branch"-that plasticity shall be confined to the softer γ -phase-is admissible. Otherwise, the whole compound is assumed to be plastic and J 2 -theory is adopted throughout 4 . The general approach for the numerical integration of the constitutive equations is hence to first compute a plastic correction of an elastic predictor stress (obtained from Hooke's law by assuming zero plastic strain) using the equations of the "TRIPbranch" and subsequently, if the resulting stress exceeds the mixed yield strength σ Y , reverting the correction and instead correcting the elastic predictor using classical projection algorithms on the "J 2 -branch" [20,Ch. 9.14]. Note that Leblond's model therefore reduces to classical J 2 -plasticity if only one phase is present (i.e. before the phase transformation, z = 0, or after its completion, z = 1) since then˙ p =˙ cp σ =˙ cp T =˙ tp = 0 in the "TRIP-branch", implying no correction of the elastic predictor at all if its corresponding equivalent stress is smaller than σ Y . 3 According to Leblond et al. [17], this threshold is set to z c = 0.03, whereas Leblond later published the model in [20,Ch. 9.14] setting z c = 0 instead. The algorithm proposed in this paper is formulated to accommodate any choice of z c . 4 This appears reasonable within the context of Leblond's modeling, since the weaker γ -phase will be entirely plastic for stresses larger than σ Y ≥ σ y γ w.r.t. the von Mises norm, regardless of Greenwood and Johnson's mechanism. However, for martensitic transformations under very large applied stresses, Magee's mechanism may become the dominant mode of TRIP, as shown by Liu et al. [22], so that applying J 2 -plasticity for σ eq ≥ σ Y may be an oversimplification in such cases and constitutive models including Magee's mechanism may be necessary, see e.g. [7].

Implicit integration scheme for Leblond's model
In order to employ Leblond's model in solid mechanics, the mechanical boundary-value problem shall be solved by application of an incremental-iterative solution scheme since the constitutive model is both nonlinear and path dependent. Focusing on the evaluation of the constitutive equations, this usually involves the computation of updated mechanical state variables for a given deformation history [31]. As an example see [4,Ch. 2] for detailed application guidelines in the context of a displacement-based FE approach. Specifically, for the plasticity model at hand, the stress must be updated based on given values of phase proportion z, temperature T and total strain t . To this end, we derive an implicit integration scheme for Leblond's model (5)-(9) and furthermore the associated consistent tangent modulus which is usually required by the global incrementaliterative method applied to the mechanical boundary-value problem.

Notation and operator matrices
The derivations in the following sections will be formulated based on the following Voigt mapping for stresses and strains, σ T = σ 11 σ 22 σ 33 σ 12 σ 13 σ 23 T and T = 11 22 33 2 12 2 13 where we write single underscores for column vectors representing second-order tensors. Correspondingly, matrices representing fourth-order tensors will be written using double underscores henceforth. The presentation of the integration scheme is thus very convenient for implementation since the mapping Eq. (14), which is commonly used in finite-element frameworks, reduces algebraic tensor operations to basic operations on matrices and vectors with the aid of properly defined operator matrices. Throughout our work, the following operators are used: They are useful to represent double contractions of deviators of second-order tensors (e.g. σ : σ → σ T Pσ ) and the stress deviator (σ → Pσ ), respectively. Finally, we declare the representations of the second-and fourth-order identity tensors in Voigt-notation as respectively, where diag(•) denotes a square matrix with the operator's arguments as diagonal elements.

Discretization of the constitutive rate equations
We start the derivation of the solution scheme by considering the timestep t → t + Δt, where Δt > 0 denotes the increment of (pseudo-)time t. To obtain the incremental form of the constitutive equations presented in Sect. 2.2, we approximate the time derivatives by finite differences. For example, the plastic strain p may be expanded into a Taylor series at the time t to that end, and by dropping the terms of higher than first order, the time derivative˙ p can be approximated aṡ which amounts to the application of an implicit backward-Euler difference scheme if applied analogously to all rates in the constitutive equations at hand. It is obvious that the time increment Δt may be canceled out if such a difference quotient is substituted for each rate and by doing so, we arrive at based on the time derivatives of Eqs. (2), (4) and (3) and likewise at based on Eqs. (5), (7) and (9). Note that in these equations and in the remainder of this work, we adopt the convention that whenever no explicit indication is given as to whether a quantity (other than an increment) refers to the increment's beginning or its end, it shall always refer to its end. Furthermore, to simplify the notation we denote the time at the beginning and the end of the n-th increment as t n and t n+1 = t n + Δt, respectively, so that we may use the indices n and n + 1 to conveniently indicate quantities referring to these points of time, e.g. z n = z(t n ) and z n+1 = z(t n+1 ) = z.

Predictor-corrector scheme
As outlined in the conclusion of Sect. 2.2, the general approach to compute the stress σ n+1 at the end of a time increment t n → t n+1 corresponds to a predictor-corrector scheme, where we start by computing the elastic predictor using Hooke's law as presuming the increment to be entirely elastic and hence setting Δ p = 0 for the predictor step. Here, C denotes the isotropic elastic modulus given by where ν denotes Poisson's ratio. Then, a plastic correction for the predictor σ t is first computed using the "TRIP-branch" of the model to determine σ n+1 , see Sect. 3.4. This procedure may be interpreted as elasticplastic operator split as discussed by Simo and Hughes [31] for the case of J 2 -plasticity. Subsequently, if the equivalent stress σ eq (σ n+1 ) exceeds σ Y , the correction is revoked and the "J 2 -branch" is used to obtain a new correction for σ t , see Sect. 3.5.

Solution of the "TRIP-branch"
We start the derivation by solving Eq. (19) for the elastic strain increment Δ e and applying Hooke's law to obtain the stress increment, Adding σ n on both sides and using Eq. (24) immediately yields the expression for the corrector step, where we define the scalar function β 2 (σ n+1 ) for brevity of notation as Note that the nonlinear equation system Eq. (27) for the stress σ n+1 is similar in form to the system obtained for ideal J 2 -plasticity by application of a backward-Euler scheme (see e.g. [4, Ch. 7.3]), but β 2 (σ n+1 ) which is a function of the unknown stress appears instead of the plastic multiplier increment encountered in J 2 -theory.

If
Substituting Eq. (32) for the derivative of β 2 (σ n+1 ), reordering and inverting then eventually yields Note that the Jacobian occurring on this branch can be inverted very efficiently as follows: Let Then, the Jacobian J r σ (σ n+1 ) can be written as a rank-1 update of the matrix A, i.e.
According to the Sherman-Morrison formulae [2], the inverse of this matrix is then given as and due to the simple structure of A, the remaining inverse A −1 occurring in this expression is obtained as where G = 1 2 E(1 + ν) −1 denotes the shear modulus. Finally, we note that the results of this section may be equally applied to the modified Leblond model proposed by Taleb and Sidoroff [36] by means of the modifications detailed in Appendix A.

Solution of the "J 2 -branch"
If the stress σ n+1 obtained on the "TRIP-branch" exceeds σ Y (z) given by Eq. (13), we revert to the trial stress σ t and compute a new correction based on J 2 -theory (Eqs. (6) and (8)) to enforce the associated flow rule and admissibility of stress. To this end, the classical radial return algorithm of Wilkins [42] with the consistent tangent given by Simo and Taylor [32] may be applied. For completeness and convenience, we briefly state the procedure as summarized by Simo and Hughes [31] with adaptation to the problem at hand and the notation chosen in this paper. Let F t denote the yield function evaluated at the predictor and F t > 0 (otherwise elasticity applies). The unit normal to the von Mises yield surface is given by where n and n correspond to the Voigt spaces for stresses and strains, respectively. The plastic multiplier and the plastic strain increment may be computed in closed form as where G again denotes the shear modulus. Hence, the stress at the increment's end follows from Finally, for the ideal plastic case at hand, the consistent tangent modulus given in [31] degenerates to where K = 1 3 E(1 − 2ν) −1 is the bulk modulus and the scalar θ is defined by . (47)

Numerical study and comparison of Leblond's model to conventional J 2 -plasticity
To demonstrate the proposed solution scheme and to discuss the role of transformation-induced plasticity as described by Leblond's model in a tractable application, we consider the steel plate shown in Fig. 4a. The plate is assumed to be homogeneous, of initial thickness 2d and infinite along its symmetry plane perpendicular to the η-axis. Starting from a uniform initial temperature, the plate is symmetrically cooled down to ambient temperature by free convection on both surfaces, causing phase transformation from austenite to martensite. In this setting, the solution to the thermo-mechanical problem is symmetric about the mid-plane and spatially dependent on the η-coordinate only. Therefore, the analysis may be reduced to the half plate thickness and a sufficient spatial discretization of the plate is given by the finite element (FE) model shown in Fig. 4b. In this model, only one column of elements based on a generalized plane strain formulation is laid out along the η-direction. Furthermore, the heat conduction problem and the thermo-mechanical problem are sequentially coupled since we are concerned with small strains and disregard influences of the mechanical state on heat conduction and phase transformation.
For convenience, the parameters described in the following sections and their corresponding values for our study are collected in Table 2. We decided to assume the material parameters to be independent of temperature to simplify the discussion, although this is no requirement for the applicability of the proposed approach. Heat conduction: Initially, the plate is at uniform temperature T 0 . We assume linear isotropic heat conduction governed by Fourier's law. The outer edge of the FE model shown in Fig. 4b is subject to heat transfer defined by a constant heat transfer coefficient h T while the remaining edges are adiabatic so that the plate eventually cools down to ambient temperature T ∞ < T 0 . For simplicity, latent heat associated with the phase transformation is neglected and both the thermal conductivity λ and the specific heat capacity c p are assumed to be constant and equal for all phases. Kinetics of phase transformation: Besides the solution of the heat equation, Leblond's model must be complemented by evolution equations describing the phase transformation kinetics, i.e. the evolution of the volume fractions of all occurring phases. In this study, we enforce sufficiently fast cooling to deliberately restrict our   MPa yield strength of the α -phase z c 0.03 threshold phase proportion for˙ tp consideration to diffusionless phase transformations. In particular, we assume purely martensitic transformation γ → α starting with a completely austenitic microstructure and adopt the classical Koistinen-Marburger [13] relation to describe the martensite phase proportion z as function of temperature, Here, M s denotes the martensite start temperature and α M the rate parameter of the model.

Numerical solution of the heat equation and transformation kinetics
The computed temperature and phase proportion solutions are shown in Fig. 5a and b, respectively. Since the temperature near the surface decreases sooner and faster than in the core of the plate, the transient temperature field and hence the thermal contraction are inhomogeneous along the plate's thickness promoting transient thermal stresses. Phase transformation-as predicted by the Koistinen-Marburger relation-begins only once the martensite start temperature M s T 0 is reached at the surface while the transformation of the plate's interior is delayed. Note that the transformation remains incomplete upon reaching the stationary state due to insufficient undercooling experienced for the chosen parameter set.
The thermo-metallurgical eigenstrain thus comprises both the thermal contraction and-by virtue of the offset between the above functions-the expansion arising from the transformation, the latter of which clearly dominates once transformation starts, see Fig. 5c.
Concluding the paragraph, we note that Fig. 5 directly represents the three "driving forces" of TRIP for our considered problem which are associated to Eqs. (5), (7) and (9): The variations of temperature and phase proportion as well as the inhomogeneous thermo-metallurgical eigenstrain which by virtue of kinematic compatibility gives rise to variations of stress.

Numerical results for Leblond's model and conventional J 2 -plasticity
Based on the results discussed in Sect. 4.1, we now turn our attention to the mechanical problem and Leblond's model. Since we particularly intend to discuss the influence of TRIP, we compare the results obtained from Leblond's model with a "conventional" constitutive modeling approach to plasticity under phase transformations based on J 2 -theory that is commonly applied in engineering applications.
Let us clarify precisely what we refer to as "conventional" J 2 -approach. We thereby denote the application of ideal J 2 -plasticity with the yield strength defined by Eq. (13) and the thermo-metallurgical strain thm given by Eq. (3) and imposed as a (pseudo) thermal eigenstrain. This approach is usually trivially implemented even in commercial FE frameworks since both yield strength and thermal strain may typically be defined by straightforward user coding without the need to provide a solution scheme for the whole plastic constitutive law. Additionally, it is recovered as special case of Leblond's model when setting the plastic strains on the "TRIP-branch" to zero, cp σ ≡ cp T ≡ tp ≡ 0, consequently providing a natural counterpart neglecting TRIP. In the plate we consider, the stress state is equal biaxial so that σ ξξ is the only independent component of stress and σ ηη = 0 as well as σ i j = 0 for i = j. Since the plastic behavior is incompressible according to both Leblond's model and classical J 2 -theory 5 , the out-of-plane and in-plane plastic normal strains are related by p ηη = −2 p ξξ while no shear strains occur. Hence, we may limit the presentation to σ ξξ and p ξξ without loss of information. Furthermore, static equilibrium requires that the in-plane normal stress averages to zero along the thickness, thus establishing an intuitive reciprocal relation between regions loaded in tension and compression.
To start the discussion, we first consider the simpler case of conventional J 2 -plasticity. In Fig. 6a, the regions of the plate where plastification in tension or compression occurs along the in-plane direction are shown throughout the cooling of the plate. From this representation, an immediate qualitative assessment of the process can be drawn: -First, the surface of the plate is rapidly cooled while the plate's interior remains at the initial temperature, leading to thermal shrinkage of the outer material only, see Fig. 5c. Since the core counteracts this shrinkage, tensile stresses develop at the surface while the core is loaded in compression due to equilibrium, Eq. (50), as shown in Fig. 7a. Hence, as seen in Fig. 6a, the surface region becomes plastic in tension and subsequently the core plastifies compressively. The size of the plastic domain reaches a maximum at t ≈ 1 s when the core starts cooling and the thermal strain gradient is at its largest. -Next, the incipient thermal shrinking of the plate's core (ref. to Fig. 5c) relieves the stresses both in the surface and the core region so that the former and subsequently the latter gradually become purely elastic again while the region in between remains plastic in tension the longest. However, since the surface layers endured excessive plastic elongation in the in-plane direction, further thermal shrinking of the core gives rise to compressive stresses in the surface zone and in return to tension in the core, effectively reversing the stress distribution as seen in detail in Fig. 8a. Therefore, secondary plastification occurs at the surface, this time in compression. -At t ≈ 4 s the surface reaches the martensite start temperature and the γ → α phase transformation sets in, see Fig. 5b. From Fig. 5c it is seen that the corresponding increase of specific volume outweighs the thermal contraction locally once the transformation started. Consequently, the already existing compressive loading of the surface region intensifies while by virtue of equilibrium the interior necessarily reverts to tensile loading, Fig. 7a. At the same time, the formation of the much harder product phase significantly increases the yield strength of the phase mixture given by Eq. (13). Nevertheless, looking at Fig. 6a, compressive plastic flow in the surface region is sustained so that the local compressive stress continues to rise and equilibrium in turn requires the untransformed bulk material to become plastic in tension along an increasing portion of the plate's depth. -At t ≈ 7 s the phase transformation reaches the bulk of the plate (0 ≤ η 0.5d). The corresponding reversal of the thermo-metallurgical strain seen in Fig. 5c thus relieves the tensile stress in the bulk material, as observed in Figs. 7a and 8a, and by equilibrium, the compressive stress of the surface region is relaxed, too. Therefore, the surface zone becomes purely elastic again and subsequently, once the transformation has begun locally, the bulk material also returns to elasticity as its yield strength still increases due to the transformation. As a result of this, the plate gradually becomes elastic again throughout. -By then, the thermo-metallurgical strain gradually starts to saturate from the surface on into the bulk.
Hence, with the temperature gradient and the transformation subsiding, the outer regions of the plate essentially become strain driven by the plate's core material from then on. Therefore, the remaining thermometallurgical expansion of the core causes another reversal of stress along the thickness, as seen in Figs. 7a and 8a at t ≈ 10 s. Close inspection of Fig. 6a reveals a small region of tertiary plasticity in tension at the surface towards the transformation's end but apart from this, the plate remains thermoelastic until the steady state is reached.
Having understood the cooling process from the viewpoint of conventional J 2 -theory, we may now focus on the comparison to the results of Leblond's model to elaborate the predicted impact of TRIP. First, we  immediately find from Figs. 6, 7 and 8 that both constitutive models indeed yield the same results up to the onset of transformation at t ≈ 4 s, as expected. At this time, when the outer region of the plate is loaded in compression, the reversal of the local thermo-metallurgical expansion (Fig. 5c) as before promotes plasticity in compression, see Fig. 6b. But as Fig. 8b indicates, the plastic flow predicted by Leblond's model exceeds the prediction from J 2 -theory. Since up to this time the stresses, strains and governing equations were equal for both models, this difference arises from the solution of the "TRIP-branch". Indeed, in our example, we found this branch to exclusively govern the plastic behavior at any point in the plate subsequent to the local onset of transformation. Resulting from the just mentioned increased compressive plastic flow for Leblond's model, a slight initial relaxation of the surface region's stress is observed in Fig. 8a at the beginning of the transformation. Consequentially, in the initial stage, where transformation is confined to the near-surface zone of the plate, the stress continues to remain smaller there according to Leblond's model as opposed to J 2 -theory. Note that for reasons of equilibrium and compatibility, the predicted stresses and plastic strains in the bulk may already deviate ahead of the local transformation even though both plasticity models still coincide there at that time. In particular, equilibrium requires the stresses in the bulk of the plate to be smaller as well when compared to the prediction from J 2 -theory, as seen in Fig. 8a. Therefore, in the bulk of the plate, the domain of plasticity in tension prior to the local transformation appears significantly smaller and delayed, compare Fig. 6b to Fig. 6a. In a certain region between the plate's surface and its core, elasticity is even maintained up until the beginning of the local transformation according to Leblond's model, whereas J 2 -theory predicted plasticity in tension.
By the time the transformation reaches the inner half of the plate at t ≈ 8 s, the surface region-as noted above-is already elastic according to J 2 -theory due to the strongly increased yield stress σ Y . But since Leblond's model permits plasticity for stresses below this lower bound of J 2 -plasticity, the surface remains plastic in compression until the transformation almost reaches the plate's core, Fig. 8b. Accordingly, relaxation of the compressive stress in the middle part of the considered model (η ≈ 0.5d) is now possible due to TRIP. During this stress relaxation (seen just before t ≈ 9 s in Fig. 8a), the initial tensile stress in the middle region Fig. 9 Distributions of residual stress and plastic out-of-plane residual strain computed from Leblond's model (L.) and the conventional J 2 -plasticity approach. Note that the small "corrugation" in the distributions computed by J 2 -theory correspond to the small portion of the plate's thickness, where according to Fig. 6a no secondary plastification occured is at first slightly overcompensated by the plastic accommodation of the thermo-metallurgical strain increase, similar as previously observed at the surface upon the onset of transformation. However, as more of the interior transforms, an almost stress-free state is reached in the middle region until the core starts to transform.
As seen in Fig. 8a, the stresses finally reverse at t ≈ 9 s similar as predicted by J 2 -theory since the transformation near the surface subsides and the thermo-metallurgical strain of the bulk predominates. At that point, the transient stresses (as opposed to the plastic strains and hence the permanent deformation) are rather similar for both constitutive models. But since from then on J 2 -theory dictates elasticity virtually throughout the remainder of the cooling while Leblond's models promotes stress relaxation by plastic accommodation, the deviations subsequently become considerable especially near the plate's surface and core, see Fig. 8a. Driven by the core's thermo-metallurgical expansion, the stress evolution follows the same general pattern as before in J 2 -theory: The still expanding bulk of the plate becomes loaded in compression and in turn, the surface stress becomes tensile.
The plastic relaxation due to TRIP, however, results in much smaller magnitudes of stress so that the remaining residual stresses shown in Fig. 9a-although mostly coincident with respect to sign-differ considerably in magnitude. Hence, as seen from this example, neglecting TRIP potentially entails a large overestimation of transient and residual stresses since any plastic accommodation associated with plasticity in the weaker mother phase is disregarded for stresses smaller than an artificially defined mixed yield stress of similar type as given by Eq. (13). Moreover, the subsurface residual strains computed here even differ in sign along most of the plate's depth as seen in Fig. 9b where only the plastic component is shown since thm I is homogeneous in the steady state. This result indicates the possibility of predicting even qualitatively different permanent deformations for components of greater geometrical complexity than our plate when neglecting TRIP in computations of thermo-mechanical processes which involve phase transformations.
Lastly, we conclude the discussion by considering the individual contributions of the three strains cp σ , cp T and tp defined on the "TRIP-branch" by Eqs. (5), (7) and (9), respectively. In Fig. 10, the in-plane components of these strains are shown alongside the plastic strain J 2 ξξ accumulated on the "J 2 -branch" of Leblond's model during the same simulation. As noted, the plastic strains on the "TRIP-branch" vanish prior to the local beginning of the phase transformation whereas-in our example-plasticity is governed exclusively by the "TRIP-branch" afterwards. Furthermore, Fig. 10 demonstrates that plastic flow past the onset of transformation is dominated by tp , while the contribution of the other two components is rather small in this situation. Note, however, that these proportions are subject both to the material and the processing parameters. For example, in case of a slower phase transformation (i.e. smallerż) and a smaller difference of specific volume between the phases (i.e. smaller Δ th γ →α ) or higher rates of temperature and equivalent stress (largerṪ andσ eq , respectively 6 ), cp σ and cp T tend to have a larger impact on the overall plastic strain. As a final remark on the convergence characteristics of the proposed algorithm and the applicability of Leblond's model, we point the reader to Appendix B, where we demonstrate the convergence behavior in the FE context with the plate problem discussed in this section, as well as Appendix C, where predictions

Conclusion
In this work, we provide an implicit solution scheme for Leblond's model of TRIP suitable for straightforward implementation in computational frameworks for solid mechanics problems. Furthermore, the model is compared to a standard J 2 -plasticity-based approach in a tractable application by considering the symmetric cooling of an infinite plate. Despite the simplicity of this plate problem, the predictions for spatial and temporal evolution, intensity and even the very occurrence of plasticity in the presence of phase transformation differ significantly between the two constitutive models. Nonetheless, in this transparent example, the smaller transient and residual stresses could be attributed to the additional potential for plastic accommodation enabled by TRIP for stresses below an artificial yield strength of any sort as commonly assumed in conventional J 2 -plasticity approaches for the entire phase compound using simple mixture rules.
In applications of greater complexity regarding geometry, boundary conditions and thermo-mechanical loading, however, the interpretation and impact of TRIP is less apparent but the considered constitutive models may still predict substantially different distributions of residual stress and strain. Indeed, given the dependency of TRIP on the loading history in relation to the progression of the phase transformation, residual stress distributions differing even in sign from predictions of plain J 2 -theory across most of the component are within the realms of possibility, as we found in further numerical investigations of quenched cylindrical models under varying cooling rates. As an extreme case, though arguably a pathological one for practical choices of material and process parameters, the conventional J 2 -plasticity approach may even predict no residual stresses at all in opposition to Leblond's model if the transient stresses never reach the "mixed" yield strength of the compound. On the other hand, it may indeed be acceptable to assume elasticity for stresses below such an artificial compound yield strength if the phases' yield strengths are rather similar. In general however, with regard to its potentially unforeseeable implications, we conclude that TRIP should not be neglected a priori in thermo-mechanical investigations if phase transformations occur. Hence, particularly obvious applications of relevance include the simulation of heat treatment and welding processes.
In terms of implementation, the scheme we provide requires only a moderate amount of additional effort compared to implementing J 2 -plasticity. In particular, only a concise set of matrix-vector-operations that we covered in detail using the appropriate Voigt mapping remain as additional effort for implementing the "TRIPbranch". Hence, employing the algorithm as a user defined material model is straightforward even when having to comply with the often somewhat limited application programming interfaces of commercial FE software. However, from the viewpoint of computational efficiency, the integration point local Newton-Raphson iteration for solving the governing equations of TRIP obviously entails additional effort compared to J 2 -plasticity.
To conclude this work, we remark that further analysis should be carried out to enhance and apply the proposed numerical scheme and to further compare Leblond's model and plain J 2 -approaches. Regarding the constitutive modeling, the most notable limitation of the proposed scheme is the restriction to ideal plastic phases. Furthermore, the comparison of the constitutive models we considered should be extended to more complex geometries and processes amenable to cross-checking with carefully designed experiments comprising e.g. residual stress or plastic deformation measurements. Finally, an application to large scale problems would allow an evaluation of the additional computational cost compared to J 2 -plasticity (and hence of potential limits of practicality) when modeling and optimizing heat-treatment and welding of complexly shaped components. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.

A Adaption to Taleb and Sidoroff's modification of Leblond's model
Taleb and Sidoroff [36] proposed the following modification of Eq. (9) to replace the cutoff z c in the original formulation: G and K denote the shear and bulk modulus, respectively. Switching over to the Voigt notation used in this work, this may be expressed analogously to Eq. (23) as where H (z, z l ) is defined using the left-continuous Heaviside function H 0 (0) = 0, H 0 (z) = H (z) for z = 0, as H (z, z l ) = 1 − H 0 (z − z l ) ln(z l ) ln(z) The procedure and tangent described in Sect. 3.4 are then equally valid for the modified model, if H (z − z c ) is simply replaced by H (z, z l ) in β 2 (σ n+1 ) while β 1 (σ n+1 ) degenerates to . (54)

B Convergence behavior
Using the consistent tangent modulus plays an important role for the rate of convergence of the iterative method employed for the solution of the nonlinear mechanical boundary value problem, as pointed out by Simo and Taylor [32]. In particular, if the Newton-Raphson scheme is applied to solve the spatially discretized boundary value problem, as in our FE example, the quadratic rate of convergence of the Newton iteration may degrade and thus harm computational efficiency if no consistent tangent is provided for the constitutive response [31,32].
To provide insight into the convergence behavior of the algorithm proposed in this work within FE analyses, Fig. 11 shows the maximum residual forces in the model during the global Newton iteration for the quenching problem considered in Sect. 4.2. It can be seen that the expected superlinear convergence is indeed attained.

Fig. 11
Maximum residual forces ||F|| ∞ in the discretized model per global equilibrium iteration of the FE solver for all time increments that required more than two iterations towards equilibrium

C Comparison of predicted TRIP with experimental results from dilatometer tests
To compare the predictions of our implementation of Leblond's model with experimental data and thus to provide validation from a practical viewpoint, we consider the dilatometry experiments conducted by Taleb and Petit [35] in their investigation of TRIP in a low alloyed steel. Such an experiment may be adequately captured by a simple FE model consisting of a single element subject to a traction load and a predefined temperature evolution, since the uniaxial stress state and the temperature are homogeneous and known beforehand. Two thermal cycles between T ∞ = 20 • C and T max = 1100 • C have been prescribed following [35] and a power law has been assumed for the description of the phase transformations both during cooling-where bainite is formed according to Taleb and Petit [35]-and during austenitization: As per the dilatometry data given in [35], the transformation start and finish temperatures during cooling are T H = 560 • C and T L = 400 • C, respectively, while T L = 760 • C and T H = 840 • C for the (re)austenitization. Likewise, the thermal expansion coefficients of austenite, α γ , and of the product phase bainite, α α , as well as the strains at ambient temperature have been extracted from [35], while the elastic properties as well as the austenite's yield strength are taken from [40]: where T is given in • C. For the product phase, we assumed σ y α = 900 MPa. The first thermal cycle is stress free, whereas in the second thermal cycle a tensile load of 23 MPa is applied during cooling shortly before the beginning of the transformation. The strains computed using our proposed algorithm for Leblond's model agree well with the measurements reported by Taleb and Petit [35], as shown in Fig. 12. Note in particular that the strain evolution during the phase transformation in the final cooling step, which is the distinct feature governed by TRIP in this experiment, is well approximated by our simulation. In contrast, this feature cannot be predicted at all by a classical J 2based approach akin to the one discussed in Sect. 4.2 of this work. We remark that the experimental data of Taleb and Petit [35] displays a slightly different thermal expansion on cooling than on reheating which additionally appears to be somewhat nonlinear in temperature below T L during the final cooling, as opposed to our assumption of linear thermal expansion of each phase. However, it is difficult to draw a conclusion on the reasons for this based solely on the information given in [35]. Total axial strain during two subsequent thermal cycles (labeled C1 and C2, respectively) in a dilatometry experiment as predicted by the proposed implementation of Leblond's model compared to the experimental data provided by Taleb and Petit [35]