Phase-field Simulation of Habit Plane Formation during Martensitic Transformation in Low-carbon Steels

Abstract

The origin of the habit plane of the martensite phase (α′) in low-carbon steels is elucidated by three-dimensional phase-field simulations. The cubic → tetragonal martensitic transformation and the evolution of dislocations with Burgers vector a_α′/2〈111〉_α′ in the evolving α′ phase are modeled simultaneously. By assuming a static defect in the undercooled parent phase (γ), we simulate the heterogeneous nucleation in the martensitic transformation. The transformation progresses with the formation of the stress-accommodating cluster composed of the three tetragonal domains of the α′ phase. With the growth of the α′ phase, the habit plane of the martensitic cluster emerges near the (111)_γ plane, whereas it is not observed in the simulation in which the slip in the α′ phase is not considered. We observed that the formation of the (111)_γ habit plane, which is characteristic of the lath martensite that contains a high dislocation density, is attributable to the slip in the α′ phase during the martensitic transformation.

1. Introduction

In ferrous alloys, various types of martensite such as thin plate, lenticular, butterfly, and lath have been observed.¹⁾ Among them, the lath martensite appears in most heat-treatable commercial steels and is therefore extremely important from an industrial viewpoint. The characteristics of lath martensite, which is often observed in low-carbon steels, can be summarized as follows:^{2,3,4,5,6,7,8)} (1) the martensite phase (hereafter denoted as α′) contains a high dislocation density; (2) habit planes are close to the {111} or the {557} crystallographic planes of the parent phase (hereafter denoted as γ); and (3) the Kurdjumov–Sachs (K–S) orientation relationship ((111)_γ//(011)_α′,[101]_γ//[1 11]_α′) is satisfied between the γ phase (with a face-centered cubic (fcc) lattice) and the α′ phase (with a body-centered tetragonal (bct) lattice). The habit plane of the martensite coincides with the invariant plane because fitting the γ and α′ phases together along the invariant plane causes no elastic strain and minimizes the elastic energy associated with the transformation. The absence of elastic strain enables us to discuss the habit plane formation purely on the basis of crystallographic theory^{9,10,11,12,13,14,15,16)} without resorting to elasticity theory.

The phase-field method is a powerful tool for simulating the evolution of a material microstructure.¹⁷⁾ Phase-field models of proper and improper martensitc transformations (MTs) have been proposed by Khachaturyan and his collaborators.^{18,19,20,21,22)} The MT in steels have been simulated by simultaneously modeling the cubic → tetragonal lattice deformation and the plastic deformation in the γ and/or α′ phases.^{23,24,25,26,27,28,29,30,31,32,33)} In most of these studies, the elastic energy calculation is based on the microelasticity theory;^14,34) hence, the long-range elastic interaction between the transformation strain (Bain strain), plastic strain, and external applied stress is automatically taken into account and its effect on the microstructure evolution has been successfully simulated. However, the formation of the {111}_γ habit plane has not yet been simulated; although Yeddu et al.²⁶⁾ predicted the habit plane of an infinitesimal α′ phase as (111)_γ, the final habit plane was determined to be (211)_γ after the growth of the α′ phase. A phase-field model that can reproduce the lath-martensite characteristics, one of which is the {111}_γ habit plane, seems to be essential for simulating microstructure evolution in practical low-carbon steels. In phase-field simulations, the microstructure evolves so as to reduce the total free energy of the microstructure; hence, the formation of the {111}_γ habit plane should be simulated because it is presumably effective for minimizing the elastic energy during the MT.

As mentioned earlier, the lath martensite contains a high dislocation density.^2,6) As an inevitable result, the formation of the lath martensite is closely related to the slip in the α′ phase.^15,16) Presumably, the elastic energy minimization during the MT plays an important role in the habit plane formation in the lath martensite, where dislocations with Burgers vector a_α′/2〈111〉_α′ have often been observed.³⁾ There are some phase-field approaches in modeling defects and deformation in materials.³⁵⁾ In the microelasticity theory of dislocations,³⁶⁾ phase fields were defined to characterize inelastic strain field of individual dislocations on slip planes. Zhou et al.³⁷⁾ extended this model to a larger length scale and introduced a new set of phase fields (density functions for dislocation), which describes the plastic deformation produced by dislocations from a specific slip system. In this study, the cubic → tetragonal MT and the evolution of dislocations in the evolving α′ phase are simultaneously simulated by the phase-field method. Notably, for the first time, the elastic energy associated with both the Bain strain and the dislocations with Burgers vectors a_α′/2〈111〉_α′ in the α′ phase is taken into account; the density functions for dislocation³⁷⁾ are employed as field variables to describe the plastic deformation in the α′ phase. From simulation results with/without consideration of slip in the α′ phase, the origin of the {111}_γ habit plane of the lath martensite is explored.

2. Calculation Method

2.1. Phase-field Model

Long-range order parameters {φ_i(r, t)} = φ₁(r, t), φ₂(r, t), φ₃(r, t) (i = 1, 2, 3 numbering crystallographically equivalent orientation domains) and density functions for dislocation { $p_{\mathrm{\alpha }}^{(i)}$ (r, t)} = $p_1^{(i)}$ (r, t), $p_2^{(i)}$ (r, t), …, $p_n^{(i)}$ (r, t) (α = 1, 2, …, n numbering slip systems in the α′ phase) are adopted as field variables, which are functions of space (r) and time (t). φ_i(r, t) represents the structural density of each of the three tetragonal martensitic domains derived from the fcc → bct MT. Therefore, φ_i(r, t) is equal to unity if the vector r is inside a martensitic domain of i and is zero otherwise. $p_{\mathrm{\alpha }}^{(i)}$ (r, t) is related to $m_{\alpha }^{(i)}$ (r, t), which is the average number of lattice planes between neighboring sheared regions by dislocations of slip system α:   

