A novel 3D mixed-mode multigrain model with efficient implementation of solute drag applied to austenite-ferrite phase transformations in Fe-C-Mn alloys

Delft University of Technology

A novel 3D mixed-mode multigrain model with efficient implementation of solute drag

applied to austenite-ferrite phase transformations in Fe-C-Mn alloys

Fang, H.; van der Zwaag, S.; van Dijk, N. H.

DOI

10.1016/j.actamat.2021.116897

Publication date

2021

Document Version

Final published version

Published in

Acta Materialia

Citation (APA)

Fang, H., van der Zwaag, S., & van Dijk, N. H. (2021). A novel 3D mixed-mode multigrain model with

efficient implementation of solute drag applied to austenite-ferrite phase transformations in Fe-C-Mn alloys.

Acta Materialia, 212, [116897]. https://doi.org/10.1016/j.actamat.2021.116897

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

Acta

Materialia

journalhomepage:www.elsevier.com/locate/actamat

Full

length

article

A

novel

3D

mixed-mode

multigrain

model

with

efficient

implementation

of

solute

drag

applied

to

austenite-ferrite

phase

transformations

in

Fe-C-Mn

alloys

H.

Fang

a,b,c,∗

_,

_S.

_van

_der

_Zwaag

b,d

_,

_N.H.

_van

_Dijk

a

a Fundamental Aspects of Materials and Energy group, Faculty of Applied Sciences, Delft University of Technology, Mekelweg 15, 2629 JB Delft, The Netherlands

b Novel Aerospace Materials group, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands c Department of Mechanical Engineering, Technical University of Denmark, DK-2800, Kgs. Lyngby, Denmark

d School of Materials Science and Engineering, Tsinghua University, Beijing 10 0 084, China

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 23 September 2020 Revised 4 March 2021 Accepted 10 April 2021 Available online 16 April 2021

Keywords:

Diffusional phase transformations Mixed-mode multigrain model Solute drag

Mn partitioning Low-alloyed steels

a

b

s

t

r

a

c

t

Acomputational3Dmodelthataccountsforbothnucleationandinterfacemigrationisaveryusefultool tomonitorand graspthecomplexity ofmicrostructureformationinlow-alloyed steels. Inthepresent studywehavedevelopeda3Dmixed-modemultigrainmodelfortheaustenite-ferriteandthe austenite-ferrite-austeniteformationcapableoffollowingdiffusionalphasetransformationsunderarbitrarythermal routes.Thisnewmodelincorporatesthesolutedrageffectofasubstitutionalelement(inthiscaseMn) andensuresanautomaticchangeintransformationdirectionwhenchangingfromheatingtocoolingand vice-versa.Ananalyticalsolutionforcalculatingtheenergydissipationofsolutedragtogetherwith mul-tipleregressionapproximationsforchemicalpotentials areproposed whichsignificantlyacceleratethe computation.ThemodellingresultsarefirstbenchmarkedforanFe-0.1C-0.5Mn(wt.%)alloyunder differ-entcontinuouscoolingandisothermalholdingconditions.Themodelrevealedrelativelylargevariations intransformationkineticsofindividualgrainsasaresult ofinteractionswithneighboringgrains.Then themodel isappliedtopredict thetransformationkineticsofaseriesofFe-C-Mnalloys duringcyclic partialphasetransformations.Thecomparisonwithexperimentaldilatometerresultsnicelyvalidatesthe predictionsofthismodelregardingthechangeinoveralltransformationkineticsoftheferrite transfor-mationasafunctionoftheMncontent.Newfeaturesofthismodelareitsefficientalgorithmtocompute energydissipationbysolutedrag,itscapabilitiesofpredictingthemicrostructuralstateforspatially re-solvedgrainsandtheminimalfinetuningofmodellingparameters.Thecodetoimplementthismodelis publiclyavailable.

© 2021TheAuthor(s).PublishedbyElsevierLtdonbehalfofActaMaterialiaInc. ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

The kineticsof theaustenite(

γ

) to ferrite (

α

) phase transfor-mation in steelshasbeen studiedfordecades usingboth experi-mentalandmodellingapproaches [1,2].Morerecently,an increas-ing interest in the kineticsofthe ferrite-to-austenite transforma-tion has beenobserved, motivated by thedesign and production of advanced Mn-rich high-strength steels [3-5]. It is well known that commonsubstitutionalalloyingelementssuchasMn,Ni and Mo strongly retard the interface migration during the

austenite-∗_{Corresponding author}

E-mail addresses: h.fang@tudelft.nl , hfang@mek.dtu.dk ,

haixingfang868@gmail.com (H. Fang).

ferrite phase transformation [6,7]. This retardation effect is re-flected most noticeably by experimental observations of a stag-nant austenite-to-ferritetransformationwhere theferrite fraction does not change for a notable period of time in Fe-C-Mn, Fe-C-Ni,Fe-C-MoandFe-C-Mn-Moalloysincaseofreversionofcooling intoheatingorvice-versa [8-14].Ascanbeseeninrecentstudies on cyclicpartial phase transformations, the stagnant stage exists in both austenite-to-ferrite and ferrite-to-austenite phase trans-formations [15,16]. The experimentally observed stagnant trans-formation is generally considered to be caused by partitioning of substitutional elements in the vicinity of the migrating inter-face. Due to the large difference indiffusivities of interstitial el-ements(e.g.C)and thatof substitutionalelements(e.g. Mn),the transformation kinetics during linear cooling, isothermal

anneal-https://doi.org/10.1016/j.actamat.2021.116897

1359-6454/© 2021 The Author(s). Published by Elsevier Ltd on behalf of Acta Materialia Inc. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

ing or more complex thermal routes generally lies in between the predictions for local equilibrium (LE) and para-equilibrium (PE). Understanding the effect of the slower diffusional substi-tutional elements on the transformation kinetics is essential to be ableto tune the microstructure. So far, it hasremained chal-lenging to predict and quantify the interface migration for an arbitrary imposed thermal profile involving both austenite and ferrite.

There aretwo main approachesestablishedto accountforthe effect of substitutional elements on the interface migration. The first approach uses the concept of an effective interface mobility that considers theeffectthe alloyingelement wouldhaveonthe characterofinterfacemobility,resultinginadecreaseinthevalue of the interface mobilitycompared to the value of intrinsic

γ

/

α

interfacemobility [17-20].Thisapproachisdemonstratedtobe ef-ficientin computations,butrequiresa goodestimationofthe in-fluence on the effectivemobility beforehand.As the effective in-terface mobility can vary with the steel compositions, the cool-ing rate and the transformation direction, the value of the pre-factor M0 in the Arrhenius expression of the effective interface

mobilitycan varyover severalordersof magnitudefrom10-10 _to

10−6 molmJ−1s−1 [21-24].Thiswiderangeofvaluesforthe effec-tive interface mobilitygenerally requiresrepeated trialsbefore it can beusedforaccurate predictions.Therefore,thegeneral appli-cation of thisapproach is limited. The alternative approach isto adoptthesolutedragtheorytoaccount fortheeffectofa substi-tutionalelement,whileusingtheintrinsicinterfacemobilityto ac-countfortheenergyconsumptionbytheinterfacefriction [25-30]. Inthesolutedragtheory,thesubstitutionalelementcansegregate at the interface, causingenergy dissipation by trans-diffusion in-side the interface. Thedriving force forthe phasetransformation is potentially consumed by both the solute trans-diffusioninside the interface and the interface friction. As the intrinsic

γ

/

α

in-terface mobilitycan be derived experimentally from themassive phase transformationinbinaryFe-X (X= Mn,Ni,Moetc.)alloys, it could be regarded as generic for austenite-ferritephase trans-formations,regardlessofthechangesinnominalcompositionsand cooling/heating rates.A recent studyon a seriesof Fe-Ni, Fe-Mn andFe-Coalloysprovidesagoodreferencefortheparametersthat determinetheintrinsicinterface mobility [31].Oneother parame-ter one needs toknowto calculatetheenergy dissipationdueto solute drag is the binding energy of a specific substitutional el-ement, which describes the affinity of the solute species to the austenite-ferrite interface. Although an accurate determinationof the binding energy is not straightforward, thisparameter can be estimatedfromdetailedmeasurementsofthe compositionprofile across the interface by atom probe tomography [32,33] or from first-principlescalculations [34].Thesolutedrageffectofdifferent alloying elementscan now be included andthus leadsto a cou-pled solute drageffect,whichhasbeen demonstratedtoperform well intherecentworkonFe-C-Mn-Moalloys [13].Therefore,the approachthatusesthesolutedragtheorytoaccountfortheeffect of substitutional elementson interface migration shows a better transferabilitythan theoneemployingtheeffectiveinterface mo-bility.

