cta mater. VoL. 45, No. 2, pp. 759-773. 1997 Pergamon Copyright c 1997 Acta Metallurgica Published by elsevier Science Printed in Great Britain. All rights reserved 1359-645497517.00+0.00 THREE-DIMENSIONAL FIELD MODEL AND COMPUTER MODELING OF MARTENSITIC TRANSFORMATIONS Y WANGt and A G KHACHATURYAN Department of Ceramics, Rutgers University, P. O. Box 909, Piscataway. New Jersey 08855-0909, U.S.A. Receied 14 December 1995; accepted 8 April 1996) Abstract-A three-dimensional (3D) continuum stochastic field kinetic model of martensitic which explicitly takes into account the transformation- induced elastic strain is developed he model is able to predict the major structural characteristics of martensite during the entire ransformation including nucleation, growth and eventually formation of internally twinned plates which are in thermoelastic equilibrium with the parent phase. No a priori constraints are made on the possible configurations and sequences of structural patterns formed by orientation variants of the martensite. 3D computer simulations are performed for a generic cubic- tetragonal martensitic transformation in a prototype crystal which is elastically isotropic and elastically homogeneous. The simulations predict that (i)nucleation of martensite in a perfect crystal occurs collectively to accommodate the coherency strain. e.g. the critical nuclei are formed by two ally twinned orientation variants: (ii) the ultimate structure consists of plate-like martensite and retaining parent phase. The martensitic plates consist of twin-related platelets of two orientation variants and the habits of the plates meet the invariant plane requirement These simulation results arc in good agrecment with expcrimental observations. Copyright@ 1997 Acta Metallurgica 1 INTRODUCTION force, IAf(TI,(a difference between the specific Martensitic transformation is still one of the least stress-frce free energies of the martensitic and parent studied solid state reactions. The transformation has specific features that distinguish itself from the conventional phase transformations described by classical Gibbs thermodynamics. These features where G is the typical shear modulus, fo is the typical include (i) the occurrence of a thermoelastic stress-free transformation strain and T is the two-phase equilibrium within a certain temperature temperature. If s>1, then the strain energy range;(i)the stress-induced transformation and contribution to the total free energy during the related hysteresis; (ii) the shape memory effect and transformation exceeds or is at least commensurate (iv) the appearance of a typical"martensitic" with the transformation driving force and hence the microstructure consisting of internally twinned martensitic features of the transformation prevai te unusual in the Due to the close relationship between the could not be explained without revisions of some microstructure and macroscopic behavior of ma basic concepts of the conventional theory of phase terials, intensive research has been made during the transformation last few decades towards understanding the origin of It is generally accepted that the above distinctions the unique structures of martensite. The crystallo- of martensitic transformation are caused by the graphic(geometrical)theories of martensitic trans- transformation-induced elastic strain involved in the formation proposed by Wechsler, Lieberman and ansformation, an effect ignored in the conventional Read(WLr)[2] and Bowles and MacKenzie(BM)[31 stress-free themodynamics. This strain is generated explain the habit plane orientation and multidomain by the crystal lattice misfit between the martensite structure of martensite reasonably well.These and its parent phase as well as between different theories are based on an assumption that the orientation variants of the martensite. The criterion austenite martensite interfaces are invariant under indicating in what situations the unique"martensitic" the transformation and thus macroscopically features are expected is given by the ratio of the perfect fit along the habit plane interfaces is obtained typical stra Gea. to the*chemical " dri The invariance is realized through the formation of a multidomain structure of martensite consisting of two twin-related orientation variants. The Wlr and +Permanent address tment of Materials Science BM theories are formulated in geometrical neering, The terms since under the invariant plane interface 43210,USA hypothesis the system is stress free and, thus. no
PI1 s1359-6454(%)00180-2 Acta mum. Vol. 45, No. 2, pp. 159-713. 1997 Copyright 0 1997 Acta Metallurgica Inc. Published by Elsevier Science Ltd Printed in Great Britain. All rights reserved 13596454197 $17.00 + 0.00 THREE-DIMENSIONAL FIELD MODEL AND COMPUTER MODELING OF MARTENSITIC TRANSFORMATIONS Y. WANG? and A. G. KHACHATURYAN Department of Ceramics, Rutgers University, P. 0. Box 909. Piscataway. New Jersey 08855-0909, U.S.A. (Received 14 December 1995; accepted 8 April 1996) Abstract-A three-dimensional (3D) continuum stochastic field kinetic mode1 of martensitic transformations which explicitly takes into account the transformation-induced elastic strain is developed. The model is able to predict the major structural characteristics of martensite during the entire transformation including nucleation, growth and eventually formation of internally twinned plates which are in thermoelastic equilibrium with the parent phase. No a priori constraints are made on the possible configurations and sequences of structural patterns formed by orientation variants of the martensite. 3D computer simulations are performed for a generic cubic-+tetragonal martensitic transformation in a prototype crystal which is elastically isotropic and elastically homogeneous. The simulations predict that (i) nucleation of martensite in a perfect crystal occurs collectively to accommodate the coherency strain, e.g. the critical nuclei are formed by two internally twinned orientation variants; (ii) the ultimate structure consists of plate-like martensite and retaining parent phase. The martensitic plates consist of twin-related platelets of two orientation variants and the habits of the plates meet the invariant plane requirement. These simulation results are in good agreement with experimental observations. Cop>vight 0 1997 Acta Metallurgica Inc. 1. INTRODUCTION Martensitic transformation is still one of the least studied solid state reactions. The transformation has specific features that distinguish itself from the conventional phase transformations described by classical Gibbs thermodynamics. These features include (i) the occurrence of a thermoelastic two-phase equilibrium within a certain temperature range; (ii) the stress-induced transformation and related hysteresis; (iii) the shape memory effect and (iv) the appearance of a typical “martensitic” microstructure consisting of internally twinned plates. They are quite unusual in the sense that they could not be explained without revisions of some basic concepts of the conventional theory of phase transformation. It is generally accepted that the above distinctions of martensitic transformation are caused by the transformation-induced elastic strain involved in the transformation, an effect ignored in the conventional stress-free themodynamics. This strain is generated by the crystal lattice misfit between the martensite and its parent phase as well as between different orientation variants of the martensite. The criterion indicating in what situations the unique “martensitic” features are expected is given by the ratio of the typical strain energy, GE,‘, to the “chemical” driving tpermanent address: Department of Materials Science and Engineering, The Ohio State University, Columbus, OH 43210, U.S.A. force, IAf(T) 1, (a difference between the specific stress-free free energies of the martensitic and parent phases)[l]: i = G&V(T) I (1) where G is the typical shear modulus, to is the typical stress-free transformation strain and T is the temperature. If [ > 1, then the strain energy contribution to the total free energy during the transformation exceeds or is at least commensurate with the transformation driving force and hence the martensitic features of the transformation prevail. Due to the close relationship between the microstructure and macroscopic behavior of materials, intensive research has been made during the last few decades towards understanding the origin of the unique structures of martensite. The crystallographic (geometrical) theories of martensitic transformation proposed by Wechsler, Lieberman and Read (WLR) [2] and Bowles and MacKenzie (BM) [3] explain the habit plane orientation and multidomain structure of martensite reasonably well. These theories are based on an assumption that the austenitejmartensite interfaces are invariant under the transformation and thus a macroscopically perfect fit along the habit plane interfaces is obtained. The invariance is realized through the formation of a multidomain structure of martensite consisting of two twin-related orientation variants. The WLR and BM theories are formulated in pure geometrical terms since under the invariant plane interface hypothesis the system is stress free and, thus, no 759
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL elastic strain analysis is required. Further develop- limited number of structural configurations for the ment of these theories where the problem of critical nucleus. They cannot provide a complete test microstructure optimization is treated under a zero of all possible heterogeneous structural configur stress condition has been proposed by ball and James ations which could be the critical nuclei. To do this 4, Kohn [5]and Bhattacharaya [6] a stochastic dynamic approach has to be developed The micromechanic approach to the problem was which should provide a random sampling of all developed by Khachaturyan and Shatalov [7]. They configurational space and eventually single out the formulated the strain energy functional of an critical nucleus heterogeneity arbitrary coherent two-phase microstructural pattern Summing up the foregoing discussion, it is quite for the homogencous modulus casc. Minimization of obvious that further progress in thcorctical character the strain energy functional for arbitrary ization of martensitic transformations requires the microstructure results in the observed morphological development of a dynamic approach to the problem pattern of martensite, i.e. a polydomain plate However, any dynamic approach has to deal with consisting of alternating twin-related lamellae of complex coherent multiparticle configurations with different orientation variants. The habit of this plate the particles interacting with each other through the is an invariant plane. Roytburd obtained the same infinitely long-ranged overlapping strain fields.It sult by minimizing the strain energy of polydomain seems that using numerical methods is the only way plate with the invariant plane habit with respect to an to advance in this area. With rapid advances in alternating period [8]. The theory of branching modern computer technology, several numerica domain walls on the martensite/austenite interface, approaches have been made. Wen et al. [11] and imilar to the lifshits theory of branching magnetic recently Xu and morris[12]employed a finite element domains, has been developed by Kohn and Muller model based on the strain energy functional proposed in Ref [7]. Their computer simulations of two-dimen- Although the strain energy minimization predicts sional (2D) martensitic transformations have pre- some important aspects of the martensitic micro- dicted many characteristics typical of martensitic structurc,it is insufficient to providc a complctc transformations and have given important insights on description of the morphology and the transform- the origin of the thermal hysteresis of the transform ation thermodynamics and kinetics. For example a ations. However, there is one critical aspect of the single-phase and single orientation variant state of transformation that a 2D model is unable to martensitic phase providing zero strain energy is reproduce. This is the formation of multidomain never observed. To fully understand the transform- martensite plates with invariant plate habits. The ation, one needs to consider the transformation reason for this shortcoming is the principle dynamics. The reason for this is as follows. If qualiLative difTerence between 2D and 3D sysLems condition(1) holds, the martensitic transformation The formation of a 3D single-variant martensitic can proceed only by developing a series of particle would always raise the volume-dependent multivariant structure patterns which provide the strain energy which is prohibitively high for the elastic strain accommodation. Escaping from this particle to be formed. To avoid this strain energy path (deviation from the strain-accommodating increase, the martensitic phase is formed as plates of microstructures)would be penalized by a macro- invariant plane habits, each plate consisting of alternating lame result, the martensitic system usually cannot reach variants. This morphology is the main feature of the thermodynamic equilibrium, i.e. a single-variant state observed martensitic structures. Unlike the 3D which has the lowest free energy. Instead, it is model, a 2D model permits formation of a arrested at the first encountered deep free energy local single-variant particle since the volume-dependent minimum along the strain-accommodating trans- strain energy contribution can always be eliminated formation path. Therefore, the martensitic systems if this particle is a"strip"along an invariant line meeting condition (1)are not ergodic and the which is not aftected by the crystal lattice ultimate microstructure of the martensitic phase is rearrangement. Such a line practically always exists in path dependent. 2D modcl of martensitic transformations, The Another important aspect for the fundamental opportunity to form a single-domain particle in 2D concerns the nucleation mechanisms of martensite. a multidomain martensitic plates. Therefore, any good and detailed review of this problem can be realistic computer modeling of martensitic transform- found in Ref. [10]. A well-accepted argument ations which is able to catch the main physics of the concerning the homogeneous nucleation is that the morphological pattern formation should be formu- strain cncrgy contribution to the nucleation barrier of lated in 3D. An important step to 3D numerical single variant martensitic particle is too high to be modeling of martensitic transformations has been overcome by thermal fluctuations. However, it made by Clapp et al. [13] who used the Molecular should be noted that the existing nucleation models, Dynamic(MD) method. The massive MD approach both classical and non-classical, assume only a is straightforward and is certainly the future for the
760 WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL elastic strain analysis is required. Further development of these theories where the problem of microstructure optimization is treated under a zero stress condition has been proposed by Ball and James [4], Kohn [5] and Bhattacharaya [6]. The micromechanic approach to the problem was developed by Khachaturyan and Shatalov [7]. They formulated the strain energy functional of an arbitrary coherent two-phase microstructural pattern for the homogeneous modulus case. Minimization of the strain energy functional for an arbitrary microstructure results in the observed morphological pattern of martensite, i.e. a polydomain plate consisting of alternating twin-related lamellae of different orientation variants. The habit of this plate is an invariant plane. Roytburd obtained the same result by minimizing the strain energy of polydomain plate with the invariant plane habit with respect to an alternating period [8]. The theory of branching domain walls on the martensite/austenite interface, similar to the Lifshits theory of branching magnetic domains, has been developed by Kohn and Miiller 191. Although the strain energy minimization predicts some important aspects of the martensitic microstructure, it is insufficient to provide a complete description of the morphology and the transformation thermodynamics and kinetics. For example, a single-phase and single orientation variant state of martensitic phase providing zero strain energy is never observed. To fully understand the transformation, one needs to consider the transformation dynamics. The reason for this is as follows. If condition (1) holds, the martensitic transformation can proceed only by developing a series of multivariant structure patterns which provide the elastic strain accommodation. Escaping from this path (deviation from the strain-accommodating microstructures) would be penalized by a macroscopically large increase in the strain energy. As a result, the martensitic system usually cannot reach thermodynamic equilibrium, i.e. a single-variant state which has the lowest free energy. Instead, it is arrested at the first encountered deep free energy local minimum along the strain-accommodating transformation path. Therefore, the martensitic systems meeting condition (1) are not ergodic and the ultimate microstructure of the martensitic phase is path dependent. Another important aspect for the fundamental understanding of the microstructure development concerns the nucleation mechanisms of martensite. A good and detailed review of this problem can be found in Ref. [lo]. A well-accepted argument concerning the homogeneous nucleation is that the strain energy contribution to the nucleation barrier of a single variant martensitic particle is too high to be overcome by thermal fluctuations. However, it should be noted that the existing nucleation models, both classical and non-classical, assume only a limited number of structural configurations for the critical nucleus. They cannot provide a complete test of all possible heterogeneous structural configurations which could be the critical nuclei. To do this, a stochastic dynamic approach has to be developed, which should provide a random sampling of all configurational space and eventually single out the critical nucleus heterogeneity. Summing up the foregoing discussion, it is quite obvious that further progress in theoretical characterization of martensitic transformations requires the development of a dynamic approach to the problem. However, any dynamic approach has to deal with complex coherent multiparticle configurations with the particles interacting with each other through the infinitely long-ranged overlapping strain fields. It seems that using numerical methods is the only way to advance in this area. With rapid advances in modern computer technology, several numerical approaches have been made. Wen et al. [ll] and recently Xu and Morris [12] employed a finite element model based on the strain energy functional proposed in Ref. [7]. Their computer simulations of two-dimensional (2D) martensitic transformations have predicted many characteristics typical of martensitic transformations and have given important insights on the origin of the thermal hysteresis of the transformations. However, there is one critical aspect of the transformation that a 2D model is unable to reproduce. This is the formation of multidomain martensite plates with invariant plate habits. The reason for this shortcoming is the principle qualitative difference between 2D and 3D systems. The formation of a 3D single-variant martensitic particle would always raise the volume-dependent strain energy which is prohibitively high for the particle to be formed. To avoid this strain energy increase, the martensitic phase is formed as plates of invariant plane habits, each plate consisting of alternating lamellae of two twin-related orientation variants. This morphology is the main feature of the observed martensitic structures. Unlike the 3D model, a 2D model permits formation of a single-variant particle since the volume-dependent strain energy contribution can always be eliminated if this particle is a “strip” along an invariant line which is not affected by the crystal lattice rearrangement. Such a line practically always exists in a 2D model of martensitic transformations. The opportunity to form a single-domain particle in 2D removes the main driving force for the formation of multidomain martensitic plates. Therefore, any realistic computer modeling of martensitic transformations which is able to catch the main physics of the morphological pattern formation should be formulated in 3D. An important step to 3D numerical modeling of martensitic transformations has been made by Clapp et al. [13] who used the Molecular Dynamic (MD) method. The massive MD approach is straightforward and is certainly the future for the
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL omputer simulation of phase transformations, but 2 FIELD KINETIC MODEL OF MARTENSITIC right now even the power of state-of-the-art parallel TRANSFORMATIONS computers is insufficient to carry out MD simulations for systems of a reasonable size (i.e. big enough to The continuum field approach to martensitic exclude or, at least to make qualitatively insignificant. ransformations is based on a simple idea that an e size effect) during a time of more than about arbitrary microstructure. consisting of all possible orientation variants of martensite coherently imbed An alternative approach is the Monte Carlo(MC) ded in its parent phase matrix, can be described imulation. It is an efficient technique for studying bv tinuous functions (field variables) lacie transformations involving simple atomic n (r)where r is the coordinate vector xchanges among rigid lattice sites. Recently, 3D Each of them describes an Iro parameter profile MC simulation has been used to model ordering the orientation variants. Thus the total number of important role [14. 15]. In this approach. the orientation variants which are determined by the ordering is described by the Ising model and the symmetry of the crystal lattice rearrangement. Each simulation has demonstrated that the transform- function assumes a non-zero value within martensitic ation strain accommodation caused by the ordering particles of the corresponding orientation variant and produces the so-called tweed mesoscopic structure. vanishes within the parent phase or other orientation As far as we know, however, this approach has not variants. Using this approach. Wang et al.have been applied to the study of martensitic transform- simulated the microstructure evolution during pre- ations cipitation of tetragonal ZrO, particles from a cubi At the present stage of computer technology, the matrix in partially stabilized zirconia [17 ontinuum field method based on the phenomeno. Below we will give a brief summary of the field gical Time-Dependent Ginzburg-Landau(TDGL) approach specially modified to describe the improper kinetic equation [16] seems to be the best choice for martensitic transformations. Particular emphasis will theoretical modeling of martensitic transformations. be given to the formulation of a free energy This approach is used below to describe the functional that incorporates the coherency strait transformation dynamics in the so-called improper energy for the transformation. A good review of the artensitic transformations. The thermodynamic field theory itself can be found in Ref. [16] parameters characterizing an ransformation arc long-range order (lro) par- 2. 1. The field kinetic equation meters which describe the softening displacive The temporal dependence of the field functions on ptical modes causing the transition. The macro- time describes the temporal evolution of the scopic transformation-induced elastic strain is a microstructure. Therefore, if the time-dependent secondary effect caused by its coupling with the Iro equations for those field variables are formulated. the parameters. The improper martensitic transform- problem of the theoretical description of microstruc ations occur mostly in ceramic materials. A good tural dynamics is simply reduced to a solution of and well-investigated example is the martensitic these equations. The simplest form transformation in partially stabilized ZrO The equation of this type is the stochastic Langevin formulation of the "coarse-grained"free equation based on the TDGL. equation which which is the basis for formulating the TDGL assumes that the rate of evolution of a field is a linear equations, is much simpler for an im martensitic transformation than for a proper forces [16] function with respect to the transformation driving nsformation. In the latter (e f c.c. -b.c.c. martensitic transformations in steels) the thermodynamic parameters characterizing the cp(,) +5n(r,D)(2) transformation are the macroscopic strains and the an,(r, t) alytical with respect to the strain gradients where the indexes P, q=l n is a matrIx The main purposes of this re (ite (operator)of the kinetic coefficient, F is the free formulate the first 3D phenomenological field kinetic driving force and 5 (r, 1) is the Langevin noise term model which is able to describe the microstructure which is taken to be Gaussian distributed and its development of martensite without any a priori constraint on the transformation path and (i) to fluctuation-dissipation theorem [18] orrelation properties meet the requirements of the carry out 3D computer simulations of an improper martensitic transformation in a model system by (5,(r, t)s,(r, t')>=2kBTLn,8p,8(r-r)8(t-r) using computationally tractable algorithms. It should be mentioned that the improper martensitic transformation has all the martensitic features where kn is the Boltzmann constant. The origin of the above loIse term is related to the microscopic degrees of
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL 761 computer simulation of phase transformations, but right now even the power of state-of-the-art parallel computers is insufficient to carry out MD simulations for systems of a reasonable size (i.e. big enough to exclude or, at least to make qualitatively insignificant. the size effect) during a time of more than about lo-* s. An alternative approach is the Monte Carlo (MC) simulation. It is an efficient technique for studying replacive transformations involving simple atomic exchanges among rigid lattice sites. Recently, 3D MC simulation has been used to model ordering transformations with a point group symmetry reduction where the transformation strain plays an important role [14. 151. In this approach, the ordering is described by the Ising model and the simulation has demonstrated that the transformation strain accommodation caused by the ordering produces the so-called tweed mesoscopic structure. As far as we know, however, this approach has not been applied to the study of martensitic transformations. At the present stage of computer technology, the continuum field method based on the phenomenological Time-Dependent Ginzburg-Landau (TDGL) kinetic equation [16] seems to be the best choice for theoretical modeling of martensitic transformations. This approach is used below to describe the transformation dynamics in the so-called improper martensitic transformations. The thermodynamic parameters characterizing an improper martensitic transformation are long-range order (lro) parameters which describe the softening displacive optical modes causing the transition. The macroscopic transformation-induced elastic strain is a secondary effect caused by its coupling with the ho parameters. The improper martensitic transformations occur mostly in ceramic materials. A good and well-investigated example is the martensitic transformation in partially stabilized ZrOl. The formulation of the “coarse-grained” free energy, which is the basis for formulating the TDGL kinetic equations, is much simpler for an improper martensitic transformation than for a proper transformation. In the latter case (e.g. the f.c.c. -tb.c.c. martensitic transformations in steels), the thermodynamic parameters characterizing the transformation are the macroscopic strains and the coarse-grained free energy is not analytical with respect to the strain gradients. The main purposes of this paper are (i) to formulate the first 3D phenomenological field kinetic model which is able to describe the microstructure development of martensite without any a priori constraint on the transformation path and (ii) to carry out 3D computer simulations of an improper martensitic transformation in a model system by using computationally tractable algorithms. It should be mentioned that the improper martensitic transformation has all the martensitic features discussed above. 2. FIELD KINETIC MODEL OF MARTENSITIC TRANSFORMATIONS The continuum field approach to martensitic transformations is based on a simple idea that an arbitrary microstructure, consisting of all possible orientation variants of martensite coherently imbedded in its parent phase matrix, can be described by continuous functions (field variables), q,(r). +(r), . q,,(r) where r is the coordinate vector. Each of them describes an lro parameter profile which characterizes the spatial distribution of one of the orientation variants. Thus the total number of functions, II, is equal to the number of all possible orientation variants which are determined by the symmetry of the crystal lattice rearrangement. Each function assumes a non-zero value within martensitic particles of the corresponding orientation variant and vanishes within the parent phase or other orientation variants. Using this approach, Wang et al. have simulated the microstructure evolution during precipitation of tetragonal ZrOz particles from a cubic matrix in partially stabilized zirconia [17]. Below we will give a brief summary of the field approach specially modified to describe the improper martensitic transformations. Particular emphasis will be given to the formulation of a free energy functional that incorporates the coherency strain energy for the transformation. A good review of the field theory itself can be found in Ref. [16]. 2.1. The ,field kinetic equation The temporal dependence of the field functions on time describes the temporal evolution of the microstructure. Therefore. if the time-dependent equations for those field variables are formulated, the problem of the theoretical description of microstructural dynamics is simply reduced to a solution of these equations. The simplest form of kinetic equation of this type is the stochastic Langevin equation based on the TDGL equation which assumes that the rate of evolution of a field is a linear function with respect to the transformation driving forces [ 161: where the indexes p, q = 1,2, , n. &, is a matrix (operator) of the kinetic coefficient, F is the free energy functional, 6F/6q,(r, t) is the thermodynamic driving force and &(r, t) is the Langevin noise term which is taken to be Gaussian distributed and its correlation properties meet the requirements of the fluctuation-dissipation theorem [ 181: (<,(r, t)&(r’, t’)) = 2ksT&&&,6(r - r’)6(t - t’) (3) where ke is the Boltzmann constant. The origin of the noise term is related to the microscopic degrees of
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL freedom, e.g. to the thermal vibrations (phonons), related to the finite-range atomic interaction, the high-order correlations with short relaxation time, strain energy contribution resulting from the etc,which are assumed to be in equilibrium at the infinite-range strain-induced intcraction has to be temperature T to which the system is quenched. considered. Below we formulate both energies The field equation(2)is valid under the condition functionals of the lro parameter fields variables with respect to all other degrees of freedom. 2.2. The strain energy functional Then the averaging with respect to the" fast"degrees It has been shown by Khachaturyan [7, 19]that the of freedom over the time-dependent ensemble gives the strain energy of an arbitrary coherent multiphase non-cquilibrium free energy whose deviation from the system can be presented in a closed form for the quilibrium value is the driving force for the ensemble homogeneous modulus case. The strain energy in evolution. This condition is met for the improper [7, 19] is expressed as a functional of the transform martensitic transformations considered in this paper. ation- induced stress-free strain e(r)which is pre The Iro parameters characterizing the improper sented through an arbitrary shape function 0, (r), e. g. ransformations are soft optical phonon modes which are usually overdamped and whose thermalization time is considerably longer than that of thefast f(r)=∑rygo(m) degrees of freedom which are the regular phonons. It is the averaging over the regular phonons that gives the where 0()equals unity if a point r, is within a ncw non-equilibrium free energy in the theory of improper phase particle of type p and zero outside it, and the transformations in ferroelastic materials tensor Ew (p) describes the stress-free transformation However, it should be mentioned that this is not strain transforming the parent phase into the pth always the case for the proper martensitic transform- orientation variant of the product phase ations, The condition formulated above which is the To incorporate the strain energy into the field basis for the kinetic equation (2)may be violated for kinetic equation (4), we have to modify the theory the strongly undercooled proper martensitic trans- [7, 19]by expressing the stress-free strain through the formations which develop thermally. The Iro continuum field ne(r). Let us consider the most parameters characterizing this type of transformation common case of improper martensitic transform- are the breaking symmetry acoustic strains which are ations where the free energy is invariant with respect not soft modes and whose thermalization time is to the ne transition. This if the ommensurate with the thermalization time of other np +-n transformation is equivalent to one of phonons. In this case, a dynamic approach rather the parent phase symmetry operations. In this case than the approach based on the thermodynamics of the first non-vanishing term of the stress-free strain the non-equilibrium state(e.g. the TDGL equation expansion is quadratic with respect to n (r). Then the and other conventional kinetic equations defined by equation for the stress-free strain becomes averaging over the time-dependent ensemble) should be employed and effects such as the mechanical inertia should be taken into account ()=∑m()(p) The Iro parameters in equation(2)are non-con- served fields because they have no kinematIc Equation (7)actually describes a situation where the constraints and their values are determined by the change ne+-ne does not affect the transformation free cnergy minimization condition only. For strain (for cxample, his change produces an served fields, the first non-vanishing term the long wave approximation for the kinetic According to Ref. [19] the total elastic strain coefficient matrix constant,i.e.Lp=LP. energy counted from the stress-free state can be diagonal, the kinetic equation(2)can be simplified as terms: an (r, n)odF (4) Ee= eyqueen + Erom squeeze"all the stress-free orientation variants of (Sp(r, t)sp(r, 1)>=2k TL(r-r)o(t-t).(s) the martensitic phase back into their parent phase structure or, in other words, it is the energy of elastic To present equation(4)in an explicit analytical form strain required to compensate the transformation-in- itable for a numerical solution, we have to find the duced stress-free strain e(r). Erea and Ereiax are the analytical equation for the thermodynamic driving encrgy reductions due to the macroscopic homo- force 8F/Sn, (r, 1). To do so, we have to formulate the geneous strain relaxation and local heterogeneous ee energy functional F in terms of the Iro strain relaxation, respectively, They are the energy parameters n, first. In martensitic transformations, in changes caused by relaxation of the system from the addition to the "chemical"free energy which is non-strained (squeezed) state into the new elastic
162 WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL freedom, e.g. to the thermal vibrations (phonons), high-order correlations with short relaxation time, etc., which are assumed to be in equilibrium at the temperature T to which the system is quenched. The field equation (2) is valid under the condition that the lro parameters are “slow” relaxing dynamic variables with respect to all other degrees of freedom. Then the averaging with respect to the “fast” degrees of freedom over the time-dependent ensemble gives the non-equilibrium free energy whose deviation from the equilibrium value is the driving force for the ensemble evolution. This condition is met for the improper martensitic transformations considered in this paper. The lro parameters characterizing the improper transformations are soft optical phonon modes which are usually overdamped and whose thermalization time is considerably longer than that of the “fast” degrees of freedom which are the regular phonons. It is the averaging over the regular phonons that gives the non-equilibrium free energy in the theory of improper transformations in ferroelastic materials. However, it should be mentioned that this is not always the case for the proper martensitic transformations. The condition formulated above which is the basis for the kinetic equation (2) may be violated for the strongly undercooled proper martensitic transformations which develop athermally. The lro parameters characterizing this type of transformation are the breaking symmetry acoustic strains which are not soft modes and whose thermalization time is commensurate with the thermalization time of other phonons. In this case, a dynamic approach rather than the approach based on the thermodynamics of the non-equilibrium state (e.g. the TDGL equation and other conventional kinetic equations defined by averaging over the time-dependent ensemble) should be employed and effects such as the mechanical inertia should be taken into account. The lro parameters in equation (2) are non-conserved fields because they have no kinematic constraints and their values are determined by the free energy minimization condition only. For non-conserved fields, the first non-vanishing term in the long wave approximation for the kinetic coefficient matrix is a constant, i.e. i,, = Lj,. Assuming the kinetic coefficient matrix to be diagonal, the kinetic equation (2) can be simplified as and (&(r, f)&,(r’, t’)) = 2keTLoG(r - r’)s(t - t’). (5) To present equation (4) in an explicit analytical form suitable for a numerical solution, we have to find the analytical equation for the thermodynamic driving force 6F/6qp(r, t). To do so, we have to formulate the free energy functional F in terms of the lro parameters qP first. In martensitic transformations, in addition to the “chemical” free energy which is related to the finite-range atomic interaction, the strain energy contribution resulting from the infinite-range strain-induced interaction has to be considered. Below we formulate both energies as functionals of the lro parameter fields. 2.2. The strain energy functional It has been shown by Khachaturyan [7, 191 that the strain energy of an arbitrary coherent multiphase system can be presented in a closed form for the homogeneous modulus case. The strain energy in [7, 191 is expressed as a functional of the transformation-induced stress-free strain 6:(r) which is presented through an arbitrary shape function O,(r), e.g. where O,(r) equals unity if a point r, is within a new phase particle of type p and zero outside it, and the tensor c;“(p) describes the stress-free transformation strain transforming the parent phase into the pth orientation variant of the product phase. To incorporate the strain energy into the field kinetic equation (4), we have to modify the theory [7, 191 by expressing the stress-free strain through the continuum field VP(r). Let us consider the most common case of improper martensitic transformations where the free energy is invariant with respect to the Q,+- nP transition. This is the case if the VP+-% transformation is equivalent to one of the parent phase symmetry operations. In this case, the first non-vanishing term of the stress-free strain expansion is quadratic with respect to qp(r). Then the equation for the stress-free strain becomes Equation (7) actually describes a situation where the change Q, + - Q does not affect the transformation strain (for example, this change produces an antiphase domain shift). According to Ref. [19], the total elastic strain energy counted from the stress-free state can be presented as a sum of three physically distinctive terms: -&I = Esqueczc + E::: + Eh& (8) where E,,, is the energy required to elastically “squeeze” all the stress-free orientation variants of the martensitic phase back into their parent phase structure or, in other words, it is the energy of elastic strain required to compensate the transformation-induced stress-free strain c;(r). Es: and EF,=, are the energy reductions due to the macroscopic homogeneous strain relaxation and local heterogeneous strain relaxation, respectively. They are the energy changes caused by relaxation of the system from the non-strained (squeezed) state into the new elastic
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL quilibrium in a structurally inhomogeneous state ments vanishes on the surfacc. In reality, the (e.g. the transformed two-phase coherent mixture) displacements oscillate around zero on the surface described by spatial distribution of the Iro parameter However, the length of the oscillation are within the same scale as that of the structural heterogeneit Keeping in mind the correspondence between na (r) Therefore the inaccuracy caused by this boundary nd 0(r)in equations(6) and (7)and following condition (e.g. the extra strain caused by the exactly the same line of reasoning as in Ref [19]. we oscillation) is localized within a surface layer whose can express the"squeeze""energy (unrelaxed strain thickness is commensurate with the scale of the )and the relaxation energies as follow structural heterogeneity. It is proportional to the surface area of the sample and thus can be neglected Egquere=22 Cuity()e(p)noe, )comparing with the volume-dependent part of the strain ener It can be seen from equations (9), (10)and (12)that are independent of the microscopic morphology of the structural pattern. They depend on the macroscopic parameter fields ne(r). It is, therefore, morphology where V is the total volume of the system, cwkr is the dependent. It is interesting that all the information related to the morphology of the mesoscopic microstructure is contained in the Iro paramete n(r)d profiles np(r). As for the information related to the intrinsic material characteristics, both elastic and crystallographic, it is imbedded in the characteristic function B(e). B(e)depends the elastic is the volume fraction of all domains of the pth properties of the material through the Green function orientation variant (nn is the equilibrium value of the $(e)and on the crystallography of the crystal lattice Iro parameter at a given temperature) and rearrangement through the transformation strain coefficients ea(p)entering a(P)[see equation(13) e=7∑e(p)o (12) Ihe back Fourier transform of Bo, k/k). the function We(r-r), is the induced pairwise interaction cnergy between finite clements of particles of type p is the macroscopic homogeneous strain describing the and q at points r and r. It should be noted that caused by the transformation. In equation(mody B (k k) has a singularity at k=0, since its limit at k+0 depends on the direction of k. This singularity Beo(e)-co(p)Q2 x(e)okr(q)e, results in a Iyr' asymptote of the back Fourier (13) transform Wp(r), providing an infinitely long-ranged where e=k/k is a unit vector in reciprocal space and dipole dipole like interaction. Such an interaction e; is its ith component, o(P)=cuE (P)and S2,(e) is cannot be reduced to the gradient terms in a Green function tensor which is inverse to the tensor continuum limit. Finally, it should be emphasized that Q(e)=calee quation(11)is obtained as a solution of the elastic equation at given Iro parameter profiles np(r). This r F"(r)expo procedure expresses the elastic strain in terms of the Iro parameter profiles and thus excludes the strain from the explicit form of the strain energy. Therefore, is the Fourier transform of mp(r)and (nd(r)* is the the strain energy in the forms of equations(8)(11) complex conjugate of n(r)Jk. In equation( 13), the becomes a functional of the Iro parameter profiles Einstein suffix notation is used, i. e. when a suffix only occurs twice in the same term, it indicates summation Below we consider a particular case of with respect to that suffix. The integration in bic-tetragonal transformation which is quite com- equation( 1)is carried out over the infinite reciprocalmon martensitic transformations. TH space and the sign indicates that a volume of(2x)/v transformation produces three orientation variants of about k =0 is excluded from the integration. when the product martensitic phase and thus can be y is asymptotically large, this exclusion defines the described by three Iro parameter fields, n,(r), m(r)and principle value"of the integral ns(r). The crystallography of the martensitic crystal The strain energies( 8)(11)are obtained by solving lattice rearrangement enters the theory through three the elastic equation for a macroscopically homo- tensor constants e(p)(p=1. 2. 3)in the elastic geneous unconstrained body using a boundary energy equations. These constants are related to the condition that the " heterogeneous"part of displace- corresponding Iro parameters and are determined by
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL 763 equilibrium in a structurally inhomogeneous state (e.g. the transformed two-phase coherent mixture) described by spatial distribution of the lro parameter fields. Keeping in mind the correspondence between $(r) and 9,,(r) in equations (6) and (7) and following exactly the same line of reasoning as in Ref. [19], we can express the “squeeze” energy (unrelaxed strain energy) and the relaxation energies as follows EL, = -i z 5% 4 (k){$ WM$ WiZ (11) where V is the total volume of the system, cllkl is the elastic modulus tensor, s vj W d3r wp = dV is the volume fraction of all domains of the pth orientation variant (Q is the equilibrium value of the lro parameter at a given temperature) and is the macroscopic homogeneous strain describing the macroscopic shape change of a non-constrained body caused by the transformation. In equation (11) &(e) = e&~)Q,4e)&q)ei (13) where e = k/k is a unit vector in reciprocal space and e, is its ith component, oz@) = Cam&@) and Q(e) is a Green function tensor which is inverse to the tensor Q(e); ’ = c,ki,eke/ , ($ (r)}k = s d3k m q:(r) exp( - ikr) is the Fourier transform of q;(r) and {qgZ(r)}f is the complex conjugate of {qi(r)JL. In equation (13), the Einstein suffix notation is used, i.e. when a suffix occurs twice in the same term, it indicates summation with respect to that suffix. The integration in equation (11) is carried out over the infinite reciprocal space and the sign j indicates that a volume of (2n)‘/ V about k = 0 is excluded from the integration. When V is asymptotically large, this exclusion defines the “principle value” of the integral. The strain energies (8x11) are obtained by solving the elastic equation for a macroscopically homogeneous unconstrained body using a boundary condition that the “heterogeneous” part of displacements vanishes on the surface. In reality, the displacements oscillate around zero on the surface. However, the length of the oscillation are within the same scale as that of the structural heterogeneity. Therefore, the inaccuracy caused by this boundary condition (e.g. the extra strain caused by the oscillation) is localized within a surface layer whose thickness is commensurate with the scale of the structural heterogeneity. It is proportional to the surface area of the sample and thus can be neglected comparing with the volume-dependent part of the strain energy. It can be seen from equations (9), (10) and (12) that the “squeeze” and homogeneous relaxation energies are independent of the microscopic morphology of the structural pattern. They depend on the macroscopic characteristics, w,, only. However, the heterogeneous strain relaxation energy (11) is a functional of the lro parameter fields q,(r). It is, therefore, morphologydependent. It is interesting that all the information related to the morphology of the mesoscopic microstructure is contained in the lro parameter profiles qp(r). As for the information related to the intrinsic material characteristics. both elastic and crystallographic, it is imbedded in the characteristic function &(e). B,,(e) depends on the elastic properties of the material through the Green function Q,(e) and on the crystallography of the crystal lattice rearrangement through the transformation strain coefficients t:(P) entering 0:@) [see equation (13)]. The back Fourier transform of B,,(k/k), the function wf(r - r’), is the strain-induced pairwise interaction energy between finite elements of particles of type p and q at points r and r’. It should be noted that B,,(k/k) has a singularity at k = 0, since its limit at k+O depends on the direction of k. This singularity results in a l/r-’ asymptote of the back Fourier transform KY(r), providing an infinitely long-ranged dipole-dipole-like interaction. Such an interaction, cannot be reduced to the gradient terms in the continuum limit. Finally, it should be emphasized that equation (11) is obtained as a solution of the elastic equation at given lro parameter profiles VP(r). This procedure expresses the elastic strain in terms of the lro parameter profiles and thus excludes the strain from the explicit form of the strain energy. Therefore, the strain energy in the forms of equations (8)( 11) becomes a functional of the lro parameter profiles only. Below we consider a particular case of cubic+tetragonal transformation which is quite common in martensitic transformations. The transformation produces three orientation variants of the product martensitic phase and thus can be described by three lro parameter fields, q,(r), q?(r) and q,(r). The crystallography of the martensitic crystal lattice rearrangement enters the theory through three tensor constants c:(p)@ = 1,2, 3) in the elastic energy equations. These constants are related to the corresponding ho parameters and are determined by
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL the crystal I parameters of the parent and The local specific free energy f(nl, nz,..., n,)in martensitIc as well as by the orientation equation(15) defines the basic bulk thermodynamic erties of the system. It should be Cartesian system axes are parallel to the cubic respect to all rotations and reflections of symmetry crystallographic axes of the parent phase, these operations of the parent phase. It can be approxi- tensors have the following forms: mated by the Landau polynomial expansion with respect to the lro parameters which form independent symmetry invariants [18]. Rich applications of the Landau theory for the theoretical characterization of phase transformations in ferroelastic materials can be found in Refs [20] and [21]. For the cubi tetragonal transformation considered in this paper, E8(3)=0610 (14)f((n, n2,,, nn)=f(n, 2, n3 ). The polynomial ap- proximation of f(n, n2, n3 )up to the sixth-order terms in n, n2 and ns forming independent symmetry variants of the appropriate order is n,n,)=2(7+n+)+4++m Cr and a, are the crystal lattice parameters of the ( ni+n2+n3)+2ninin3 equilibriun martensitic phase at the stress-free state and a is the crystal lattice parameter of the equilibrium parent phase at the stress-free state. The (+n+ stress-free transformation strain is described by uation(7) with p= 1, 2, 3 and e (p)given by +闭n+n+n)》+ equation (14). 23. The chemical free energy The elastic strain energy [equation( 8)] introduced n2 -ne symmetry(this is usually the case for most and hence it increases the total free energy in an absence of all odd terms in equation(16)since driving force that makes the transformation occur. this transition. To describe a first-order transition The non-equilibrium chemical free energy can be from a metastable parent phase to a stable tetragonal approximated by the conventional Ginzburg-Landau phase, a choice of the coefficients A, in equation (16) phenomenological coarse-grained functional of the should provide a local minimum at Iro parameter fields. It contains terms depending on the local specific free energy f(n, nz,..., ne) and (17a) non-local gradient terms, i.e corresponding to the metastable parent phase and sixfold degenerated global minima at F-2 i A(p)v m, VMn η=土n,n2=n=0 72=士7,==0 +f(m,n,…,n)dr(15) ±阳,η=2 where not is the equilibrium Iro parameter. The Bu(p)are coefficients at the gradient terms ter minima correspond to three orientation variants of the stable martensitic phase as well as two arc positively dcfined sccond-rank tensors, i, possible antiphase states for each of the orientation are indexes of Cartesian coordinates and vi= d/ ar is a differential operator. The integration in equation variants distinguished by the signs at the Iro al, the gradient terms in equation( 15) provide stability of the metastable parent phase against min( f(n, nz,n3))=f(+ o, 0, 0)=f(o, +no, 0) infinitesimal Auctuations of Iro parameters. The ymmetry of the gradient coefficient tensors B, is f(,0,±n) determined by the symmetry change associated with Substituting equations (14)into equation(7)and the crystal lattice rearrangement during the phase using the minimizing Iro parameters [equations transformations. It characterizes the anisotropy of (17b d)], one can find that the stress-free transform the interfacial and or domain boundary energie ation strain tensors corresponding to these par
764 WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL the crystal lattice parameters of the parent and martensitic phases as well as by the orientation relationship between them. For example, if the Cartesian system axes are parallel to the cubic crystallographic axes of the parent phase, these tensors have the following forms: El 0 0 c?(3) = 0 tl 0 ( i (14) 0 0 63 where 63 = ct - a, 2 wk ct and a, are the crystal lattice parameters of the equilibrium martensitic phase at the stress-free state and a, is the crystal lattice parameter of the equilibrium parent phase at the stress-free state. The stress-free transformation strain is described by equation (7) with p = 1, 2, 3 and t;(p) given by equation (14). 2.3. The chemical free energy The elastic strain energy [equation (8)] introduced by a martensitic transformation is always positive and hence it increases the total free energy suppressing the transformation. It is the “chemical” driving force that makes the transformation occur. The non-equilibrium chemical free energy can be approximated by the conventional Ginzburg-Landau phenomenological coarse-grained functional of the lro parameter fields. It contains terms depending on the local specific free energy f(r,, Q, . . , a), and non-local gradient terms, i.e. F= S[ ” ffh, m, . . . , m) d3r (15) 1 where p,-(p) are coefficients at the gradient terms which are positively defined second-rank tensors, i, j are indexes of Cartesian coordinates and Vi 3 d/ar, is a differential operator. The integration in equation (15) is carried out over the entire system volume V. As usual, the gradient terms in equation (15) provide stability of the metastable parent phase against infinitesimal fluctuations of ho parameters. The symmetry of the gradient coefficient tensors /Iti is determined by the symmetry change associated with the crystal lattice rearrangement during the phase transformations. It characterizes the anisotropy of the interfacial and/or domain boundary energies. The local specific free energy f(~, Q, . . . , qn) in equation (15) defines the basic bulk thermodynamic properties of the system. It should be invariant with respect to all rotations and reflections of symmetry operations of the parent phase. It can be approximated by the Landau polynomial expansion with respect to the lro parameters which form independent symmetry invariants [18]. Rich applications of the Landau theory for the theoretical characterization of phase transformations in ferroelastic materials can be found in Refs [20] and [21]. For the cubic+ tetragonal transformation considered in this paper, &I, v2, . . . , q.) =fh,q2,yl3). The polynomial approximation of f(q,, q2, n3) up to the sixth-order terms in ~1, q2 and n3 forming independent symmetry invariants of the appropriate order is + +4: + rl: + vi)' f.. . (16) where Ai(i = 1,2, . . ) are the expansion coefficients at fixed temperature and pressure. The assumed Q-+ -r?, symmetry (this is usually the case for most of the improper martensitic transformations) results in an absence of all odd terms in equation (16) since the free energy should be invariant with respect to this transition. To describe a first-order transition from a metastable parent phase to a stable tetragonal phase, a choice of the coefficients Ai in equation (16) should provide a local minimum at ?I = r/z = qs = 0 (17a) corresponding to the metastable parent phase and sixfold degenerated global minima at 91 = *no, @ = 'I3 = 0 (17b) q2 = *f/o, VI = q3 = 0 (174 q3 = +ro, Q=?z=O (174 where no # 0 is the equilibrium lro parameter. The latter minima correspond to three orientation variants of the stable martensitic phase as well as two possible antiphase states for each of the orientation variants distinguished by the f signs at the lro parameter. This sixfold degeneration is reflected by the following equalities: min{f(rll,q2, 13>> =f(+rbO,O) =f@, +VO, 0) =m 0, _frlo). (18) Substituting equations (14) into equation (7) and using the minimizing lro parameters [equations (17&d)], one can find that the stress-free transformation strain tensors corresponding to these par-
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL ameters are tetragonal and provide all three possible rientation variants. It should be noted that the e8(3)=60-10 (20) prcscntation of free energy [equation(16)]as series expansion with respect to the cubic symmetry invariants of the Iro parameters n. n and ns automatically provides the correct degener- They describe a simple shear lattice distortion ation of the global minimum and thus the involving expansion along one of the three cubic axes correct number of orientation variants of the and compensating co actions along the other two martensite ents the magnitude of the shear The simplest form of the polynomial free It is well documented that martensitic transform- energy approximation meeting all the above require- ations in bulk materials always produce all orien- tation variants of the low-symmetry martensitic phase so that the average macrosco try of f(n.z, n3)=5A(ni+m2+ n:) the Am+n+m)+A团+n+n)(19) cubic-tetragonal transformation, this means that where At, A2 and A, are positive constants depending the macroscopic retry of the system after the on temperature and composition of the syster transformation remains cubic. In this case, the Substituting the above formulated elastic strain macroscopic strain (0 m(P )5, since any second- nergy and chemical free energy into kinetic equation rank tensor in a cubic lattice is proportional to the (4)and taking the variational derivative with respect Kronecker symbol 8y. For the pure shear type of to np(r) gives an explicit general form of the field stress-free crystal lattice rearrangement considered tic equation for the cubic+tetragonal martensitic here the trace of the stress-free transformation strain transformation, The model is in principle applicabl Emn (P) is zero. Therefore the macroscopic strain is for any transformations with any arbitrary crystal zero as well. The same argument and conclusion are lattice rearrangement in elastically anisotropic crys- applicable to a polycrystalline material without texture. In this case, the macroscopic symmetry of the system is spherical, which results in a pure 3. APPLICATION TO A MODEL SYSTEM dilatational macroscopic deformation Ey- du. Since he assumed transformation model does not consider In the field approach formulated above, the the dilatational effect [Emm (p)=0]. 60=0. Therefore problem of theoretical characterization of martensitic the term(10)in the free energy can be neglected transformations has heen reduced to finding solutions Neglecting the anisotropy of the interfacial energy of the corresponding field kinetic equations for all we may assume that the gradient coefficient tensor, possible orientation variants of the martensitic phase. Bu, in equation (15)is isotropic, e.g. B,=B8 ntegro-difference equations. Below we show how to and chemical free energy equations and taking the solve those equations numerically using a Fourier variational derivative of the free energies with respect transformation technique for a 3D prototype model to the Iro parameter fields, the field kinetic equation system. The system undergoes a cubic-tetragonal for the cubic, tetragonal transformation becomes martensitic transformation below the transition an,(r, 1) It is typical for all structural transformations hat the shear part of the stress-free transformation strain e (p) is much greater than its dilatational rt that is associated with the volume change of =-L0-B△n2(r,1)+ af--2np(r, Onp(r, I) dilatational part is zero and thus the martensitic transformation develops without volume chan This assumption is equivalent to a requirement that Trace(eg (p))=e m (p)=0. Then the stress-free strains lequations(14)become (21) 0 e"(1)=0 0 wheref is given by equation (19) We can further simplify the numerical calculation 100 should be, however, emphasized that this assumption as well as the assumption that the trace of the (2)=6020 tress-free transformation strain is zero results in no loss of the essential characteristics of the
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL 765 ameters are tetragonal and provide all three possible orientation variants. It should be noted that the presentation of free energy [equation (16)] as a series expansion with respect to the cubic symmetry invariants of the lro parameters nI, Q and q3 automatically provides the correct degeneration of the global minimum and thus the correct number of orientation variants of the martensite. The simplest form of the polynomial free energy approximation meeting all the above requirements is - $ A,($ + r/4 + q;) + ; A,($ + r/t + SS)’ (19) where A,, AZ and A, are positive constants depending on temperature and composition of the system. Substituting the above formulated elastic strain energy and chemical free energy into kinetic equation (4) and taking the variational derivative with respect to I gives an explicit general form of the field kinetic equation for the cubic-+tetragonal martensitic transformation. The mode1 is in principle applicable for any transformations with any arbitrary crystal lattice rearrangement in elastically anisotropic crystals. 3. APPLICATION TO A MODEL SYSTEM In the field approach formulated above, the problem of theoretical characterization of martensitic transformations has been reduced to finding solutions of the corresponding field kinetic equations for all possible orientation variants of the martensitic phase. These equations form a set of coupled nonlinear integro-difference equations. Below we show how to solve those equations numerically using a Fourier transformation technique for a 3D prototype model system. The system undergoes a cubic+tetragonal martensitic transformation below the transition temperature. It is typical for all structural transformations that the shear part of the stress-free transformation strain $‘@) is much greater than its dilatational part that is associated with the volume change of the system. For simplicity we assume that the dilatational part is zero and thus the martensitic transformation develops without volume change. This assumption is equivalent to a requirement that Trace{ci!(p)} = E,$,@) = 0. Then the stress-free strains [equations (14)] become 20 0 t?(l) = 6” 0 -1 0 i i ) 0 0 -1 -1 0 0 CF(2)=Ea 0 2 0 ( i i 0 0 -1 -1 0 0 tY(3)=to t 0 -1 0. i (20) 0 0 2 They describe a simple shear lattice distortion involving expansion along one of the three cubic axes and compensating contractions along the other two axes; c0 represents the magnitude of the shear. It is well documented that martensitic transformations in bulk materials always produce all orientation variants of the low-symmetry martensitic phase so that the average macroscopic symmetry of the system coincides with the point group symmetry of the parent phase. In the specific case of a cubic+tetragonal transformation, this means that the macroscopic symmetry of the system after the transformation remains cubic. In this case, the macroscopic strain Ci - ~:~@)6, since any secondrank tensor in a cubic lattice is proportional to the Kronecker symbol 6,. For the pure shear type of stress-free crystal lattice rearrangement considered here, the trace of the stress-free transformation strain C::,“,(P) is zero. Therefore the macroscopic strain is zero as well. The same argument and conclusion are applicable to a polycrystalline material without texture. In this case, the macroscopic symmetry of the system is spherical, which results in a pure dilatational macroscopic deformation Fi - 6,, Since the assumed transformation model does not consider the dilatational effect [E:~@) = 01, FP, = 0. Therefore, the term (10) in the free energy can be neglected. Neglecting the anisotropy of the interfacial energy, we may assume that the gradient coefficient tensor, /r ‘i 3 in equation (15) is isotropic, e.g. /I, = PS,,. Substituting $ = 0 and fi, = /-?S, into the strain energy and chemical free energy equations and taking the variational derivative of the free energies with respect to the lro parameter fields. the field kinetic equation for the cubic-+tetragonal transformation becomes: = -Lo af -pAqP(r, t) + ~ a%(r, t) - 2gp(r, t) x;j&m(;) {f?,‘(r, t))k exptikr) + &(r, t) (21) where f is given by equation (19). We can further simplify the numerical calculation by assuming that the system is elastically isotropic. It should be, however, emphasized that this assumption as well as the assumption that the trace of the stress-free transformation strain is zero results in no loss of the essential characteristics of the martensitic
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL transformation. Besides, the numerical simulation and without these assumptions does not pose any principal difficulties <Sp(r", t)sp(r*,t)=2 o(r-r")(t-r) For an elastically ic medium. the green function in equation(13)bece Since, according to equations(20) and(23), B(e)is proportional to GEo, it turns out that all the (22) dimensional material constants entering the simu- lation are included in just two dimensionless where G is the shear modulus and y is Poissons ratio arameters: 5=GE/A/ and =B/P/4n. The parameter s discussed in the beginning of the paper Substituting equation (22) into (13)and using characterizes the""strength"of the martensitic equations(20) in the definition of oi(p)=coek(p) transformation while the parameter By(e)=4e"(p)W(q)2-2G1 X (ee (P)e, ee(p)er ].(23) characterizes the contribution of the interfaci energy between orientation variants and between It follows from equations (20)and (23)that the phases to the system thermodynamics and kinetics dimensional part of the strain energy is Ge It is convenient to present the kinetic equation simulated microstructure. For example, if we know in a dimensionless form for the numerical solution This will also allow us to reveal the key physical transformation, the interfacial and twin boundary do this, we have introduced a reduced time, defined energies for a particular system, then the physical length of the computational grid increment I as t= LolA/lt, and reduced spatial coordinates, equation(26)can be directly calculated from the is the transformation driving force [Af Imin f(ni, nz, 3)l=U(no, 0, O)1] and I is the length n our computer simulation, equation(24)is solved numerically in reciprocal space using a periodical boundary condition along all three dimensions. It has Substituting these reduced variables into equation been assumed that the size and macroscopic shape of (21), we get the dimensionless form of the kinetic the system do not change during the transformation This is actually equivalent to a boundary condition in(r*,t) hat the displac This boundary condition is met if the transformation occurs in single grains of a polycrystalline body or in a unconstrained single-crystalline system where the --φ△n(°,)+ transformation does not produce a volume effect and maintains the macroscopic homogeneity and sym metry of the body(this is the case if the sizes of the d'K martensitic particles are much smaller than the size of J(2n)B 2* n(r* t)]k exp(ik*r") the system and if all orientation variants are equally present) +5(r*,T) (24) To provide better visualization of the multidomain and multi-variant microstructure through black-and we have made a simplification by permitting the appearance of only two orientation r+=(p1,p2,p3),k=(K,K2,K), variants, i.e. we maintain na(r, t)=0 during the nulation. In this case, the macroscopic shape of the orientation variants. Therefore, our boundary con A, is the Laplacian in p dition assuming a lack of external boundary fa=ja(ni+n2+n3)-ta20ni+n2+n3) displacements actually means that the system is constrained. The transformation raises an additional +5a3(ni+2+ ni)(25) clastic strain energy resulting from the constrained boundary and thus a larger chemical d where following values of the dimensionless constants have been chosen:{=0.056,中=0.01,a1=0.02 and ay=0. 1. One may note that the value of s chosen
766 WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL transformation. Besides, the numerical simulation without these assumptions does not pose any principal difficulties. For an elastically isotropic medium, the Green function in equation (13) becomes (22) where G is the shear modulus and v is Poisson’s ratio. Substituting equation (22) into (13) and using equations (20) in the definition of c@) = cijk&‘(P) gives 1 B,(e) = 4Gei$‘(p)~,?(q)e, - 2G - l-v x [ei2”(p)ejl[ek~kg,(p)e,l. (23) It follows from equations (20) and (23) that the dimensional part of the strain energy is Gci. It is convenient to present the kinetic equation in a dimensionless form for the numerical solution. This will also allow us to reveal the key physical parameters that control the system behavior. To do this, we have introduced a reduced time, defined as z = LolAAt, and reduced spatial coordinates, defined as pi = xi/l. In these definitions, [Afl is the transformation driving force [IAfl = lminf(q, ~2, ~41 = Ifh, 0,O)ll and 1 is the length scale assigned to the computational grid increment. Substituting these reduced variables into equation (21), we get the dimensionless form of the kinetic equation: + tXr*, 7) (24) where r* = (~1, ~2, p3), k* = (ICI, ~2, ~31, K=kl, 4=f$ Ap is the Laplacian in pi, + ; a&: + q: + $)’ (25) where L?=&, a2=&. a3 =& and (5;(r*, z)T;(r*‘, 7’)) = 2 $$ S(r* - r*‘)6(z - z’). Since, according to equations (20) and (23), B,,(e) is proportional to GE:, it turns out that all the dimensional material constants entering the simulation are included in just two dimensionless parameters: 5 = Gti/lAfl and 4 = p/l’/IAfl. The parameter i discussed in the beginning of the paper characterizes the “strength” of the martensitic transformation while the parameter (26) (where lcapillar = JpI14fl is the capillarity length) characterizes the contribution of the interfacial energy between orientation variants and between phases to the system thermodynamics and kinetics. The value of 4 determines the length scale of the simulated microstructure. For example, if we know the chemical driving force of the martensitic transformation, the interfacial and twin boundary energies for a particular system, then the physical length of the computational grid increment 1 in equation (26) can be directly calculated from the coefficient 4 used in the simulation. In our computer simulation, equation (24) is solved numerically in reciprocal space using a periodical boundary condition along all three dimensions. It has been assumed that the size and macroscopic shape of the system do not change during the transformation. This is actually equivalent to a boundary condition that the displacements at the body surface are zero. This boundary condition is met if the transformation occurs in single grains of a polycrystalline body or in a unconstrained single-crystalline system where the transformation does not produce a volume effect and maintains the macroscopic homogeneity and symmetry of the body (this is the case if the sizes of the martensitic particles are much smaller than the size of the system and if all orientation variants are equally present). To provide better visualization of the multidomain and multi-variant microstructure through black-andwhite pictures, we have made a simplification by permitting the appearance of only two orientation variants, i.e. we maintain q3(r, 2) = 0 during the simulation. In this case, the macroscopic shape of the system changes due to the presence of only two orientation variants. Therefore, our boundary condition assuming a lack of external boundary displacements actually means that the system is constrained. The transformation raises an additional elastic strain energy resulting from the constrained boundary and thus a larger chemical driving force has to be used. For the computer simulation, the following values of the dimensionless constants have been chosen: [ = 0.056, 4 = 0.01, aI = 0.02, a2 = 0.5 and a3 = 0.1. One may note that the value of i chosen
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL here is relatively small compared to the criterion switching off the noise term. The nine small squares discussed at the beginning of the paper. This is presented in each micrograph of Fig. I are the because a greater chemical driving force lAf in consecutive equally spaced (010)cross-sections of the equation(1)is required for the transformation under 3D cubic computational cell with 64 x 64 x 64 mesh the constrained condition points of the cubic grid. They are arranged in such an The size of our computational cell IS 64 x 64 x 64. orientation that the edges of the squares are parallel It turns out that the size effect and the inaccuracy to the [100] and [00l] directions, respectively introduced by the assumption of macroscopic Since the strain accommodation is the dominant homogeneity, which are both associated with the factor controlling a martensitic tra ansformation, the comparatively small system volume, are not too strain-induced interaction between structural hetero gnificant to prevent reproducing the major features geneities should drastically affect the nucleation path of the martensitic structure The simulation result presented in Fig. I(a)sh that this is exactly the case. The critical nuclei generated by the stochastic random noise are 4. RESULTS AND DISCUSSION multidomain particles consisting of two twin-related The simulation is started with an initial conditio orientation variants, Such a transform path is ery different from that assumed by the exis Sung to a homogeneous cubic solid solution which, nucleation theories of martensitic transformation according to equation (25) with the chosen pa ameters for a1, ay and a,, is a metastable state particle event. I his simulation result is actually not corresponding to a local free energy minimum. The surprising because the internally twinned structure of the embryos considerably reduces the strain energy simulation results, i.e. the solutions of equation(24), Formation of the polytwin structure shown in are presented in Figs 1-6. They show the time Fig. l(a) is a collective phenomenon. This has been volution of the Iro parameter profiles, which describes the microstructure development of marten. Confirmed by a specially designed computer exper site during the transformation. In those 2D iment. For example, we have performed a parallel simulation using exactly the same parameters as those black-and-white micrographs, the shades of gray employed in the above simulation but with only one the higher the value, the brighter the shade. orientation variant. No nucleation event takes place Therefore, the black and white domains represent even with extended fluctuation time. It is therefore two different orientation va of the tetrag necessary to have at least two orientation variants martensitic phase and the gray background rep. involved in a single nucleation event. We may call resents the cubic matrix. In those 3D pictures the t such a nucleation process a collective nucleation orientation variants are represented by either However, if the martensitic transformation IS not two different colors, e.g. blue and green(Fig. 3)or ery "strong"(e.g. the parameter c is small). two difterent gray levels e.g. black and gray(Fig. 5). by a step-by-step sympathetic nucleation of orientation variants which are twin- relate other 4. 1. Nucleation of martensite The above result is interesting since the mechan Since the initial cubic phase is in a metastable state, isms of nucleation and the structures of critical nuclei the martensitic transformation occurs through a of martensite are still a subject of intensive nucleation process driven by the Langevin noise term discussions. Several estimates(see, e.g., Cohen and in the kinetic equation which simulates the thermal Wayman[22])have resulted in a conclusion that the fluctuations. Because this stochastic approach does activation energy required for homogeneous nucle not impose any constraint on the transformation ation of a single variant martensitic particle in perfect path, it allows the system tochoose"the optimal bulk crystal is far too high (- 10'kBT)to be overcom evolution path itself and thus gives the most general by thermal fluctuations. This conclusion seems to rule realization of a critical structural heterogeneity out any chance of homogeneous nucleation of a rresponding to a critical nucleus state which is single-variant particle. However, at least one reliable consistent with the characteristics of the free observation of the homogeneous nucleation To distinguish the critical"and "non-critical" martensite within Fe-Co precipitates formed in a heterogeneities, we have artificially"turned off"the Cu-based matrix has been reported by Lin et al.[23) noise term in equation (24)after a certain period of This contradiction could be resolved if the adaptive ime. Then, all the lro parameter heterogeneities phase nucleation model [1] is assumed. According to corresponding to the "sub-critical"fluctuations Ref [1]. the strain energy is drastically reduced if the would eventually disappear and the survived" nucleus is a particle of a transient adaptive phase heterogeneities that are able to grow are actually the whose atomic structure is formed by periodically critical or operational nuclei. Figure I(a)shows the alternating atomic scale microtwins of the martensitic morphology of the survived nuclei at t =2 after phase lattice. The thicknesses of the microtwins in the
WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL 161 here is relatively small compared to the criterion discussed at the beginning of the paper. This is because a greater chemical driving force ]Afl in equation (1) is required for the transformation under the constrained condition. The size of our computational cell is 64 x 64 x 64. It turns out that the size effect and the inaccuracy introduced by the assumption of macroscopic homogeneity, which are both associated with the comparatively small system volume, are not too significant to prevent reproducing the major features of the martensitic structure. 4. RESULTS AND DISCUSSION The simulation is started with an initial condition q,(r*, 0) = q2(r*, 0) = q3(r*, 0) = 0. This corresponds to a homogeneous cubic solid solution which, according to equation (25) with the chosen parameters for a,, a2 and u3, is a metastable state corresponding to a local free energy minimum. The simulation results, i.e. the solutions of equation (24), are presented in Figs l-6. They show the time evolution of the lro parameter profiles, which describes the microstructure development of martensite during the transformation. In those 2D black-and-white micrographs, the shades of gray represent the values of ]q,(r*, r)] - Iq2(r*, z)l, e.g. the higher the value, the brighter the shade. Therefore, the black and white domains represent two different orientation variants of the tetragonal martensitic phase and the gray background represents the cubic matrix. In those 3D pictures, the two orientation variants are represented by either two different colors, e.g. blue and green (Fig. 3) or two different gray levels, e.g. black and gray (Fig. 5). The retaining cubic phase is transparent. 4.1. Nucleation of martensite Since the initial cubic phase is in a metastable state, the martensitic transformation occurs through a nucleation process driven by the Langevin noise term in the kinetic equation which simulates the thermal fluctuations. Because this stochastic approach does not impose any constraint on the transformation path, it allows the system to “choose” the optimal evolution path itself and thus gives the most general realization of a critical structural heterogeneity corresponding to a critical nucleus state which is consistent with the characteristics of the free energy. To distinguish the “critical” and “non-critical” heterogeneities, we have artificially “turned off” the noise term in equation (24) after a certain period of time. Then, all the lro parameter heterogeneities corresponding to the “sub-critical” fluctuations would eventually disappear and the “survived” heterogeneities that are able to grow are actually the critical or operational nuclei. Figure l(a) shows the morphology of the survived nuclei at 7 = 2 after switching off the noise term. The nine small squares presented in each micrograph of Fig. 1 are the consecutive equally spaced (010) cross-sections of the 3D cubic computational cell with 64 x 64 x 64 mesh points of the cubic grid. They are arranged in such an orientation that the edges of the squares are parallel to the [loo] and [OOl] directions, respectively. Since the strain accommodation is the dominant factor controlling a martensitic transformation, the strain-induced interaction between structural heterogeneities should drastically affect the nucleation path. The simulation result presented in Fig. I(a) shows that this is exactly the case. The critical nuclei generated by the stochastic random noise are multidomain particles consisting of two twin-related orientation variants. Such a transformation path is very different from that assumed by the existing nucleation theories of martensitic transformation where the nucleation is considered as a single variant particle event. This simulation result is actually not surprising because the internally twinned structure of the embryos considerably reduces the strain energy. Formation of the polytwin structure shown in Fig. l(a) is a collective phenomenon. This has been confirmed by a specially designed computer experiment. For example, we have performed a parallel simulation using exactly the same parameters as those employed in the above simulation but with only one orientation variant. No nucleation event takes place even with extended fluctuation time. It is therefore necessary to have at least two orientation variants involved in a single nucleation event. We may call such a nucleation process a collective nucleation. However, if the martensitic transformation is not very “strong” (e.g. the parameter [ is small), formation of the polytwin structure may also occur by a step-by-step sympathetic nucleation of adjacent orientation variants which are twin-related to each other. The above result is interesting since the mechanisms of nucleation and the structures of critical nuclei of martensite are still a subject of intensive discussions. Several estimates (see, e.g., Cohen and Wayman [22]) have resulted in a conclusion that the activation energy required for homogeneous nucleation of a single variant martensitic particle in perfect bulk crystal is far too high (N 105kB 7) to be overcome by thermal fluctuations. This conclusion seems to rule out any chance of homogeneous nucleation of a single-variant particle. However, at least one reliable observation of the homogeneous nucleation of martensite within Fe-Co precipitates formed in a Cu-based matrix has been reported by Lin et al. [23]. This contradiction could be resolved if the adaptive phase nucleation model [l] is assumed. According to Ref. [ 11, the strain energy is drastically reduced if the nucleus is a particle of a transient adaptive phase whose atomic structure is formed by periodically alternating atomic scale microtwins of the martensitic phase lattice. The thicknesses of the microtwins in the
ANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL adaptive phase are of several interplanar distances. adaptive phase. Of course, being constrained by the The ratio of the volume fraction of the twin continuum model, we are unable to reproduce the mponents is self-adjusted to provide the stress-frce real atomic structure of the martensitic phase invariant plane strain. This transient phase structure However, we can show that the scale of the transforms into the conventional martensite structure microtwins is commensurate with the(110)interpla during growth. The adaptive phase structure, in fact, nar distances. For example, if we assign a value of is a miniaturized version of the conventional the order of I nm to the grid increment length I in train-accommodating polytwinned martensite whose the definition of in equation (26), then the structure consists of periodically alternating twin-re- chosen parameter =0.01 is obtained lated lamellae of two orientation variants. The 7R materials parameter B/An=apiary 1A[see and 9R martensites are examples of such adaptive definition (26)). This is a physically reasonable value for a typical martensitic phase since VB/ ur computer simulation results lend support to 4f 7/4f where y is the surface energy of the nuclei obtained in the simulation and shown in is of the order of 10? cal/mole. In futuro the adaptive phase nucleation model, The twinned twins, which is of the order of 10 erg/ cm?' and An ig. l(a) can always be interpreted as nuclei of an going to investigate the heterogeneous nucleation Fig. 1. Simulated microstructure development of martensite during isothermal transformation through acleation. The nine small squares in each micrograph are the consecutive 2D cross-sections of a 3D cube along the [010]axis. The black and white dumains represent two different orientation variants of the tetragonal martensitic phase and the gray background represents the cubic matrix;(aH(d correspond to the reduced times t= 2, 6, 10 and 20
168 WANG and KHACHATURYAN: THREE-DIMENSIONAL FIELD MODEL adaptive phase are of several interplanar distances. The ratio of the volume fraction of the twin components is self-adjusted to provide the stress-free invariant plane strain. This transient phase structure transforms into the conventional martensite structure during growth. The adaptive phase structure, in fact, is a miniaturized version of the conventional strain-accommodating polytwinned martensite whose structure consists of periodically alternating twin-related lamellae of two orientation variants. The 7R and 9R martensites are examples of such adaptive phases. Our computer simulation results lend support to the adaptive phase nucleation model. The twinned nuclei obtained in the simulation and shown in Fig. l(a) can always be interpreted as nuclei of an adaptive phase. Of course, being constrained by the continuum model, we are unable to reproduce the real atomic structure of the martensitic phase. However. we can show that the scale of the microtwins is commensurate with the (110) interplanar distances. For example, if we assign a value of the order of 1 nm to the grid increment length 1 in the definition of 4 in equation (26) then the chosen parameter C$ = 0.01 is obtained if the materials parameter m = lcapillary - 18, [see definition (26)]. This is a physically reasonable value for a typical martensitic phase since J/I/ ]Af] - y/]Af] where y is the surface energy of the twins, which is of the order of 10 erg/cm2 and ]Affl is of the order of 102cal/mole. In future, we are going to investigate the heterogeneous nucleation Fig. 1. Simulated microstructure development of martensite during isothermal transformation through homogeneous nucleation. The nine small squares in each micrograph are the consecutive 2D cross-sections of a 3D cube along the [OlO] axis. The black and white domains represent two different orientation variants of the tetragonal martensitic phase and the gray background represents the cubic matrix; (ak(d) correspond to the reduced times r = 2, 6, 10 and 20