$$p_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)=\frac{\left|{\mathbf{b}}_{\alpha }^{(i)}\right|}{D_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)}=\frac{\left|{\mathbf{b}}_{\alpha }^{(i)}\right|}{d_{\alpha }^{(i)}{m}_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)},$$(1)

where $d_{\alpha }^{(i)}$ is the lattice plane spacing and ${\mathbf{b}}_{\alpha }^{(i)}$ is the Burgers vector for slip system α. In Eq. (1), $D_{\alpha }^{(i)}$ (r, t)= $d_{\alpha }^{(i)}$ $m_{\alpha }^{(i)}$ (r, t) represents the average spacing between neighboring sheared regions by dislocations of slip system α.³⁷⁾ The temporal evolution of the field variables is given by the time-dependent Ginzburg–Landau (TDGL) equations:   

$$\frac{\partial {\phi }_i(\mathbf{r},\mathrm{\hspace{0.17em} }t)}{\partial t}=-L_{\phi }\frac{\delta G}{\delta {\phi }_i(\mathbf{r},\mathrm{\hspace{0.17em} }t)},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }i=1,\mathrm{\hspace{0.17em} }2,\mathrm{\hspace{0.17em} }3,$$(2)

  

$$\frac{\partial p_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)}{\partial t}=-L_p\frac{\delta G}{\delta p_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }i=1,\mathrm{\hspace{0.17em} }2,\mathrm{\hspace{0.17em} }3,\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\alpha =1,\mathrm{\hspace{0.17em} }2,\mathrm{\hspace{0.17em} }\mathrm{...},\mathrm{\hspace{0.17em} }n,$$(3)

where L_φ and L_p are kinetic coefficients and G is the total free energy of the microstructure. The G is a functional of field variables:   

$$\begin{array}{l}G={\int }_{\mathbf{r}}[f_{\text{L}}\left(\left\{{\phi }_i\right\}\right)+\frac{{\kappa }_{\phi }}{2}{\sum }_{i=1}^3{\left(\nabla {\phi }_i\right)}^2\\ \hspace{1em}+\frac{{\kappa }_p}{2}{\sum }_{i=1}^3{\sum }_{\alpha =1}^n{\left({\mathbf{n}}_{\alpha }^{(i)}\times \nabla p_{\alpha }^{(i)}\right)}^2]d\mathbf{r}+E_{\text{el}}\left(\left\{{\phi }_i\right\},\left\{p_{\alpha }^{(i)}\right\}\right),\end{array}$$(4)

where f_L is the Landau energy density, κ_φ and κ_p are gradient energy coefficients, ${\mathbf{n}}_{\alpha }^{(i)}$ is the unit vector of the slip plane normal,^36,37) and E_el is the elastic energy. Note that the “crystalline” energy³⁶⁾ is not included in Eq. (4), because $p_{\mathrm{\alpha }}^{(i)}$ (r, t) defined by Eq. (1) does not assume integer values.³⁷⁾ The f_L({φ_i}) is approximated by the Landau polynomial expansion with respect to {φ_i(r, t)}:^18,19,20)   

$$f_{\text{L}}=\mathrm{\Delta }f\left\{\frac{a}{2}{\sum }_{i=1}^3{\phi }_i^2-\frac{b}{3}{\sum }_{i=1}^3{\phi }_i^3+\frac{c}{4}{\left({\sum }_{i=1}^3{\phi }_i^2\right)}^2\right\}.$$(5)

With the conditions a − b + c = 0 and 6a − 4b + 3c + 12 = 0, Δf represents the driving force for the MT. In this study, a = 1, b = 15, and c = 14 are chosen.²⁰⁾ The E_el({φ_i},{ $p_{\mathrm{\alpha }}^{(i)}$ }) is given by^14,34)   

$$E_{\text{el}}={\int }_{\mathbf{r}}\left[\frac{1}{2}{C}_{ijkl}\left\{{\epsilon }_{ij}(\mathbf{r},\mathrm{\hspace{0.17em} }t)-{\epsilon }_{ij}^0(\mathbf{r},\mathrm{\hspace{0.17em} }t)\right\}\left\{{\epsilon }_{kl}(\mathbf{r},\mathrm{\hspace{0.17em} }t)-{\epsilon }_{kl}^0(\mathbf{r},\mathrm{\hspace{0.17em} }t)\right\}\right]d\mathbf{r},$$(6)

where C_ijkl is the elastic constant, ε_ij is the total strain, and ${\epsilon }_{ij}^0$ is the eigenstrain. Note that the Einstein summation convention is applied in Eq. (6). For simplicity, the material is assumed to be isotropic, and the tensor C_ijkl is expressed in terms of Lamé’s constants: C_ijkl = λδ_ijδ_kl + μ(δ_ikδ_jl + δ_ilδ_jk), where δ_ij is the Kronecker delta function. By solving mechanical equilibrium equation (∂σ_ij/∂r_j = 0, σ_ij = C_ijkl{ε_kl(r, t) − ${\epsilon }_{kl}^0$ (r, t)}) in Fourier space, we calculate the total strain ε_ij from the spatial distribution of the eigenstrain ${\epsilon }_{ij}^0$ . Here, the body is assumed to be fully clamped (strain-controlled boundary condition) and the macroscopic strain ε_ij≡∫_r ε_ijdr/V, where V is the system volume, is assumed to be zero. In this study, following the model of Zhang et al.,²²⁾ we model heterogeneous nucleation in the MT by assuming a static defect in the undercooled γ phase. The eigenstrain ${\epsilon }_{mn}^0$ is defined by   