As the transformationkinetics isnot only controlled by inter-face migration,butalsoby nucleationofthenewphase,the evo-lution of the newand parent grain structure is a resultof both effects. Thus,the effectof alloyingelementsonthe interface mi-gration cannotsolelybe studiedbyconventionalcontinuous cool-ing or heating, or by isothermal holding experiments in which nucleation and growthof the ferrite phase occur simultaneously. An effectiveapproach tocircumventthisissueisto adoptthe so-calledcyclicpartialphasetransformation,wherethesteelis ther-mallycycledbetweentwo temperatures,that bothlieinthe two-phase

γ

-

α

region,sothat onlythemoreenergeticallyfavored

in-terfacemigrationtakesplace,ratherthanthesimultaneous occur-rence of interface migration and nucleation of new grains [11]. This behavior hasbeen observed by our recent neutron depolar-izationmeasurements,whichdemonstratethat newnucleationof ferritegrainsisindeedabsentandthatthetransformationkinetics is only controlled by the interface migration duringcycling [35]. Althoughall cyclingexperiments ofFe-C-Mnalloys(compositions with0.023-0.25wt.%Cand0.17-2.1wt.%Mn)showstagnantstages inthetransformationofaustenitetoferriteandviceversa,the cy-cling behavior is distinctly different. Forexample, in a lean Mn-alloyedFe-0.023C-0.17Mn(inwt.%)steelanopenloopoftheferrite fractionisobserved [11],whereasagradualnetincreaseofthe fer-ritefractionispresentinanFe-0.25C-2.1Mn(inwt.%)steelduring cycling [35].Anotherstudyonausteniteformationduringthermal cyclinginamedium-Mnsteel(Fe-0.2C-4.5Mninwt.%)showedthat thefractionofthenewlyformedaustenitedecreasesineach suc-cessivecycle,whilstthetotalamountofausteniteincreasesduring cycling [36].Thesedifferencesinbehaviorareduetothedynamics ofthemodificationsofthelocalchemicalcompositionsacrossthe interface,theelement chemicalpotentialsandtherelative veloci-tiesoftheMndiffusionandthemovinginterface.Oneofthebest waysto clarifythis complexityin dynamicsisto develop models thatrevealtheeffectsoflocalchemistry(includingCandMn)and phase structure at the same time under a mixed-modeinterface condition. To do this, in recent years extensive modelling stud-ies havebeen conductedto investigatethekinetics ofcyclic par-tial austenite-ferrite phase transformations using the solute drag theory [37-42].However,giventheheavycomputationaldemands afterincorporating solute drag, mostestablished models are lim-itedto1D,withtheexceptionofone2Dcellularautomatonmodel

[41], andall of these models only focus on the interface migra-tion itself. Althoughthese studies provide a detailedinsight into thetransformation behavior duringlinearcooling, isothermal an-nealingandthermalcycling, twokey parameters that needto be includedtomakethesteptomodeltransformationkineticsinreal (3D)steelsarestillmissing. Firstly,the1Dsimulationprovidesno spatialinformationand givesno insightinto thevariation in be-haviorofdifferent ferrite grains. Secondly,thevarying degreesof nucleation(i.e.degreeofundercooling)ofindividualferritegrains before entering thermalcycling hasbeen discarded,and asa re-sultthepreviousmodelscontainnoinformationonthegrainsize andthegrainsizedistribution.However,aspointedoutinour pre-vious study [43], theaverage grainsize andthe grainsize distri-bution are essential information for a model to predict the mi-crostructure and shedlight onthe underlying physics. Therefore, itisverydesirable todevelopa modelthatnot onlyincorporates solutedrag,butalsoincludesthespatialinformationin3Dto gen-erate a statistically relevant number of events for ferrite forma-tion.Amodelofthiskindwillalsomakeitpossibletostudyhow thesizedistributionofferritegrainsformedbeforecyclingaffects the interface migrating back and forth in the presence of solute drag.

Inthepresentwork,we haveextendedourprevious3D multi-grainmixed-modemodelbycouplingitwiththesolutedrag the-ory.Thisnewmodel accountsforferritenucleation basedon the classical nucleation theory and interface migration based on the balanceinGibbsfreeenergiesbetweenthechemicaldrivingforce and the energy dissipations due to interface friction and solute drag, all calculated locally for each moving interface. The mod-ellingresults are firstbenchmarkedfor an Fe-0.1C-0.5Mnternary alloyundercontinuouscooling andisothermalholding.Next,this model is applied to describe the transformation kinetics during cyclicpartialphasetransformations ofthesamesteel.Finally,the model is applied to calculate the cyclic transformation behav-ior asa function of C and Mn concentrations in the steels. The aim is to develop a versatile and flexible modelling tool to

pre-dict thekineticsofaustenite-ferritephase transformationswitha minimal numberoffittingparameters.The moreefficientmethod to compute the local energy dissipation due to solute drag pro-posed in the present work is also expected to provide inspira-tion for further 3D simulations of multigrain solid-state phase transformations.

2. Model

The startingstructurein themodelis afullyaustenitic multi-grain structure produced by constructing Voronoi cells in a cu-bic box with a length of L_b. The number density of austenite grains (

ρ

_γ) ispreset, resulting inan average austenitegrain size ¯d_{γ =}

₍

6_/

πρ

_γ

)

1/3_{. The centers of the Voronoi cells are randomly}

produced with an imposed minimum distance (dmin) to control

the austenite grain size distribution. The corners of the Voronoi cellsareregardedasthemainpotential nucleationsitesforferrite grains, asferrite nucleation is foundto predominantlytake place at grain corners. Besides the grain corners, also edges and faces canact aspotentialnucleationsiteswithalower probability [44]. Onceastableferritenucleusformsatoneofthegraincorners(or alternatively at a grain edge orat a grain face), the ferrite grain is assumed to grow isotropically into the austenite as a sphere. Asshowninthemicro-beamX-raydiffractionduring austenite-to-ferrite andferrite-to-austenite transformations [45,46], the nucle-ationbarrierisapparentlyhigherthanthekineticenergies. There-fore, it is reasonable to use classical nucleation theory (CNT) to compute thenumberdensityofferritegrains [47-49].The migra-tionoftheaustenite-ferriteinterfaceiscontrolledbythediffusion ofbothcarbonandsubstitutionalelementsandbyalattice recon-struction.Thisisaso-calledmixed-modeapproach [1,3,18]andthe relevantassumptions ofthemodelasusedaredescribedindetail later.Thesolutedragtheoryisappliedtodescribetheenergy dissi-pationduetothetrans-diffusionofsubstitutionalsoluteinsidethe interface. Aquasi-steady state is reached when the pressure im-posedonthemovinginterfaceiszero,indicatingthatthechemical drivingforceisbalancedbytheenergyconsumptionbyboth inter-facefriction andsolutedrag.The aboveprinciplehasbeenwidely adopted in various modelling approaches and currently sets the standard. However, the accuracy ofthe model depends on being abletocalculatetheconcentrationsofthekeyinterstitialand sub-stitutionalalloyingelementsatmovinginterfacesaswellasinthe parentphasestheinterfaceismovingtowards.

In a previous work we derived analytical expressions for the carbon(i.e.theinterstitial element)concentrationattheinterface andfaraway fromtheinterfacefornon-overlappingand overlap-ping of diffusion fields (soft impingement).The base modelalso accounted for the hard impingement of growing ferrite grains. Theseexpressionswereshowntobecorrectregardlessof transfor-mationdirections(i.e.ferritegrainsaregrowingoraustenitegrains aregrowing).

2.1. Ferritenucleation

In the simulations to be presented,ferrite nucleation was as-sumed to take placeonly at the corners of austenitegrains. The code hashowever been generalized such that it can handle dif-ferent nucleation sites, such as grain corners, grain edges and grain faces or mixtures of these. Dedicated simulations showed thatthechoiceofthedominantnucleationsitehadamarginal ef-fect on thetransformation and doesnot to affect the main con-clusions ofthe presentwork(minor changes were onlyobserved inthe laterstagesofthetransformation).The classicalnucleation theory is used to describe the nucleation rate dN/dt in a se-lected volume. According to the CNT, the nucleation rate can be

writtenas: dN dt =AZN0

(

1− fα

(

T

)

)



kBT h



exp



−QD kBT



exp



−



kBT

(



GV

(

T

)

−



GS

)

2



, (1)

whereAisa pre-factor,Zis theZeldovich factorandnearly con-stant(Z_{≈ 0.05),}N0representsthenumberofpotentialnucleation

siteswhentransformation starts,f_α isthevolume fractionof fer-rite phase, kB is the Bolzmann constant, T is the temperaturein