$$\begin{array}{l}{\epsilon }_{mn}^0(\mathbf{r},\mathrm{\hspace{0.17em} }t)\\ ={\sum }_{i=1}^3\left\{{\epsilon }_{mn}^{(i)}{\phi }_i(\mathbf{r},\mathrm{\hspace{0.17em} }t)+{\sum }_{\alpha =1}^n\frac{{\mathbf{n}}_{\alpha }^{(i)}\otimes {\mathbf{b}}_{\alpha }^{(i)}+{\mathbf{b}}_{\alpha }^{(i)}\otimes {\mathbf{n}}_{\alpha }^{(i)}}{2\left|{\mathbf{b}}_{\alpha }^{(i)}\right|}{p}_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)\right\}\\ +{\epsilon }_{mn}^{(\text{d})}\zeta (\mathbf{r}).\end{array}$$(7)

Here, ${\epsilon }_{mn}^{(i)}$ is the transformation strain (Bain strain) derived from the fcc → bct MT, ${\epsilon }_{mn}^{(\text{d})}$ is the stress-free transformation strain induced by the static defect in the γ phase, and ζ(r) is a shape function that is equal either to unity if r is inside the static defect or to zero otherwise. In the undercooled γ phase, the MT is triggered by the static defect. As the α′ phase grows, the eigenstrain field of the static defect, which is described by ${\epsilon }_{mn}^{(\text{d})}$ ζ(r) in Eq. (7), is assumed to be inherited into the α′ phase. The ${\epsilon }_{mn}^{(i)}$ is given by   

$${\epsilon }_{mn}^{(1)}=\left(\begin{array}{ccc}{\epsilon }_a& 0& 0\\ 0& {\epsilon }_a& 0\\ 0& 0& {\epsilon }_c\end{array}\right),\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\epsilon }_{mn}^{(2)}=\left(\begin{array}{ccc}{\epsilon }_c& 0& 0\\ 0& {\epsilon }_a& 0\\ 0& 0& {\epsilon }_a\end{array}\right),\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{\epsilon }_{mn}^{(3)}=\left(\begin{array}{ccc}{\epsilon }_a& 0& 0\\ 0& {\epsilon }_c& 0\\ 0& 0& {\epsilon }_a\end{array}\right),$$(8)

where the coordinates are along the 〈100〉 axes of the γ phase and ε_a = $\sqrt{2}$ a_α′/a_γ−1 and ε_c = c_α′/a_γ−1 are the tetragonal transformation strains along the a and c axes, respectively. In the martensitic domains of φ₁ = 1, φ₂ = 1, and φ₃ = 1, the [001]_α′ coincides with the [001]_γ, [100]_γ, and [010]_γ, respectively. The Bain crystal lattice correspondence in each of the three tetragonal martensitic domains¹⁵⁾ is summarized in the Appendix.

2.2. Simulation Conditions

We used 64 × 64 × 64 computational grids for three-dimensional simulations. The periodic boundary conditions are assumed along all three dimensions. Equations (2) and (3) were solved by the finite difference method. In numerical analysis, all the physical parameters were reduced to dimensionless quantities using the scaling factors of RT for energy, l₀ for length, and (L_φRT)⁻¹ for time, where R is the gas constant, T is the absolute temperature, and l₀ is the unit grid size. The unit time step Δt*, where the superscript * denotes a dimensionless quantity, was selected as Δt* = 0.001 to maintain numerical accuracy and stability. The parameters used in the simulation are summarized in Table 1. The difference between the free energies of the austenite and ferrite phases in Fe–0.1mass%C at 300 K was assumed to be the driving force for the MT and was estimated from the Thermo-Calc TCFE7 database. Zhang et al.²²⁾ solved the PFM equation for two (110) twin-related martensitic domains and calculated the reduced domain wall tension coefficient. Through the use of their calculation result as a reference, the coefficient κ_φ was estimated by assuming the interphase boundary to be semicoherent (the interface energy γ^semicoh~0.24 J m^{−2 38)}); the estimated κ_φ gave a unit grid size l₀~8 nm (the system size was 522×522×522 nm³). The coefficient κ_p was determined so as to guarantee a smooth transition of $p_{\mathrm{\alpha }}^{(i)}$ (r, t) across the boundary between the sheared and unsheared regions.³⁷⁾ Lattice parameters a_γ, a_α′, and c_α′ in Fe–0.1mass%C at 300 K were estimated from the literature³⁹⁾ and from the empirical formula (c/a)_α′ = 1 + 0.045(mass%C).⁴⁰⁾ Elastic constants of pure iron were used.³⁹⁾ Because of the lack of experimental data, the kinetic coefficients L_φ and L_p were not well determined and were assumed to be equal to each other.

Table 1. Parameters used in the simulation.

Driving force for the MT (J m−3)	Δf = 7.25 × 108
Gradient energy coefficients (J m−1)	κφ = 1.90 × 10−8, κp = 1.78 × 10−7
Lattice parameters (10−10 m)	aγ = 3.60, aα′ = 2.87, cα′ = 2.88
Elastic constants (GPa)	λ = 123, μ = 72

The Landau energy density f_L for φ₂ = φ₃ = 0 is shown in Fig. 1. As shown in the inset in Fig. 1, a local maximum describes the chemical energy barrier of the transformation. With respect to the slip in the α′ phase, the time evolution of $p_{\mathrm{\alpha }}^{(i)}$ (r, t) was solved for four specific slip systems (α = 1,2,3,4) listed in Table 2, whereas, in practice, a number of slip systems would be activated during the MT. This assumption is based on the experimental result that dislocations with Burgers vector a_α′/2〈111〉_α′ occur in the lath martensite.³⁾ Note that the slip is confined to the α′ phase and that the time evolution of $p_{\mathrm{\alpha }}^{(i)}$ (r, t) is solved in the region where φ_i(r, t) ≥ 0.4; plastic yield criterion of the α′ phase is not considered. Figure 2 shows the initial microstructure of the simulation; the homogeneous γ phase contains a static defect (dislocation loop) on the (111)_γ plane at the center of the computational cell. The coordinate axes of the computational cell correspond to the [100]_γ, [010]_γ, and [001]_γ axes. The dislocation loop is assumed to have the Burgers vector b = M $\hat{\mathbf{b}}$ in the (111)_γ slip plane (n = 〈1/ $\sqrt{3}$ , 1/ $\sqrt{3}$ , 1/ $\sqrt{3}$ 〉), where $\hat{\mathbf{b}}$ = a_γ/2[110]_γ is an elementary Burgers vector, and M is a multiplicity factor characterizing the number of dislocation loops located within the same slip plane.²²⁾ Then, the stress-free transformation strain ${\epsilon }_{mn}^{(\text{d})}$ in Eq. (7) is given by ${\epsilon }_{mn}^{(\text{d})}$ = (b⊗n + n⊗b)/2l₀.³⁶⁾ The multiplicity factor is set to be M = 30, which is sufficiently large to describe the athermal MT triggered by the dislocation loop.²²⁾ The effect of the magnitude of M on the MT is beyond the scope of this paper and will be reported elsewhere.

Fig. 1.

Landau energy density (φ₂ = φ₃ = 0).

Table 2. Slip systems in the α′ phase assumed in the simulation.

α	Slip system
1	(101)α′[111]α′
2	(101)α′[111]α′
3	(101)α′[111]α′
4	(101)α′[1 11]α′

Fig. 2.

Initial microstructure for simulation. The γ phase is transparent. A static defect (dislocation loop) on the (111)_γ plane is located at the center of the computational cell.

3. Results

3.1. Formation of (111)_γ Habit Plane

First, simulation results obtained without considering slip in the α′ phase are shown in Fig. 3. Here, instead of Eq. (7), the eigenstrain ${\epsilon }_{mn}^0$ is defined by   

$${\epsilon }_{mn}^0(\mathbf{r},\mathrm{\hspace{0.17em} }t)={\sum }_{i=1}^3{\epsilon }_{mn}^{(i)}{\phi }_i(\mathbf{r},\mathrm{\hspace{0.17em} }t)+{\epsilon }_{mn}^{(\text{d})}\zeta (\mathbf{r}).$$(9)

In Fig. 3, the γ phase is transparent; three different colors represent three tetragonal domains (Bain variants) of the α′ phase that are twin-related with each other, and t* represents the dimensionless time step in the numerical analysis. The parent phase inside the dislocation loop first transforms to a martensitic domain (φ₃ ≥ 0.7 in red), which is subsequently surrounded by the other two martensitic domains (φ₁ ≥ 0.7 in green and φ₂ ≥ 0.7 in greenish yellow). This result is consistent with a previous simulation result reported by Zhang et al.²²⁾ that the martensitic nucleus prefers a multidomain, stress-accommodating cluster. As evident in Fig. 3, the MT progresses with the growth of the multidomain complex composed of the three tetragonal domains. The twin boundary between the two adjacent martensitic domains is on the {110}_γ plane. However, specific habit planes are not observed between the γ and α′ phases. In contrast, Fig. 4 shows the simulation results obtained with consideration of slip in the α′ phase. In this case, the eigenstrain ${\epsilon }_{ij}^0$ contains the term associated with $p_{\mathrm{\alpha }}^{(i)}$ (r, t), as given by Eq. (7). The growth of the martensitic cluster composed of the three tetragonal domains is observed; notably, as shown in Fig. 5, the habit plane of the multidomain cluster emerges near the (111)_γ plane (the blue-colored (111)_γ plane is inserted in Fig. 5 for reference). As evident in Figs. 3, 4, 5, the habit plane close to the (111)_γ plane is formed when the slip in the α′ phase is considered in the simulation. This result suggests that the slip in the α′ phase plays a prominent role in the formation of the (111)_γ habit plane, which is characteristic of the lath martensite with a high dislocation density.

Fig. 3.

Simulation results obtained without consideration of slip in the α′ phase. The γ phase is transparent, and the three different colors represent the three tetragonal domains of the α′ phase (φ₁ ≥ 0.7 in green, φ₂ ≥ 0.7 in greenish yellow, and φ₃ ≥ 0.7 in red). (Online version in color.)

Fig. 4.

Simulation results obtained with consideration of slip in the α′ phase. The γ phase is transparent, and the three different colors represent the three tetragonal domains of the α′ phase (φ₁ ≥ 0.7 in green, φ₂ ≥ 0.7 in greenish yellow, and φ₃ ≥ 0.7 in red). (Online version in color.)

Fig. 5.

Habit plane of the multidomain martensitic cluster shown in Fig. 4 at t* = 6. The blue-colored plane inserted in the left-hand figure is the (111)_γ plane. (Online version in color.)

3.2. Slip Distribution in the α′ Phase

In the simulation shown in Fig. 4, the evolution of $p_{\mathrm{\alpha }}^{(i)}$ (r, t) is solved for four slip systems (α = 1,2,3,4) listed in Table 2. Figure 6(a) shows the multidomain martensitic cluster at t* = 6, whereas Fig. 6(b) shows the domain of φ₁ ≥ 0.7 on a cross-section of the martensitic cluster. The spatial distributions of $p_1^{(1)}$ , $p_2^{(1)}$ , $p_3^{(1)}$ , and $p_4^{(1)}$ in the domain of φ₁ ≥ 0.7 on the corresponding cross-section are shown in Figs. 6(c), 6(d), 6(e), and 6(f), respectively. As evident in Fig. 6, all four of the slip systems are activated; however, the magnitude of $p_{\alpha }^{(1)}$ (α = 1,2,3,4) differs among the four slip systems. In the region marked by the circle in Fig. 6, the magnitudes of $p_3^{(1)}$ and $p_4^{(1)}$ (corresponding to the (101)_α′ [111]_α′ and (101)_α′ [1 11]_α′ slip systems, respectively) are larger than that of $p_1^{(1)}$ and $p_2^{(1)}$ (corresponding to the (101)_α′ [111]_α′ and (101)_α′ [111]_α′ slip systems, respectively). This result is reasonable because the α′ phase is under tensile stress along the c axis of the tetragonal lattice after the fcc → bct MT; given the Burgers vectors of the four slip systems, the combination of $p_3^{(1)}$ and $p_4^{(1)}$ (or $p_1^{(1)}$ and $p_2^{(1)}$ ) can obviously effectively reduce the elastic energy associated with the Bain strain. Notably, however, the magnitude of $p_3^{(1)}$ is slightly larger than that of $p_4^{(1)}$ in the region marked by the circle in Fig. 6. This result presumably suggests that the equivalent activation of the (101)_α′ [111]_α′ and (101)_α′ [1 11]_α′ slip systems is not necessarily effective for elastic energy minimization.