Kelvin,hisPlanck’sconstant,QDistheenergybarrierfordiffusion,



isaconstantthatcomprisesallthecontributionsoftheshapeof thecriticalnucleusandinterfacialenergybetweenthenucleusand thesurroundingparentgrains,



GV isthenetdifferenceinGibbs

freeenergyperunitvolumebetweenferriteandaustenite,



GSis

themisfitstrainenergyperunitvolume.Thecriticalnucleussizeis r∗=2

σ

_αγ/

(



GV

(

T

)

−



GS

)

[49],where

σ

αγ isthe interfacial

en-ergyperunitareafor

γ

/

α

boundariesandamountsto

σ

_αγ =0.62 Jm−2[50].In Eq.1theconstant



isthemostchallenging parame-tertobedeterminedprecisely.However,highenergyX-ray diffrac-tionstudiesmonitoring thenucleationofferriteoraustenite sug-gestthat



inrealityhasavalueoftheorderof10-8_J3_m−6 _[45].

Therefore,a value of



≈ 5× 10-8 _J3_m−6 _{isused inthepresent}

work.The



GSiscalculatedby



Gs= ₉2(₍1₁₊−ν)_ν)

μ

(

2a3_α/a3_γ− 1

)

2[51],

where

ν

isthe Poisson’sratio,

μ

isthe shearmodulus(

ν

= 0.33 and

μ

= 60 GPa for pure iron). The lattice parameter of ferrite a_α andthelatticeparameterofaustenitea_γ bothdependon tem-peratureandcarbonconcentration,andarecalculatedaccordingto

[52].

2.2.

γ

/

α

interfacemigration

The chemical drivingforce per mole ofatoms



Gchem

m forthe

γ

/

α

interface migration after ferrite nucleation can be expressed as:



Gchem m = n  i x0 i



μ

γ α i



xγ α_i

−

μ

αγ i



xαγ_i

, (2)

wheresubscriptiistheelementinthealloy,nisthetotalnumber ofelements,x0

i isthecomposition ofelementitransferred across

the interface, the superscripts of

γ α

and

αγ

correspond to the austeniteside andferrite side on the interface, respectively,

μ

is the chemical potential and x is the mole fraction. For an inter-face moving in aquasi-steady state the chemical driving force is balanced bytheinterface friction



Gfriction

m andthedissipation of

substitutionalelementsinsidetheinterface



Gdiff m :



Gchem

m =



Gfrictionm +



Gdiffm . (3)

Suchan approachtosplitthedissipationofthechemical driv-ingforce intotwocontributionswasfirstproposed byHillertand coworkers [25,53,54] and is now a common approach (e.g. refs.

[13,30,37]). Eq.3alsoindicates thatlong-rangediffusionofa sub-stitutionalelement(Mn inthiscase)inferrite orausteniteis not considered, although it may contribute to the Gibbs free energy dissipation. According to [55], it is not evident how to separate the energydissipation dueto trans-diffusion ofsubstitutional el-ementinside theinterfacefromthedissipation bythelong-range diffusion.Sinceinterstitialatomsdiffusemuchfasterthan substitu-tionalatoms,itisreasonabletoassumethatthereisalong-range diffusionoftheinterstitial element (Cinthiscase)anda steady-state conditionofthe substitutionalelement (Mnin thiscase)in andacrossthethininterfaceregion [56],asisdoneinthecurrent model.

Theinterfacefrictiondependsontheinterfacevelocityv_int and theintrinsicmobilityofthe

γ

/

α

interfaceMint:



Gfrictionm =vintVm/Mint, (4)

where Vm is the molar volume of Fe and assumed to be the

same for austenite and ferrite, Mint is expressed by Mint=

M_0,intexp

(

_−Qint/RT

)

whereM_0,int is a pre-factor in molmJ−1s−1, Qint istheactivationenergyandR isthegasconstant.Theuseof

the intrinsicinterface mobility,rather thanthe effectiveinterface mobility,provides asubstantialadvantage intransferabilityofthe modelfromonesteelcomposition toanother.Thedissipation en-ergy duetotrans-diffusionofsolute insidetheinterface



Gdiff

m is

calculated usingthe methodproposed byPurdy andBréchet [57]. As



Gdiff

m dependsonvint,whichisrelatedtotheinterfacial

com-position,the conventionalwaytoderive thesolutionofvint isby

calculating the values of the three energies (



Gchem

m ,



Gdiffm and



Gfriction

m ) as a function of vint inthe expected range. The

inter-sectionpointbetween



Gchem

m and



Gdiffm +



Gfrictionm isregardedas

thesolutionofvint (whenmultipleintersectionpointsarepresent,

the minimumv_int is adopted).Although thismethodis robust,it is computationally heavy asit needsto be calculated via numer-ical procedures for every moving interface while takinginto ac-countthelocal(thermalhistorydependent)chemicalcomposition. In this work, we have derived an analytical solution for



Gdiff

m ,

which avoids the need to compute the Gibbs free energies as a function of v_int and allows one to directly compute the solution. Thereby,thecomputationcanbesignificantlyacceleratedenabling themodeltobeusedinamultigrainsetting.Adetaileddescription onhowtocalculate



Gdiff

m ispresentedin AppendixA.

The interfacial carbon concentrations and diffusion profiles need to be solved to derive vint. As the diffusivity of carbon in

ferrite Dα_C ismuchlargerthan thatinausteniteDγ_C,we only con-siderthebulkdiffusionofcarboninaustenite.Thisway,thereisno carbonconcentrationgradientinferriteandtherebytheinterfacial carbonconcentrationCα amountstotheequilibriumconcentration Cαγ_eq.Thediffusionprofileofcarboninausteniteisapproximatedby a second-orderpolynomial.Thisapproximationyields mathemati-cal simplicityandwasfoundtogiveagoodmatchwiththe limit-ing diffusionalcasesforwhichexactsolutionsexist.Detailsofthe mathematical treatment of the carbon diffusion are presentedin

Appendix B.Fornon-overlappingdiffusionfieldstheruleofmass conservationofCresultsin:



C0− Cαγeq

V_α =4

π

 L 0

(

C

(

r

)

− C0

)

(

r+Rα

)

2dr =−2

π

(

C0− Cγ

)



L3_+5L2_R α+10LR2 α

15 , (5)

whereC0 isthenominalcarbonconcentration,Cαγeq isthe

equilib-riumconcentrationofcarboninferrite,R_α istheradiusandV_α is thevolumeoftheferritegrain,Listhediffusionlengthandristhe distance,Cγ istheinterfacialcarbonconcentrationattheaustenite side.Asthereisnocarbonaccumulatedattheinterface itself,the followingequationcanbederived:

vint



Cγ − Cαγeq

=2DγC L

(

Cγ− C0

)

. (6)

Nowwecanderivethesolutionsfortheunknownparametersvint,

Cγ and the diffusion length L by solving the set of Eqs. 2-6. In

Section2.3wepresenthowwesolvethisproblem.

Fortheoverlappingdiffusionfields(soft-impingement)the car-bonconcentrationatthesoftimpingementpointCγmincreases.The

massconservation ofcarbon andtheabsence ofaccumulation of carbon at the interface still applies in this stage. Therefore, we have:



C0− Cαγeq

V_α =4

π

 L 0



(

Cm− C0

)

+

(

Cγ− Cm

)



1−r L

2

(

r+R_α

)

2dr =−2

π

L



9L2_C m− 10L2C0+L2Cγ +20R2_αCm− 30R2_αC0+10R2_αCγ+25LRαCm− 30LRαC0+5LRαCγ

15 , (7) vint



Cγ− Cαγeq

= 2DγC L

(

Cγ − Cm

)

, (8) L=L0− Rα, (9)

whereL0isthedistancewherethecarbondiffusionfieldofa

fer-ritegrainstartstocontactwiththediffusionfieldsofitsneighbors. Nowthesolutionsforvint,Cγ,Cm andLcanbederivedbysolving

thesetof Eqs.2-4and Eqs.7-9.

Todistinguishwhethersoftimpingement occurs,we firsttreat thediffusion ofcarbonwithout softimpingement andderive the diffusion length L.When theprofile of grain jstarts to intersect with that of any neighboring grain k, the following condition is fulfilled:

R_α,j+Lj+Rα,k+Lk=rjk

(

j=k

)

, (10)

whererjkisthedistancebetweenthecentersofferritegrainsjand