Fig. 6.

(a) Multidomain martensitic cluster shown in Fig. 4 at t* = 6, and (b) domain of φ₁ ≥ 0.7 on a blue-colored cross-section of the martensitic cluster. (c), (d), (e), and (f) show the spatial distribution of $p_1^{(1)}$ , $p_2^{(1)}$ , $p_3^{(1)}$ , and $p_4^{(1)}$ in the domain of φ₁ ≥ 0.7 on the cross-section, respectively. (Online version in color.)

There are twenty-four α′ variants for the K–S orientation relationship (OR). These variants are divided into four groups (CP groups); each of the group contains the six K–S variants whose close-packed planes ({011}_α′) are parallel to one of the four close-packed planes of the γ phase ({111}_γ). In each of the CP groups, the {111}_γ habit plane is shared by three Bain variants, and each of the Bain variants is shared by two K–S variants.^15,16) In this study, as we assumed the Bain crystal lattice correspondence shown in the Appendix, six K–S variants can derive from the three Bain variants, which share the (111)_γ habit plane, as a result of the slip in the α′ phase; the activation of the slip systems shown in Figs. 6(c)–6(f) would give rise to two K–S variants in the Bain variant domain of φ₁ ≥ 0.7, and in the same way, the other four K–S variants would emerge in the Bain variant domains of φ₂ ≥ 0.7 and φ₃ ≥ 0.7. As the phase-field model in this study is based on the microelasticity theory,¹⁴⁾ which is formulated in terms of strain that is invariant with respect to rigid body rotation, the OR between the γ and α′ phases cannot be predicted directly by the simulation. However, as discussed by Gao et al.,⁴¹⁾ there is one-to-one correspondence between the OR and invariant plane orientation. Hence, the OR can be determined if the invariant plane normal is predicted by the phase-field simulation based on the microelasticity theory.

3.3. Effect of the Initial Microstructure on the Habit Plane Formation

As shown in Figs. 4 and 5, the formation of the (111)_γ habit plane during the MT was simulated by the phase-field simulation. However, the configuration of the static defect (dislocation loop) shown in Fig. 2 could influence the simulation results; if the dislocation loop is located on a plane other than the (111)_γ plane, another simulation result would be obtained. Accordingly, a simulation was performed where the dislocation loop with the Burgers vector b = M $\hat{\mathbf{b}}$ (M = 30 and $\hat{\mathbf{b}}$ = a_γ/2[110]_γ) was located on the (111)_γ plane (n = 〈1/ $\sqrt{3}$ , 1/ $\sqrt{3}$ , −1/ $\sqrt{3}$ 〉) at the center of the computational cell, as shown in Fig. 7. The stress-free transformation strain ${\epsilon }_{mn}^{(\text{d})}$ in Eq. (7) is given by ${\epsilon }_{mn}^{(\text{d})}$ = (b⊗n + n⊗b)/2l₀. The simulated microstructure evolution is shown in Fig. 8, where the MT is observed to progress with the growth of the martensitic cluster composed of three tetragonal domains. The habit plane is not clearly observed until t* = 2; however, it emerges near the (111)_γ plane as the α′ phase grows. Figure 9 shows the simulated microstructure at t* = 6. A portion of the α′ phase grows out of the (111)_γ habit plane, as indicated by the arrows in Fig. 9; however, its effect on the habit plane formation appears to be negligible. Although this search was not exhaustive and more simulations with different initial microstructure are needed, we presumed that the formation of the (111)_γ habit plane is energetically favorable, irrespective of the configuration of the static dislocation loop.

Fig. 7.

Initial microstructure for simulation. The γ phase is transparent. A static defect (dislocation loop) on the (111)_γ plane is located at the center of the computational cell.

Fig. 8.

Results for the simulation where the initial microstructure shown in Fig. 7 was used. Slip in the α′ phase is considered. The γ phase is transparent, and the three different colors represent the three tetragonal domains of the α′ phase (φ₁ ≥ 0.7 in green, φ₂ ≥ 0.7 in greenish yellow, and φ₃ ≥ 0.7 in red). (Online version in color.)

Fig. 9.

Habit plane of the multidomain martensitic cluster shown in Fig. 8 at t* = 6. The blue-colored plane inserted is the (111)_γ plane. (Online version in color.)

4. Discussion

4.1. Slip Condition for an Invariant Plane to Exist

As previously mentioned, no elastic strain occurs along an invariant plane; hence, the invariant plane formation is effective for minimizing the elastic energy associated with the MT. The slip condition for an invariant plane to exist is considered on the basis of crystallographic theory, as shown below.¹⁴⁾ Let the atom positions of the γ phase (r_γ) shift to atom positions of the α′ phase ( ${r^{\prime }}_{{\alpha }^{\prime }}$ ) under the MT. For the martensitic domain of φ₁ = 1, the transformation can be written in the fcc crystal lattice basis as:   

$${{\mathbf{r}}^{\prime }}_{{\alpha }^{\prime }}=\mathbf{A}{\mathbf{r}}_{\gamma },$$(10)

where the transformation matrix A is given by   

$$\mathbf{A}=\left(\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ A_{31}& A_{32}& A_{33}\end{array}\right)$$