k.AtthismomentthediffusionlengthListakenasthemaximum distancethat carboncan diffuseforthat grain.Withthe advance ofthe interface migration, the ferrite grainscould be in physical contactwithneighboringgrains,resultinginhardimpingement.To accountforthishardimpingementeffect,weusethesame correc-tionasthe one appliedin ourprevious work [43],for whichthe sourcecodenowcanbefoundat(https://github.com/haixingfang/ 3D-mixed-mode-model).

2.3. Simulationconditions

The model is used to compute transformation kinetics under various thermalconditions,includingcontinuous cooling, isother-mal holding and thermal cycling for a series of Fe-C-Mn alloys. Thismodelisfirstbenchmarkedforcontinuouscoolingatdifferent ratesandisothermalholdingataselectedtemperatureforthe Fe-0.1C-0.5Mnalloy(seetemperatureprofilein Fig.1aand 1b). Then wefocus onusingthismodelforcyclicpartial phase transforma-tions inthe

γ

/

α

two-phase regionbetween temperaturesT1 and

T2 (see Fig. 1c). Table 1 givesthese temperaturesand the

corre-spondingequilibriumferritefractionscalculatedwithThermo-Calc softwareusingtheTCFE8database. Thetwo cyclingtemperatures werechosentobethesameasinthereportedexperimental stud-ies (for Fe-0.023C-0.17Mn [11] and Fe-0.1C-0.5Mn alloys [38]) or chosentoensurethat theferritefractionisnottoosmallandalso nottoolargeattheendofisothermalholding(fortheotheralloys examined here). For those alloys the cycling temperatures were chosentoshowortho-equilibriumferritefractionsoffOE

α (T1)≈ 0.70

andfOE

α (T2) ≈ 0.20. For alloyswith 1.5wt.% Mn or morethe T1

wassetasf_α(T1) ≈ 0.50for theendofthe isothermal stageand

T2 fulfills T2 - T1 ≈ 50 K. The reasonfor setting up a different

criterionforhigherMnalloysisthatasubstantialferrite transfor-mationonlyoccurswellbelowthetransitiontemperaturebetween

Fig. 1. Schematic temperature profiles for (a) continuous cooling, (b) isothermal and (c) cyclic partial phase transformations. All simulations start at A 3 (PE) and end at A 1e

for continuous cooling and cyclic transformations, while the simulation does not end until there is negligible change of f αfor the isothermal transformation. For simulations

of the cyclic phase transformation, the primary cooling rate is set as 20 K/s, the isothermal time at T 1 is set as 30 min and the number of cycles is chosen to be 3. Table 1

Compositions of the Fe-C-Mn alloys and the cycling temperatures T 1 and T 2 for each simulation. The ferrite fractions under

ortho-equilibrium f OE

α, the A 3 and A 1 temperatures under ortho-equilibrium (A 3e and A 1e , respectively), para-equilibrium

temperature A 3 (PE) and PLE/NPLE temperature are also listed.

Alloy (wt.%) T 1 (K) T 2 (K) fαOE(T 1 ) (-) fαOE(T 2 ) (-) A 3e (K) A 3 (PE) (K) A 1e (K) PLE/NPLE (K)

Fe-0.1C-0.17Mn 1073 1127 0.70 0.20 1136 1132 994 1131 Fe-0.1C-0.5Mn 1058 1115 0.71 0.20 1124 1118 983 1108 Fe-0.1C-1.0Mn 1040 1095 0.70 0.24 1107 1096 966 1074 Fe-0.1C-1.5Mn 993 1043 0.86 0.71 1092 1073 947 1039 Fe-0.1C-2.0Mn 953 1003 0.91 0.82 1077 1051 926 995 Fe-0.1C-2.5Mn 938 983 0.92 0.85 1063 1033 903 988 Fe-0.25C-0.17Mn 998 1075 0.68 0.20 1091 1089 1091 995 Fe-0.25C-2.1Mn 939 983 0.77 0.63 1045 1019 1045 938 Fe-0.023C-0.17Mn 1133 1158 0.80 0.42 1166 1163 1166 994 Fe-0.05C-2.0Mn 998 1062 0.82 0.45 1090 1064 1090 920

partitioninglocalequilibriumandnegligibleportioninglocal equi-librium(PLE/NPLE)accordingtopreliminarysimulationresults.The PLE/NPLEtemperaturesarealsolistedin Table1.

The length of the cubic sample box is set as Lb = 70

μ

m.

The average austenite grain size is chosen to be d¯_{γ =} 20

μ

m (

ρ

_γ = 2.4× 1014 _m−3_{) withdmin = 12}

_μ

_m.Giventhese

dimen-sions andtwoimposed coolingratesthemodeldescribesthe col-lectiveandindividualbehaviorofabout151parentaustenitegrains andabout58-204emergingferrite grainsdependingonthe cool-ingrate.Asthechemicalpotentialsofallelementsdependonboth temperatureandcompositions,theymustbecalculatedseparately for each grain at each time step. To accelerate the computation, we usemultiple nonlinear regression to approximatethe depen-dency of the chemical potentials of Fe, C andMn beforehandto avoidcalls to theThermo-Calcdatabase. Adetaileddescriptionis presentedin AppendixB.Usingthisapproximation,combinedwith the analytical expression derived for



Gdi f fm (shown in Eq.A4 in

Appendix C), the calculation ofthe transformationkineticsin 3D foramultigrainsettingcanbesignificantlyacceleratedandbe re-alizedinrealisticcomputingtimes.

The phaseboundarylinesof(

α

₊

γ

)/

γ

and

α

/(

α

₊

γ

) are calcu-lated with Thermo-Calc and fitted to a second-order polynomial forthetemperaturerangeofinterest.Theparametersforintrinsic interfacemobilityproposedbyZhuandcoworkers [31]areusedin the present work. The diffusivityof Mn in the interface is taken asDMn

int =



Dγ_MnDα_Mn.Except forthethermodynamic data,the val-ues ofall other parameters arekept constant independentofthe alloycomposition andthermaltreatment.The modelling parame-tersarelistedin Table2.ThebindingenergyE0 forMnonthe

α

/

γ

interfaceisaveryimportantparameterasitgovernstheMn parti-tioninginsidetheinterface,therebystronglyinfluencingthe trans-formation kinetics. Slightlydifferent values forE0 havebeen

de-rivedusingdifferentapproaches:scanningAugermicroprobe stud-iesonaustenitegrainboundariesyieldedE0=8± 3kJmol−1[58];

fittings to experimental transformation kineticsyielded E0 = 5 ~

9.9 kJmol−1 [37]; calculationsfrom theMn profileacross the

in-terface measured by atom probe tomography gave E0 = 6.0 ±

1.4kJmol-1 _{[32] andfirst principlescalculationsprovided avalue}

of E0 ≈ 12.5 kJmol−1 [34]. In the presentwork we havechosen

an intermediatevalue ofE0 =7kJmol−1.The choiceofaslightly

differentvalue forE0 wouldhave affectedto overallkinetics,but

wouldnothaveaffectedthekeyfindingsinthesimulations. In this model the only adjustable parameter is the pre-factor Ain Eq.1.Thispre-factorAaffectsthemaximumnumberof nu-clei and the effectivetime (temperature) window for nucleation. The parameterdepends onthe rateof transformation,which de-pends onthe imposedtemperature path.It shouldbe noted that any small uncertainty in the modelled energy barrier for nucle-ation (controlledby the parameter



) willdirectly translate into a significant change inthe value of the required pre-factorA. In the currentsimulations forlinear cooling, the value for the pre-factorAisestimatedfromthemaximumnumberdensityofferrite grains,whichwasestimatedfrommetallographicobservations.For isothermalholdingandthermalcyclingsimulations,thepre-factor Aissetintherangeof0.001 ~ 0.0033(i.e.1/1000~ 1/300).This choicewasbasedonreportedvaluesinsynchrotronX-ray diffrac-tionstudiesonisothermalferrite-to-austenitetransformations[46]. Foreachtimestep



tduringthecalculationofthephase trans-formation thenumber offerrite grainsis calculatedaccording to theCNT (see equation 1).The interface velocity foreach grain is then calculated as described in Section 2.2. The effective radius R_α,j (aftercorrectionsofhardimpingementifnecessary)offerrite grainjatatimetiscalculatedusingthevelocityderivedfromthe lasttimestep

v

int, j

(

t−



t

)

:

R_α,j

(

t

)

=Rα,j

(

t−



t

)

+

v

int, j

(

t−



t

)



t. (11) Theoverallmicrostructuralcharacteristics,includingtheferrite volume fractionf_α, theaverage grain radius

δ

_α andthe standard deviation

σ

_α fortheradiusdistributionofferritegrains,are com-putedforeachtimestep.Theprogrammingcodesfor implement-ingthemodelpresentedherearepubliclyavailable(https://github. com/haixingfang/3D-GEB-mixed-mode-model).

Table 2

Modelling parameters and the data references.

Parameter Value or equation Unit Ref.

Density of potential nucleation site, ρ0 1 . 14 × 10 15 m −3 -

Pre-factor A in Eq. 1 1/1000 ~ 1/300 for simulations on cyclic phase transformations - -

Activation energy for Fe self diffusion, Q D 3 . 93 × 10 −19 J [45]

Molar volume of Fe, V m 7 . 1 × 10 −6 m 3 mol −1 -

Interface thickness, δint 0 . 5 × 10 −9 m [37]

Intrinsic mobility of γ/ αinterface, M int 2 . 7 × 10 −6 exp ( −145kJmol

−1

RT ) m 4 J −1 s −1 [31]

Binding energy of Mn, E 0 7 kJmol −1 -

Diffusivity of Mn in austenite, D γMn 0 . 178 × 10 −4 exp ( −264kJmol −1 RT ) m 2 s −1 [59] Diffusivity of Mn in ferrite, D α Mn 0 . 756 × 10 −4 exp ( − 224.5kJmol−1 RT ) m 2 s −1 [59]

Diffusivity of C in austenite, D Cγ 4 . 53 × 10 −7( 1 + y C( 1 − y C)8339T.9K) exp { −(1TK − 2 . 221 × 10 −4)( 17767 − 26436 y C)} where y C = 1−xxCC m

2 s −1 _[60] Diffusivity of C in ferrite, D α C 0 . 02 × 10 −4 exp ( − 10115K T ) exp { 0 . 5898[ 1 + 2 πarctan( 14 . 985 −15309T K) ] } m 2 s −1 [61]

Fig. 2. Kinetics of austenite-to-ferrite phase transformation in Fe-0.1C-0.5Mn alloy during continuous cooling at constant rates of 0.4 and 10 K/s predicted by the present model. (a) Nucleated ferrite number density ρα, (b) ferrite fraction f α, (c) average ferrite grain size δαand standard deviation σα, and (d) final grain size distribution of

ferrite for 10 K/s in comparison with metallographic measured results from [19] . Lines in (d) are fitting curves of lognormal distribution.

3. Resultsanddiscussion

3.1. Transformationsduringcontinuouscooling

Fig.2 showsthe ferritetransformationkineticsinthe Fe-0.1C-0.5Mnalloyduringcontinuouscoolingatarateof0.4and10K/s. Inthesimulation,thepre-factorAin Eq.1wassetatA=1/3000 andatA=1/35for0.4and10K/s,respectively.The pre-factoris the onlyparameterrequiredto beadjusted anditischosen such tomakesurethatthemaximumvaluefortheferritenumber den-sitywascomparabletothatofmetallographicobservationsonthe fully transformed sample [62]. As a result, the maximum ferrite numberdensityfor10K/s isabout3.5timeshigherthanthat for 0.4K/s,eventhough thenucleationstartsataboutthesame tem-perature of 1095 K,as shown in Fig. 2a. Fig. 2b shows that the model predictedferrite fractions (f_α) that are in goodagreement withthedilatometerresultsforthetwocoolingrates. Fig.2cplots the evolutionofthecalculated averageferrite grainsize (

δ

_α) and

standard deviationofthe ferrite grain size distribution(

σ

_α) asa functionoftemperature.Thefinal valuesfor

δ

_α (10.7and6.9

μ

m for0.4and10K/s,respectively)areasintendedclosetothe exper-imental data determined on metallographicimages (10.7 and6.5

μ

m,respectively) [61]. Fig.2d showsa comparison ofthe ferrite grainsizedistributionafteracompletetransformationbetweenthe modelpredictionandtheexperimentalresultsfor10K/s.The orig-inalexperimentaldatafortheferritegrainsizewere presentedas theequivalentcirculardiameterin2D [19].Tomakeacomparison withour3Ddata,weconvertedthepublished2Ddatainto3D by R_α,3D= 4Rα,2D/

π

for each ferrite grain. The figure showsa

rea-sonable agreementbetweenthe modellingresultsandthe exper-imentaldata,both ofwhichcanbe approximatedby alognormal distribution.Itcanalsobeseenin Fig.2dthatthemodelpredictsa narrowersizedistributionandaslightlylargeraverage grainsize, whichcouldpartlybeattributedtoanearliertransformationonset forthemodelling (see Fig. 2b).It isalso worth pointingout that theactualgrain sizeandshapedistributionofthefullyaustenitic

Fig. 3. Visualization of the formed grains in parent austenitic structure at constant cooling rates of (a, b) 0.4 K/s and (e, f) 10 K/s for two different f α, and developments of two grains (marked in αA and αB ) as a function of temperature at cooling rates of (c, d) 0.4 K/s and (g, h) 10 K/s. (a, e) f α= 0.10, (b, f) f α= 0.25, (c, g) ferrite radius

and interfacial migration rate, and (d, h) chemical driving force Gchem

m , interfacial friction energy G f riction

m and energy dissipation due to trans-interfacial diffusion of Mn,

Gdi f fm . In (a, b, e, f) spherical surfaces of ferrite grains are shown in red and their cut-off planes on the edge of the cubic box are shown blue. Austenitic grains are made

semi-transparent and their grain boundaries are shown in black lines. In each of (c, g) five circle points are shown to mark the radius of the ferrite grain αA when f α= 0.01,

0.05, 0.10, 0.25 and 0.50. In (d, h) the temperatures for starting growth of αA and αB are marked as T A and T B , respectively.

startingstructureinfluencesthefinalferritegrainsizedistribution. Givenarelativelyuniformandisotropicstartingmicrostructure(as presented here), the final grain structure is mainly controlled by theaveragegrainsizeoftheaustenite.Themodelitselfcaneasily beadaptedtoincludenon-equiaxedstartingmicrostructures.

Wewillnowtryanddemonstratethecapabilitiesofthemodel inpredictingthetransformationkineticsbothattheoverallsample level andatthe levelof individual emergingferrite grains. Fig. 3

showsthedevelopmentoftwoselectedferritegrains(

α

A and

α

B)

nucleated at the same position but under differentcooling rates of 0.4 and 10 K/s, respectively. Snapshots of Fig. 3a-b and 3e-f show their positionsandtheoverall microstructurechanges from f_α = 0.10to0.25. Fig.3cand 3gshow theevolutions oftheir ra-diusandinterfacevelocities.Toanalyzethecontrollingmechanism for thegrowthof

α

A and

α

B, we plotthechanges ofthe

chemi-cal driving force(



Gchem

m ) andtheenergydissipations(



Gmf riction and



Gdi f f_m )in Fig.3dand 3h.Forcoolingat0.4K/s,

α

Astartsto

grow at1095K,abovewhich



Gchem

m <



Gdi f fm andtherefore,the interface cannotmove. Below1095K,theinterface velocityof

α

A

increases,peaksat1072Kandfallsbacktonearlyzeroas



Gdi f fm becomes very close to



Gchem

m .

α

B startsto grow later than

α

A.

Butonceitstartsat1081K,ithasaveryhighinitialvelocitydue to a large undercooling, witha larger driving force compared to

α

Aatthestartingpoint.Differentfrom

α

A,theinterfacevelocities

of

α

B continuously decreasetozerowitha smallbump occurring

at 1072 K.Forcooling ata rate of10 K/s,

α

A and

α

B show very

similar behaviorinterms ofthestartinggrowthtemperature,the energychangesandtheinterfacialvelocitiesduringthewhole pro-cess,whichisverydifferentfromthebehaviorobservedata cool-ingrateof0.4K/s.

The averageandlocalconcentration profiles ofC andMnalso evolve distinctlydifferentlyforthetwograinsatdifferentcooling rates.Itcanbeseenfrom Fig.4bthatMnconcentrationatthe in-terfacebuildsuptoahighlevel,especiallywhenf_α =0.50, while theCdiffusiongradientissmallandevensmearedoutatf_α =0.50 for

α

A upon cooling at 0.4 K/s (see Fig. 4a). This indicates that

trans-diffusionofMn in theinterface dominates thekinetics.For thesame grainnucleated atthesame temperatureandthesame positionbutexposed to a highercoolingrate, the Mn concentra-tionacrosstheinterface ismuchlowerandtheCdiffusion gradi-entismuchhigher(see Fig.4cand 4d),indicatingthatinthiscase (acoolingrateof10K/s)thetransformationkineticsismainly con-trolledbyCdiffusion.

Themodellingresultsofthetransformationkineticsduring con-tinuouscooling demonstratethat: 1) the modelperforms well in predicting the overall transformation kinetics during continuous cooling; 2) the development of physical and chemical character-isticscanbemonitoredforeachindividualgrainand3)the contri-butionanddissipationsoftheGibbsfreeenergiescanbecalculated tounderstandthemechanismforthechangeininterfacialvelocity.

3.2. Isothermaltransformations

A further modelling test wasmade on the same alloy of Fe-0.1C-0.5Mn but imposed to isothermal phase transformations at 1058 and 1048 K to enable model validation against phase field simulations and dilatometry measurements, respectively. In both simulations thepre-factor in Eq. 1wasset atA = 1/800. Fig. 5a showsthepredictedf_α asa function ofisothermal timetogether withthe resultobtainedfrom1D phasefield simulation.In both models,Mnismodelledsuch thattrans-diffusionintheinterface, i.e.solutedrag, cantake place.Thefigure showsthatinboth ap-proaches f_α evolves in a similar manner and ends at about the same level of 0.65, which is smaller than the para-equilibrium fraction fPE

α (T= 1058 K) = 0.72. Ashas beenreported in

litera-ture[e.g.12, 13,37, 63],thisis a caseofincomplete transforma-tionbecausetheenergydissipationofsolute dragbecomeslarger thanthedrivingforceatlongerisothermaltimes.Notably,thereis an obviouslatenttransformationintheearlystage inthecurrent modellingresults(markedbytheredarrowsin Fig.5a),whilstthis effecthasnotbeenobservedforthe1Dphasefieldmodelling [39]. Thisis dueto differencesin nucleation modesandgrain geome-tries:thepresentmodelimplementsacontinuousnucleationmode