$A_{11}=\left(1-\frac{1}{2m_3^{(1)}}-\frac{1}{2m_4^{(1)}}\right){\eta }_1$ , $A_{12}=\left(\frac{1}{2m_3^{(1)}}+\frac{1}{2m_4^{(1)}}\right){\eta }_1$ , $A_{13}=\left(\frac{1}{2m_3^{(1)}}-\frac{1}{2m_4^{(1)}}\right){\eta }_1$ , $A_{21}=\left(\frac{1}{2m_1^{(1)}}+\frac{1}{2m_2^{(1)}}\right){\eta }_1$ , $A_{22}=\left(1-\frac{1}{2m_1^{(1)}}-\frac{1}{2m_2^{(1)}}\right){\eta }_1$ , $A_{23}=\left(-\frac{1}{2m_1^{(1)}}+\frac{1}{2m_2^{(1)}}\right){\eta }_1$ , $A_{31}=\left(-\frac{1}{2m_1^{(1)}}+\frac{1}{2m_2^{(1)}}-\frac{1}{2m_3^{(1)}}+\frac{1}{2m_4^{(1)}}\right){\eta }_3$ , $A_{32}=\left(\frac{1}{2m_1^{(1)}}-\frac{1}{2m_2^{(1)}}+\frac{1}{2m_3^{(1)}}-\frac{1}{2m_4^{(1)}}\right){\eta }_3$ , $A_{33}=\left(1+\frac{1}{2m_1^{(1)}}+\frac{1}{2m_2^{(1)}}+\frac{1}{2m_3^{(1)}}+\frac{1}{2m_4^{(1)}}\right){\eta }_3$ .Here, η₁ = $\sqrt{2}$ a_α′/a_γ and η₃ = c_α′/a_γ. $m_{\alpha }^{(i)}$ are the average number of lattice planes between neighboring regions sheared by dislocations of slip system α and are related to $p_{\mathrm{\alpha }}^{(i)}$ through Eq. (1). In general, the matrix A can be represented as the product of symmetric matrix F and rigid body rotation matrix R as   

$$\mathbf{A}=\mathbf{R}\mathbf{F}.$$(11)

If A⁺ is the transposed matrix of A, then   

$${\mathbf{A}}^{+}\mathbf{A}={(\mathbf{R}\mathbf{F})}^{+}(\mathbf{R}\mathbf{F})={\mathbf{F}}^{+}({\mathbf{R}}^{+}\mathbf{R})\mathbf{F}={\mathbf{F}}^2.$$(12)

If an invariant plane is present, one of the eigenvalues of F² will be less than unity, the other will be equal to unity, and the third will be greater than unity.¹⁴⁾ Figures 10(a) and 10(b) show the slip conditions for an invariant plane to exist when $m_1^{(1)}$ = 20 and $m_1^{(1)}$ = 80, respectively; an invariant plane exists if the slip condition is on the curvilinear surface shown in Figs. 10(a) and 10(b). We see that the geometry of the curvilinear surface has a small dependence on the magnitude of $m_1^{(1)}$ . The values of $m_2^{(1)}$ , $m_3^{(1)}$ , and $m_4^{(1)}$ in the simulation shown in Fig. 4 (at t* = 6) are extracted by Eq. (1) and are plotted in Figs. 10(c) and 10(d). Here, $m_1^{(1)}$ values are outside our scope of concern and are not shown in the figure. The simulated slip conditions shown in Figs. 10(c) and 10(d) distribute near the curvilinear surface shown in Figs. 10(a) and 10(b). This result suggests that the slip systems assumed in the simulation have been activated so as to form an invariant plane and to reduce the elastic energy during the MT. Interestingly, as a result of the slip in the α′ phase, the (111)_γ habit plane of a multidomain martensitic cluster was simulated (see Figs. 5 and 9) without any a prior assumption about the habit plane formation.

Fig. 10.

Slip conditions for an invariant plane to exist when (a) $m_1^{(1)}$ = 20 and (b) $m_1^{(1)}$ = 80. (c) and (d) show the values of $m_2^{(1)}$ , $m_3^{(1)}$ , and $m_4^{(1)}$ obtained from the simulation shown in Fig. 4 at t* = 6. (Online version in color.)

4.2. Plasticity during the MT

According to the phase-field microelasticity theory of dislocations,³⁶⁾ dislocations exist where the value of $p_{\mathrm{\alpha }}^{(i)}$ (r, t) varies. From the simulated spatial distribution of $p_{\mathrm{\alpha }}^{(i)}$ (r, t), the local dislocation density can be estimated by   

$$\rho (\mathbf{r},\mathrm{\hspace{0.17em} }t)={\sum }_{i=1}^3{\sum }_{\alpha =1}^4\frac{\left|{\mathbf{n}}_{\alpha }^{(i)}\times \nabla p_{\alpha }^{(i)}(\mathbf{r},\mathrm{\hspace{0.17em} }t)\right|}{\left|{\mathbf{b}}_{\alpha }^{(i)}\right|}.$$(13)

Figure 11 shows the change in the spatially-averaged dislocation density during the MT. The dislocation density increases as the MT progresses, and, when the transformation completes, it reaches a value of approximately ρ = 3.7 × 10¹⁶ m⁻² (solid line in Fig. 11); this value is an order of magnitude greater than the experimentally measured value.^2,6) In Eq. (13), by use of the central difference formula, the value of ∇ $p_{\mathrm{\alpha }}^{(i)}$ (r, t) is calculated from the values of $p_{\mathrm{\alpha }}^{(i)}$ (r, t) on the neighboring computational grids and the unit grid size l₀, which we estimated to be l₀~8 nm by assuming that the interphase boundary is semicoherent. However, if the interphase boundary is assumed to be incoherent (the interface energy γ^incoh = 0.5 – 1.0 J m^{−2 42)}), l₀ should be much larger and the value of ρ would therefore be smaller. As an example, if the interface energy is assumed to be γ^incoh = 1 J m⁻², then the l₀ and ρ are estimated to be l₀~34 nm and ρ~8.9 × 10¹⁵ m⁻², respectively (see the dashed line in Fig. 11). Furthermore, the values of kinetic coefficients L_φ and L_p in Eqs. (2) and (3) could influence the value of ρ. If the value of L_φ is significantly larger than that of L_p, the value of ρ will decrease. A detailed discussion on the effect of these modeling parameters on the value of ρ will be presented in a separate paper.

Fig. 11.

Change in the spatially-averaged dislocation density obtained from the simulation shown in Fig. 4; l₀ is the unit grid size.

On the basis of the simulation results shown in Figs. 3 and 4, formation of the (111)_γ habit plane can be attributed to the slip in the α′ phase. This result is consistent with the fact that the {111}_γ habit plane is characteristic of the lath martensite that contains a high dislocation density. However, the plastic deformation has also been reported to occur in the γ phase surrounding the lenticular and lath martensite in ferrous alloys.^43,44,45) Furthermore, it has been noted that the dislocation in the γ phase is inherited by the α′ phase near the interface region.⁴⁴⁾ Although the plastic strain in the γ phase would be smaller than that in the α′ phase,^23,24,26) the plastic deformation in the γ phase could reduce the elastic energy associated with the MT and hence could influence the transformation kinetics and martensitic morphology. In future works, the plasticity in the γ phase should also be incorporated into the model and its effect on the habit plane formation should be explored.

The MT is influenced by materials parameters such as kinetic coefficients L_φ and L_p, driving force for the MT, Bain strain, elastic constants, and plastic yield criteria of the γ and α′ phases. As the phase-field model in this study allows both twinning and slip modes, some parametric simulations, where each of the materials parameters is systematically changed, should be effective for elucidating the key materials parameters that determine whether to twin or to slip. In other words, in the near future, parametric simulation study would give an appropriate explanation of the experimental observations that twinning happens mainly in high-carbon steels at relatively low temperature while slip happens mainly in low-carbon steels at relatively high temperature.

5. Conclusions

The fcc → bct MT and the evolution of dislocations with Burgers vector a_α′/2〈111〉_α′ in the α′ phase were modeled simultaneously using the phase-field method. By assuming a static defect (dislocation loop) in the undercooled γ phase, we simulated heterogeneous nucleation in the MT. The transformation progressed with formation of the stress-accommodating cluster composed of the three Bain variant domains of the α′ phase. As the α′ phase grew, the habit plane of the multidomain martensitic cluster emerged near the (111)_γ plane, whereas it was not observed in the simulation where the slip in the α′ phase was not considered. The habit plane formation was also simulated in the case where the configuration of the static dislocation loop was changed, which presumably implies that the formation of the (111)_γ habit plane is energetically favorable, irrespective of the configuration of the static dislocation loop. We observed that the formation of a habit plane close to the (111)_γ plane, which is characteristic of the lath martensite with a high dislocation density, can be attributed to the slip in the α′ phase during the MT.

Acknowledgements

This work was supported by Grant-in-Aid for Young Scientists (B) (Grant Number 26820282) of the Japan Society for the Promotion of Science (JSPS), and was also supported by ISIJ Research Promotion Grant.

Appendix

The Bain crystal lattice correspondence in each of the three tetragonal martensitic domains is summarized.

“Martensitic domain of φ₁ = 1”   

$${\left[\begin{array}{c}a\\ b\\ c\end{array}\right]}_{\gamma }=\left(\begin{array}{ccc}1/2& 1/2& 0\\ -1/2& 1/2& 0\\ 0& 0& 1\end{array}\right){\left[\begin{array}{c}u\\ v\\ w\end{array}\right]}_{{\alpha }^{\prime }}$$

  

$${(hkl)}_{\gamma }={(pqr)}_{{\alpha }^{\prime }}\left(\begin{array}{ccc}1& -1& 0\\ 1& 1& 0\\ 0& 0& 1\end{array}\right)$$

“Martensitic domain of φ₂ = 1”   

$${\left[\begin{array}{c}a\\ b\\ c\end{array}\right]}_{\gamma }=\left(\begin{array}{ccc}0& 0& 1\\ 1/2& 1/2& 0\\ -1/2& 1/2& 0\end{array}\right){\left[\begin{array}{c}u\\ v\\ w\end{array}\right]}_{{\alpha }^{\prime }}$$

  

$${(hkl)}_{\gamma }={(pqr)}_{{\alpha }^{\prime }}\left(\begin{array}{ccc}0& 1& -1\\ 0& 1& 1\\ 1& 0& 0\end{array}\right)$$

“Martensitic domain of φ₃ = 1”   

$${\left[\begin{array}{c}a\\ b\\ c\end{array}\right]}_{\gamma }=\left(\begin{array}{ccc}-1/2& 1/2& 0\\ 0& 0& 1\\ 1/2& 1/2& 0\end{array}\right){\left[\begin{array}{c}u\\ v\\ w\end{array}\right]}_{{\alpha }^{\prime }}$$

  

$${(hkl)}_{\gamma }={(pqr)}_{{\alpha }^{\prime }}\left(\begin{array}{ccc}-1& 0& 1\\ 1& 0& 1\\ 0& 1& 0\end{array}\right)$$

References