Fig. 4. Carbon (a, c) and manganese (b, d) profiles across the interface for grain αA (as visualized in Fig. 3 ) at cooling rates of (a, b) 0.4 K/s and (c, d) 10 K/s. r / δis a

dimensionless quantity where r is distance and δis the half of the interface thickness.

Fig. 5. Ferrite fraction f αas a function of time for isothermal holding of Fe-0.1C-0.5Mn (a) at 785 °C in comparison with the simulation result of the 1D phase field modelling [39] considering the Gibbs free energy dissipation due to Mn diffusion inside the interface and (b) at 775 °C against experimentally derived f αfrom dilatometer measurements

with a primary cooling rate of 1 K/s. The time is expressed relative to the starting moment of the isothermal holding. The experimental f αin (b) is derived from average of

two separate dilatometry measurements (each starting with a fresh sample).

and operatesin 3D,leading toa scaling off_α withR_α3_,whereas

an instantaneous nucleation and a planar geometry are assumed inthe phasefield andtherefore,f_α scales withR_α.Differences in themedianisothermalstage,wheref_α =0.25~ 0.55,arealsoseen, whichcould beduetothefact thatourmodeltakesintoaccount both nucleationandspatialcorrelationsbetweenneighboring fer-ritegrains.

Althoughthedifferencesintheoverallf_α betweenthepresent 3D multigrain model and the 1D phase field model predictions presentedin Fig.5aarenotlarge,the3Dmultigrainmodelshows asignificantvariationintransformationbehaviorforindividual fer-rite grains due to grain interactions (related to nucleation time

andlocalenvironment), whereasthe 1D modelbydefaultcannot give such information. Other experiments andmodelling studies

[45,64]alsodemonstratedthat thetransformationbehaviorof in-dividualgrainscanvarysignificantly.Thereareanumberof signif-icantdifferencesbetween1D and3Dmodels.Astothe3Dmodel, the transformation behavior for each grain is determined by the nucleationandgrowthbehaviorofthesurroundingmaterial,in ad-ditiontowhathasbeenprescribed forthatparticulargrain. Such behaviorcanneverbecapturedaprioriina1Dmodel.Onlywhen oneknowsthefull3Dbehaviorthenthesettingsofthe1D simu-lationcan besetproperly. Thedimensionalitywill alsoaffectthe waythat soft andhard impingement play a role inthe

transfor-Fig. 6. Modelling predictions and experimental dilatometry results [11,38] determined for the ferrite fraction as a function of temperature during cyclic partial transformation in (a) Fe-0.023C-0.17Mn and (b) Fe-0.1C-0.5Mn alloys. S: stagnant transformation where transformation kinetics is negligible; N: transformation proceeds in the direction in accordance with the temperature change; I: transformation proceeds in an opposite direction as would be expected from the temperature change.

mation kinetics.Thede-activationofnucleationsites bysoft/hard impingement andthe evolutionof thegrain size distributioncan onlybeobtainedina3Dmodel.Giventheirgeometricalsimplicity, 1D modelsare moresuitabletodemonstratetheeffectofvarious assumptions regardingsolute behavior near the moving interface onthetransformationkinetics.

Fig.5bcompares thepredictedf_α withtheexperimentally de-terminedf_α fromdilatometrymeasurementsbyapplyingthelever rule [65]. By setting the average austenite grain size to < d_γ

> =50

μ

minaccordancewithexperiments [66]andthesample box Lb = 175

μ

m (using thesame ratioofLb to< dγ >), while

keepingthe pre-factorAthe sameasinthesimulation shownin

Fig.5a, wefound thatthepredictedf_α matches theexperimental resultcloselyasshownin Fig.5b.Themodelandtheexperiment showalatenttransformationintheearlystageandasimilarvalue off_α atthefinalstagethatisclosetothepara-equilibriumfraction

fPE

α (T=1048K)=0.76).Thegoodagreementbetweenthe

simula-tions (usingonlyonefreely adjustableparameterand11 indepen-dentlydeterminedparameters)andthedilatometermeasurements validatesthegoodtransferabilityofthecurrentmodel.

It is worth to note that the present model is aimed to de-scribebulkbehaviorandthereforeassumesperiodicboundary con-ditions.Inthecurrentformat,themodelisnotappropriateto de-scribeandpredictthetransformationkineticsforisothermal trans-formationexperimentsthatarebasedontheremovalofcarbonvia gas/samplereactions atthesurface(includingboth masslossand macroscopiccarbondiffusiontowardstheoutersurfaceofthe ma-terial),e.g.thedecarburizationexperimentsasconductedin [14].

3.3. Cyclicpartialphasetransformations

When analyzing a cyclic partial phase transformation the austenite-to-ferritetransformationisclassifiedasleadingtoa pos-itive interfacialmigration (vint > 0)while the ferrite-to-austenite

transformationleadstoanegativeinterfacialmigration (v_{int <}0). Similarly, thedrivingforce fortheaustenite-to-ferrite transforma-tion is connected to positive changes in the Gibbs free energy whiletheferrite-to-austenitetransformationisconnectedto nega-tivechangesintheGibbfreeenergy.

3.3.1. Comparisonwithdilatometerexperiments

Wefirstverifyoursimulationresultsagainstdilatometer exper-imentsfortwodifferentalloysduringcyclingatarateof10K/min.

Fig.6 showsthat ourmodelcanadequately reproducethe exper-imental cyclic behavior of f_α forboth Fe-0.023C-0.17Mnand Fe-0.1C-0.5Mnalloys.Thestagnantandinversetransformationstages,

whichare two distinct features of thecyclicbehavior in lowMn steels,arecapturedinthesimulations. Whilstitcan beseen that oursimulations predicta longerstagnant stage(especiallyduring cooling) and a shorter inverse transformation stage, the simula-tions predict quite a good range of f_α where the transformation could span during each cycle and how f_α evolves in each cycle comparedtotheexperimentaldata.

Interestingly,inbothalloysf_αstartsatahigherlevelatthe be-ginningofthefirstcycle, whilstf_α neverreturnscompletelytoits original level in the following cycles, which themselves are per-fectly reproducible. This suggests that the transformation at the endofisothermalholdinghasreachedanequilibriumstage,which is not reached again during cycling due to the prevailing non-equilibriuminterfacialconditionscausedbydynamical changesin temperature.Itshouldbenotedthatthesetwotransformation fea-tures not always occur. In our previous experimental study on a steelwithahigherMnconcentration (Fe-0.25C-2.1Mn alloy [35]), f_α was found to be far from equilibrium at the end of isother-mal holding. Consequently, f_α progressively increased after each thermalcycle,showinganoppositebehaviorfromwhatisseenin

Fig.6.

3.3.2. Cyclicbehaviorofindividualgrains

Fig. 7 shows the behavior of the first nucleated grain during cyclingfromthesimulations ofFe-0.1C-0.5Mn alloy.Similar tof_α, R_α alsoexperiencesaloopconsistingoftwostagnant,twonormal andtwo inverse stages duringeach cycleshown in Fig.7a. Each stageischaracterizedbydifferentfeaturesoftheinterfacial veloc-itiesasplottedin Fig.7bandthefeaturescanbeunderstoodfrom thedifferencebetweenthechemicaldrivingforce



Gchem

m andthe energydissipation ofsolute drag



Gdi_mf f (see Fig.7c-e). Toassist theanalysis, wemarked pointsof interestby A-G ineach graph. Thecharacterofthetransformationstage markedby thesepoints issummarizedasfollows:

1) Stagnant stage (frompoint A to C). The interface hasreached the stasis state at the beginning of the cycling (point A) be-cause



Gchem

m <



Gdi f fm .Whentemperatureincreases,



Gchemm decreasesandreacheszeroatpointB.Afterthat,both



Gchem

m and



Gdi f f_m becomenegative,indicatingatendencyfor

α

_→

γ

. However,

α

→

γ

cannot occur because

|



Gchem

m

|

<

|



G di f f

m

|

,

untilpointCwhere



Gchem

m =



Gdi f fm .

2) Normal stage where

α

→

γ

proceeds during heating (from pointCtoD).AnincreaseintemperatureabovepointCmakes thedriving force sufficientlylarger thanthe solute dragforce,

Fig. 7. (a) Radius, (b) velocity and (c) Gibbs free energy of the chemical driving force Gchem

m and the solute drag Gdi f fm for the first nucleated ferrite grain as a function

of temperature during three successive cyclic phase transformations in Fe-0.1C-0.5Mn alloy. (d) and (e) are zoom-ins of the regions marked by rectangular boxes in (c). The data in (c)-(e) are for the first cycle. The arrows illustrate the direction of the transformation during the cycle.

i.e.

|



Gchem

m

|

>

|



Gdi f fm

|

. Thereby,

α

→

γ

proceedsuntil the temperaturereachesitsmaximumatpointD.

3) Inverse stage where

α

→

γ

proceeds during cooling (from point D to E). Although

|



Gchem

m

|

decreases with decreasing temperature after point D,

|



Gchem

m

|

>

|



Gmdi f f

|

still holds. Therefore,

α

_→

γ

continues.

4) Stagnant stage (frompoint Eto G).The magnitudeof



Gchem m decreases with decreasing temperature and remains smaller than



Gdi f f_m .



Gchem

m crosses zero at point F, after which



Gchem

m increases with decreasing temperature, but remains smaller than



Gmdi f f. Therefore, the interface stays immobile untilpointGwhere



Gchem

m =



Gmdi f f.

5) Normal stage where

γ

_→

α

proceeds during cooling (from pointGtothestartofnextcycle)because



Gchem

m >



Gmdi f f. Toinvestigatetheevolution ofthecyclicbehaviorover succes-sive cycles,we plot the radius andthe interfacial velocities asa function oftemperature in Fig. 8for sixferrite grains(including the one shown in Fig. 7). Five of the grains behave similarly in termsofbothR_α(formingaloop)andv_int,whereasonegrainthat nucleatedthelatestshrinksduringheatingafterastagnantperiod andfinallydisappears,showinganabnormalbehavior. This abnor-malbehaviorisalsofoundforsomeothergrainsthatnucleatelate and are relatively smallbefore cyclingbegins. The reasonis that the overlapping of carbon concentration fields from their neigh-boringlargergrains,togetherwiththetemperatureincreaseduring heating, decreasesthechemical drivingforceforferriteformation tostopthecontinuousshrinkage.Inprinciple,thecurvatureisalso expected to play a role here, but thiseffect isnot considered in thecurrentmodelling.

Forthe fivegrainsshowing anormal transformationbehavior, their radii span different ranges and their velocity peaks at dif-ferent temperatures, showing a significant variation in behavior, as hasalso been observedin in-situ EBSD andTEMexperiments

[67,68].Twoofthefivegrains(



tnuc =21.3and106.1sin Fig.8a)

show a considerablylongerinversetransformationstage thanthe rest.Forone ofthem(



tnuc = 106.1s)theinterface doesnotgo

Fig. 8. (a) Radius and (b) interfacial velocity of six ferrite grains as a function of temperature during the first cycle. The nucleation time relative to the earliest nucle- ation event, denoted as t nuc , is marked for each grain, among which t nuc = 113.4

s corresponds to the latest nucleated grain. Arrows illustrate the direction of the transformation during the cycle.

backtotheoriginalposition,whichcoincideswiththeobservation thatf_α doesnotrecovertothestartingvalueuponthecompletion ofthefirstcycle(see Fig.6b).Thereasonisthattheferritegrains whichnucleatedlateinthecyclehavehigherremotecarbon con-centrationintheirparentaustenitegrainsbecauseofimpingement ofthecarbondiffusionfieldsfromearliernucleatedferritegrains. Asaresult,

α

_→

γ

ismorefavoredthan

γ

_→

α

fortheselate nu-cleatedgrains. Itcanbeseen in Fig.8bthat thevint forthegrain

with



tnuc = 106.1s becomesnegativeearlier(i.e.

α

→

γ

starts

earlier)andpositivelater (i.e.

γ

_→

α

startslater)thanother ear-lier nucleatedferrite grains.Forthegrain with



tnuc =113.4 sit

evendisappearsduringthefirstcycle.

Sincethese6grainsselectedincludetheearliestnucleatedone, the latest nucleatedone andthe some nucleatedin between, to-gether they are expectedto representthe spectrumof cyclic be-haviorforallgrains.Itcanbeseenfrom Fig.8bthattheinterface velocityalwaysshowsapeakanditsbehavioragreeswellwiththe recentin-situ TEMobservationontwosingle

γ

/

α

interfaces [68]. The peaks of interface velocities for the normal behaved grains are locatedinthe range803-808°Cfor

α

_→

γ

and~ 840°Cfor

γ

→

α

,whichisincloseagreementwiththeTEMresults. It shouldbe pointedout that afullgrain shrinkageduring cy-clingcannotbe predictedinanypreviousmodelsforcyclic trans-formations becausethey donotconsiderferrite nucleationbefore the cyclingstarts. However, in-situ neutrondepolarization exper-iments showed that the actual number density of ferrite grains during cycling over certain temperature windows can decrease slightly,duetotheoccurrenceofcoarsening(resulting inthe dis-appearanceofsomesmallgrains) [35].Tounambiguouslytrackthe behavior of individual grains duringcyclic transformations, grain resolved 3D/4D techniques based on synchrotron high-energy X-ray diffraction using micro beam [69,70] seems to be the most promisingapproach.

3.3.3. CyclicbehaviorasafunctionofCandMnconcentrations

The modelcanalsobe extendedtopredict thecyclic transfor-mation kineticsasa function of wide rangeconcentrations forC and Mn, which will allow us to elucidate their interplay. Fig. 9

showstheevolutionoff_α foraseriesofFe-C-Mnalloyswith tem-peratureincomparisonwiththeequilibriumfraction.Different be-haviors for f_α are clearly observed. As the Mn concentration in-creases from0.17to 1.0wt.%forFe-0.1C-Mn alloys(Fig.9a-f), an open loopforf_α canclearlybeobserved,whilethestagnantstage (measuredbythetemperaturerangewhereinterfaceisimmobile) becomes longer. It is interesting to note that when theMn con-centration increases to 1.0 wt.%, f_α spans a much smaller range comparedtothelower-Mnalloys.WithfurtherincreasingMn con-centrations,theamplitudeofthechangesinf_α ineachcycleisso smallthattheopenloopishardlyseen.Thismeanstheinterfacial Mn partitioninglargely controlsthecyclickineticsandlimitsany substantialinterfacialmovementsinthesealloys.However,carbon partitioningalsoplaysarolehere.EvenforarelativelyhighMn al-loy,whenCconcentrationdecreases,thetransformationkineticsof both

γ

→

α

and

α

→

γ

couldbeaccelerated.Thisaccelerationis observable whencomparing Fig.9eand 9g. Conversely,whenthe C concentration isincreasedinalloyswitha fixedMn concentra-tion, thenthe kineticsis gradually suppressedandf_α cyclesover a smallerrange.Thisiswell illustratedwhen Fig.9aiscompared with Fig. 9h. Similar behavior is seen when Fig. 9e is compared with Fig.9i.

Previous studies have focused on the effect of a single alloy-ing element(either CorMn)on thecyclictransformation behav-ior [66,71,72].Tocharacterizethestagnantbehaviorduringcycling, the lengthofthestagnantstage (the temperaturerangewheref_α does not change) duringcycling or the length of the post cyclic stagnant stage (the temperaturerangewhere f_α does notchange

afterthethermalcyclingisfinished, whichisduetoresidualMn spikes formed duringcyclinginfront of theinterfaces) are used. It is shown that increasing the C or Mn concentration leads to morestagnant movementofinterfacesduringcycling [71,72].The postcyclicstagnant stage,presentrightafterthe cyclingfinishes, becameevident (> 2.5K) when the Mn concentration increased above 1.0 wt.% Mn [66]. It can be seen from Fig. 9 that C and Mnhaveacoupledeffectonthecyclictransformationbehavior.To quantifythiseffect,itisnotverysuitabletousethelengthof stag-nantstageassuch,ashasbeendone inpreviousstudies,because theequilibriumstatesatthebeginningofthecyclingandthe inter-valofequilibriumfractionsbetweenT1 andT2 aredifferent.Here,