• 1)   M.  Umemoto,  E.  Yoshitake and  I.  Tamura: J. Mater. Sci., 18 (1983), 2893.
• 2)   M.  Kehoe and  P. M.  Kelly: Scr. Metall., 4 (1970), 473.
• 3)   B. P. J.  Sandvik and  C. M.  Wayman: Metall. Trans. A, 14 (1983), 809.
• 4)   P. M.  Kelly,  A.  Jostsons and  R. G.  Blake: Acta Metall. Mater., 38 (1990), 1075.
• 5)   S.  Morito,  H.  Tanaka,  R.  Konishi,  T.  Furuhara and  T.  Maki: Acta Mater., 51 (2003), 1789.
• 6)   S.  Morito,  J.  Nishikawa and  T.  Maki: ISIJ Int., 43 (2003), 1475.
• 7)   H.  Kitahara,  R.  Ueji,  N.  Tsuji and  Y.  Minamino: Acta Mater., 54 (2006), 1279.
• 8)   S.  Morito,  X.  Huang,  T.  Furuhara,  T.  Maki and  N.  Hansen: Acta Mater., 54 (2006), 5323.
• 9)   J. S.  Bowles and  J. K.  Mackenzie: Acta Metall., 2 (1954), 129.
• 10)   J. K.  Mackenzie and  J. S.  Bowles: Acta Metall., 2 (1954), 138.
• 11)   J. S.  Bowles and  J. K.  Mackenzie: Acta Metall., 2 (1954), 224.
• 12)   B. P. J.  Sandvik and  C. M.  Wayman: Metall. Trans. A, 14 (1983), 835.
• 13)   P. M.  Kelly: Mater. Trans. JIM, 33 (1992), 235.
• 14)   A. G.  Khachaturyan: Theory of Structural Transformations in Solids, Dover, New York, (2008).
• 15)   K.  Iwashita,  Y.  Murata,  Y.  Tsukada and  T.  Koyama: Philos. Mag., 91 (2011), 4495.
• 16)   L.  Qi,  A. G.  Khachaturyan and  J. W.  Morris, Jr.: Acta Mater., 76 (2014), 23.
• 17)   L. Q.  Chen: Ann. Rev. Mater. Res., 32 (2002), 113.
• 18)   Y.  Wang and  A. G.  Khachaturyan: Acta Mater., 45 (1997), 759.
• 19)   A.  Artemev,  Y.  Jin and  A. G.  Khachaturyan: Acta Mater., 49 (2001), 1165.
• 20)   Y. M.  Jin,  A.  Artemev and  A. G.  Khachaturyan: Acta Mater., 49 (2001), 2309.
• 21)   Y.  Wang and  A. G.  Khachaturyan: Mater. Sci. Eng. A, 438–440 (2006), 55.
• 22)   W.  Zhang,  Y. M.  Jin and  A. G.  Khachaturyan: Acta Mater., 55 (2007), 565.
• 23)   V. I.  Levitas and  A. V.  Idesman: Philos. Mag. A, 82 (2002), 429.
• 24)   A.  Yamanaka,  T.  Takaki and  Y.  Tomita: Mater. Sci. Eng. A, 491 (2008), 378.
• 25)   H. K.  Yeddu,  A.  Malik,  J.  Ågren,  G.  Amberg and  A.  Borgenstam: Acta Mater., 60 (2012), 1538.
• 26)   H. K.  Yeddu,  A.  Borgenstam,  P.  Hedström and  J.  Ågren: Mater. Sci. Eng. A, 538 (2012), 173.
• 27)   A.  Malik,  H. K.  Yeddu,  G.  Amberg,  A.  Borgenstam and  J.  Ågren: Mater. Sci. Eng. A, 556 (2012), 221.
• 28)   H. K.  Yeddu,  V. I.  Razumovskiy,  A.  Borgenstam,  P. A.  Korzhavyi,  A. V.  Ruban and  J.  Ågren: Acta Mater., 60 (2012), 6508.
• 29)   Z.  Cong,  Y.  Murata,  Y.  Tsukada and  T.  Koyama: Mater. Trans., 53 (2012), 1598.
• 30)   Z.  Cong,  Y.  Murata,  Y.  Tsukada and  T.  Koyama: Mater. Trans., 53 (2012), 1822.
• 31)   Z.  Cong,  Y.  Murata,  Y.  Tsukada and  T.  Koyama: Philos. Mag., 93 (2013), 1739.
• 32)   H. K.  Yeddu,  A.  Borgenstam and  J.  Ågren: Acta Mater., 61 (2013), 2595.
• 33)   A.  Malik,  G.  Amberg,  A.  Borgenstam and  J.  Ågren: Acta Mater., 61 (2013), 7868.
• 34)   T.  Mura: Micromechanics of Defects in Solids, 2nd, rev. ed., Kluwer Academic, Dordrecht, (1987).
• 35)   Y.  Wang and  J.  Li: Acta Mater., 58 (2010), 1212.
• 36)   Y. U.  Wang,  Y. M.  Jin,  A. M.  Cuitiño and  A. G.  Khachaturyan: Acta Mater., 49 (2001), 1847.
• 37)   N.  Zhou,  C.  Shen,  M.  Mills and  Y.  Wang: Philos. Mag., 90 (2010), 405.
• 38)   Z.  Yang and  R. A.  Johnson: Model. Simul. Mater. Sci. Eng., 1 (1993), 707.
• 39)  Kinzoku Data Book, 4th rev. ed., ed. by The Japan Institute of Metals, Maruzen, Tokyo, (2004), (in Japanese).
• 40)   C. S.  Roberts: Trans. AIME, 197 (1953), 203.
• 41)   Y.  Gao,  N.  Zhou,  D.  Wang and  Y.  Wang: Acta Mater., 68 (2014), 93.
• 42)   D. A.  Porter,  K. E.  Easterling and  M. Y.  Sherif: Phase Transformations in Metals and Alloys, 3rd ed., CRC Press, Boca Raton, (2009).
• 43)   A.  Shibata,  S.  Morito,  T.  Furuhara and  T.  Maki: Scr. Mater., 53 (2005), 597.
• 44)   A.  Shibata,  S.  Morito,  T.  Furuhara and  T.  Maki: Acta Mater., 57 (2009), 483.
• 45)   G.  Miyamoto,  A.  Shibata,  T.  Maki and  T.  Furuhara: Acta Mater., 57 (2009), 1120.