weintegratef_α overthetemperaturenormalizedbyintegrationof thenetpara-equilibriumfactionf_αPE _{surplusfromthatofT}

2 over

thesametemperaturerangeto describethestagnant transforma-tionbehavior(denotedas

η

stag):

η

stag=1−  f_αdT/  T2 T1



fPE α − fαPE

(

T2

)

dT. (12)

Ahighvalueof

η

stag correspondsto amorestagnant behavior

duringthecyclicphase transformation.Intuitively, thenumerator in Eq.12representsthearea ofthecyclicloopin thef_α – Tplot andthedenominatorrepresentsthearea off_α overTasifit was formedbyquenchinganalloyfromT2toT1,thenholdingatT1for

alongtimeandfinallyheatingbacktoT2underpara-equilibrium

conditions.Sincethetransformationbehavior becomesrepeatable fromthesecondcycleon,weusedthedataforthesecondcycleto calculate

η

stagforallalloys.

Fig.10showsacontourplotof

η

stag asafunctionoftheCand

Mn concentrations. The figure clearlyshows that the interface is leaststagnant atthecorneroflow-C andlow-Mnconcentrations, while it is more stagnant when either C or Mn increases. Two dashedlinesdepictedin Fig.10markthetransitionregionbetween a pronounced open loop (the interface movessubstantially when

η

stag<0.85)andanearlyinvisible loopbehavior(the interfaceis

nearlyimmobilewhen

η

stag>0.90)fordifferentalloys,inbetween

theinterface movesgently(0.85≤

η

stag ≤ 0.90).For

η

stag >0.85

thefigureshowsthat thestagnantbehavior ismoreenhanced by increasingtheMnconcentrationandonlyhasaweakdependence ontheCconcentration.Theresultsareinexcellentagreementwith theexperimentalresultsobtainedbyFarahaniandcoworkers [66]. Themetric

η

stagcanbelinkedtotheevolutionoftheCandMn

concentrationprofilesduringcyclicphasetransformation.In Fig.11

weplottheradiusofthefirstnucleatedferritegrainsandtheir C andMn concentrationprofilesacrosstheinterfaceatfiveselected positionsinthefirstcyclefortheFe-0.1C-0.1Mnalloywitha rela-tivelylow

η

stagandfortheFe-0.1C-1Mnalloywitharelativelyhigh

η

stag(asmarked in Fig.10),respectively. Althoughthetwograins

bothshowcharacteristiccyclicbehavior(stagnant,inverseand nor-maltransformations)asdescribedin Section3.3.2,theamplitudes ofthechangeinR_α andtheresultingloopoutlinedbyR_α are sig-nificantlydifferent(see Fig.11aand 11d).PlotsoftheCprofilesat thefiveselectedpositionsindicatethattheevolutionoftheC pro-file is much more pronounced forthe Fe-0.1C-0.17Mn alloy than fortheFe-0.1C-1Mnalloy(see Fig.11band 11e).However,inboth alloystheCprofileinaustenitehasaverysmallconcentration gra-dient, suggestingthat Cdiffusionis not theratelimitingstepfor both alloys. UnlikeC, Mn spikes appear atall positions for both alloys.When

γ

→

α

isfavored(positionsI,IVandV),theMn con-centrationbuildsupattheinterfaceformingpositivespikes,where thepositions thatcorrespond tothestagnant transformation(e.g. positionI andIV)havethe highestmagnitudes.When

α

→

γ

is favored(positionsIIandIII),negativeMnspikesareseenatthe in-terfaceforbothalloys. Themagnitudesofall thespikes aremuch larger inthe Fe-0.1C-1Mn alloythan in theFe-0.1C-0.17Mn alloy,

Fig. 9. Evolution of the ferrite fraction f αduring cyclic partial phase transformation as a function of temperature at a rate of 10 K/min in different Fe-C-Mn alloys. The ferrite

fractions under ortho-equilibrium (Ortho-E), para-equilibrium (PE) and negligible partitioning local equilibrium (NPLE) are also plotted.

Fig. 10. Contour plot of ηstag as a function of C and Mn concentrations. Two dashes lines with ηstag = 0.85 and 0.90, respectively, mark a transition region for distinguishing

Fig. 11. Radius and C and Mn concentration profiles across the interface for the first nucleated ferrite grain during the cyclic phase transformations for the (a-c) Fe-0.1C- 0.17Mn and (d-f) Fe-0.1C-1Mn alloys. Positions I – V correspond to transformation behaviors of stagnant, normal α→ γ, inverse α→ γ, stagnant and normal γ → α. The values of ηstag for the two alloys are marked in Fig. 10 . The Mn profiles are plotted as a function of the normalized quantity r / δ. The inserts in (c) and (e) show

enlargements of the regions marked by the two boxes, respectively.

consistent withthe observationthat

η

stagfortheFe-0.1C-1Mn

al-loyislargerthanthatfortheFe-0.1C-0.17Mnalloy.

4. Conclusions

A novel 3D mixed-model has been developed by considering ferrite nucleationandgrowthgovernedbya balanceintheGibbs free energycontributionsbetweenchemicaldrivingforceand dis-sipation by interface friction and by solute drag for austenite-ferrite transformations inFe-C-Mn alloys. A newimplementation ofthe solutedragdissipationcalculationsisproposed andshown tobecomputationallyveryefficient.Thecapabilitiesofthismodel havebeenillustrated.Basedonthemodelapplicationto austenite-ferritetransformationsinFe-C-Mnalloys,thefollowingconclusions canbedrawn:

1) Thederivedanalyticalsolutionforcalculatingtheenergy dissi-pation due tosolute dragiscomputationally much more effi-cient thantheconventional approachthathasto computethe Gibbsfreeenergiesasafunctionofinterface velocityina pre-definedrange.Thisnewmethodproposed hereisgeneraland canbeusedinother3Dsimulations.

2) Themodelperforms wellinpredictingtheoverall transforma-tionkineticsduringcontinuouscooling,isothermalholdingand inter-criticalthermalcycling.Thegrainsizedistributioncanbe predictedreasonablywellforcontinuouscooling.

3) The transformation behavior of each individual grain can be monitored. Awide variation in behavior wasfound. Thelocal Gibbs free energyofthe chemicaldriving force andthe dissi-pationscanbecalculatedtounderstandthemechanismforthe change ininterfacialvelocity.The carbon partitioningandthe solutedrageffectofMncomplementeachotherdependingon thespecificconditions.

4) Small ferrite grains that nucleated relatively late may fully shrinkuponheatingduringcyclicpartialphasetransformations. This phenomenon has not been predicted by models that do notconsiderferritenucleationbeforethebeginningofthermal cycling.

5) Anew metricis proposed to describe the stagnant interfacial behaviorduringcyclicphase transformations.It isshownthat the coupled effect of C and Mn on the cyclic transformation behaviorcanadequately bedescribed. Atransitionregionasa functionofCandMnconcentrationsisidentifiedtodistinguish differenttypesofcyclicbehavior.

Although we only demonstrate the application of the model for Fe-C-Mn alloys, the model can also be used for calculating thetransformation kineticsfor other ternary alloysonce the cor-respondingthermodynamic informationisavailable.Sincea more efficient implementation of computing solute drag dissipation is proposed, the present work is expected to encourage more ad-vanced 3D simulations using the principle ofthe Gibbs free en-ergy balance. To this aim the code is made publicly available at

https://github.com/haixingfang/3D-GEB-mixed-mode-model.

DeclarationofCompetingInterest

Wedeclarenoconflictofinterest.

Acknowledgements

The authors thank Dr. Hussein Farahani for providing the dilatometrydatafortheisothermaltransformationofthe Fe-0.1C-0.5Mn alloy.H. Fangis grateful to the China Scholarship Council (CSC)forprovidingfinancialsupportforhisPhDstudyatDelft Uni-versityofTechnology.HeisalsogratefultoProf.DorteJuulJensen for providing a postdoc position that is financially supported by theEuropean ResearchCouncil undertheEuropean Union’s Hori-zon2020researchandinnovationprogram(M4Dgrantagreement No.788567).

AppendixA:Solutedrageffectofsubstitutionalelements

Duringthediffusionalphase transformations,substitutional el-ementscansegregateattheinterfacebecausetheinterfacecanbe