Phase Field Modelling of Grain Growth withParticle Pinning and Solute DragbySina ShahandehMaster of Science, Sharif University of Technology, 2007Bachelor of Science, Ferdowsi University of Mashhad, 2004a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoralstudies(Materials Engineering)The University of British Columbia(Vancouver)September 2013c? Sina Shahandeh, 2013AbstractSecond phase particles and solute atoms have been used as an important constituentin the design of materials and processes due to their ability to restrain the motionof grain boundaries. The drag effect occurs on a scale comparable to the particlediameter and interface thickness. However, to simulate grain growth with numericalefficiency one requires a model that captures the drag pressure on the interfacesbut does not resolve the particles or solute segregation spike. In this work, amulti-scale modelling scheme is proposed to simulate grain growth with particlepinning and solute drag. The interaction of a grain boundary with an ensembleof particles is simulated to obtain the pinning pressure. A phase field model isthen developed that incorporates the drag pressure in the meso-scale and simulatesgrain growth. The accuracy of the model is confirmed in comparison with analyticalexpressions. The application of the model is then presented for grain growth in two-and three-dimensional systems under the influence of particle pinning. Measuringthe curvature of the grain boundary network reveals that in the completely pinnedstructure, the average driving pressure is not equal to but lower than the pinningpressure. The results of the nano-scale simulations for pinning pressure is combinedwith the results from the meso-scale to produce a limiting grain size that coincideswith the experiments. This curvature analysis provides a kinetic model thatdescribes the evolution of the structure more accurately than that of the meanfield theories.The proposed phase field formulation is also applied to simulate grain growth inthe presence of solute drag. The grain growth kinetics follows a phenomenologicalrelationship that is described using a power law, with a time exponent in the rangeof 0.35 to 0.50. The deviation from ideal grain growth, associated with a timeexponent lower than 0.5, and its correlation with the solute drag parameters isinvestigated.iiPrefaceThis dissertation is written based on original research conducted by theauthor, S. Shahandeh. All of the work presented henceforth was conductedin the Materials Engineering Department of the University of BritishColumbia, at the Point Grey Campus.Figures 2.3, 2.5, 2.7, 2.8, 2.9, 2.10, 2.11, 2.12, 2.13, 2.14, 3 3.5 and 3.6in the literature review chapters have been taken with permission from thecited sources.The contents of Chapters 6, 8 and 9 are based on the author?s previouslypublished papers. The paper, S. Shahandeh, M. Greenwood, and M.Militzer, Modelling and Simulation in Materials Science and Engineering,vol. 20, no. 6, p. 065008, 2012, has been chosen by the Institute of PhysicsPublishing as a ?Featured Article? in 2012. As the lead investigator, Iam responsible for formation of the concept and for performing all of thesimulations and preparation of the manuscript. M. Greenwood contributedto the formulation of the phase field model and edited the manuscript.Another paper, S. Shahandeh and M. Militzer, Philosophical Magazine,2013, DOI:10.1080/14786435.2013.805277 is to appear in print. M. Militzerwas the supervisory author on this project and was involved throughout theproject in its concept formation and manuscript editing.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xxDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 An Overview of Grain Growth and Drag Phenomena . . . 52.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Ideal Grain Growth . . . . . . . . . . . . . . . . . . . . . . . . 62.2.1 The Kinetics of Ideal Grain Growth . . . . . . . . . . 62.2.2 The Topological Aspects of Ideal Grain Growth . . . . 112.2.3 The Computer Simulation of Ideal Grain Growth . . . 152.3 Grain Growth with Drag Pressure . . . . . . . . . . . . . . . 192.3.1 Grain Boundary Drag by Second Phase Particles . . . 192.3.2 The Kinetics of Grain Growth with Pinning Pressure . 24iv2.3.3 Computer Simulations of Particle Pinning . . . . . . . 262.4 Grain Boundary Drag by Solute Atoms . . . . . . . . . . . . 333 Phase Field Models . . . . . . . . . . . . . . . . . . . . . . . . 393.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2 A Historical Overview . . . . . . . . . . . . . . . . . . . . . . 403.3 Describing Microstructures in Phase Field Models . . . . . . 423.4 The Phase Field Model Construction . . . . . . . . . . . . . . 443.4.1 Thermodynamics . . . . . . . . . . . . . . . . . . . . . 443.4.2 The Kinetics of Evolution . . . . . . . . . . . . . . . . 483.5 The Simulation of Microstructural Evolution in Polycrys-talline Materials . . . . . . . . . . . . . . . . . . . . . . . . . 503.6 The Phase Field Modelling of Particle Pinning . . . . . . . . 533.7 Phase Field Modelling of Solute Drag . . . . . . . . . . . . . 563.8 Conclusions of the Literature Review . . . . . . . . . . . . . . 584 Scope and Objectives . . . . . . . . . . . . . . . . . . . . . . 605 The Methodology and Numerical Techniques . . . . . . . . 615.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615.2 The Numerical Solution of the Phase Field Equation . . . . . 625.2.1 The Finite Difference Method . . . . . . . . . . . . . . 625.2.2 Quantitative Phase Field Model . . . . . . . . . . . . 645.2.3 Numerical Optimization Techniques . . . . . . . . . . 665.3 Post-Processing Algorithms . . . . . . . . . . . . . . . . . . . 685.3.1 Visualization of the Structure . . . . . . . . . . . . . . 685.3.2 The Statistics of Grain Size . . . . . . . . . . . . . . . 695.3.3 The Topology of Grains . . . . . . . . . . . . . . . . . 705.3.4 The Curvature of Grain Boundaries . . . . . . . . . . 716 Ideal Grain Growth Simulations . . . . . . . . . . . . . . . . 756.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.2 Two-Dimensional Systems . . . . . . . . . . . . . . . . . . . . 756.3 Three-dimensional Systems . . . . . . . . . . . . . . . . . . . 806.4 Discussions and Conclusions . . . . . . . . . . . . . . . . . . . 82v7 Modelling the Particle-Grain Boundary Interaction . . . . 857.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 857.2 The Pinning Force of a Particle . . . . . . . . . . . . . . . . . 867.2.1 Energetics Approach . . . . . . . . . . . . . . . . . . . 867.2.2 The Kinetic Approach . . . . . . . . . . . . . . . . . . 897.3 The Pinning Pressure of a Distribution of Particles . . . . . . 937.3.1 The Planar Grain Boundary . . . . . . . . . . . . . . . 937.3.2 Curved Grain Boundary . . . . . . . . . . . . . . . . . 967.4 The Interaction of Grain Boundary Junctions with Particles . 997.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1018 The Friction Pressure Model for Grain Boundary Drag . . 1038.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1038.2 Phase Field Formulation . . . . . . . . . . . . . . . . . . . . . 1048.2.1 Applying Pressure on the Interface . . . . . . . . . . . 1048.2.2 Applying Drag Pressure . . . . . . . . . . . . . . . . . 1068.3 Numerical Testing and Benchmarking . . . . . . . . . . . . . 1098.3.1 Constant Pinning Pressure . . . . . . . . . . . . . . . 1098.3.2 Velocity-Dependent Solute Drag Pressure . . . . . . . 1118.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1129 Grain Growth With Drag Pressure . . . . . . . . . . . . . . 1149.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1149.2 Grain Growth with Pinning Pressure . . . . . . . . . . . . . . 1159.2.1 Grain Structure and Limiting Grain Size . . . . . . . . 1159.2.2 The Evolution of Curvature . . . . . . . . . . . . . . . 1199.2.3 The Kinetics of Grain Growth . . . . . . . . . . . . . 1239.2.4 The Topological Aspects of Grain Growth . . . . . . . 1279.3 Grain Growth with Solute Drag . . . . . . . . . . . . . . . . . 1319.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13510 Overall Conclusions . . . . . . . . . . . . . . . . . . . . . . . 13810.1 Achievements and Contributions . . . . . . . . . . . . . . . . 13810.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 139Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142viAppendix A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157viiList of TablesTable 6.1 Phase field model parameters for ideal grain growth simu-lations in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 76Table 6.2 The kinetics constants of ideal grain growth. . . . . . . . . 83viiiList of FiguresFigure 2.1 The geometrical and statistical aspects of a grain struc-ture are categorized into metric and topological measureswith the goal of attaining a kinetic relationship. . . . . . 8Figure 2.2 The schematic representation of grain boundary curvatureas a function of number of sides. . . . . . . . . . . . . . . 12Figure 2.3 Simulations of grain growth using front tracking method.The grain boundary segments are discretized with smallflat elements. . . . . . . . . . . . . . . . . . . . . . . . . . 17Figure 2.4 Grain size distributions comparing simulation resultsfrom [16] with analytical models. . . . . . . . . . . . . . . 18Figure 2.5 The topological distribution of 3D ideal grain structure. . 19Figure 2.6 The balance of interface tension force for a system withone spherical particle and a curved interface. . . . . . . . 20Figure 2.7 Pinning force of a spherical particle on a grain boundaryas a function of its position. . . . . . . . . . . . . . . . . . 22Figure 2.8 Ratio of limiting grain size to the particle radius vs.volume fraction of second phase. . . . . . . . . . . . . . . 23Figure 2.9 Kinetics of grain growth and grain size distributionpredicted by model of Hunderi and Ryum [70] . . . . . . 25Figure 2.10 Phase field modelling of 2D polycrystalline system con-taining second-phase particles. . . . . . . . . . . . . . . . 27Figure 2.11 Monte Carlo [11] and phase field [12] 3D simulations ofgrain growth with particle pinning. . . . . . . . . . . . . . 28ixFigure 2.12 The simulation of a grain boundary segment interactionwith an ensemble of particles using finite element method[8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 2.13 The velocity of a grain boundary interacting with aparticle and with a distribution of particles modelled withthe finite element method. . . . . . . . . . . . . . . . . . . 30Figure 2.14 Phase field simulation results for grain boundary velocityand grain structures by Apel et al. [79] . . . . . . . . . . 32Figure 2.15 Schematic representation of a grain boundary interactionwith solute atoms. . . . . . . . . . . . . . . . . . . . . . . 33Figure 2.16 The Effect of the grain boundary velocity on soluteconcentration profile and the consequent drag pressure. . 34Figure 2.17 A schematic representation for the variation of the in-terface velocity versus the applied driving pressure forsystems with different solute concentrations. . . . . . . . . 36Figure 3.1 Schematic description of microstructures with field vari-ables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.2 A polycrystalline structure defined by a set of orderparameters, ?i, associated with each grain. . . . . . . . . 44Figure 3.3 The local energy density of a system that decomposesinto two regions as a function of composition and orderparameter. . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 3.4 Homogeneous energy density and order parameters for asystem with two phase fields. . . . . . . . . . . . . . . . . 48Figure 3.5 The evolution of grain structure simulated with 2D phasefield. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 3.6 The grain structures from a phase field model in 3D. . . . 52Figure 3.7 The topological properties of the 3D ideal grain structureas simulated by a phase field model. . . . . . . . . . . . . 53Figure 3.8 Cross-sectional snapshots of the interaction of a sphericalgrain boundary with six particles at different times [9]. . . 54xFigure 3.9 Interaction of a grain boundary with a particle and thechange in total energy of the system. . . . . . . . . . . . . 55Figure 3.10 Concentration spike at an interface for three differentinterface velocities modelled with phase field [87]. . . . . . 58Figure 5.1 A 2D and a 3D stencil for calculation of the Laplacianterm in the phase field equation. . . . . . . . . . . . . . . 63Figure 5.2 Schematic representation of the sparce data structure. . . 68Figure 5.3 Images of a 2D grain structure constructed by ? andMaxind arrays. . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 5.4 The Maxind array of a 3D structure is analyzed tocalculate the number of faces of a grain. . . . . . . . . . . 71Figure 5.5 A circular order parameter and its corresponding curva-ture field. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 5.6 Colour map of positive and negative values of the curva-ture fields for all order parameters. . . . . . . . . . . . . . 73Figure 6.1 The grain boundary network in 2D ideal grain growthsimulations. . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 6.2 Normalized growth rate of each grain as a function ofequivalent radius for ideal grain growth in 2D. . . . . . . 77Figure 6.3 The self-similar grain boundary curvature distribution for2D ideal grain growth simulations. . . . . . . . . . . . . . 78Figure 6.4 Equivalent grain radius as a function of normalized timefor different simulations in 2D. . . . . . . . . . . . . . . . 80Figure 6.5 Curvature map of the grain boundary network for idealgrain growth in 3D. . . . . . . . . . . . . . . . . . . . . . 80Figure 6.6 Normalized growth rate of individual grains as a functionof size for ideal grain growth in 3D. . . . . . . . . . . . . 81Figure 6.7 Distribution of the grain boundary curvature at differenttime steps in 3D ideal grain growth. . . . . . . . . . . . . 82Figure 6.8 Equivalent grain radius as a function of normalized timein 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82xiFigure 7.1 Energy density of a system with a particle and a grainboundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . 87Figure 7.2 Cross-section snapshots for a 3D system with a movinggrain boundary and a particle. . . . . . . . . . . . . . . . 87Figure 7.3 The total energy of the system and its derivative as afunction of the grain boundary position. . . . . . . . . . . 89Figure 7.4 The iso-surfaces of ? illustrating the interaction of fourparticles at the centre of curved grain boundaries. . . . . 90Figure 7.5 Velocity of grain boundary and the corresponding pinningforce for interaction of a particle with a curved grainboundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 7.6 Pinning force as a function of the grain boundary positionwith respect to the particle. . . . . . . . . . . . . . . . . . 92Figure 7.7 A planar grain boundary moves through a distribution ofparticles. . . . . . . . . . . . . . . . . . . . . . . . . . . . 94Figure 7.8 The instantaneous grain boundary velocity and pinningpressure for a planar grain boundary. . . . . . . . . . . . . 95Figure 7.9 Motion of a curved grain boundary through a dispersionof particles. . . . . . . . . . . . . . . . . . . . . . . . . . . 96Figure 7.10 Velocity of the grain boundary and the normalized pin-ning pressure for a curved grain boundary. . . . . . . . . . 97Figure 7.11 Interaction of triple junctions and quadruple points withparticles. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100Figure 7.12 The velocity of the grain boundary and the pinning forcefor the interaction of junctions with particles. . . . . . . . 101Figure 8.1 Local energy density as a function of order parameters ina system with two grains. . . . . . . . . . . . . . . . . . . 105Figure 8.2 Schematic representation of the friction pressure on aninterface. . . . . . . . . . . . . . . . . . . . . . . . . . . . 107Figure 8.3 The grain boundary velocity of a shrinking circular grainas a function of driving pressure in phase field simulationswith pinning pressure. . . . . . . . . . . . . . . . . . . . . 110xiiFigure 8.4 The normalized grain boundary velocity of a shrinkingcircular grain as a function of normalized driving pressurein the presence of solute drag pressure. . . . . . . . . . . . 111Figure 8.5 A four-sided grain shrinks in a symmetric boundarycondition. . . . . . . . . . . . . . . . . . . . . . . . . . . . 112Figure 9.1 Snapshots of a 2D grain growth simulation with pinningpressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115Figure 9.2 Equivalent grain radius and limiting grain size for 2Dgrain growth simulations with pinning pressure. . . . . . . 116Figure 9.3 Grain structure of a 3D simulation with the pinningpressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117Figure 9.4 Equivalent grain radius and limiting grain size for 3Dgrain growth simulations with pinning pressure. . . . . . . 117Figure 9.5 The experimental data on limiting grain size as a functionof particle size and volume fraction is compared with thesimulations. . . . . . . . . . . . . . . . . . . . . . . . . . . 118Figure 9.6 The Grain structure and associated curvature distribu-tions for a 2D system with pinning pressure. . . . . . . . 120Figure 9.7 The grain structure and associated curvature distribu-tions for 3D system with pinning pressure. . . . . . . . . . 121Figure 9.8 ?1 as a function of R? during 3D grain growth for threesystems with different pinning pressures. . . . . . . . . . . 122Figure 9.9 Kinetics of grain growth with pinning compared to ana-lytical models. . . . . . . . . . . . . . . . . . . . . . . . . 124Figure 9.10 Fitting Andersen?s equation on the simulation data for a3D system. . . . . . . . . . . . . . . . . . . . . . . . . . . 125Figure 9.11 The normalized growth rate of individual grains as afunction of relative size for 3D simulations with pinningpressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127Figure 9.12 The normalized growth rate as a function of the numberof faces of a grain for a 3D system. . . . . . . . . . . . . . 129xiiiFigure 9.13 The normalized growth rate of individual grains as afunction of the number of faces for the grain growth withpinning pressure. . . . . . . . . . . . . . . . . . . . . . . . 131Figure 9.14 Snapshots of 2D grain growth with solute drag pressure. . 132Figure 9.15 The effect of solute drag pressure on the kinetics of graingrowth for simulations with different b values in Cahn?sequation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 133Figure 9.16 The time exponent of grain growth (?) as a function ofthe solute drag parameter, b. . . . . . . . . . . . . . . . . 134xivNomenclature? Constant in Hillert?s grain growth kinetics law?1 Proportionality constant between average curvature and aver-age grain size?i External angle between two faces of a grain intersecting at atriple junctionK? Average curvature of the grain boundary network, [m?1]?N Area fraction of N -sided grains in 2D?t Time step duration, [s]?x Numerical grid spacing, [m] Coefficient in local energy density of a system with particle?i Phase field order parameter i? Concentration of solute atoms at the grain boundary? Energy gradient coefficient, [Jm?3]? Time exponent of grain growth kinetics?d Distance between nearest vertices of a grain, [m]r? Axis of the curvilinear coordinate system normal to iso-?isurfacesxvr Position vector of a grid point in the domainF Total energy of the system, [J ]S Surface area of the grain boundary network, [m2]?R Ratio of surface area of a grain to the cubic root of its volume,[m]? Radius of curvature of grain boundary, [m]? Interface or grain boundary energy, [Jm?2]? Contact angle between interface and particle? Exponent of fv in limiting grain size relationship% Grain size normalized by Rcr? Proportionality constant for limiting grain sizeA Area of a grain in 2D, [m2]a, b Constant in Cahn?s solute drag pressure relationshipa0 Lattice constant, [m]Ap/gb Intersection area between a particle and a grain boundary, [m2]C0 Concentration of solute atoms in the lattice, [m?3]C1, C2, C3 Phenomenological constants in Louat?s grain size distributionD Diffusivity of solute atomsE Interaction energy of solute atoms and grain boundaryei Length of a triple line, [m]f0 Local energy density, [Jm?3]fa Area fraction of particles in 2D systemsxvifK Probability density function of the grain boundary networkcurvature distributionfN Probability density function of the grain face number distribu-tionfR Probability density function of the grain equivalent radiusdistributionfv Volume fraction of particlesFZ Zener pinning force, [Jm?1]G(u) Scaling function that relates growth rate of a grain to theaverage grain growth rateG1, G2 Scaling function in Mullins topological relationshipGH Scaling function in Hilgenfeldt topological relationshipK Local grain boundary curvature, [m?1]k Grain growth kinetic constantK1, K2 Principal curvatures of a surfacekB Boltzmann constants, [JK?1]ks Fitting parameter for the kinetics of grain growth with solutedragKetai Curvature field associated with the order parameter ?iL Phase field kinetic coefficient, [J?1s?1]lx, ly, lz Length of the phase field system in each direction,[m]lgb Thickness of a grain boundary in phase field model, [m]M Interface or grain boundary mobility,[m4J?1s?1]xviim local energy density coefficient, [Jm?3]m1, m2 Constants in the local energy density for a system with oneorder parameter, [Jm?3]Meff Effective (observed) mobility obtained phenomenologically,[m4J?1s?1]N Number of faces of a grain in 3D and number of sides of a grainin 2Dn Dimensions of the systemNp Number density of particles intersecting with a grain boundary,[m?2]Nv Number of atoms per unit volumePd Driving pressure on an interface, [Jm?3]Pf Drag (friction) pressure, [Jm?3]Ps Solute drag pressure, [Jm?3]Pz Pinning pressure on a grain boundary, [Jm?3]Pz Pinning pressure, [Jm?3]R Grain equivalent radius, [m]R? Radius of grains with maximum frequency in the grain sizedistributionrp Particle radius,[m]Rcr Equivalent radius of grains with the critical size that do notshrink or growS Surface area of all faces of a grain in 3D, [m2]s Perimeter length of a grain in 2D, [m]xviiiT Temperature, [K]t time, [s]tn Time step numberu Grain size normalized by R?ucr Relative size of the critical grain size (Rcr/R?)V Volume of a grain, [m3]v Velocity of interface or grain boundary, [ms?1]x Average position of a grain boundary, [m]z Zener factor; Pinning pressure in unit of ?fv/rpPDE Partial Differential EquationL Mean Width measure of a grain, [m]R? Average grain equivalent radius, [m]R?lim Limiting grain size; average grain size of a fully pinnedstructure, [m]R?0 Initial average grain radius, [m]xixAcknowledgementsFirst of all, I would like to thank my supervisor, Prof. Matthias Militzerfor providing continuous support and wisdom throughout my studies atUBC. Secondly, I had a pleasure of spending wonderful time with colleaguesand friends at the Microstructural Engineering group at the Department ofMaterials Engineering. In particular, I am thankful to Michael Greenwoodfor mentoring me and stimulating discussions. I will never forget ourjoyous discussions over drinks with Michael, David, Abhijit, Arnaud, Beth,Guillaume, Jennifer, Phil, Tegar, Thomas and many more friends...I am immensely grateful for the financial support of Natural Scienceand Engineering Research Council of Canada and for my scholarshipsfunded by the government of British Columbia and the University of BritishColumbia. This project was possible by utilizing computational resourcesof the Compute Canada, the West Grid clusters.I am thankful to Olga for her friendship and support during the lastfour years. My life in Vancouver was truly enriched through her enterprisingspirit, by her music and kindness.Above all, I would like to thank my parents for all their love and care. Iam forever indebted to them for their many sacrifices.xxDedicationThis work is dedicated to the countless intelligent children indeveloping countries whose lack of wealth and privilege havedeprived them of pursuing their natural curiosity. Your sadlymissed contribution toward the fields of science reminds us ofour own responsibility to create a world in which the science andtechnology is in the service of all humanity...xxiChapter 1IntroductionThe evolution of cellular networked structures occurs in many physicalsystems from soap froths to grain growth in polycrystalline materials. Inmaterials, this phenomenon has important technological applications inprocesses with high homologous temperatures, where grain boundaries aremobile and structural coarsening occurs. The ability to understand andcontrol materials? polycrystalline structure is the key in obtaining the desiredphysical and mechanical properties in materials ranging from structural steelto copper thin films in microelectronics. To develop a desired microstructure,a material may be heat treated starting from different initial structures, e.g.deformed and recrystallized or as electrodeposited structures, allowing graingrowth to occur. In some materials, e.g. in high strength and fracturetoughness linepipe steels, fine grains are desirable, whereas in others,e.g. in electrical steels, coarse grains with a particular crystallographicorientation, are preferred. The variety of processes and alloying techniquesused for developing a microstructure is virtually unlimited, and in all cases,controlling the motion of grain boundaries is paramount.Solute atoms and/or second phase particles can be used to reduce therate of the grain boundary migration or to halt its motion [1, 2]. As anexample, in high strength low alloy (HSLA) steels, the addition of smallamounts of micro-alloying elements creates fine carbide precipitates [3].These particles limit the motion of grain boundaries by a particle pinning1mechanism. In addition, solute atoms in the matrix segregate to grainboundaries and impede their movement by solute drag [4].Employing particles to restrain grain growth kinetics has been studied forseveral decades since the seminal work of Zener and Smith [5]. The effectof particles is explained by considering a pinning pressure which opposesthe motion of the grain boundaries. The goal of the pinning theories is todetermine the kinetics and morphology of the grain structure with pinningand to estimate the limiting grain size beyond which no grain growth canoccur. The latest analytical theories have used complex geometrical modelsfor the interaction of the grain boundaries with particles, to ascertain thepinning force [6, 7]. Numerical simulation tools such as the finite element[8] and phase field models [9], have been developed recently to study thepinning process in more detail. However, in spite of decades of research ondeveloping a particle pinning theory [10] and the recent advancement of thecomputer models e.g., the Monte Carlo [11] and phase field [12] simulations,a theory to accurately describe all aspects of grain growth with particlepinning has not yet been developed.In the solute drag phenomenon, solute atoms are the obstacles to grainboundary motion. The atomic structure of grain boundaries causes thesegregation of solute atoms to occur. The solute atmosphere interacts withthe moving grain boundary and causes a drag pressure. A continuum theoryfor solute drag, proposed by Cahn [13] and then by Hillert [14], forms thebasis of the current understanding of the phenomenon. The focus of thesetheories has been to predict the level of segregation and the drag pressureon the interface. However, works studying the effect of solute drag on thegrain microstructure has been limited.An extensive body of literature has been dedicated to describing, the-orizing and simulating microstructural changes in polycrystalline systems.Traditionally, phenomenological observation and mean field theories havebeen utilized. In the past two decades, numerical modelling techniques haveemerged as powerful tools for investigating many details about the kineticsof grain growth and the morphology of grain structure [15]. Currently, large-scale computer simulations [16, 17] can be performed to refine the theories2and explain more details about the ideal grain growth dynamics. However,the application of such models to study grain growth with drag is limited.The challenge facing computer models is that the interface drag phenomenonis a multi-scale problem. The interaction of grain boundaries with particlesor solutes occurs on the scale of the interface thickness (? 1 nm) or theparticle diameter (1-100nm). On the other hand, the grain structure evolveson the length scale of the grain size (1-100 ?m). It is computationallyunfeasible to create a model that resolves a particle or solute atmosphereand at the same time simulates the evolution of microstructure. Therefore,a proper description of the phenomenon requires a microstructural evolutionmodel that has been formulated on the length scale of the microstructure,while appropriately including the effect of interface drag.The phase field method is a versatile modelling technique and has beenutilized in the simulation of many aspects of phase transformations. Sincethis model utilizes a continuum field, there is no need to track the positionof the interfaces. This technique is particularly suitable for grain growthsimulations where the complex motion of the grain boundary network canbe implemented with ease. Diffusion of solutes, external field e.g. stress andother meso-scale variables, e.g. thermal gradient can be readily included inthe model.Outline of the TextThis thesis uses the phase field modelling technique to study the effect ofparticle pinning and solute drag on grain growth kinetics. In the first twochapters of this thesis, the current state of knowledge about grain growthand drag phenomena including particle pinning and solute drag is reviewed.Then, an overview on the phase field modelling, is presented. The focus is onthe application of the phase field method to the simulation of grain growthwith drag. The numerical techniques for implementation of the model inlarge-scale computer simulations and post-processing algorithms is describedin Chapter 4. A new method for measuring grain boundary curvature inphase field models is presented. The curvature analysis technique is then3applied to ideal grain growth structures where fundamental quantities suchas the distribution of driving pressure and geometrical constants of graingrowth is are obtained. Chapter 6 deals with the simulation of a grainboundary interaction with particles. The length-scale of these calculationsis small enough to resolve the particles and the pinning force of individualparticles on the grain boundary. In Chapter 7, the phase field model capableof simulating drag pressure on the meso-scale is introduced and its accuracyis demonstrated. We then use this formulation to simulate grain growthwith the presence of pinning pressure and solute drag. The last chapterconcludes the work and suggest directions for further research.4Chapter 2An Overview of GrainGrowth and DragPhenomena2.1 IntroductionResearch on the physical process of grain boundary migration and graingrowth has a long history in metallurgy and materials science. It waswell known in the 17th century that some materials form crystals. Forinstance, Kepler [18] hypothesized that the six-fold symmetry of ice crystalsis due to the fact that the ?water particles? solidify in a hexagonalarrangement. The observation of small grains in metallic materials wasnot possible until optical microscopes and metallographic techniques weredeveloped during the period of 1880 to 1910 [19]. In metallurgy, theinterest in studying grain structures has emerged through the observationof the correlation between the softening of as-deformed metals after hightemperature annealing and their microstructure. This softening is followedby the formation of small crystals and their subsequent coarsening. Theemergence of crystallography and x-ray diffraction techniques then paved theway for the modern understanding of the grain structure in polycrystalline5materials. As an alternative view, the grain structure can be considereda network of connected two-dimensional (2D) crystalline defects, i.e. thegrain boundaries. In this overview, the focus is on the analytical theoryand computer simulation of grain growth kinetics. First, the evolution ofan isotropic grain boundary network in pure materials is considered. Inthe next section, theories on grain growth with the presence of obstaclesincluding second phase particle and solute atoms is reviewed.2.2 Ideal Grain Growth2.2.1 The Kinetics of Ideal Grain GrowthThe concept of ideal grain growth is an approximate model to describe thecoarsening of the grain structure of polycrystalline materials. It assumesthat all grain boundaries have the same interface energy and mobility,and that the velocity of a grain boundary?s motion is only proportionalto its local curvature [2]. Ultra-pure single phase materials with close-to-isotropic interface properties exhibit ideal grain growth behaviour at hightemperatures. In this condition, the average equivalent radius of grains, R?,increases parabolically with time [20]:R?2 ? R?20 = M?kt, (2.1)where R?0 is the initial grain size, ? is the grain boundary energy, M is thegrain boundary mobility and k is a constant. To explain this observation,Burke and Turnbull assumed that the rate of migration of a grain boundaryis proportional to its capillarity driving pressure Pd or inversely proportionalto its radius of curvature, ?:v = MPd = M??, (2.2)where v is the grain boundary velocity. They assumed that the curvature is6proportional to the average grain size (? ? R?) and that the rate of growthof the average grain size is proportional to its migration rate (dR?/dt ? v).As a result, the rate of grain growth can be written as:dR?dt= M?k21R?. (2.3)Solving Eq. 2.3 for R? produces the parabolic grain growth law (Eq. 2.1)[20]. Since the 1950s, many studies have been dedicated to finding ananalytical relationship for the rate constant k, as well as for the grain sizedistribution. The main factor that makes this problem more complicatedcompared with that of processes such as particle coarsening is the complexconnectivity of the grain boundary network. The growth or shrinkage rateof a grain not only depends on its size but is also function of the size of itsneighbours. This connectivity between different elements of the structureconstitutes the topological aspects of the grain structure. For example, thenumber of faces of a grain, N , identifies the growth-shrinkage pairs betweena grain and its neighbours. Fig. 2.1 schematically illustrates the geometricand statistic aspects of a grain structure.The theoretical study of grain growth is conducted to establish rela-tionships between the different aspects of grain growth shown in Fig. 2.1.Feltham [21] was first to obtain a kinetic relationship for individual grains aswell as the parabolic grain growth law using an experimentally determinedgrain size and grain?s face number distributions. He suggested a relationshipfor the growth rate of individual grains using a phenomenological log-normaldistribution for the grain size:dRdt=M?Rln(RR?), (2.4)Here R? is the radius with the maximum frequency in the grain sizedistribution. Hillert [22] later postulated a relationship for the grain growthrate as a function of the grain size. His theory is based on the assumptionthat on average, a grain is surrounded by grains of a critical radius Rcr. Ifthe grain is smaller than Rcr it shrinks and if it is larger it grows. The rate7Figure 2.1: The geometrical and statistical aspects of the grainstructure (d) are categorized into metric (a-c) and topological (f-h) measures with the goal of attaining a kinetic relationship (e).Grain size distribution and face number distribution are shownin (a) and (f), respectively. The number of faces of grains (N)is related to their size (g). The individual growth rate of a grainis both function of its size (c) and its topological class (h). Thedriving pressure for the evolution of the structure is determinedby distribution of the grain boundary network curvature, K (b).of growth is given by:dRdt= M??(1Rcr?1R). (2.5)Here, ? is a proportionality constant, with ? = 1 for 3D systems8and ? = 0.50 for 2D systems. The relationship between the grain sizedistribution, fR(R), and the individual grain growth rate can be obtainedfor any n?dimensional system under two conditions: (i) the conservation ofvolume should be maintained, i.e.? ?0RndRdtfR dR = 0, (2.6)and (ii) the continuity equation for the grain size distribution should besatisfied:?fR?t+??R(fRdRdt)= 0. (2.7)The grain size distribution in ideal grain growth remains self-similar, i.e.at any time, the distribution function of R/R? remains the same. For thiscondition to be satisfied, Hillert used a stability criterion similar to that ofthe condition used for particle coarsening by Lifshitz and Slyozov [23] andobtained the self-similar grain size distribution given by [22]:fR = (2 exp(1))n ?n%(2? %)2+n? exp(?2n2? %), (2.8)where % = R/Rcr, and n = 2 and 3 for 2D and 3D systems, respectively.Eq. 2.8 can be averaged to obtain R? = Rcr for 2D and R? = 8/9Rcr for 3Dsystems. The parabolic grain growth kinetics law was also derived as:dR2crdt=12M??t. (2.9)In consequence, Hillert?s theory predicts rate constants of k = 0.25 for2D and k = 0.39 in 3D for the parabolic grain growth kinetics. This theorydoes not consider any topological aspect of the grain structure. Althoughit has become the cornerstone of further investigations, Hillert?s grain sizedistribution and the kinetic constant do not match the results obtained fromnumerical simulation of ideal grain growth, nor the experimental results.Louat [24] proposed a fundamentally different approach. He argued that9Hillert?s assumption about the growth rate of grains (Eq. 2.5) is subjectto a substantial scatter and can only be applied to the average of theensemble. He then postulated that the evolution of the structure is aresult of the statistical movement of different grain boundary segments.This stochastic process influences the grain size distribution function viaa ?diffusive? differential equation:?fRdt= C1?2fR?R2, (2.10)where C1 is a rate constant. The solution for this equation with theboundary condition fR(R = 0) = 0 is:fR(R, t) = C2Rexp(?R2/4C1t)C3t3/2, (2.11)where C2 and C3 are constants that are adjusted based on the initial grainsize distribution. The grain size distribution given by Eq. 2.11, also calledthe Rayleigh distribution in statistics, is very similar to the log-normaldistribution observed in experimental results.Efforts to refine Hillert?s and Louat?s theories continued by authorssuch as Aboav and Langdon [25], Hunderi and Ryum [26] and Rios [27?29]proposing other grain size distributions. The problem is that the proposedgrain size distributions eventually require an adjustable (non-physical)parameter to match the experimental distribution [29]. This is becausethe growth rate of individual grains and the grain size distribution arerelated through the continuity equation. For any self-similar distribution, aparticular growth kinetic law can be achieved. Therefore, to specify a uniquesolution that is consistent with the physical nature of the grain boundarynetwork?s connectivity, the topological aspect of grain growth should beconsidered.102.2.2 The Topological Aspects of Ideal Grain GrowthThe kinetics of grain growth depends on the structural connectivity, i.e.topology of the grain boundary network (Fig. 2.1). For example, thenumber of faces of a grain, N , determines the number of neighbours thegrain possesses. A grain with a greater number of neighbours is larger andincreases in size whereas a small grain has fewer neighbours and shrinks.A topological relationship for ideal grain growth in 2D has been proposedby von Neumann and Mullins [30, 31]. This relationship is directly derivedfrom the fact that the velocity of a grain boundary is v = M?K, where Kis the grain boundary curvature. Therefore, the rate of change in the areaof a grain, A, is:dAdt= ?M??sKds = ?M?(2pi ?pi3N), (2.12)where S defines the grain boundary length surrounding the grain. Theintegral of curvature over a closed curve is always 2pi, except for that ofa polygonal grain, where abrupt changes occur at each junction (Fig. 2.2).The equilibrium angular difference at the triple points is pi/3, and there areN junctions for an N -sided grain. Therefore, for each junction, a value ofpi/3 is subtracted from the integral:dAdt=M?pi3(N ? 6). (2.13)The von Neumann-Mullins relationship is applicable to a structure atfull equilibrium, where the grain boundary contact angle at all the triplejunctions is 120?. Grains with six sides have flat grain boundary segmentsand do not change their area. Five- and seven- sided grains have concaveand convex segments, respectively, and shrink or grow at a constant rate(Fig. 2.2).The von Neumann-Mullins relationship (Eq. 2.13) enables one to obtaina link between the topological distribution of the grain structure and theaverage kinetics. The average grain growth rate is related to the individualgrowth rate of grains through a scaling function G(u):11?/3 ?/3?/3(a) (b) (c)Figure 2.2: The schematic representation of grain boundary curvatureas a function of number of sides. A five-sided grain (a) hasconcave grain boundary segments and shrinks, while grains withmore than 6 sides (c) have convex segments and grow. The flatsegments of a grain with six sides (b) do not move.dR/dt = G(u)dR?/dt, (2.14)where u = R/R?. The function G(u) does not change with time for aself-similar structure. The topological distribution of the grains can becharacterized with a histogram as shown in Fig. 2.1f. It can be perceived thatthe size of a grain is directly proportional to its number of faces (Fig. 2.1g).It was shown by Mullins [32] that an invariant scaling function, s?(u), can beintroduced that gives the number of sides of a grain with the relative sizeu. He then derived a relationship between the growth rate of a grain, theparabolic grain growth constant, the grain size distribution and s?(u) [32]:G(u) = u?2fR(u)? umaxufR(u?)du? =16ku(s?(u)? 6) . (2.15)Eq. 2.15 defines the relationship between k, s?(x) and fR where, if thetwo are given by experiments, the third can be derived. The rate constantcan also be determined from a structure with:12k =23R?2R2R4?N?N (N ? 6), (2.16)where, ?N is the fraction of the total area of the structure occupied byN -sided grains. This still does not provide an analytical theory to directlyobtain the kinetics of grain growth, and requires experimental data [32] orgeometrical assumptions [33?36]. To date, there is no analytical theory forideal grain growth in 2D that can provide accurate results for the kineticsand grain size distribution.The complexity of grain growth in 3D is even more pronounced. Mullinsinitially proposed an approximate relationship between the number of facesof a grain and its volumetric growth rate [37]:V ?13dVdt= M?G1(N)G2(N), (2.17)where V is the volume of a grain and,G1(N) =[pi3? 2 tan?1(1.86?N ? 1N ? 2)], (2.18)and,G2(N) =(34pi) 13[5.35N23(N ? 22?N ? 1?38G1(N))?13]. (2.19)Mullins? topological expression for the rate of growth of a grain wasobtained by assuming that the grains are regular polyhedrons. Then, it wasused to derive a relationship for the parabolic grain growth constant in 3D:k =u?2piu6? [fR(u)u4?NfNG1G2]du. (2.20)Here, fN is the distribution function for grains of topological class Nshown in Fig. 2.1f. This relationship is analogous to Eq. 2.16 in whichG1G2 plays the role of the von Neumann-Mullins equation. The need for13experimental data still remains. Mullins proposed k = 0.50 for 3D graingrowth using an experimentally determined probability distribution functionfor the grain size, fR, and the grain face numbers, fN [37].More recently, Hilgenfeldt et al. [38] proposed another analyticalrelationship for the volumetric growth rate:V ?13dVdt= M?GH(N) (2.21)where,GH(N) =12223(34pi) 13tan13(?N2)(pi3? ?N)((N ? 2) tan(pixN)) 23,(2.22)?N = 2 tan?1(?4 sin2(pixN)? 1), (2.23)xN = 6(1?2N). (2.24)The functions G1G2 and GH are positive for grains with equal or morethan 14 faces and negative for grains with 13 or less faces. This shows that,unlike six-sided grains in 2D, no topological class in 3D exists having flatfaces and lacking growth.A fundamental relationship for the growth rate of a 3D grain that doesnot include any approximations was derived by MacPherson and Srolovitz[39]:dVdt= M??SKdS = ?2piM?(L(D)?16?eei(D))(2.25)Here, L(D) is a one-dimensional measure for the grain D, ei is the lengthof triple line i and the sum is over all the triple lines surrounding the grain.The measure, L, is called the mean width (not to be mistaken for the graindiameter or linear intercept length) and can be simplified for a polyhedron,14as in [39]:L =12pi?e?iei, (2.26)where ?i is the external angle between the two faces of a grain intersecting ata triple junction. Although Eq. 2.25 is general and describes the governingdynamics of grain growth in 3D, it is not a purely topological relation.However, Rios et al. [40] recently suggested an approximate relationship forthe mean width as a function of the number of sides of a regular polyhedron:L?d? 2.07(?N ? 0.96), (2.27)where ?d is the distance between the nearest vertices. Recently, Lazar etal. [41] used a more comprehensive description for the topology of thecellular structures given by Matzke [42] to differentiate between an idealgrain growth structure and a random Voronoi tessellation network. Theyshowed that the number of faces of a grain in 3D is a poor method fortopological characterization. It has been confirmed by computer simulationsthat 3D grains with the same number of faces can have different normalizedgrowth rates [43].A comprehensive theory of grain growth requires both a kinetic re-lationship, as expressed in von Neumann-Mullins equation and its 3Dequivalents, and a statistical distribution of topological classes. It hasbeen proven challenging to gain insight through analytical methods into thecomplex arrangement of cellular structures and their statistical attributes.Alternatively, a better picture for the process of ideal grain growth can berevealed by simulating the evolution of the microstructure.2.2.3 The Computer Simulation of Ideal Grain GrowthA natural way to numerically represent a grain microstructure is to discretizethe structure into pixels (e.g. similar to that performed for a digital image)15and to associate a grain index to each pixels. In this method, each area (orvolume) associated with a grain have a unique index. Grain boundariesare distinguished when neighbouring pixels have unequal indices. Theevolution of this structure is facilitated by changing the indices of the pixelsat the grain boundaries. The governing dynamics of the change can bederived from probability rules using the Monte Carlo algorithm [44]. Inthis method, the probability of an exchange of an index is a function ofits neighbours? configurations to decrease the energy of the system [45, 46].Using this simulation approach, the parabolic grain growth kinetics wereobtained, along with other statistical aspects of the grain structure such asthe grain size and its topological distribution [44, 47, 48]. In this method,it is convenient to construct a structure; however, tracking the temporalprogress of a microstructure is not straight-forward since there is no physicalconnection between each Monte Carlo time step and the microstructureevolution time.Grain growth can be simulated by modelling the grain boundaries ascurved surfaces and tracking the position of each surface. At each timestep, the velocity of a grain boundary is proportional to its curvature.This method is used in vertex models [49?51] where only the position ofthe junctions is stored. The curvature of the grain boundary segments iscalculated based on the assumption that the junctions are at equilibrium;e.g. in 2D, there is a 120? angle between the segments. In vertex models,the topological changes of grains, i.e. the disappearance or creation of ajunction, have to be accounted by explicit instructions in the algorithm. Asa result, this technique has a limited application to 3D systems where it isdifficult to account for all topological variations.Another front tracking technique, the surface evolver, has been intro-duced [52, 53] that models a grain boundary segment with a discretizednumerical mesh. Resembling finite element models, each element on thesurface of a grain boundary evolves based on the local curvature andinteracts accordingly with neighbouring elements (Fig. 2.3a). Topologicalvariations can be more easily tracked by adjusting the elements in a casewhere a face or segment disappears or is created. Lazar et al. [17, 54]16recently developed a modified front tracking scheme in which the volume(area) change of a grain is subject to the condition that it obey theMacPherson-Srolovitz (or von Neumann-Mullins) relation. As shown inFig. 2.3b, in this method, the variation of curvature on the grain faces islimited to fully convex or concave due to its radial discretization.(a) (b) (c)Figure 2.3: Simulations of grain growth using front tracking method.The grain boundary segments are discretized with small flatelements in a simulation by (a) Wakai et al. [53] and (b) Lazaret al. [17](c) The grain boundary network is an aggregate of theinterconnected elements [17].Another major technique for simulating grain growth utilizes the phasefield models that are reviewed in detail in Sec. 3.5. Level set methods havealso been also employed to describe the evolution of grain structure. Inthis model an iso-surface of a smooth function (or a set of kernel functions)identifies a grain boundary segment. A junction can be distinguished bythe intersection of two or more iso-surfaces [55]. The motion of the grainboundaries are modelled by solving a stable linear differential equation thatreproduces a curvature driven motion for the iso-surface. Similar to phasefield models, the topological transition of the grains is implicitly modelledwith the interaction of the kernel functions. A large scale simulation ofideal grain growth was performed by Elsey et al. [16, 56]. Fig. 2.4 comparesnormalized grain size distribution from a simulation with analytical models.Different simulation methods provide relatively close values for the graingrowth rate constant k in Eq. 2.1. In 2D simulations, the following valueswere obtained: 0.27 with the phase field [57], 0.27 with the vertex model [50],170 0.5 1.0 1.5 2.0 2.5 3.000.20.40.60.81R /RfrequencyLouat (Eq. 2.11)LognormalWiebullRios [28]Hillert (Eq. 2.8)Figure 2.4: Grain size distributions comparing simulation results from[16] with analytical models.0.29 with te Monte Carlo simulation [58] and 0.26 with the front-trackingmodel [54]. The k values from 3D simulations were reported to be: 0.5 [59]and 0.49 [60] for phase field simulations, 0.43 for Monte Carlo simulations[61] and 0.42 for the level set method [16]. Lazar et al. [17] reported k = 0.39with the front tracking method which is smaller than the value of othersimulation techniques.Other metric and topological aspects of the simulated grain structuresare also reported to be similar for different models. Fig. 2.5 compares theresults of several simulations methods for the distribution of the number offaces of grains. Evaluating the advantages or disadvantages of a particularmodelling method is a rather subjective endeavour. It is not conclusivelyclear which method is faster or more efficient since computational time isusually not reported and depends significantly on the nature of the computerhardware, software and numerical optimization algorithms employed. Theversatility of a simulation technique is also an important factor since,eventually, the technological application of studies of structural evolutionof polycrystalline materials does not remain limited to ideal grain growth,but is extended into more complicated systems with anisotropy, inclusion ofimpurities, solutes or second phase particles.18Level setPhase fieldMC1MC2Front trackingVertexFigure 2.5: The topological distribution of 3D ideal grain structureobtained from level set [17], phase field [62], Monte Carlo [48,63], front tracking [53] and vertex [64] simulation techniques.2.3 Grain Growth with Drag Pressure2.3.1 Grain Boundary Drag by Second Phase ParticlesSecond phase particles and solute atoms have been used as importantconstituents to design materials and processes due to their ability to controlthe motion of grain boundaries [1, 2]. As an example, in high strengthlow alloy (HSLA) steels, the addition of small amounts of micro-alloyingelements creates fine carbide precipitates [65]. These particles limit themotion of grain boundaries by the particle pinning mechanism. Particlesexert an opposing force on the interface motion. This drag force reducesthe velocity of the interface or, if it exceeds the driving pressure, it cancompletely pin the interface.The first theory of particle pinning was proposed by Zener and Smithin 1948 [5]. It consist of three steps: (i) describing the particle-interfaceinteraction and calculating the drag force of an individual particle on theinterface, (ii) establishing the cumulative effect of a dispersion of particles onthe interface and finding the pinning pressure (iii) explaining the effects of19the pinning pressure on the grain structure and its evolution kinetics. Latertheories were built on the same principles with a number of refinements ofthe Zener-Smith theory.Fig. 2.6 shows an interface in contact with a particle with radius rp. Theforce on the interface can be obtained by balancing the interface tension forcewith the drag force, assuming a stationary particle. When the angle betweeninterface and particle is ? = pi/4, the force is at maximum[2], i.e:FZ = pirp? (2.28)Here, FZ is the detachment force or the force required to push theinterface away from the particle. Ashby [66] obtained a formula similarto Eq. 2.28 by considering the energy of the system before and afterdetachment. When the grain boundary is in contact with the particle, anarea of the grain boundary is removed; therefore, detachment of the interfaceadds an amount of energy of 2pir2p? to the system. This is equal to the workdone by the pinning force of FZ = pirp? over the length of 2r.Figure 2.6: The balance of the interface tension force of a systemcomprising one spherical particle and a curved interface. Thedrag force reaches maximum at ? = pi/4.20The second part of the analysis is related to how the distribution ofparticles affects the interface motion. In the most simplistic case, one canassume a planar interface that intersects with a uniform distribution ofparticles. The number of particles per unit volume (the number density)is:Np =3fv4pir3p(2.29)where fv is the volume fraction of the particles. If the interface remainsplanar, it interacts with Np ? 2rp particles per unit area. Therefore, thedrag pressure from the particles, Pz, is given by:Pz = 2rpNp ? Fz =3fv?2rp= z?fvrp(2.30)The z = 3/2 parameter is known as the Zener factor. It is assumed thatif the drag pressure reaches the driving pressure for the interface movement,the interface then stops and the structure reaches a limiting grain size. Byassuming a spherical shape for a grain, the driving pressure is Pd = 2?/R?,where R? is the average grain radius. The limiting grain size is obtained as:R?lim =4rp3fv(2.31)This value for the limiting grain size is six to eight times larger than thatof the experimental observations [1, 10]. This is due to crude assumptionsmade regarding the interaction of the particle distribution with the grainboundaries and the estimation of grain growth driving pressure. Forinstance, the grain boundary does not remain planar during the interactionand thus the number of particles per unit area of grain boundary is higherthan 2rpNp. Several modifications have been made to Zener?s model.Hellman and Hillert [67] considered that the grain boundary assumes anaxisymmetric catenoid shape during anchorage and calculated the pinningpressure as a function of the position of the grain boundary. Fig. 2.7 shows21that the pinning force is a function of the radius of curvature of the grainboundary, ?. When the ratio of ?/rp increases, the grain boundary remainsin further contact with the particle, but the maximum detachment forceremains at FZ .Figure 2.7: The pinning force of a spherical particle on a grainboundary as a function of its position, x. The centre of theparticle is at x = 0 [6].Louat [24] and Hillert [6] assumed that particles in front of the movinginterface actually attract the interface because it reduces the energy of thesystem; i.e., the pinning force is applied only during the detachment. Inaddition, rather than assuming a planar interface, a more realistic shape forthe interface was considered [7, 68]. Such modifications lead to a similarformula for the limiting grain size:R?lim = ?rpf ?v, (2.32)where ? and ? are parameters depending on the underlying assumptions ofthe model. For example, Hillert [6] proposed ? = 4/9 and ? = 0.93 andWo?rner [7, 68] suggested ? = 0.31 and ? = 0.87. Fig. 2.8 summarizes someof the experimental results for the limiting grain size in metals containingsecond phase precipitates. For most systems with low particle volumefractions, the values of ? = 0.17 and ? = 1 [10] create a good fit.Particle pinning in 2D has very limited practical applications. Evenin polycrystalline thin films, second phase particles rarely occupy the entirethin film thickness. Nonetheless, corresponding relationships for the pinning22Figure 2.8: The ratio of the limiting grain size to the particle radiusversus the volume fraction of the second phase particles. Mostof the experimental data for systems with low volume fractionscoincide with the generalized formula (Eq. 2.32) with ? = 0.17and ? = 1 [10].pressure and the limiting grain size can be simply derived as [69]:Fz(2D) = 2?, (2.33)R?lim(2D) =pi2rpfa, (2.34)where fa is the area fraction of particles. The latter relationship for thelimiting grain size has not been observed in 2D computer simulations;instead, R?lim ? 1/?fa [64].To obtain the pinning pressure on a moving interface, one needs toaddress the complex geometrical and stochastic problem of the interfaceinteraction with an ensemble of particles. Furthermore, the influence ofpinning pressure on the evolution of grain structure requires an accurate23description of the grain growth driving pressure and the balance betweenthe driving and drag pressures. In particular, describing the kinetics ofgrain growth has important technological applications.2.3.2 The Kinetics of Grain Growth with Pinning PressureAlthough ideal grain growth has been studied extensively, the kinetics andtopological aspects of grain growth under the influence of drag pressuresis less known. Initially, Hillert introduced the drag pressure of particlesinto the kinetic relationship of a grain by modifying Eq. 2.5. The pinningpressure is included as a parameter that reduces the growth or shrinkagerate of individual grains [22]:dRdt= M[??(1Rcr?1R)? Pz](2.35)The sign of ? depends on the radius of each grain. It is negative forlarge grains when ??/R < ??/Rcr ? Pz and positive for small grains when??/R > ??/Rcr + Pz. Between this range dR/dt = 0. Assuming thatthe grain size distribution is self-similar, the following relationship for thegrowth of Rcr can be obtained [22]:dR2crdt=M2??(1?PzRcr??)2(2.36)Hunderi and Ryum [70] later solved numerically the continuity differen-tial equation for grain size distribution (Eq. 2.7). They used Hillert?s grainsize distribution (Eq. 2.8) as an initial condition. Fig. 2.9 shows the kineticsof grain growth and the normalized grain size distribution at various timesuntil the limiting grain size is achieved. The numerical results differ fromEq. 2.36 except for small pinning pressures and later stages of growth beforereaching the limiting grain size. In addition, when the pinning pressure islarge, the initial grain size distribution affects the value of the limiting grainsize and the final grain size distribution. The grain size distributions shownin Fig. 2.9b has not been observed experimentally.24(a) (b)Figure 2.9: (a) The kinetics of grain growth with the normalizedpinning pressure of Pz/? = 0.10. (b) The evolution at differenttimes of the normalized grain size distribution until the limitinggrain size is reached [70].To replicate the experimental data, Andersen and Grong [71] developeda phenomenological model in which the relationship between grain boundarycurvature and the average driving pressure for grain growth is assumed:P?d =?1?R?. (2.37)The constant, ?1, can be chosen to fit the experimental data. Forexample, based on measurements by Patterson and Liu [72], Rios [73]suggested that ?1 = 0.47. The following relationship was then proposedfor the average grain growth rate by subtracting the pinning pressure fromthe driving pressure:dR?dt= M(?1?R?? Pz). (2.38)This rate equation can be integrated from an initial grain size to obtainthe average grain radius as a function of time. The limiting grain size isreached when dR/dt = 0, i.e:R?lim =?1?Pz. (2.39)25Eq. 2.38 has been used to create consistent fits with experimental datafor systems with the presence of stable or dissolving precipitates [71, 74].However, there have always been adjustable parameters in the model, such as?1 or the Zener factor, z, in Eq. 2.30 to tune the model to fit the experimentaldata. An analytical model can than capture the evolution of the grainstructure and its kinetics has not yet been developed. The challenge in thisproblem is that the drag pressure adds yet another layer of complexity tothe already unresolved problem of ideal grain growth. The geometric andstatistical complexity of the grain growth with pinning can be addressed bysimulating the microstructural evolution.2.3.3 Computer Simulations of Particle PinningA variety of microstructural evolution models have been utilized to simulatethe effect of second phase particles on grain growth. Initially, Srolovitz etal. [75] proposed a Monte Carlo simulation method in which inert particleswere considered as areas with a specific index and with zero probabilityfor changing to other grain indices [75]. This creates an obstacle for grainboundary motion in which the pinning force is a function of the simulationtemperature. With low thermal activation, grain boundaries can not detachfrom particles. In 2D simulations, values of ? = 1.18 and ? = 0.52 (Eq. 2.32)were obtained [76]. Weygand et al. [64] extended their 2D vertex model toinclude particles with a designated detachment force and predicted ? = 0.5similar to the Monte Carlo simulations. They have not reported a value for?.For 2D systems, Moelans et al. [77] obtained a similar result (? = 1.28and ? = 0.5) using a phase field model. An example of the simulationsnapshots is shown in Fig. 2.10 where a dispersion of particles with a volumefraction of 0.015 is considered. Grain boundaries interact with particles andif their curvature is sufficiently large they can pass the particles. As grainsgrow and the curvature of the grain boundary segments decreases, the grainboundaries are pinned by the particles until none of the segments can move.In 2D simulations, the limiting grain size is much smaller than that26Figure 2.10: Microstructural evolution in a 2D polycrystalline systemcontaining second-phase particles for rp = 3 and fv = 0.015.Snapshots at time steps (a) 200, (b) 3000 and (c) 30,000 areshown. No further evolution was observed after approximately30,000 time steps, where Rlim = 29.5 [77].of the Zener model prediction (Eq. 2.34). As can be seen in Fig. 2.10,a large fraction of particles is located on the grain boundaries and triplepoints. This positioning of particles at the junctions creates a much largerpinning pressure compared to that of a random distribution. In addition,the interaction of a particle with a grain boundary segment can pin theentire segment by making it flat and eliminating its driving pressure. Basedon the observation that three particles are sufficient to pin all segments ofa grain, Srolovitz et al. [47] suggested ? =?3 and ? = 0.5 for the limitinggrain size in 2D which is more consistent with the computer simulations.This relationship can be further refined by considering the fraction of theparticles that are not positioned on grain boundaries and the particles thatare located at triple points [6].A simulation of particle pinning in 3D requires much larger computa-tional resources. In a large-scale Monte Carlo simulation, Miodownik etal. [11] obtained values of ? = 0.728 and ? = 1.02 for a range of volumefractions 0.025 < fv < 0.15. A simulated structure is shown in Fig. 2.11a atlimiting grain size. The ? value is much higher than the experimental valueof ? = 0.17. In contrast to experimental observations (Fig. 2.8), ? does notchange when the particle fraction increases above 0.05.27Suwa et al. [12] used the phase field model proposed by Moelans etal. [9] to simulate grain growth in 3D systems with 0.04 < fv < 0.12.They obtained values of ? = 1.42 and ? = 0.870 that differs from those ofthe Monte Carlo simulations and are not compatible with the experiments,either.Figure 2.11: 3D simulations of grain growth with particle pinning.(a) A cross-section of a structure modelled with the MonteCarlo method (fv = 0.025 with particles as cubes of 3? 3 gridpoints). The circle exhibits the dimpling of the grain boundaryin contact with the particles [11]. (b) Phase field simulationwith fv = 0.04 and particle size of 9? 9 grid points [12].Overall, the results of the microstructural modelling of grain growthwith pinning have not been consistent with experimental data despite theircomplexity and computational cost. A fundamental source of discrepancy isthe difference in the length scale of the grain boundary-particle interactioncompared to that of the grain size. In most materials, particle size is inthe range of 10 nm to 500 nm whereas the grain size is 1 ?m to 100 ?m. It isthus unfeasible to simulate microstructural evolution with such variations inscale. The multi-scale nature of the particle pinning process thus requires amodel that can resolve grain boundary-particle interaction on the nano scaleand yet be large enough to include hundreds of grains on the meso-scale.As an alternative, first, the interaction of particles with a grain boundary28can be modelled on a small scale to obtain the pinning pressure. Then,the pinning pressure can be included into another meso-scale grain growthmodel to study the kinetics of grain growth.Couturier et al. [8, 78] developed a finite element model for theinteraction of a moving grain boundary segment with an ensemble ofparticles. In this simulation, the grain boundary is modelled analogousto an elastic membrane that can be merged with and detached from thestationary particles. Fig. 2.12 shows snapshots of the system?s setup.Figure 2.12: The simulation of a grain boundary segment interactionwith an ensemble of particles using finite element method [8].The model consisted of a canonical domain in which a grain boundarymoves from the top to the bottom of a domain while concurrently reducingits curvature and driving pressure. This geometry can be transformed toa cylindrical system, as shown in Fig. 2.12, where the boundary maintainsits curvature but the particles on the top become larger and more distantcompared with those on the bottom. This allowed the authors to virtuallysimulate the process of the growth of a grain and its final stagnation ina single run. The particles radii were 0.005 times the radius of the grainboundary curvature and fv = 5? 10?4.To determine the pinning pressure, the velocity of the grain boundarywas measured. Fig. 2.13a shows the grain boundary velocity as a function29of the grain boundary position for a system with one particle. First,the grain boundary is attracted to the particle and the velocity increases.This is followed by the anchoring of the grain boundary with the particle,which slows down its motion over a longer distance. The interaction ofmany particles with a grain boundary can be interpreted as a randomsuperposition of many particle-grain boundary interactions and is shownin Fig. 2.13b. For the given particle distribution, the motion of the grainboundary finally halted when its declining driving pressure finally matchedthe pinning pressure of the particle distribution. The pinning pressure atthe time of the anchoring was given as a function of particle size and volumefraction [8]:Pz = 3.3?fvrp(2.40)Figure 2.13: The velocity of a grain boundary interacting (a) with aparticle and (b) with a distribution of particles modelled withthe finite element method [8].Compared to Eq. 2.30, this pinning pressure is 2.2 times larger than thatof the pinning pressure assumed by Zener. A higher number of particlesinteract with a grain boundary with a flexible geometry than with onepossessing a flat surface. Furthermore, not all the particles in contact withthe grain boundary at a moment exert the maximum force on the grainboundary. These circumstances are both accurately captured by the finite30element model. In addition, due to the sharp nature of the grain boundary,this model is more suitable for taking into account the difference in scalebetween the interface thickness (? 1 nm) and the particle diameter comparedwith those of the phase field models. However, a large scale simulation ofgrain growth is impractical with this technique due to its computationalexpense [8].To include the effect of particle pinning in a meso-scale model withoutresolving the fine particles, one can employ a conventional simulationtechnique for grain growth but use a reduced (effective) mobility for the grainboundaries. Apel et al. [79] used a phase field model to simulate interactionof a planar grain boundary pushed through an ensemble of particles withdriving pressure Pd. This method requires a smaller simulation domaincompared to a geometry with a curved grain boundary. Fig. 2.14a shows thevelocity of the grain boundary as a function of Pd for systems having differentvolume fractions of particles. Intrinsic mobility of the grain boundarywas M =5? 10?3 cm4/Js. With the presence of particles, the velocityof the grain boundary decreases. It has a linear relationship with Pd athigh velocity regimes as predicted by Zener?s theory (v = M(Pd ? Pz)).However, at low velocities, a driving pressure dependent reduction invelocity was observed. The effective mobility can be deduced by fittinga phenomenological relationship on the velocity as a function of the drivingpressure [79]:Meff = M exp(?0.12PzPd ? Pz), (2.41)such that v = MeffPd. The mobility given by Eq. 2.41 was used tomodel grain growth on the meso-scale starting from different initial grainsizes and with the presence of different pinning pressures. Three regionswere identified: (i) grain growth reaches a stable limiting grain size, (ii)large grains appear in the structure, i.e. abnormal grain growth and (iii)no significant grain growth occurs since the driving pressure of the graingrowth is smaller than the pinning pressure (Fig. 2.14b). Unfortunately,31the authors do not explain the non-linear behaviour of velocity versus thedriving pressure as shown in Fig. 2.14a. It is not expected that a systemwith stable particles would have a velocity dependent pinning pressure orabnormally large grains would appear.Figure 2.14: (a)The velocity of a grain boundary interacting with anensemble of particles with a radius of 220 nm and differentvolume fractions. The dotted line represents a system with noparticles [79]. (b) A 2D meso-scale simulation of grain growthfor different pinning pressures and initial grain sizes [79].This problem can be investigated using the vertex models where thedriving pressure on each segment is calculated according to the localcurvature at each time step. This allows the pinning pressure to be directlysubtracted during the calculations. In a recent 2D simulation of grain growthwith pinning pressure, limiting grain sizes without any abnormal growthwere observed [80]. The grain size distribution changes from the self-similarstate and widens as the structure approaches the limiting grain size.Recent studies to devise a multi-scale approach to model grain growthwith particle pinning are promising. Nevertheless, a simulation resultconsistent with experimental data in 3D has not yet been reported.Constructing an accurate simulation technique for grain growth with particlepinning paves the way for studying microstructural evolution in importanttechnological materials and processes.322.4 Grain Boundary Drag by Solute AtomsThe atomic arrangement of a crystal is disrupted at the grain boundaryregions. Therefore, solute atoms have different energies when they are posi-tioned at grain boundaries. This energy difference creates a thermodynamicdriving force that causes solute atoms to segregate to grain boundaries (Fig.2.15). This concentration spike is the reason for solute drag [1]. When theinterface moves under a driving pressure, the solute atoms need to diffuseat the same speed in order to maintain their equilibrium concentrationprofile. Since the grain boundary moves and the solute spike follows, theconcentration profile becomes asymmetric and exerts an opposing pressure(i.e. solute drag pressure) on the grain boundary.Figure 2.15: The characteristics of a grain boundary interaction withsolute atoms. a) Energy of solute atoms, E, as a functionof distance from centre of an interface. Atoms have lowerenergy when they are within the grain boundary. b) Theinteraction force (F = ?dEdx ) between the solute atoms andthe interface. c) The equilibrium concentration profile of soluteatoms. d) The variation of the diffusion coefficient for soluteatoms around a grain boundary.Lu?cke and Detert [81] were the first to propose an energetics approachto describe this effect. The segregation level of atoms depends on their33interaction energy, E, with the grain boundary. The segregation can bedescribed by:? =4?2C0a20exp(?EkBT), (2.42)where ? is the number of solute atoms per unit area of the grain boundary,C0 is the bulk concentration of the solute atoms, a0 is the lattice constant,kB is the Boltzmann constant and T is the temperature.Cahn later [13] developed a widely used theory of solute drag. Hiscontinuum theory assumes an interaction energy profile and a continuouschange in diffusivity across the grain boundary. The diffusion equationcombined with the interaction energy can be solved to obtain a steadystate composition profile. When the grain boundary moves, the soluteconcentration profile changes from its equilibrium shape to an asymmetricalprofile, and creates a pressure in the opposite direction of the grain boundarymovement. Fig. 2.16a illustrates the concentration profile of solutes arounda grain boundary as a function of its velocity. At higher velocities, soluteatoms can not diffuse as fast as the movement of the grain boundary and,consequently, the concentration spike, as well as the solute drag pressuredecreases [13]. The change in the drag pressure as a function of the grainboundary velocity is shown schematically in Fig. 2.16b.Figure 2.16: The Effect of the grain boundary velocity on (a) soluteconcentration profile and (b) the consequent drag pressure.Cahn proposed the following approximate relationship for the drag34pressure, Ps in steady state condition [13]:Ps =av1 + bv2, (2.43)where v is the interface velocity and a and b are constants depending on thesolutes diffusivity and interaction energy:a = 4NvkBTC0? +???sinh2 E(x)2kBTD(x)dx, (2.44)a?b=NvkBT? +???(dEdx)2D(x) dx. (2.45)Here, x is a distance measured on an axis normal to the interface withan origin at the centre of the interface, D(x) and E(x) are the profilesfor the solutes? diffusivity and interaction energy with the interface andNv is the number of atoms per unit of volume. Since the drag pressureis velocity dependent, the interface assumes a velocity dependent on thedriving pressure that is obtained by solving v = M(Pd ? Ps). The velocityof the interface as a function of the driving pressure is schematically shownin Fig. 2.17. For a given driving pressure, the velocity decreases with theincrease in the solute concentration. At sufficiently high driving pressures,the interface velocity approaches the intrinsic mobility of the system withouta solute spike.At higher solute concentrations the curve becomes unstable and there isa jump from high to low interface velocities. There are two paths for thistransition. First, when the driving pressure decreases from a value in thehigh velocity regime, at a critical velocity, the solute spike can build up atthe interface and significantly reduce the velocity of the interface. Second,when the driving pressure increases from a value in the low velocity regime,the interface breaks up from the solute atoms at a critical velocity andtransition from the low to the high velocity regime occurs. The transitionvelocity is not the same for the two conditions as is shown in Fig. 2.17 by35Figure 2.17: A schematic representation for the variation of theinterface velocity versus the applied driving pressure forsystems with different solute concentrations.arrows.Hillert and Sundman [14] extended the solute drag theory further. Theyargued that the dissipation of Gibbs energy by the interface motion shouldbe equal to the same amount dissipated by the diffusion of solute atoms. Inthis theory, the treatment of phase interfaces is also possible and one canapply this model for phase transformations [82]. Most of the recent studieson solute drag is dedicated to improving Cahn?s assumptions. For instance,Mendelev and Srolovitz [83] expanded the theory by assuming that phaseshave the thermodynamic behaviour of non-ideal solutions. Bre?chet andPurdy [84] considered ternary systems with interaction coefficients betweenthe solutes and the resulting effects on the segregation of elements.The complexities of the analytical models for solute drag increasesignificantly when one considers a general non-steady state systems, seenin processes with rapid cooling and heating, and when dissolution orprecipitation occurs. Furthermore, the atomistic length scale of the interfacesegregation phenomenon limits the application of continuum models. On theother hand, the time scale of the problem is relatively long to be simulatedwith such atomistic models as those for molecular dynamics, since thisprocess involves the diffusion of atoms and the motion of the interface.36Simulation techniques, e.g. kinetic Monte Carlo [85] and phase fieldcrystal [86], that have atomistic length scales, but diffusive time scales, haverecently been devised for simulating the solute drag phenomenon. While theresults follows the overall trend of Cahn?s model, such simulations revealthe discreet nature of the interface-solute atoms interaction. For example,contrary to the continuum models, it was observed that the interfacedoes not maintain its flat shape during the migration, and atomistic scaleroughening occurs [85]. In addition, at low velocity regimes, the interfacerequires a critical driving pressure to overcome the initial barrier and beginto move [86].Phase field models have also been used to model solute drag. Since phasefield has a continuum nature, it can be considered a numerical proxy forCahn?s model. The interaction energy of the solute concentration spike withthe diffuse interface can be implemented in the phase field energy functional[87] (see Sec. 3.7). This model has been used for modelling 2D grain growthon the nano-scale [88, 89] and has demonstrated the non-steady state natureof grain boundary-solute spike interaction. For example, the attachment ofsome grain boundaries to the solute atmosphere and the detachment of a fewgrain boundary segments can lead to abnormal grain growth [88]. Simulatinga truly meso-scale microstructure is not feasible with this model since thescale of the concentration spike is orders of magnitude smaller than that ofthe meso-scale features, such as grain diameter.Setting aside the nano-scale interaction of atoms with the grain bound-ary, experimental results have shown that solutes not only slow down thekinetics of grain growth but also alter the parabolic grain growth law.Experimental data suggest a general relationship as [2]:R?? ? R??0 = kt, (2.46)where ? is a constant greater than 2, and up to 5 [2]. The effect of solutedrag on the kinetics of grain growth can be modelled by considering itsinfluence on the mobility of the grain boundaries. Strandlund et al. [90]37used a velocity dependent mobility function to simulate grain growth inbinary and ternary systems with a phase field model. The effective mobilityfunction was obtained by considering the energy dissipation required for thediffusion of solutes across the grain boundary [14, 91]. The grain growthexponent was observed to be ? = 2.0-3.5 for Fe-Ni systems and ? = 3.3-5.7for Ni-Cr-Fe systems, depending on the simulation temperature and soluteconcentration.One can reconcile the kinetics of grain growth observed experimentallywith the effective mobility parameter in a microstructural evolution model.However, this does not elucidate the physical process of solute drag, nordoes it explain the balance between the driving and drag pressures. Forany system with a different chemistry or temperature, a new experiment isrequired to ascertain the fitting parameters used in the effective mobility.This problem is more pronounced for systems with soluble second phaseparticles where solutes form intermetallic phases and therefore the effectof particle pinning and solute drag are combined and interdependent.Furthermore, if the solute concentration on the two sides of the interfaceis not the same, e.g. in phase transformations, the thermodynamic effectof the chemical driving pressure alters the solute segregation and the solutedrag pressure. To resolve these convolutions and quantitatively determinethe contribution of solute drag, further studies are required to increase ourunderstanding of the solute drag phenomenon both on the atomic scale andon its influence on the overall evolution of the microstructure.38Chapter 3Phase Field Models3.1 IntroductionThe phase field method is a simulation technique to model the motion ofinterfaces and the evolution of microstructures. In the phase field method,the constituents of a microstructure are modelled with continuum fieldsand interfaces are defined as regions where the field values varies. Thiscoarse graining of discreet aspects of materials, such as arrangements ofatoms in a lattice or interface, enables the phase field method to simulatephysical processes such as interface migration on the meso-scale and overlong diffusive times. The phase field model can be readily combinedwith thermodynamic properties of the system and diffusion equations tomodel the migration of interfaces. They also possess a methodologicaladvantage where the evolution of a structure is simulated by solvingdifferential equations for the fields, and the interfaces do not need to beexplicitly tracked. This is ideal for grain growth simulations where thetopological transitions of the grain boundary network is difficult to managein other models. The flexibility of the phase field models to allow forthe incorporation of such physical phenomena as diffusion, elasticity, andanisotropy gives a compelling reason to utilize and further develop them forother applications.393.2 A Historical OverviewDescribing the motion of phase boundaries has been a subject of interestsince the early nineteenth century. For instance, the speed at whicha layer of ice forms was discussed by Clapeyron in 1831 [92]. Suchproblems in mathematics are classified as heat transfer problems subjectto a moving boundary condition [93]. Later, Stefan proposed a solution,which is now known as the sharp interface solution to the moving boundaryproblems [94]. Solutions to these problems were also used to explain thekinetics of diffusional phase transformations. The sharp interface modelsof phase transformations rely on such macroscopic phenomena as heattransfer or diffusion in phases to describe the evolution of microstructures.Therefore, boundary conditions for the solution of differential equations aredefined explicitly at the interfaces. As a result, such models are incapableof predicting phenomena such as interface instability that occur duringsolidification and dendrite formation.An alternative view to Stefan problems is that the nature of the interfaceis diffuse. An interface occupies a region in which the atomic structuregradually changes from one phase to another. Although this representationwas not new to the early chemists who observed interface layers between twoimmiscible liquids, it was not until Van Der Waals [95] and Gibbs [96] whodiscussed thermodynamic aspects such as surface energy and adsorption inthe late 19th century. The diffuse approach to the properties of interfaceswas studied by other authors [97?99] during the early to mid 20th century,until a seminal paper by Cahn and Hilliard [100, 101], in 1958, broughtup a new thermodynamic concept based on the construction of energyfrom gradient of atomic concentration across an interface. This theory, forinstance, explained the surface energy and thickness of liquid nitrogen as afunction of temperature.Advances in other fields of condensed matter physics, such as theGinzburg-Landau theory of superconductivity [102], inspired the develop-ment of phase field models by introducing the concept of spin fields or orderparameters to describe states of the microstructure and its evolution. This40idea was later used by Allen and Cahn [103] to develop a model for theevolution of anti-phase domain boundaries in intermetallics.During the 1980?s, two schools of thought emerged. One was based onthe microscopic theory of Khachaturyan [104], in which the transition atthe interface structure is expressed as a field variable closely related to aphysical property of the system. For example, a variable representing thecrystallographic orientation smoothly changes from one grain to another ata grain boundary. The thickness of the interface in this model matches itsphysical value.The second approach employs a field variable that is phenomenologicallyrelated to the energy of the system. This functionality can be adjusted fordesired interface thickness and energy. This theory is based on the treatmentof phase transition and critical phenomena explained by Hohenberg andHalperin [105]. Although the second approach has phenomenologicalrestrictions, it allows for the numerical solution of structural evolution inrealistic time and length scales.In the last two decades, there has been significant advances in the theoryand application of phase field models. Notable theories include phase fieldmodels for solidification [106, 107], the multi-phase field models of Chenand Wang [108] and Steinbach et al. [109]. Phase field models have beenapplied to a wide range of phenomena such as solidification [110, 111],solid state phase transformations [112] grain growth and recrystallization[62], microstructural evolution under stress [113, 114] and martensitictransformations [115]. Interested readers are referred to references [116?119] for a further review of the phase field models.In the next sections of this chapter, a heuristic rather than a historicapproach is adopted in reviewing phase field models. Physical principlesbehind the models are discussed with a focus on the area of modeldevelopment that is of interest for understanding particle pinning and thepresence of drag pressure.413.3 Describing Microstructures in Phase FieldModelsOne can recognize the microstructure of a material by the variation ofits local properties in space. A microstructure is defined by the spatialdependence of such attributes as crystallographic structure and orientation,chemical composition, density and spin. In that regard, phase field modelsare closely related to the definition of the structure. A phase field modelis constructed by considering a set of continuous functions across themicrostructure representing structural features or atomic concentrations.In the phase field community, these functions are known as an ?orderparameters? (?).In many systems with second order phase transitions, e.g., spinodaldecompositions or order-disorder transitions, one scalar function (field)can convey dependence of concentration or spin. For instance, in orderedintermetallics, two regions exists with similar lattice structures but differentatomic arrangements. As an example, Ni3Al has an FCC lattice with anordering of Al atoms in two configurations, one at the corner of the unitcell and one at the centre of the two opposite faces. These two variationscan be represented with a scalar field with the value of 1 for one atomicarrangement and -1 for the other arrangement (a value of 0 can be perceivedas a disordered phase). A variation of the order parameter between 1 and-1 is considered an interface or an antiphase boundary. Fig. 3.1a shows aschematic representation of such a system with one order parameter.The microstructure of alloy systems that undergo first order phasetransition consists of two or more phases. In addition to the variation ofthe order parameter, the concentration field may also have to be considered.For example, in the solidification of an alloy, a solid dendrite is formed fromliquid and the concentration of solute atoms is also changing in the system.Such a structure (Fig. 3.1b) can be described as possessing two fields: onerepresenting its constituent phases and the other one for its concentration.Order parameters can be a conserved quantity (e.g. concentration) ornon-conserved quantity. For instance, order parameters are not conserved in42(a) x?(b)?C ?cs(eq)cl(eq)xFigure 3.1: A schematic description of microstructures with fieldvariables for (a) a system with variation in spin orcrystallographic ordering, (b) a two phase system with variationof composition and phase, e.g. the solidification of a dendrite.solidification of a liquid or in the grain growth as the liquid or some grainsshrink and disappear during the microstructure?s evolution.An example of a more complex structure is when the distinction betweena variety of regions is required. For example, in a polycrystalline structure,each grain has a specific crystallographic orientation. As shown in Fig. 3.2,each grain is assigned an order parameter which is 1 inside the grain andsmoothly changes across the grain boundary to 0 outside the grain. Withthis method one can distinguish each grain as it grows or shrink. Grainboundaries and junctions are naturally identified as regions where two ormore order parameters have a non-zero value.The advantage of describing microstructures with field variables is thatdifferent phases or microstructural constituents (e.g. grains) can be easilydetermined and there is no need to explicitly track the position of theinterfaces. An interface region can be identified from a variation in theorder parameters. Relating the order parameter fields to the thermody-namic quantities is straightforward. Differential equations that govern theevolution of order parameters can be derived from the thermodynamics andevolution kinetics of the system.43x? ?5?3 ?4?1 ?2Figure 3.2: A polycrystalline structure defined by a set of orderparameters, ?i, associated with each grain.3.4 The Phase Field Model Construction3.4.1 ThermodynamicsThe notion of a diffuse interface was first introduced for systems with amiscibility gap in their phase diagrams. In these systems, a homogeneoussingle phase (e.g. a liquid or solid solution) can separate into two phases bythe process of spinodal decomposition [120]. The Helmholtz free energy,f0, of such a system is shown in Fig. 3.3a. The free energy can beinterpreted as the local energy density (energy per unit volume) at eachpoint in the system. In other words, f0 at any concentration is theenergy of a unit volume regardless of any spatial variation in the soluteconcentration. A system with a composition between the two inflectionpoints of the energy curve instantaneously decomposes into two regionswith higher and lower concentrations. As decomposition continues towarda structure with two phases with equilibrium concentrations, regions withgradients of concentration appears as the diffuse interface between the two44phases. The interface region has higher local energy density than theequilibrium phases and this contributes to the interfacial energy associatedwith the boundary.Figure 3.3: The local energy density of a system that decomposes intotwo regions as a function of a) composition, b) order parameter.One can extend the concept of immiscibility in systems with spinodaldecomposition to systems with two states separated by an interface. Inan ordered intermetallic with two ordering modes, for instance, bothatomic arrangements have similar energies but the boundary between theregions (the antiphase boundaries) have higher energy. In this case, aphenomenological field variable, ?, can be assumed to relate the structuralstate to the local energy density. The energy of the system has two minimaat +?e and ??e that represent two equilibrium atomic orderings (Fig. 3.3b).The local energy density of the system is given by [103]:f0(?) = ?m12?2 +m24?4 (3.1)where m1 and m2 are phenomenological constants. This function has twominima at ? = ??m1/m2. For instance, if m1 = m2, the phase fieldparameter is +1 and -1 inside the equilibrium domains. The disorderedphase has energy associated with ? = 0. For a system that is separated45into two ordered regions with a diffuse boundary in between, energy densityof each point is given by f0(?). The local energy density is higher for theinterface region that has ? between the equilibrium points. However, this isnot the only contribution to the total free energy since the spatial variationof ordering increases the energy of the system to a greater level than that ofa uniform system. For a non-uniform system, the local energy density canbe expanded to second order terms by a Taylor series, such as:f(?,??,?2?, ? ? ? ) = f0(?) +???f???+?2??f??2?+ (??)212?2f?(??)2+ ? ? ? (3.2)The second term in the expansion should be zero because the value of theenergy does not depend on the direction of the coordinate axis (reversingthe direction of x axis makes ?? negative). The total free energy of thenon-uniform system, F , can be written as an integral over space:F =?[f0(?) + ?1?2? + ?2(??)2 + ? ? ?]d3r (3.3)where ?1, ?2, ? ? ? are the values of the corresponding derivatives in Eq. 3.2.The integral of ?2? can be converted to the integral of (??)2 by integrationby parts. Ignoring the higher derivative terms, the total energy of the non-uniform system becomes [100]:F =?[f0(?) + ?(??)2] d3r, (3.4)where? = ?d?1d?+ ?2 = ?[?2f????2?]0+[?2f?(??)2]0(3.5)Here, the subscript 0 represents the value of the derivatives at the Taylorseries expansion point or ?? = 0, i.e. in the uniform disordered system. ? isknown as the gradient energy coefficient and has a positive value. It indicatesthat the variation of the order parameter at the diffuse interface increases46the energy of the system. Thus the value of ? is linked to the interfaceenergy. Eq. 3.4 is a functional with respect to the field ?(r) described bythe position vector r.The above energetic description can be extended to systems withmultiple order parameters. For a polycrystalline material, each grainrequires its own order parameter (Fig. 3.2). Since each grain is stable, thelocal energy density function should have a minimum energy associated witheach grain. An extension of Eq. 3.1 to this system was proposed by Chenand Yang [108]. A simplified form of the function is given by [121]:f0(?1, ?2, ? ? ? , ?p) = m??p?i=1(??2i2+?4i4)+32p?i=1p?j>i?2i ?2j +14?? . (3.6)Here, p is the number of phase field parameters and m is a constant thatdetermines the depth of the energy wells. Similar to Eq.3.1, f0 is a multi-wellfunction with a = b = m and a cross interaction term ?2i ?2j . This functionhas 2p minima located at ? =(1,0,...,0), (-1,0,...,0), (0,1,0,...,0), etc. Fig.3.4a shows the energy density for a system with two phase field parameters(i.e. two neighbouring grains).The total energy of a multi-phase system can be derived similarly to thatof the energy functional for the system with one order parameter (Eq. 3.4):F =? [floc + ?p?i=1|~??i|2]d3r (3.7)The constants m and ? determine the thickness and energy of theinterfaces between the order parameters which in this formulation are thesame for all boundaries (i.e. for a system with isotropic interface properties).The equilibrium order parameter profiles across an interface can beobtained by minimizing the total energy for a system consisting of two orderparameters (Fig. 3.4b). The profile of ?1 and ?2 defines a path on the energysurface from a minimum point to an adjacent one, represented as a solid line47in Fig. 3.4.Figure 3.4: (a) Homogeneous energy density as a function of orderparameters for a system with two phase fields. The black solidline represents a path where two phase field parameters changefrom inside one grain to a point inside another grain. Theequilibrium ? profile normal to the interface as a function ofdistance is shown in (b).In one dimension the order parameter profile through an interface isgiven by [122]:?i =12[1? tanh(?2?mr)], (3.8)where the ? sign refers to two symmetrical order parameter profiles. Theparameter lgb represents the thickness of the interface and is shown inFig. 3.4b:lgb =?8?m(3.9)3.4.2 The Kinetics of EvolutionThe non-uniform system with the energy described above evolves bydecreasing its energy towards the local gradient of the energy landscape. It48can be postulated that the rate of evolution for an order parameter (??i/?t)is proportional to the amount of energy dissipated (?F/??) associated withthat evolution [103]. For the multi-phase system, this statement is describedas:??i(r, t)?t= ?L?F??i(r, t)(3.10)Here L is a kinetic constant that can be related to the interface mobility.A set of partial differential equations that govern the dynamics of the systemis obtained by substituting the total energy functional (Eq. 3.7) into theevolution equation (Eq. 3.10), i.e.??i?t= ?L??m???3i ? ?i + 3?ip?j 6=i?2j??? ??2?i?? . (3.11)Solving this set of differential equations numerically describes theevolution of the microstructure. Numerical approaches for solving phasefield PDE will be discussed in Sec. 5.2.Model parameters in this formulation are related to the physicalquantities. The interfacial energy, ? can be calculated by substitutingEq. 3.8 into Eq. 3.7 and integrating it across the grain boundary:? =? ???[mf0(?i, ?j) + 2?(??i)2] dx = 1/3?2m?. (3.12)It can be proven [123] that a curved interface with the radius of curvatureR moves with a speed proportional to ?/R. This resembles the assumptionfor a sharp interface that the speed of movement is v = M?/R, where M isthe interface mobility and for the phase field model [122]:M =3L2?2?m. (3.13)Another multi-phase field formulation was introduced by Steinbach etal. [109]. In this model, one can use the physical parameters directly. Thefollowing equation can be obtained using a double well potential for two49order parameters:??i?t=p?j 6=iM{?[(?i?2?j ? ?j?2?i)?18?2?i?j(?j ? ?i)]?6?Gij??i?j}(3.14)Here, ?Gij is the free energy difference of the phases (grains) i and jand ? is the interface thickness. This model has been implemented into thecommercial phase field modelling software MICRESS [124].3.5 The Simulation of Microstructural Evolutionin Polycrystalline MaterialsAn immediate application of phase field models is the simulation ofmicrostructural evolution in single phase polycrystalline materials. Forinstance, a large body of literature is dedicated to the simulation ofgrain growth. Initially, Fan and Chen solved a set of phase field partialdifferential equations (Eq. 3.11) in a 2D domain with periodic boundaryconditions. Order parameters were initialized by random fluctuations inthe range of ?0.001 over the entire domain. With time, larger regionswith ?i = 1 emerge from this noisy domain that resembles the structureof a polycrystalline material. The coarsening of these domains basicallysimulates the curvature driven grain growth process. They only used alimited number of order parameters since it is computationally costly tosolve one differential equation for each grain on the entire domain. Fig. 3.5shows simulation snapshots for two systems with 4 and 36 order parameters.To reduce the computational expense of the simulations, one can usethe same order parameter to resolve two grains that are far from eachother. However, as it is evident in Fig. 3.5, when two grains with the sameorder parameter come into contact, they coalesce and cause unrealistic grainshapes. Grain coalescence is more predominant in Fig. 3.5a-d where only4 order parameters were used. To eliminate this artifact, Krill and Chen[62] proposed a numerical method that swaps the order parameter of a grain50Figure 3.5: The evolution of grain structure simulated with 2D phasefield for (a)-(d) with 4 order parameters and for (e)-(f) with36 order parameters. The timesteps are for (a)200, (b)1000,(c)2000, (d)4000; (e)200, (f)1000, (g)4000 and (h)10000 [108].when it approaches a grain with the same order parameter. A 3D simulationof structures with ?grain orientation reassignment? is shown in Fig. 3.6.Large scale simulations of grain growth using the phase field methodhave been performed. Statistical results such as grain size distribution andthe topological aspects of grain growth are shown to be consistent withother simulation techniques and analytical models [59, 60]. For instance,Fig. 3.7 shows simulation results of Darvishi et al. [59] on the distributionof the grain?s face number and the growth rate of the individual grainsas a function of their topological class. Phase field simulation results areconsistent with the large scale front tracking simulation performed by Lazaret al. [17] and the analytical models of Mullins [37] and Hilgenfeldt et al.[38]. In addition, it has been shown that the structures obtained from twophase field formulations (Eq. 3.11 and Eq. 3.14) coincide when starting froma same initial condition [57].The presented phase field method can also be used to model recrys-tallization, considering the uniform and non-uniform distribution of storedenergy [125, 126], the sub-grain boundary motion [127] and dynamic51Figure 3.6: The grain structures from a phase field model in 3D.Snapshots are taken at different times. The value of N indicatesthe number of grains in the simulation box. At time 10, onlythe initial random value of noise is shown, represented by smallnuclei in a space with all ?i = 0 [62].recrystallization [128]. In addition, the anisotropy of grain boundaryproperties has been included in the phase field model by making the modelparameters ? and m orientation dependent [122, 129, 130].52Figure 3.7: The topological properties of the 3D ideal grain structureis self-similar as simulated by a phase field model [59]. (a)The distribution of the number of faces of grains at differenttimesteps is compared with front tracking simulations by Lazaret al. [17]. (b)The normalized growth rate of grains as a functionof number of faces is compared to analytical models of Mullins(Eq. 2.17) and Hilgenfeldt et al. (Eq. 2.21).3.6 The Phase Field Modelling of ParticlePinningSince ideal grain growth rarely occurs in technologically relevant materials, itis important to incorporate second phase particles into phase field models inorder to simulate structural evolution in materials with major applications.Inert particles can be defined in the model with a separate stationary orderparameter. The particle order parameter (?) is 1 within the particle regionand 0 within the matrix. Since particles remain stationary during graingrowth, no evolution equation is required to be solved for them. However onecan manually change the ? at certain timesteps to mimic the microstructuralevolution of particles. Moelans et al. [9, 131] proposed an energy functionalfor a system containing particles:F =? [mf0 + ?p?1?2i +p?1?2(??i)2]d3r, (3.15)53where f0 is given in Eq. 3.6 and is a phenomenological constant thatdetermines the particle-matrix interface energy. Using Eq. 3.10, the partialdifferential equation for a system with inert particles is given by:??i?t= ?L??m???3i ? ?i + 3?ip?j 6=i?2j??+ 2??i ? ??2?i?? . (3.16)The dispersion of all particles of a given radius and volume fractionis specified by placing one ? field over the entire domain. Fig. 3.8 shows across-section, in 3D, of the interaction of six particles with a curved interface.The anchoring of the particle is evident in Fig. 3.8b. In this geometry, thecurvature of the grain creates enough driving pressure to detach the interfacefrom the particles.Figure 3.8: Cross-sectional snapshots of the interaction of a sphericalgrain boundary with six particles at different times [9].As proposed by Ashby [132], the pinning force of a particle can beanalyzed by comparing the energy of the system before and during theinterface-particle interaction (Fig. 3.9a-b). The energy difference, associatedwith the removal of some area of the interface, is shown in Fig. 3.9c. Fora 3D system, this energy should be 2pir2?, where r is the particle radius.However, the phase field model persistently predicts a higher decreases inthe energy [9]. This behaviour is associated with the unrealistic thicknessof the interface in the model. As shown in Fig. 3.9b, the thick interfaceremoves some of the particle-matrix interface. In reality, the thickness of54the interface is much smaller (compared with that of the particle diameter)than it is modelled in the phase field.Figure 3.9: Interaction of a grain boundary with a particle and thechange in total energy of the system. A flat interface andparticle in a phase field model shown in two configurations. Thedifference between (a) and (b) in the total energy of the systemis shown in (c) [9].When the interaction of one particle with a grain boundary is concerned,it has been shown [133] that the phase field produces comparable resultsto other simulation techniques such as the Monte Carlo or front trackingmethods for the pinning force and detachment angle. However, the interfaceremains attached to the particle for a longer distance and the pinning forceis slightly higher. Both artifacts are due to the high interface thickness inthe phase field models. Reducing the grain boundary thickness requires afiner numerical mesh and thus increases the computational cost.Phase field models can be used to study particle pinning in systemswith a coherent particle-matrix interface [134] and non-spherical (elongated)particle shapes [135]. These factors are shown to increase the pinning forceof a particle on an interface and thus create a stronger pressure for a givenparticle volume fraction compared with that of the Zener theory.Recent model development studies have been dedicated to capturing theevolution of particles. In many alloys, particle coarsening occurs during theisothermal holding of the material at high temperatures. In these conditions,55it was shown that the parabolic grain growth kinetics is replaced by a kineticswith rate constant of 3, coinciding with the coarsening time exponent [136].Small particles can also move under the force excreted by the moving grainboundary [132]. This process is facilitated by the diffusion of atoms at theparticle-matrix interface. To model the mobile particle drag, the particleswere allowed to move at a speed proportional to the force exerted by theinterface using a phenomenological parameter interpreted as the particlemobility constant [137]. It was observed that the pinning pressure is reducedas particle mobility increases.Phase field models are flexible tools employed to simulate particlepinning due to their ability to include various geometrical and physicalproperties. However, the interaction of an interface with a thicknesscomparable to a particle?s diameter produces artifacts. For this reason,the application of current phase field models is limited to structures ofsmall grain size to particle size ratio (e.g. nano crystalline materials). Thisproblem exists in many technologically important systems with grain sizesat the micron scale and fine particles at the nano scale.3.7 Phase Field Modelling of Solute DragThe diffuse nature of interfaces in phase field models makes them idealfor capturing variation in the interface physical properties that causes thesegregation of solutes at the interfacial region. Different techniques areutilized to incorporate this effect into the phase field simulations. In thefirst approach, one can define the interface as another phase with a separateorder parameter [138]. Secondly, in a more robust way, one can change theproperties of the system spatially, depending on the position of the lattice orthe grain boundary. Ordinary phase field formulations can be used with aneffective (reduced) mobility for the interface to mimic the solute drag [90].The first attempt to include the interaction of solutes and an interfacewas perfomed with the simulation of solute trapping in solidification [139,140]. Wheeler et al. used the phase field model of Wheeler-Boettinger-McFadden [106] to show that, for high solidification rates, solute atoms56begin to segregate in the front of the moving interface, causing a decreasein the driving force for the solid-liquid interface motion.For grain boundaries, Fan and Chen [141] modified their previous phasefield model [108, 142] by introducing a concentration dependent energyfunctional. They showed that, in the low velocity regime, solute dragincreases the time-dependent exponent of grain growth to values greaterthan 2. Similar to the approach of Fan and Chen, Gro?nhagen and A?gren[87] modified the height of the double-well energy density curve (Fig. 3.3b)as a function of local concentration. The total energy functional of thesystem is written as:F =?[f0(?, CB) + ?(??)2] d3r, (3.17)where f0(?, CB) can be phenomenologically related to the concentration.For instance, a linear relationship between the height of the energy curveand the concentration yields:f0 = ?2(1? ?2)(WBXB +WA(1?XB)). (3.18)Here, WA and WB are constants that can be determined from the grainboundary energy and the segregation level at the grain boundaries. In thismodel, increasing the solute concentration decreases the local energy densityof the grain boundary region. As a consequence, the solute segregates tothe grain boundaries and the total energy of the system decreases. Whenthe grain boundary moves, the solute atoms tend to diffuse as well. Thus,the diffusion equation needs to be solved simultaneously with the phase fieldequation. The diffusion equation coupled with the order parameter can bewritten as [87]:dCBdt= ?(MB(1? CB)CB(?2f0?C2B?CB +?2f0?CB????))(3.19)Therefore, by solving Eqs. 3.10 and 3.19, the solute drag at the grainboundaries can be modelled. Fig. 3.10 shows concentration profiles for threedifferent interface velocities. As expected, by increasing the velocity of the57Figure 3.10: The concentration spike at an interface for three differentinterface velocities. From left to right, the interface migratesmore quickly. For higher velocities, the concentration at theinterface decreases and, consequently, the solute drag force alsodecreases [87].interface the solute atoms can not diffuse as quickly, and the concentrationspike and therefore the solute drag force decreases.With the phase field method, one can model variable diffusivity acrossthe interface and the non-steady state conditions [89]. Current models are inthe initial stages of development and consider only single phase and binarysystems.3.8 Conclusions of the Literature ReviewIdeal grain growth is a very well studied subject in materials science.Although a comprehensive analytical theory has not hitherto been devel-oped, numerical simulations have assisted the theoretical investigations bycapturing complex geometrical and statistical aspects of grain structure. Onthe other hand, grain growth with the presence of obstacles has been lessunderstood. Zener?s theory and its subsequent modifications provide a roughapproximation for the value of the limiting grain size, while Cahn?s theoryprovides an outline for solute drag pressure. However, major questions aboutthe kinetics of grain growth and its geometrical and topological aspectsremain unanswered. For instance, the experimental results for the limitinggrain size are not compatible with those of analytical models and in manyapplications, kinetic models with fitting parameters are utilized to replicatethe observed kinetics.58Computer simulation techniques such as Monte Carlo and Phase Fieldhave been developed to model grain growth with the presence of adistribution of particles. However, the results of such simulations can notbe applied to technologically important materials where there are orders ofmagnitude of difference between the length scale of the particles and thegrain size.Phase field models have been chosen in many grain growth studiesdue to their ability to capture with ease the topological transitions of thegrain boundary network. With the recent developments, one can obtaindrag pressure by simulating the interactions of a grain boundary with theparticles or the solute concentration spike. These simulations are capableof capturing the dynamics of the drag phenomenon in nano-scale whereparticles and solute segregation spike are resolved. However, they can notprovide appropriate representation of grain growth dynamics in meso-scale.59Chapter 4Scope and ObjectivesThe extensive studies of the inhibition of grain growth by particles andsolutes over the past seven decades have not been able to accurately captureimportant aspects of grain growth kinetics and the limiting grain size. Themajor source of complication has been the gap between the scale of theinteraction of the grain boundaries with obstacles and the scale of graingrowth.The goal of this project is to develop a multi-scale approach to bridgethe gap between the length scale of particle/solute interactions with thegrain boundary and the length scale of the polycrystalline structure. Forthis purpose, the phase field modelling technique is used in two scales.On the nano scale, the pinning pressure of a particle distribution with agiven volume fraction and particle size is determined. On the next stepan opposing driving pressure model is developed to account for the pinningpressure on the migrating boundaries in meso-scale without resolving theparticles explicitly.The friction pressure model, developed in this project, will be used in 2Dand 3D grain growth simulations with the presence of a drag pressure. Forthe particle pinning, the drag pressure is a constant pinning pressure. Forthe case of solute drag, the drag pressure is velocity dependent. Simulationresults is then analyzed and compared with the available experimentalobservations.60Chapter 5The Methodology andNumerical Techniques5.1 IntroductionSolving the multi-phase field equations on the large scale is computationallyintensive, making the development of efficient numerical procedures critical.To simulate the evolution of a grain structure, hundreds of non-lineardifferential equations, one for each order parameter, need to be solved. Asa consequence, a variety of numerical optimization techniques have beendeveloped to make such calculations feasible.Post-processing algorithms play an important role in obtaining accurateand meaningful results. During simulations of the grain structure, snapshotsof the simulation domain containing order parameters are stored on a disk.The phase field data is then analyzed to obtain images of the structure, aswell as statistical information about its grain size, grain growth rate andtopological features. A novel method is also developed for measuring thecurvature of grain boundaries.The finite difference method has been widely employed for solving multi-phase field equations [62, 129, 143] in which the forward explicit Eulermethod was applied for numerical differentiations. Despite the fact thatusing the implicit Euler method increases the time step size, it has only61a small advantage in reducing the computational cost, since the requirediterations for convergence renders the larger time step size ineffective.5.2 The Numerical Solution of the Phase FieldEquation5.2.1 The Finite Difference MethodThe multi-phase field equations (Eq. 3.11) are highly non-linear partialdifferential equations (PDE) that are coupled at interfaces and junctions.Since the value of d?i/dt is explicitly defined by the phase field PDE at timestep tn, the value of ?i at the next time step is given by the forward Eulermethod:?i(tn+1) = ?i(tn) +??i?t(tn)?t, (5.1)where ?t is the duration of each time step. To find d?i/dt, the right handside of the phase field equation (Eq. 3.11) should be evaluated numerically.This requires the spatial discretization of the domain into grid points wherethe non-linear term and the Laplacian can be evaluated at each point. Ahigh order 9-point stencil in 2D and a 19-point stencil in 3D is used [144].Fig. 5.1 shows the configuration of the stencils in the 2D and 3D squaregrids. Using both the first and second nearest neighbours in the high orderstencils increases the accuracy and reduces the numerical anisotropy of thecalculations.The Laplacian in 2D is calculated using the 9-point stencil as [144]:62Figure 5.1: (a) A 2D 9-point stencil (b) A 3D 19-point stencil forcalculation of the Laplacian term in the phase field equation.?2?i(r) =?2?i?2x+?2?i?2y=16(?x)2??4?r+?x?i(r+?x) +?r+?2?x?i(r+?2?x)? 20?i(r)?? ,(5.2)and in 3D, with:?2?(r) =?2??2x+?2??2y+?2??2z=16(?x)2??2?r+?x?i(r+?x) +?r+?2?x?i(r+?2?x)? 24?i(r)?? . (5.3)Here, r+?x and r+?2?x refer to the position of the first and secondnearest neighbours of the grid point positioned at r. The summations areover all the neighbours shown in Fig. 5.1.The coupling term (?i??2j ) in Eq. 3.11 is the only term that requires63the value of other phase order parameters at the point located at r. Otherterms can be simply calculated by knowing ?i(r). The computing processfor each time step is performed by nested loops sweeping through all thegrid points. At each grid point the following procedure is performed:i The value of ?i and the value of its first and second nearest neighboursare taken from the data structure and sent to the finite differencingalgorithm.ii The Laplacian term is calculated using Eq. 5.2 or Eq. 5.3,iii The non-linear and coupling terms are calculated,iv The value of ?i(tn+1) is calculated by Eq. 5.1 and is stored in anauxiliary variable.This process continues for all order parameters until the values of ?i(tn+1)are calculated for the entire domain. This is the end of the time step tn.The values of ?i(tn+1) are replaced with ?i(tn) and the next time step isstarted.5.2.2 Quantitative Phase Field ModelA solution of the phase field differential equations is acceptable when itmatches that of the curvature driven motion of a grain boundary. It hasbeen shown that the so called sharp interface condition is achieved for apair of order parameters at an interface using Eq. 3.11 [122, 123]. To have aproper numerical accuracy, the diffuse interface should be discretized withsufficient grid points along its thickness. The accuracy of the simulation isalso subject to the ratio of the interface curvature to the interface thickness.The average grain size should be at least 5 times larger than the thicknessof the grain boundary [122, 123].The parameters of the phase field model can be adjusted to representa particular physical system. The interface thickness, energy and mobilityare given in Eqs. 3.9, 3.12 and 3.13, respectively, as a function of the modelparameters ? and m and L. Each parameter, along with the numericalparameters ?x and ?t, thus assumes the physical units of the system.For instance, a simulation with m =1 J m?3 and ? =2 J m?3 has the64grain boundary energy of ? =2/3 Jm?2. In this system, a curved grainboundary with a radius of curvature of ? =100 ?m has the driving pressureof Pd =0.01 MJ m?3= 10?14 J(?m)?3. To simulate the motion of the givengrain boundary under this driving pressure, ?x =1 ?m and ?t =1 ?s canbe used. By considering L =104 (?s)?1 J?1 which is equivalent to themobility of M =3?104 (?m)4 (?s)?1 J?1, the grain boundary velocity ofv =3 ?10?4 ?ms?1 is observed in the sharp interface as well as in the phasefield model. Other physical properties with their corresponding units canbe similarly derived.Therefore, the choice of the system of units depends on the applicationand follows the modelled physical system. In this thesis the simulations areperformed on model systems and the results are reported as dimensionlessquantities were possible. In cases where the results are not normalized, theunits are omitted for the sake of brevity and the fact that they are notconnected to any particular physical system. Nonetheless, the list of unitsfor all parameters and quantities are given in the Nomenclature of the thesis.These model parameters in this work have been chosen in a way to beproduce reasonably close results in comparison to the sharp interface model.In a numerical test, a circular or spherical grain with a radius of 150 gridpoints shrinks with a driving pressure due to its curvature. The velocityof the sharp interface is v = M?/?. The velocity of the diffuse interfaceis calculated in the test from the area/volume of the circle/sphere. In allsimulations, the phase field and numerical parameters have been chosen toproduce less than 3% error in the velocity of the interface for a circle/sphereradius greater than 5. Below this value, the grain shrinks faster than thesharp interface benchmark and disappears. However, the size of the grain isvery small compared with the simulation domain and its effect on the overallmicrostructure is negligible.In a system with a few order parameters, it is possible to store the valueof each order parameter over the entire domain in the computer memory andsolve Eq. 3.11 at the interface region where 1?10?5 < ? < (1?1?10?5). Forinstance, in simulations where a grain boundary interacts with an ensembleof particles, only two order parameters are required in a domain with the65symmetric boundary condition. However, in meso-scale simulations wherehundreds of grains are present, a special numerical technique is utilized toreduce the amount of memory and processing power required to solve thephase field equations.5.2.3 Numerical Optimization TechniquesTheoretically, the finite difference process should be performed nx?ny?nz?p times, where nx, ny and nz are the number of grid points in each directionof the system and p is the number of order parameters, each representingone grain in the structure. This makes it computationally unfeasible tomodel large scale structures with the phase field method. The calculationspeed, however, can be considerably enhanced by numerical optimizationtechniques.There is the possibility of reducing the number of grid points in thesystem. Order parameters are constant within a phase and only vary acrossthe interfaces. Therefore, there exists only the need to resolve the interfaceswith a fine mesh and to discretize within the phases with a coarse grid.Such adaptive meshing techniques have been used in modelling solidificationand solid state dendrites [111, 114, 145, 146]. There is, however, anoverhead computational cost for maintaining and re-structuring the adaptivemesh during the calculations. This renders this method less favourable forstructures with large densities of interfaces, e.g. grain structures.A grain occupies a small volume of the entire domain and only withinthat area does the order parameter have a non-zero value. Therefore, asparse data structure can be used to reduce the number of order parametersthat need to be stored at any point. The number of order parameters of non-zero value at any point depends on the topological arrangement of the grains.The maximum number in 2D occurs when a four-sided grain shrinks anddisappears. In this situation, five order parameters are required to capturethe interactions of the grains with diffuse grain boundaries. In 3D, a smallshrinking grain usually has a tetragonal shape. With much less frequency,grains with six to seven sides shrink and disappear. Therefore, eight order66parameters are sufficient to account for all topological transitions. Thisallows the merging of up to two quadruple points to be considered withouterror.The data structure consists of two arrays with the size nx?ny?nz? pnwhere pn = 8 for 3D systems and pn = 5 for 2D systems (with nz = 1). Onearray, eta, stores the value of order parameters ?i and another array, ind,stores unique indices of the grains corresponding to the order parameter i.Fig. 5.2 shows a schematic representation of the data structure for threegrains. Each layer of the eta and ind arrays are a space holder to store apair of index i and the value of ?i for any given point in the grid. If the pointcontains two non-zero order parameters (i.e. at the interface), the secondlayer is used to store the additional pair and, at triple junctions, three layersare used. A maximum of pn pairs can exist at a point to accommodate ajunction where all pn order parameters are not equal to zero.This is an unstructured data set and there is no need to sort the orderparameters from high to low at each time step. To calculate the Laplacianof ?i at any grid point, it is required to search for neighbours with index iand to send the sum of first the and second nearest neighbours to the ?2processing algorithm (Eq. 5.3). At the end of each time step, the list ofactive order parameters is maintained in the ind array. If ?i > 10?5 at anypoint, the values of other order parameters at that point are checked andare replaced with the pair of ?i and i, if ?j has a smaller value than that ofthe threshold. The algorithm also adds the index of i to the neighbours inthe ind array, so that there is always a padding of one grid point ahead ofthe growing order parameter.The initial grain structure is created by setting 0.5% of the grid points inthe sparse data structure with the value of 1, each associated with an orderparameter ?i. These seeds are positioned randomly in a spatially uniformdistribution. The remainder of the unselected grid points are initializedto 0. The seeds grow due to the higher energy of the regions outside ?iand subsequently impinge on each other, creating a fine grain structure.Subsequently, this structure is allowed to coarsen until a self-similar grainstructure, similar to that illustrated in Fig. 3.6, is obtained.67Figure 5.2: Schematic representation of the sparce data structure. Anarea of the structure having three grains is discretized with agrid and the values of the order parameters and the grain indicesare stored in a sparse data structure.5.3 Post-Processing Algorithms5.3.1 Visualization of the StructureA phase field simulation generates large data sets. For example, for a systemwith a size of 300?300?300?8, the order parameter array eta with doubleprecision floating point variables requires 1.61 GB of memory and the indarray with integer precision variable takes 820 MB of storage. This dataneeds to be processed further to obtain meaningful information about thestructures.The grain structure can be revealed as an image with two methods.68First, the position of the grain boundaries can be displayed by creating agrayscale contrast from within the grain to the core of the grain boundary.The grayscale image is created by:? =p?i=1?2i =pn?n=1(eta(:,n))2 . (5.4)Here, eta(:,n) refers to the n-th layer of the eta array. Inside thegrains, the value of ?i = 1 and at the grain boundaries both ?i and ?j aresmaller than 1. Therefore, ? = 1 inside the grain (i.e. white region in theimage) and ? < 1 at the grain boundaries (i.e. darker shades of grey in theimage).Secondly, the index of order parameters can be used to distinguish eachgrain in the structure with a random colour associated with the index. Sincethe grain boundaries are diffuse, one needs to determine to which grain agrid point belongs. To do so, the index i is selected from the ind array ifits corresponding ?i is maximum at any point of the domain:Maxind = {ind(:,n?) | eta(:,n?) = max (eta(:,n))n=1,??? ,pn} (5.5)Here, Maxind is an integer matrix representing a discretized grainstructure. Fig. 5.3 shows images of a 2D grain structure created by ? andMaxind arrays.5.3.2 The Statistics of Grain SizeFor small scale structures with few grains, the volume of a grain, V , iscalculated by:V =? lx0? ly0? lz0?i dx dy dz (5.6)This relationship naturally takes the volume of ?i at the diffuse interfaceinto account. Numerically, it is evaluated with the rectangular integrationtechnique, i.e. by the sum of all the values of ?i in the domain. This accuracy69Figure 5.3: Images of a 2D grain structure constructed by (a) ? and(b) Maxind arrays.suffices, since the grain boundary thickness is much smaller than the areaof the grains. As an alternative, for large structures with many grains,the volume of grain i is calculated by counting the number of elements inthe Maxind array. The value obtained from this method is close to thatof Eq. 5.6, since half of the thickness of the grain boundary is indexed tobelong to the grains on each side.During a grain growth simulation, the volume of all the grains is surveyedat 100 timesteps intervals and then stored in a file. These files are later usedto produce statistical data on the grain size distribution and calculating thegrain growth rate, dR/dt.5.3.3 The Topology of GrainsBy specifying an index i in the Maxind array, a grain can be selected as abinary image where the pixels have the value of 1 inside the grain and of0 outside. Fig. 5.4 shows an example for a 3D structure. The algorithmfor calculating the number of faces of the grain i uses an image processingtechnique to ?dilate? the shape of the grain in the binary image by one pixeland intersect it with the Maxind array. It then counts the number of uniqueindices in the image that contain the grain i and a one-pixel margin fromother grains. This process was performed with the image processing toolboxof MATLAB software.70Figure 5.4: The Maxind array of a 3D structure is analyzed tocalculate the number of faces of a grain. (a) A cross sectionof the Maxind array. (b) A binary image of the grain followingdilation. (c) Multiplication of the two images shown in (a)and (b) reveals only the grain and a one-pixel margin of itsneighbours.5.3.4 The Curvature of Grain BoundariesMeasuring grain boundary curvature is challenging both experimentally andin mathematical models. Indirect measurements of curvature have beenproposed using stereological relations [72, 147, 148]. The interface curvaturehas also been calculated for discrete data sets using tessellation [149] andfor structures with diffuse interfaces using differential calculus [150]. A newtechnique is used to measure the curvature of grain boundaries in phasefield simulations. This method creates less noisy results compared withthose employing the method given in [150]. The limitation of this method,however, is that it can only be applied when an analytical relationship forthe diffuse interface profile is known in the curvilinear coordinate system.Consider a parallel iso-?i surface for a boundary between grains i and jthat has a small thickness compared with its radius of curvature. The normaldirection r? is defined as a unit vector normal to an iso-?i surface, pointingfrom inside grain i to the outside. At any point on the grain boundary [103]:?2?i =?2?i?r?2+??i?r?K?i . (5.7)71Here, r? is the distance in the normal direction and K?i is a scalarfield equal to the sum of the principal curvatures (K1 + K2) on an iso-?isurface. For calculating the curvature at each grid point, the value of ?2?iis determined numerically using Eq. 5.2 or Eq. 5.3. The partial derivativeterms in Eq. 5.7 are known analytically from the ?i profile (Eq. 3.8):??i?r?=12?2?m[1? (2?i ? 1)2], (5.8)?2?i?r?2= ?2?m(2?i ? 1)[1? (2?i ? 1)2]. (5.9)To test the numerical accuracy of this method, a shrinkage of a circulargrain in 2D was modelled with m = 2 and ? = 4. The system consists of twoorder parameters that form a circle in the centre of the domain. Fig. 5.5aillustrates a colour map of ?1. The circular grain has the radius of ? = 23.61.The curvature field is calculated for ?1 from Eq. 5.7. A colour map of K?iis shown in Fig. 5.5b for grid points where 0.1 < ?i < 0.9. The curvaturevalues vary in the range of 0.040? 0.044, with the average of K??1 = 0.0418.This makes the radius of curvature 1/K??1 = 23.92, which is close to ? witha 1.3% error. The symmetrical pattern in the curvature field is due to thenegligible degree of error introduced by the numerical anisotropy associatedwith the solution of the phase field equations and evaluating the Laplacian.Figure 5.5: (a) A circular order parameter with a radius of ? = 23.61and (b) the corresponding curvature field.The above procedure can be applied to each order parameter in the grain72structure. At a grain boundary face, two order parameters have profiles withmirror symmetry; therefore, for any point, a pair of values is obtained whereK?i = ?K?j . However, Eq. 5.7, gives inaccurate curvature values close tothe junctions where more than two order parameters have non-zero values.The error is due to the fact that the analytical profile obtained for theequilibrium boundary profile (Eq. 3.8) is not valid at the junctions. Thearea of the grain boundary network close to the junctions can be neglectedif the grain size is much larger than the interface thickness, which is also anecessary condition for the validity of the phase field simulations in general.For statistical sampling, K?i is determined for all order parameters atall grid points falling within the range of 0.2 < ?i < 0.8. Each grid pointhas two values of positive and negative curvature. Fig. 5.6 shows the twocurvature values for a 2D ideal grain growth structure in two separate colourmaps.Figure 5.6: A colour map of (a) positive and (b) negative values ofthe curvature fields for all order parameters.At each grid point on the grain boundary, the magnitude of curvature,K is stored in a list provided that the difference between the curvature ofthe two order parameters is less than 10%, i.e. |K?i | ? |K?j | < 0.1|K?i |.This curvature list is then used to determine the curvature distribution ofthe grain boundary network. The average of the curvature distribution, K?,is defined by:73K? =?K dS?dS, (5.10)where S is the surface area of the entire grain boundary network. Sincethe interface thickness is the same for all grain boundaries, the frequencyof a value in the curvature list is proportional to the length of the grainboundary network having that value. As a result, the arithmetic average ofthe discretized data in the curvature list equals the average curvature (K?)of the grain boundary network.Applications of the aforementioned analytical techniques is presented inthe next chapter, where they are applied to ideal grain growth structures.74Chapter 6Ideal Grain GrowthSimulations6.1 IntroductionStudying ideal grain growth is the first step toward understanding structuralevolution in polycrystalline materials. In this chapter, ideal grain growthin 2D and 3D is modelled with large-scale phase field simulations. Thekinetics of grain growth and the statistical and topological attributes ofthe structure are analyzed and compared with the results of earlier studiesto verify the reliability of the present simulation code. The curvatureof the grain boundary network is measured to reveal the distribution ofthe driving pressure in the structure. A relationship between the averagedriving pressure and the kinetic constant of the parabolic grain growth lawis derived.6.2 Two-Dimensional SystemsIdeal grain growth simulations in 2D were performed with four sets of modelparameters. Tab. 6.1 lists the parameters used in the model along with thegrain boundary energy, mobility and thickness.A system size of 2000 ? 2000 grid points was used for the simulations.75Table 6.1: Phase field model parameters for ideal grain growthsimulations in 2DL m ? ? M lgbSet 1 1 1 2 0.667 3.000 4.000Set 2 1 1 4 0.943 4.243 5.657Set 3 1 2 3 1.155 2.598 3.464Set 4 1 2 4 1.333 3.000 4.000Fig. 6.1 shows two snapshots of the grain structure. The value of the grainboundary curvature is mapped with coulour on the structures.Figure 6.1: The grain boundary network in 2D ideal grain growthsimulations. The values of the curvature field are mapped on thestructure after (a) 20,000 timesteps and (b) 50,000 timesteps.The local curvature of the grain boundary network is not necessarilyconstant across a given grain boundary segment. In addition, a graincan have segments with both convex and concave curvatures. This meansthat, unlike the assumption of the conventional mean field theories, grainboundary curvature is influenced not only by the grain size and the numberof faces, but also by the configuration of the neighbours of a grain. Thisfact influences the growth kinetics of grains.Fig. 6.2 shows the normalized grain growth rate of individual grains(RdR/dt/M/?) for Simulation Set 4, as a function of normalized radius,76u = R/R?. Hillert?s kinetic relationship (Eq. 2.5) is plotted with ? = 0.5and Rcr = R?. The simulation data shows a considerable amount of scatteraround Hillert?s relationship.The colour of the markers in Fig. 6.2 indicates the topological class ofeach grain. The number of faces of a grain is related to its growth rateaccording to the von Neumann-Mullins relationship (Eq. 2.13), re-writtenhere as:RdRdt=M?6(N ? 6) . (6.1)Based on this equation, the values of normalized grain growth on Fig. 6.2should be in multiples of 1/6. Better agreement between the von Neumann-Mullins relationship and the simulation data is observed. However, thescatter shows that, at the level of an individual grain, the grain boundarycurvature does not strictly follow the von Neumann-Mullins assumption andsmall deviations exist depending on the local arrangement of the grains andthe topological aspect of the grain boundary network.Figure 6.2: Normalized growth rate of each grain as a function ofthe normalized equivalent radius for ideal grain growth in 2D.Different colours refer to topological class. The solid line isHillert?s equation (Eq. 2.5). Horizontal lines indicate the graingrowth rates according to the von Neumann-Mullins law.77A holistic approach for the grain growth dynamics is now considered byconsidering the distribution of the entire grain boundary network. Fig. 6.3shows the distribution of the grain boundary curvature at different times.A self-similar distribution is observed when K is normalized by K?; i.e.,when all curves at different times collapse on a master curve. The averagecurvature, K?, is related to the average equivalent grain radius, R?, using aproportionality constant, ?1:K? =?1R?. (6.2)Here, ?1 = 0.20 ? 0.02 is obtained for the self-similar 2D ideal graingrowth structure.0 1 2 3 4 500.10.20.30.40.50.60.70.8K/KLength Fractiontn=20000tn=30000tn=40000Figure 6.3: The self-similar grain boundary curvature distribution for2D ideal grain growth simulations.The driving pressure on a grain boundary is the result of its curvature(Pd = ?K). Similar to the grain boundary curvature, the driving pressurecan be averaged over the entire grain boundary network as:P?d =?Pd dS?dS= ?K?. (6.3)The value of P?d can be viewed as the average driving pressure of grain78growth that is inversely proportional to the grain size. Combining Eqs. 6.2and 6.3 gives:P?d = ??1R?. (6.4)A same relationship as Anderson and Grong phenomenological equation(Eq. 2.37) is obtained here, but the ?1 can be physically quantified withthe curvature distribution. The average grain size is then supposed to beproportional to the average driving pressure, such that:dR?dt= ?2MP?d = M?2?11R?(6.5)where ?2 is a proportionality constant. Integrating Eq. 6.5 gives theparabolic grain growth relationship:R?2 ? R?20 = 2?1?2M?t = kM?t. (6.6)Here, k = 2?1?2 is introduced as the grain growth rate constant(Eq. 2.1). Fig. 6.4 shows the square of the average equivalent grain radius asa function of normalized time ? = M?t. The simulation results, performedwith different grain boundary mobility, energy and thickness, confirm theparabolic grain growth law with k = 0.29? 0.02.The curvature analysis offers further information about the kinetics ofgrain growth. The relationship between the driving pressure and grain sizeis known by the constant ?1. The value of k is measured independentlyfrom grain growth kinetics; thus, according to Eq. 6.6, ?2 = 0.69 ? 0.10 isobtained.79Figure 6.4: Equivalent grain radius as a function of normalized time(? = M?t) for different simulations in 2D.6.3 Three-dimensional SystemsA similar approach is applied to 3D ideal grain growth to obtain the kineticconstants. Simulations were performed on a 300 ? 300 ? 300 domain with?x = 1, ?t = 0.1, L = 1, m = 1 and ? = 2. Fig. 6.5 shows two snapshots ofthe grain boundary structure by mapping the curvature values on the grainboundary position.Figure 6.5: (a) A curvature map of the grain boundary network forideal grain growth at (a) tn = 15000 and (b) tn = 40000.80Fig. 6.6a compares Hillert?s theory with the simulation data where thenormalized growth rate of grains is plotted as a function of their normalizedequivalent radius. Since grain growth is self-similar, the overall position ofthe scatter points does not change with time. The simulation data followsHillert?s relationship (Eq. 2.5) with ? = 1 and Rcr = 9/8R?. In Fig. 6.6bthe normalized volumetric growth rates are plotted against the number offaces of grains and compared with the topological relationships of Mullins(Eq. 2.17) and Hilgenfeldt et al. (Eq. 2.21).Figure 6.6: (a) The normalized growth rate of individual grains as afunction of size for ideal grain growth in 3D. Different coloursrefer to different topological classes. The solid lines representHillert?s equation (Eq. 2.5). (b) The normalized volume changeof grains as a function of the number of faces of the grainsis compared with Mullins (Eq. 2.17) and Hilgenfeldt et al.(Eq. 2.21).To address the overall grain growth kinetics, the distribution of curvatureof the grain boundary network is considered. The self-similar curvaturedistribution in 3D is shown in Fig. 6.7. The value of ?1 = K?R? = 0.27?0.03is obtained.The kinetic constant k = 0.45? 0.04 is determined from the slope of R?2versus ? from three simulation runs (Fig. 6.8). Consequently, for the givenk and ?1, Eq. 6.6 results in ?2 = 0.81? 0.16 for 3D grain growth.81Figure 6.7: Distribution of the grain boundary curvature at differenttimesteps in 3D ideal grain growth.Figure 6.8: Equivalent grain radius as a function of normalized time(? = M?t) in 3D. The variation of the kinetics constant fordifferent simulations is k = 0.45? 0.04.6.4 Discussions and ConclusionsThe results for the kinetics of grain growth are in good agreement with earliersimulations. In 2D systems, k = 0.29 compares well with the followingvalues: 0.27 from phase field [57], 0.27 from vertex model [50], 0.29 fromMonte Carlo simulation [58] and 0.26 from a front-tracking model [54]. In3D systems, the value of k = 0.45 in this work is comparable with the resultof earlier studies, e.g. k = 0.5 [59] and 0.49 [60] for phase field simulations,k = 0.43 for Monte Carlo simulations [61] and k = 0.42 for a modified level-82set model [16]. Lazar et. al [17] obtained k = 0.39 which is smaller than thevalue obtained in this work.The consistency of the results of this study with those of othersimulations is a benchmark for the accuracy of the model. Furthermore,the curvature analysis provides precise information on the average drivingpressure and its relationship with grain growth kinetics that is not evident inother studies. The parameters ?1 and ?2 are kinetic constants of ideal graingrowth and can be applied to any self-similar cellular system with isotropiccurvature-driven dynamics. Table 6.2 summarizes the values for 2D and 3Dsystems.Table 6.2: The kinetics constants of ideal grain growth demonstratedin Eq. 6.6?1 ?2 k = 2?1?22D 0.20? 0.02 0.69? 0.10 0.29? 0.023D 0.27? 0.03 0.81? 0.16 0.45? 0.04The value of ?1 is related to both the geometric factors of the grainshape and the statistical factors of the grain size distribution. To illustratethis fact, assume a large ensemble of grains, each with an equivalent radiusRi and a grain boundary area Si. The average curvature of each grain, Ki,can be measured using Eq. 5.10 by integrating over the area of the grain.Considering a grain as a polyhedra; Si = siR2i and Ki = ki/Ri where siand ki are solely geometrical constants, the average curvature of the entirenetwork can be written as:K? =?SiKi?Si=?sikiRi?siR2i(6.7)One can use simplifying assumptions for the shapes of the grains, e.g, aregular polyhedra or an average N-hedra as proposed by Glicksman [151],to evaluate si and ki. For instance, compared with that of a sphere, si isalways greater than 4pi and ki is smaller than 2, depending on the numberof faces and the shape of a grain. However, as is evident on the right handside of Eq. 6.7, a statistical factor noted by the ratio?Ri/?R2i enters83into the relationship for the average curvature which is a function of thegrain size distribution. The abovementioned analysis allows one to derive?1 analytically for the scaling distribution.The value of ?2 is related to the nature of the averaging process fora given distribution of grains. Note that in Eq. 6.5, the growth rate ofthe average grain size, dR?/dt, is obtained by the arithmetic average of R,whereas the average driving pressure was obtained by an integral over thesurface of the grains. This constant is usually estimated to be 1 in theliterature [2].Hillert?s model (Eq. 2.5) uses ? = 0.5 for 2D and ? = 1 for 3D. Thisagrees with the ideal grain growth simulation results, albeit with scatter,shown in Figs. 6.2 and 6.6b. For 2D, Rcr = R? and for 3D, Rcr = 9/8R? [22].As a result, according to Eq. 2.9, Hillert?s model predicts the rate constantof k = 0.25 for 2D and k = 0.39 for 3D which is lower than the observedvalues in numerical simulations.Phase field simulations produce consistent results with the topologicalrelationship of grain growth. The average grain growth rate of eachtopological class is in good agreement with the von Neumann-Mullins in2D and the Mullins and Hilgenfeldt et al. relationships in 3D. However,considerable scatter is observed for the growth of the grains in each class.Numerical inaccuracies in measuring grain volume and growth rate are muchsmaller compared with the observed scatter. In the simulated structures,grain shapes are observed that are different from those considered in meanfield theories. A grain can have both convex and concave faces and somegrain boundary faces have non-uniform curvature (Fig. 6.1 and Fig. 6.5.Another reason is that at any time step, when the growth rate of a grainis measured, some grains can be in transition from one topological classto another. In addition, in a microstructural evolution model such as phasefield, it takes time for a topological transition to propagate through the grainboundary network and for the junctions to reach the equilibrium angle.84Chapter 7Modelling the Particle-GrainBoundary Interaction7.1 IntroductionThe presence of second phase particles in a structure disrupts motion ofthe grain boundaries. The ideal grain growth kinetic curve slows down andreaches the limiting grain size when all the grain boundaries are pinned.To study the effect of second phase particles on grain growth, first theinteraction of a grain boundary with particles is considered on the smallscale. At this scale, only a fraction of the grain boundary network ismodelled and the grains are considered as semi-infinite domains. This allowsthe resolution of the nano-scale particles with a proper accuracy. In thischapter, the phase field model of Moelans et al. (Eq. 3.16) is used. Theobjective is to investigate the physical aspects of particle-grain boundaryinteractions that have not been reported previously and to discern thepinning pressure that particles exert on a moving grain boundary.857.2 The Pinning Force of a Particle7.2.1 Energetics ApproachThe second phase particles in this work are non-coherent (i.e. no specificcrystallographic orientation) and stable. Therefore, there is no need todistinguish between individual particles and the position of all of theparticles in the domain can be determined by one ? field. Since no evolutionequation is solved for ?, it does not need to have a diffuse interface fromwithin the particle to the matrix. However, the order parameters of thegrains form a diffuse profile around the particles. This diffuse regionconstitutes the particle-matrix interface and has a higher energy densitythan that within the grain (matrix). Fig. 7.1 illustrates energy density ofthe system consisting of a particle and a grain boundary.In the phase field model, the forces between a particle and a grainboundary can not be explicitly obtained since no mechanical interactionsare defined. In order to obtain the pinning force of a particle on a grainboundary from the phase field simulation, an energetics approach inspiredby Ashby is used [66]. When a grain boundary is in contact with a particlesome area of the grain boundary is removed and the energy of the systemdecreases; i.e., ?F = Ap/gb ? ?, where Ap/gb is the area of intersectionbetween the grain boundary and the particle. Fig. 7.1 shows the energydensity of the system for two states: first, when a grain boundary is incontact with a particle and then, when it detaches from the particle. Theenergy difference between the two states is equal to the work that is doneby the pinning force, i.e. Fz = ?F/?x, where x is the position of the grainboundary. An energy well can be introduced for the interaction of a particleand a grain boundary. The instantaneous pinning force is then obtained bydifferentiating F with respect to the grain boundary displacement.The energetics approach to this problem is very suitable for the phasefield method since it is expedient to calculate the total energy of a systemusing the energy functional. To find the pinning force of a particle on a grainboundary, a simple geometry is used where a planar interface moves with a86Figure 7.1: The energy density of a system with a particle and a grainboundary while they interact in (a) and are detached in (b).driving pressure of Pd = 0.050 and interacts with a particle at the centre ofa 3D domain. The system size of lx ? ly ? lz = 80? 80? 80 was used withsymmetric boundary conditions on all of its sides. The model parameters arem = 1, ? = 2 and L = 1, translating into ? = 0.667 and M = 3 (Eqs. 3.12and 3.13). Fig. 7.2 shows snapshots taken from the central cross-section ofthe system. Here, the particle radius, rp = 10.6.Figure 7.2: Cross-section snapshots for a 3D system with a movinggrain boundary and a particle.The mean position of the grain boundary, x, is calculated using the87volume of the grain on the left side:x =1lylz??1(r?) d3r? (7.1)The energy of the system is calculated using Eq. 3.15, where the integralis evaluated using the rectangular integration technique. The total energyof the system and its derivative, F and dF/dx, are shown in Fig. 7.3. Thepinning force is negative at the beginning of the interaction, indicating theattraction of the grain boundary to the particle. This coincides with thedecreasing energy of the system when an area of the grain boundary isremoved by the particle. As the energy of the system increases duringdetachment, the pinning pressure becomes positive and reaches the valueof Fz = 23.20. After detachment, the energy of the system slowly declinesto the same level as before the interaction. The slightly negative force afterthe detachment is due to the assimilation of the anchoring protrusion withthe grain boundary that increases the driving pressure on the grain boundarydue to its curvature.The pinning pressure obtained from analyzing the energy of the systemis close to Zener?s value of FZ = pirp? = 22.34 (Eq.2.28). The deviationincreases as the ratio of the particle size to the grain boundary thicknessdecreases. For the given model parameters, the grain boundary thickness islgb = 4. The deviation in the maximum of the pinning force remains below5% for particle sizes greater than rp = 5.6.To minimize this deviation, the grain boundary thickness should beselected in such a way that it is small compared to the particle diameter. Agreater grain boundary thickness creates artifacts due to the wide range ofits interaction with the particle. In addition, a thick grain boundary formsa larger triple junction between the grain boundary and the particle-matrixinterface that is not considered in sharp interface models.88Figure 7.3: The total energy of the system and its derivative as afunction of the grain boundary position. (a)-(f) denotes theposition of the structure snapshots shown in Fig. 7.2.7.2.2 The Kinetic ApproachThe motion of a grain boundary that is driven by its capillarity is morerelevant to the process of grain growth. In this case, a more complexgeometry and larger simulation box are required to obtain the pinningforce. Here, a geometry proposed by Harun et al. [133] is used. Thisstructure consists of five grains: a semi-infinite grain to create a curvedfront and four columnar hexagonal grains to stabilize the structure. Fig. 7.4shows the interaction of the curved grain boundary with four particles atthe centre of each grain. The system properties are ? = 0.667, M = 3,lx ? ly ? lz = 380 ? 300 ? 264 and rp = 15.6. The periodic boundarycondition is used on the sides of the domain and the symmetric boundarycondition for the faces normal to the x-axis.Columnar grains maintain their shape, creating a suitable curved front89Figure 7.4: The iso-surfaces of ? = 0.65 illustrating the interaction offour particles at the centre of four curved grain boundaries attimes (a)100, (b) 600 and (c) 750.with constant driving pressure. To measure the grain boundary velocity,Eq. 7.1 was used where ?1 was the order parameter of the semi-infinite grainto obtain the mean position of the grain boundary. A smoothing spline curvewas fitted onto the x versus time data using MATLAB software, and thederivative of the splines were calculated at each time step. Fig. 7.5 shows theinstantaneous grain boundary velocity as a function of the grain boundaryposition. The pinning force at each position can be deduced by comparingthe grain boundary velocity with or without its interaction with the particle.The overall pinning pressure on the grain boundary is calculated by:Pz =v0 ? vM, (7.2)where v0 is the velocity of the grain boundary without interaction with theparticle and v is the instantaneous grain boundary velocity. This pressureis created by four particles exerting the same pinning force; i.e:Fz =14Pzlylz. (7.3)The pinning force is plotted in Fig. 7.5. The maximum pinning force,Fz = 30.7, is close to the Zener force of FZ = pi?rp = 32.7.The curve in Fig. 7.5 for the velocity of the grain boundary resembles theresults by Couturier et al. [8] (Fig. 2.13), except that, a narrower attractionforce region is observed in the finite element model. A grain boundary90Figure 7.5: The velocity of a grain boundary and the correspondingpinning force for interaction of a particle with a curved grainboundary, rp = 15.6.with exaggerated thickness in the phase field models can interact with theparticle sooner, thus extending and amplifying the negative pinning force.However, it is expected that the phase field model converges to the behaviourof sharp interface models when the interface thickness becomes less thanthat of the particle radius. Phase field simulations were performed with thegrain boundary thickness of lgb = 4 and with particle sizes of rp = 5.6 andrp = 15.6. Fig. 7.6 compares the normalized pinning force, Pz/(pi?rp), forthe two phase field simulations and the analytical model of Hellman andHillert [67] as well as finite element [8] and front tracking [133] simulations.Hellman and Hillert only considered the anchoring of the grain boundary bythe particle and also in the front tracking simulations only the attractionpart of the curve was reported.Front tracking and phase field simulations with ?/rp = 10 producesimilar results. The Hellman and Hillert analytical solution, however,underestimates the extent of grain boundary anchoring. The negative forceregion was only reported for the finite element model and it has a differentgrain boundary curvature to particle radius ratio. Nonetheless, it is expectedthat the effect of ?/rp ratio to be much less pronounced for the attractionof the grain boundary to the particle than it is for the anchoring process.91Figure 7.6: The normalized pinning force, Pz/(pi?rp) as a functionof the normalized position of the grain boundary for severalsimulation methods.The pinning force varies strongly as the grain boundary moves throughthe particle. Notably, the acceleration of grain boundary motion by particlesat the beginning of the interaction is considerable. Unlike Zener?s theory, theinteraction of the grain boundary and particles is not limited to a range from?rp to rp but can extend from ?rp up to four times the radius of the particleduring anchoring. As a result, the number of particles interacting with agrain boundary substantially increases compared with Zener?s assumption.It is important in phase field simulations to use a large particle sizecompared with the grain boundary thickness. The negative portion of thepinning force curve is considerably widened when lgb/rp < 2 and can resultin an enlarged negative pinning pressure when a distribution of particlesinteracts with a grain boundary. This fact has been ignored in previousphase field simulations of grain growth with particle pinning [12, 131] toreduce computational costs.927.3 The Pinning Pressure of a Distribution ofParticles7.3.1 The Planar Grain BoundaryThe interaction of a moving grain boundary with a distribution of particlesis a complex process. The anchoring process changes the shape of the grainboundary locally and, as a consequence, the changes in the local curvaturealter the pinning force on the particles. The phase field method is used toinvestigate this process and obtain the pinning pressure. In this approach, adriving pressure of Pd is applied on a planar grain boundary. The assumptionof a planar grain boundary interacting with particles is not far fetched since,in a grain structure with very fine particles, the scale of the grain boundary-particle interaction is much smaller than the radius of curvature of the grainboundaries. In addition, a driving pressure exists on planar grain boundariesin systems with recrystallization or phase transformation.The simulations are similar to those of single particle-grain boundaryinteractions (Sec. 7.2.1) but with a much larger domain size. Here, a domainsize of lx ? ly ? lz = 300 ? 300 ? 300 with symmetric boundary conditionis used. Different simulations were performed with a driving pressure ofPd = 0.05. The simulation parameters are ? = 0.667, M = 3, ?x = 1 and?t = 0.05. Particle distribution is specified by the ? field (Eq. 3.16) in whichspherical regions with a radius of 10 grid points are randomly distributed.The volume fraction of the particles is 0.047. The planar grain boundary isinitially situated at x = 30 and moves to x = 270. Fig. 7.7 shows a snapshotof the system after 5000 timesteps.At any moment, the grain boundary encounters different local arrange-ment of particles and is in contact with several particles. As the grainboundary anchors on each particle, the shape of the grain boundary changes.The ?average? position of the grain boundary, x, is inferred from the volumeoccupied by ?1, taking into account the particles volume fraction. Thevolumetric integral of ?1 does not include the volume of the particles, since?1 is zero within the particles, but the volume swept by the grain boundary93Figure 7.7: A planar grain boundary moves through a distribution ofparticles with fv = 0.047 and rp = 10.6, with a driving pressureof Pd = 0.050.should include the entire domain, i.e:x =1lylz1(1? fv)??1(r) d3r (7.4)The velocity of the grain boundary at each position can be obtainedby fitting a smoothing spline curve on x versus time and calculating thederivative of the spline curve. The pinning force exerted by each particlecould be at a level ranging from negative to positive depending on theposition of the grain boundary relative to the particle. The pinning pressurecan be perceived as the collective force of all of the particles on a grainboundary that opposes the driving pressure; i.e:Pz = Pd ?vM. (7.5)Fig. 7.8 shows the instantaneous grain boundary velocity and thenormalized pinning pressure, Pz/(?fv/rp).The average pinning pressure is P?z = 1.66?fv/rp. The fluctuations in the94Figure 7.8: The instantaneous grain boundary velocity and pinningpressure for a planar grain boundary moving under a the drivingpressure of Pd = 0.05 and interacting with a particle distribution(fv = 0.047 and rp = 10.6). The horizontal solid line is thevelocity of the grain boundary without particles, and the dashedline is the average pinning pressure.pinning pressure can be explained by the variation in the local configurationof particles in contact with the grain boundary. The amplitude of thefluctuations is a function of the number density of the particles on the grainboundary. This system has N = 322 particles in the domain. Consideringthe particle-grain boundary interaction range to be 3r, this means that onaverage, the grain boundary interacts with Nv ? 3rlylz = 32 particles. Ahigher number density of particles or a larger domain size reduces the levelof fluctuations. Therefore, in a grain structure, the variation of the pinningpressure is a function of the ratio of the grain size to the particle inter-spacing, lp ? rp/ 3?fv. For fine particles where lp << R, the uniform pinningpressure is a good approximation, whereas for a coarse distribution ofparticles where lp is in the same order of magnitude as R, higher fluctuationscan be observed. In these situations, some grain boundary segments maybe completely pinned while others are free to move.957.3.2 Curved Grain BoundaryIn this section, the interaction of a curved grain boundary with a distributionof particle is examined. To construct a grain boundary with a stablecurvature, a geometry of columnar grains with the shape of hexagonal prismswas used. The system setup resembles the structure used in Sec. 7.2.2, buton a larger scale to allow for the interaction of a statistically significantnumber of particles with the grain boundary. The simulation parametersare ? = 0.667, M = 3 and lx ? ly ? lz = 380 ? 300 ? 264 with fv = 0.048and rp = 10.6. Fig. 7.9 shows a snapshot of the structure following 5000timesteps.Figure 7.9: The motion of a curved grain boundary through adispersion of particles with fv = 0.048 and rp = 10.6.The flat grain boundary segments along the columnar grains are slightlydistorted by the presence of particles. However, they remain stationaryduring the simulation and have a negligible effect on the curvature of themoving front.The position of the grain boundary is calculated using Eq. 7.4. Asmoothing spline curve is fitted on the x versus time data and the grain96boundary velocity is calculated from the derivative of the spline curve. Thepinning pressure is obtained from Eq. 7.2, in which v0 was calculated froma simulation employing the same geometry, but in the absence of particles.Fig. 7.10 shows the grain boundary velocity and the normalized pinningpressure, z = Pz/(?fv/rp) (Eq. 2.30) as a function of the grain boundaryposition.Figure 7.10: The velocity of the grain boundary and the normalizedpinning pressure as a function of the grain boundary positionfor the system shown in Fig 7.9.The average pinning pressure for the motion of a curved grain boundaryis:P?z = 2.02?fv/rp. (7.6)The value of z = 2.02 here is higher than the value for the planargrain boundary, which is close to z = 3/2, as predicted by Zener?s theory.The agreement here is rather fortuitous. There are two main simplifyingassumptions in the theory that are not valid in the computer simulations.First, the particles in contact with the grain boundary do not exert themaximum pinning force, but depending on the position with respect to thegrain boundary, the pinning force exhibits a positive or negative value. Thissubstantially lowers the total force on the grain boundary and thus decreases97the pinning pressure compared with Zener?s prediction. Secondly, in Zener?stheory, a grain boundary remains completely flat, while in computer modelsthe grain boundary is flexible and can interact with more particles per cross-section area of the domain and thus it increases the pinning pressure. Thesefactors influence z in two opposite ways.The higher pinning pressures for the curved grain boundary demonstratethat the pinning pressure depends on the source of the driving pressure andthe geometry of the grain boundary. If a grain boundary is driven by itscurvature, the driving pressure in the vicinity of the particle is altered locallydue to the negative curvature caused by the anchoring shape. This effectivelyreduces the driving pressure and increases the pinning pressure.The pinning pressure from the phase field simulations obtained here issmaller than the value of Pz = 3.3?fv/rp from the finite element models [8].The main difference in the simulation of Couturier et al. is the fact that thepinning pressure is reported when the grain boundary motion has completelyhalted. It has been postulated by Hillert [6] that a larger driving pressureis required to unpin a grain boundary from a distribution of particles ascompared with the driving pressure required to sustain the grain boundarymotion through it. A moving grain boundary interacts with particles bothin the positive and negative force positions. However, a stagnated grainboundary is only in contact with particles on the anchoring or the zeroforce position. This configuration is formed during the moments prior to itscomplete halt, when the grain boundary stops next to a few particles, ata local region, and then surpasses other particles in the front due to theirattraction (negative pinning force), but is not able to move forward further.Without the attraction force, higher driving pressure is needed to detachthe grain boundary and resume the movement.The 30% increase in the pinning pressure for the anchored grainboundary compared with that of the mobile grain boundary is consistentwith the relative area under the pinning force curve for the negative andpositive parts shown in Fig. 7.6. In addition, the grain boundary had alarger ?/rp in the finite element model compared with that of the phasefield, and thus had a larger interaction range.98Considering all factors influencing the particle-grain boundary interac-tion, it is suggested that the pinning pressure on a grain boundary drivenby its curvature is in the range of 2 to 3 times the value of ?fv/rp,depending on the ratio of ?/r, and can increase to 3.3 times when thegrain boundary motion is halted. Note that these values are obtained withthe assumption of stable, non-coherent and spherical particles with randomspatial distributions that may not be applicable to particle pinning in allmaterials.The fluctuation of the pinning pressure observed in Fig. 7.10 is due to thevariation of the instantaneous velocity of the grain boundary. The amplitudeof this fluctuation depends on the size of the simulation domain since itis the function of how many particles at any moment interact with thegrain boundary. In a suitably large domain, the average number of particlesinteracting with the grain boundary remains constant and thus less variationin velocity is observed. The pinning pressure obtained in Eq. 7.6 is presentedhere as a general trend and as a demonstration of the simulations in nano-scale. More simulations with variety of particle sizes and volume fractionsare required to obtain an exact value and its standard deviation.7.4 The Interaction of Grain Boundary Junctionswith ParticlesOne may consider junctions as topological regions of the grain boundarynetwork where more than two grains meet. Triple lines are 1D defectscreated by the intersection of three grain boundary faces, and quadruplepoints are formed at the intersection of four grains. It is expected thatthe pinning force of a particle on a junction is stronger than on a grainboundary face, since the positioning of a particle at a junction removesa greater amount of grain boundary surface area. Based on the sharpinterface assumptions, the maximum intersection area for triple lines isAp/gb = 3/2pir2p and for quadruple points is Ap/gb = 2pir2p [152].The system with four hexagonal prismatic grains is used here, with sixparticles positioned on a plane at the centre of the domain. In Fig. 7.11a,99when the curved grain boundary front moves forward, the particles areplaced at points where they interact with the triple lines and then detachto the matrix of the semi-infinite grain. Similarly, in Fig. 7.11b, theparticles are positioned where they can interact with the quadruple pointsbefore detachment. Simulation parameters are ? = 0.667, M = 3 andlx ? ly ? lz = 150? 100? 80, with rp = 5.67.Figure 7.11: The interaction of four columnar grains with six particlespositioned at (a) triple lines and (b) quadruple points. The iso-?1 = 0.5 surface is plotted in (a) and the iso-? = 0.6 surface isplotted in (b).The velocity of the grain boundary and the pinning force for the twocases are plotted in Fig. 7.12. A smaller attraction force is seen since theparticles are already positioned on the flat segments of the grain boundarynetwork. On the other hand, a large pinning force is associated with thedetachment of the particles from the junctions to the matrix.The pinning force for the triple junction is Fz = 1.41pi?rp and for thequadruple point is Fz = 1.69pi?rp. These forces can be compared withFz = 0.93pi?rp when a particle interacts with the grain boundary face. Thehigher force exerted by a particle positioned at a junction can change theoverall pinning pressure on a grain boundary. This effect becomes importantwhen the number of particles on a grain boundary face is comparable to thenumber of particles at the junctions. For example, in a nano-crystalline100Figure 7.12: The velocity of the grain boundary and the pinning forceexerted from each individual particle on the grain boundaryfor (a) particle interaction with triple lines and (b) particleinteraction with quadruple points.material with a grain size of 100 nm, rp =10 nm and fv = 0.5%, the numberdensity of the quadruple points is almost equal to the particle numberdensity. However, this number significantly decreases for larger grain sizesbeyond 1 ?m where there is only 1 quadruple point per 1000 particles.7.5 ConclusionsThe physical phenomenon of particle-grain boundary interaction has re-sounding similarities with the concept of friction in mechanics. Each particleacts as a small energy well that the driving pressure on the grain boundaryneeds to overcome. The particles have a spatial distribution; thus, thevelocity of the grain boundary is subject to fluctuations both spatially andtemporally on the small length-scale. What complicates this process furtheris that, since the grain boundary is flexible, it can deform and extend itsinteraction range depending on the local arrangement of the particles. Thecontribution of all forces during the attraction of the grain boundary to theparticle or during the anchoring, amounts to a pinning pressure on the grainboundary.It is shown that the pinning pressure on a moving grain boundary, i.e.101the dynamic pinning pressure, is lower than the pinning pressure on ananchored grain boundary, i.e. the static pinning pressure. This fact and thestochastic nature of the particle-grain boundary interaction has importantphysical consequences, when such fluctuations can trigger the motion ofa specific grain boundary segment while other segments are pinned. Thisprocess may be responsible for the nucleation of abnormally large grainsduring heating processes, where the partial dissolution of particles occurs,and should be investigated further.The pinning pressure also depends on how a grain boundary is driven.Smaller pinning pressures are observed when a planar grain boundary ismoving by the application of a driving pressure (e.g. in recrystallization)compared with curvature driven grain boundary motion. The differenceof 25% observed here is, however, difficult to observe in experiments sincethere are usually larger uncertainties in the measurement of stored energyand grain boundary energy.It has been shown that the interaction of a particle with a grain boundaryjunction can create a pinning pressure 1.5 times larger than that for a grainboundary face. This additional pressure is considerable when grain sizeand the distance between particles are comparable. Unfortunately, this hasbeen the situation for most numerical simulations of particle pinning whereparticles were resolved in the same domain as the grains [11, 12]. It rendersresults of such simulations inapplicable to studying particle pinning in manytechnologically important materials with nano-scale particles and meso-scalegrains. The true physical interpretation of the simulated structures shownin Fig. 2.11 is the grain growth of a nano-crystalline material.This amplifies the need for a model that can include the effects of pinningpressure in the structure without resolving the particles. The inspiration forsuch a model can be drawn from the friction phenomenon where a uniformpinning pressure acts on grain boundaries.102Chapter 8The Friction Pressure Modelfor Grain Boundary Drag8.1 IntroductionThe collective force of the fine particles or solute atoms on a grain boundaryin meso-scale can be viewed as an application of a uniform pressure onthe grain boundary, opposing the motion of the interface. On this scale,the underlying mechanism of drag is immaterial, but the effect of thedrag pressure on the evolution of the microstructure is of significance. Inthis chapter, a phase field model is developed that applies a constant orvelocity dependent drag pressure on the grain boundaries. First, a methodis introduced to exert driving pressure on the interface in the the Chen-Yang phase field model. The applied pressure is then made dependent onthe direction and velocity of the grain boundary. The ?friction pressure?formulation is tested numerically to reproduce the expected behaviour of asharp interface under pinning and solute drag pressure.1038.2 Phase Field Formulation8.2.1 Applying Pressure on the InterfaceThe pressure on an interface can be seen as an energy difference betweenthe two sides of the interface. The energy difference can be modelled bychanging the internal energy of each order parameter with respect to thestandard state. Thus, the total energy functional of the system is writtenwith a modified local energy density, floc:F =? [floc + ?p?i=1|~??i|2]d3r, (8.1)where, floc consists of two parts: (i) the variation of the local energy densityacross the interface (f0) normalized by the standard state as given by Eq. 3.6and (ii) a contribution due to the internal bulk energy of the grains (fd), i.e:floc = f0 + fd. (8.2)The standard local energy density, f0, was discussed in Sec. 3.4.1. It ischosen to have a minimum of 0 inside each grain and a maximum at theinterface resembling a double-well shape. Here, the focus is on defining abulk energy, Gi, for a grain i, in comparison to the standard state. In thiscase, the value of the local energy density, floc, inside each grain is Gi and therole of the fd function is to interpolate the bulk energy of each grain acrossthe grain boundary. For two order parameters, the contribution toward thelocal energy density due to the bulk energy is written as:fd = 32?i=12?j 6=i(?2ij2??3ij3)?Gij +122?i=1Gi (8.3)where, ?ij is 0 inside grain i and 1 inside grain j, i.e:?ij =?j ? ?i + 12(8.4)Fig. 8.1 shows the local energy density, floc of a system with G1 = 0.3104and G2 = 0. The value of floc at ?i = 1 and ?j = 0 is equal to the Gi values.Across the interface, the local energy density follows a curve of the doublewell shape.Figure 8.1: The local energy density, floc, as a function of the orderparameters in a system with two grains. The solid line is theenergy density associated with the equilibrium value of orderparameters across the interface for G1 = 0.3 and G2 = 0In a case where the bulk energy of all grains is equal (Gi = Gj), fd =const. and can be omitted. The driving pressure of the grain boundarymovement is only due to its curvature. However, in a case of different bulkenergies, Gi 6= Gj , there is an additional driving pressure of ?Gij = Gj?Gion the grain boundaries.According to Eq. 3.10, the evolution of the system follows the Allen-Cahn relaxation dynamics, resulting in a set of partial differential equations(PDEs):??i?t= ?L(?f0??i+?fd??i? ??2?i). (8.5)The condition ?j = 1 ? ?i is met at the grain boundaries [122]. As aresult, the derivative of fd can be simplified to:?fd??i= 3?i?j?Gij . (8.6)105It is shown in Appendix A that Eq. 8.6 satisfies the sharp interfacelimit assumptions. Each pair of order parameters between ?i and ?j hasa contribution to the ?fd/??i term. Therefore, a set of evolution PDEs isobtained by generalizing Eq. 8.6 over all the order parameter pairs:??i?t= ?L??m???3i ? ?i + 3?ip?j 6=i?2j??? 3?ip?j 6=i?j?Gij ? ??2?i?? .(8.7)This is a more general form of the standard phase field equation(Eq. 3.11), where, in addition to the interface curvature, the driving pressureof ?Gij exists between each pair.8.2.2 Applying Drag PressureIn the previous section, it was described how to apply a driving pressure onan interface by shifting the local energy density. In this section, the samemethod is used to apply a drag pressure on the interface in meso-scale whilethe particles or solute atmosphere are not explicitly resolved. For the sakeof simplicity, we assume that the driving pressure Pd on an interface is onlydue to its curvature. A drag (or friction) pressure, Pf , can be consideredas a negative driving pressure or a pressure that opposes movement of theinterface. If the driving pressure is smaller than the drag pressure, theinterface stops moving. To create the effect of drag, the local energy densitycan be shifted depending on the direction of the interface velocity as it isshown schematically in Fig. 8.2. In general, the drag pressure can be afunction of the normal interface velocity vij of grain i moving into grain j.In such a formalism, the energy shift between grains i and j can be writtenas:?Gij(r) ={?Pf (vij) Pd > Pf?Pd Pd < Pf(8.8)Applying this approach to the phase field method, we define the normalvelocity in the curvilinear coordinate system as a scalar velocity field for106Figure 8.2: A schematic representation of the friction pressure on aninterface. The local energy density of two neighbouring grainsin (a) and (c) is shifted in order to create an opposing pressureto the interface movement. (b) When the driving pressure issmaller than the pinning pressure, the interface stops.each order parameter ?i using the chain rule, i.e.,v?i(r) =?ni?t=??i?t?ni??i. (8.9)Here ni is the distance along the curvilinear coordinate ni normal to theinterface pointing from within grain i, into grain j. The order parameterprofile across the interface is known analytically, i.e. [122]:??i?ni=12?m2?[1? tanh2(?m2?ni)]=12?m2?[1? (2?i ? 1)2] . (8.10)At each point on the ni axis through the boundary between grains iand j, v?i has the same value and is equal to the local grain boundaryvelocity v?i = |vij | and v?i = ?v?j . According to Eq. 8.8, the value of thelocal energy shift ?Gij in Eq. 8.7 can be obtained with the condition thatPf (?v?i) = ?Pf (v?i), since ?Gij = ??Gji. The curvature driven pressurein the phase field is distributed across the interface. Similar to Eq. 8.6,f ? = 3?i?jPd where f ? is:107f ?i =?f0??i? ??2?i. (8.11)Consequently, by substituting appropriate values for the energy shift inEq. 8.7, the phase field equation can be obtained:??i?t= ?L[?f0??i? ??2?i ? P?f]. (8.12)Here P ?f is given by:P ?f ={3?i(1? ?i)Pf (v?i) |f?i |>|3?i(1? ?i)Pf |f ?i |f?i |?|3?i(1? ?i)Pf |(8.13)such that the boundary becomes immobile when the drag pressure is largerthan the driving pressure due to its curvature.Eq. 8.12 is applicable to any general form of grain boundary drag. Inthe special case of pinning by a stationary particle distribution, the pinningpressure is a constant opposing the motion of the interface. Consequently,considering a constant and positive value for the pinning pressure, Pz, thefriction pressure in Eq. 8.13 (for Pd > Pf ) is written as:Pf = ?sgn(v?i)Pz. (8.14)According to Eq. 8.5, f ?i and d?i/dt have opposite signs. Therefore, thedrag pressure term (Eq. 8.13) in the phase field equation simplifies to:P ?f ={3?i(1? ?i)sgn(f ?i)Pz |f?i |>|3?i(1? ?i)Pz|f ?i |f?i |?|3?i(1? ?i)Pz|(8.15)In another special case, i.e. sole solute drag, the drag pressure, Pf = ?Psis never larger than the driving pressure. In other words, the grain boundarynever stops unless the driving pressure goes to zero. Since the case of Fig.8.2b never occurs, the resulting drag pressure term is given by:P ?f = ?3?i(1? ?i)Ps(v?i). (8.16)108An accurate solution of Eq. 8.12 with the drag term given by Eq. 8.16requires an iterative process, since velocity is also a function of ??i/?t.A typical algorithm should, at a time step tn, first calculate v?i from theprevious value of ??i/?t at a grid point, then solve the PDE for the newvalue of ??i/?t and repeat the process until convergence of v?i is obtained.Using a sufficient number of grid points through the interface (lgb/?x >4), convergence is attained within four to five iterations. Numerically, thevelocity field can be obtained only at the interface where d?i/dni is not zero.8.3 Numerical Testing and Benchmarking8.3.1 Constant Pinning PressureA simple geometry of a circular grain shrinking inside a larger grain isconsidered to test the accuracy of the proposed phase field model. Thephase field equation (Eq. 8.12) is solved with the drag pressure term givenby Eqs. 8.15 and 8.16, respectively. The grid spacing of ?x = 1, time step of?t = 0.5, m = 2 and ? = 4 are used with the forward Euler finite differencemethod. Only two order parameters are considered. The initial conditionis created by setting ?1 to the value of 1 inside a circle with a prescribedradius and the value of 0 outside. ?2 is initialized to 1 ? ?1 in the entiredomain to create the outer grain. The phase field equation is solved withthe absence of drag pressure for 300 time steps in order to reach the steadystate profile shape for the order parameter across the grain boundary.The resulting circular grain is set to shrink with the drag pressure. Inthis case, the driving pressure on the interface (Pd) is only due to the grainboundary curvature. For a circular grain, Pd = ?/?, where ? is the radiusof the circle. Since the grain always remains a circle, ? is calculated fromthe area of the domain under ?1. Interface energy and mobility are givenby Eqs. 3.12 and 3.13. The driving pressure increases monotonically asthe circle shrinks. Thus, the grain boundary velocity can be obtained forlow to high driving pressures in one simulation. The velocity of the grainboundary is calculated from the radius of the circle as a function of time.109Figure 8.3: The grain boundary velocity of a shrinking circular grainas a function of driving pressure in phase field simulations withpinning pressure.The velocity is then compared to the corresponding analytical expressions forthe shrinkage of a circular grain. Based on the sharp interface assumption,the following relationship is expected to be observed in the model.v ={M(Pd ? Pf ) Pd > Pf0 Pd < Pf(8.17)For the case of ?Zener pinning?, a constant pinning pressure is consid-ered. Eq. 8.12 is solved for the constant drag pressure Pz in the range of 0 to0.0354. In Fig. 8.3, the normalized velocity of the interface is plotted againstthe normalized driving pressure. Each line represents a simulation with adifferent Pz. The solid line with slope 1 represents a shrinking circle withoutdrag pressure. As the pinning pressure increases, the curves shift towardshigher driving pressures but the slope remains unchanged. The velocity ofan interface with a driving pressure below Pz is zero as implemented byEq. 8.13. These phase field simulations reproduce the results according toEq. 8.17.1108.3.2 Velocity-Dependent Solute Drag PressureFor a velocity dependent drag pressure caused by solute drag, a generalrelationship such as Cahn?s (Eq. 2.43) can be used, i.e. Ps = av/(1 + bv2).The phase field equation, Eq. 8.12, is solved for the shrinking circulargrain with a = 1 and b = 200 to 50, 000. The normalized velocity ofthe grain boundary is plotted in Fig. 8.4 as a function of the normalizeddriving pressure and compared with the analytical solutions. The analyticalsolution for a grain boundary velocity as a function of driving pressure canbe obtained by inserting Cahn?s equation, Eq. 2.43 into Eq. 8.17 (wherePf = Ps) and solving for v with respect to Pd.Figure 8.4: The normalized grain boundary velocity of a shrinkingcircular grain as a function of normalized driving pressure inthe presence of solute drag pressure. Simulated curves (dotmarkers) are compared with the analytical solution (solid lines)for different values of b.The parameter a determines the slope of the low velocity regime and bdetermines the driving pressure where the transition from the low to the highvelocity regime occurs. For all calculations, there is good agreement betweenthe solution of the phase field equation and the corresponding analyticalsolution. These two examples demonstrate that the phase field equationssatisfy the presumed sharp interface limit assumptions.111At grain junctions, more than two order parameters have a non-zerovalue. The application of Eq. 8.12 at the junctions distributes the dragpressure to the existing order parameters however, since??i 6= 1 [122], thecombined drag force is slightly higher than the expected value. The effect ofthis error on the overall grain boundary movement is small if the diametersof the grains are much larger in comparison to the interface thickness andthe area under the junctions. To test this effect, the shrinkage of a four-sided grain was simulated with symmetric boundary condition (Fig. 8.5)using Eq. 8.15. The velocity of the shrinking grain?s interface was comparedwith a corresponding sharp interface assumption in Eq. 8.17. There is lessthan 5% error in the observed value of Pf and less than 1% for M for grainswith diameters as small as 5 times the interface thickness.Figure 8.5: A four sided grain shrinks in a symmetric boundarycondition from an equivalent radius of (a) 37 grid points to(b) 13 grid points. The interface thickness is approximately 5grid points.8.4 ConclusionsA phase field formulation was presented here that applies a friction pressureon grain boundaries and models the presence of drag pressure. The modelis specifically evaluated for particle pinning with constant drag pressure andsolute drag with a velocity dependent drag pressure. The intrinsic nature ofthe interface drag phenomenon is not taken into account, but the resultingdrag pressure is an input to the model. One can choose a constant dragpressure, any velocity dependent drag pressure or even a combination of112both to simulate the evolution of the grain structure.This approach can be readily applied to grain growth simulationsand when an additional driving pressure ?Gij 6= 0 is present, e.g. forrecrystallization and phase transformations. To apply the proposed methodto a specific material, it will be necessary to find the relevant friction pressure(Pz, a and b, etc.) from independent experimental or simulation studies.The drag pressure in this model is applied uniformly in the domain. Inprinciple, this is true for cases where the particles are very fine compared tothe grains and each grain boundary segment interacts with a large numberof particles. In addition, placing particles at grain boundary junctions mayinsert different drag pressures than what is calculated using this method.However, if the size of grains are much larger than the particles or the solutesegregation spike and the resulting number density of the triple junctionsare small, the error is negligible on the overall grain growth kinetics.113Chapter 9Grain Growth With DragPressure9.1 IntroductionThe friction pressure model allows the evolution of a grain structure to bemodelled with a uniform drag pressure in meso-scale. Since second phaseparticles and/or solute segregation spikes are not resolved, it is feasible toperform a large scale simulation of grain growth with the presence of drag.In Sec. 9.2, the result of grain growth in 2D and 3D systems with constantpinning pressure is presented. The curvature of grain boundaries is analyzedto obtain the distribution of driving pressure. The limiting grain size andthe kinetics of grain growth are analyzed and compared with the availableexperimental data. In Sec. 9.3, grain growth in 2D systems with solute dragpressure according to Cahn?s model is considered and the kinetics of graingrowth is discussed.1149.2 Grain Growth with Pinning Pressure9.2.1 Grain Structure and Limiting Grain SizeGrain growth in 2D is simulated in systems with 2000 ? 2000 or 2400 ?2400 grid points, depending on the grain size, to ensure the presence of astatistically significant number of grains. A structure with self-similar grainsize distribution is used as an initial condition. The pinning pressure isapplied in the range of Pz = 0.01 to 0.05. Model parameters are m = 2,? = 3 and L = 1 with ?x = 1 and ?t = 0.1. In all simulations, the appliedpinning pressure is selected to be lower than the initial driving pressure,allowing grain growth to occur.Fig. 9.1 shows 2D grain structures for different times in a simulation witha pinning pressure of Pz = 0.02. At the beginning of the simulation, graingrowth occurs in all regions of the domain. As the simulation progresses,some grain boundary segments with lower curvature and stable topologicalneighbourhoods (grains with 6 sides) freeze. An example of such regionsis shown with the box in Figs. 9.1b and 9.1c. The fraction of the frozenregions increases with time and eventually results in a total halt in the grainboundary motion. The final structure consists of grains of different sizesand topological classes. The frozen grain boundary segments can be curvedor flat as shown in Fig. 9.1c with arrows I and II, respectively.Figure 9.1: Snapshots of a 2D grain growth simulation with Pz =0.020, ? = 1.1547 and M = 2.598 at times of (a) 200 (initialstructure), (b) 2000 and (c) 20000 (final frozen structure).115The average grain radius versus time is plotted in Fig. 9.2a for differentpinning pressures. All systems start from the same initial condition. Thegrain growth curve soon deviates from the ideal parabolic grain growthlaw and approaches a limit. The limiting grain size, R?lim, is inverselyproportional to the pinning pressure Pz, as shown in Fig. 9.2b:R?lim = 0.73?Pz? ? ? ? ? ? (for 2D). (9.1)Figure 9.2: (a) Equivalent grain radius as a function of time for 2Dgrain growth simulations with different pinning pressures. (b)Limiting grain size as a function of the inverse pinning pressure.Grain growth simulations in 3D are performed in systems with dimen-sions of 300 ? 300 ? 300. Model parameters are L = 1, m = 1 and ? = 2with ?x = 1, ?t = 0.05. Two initial structures with different grain sizes areused with pinning pressures of Pz = 0.010?0.030. The initial structures aretaken from an ideal grain growth simulation in two time steps. This allowsone to choose a desired initial grain size with self-similar distribution. The3D simulations with pinning are the most expensive in this study due totheir size and the long periods of time required to reach the limiting grainsize. Fig. 9.3 shows two snapshots of the structure with Pz = 0.010.The average grain radius R? as a function of time is shown in Fig. 9.4a for3D simulations with different pinning pressures. As shown in Fig. 9.4b, the116Figure 9.3: The grain structure of a 3D simulation with the pinningpressure of Pz = 0.010 at (a) timesteps tn = 2000 and (b) atthe limiting grain size when tn = 118000.Figure 9.4: (a) The equivalent grain radius as a function of time for3D grain growth simulations with different pinning pressures.(b) Limiting grain size as a function of the inverse pinningpressure.limiting grain size, R?lim is inversely proportional to the pinning pressure:R?lim = 0.60?Pz? ? ? ? ? ? (for 3D). (9.2)The experimental results for the limiting grain size have been compiled117by Manohar [10] for a variety of systems (Fig. 2.8). The relationshipR?lim = 0.17?rp/fv was fitted on the data for systems with fine particlesand small volume fractions. The same data is replotted in Fig. 9.5 and iscompared with the results of limiting grain size from different simulationsand analytical models.The limiting grain size depends on the pinning pressure as shown withthe friction pressure model in meso scale (M-PFM). In order to comparethis result with the experimental data, a relationship between the particledistribution and the pinning pressure is required. As discussed in Sec. 7.3,this relationship can be obtained from nano-scale simulations where particlesare resolved. From the nano-scale phase field calculations (N-PFM), Pz =2.02?fv/rp was obtained for a curved interface interacting with an ensembleof particles (Eq. 7.6). Combining this Pz with the results from the meso-scale friction pressure model (Eq. 9.2) establishes the prediction R?lim =0.30?rp/fp. Similarly, the pinning pressure of Pz = 3.3?fp/rp (Eq. 2.40)was obtained from finite element simulations by Couturier et al. [8]. Usingthis value in Eq. 9.2 results in R?lim = 0.18?rp/fp.10-4 10-3 10-2 10-1101102103fvRlim/rpManohar [10]Zener [5]Miodownik [11]Suwa [12]M-PFM+N-PFMM-PFM+Couturier? = 1.33, ? = 1? = 0.73, ? = 1.02? = 1.42, ? = 0.87? = 0.30, ? = 1? = 0.18, ? = 1R?lim = ?rpf?vFigure 9.5: The experimental data on limiting grain size as a functionof particle size and volume fraction is compared with thesimulations. M-PFM+N-PFM refers to the combination ofEqs. 9.2 and 7.6. M-PFM+Couturier refers to the combinationof Eqs. 9.2 and 2.40.118The best agreement with the experiments is observed for the multi-scaleapproach where the friction pressure phase field model is combined with thevalue of the pinning pressure from obtained from smaller scale simulations,one employing the finite element model, and the other the phase field model.9.2.2 The Evolution of CurvatureThe pinning pressure affects the dynamics of the grain boundary networkevolution and alters the kinetics of grain growth. Further, the value of thelimiting grain size is established by the balance between the driving andthe pinning pressures. The analysis the curvature of the grain boundarynetwork reveals quantitative information about the driving pressure in thesystem.Fig. 9.6 shows the curvature distribution and the corresponding mi-crostructures during 2D grain growth under a pinning pressure of Pz/? =0.0087. Initially, in this structure, only 32% of the length of the grainboundary network have a higher curvature than the Pz/? value and is thusable to move. Growth of the grain structure continues until 75% of theinitial grains disappear and the movement of all segments comes to a halt.During this process, two distinct classes appear in the curvature distribution.A large fraction of the grain boundary network?s length possesses zerocurvature as shown in Fig. 9.6c. A smaller peak appears at the limitingcurvature of Pz/?.Fig. 9.7 shows the grain boundary network curvature and curvaturedistribution histograms for three time steps during a 3D grain growthwith Pz/? = 0.015. There is a stark transition from the initial curvaturedistribution in Fig. 9.7a to the distribution with pinning pressure in Fig. 9.7bthat occurs relatively quickly compared with the time required to completelyhalt grain growth. Points on the grain boundary network with highercurvature than Pz/? move; however, they also affect areas with lowercurvature. Figs. 9.7a-b show two points on the grain boundary networkwith arrows (I) and (II), that have higher and lower driving pressure thanPz, respectively. Movement of the grain boundary at point (I) continues119(a) (b) (c)0 0.005 0.01 0.015 0.02 0.025K0 0.005 0.01 0.015 0.02 0.025K 0 0.005 0.01 0.015 0.02 0.025KLength CountLength CountLength CountFigure 9.6: The Grain structure and associated curvature distribu-tions for a 2D system with pinning pressure of Pz = 0.010. (a)The evolving structure at tn = 20000, (b) the partially pinnedstructure at time step 100000 and (c) the final frozen structure.The green line represents the average curvature (Pd/?) and thered line is at Pz/?.due to the high curvature. In comparison, the curvature at point (II)increases as the neighbouring grains disappear and the adjacent triple linemoves, demonstrating that topological transitions propagate through thegrain boundary network.The grain boundary curvature distribution transforms from the self-similar state of ideal grain growth. Therefore, this changes the averagedriving pressure for a given grain size. According to Eqs. 6.3 and 6.4, theaverage curvature of the system and the driving pressure are related through?1 = K?R?. To obtain ?1, first the grain boundary curvature at each pointis calculated. Then, the average curvature of the grain boundary network,K?, is obtained by evaluating Eq. 5.10. For ideal grain growth, ?1 is aconstant (Tab. 6.2); however, with pinning pressure, ?1 changes since thegrain structure deviates from the self-similar structure. Fig. 9.8 shows thevalue of ?1 as a function of R?, measured at different times during the 3Dgrain growth simulations until stagnation is reached.120(a) (b) (c)Figure 9.7: The grain structure and associated curvature distributionsfor 3D system with pinning pressure of Pz = 0.01. (a) the initialstructure, (b) the partially pinned structure at time step 2000and (c) the final frozen structure. The dashed green lines areplaced at K? and the solid red lines are superimposed at thevalue of Pz/? on the curvature distributions.For 3D systems, the value of ?1 changes from 0.27 for the initial idealgrain structure to ?1,lim = 0.5 for the frozen structure. In 2D, ?1 increasesfrom 0.20 for the ideal structure to ?1,lim = 0.26 at the limiting grain size.This increase in ?1 is accompanied by the formation of a peak at Pz/? in thecurvature distribution. The higher increase in 3D is due to the sharper peakin the curvature distribution compared with that found in 2D. In addition,there is no accumulation of flat boundaries in 3D (Figs. 9.6c and 9.7c).Another observation concerns the relationship between the pinning anddriving pressures for the limiting structure. The conventional theories ofZener pinning suggest that the limiting grain size occurs when the pinningpressure is equal to the average driving pressure, i.e. when Pz = P?d,lim =?1?/R?lim [73]. This relationship, however, seems not to be accurate sincethe motion of the entire grain boundary network ceases when the most121Figure 9.8: ?1 as a function of R? during 3D grain growth for threesystems with different pinning pressures. The markers areplotted in the interval of 5000 time steps.curved segment in the system, with the maximum value of Pd, reachesthe Pz value. From Figs. 9.6c and 9.7c, it is evident that for the limitinggrain structure, max(Pd,lim) = Pz and P?d < Pz. Therefore, to quantify therelationship between the driving and pinning pressures, a scaling factor, ?,can be introduced as:P?d,lim = ?Pz. (9.3)The constant ? is calculated from the curvature distributions usingEq. 9.3. For different simulations, the value of ? = 0.34 ? 0.02 is obtainedfor the 2D systems and ? = 0.85 ? 0.03 for the 3D systems. Eqs. 6.4 and9.3 are combined to obtain the limiting grain size as:R?lim =?1,lim??Pz. (9.4)As a consequence, the proportionality of the limiting grain size to thepinning pressure is composed of two constants, ?1,lim and ?, that canbe calculated from the grain boundary curvature distributions. Thesetwo parameters are related to two independent phenomena. ?1,lim is ageometrical/statistical measure that correlates the curvature of the grain122boundary segments with the size of the grains. ? on the other hand, is ameasure of the width of the grain boundary curvature distribution at thelimiting structure. Since a fraction of the grain boundary network alwayspossesses curvatures less than Pz/? at the limiting size, the value of ? isalways less than 1.The R?lim given in Eqs. 9.1 and 9.2 for 2D and 3D, respectively, can becompared with Eq. 9.4. For 2D systems, ?1,lim/? = 0.76 ? 0.16 is close tothe value of 0.73 in Fig. 9.2b. And for 3D systems, ?1,lim/? = 0.56? 0.1 iscomparable with the slope of 0.60 observed in Fig. 9.4. This confirms thatthe ?1 and ? obtained from the curvature measurements produce the samelimiting grain size as is measured by the equivalent radius of the grains.9.2.3 The Kinetics of Grain GrowthAverage Grain Growth RateBased on the above curvature analysis, a physical model can be presentedto describe the overall kinetics of grain growth and the limiting grain size.In accordance with Eq. 6.5 (dR?/dt = ?2MP?d), the rate of growth of R? isproportional to the available driving pressure when the pinning pressure issubtracted:dR?dt= M?2(P?d ? ?Pz)= ?2M[?1?R?? ?Pz](9.5)The scaling constant ? is multiplied by the pinning pressure since at thelimiting size, dR?/dt = 0 and Eq. 9.3 has to hold. The kinetics of graingrowth is obtained by integrating Eq. 9.5:[??1?ln(?1? ? ?PzR?) + ?PzR?2P 2z]R?R?0= ?2Mt. (9.6)R? starts at an initial grain size R?0 and asymptotically reaches the limitinggrain size. Fig. 9.9 shows the average equivalent grain radius as a function oftime for 2D and 3D systems with Pz = 0.010. Simulation results are plottedin comparison with different models. The only model that can describe the123kinetics and limiting grain size accurately is Eq. 9.6, where the constants areobtained from the curvature distribution. For 2D, ?1 = 0.26, ? = 0.34 and?2 = 0.70 were used, while for 3D, ?1 = 0.50, ? = 0.85 and ?2 = 0.80 wereemployed. Although the value of ?1 changes from the initial structure to thelimiting structure, using the constant ?1,lim for the entire curve producessatisfying results. The reason is that, the transition of ?1 from the valuefor the ideal grain growth structure to the value for the structure underpinning pressure, occurs quickly compared with the time required to reachthe limiting grain size (Fig. 9.8).Hunderi-Ryum [70]Andersen-Grong (Eq. 2.38)Hillert(Eq. 2.36)Eq. 9.6Hunderi-Ryum [70]Andersen-Grong (Eq. 2.38)Hillert(Eq. 2.36)Eq. 9.6Figure 9.9: The kinetics of grain growth with pinning compared withthe analytical models. (a) The 2D system with Pz = 0.010,? = 1.155 and M = 2.598. (b) The 3D system with Pz = 0.010,? = 0.667 and M = 3.None of the conventional models is consistent with the simulation data.One reason for the discrepancy is the fact that, to stop grain growth, theaverage driving pressure is not equal to the pinning pressure. But thelimiting grain size is reached when the maximum of the driving pressuredistribution reaches Pz. Failing to account for the distribution of the drivingpressure causes an overestimation of the effect of the pinning pressure ongrain growth. Other reasons for the discrepancies are simplified assumptionsfor the evolution of the grain size distribution or reliance on ideal graingrowth to explain the evolution of a system with pinning pressure.In most experimental studies, grain size can be readily measured.124However, the pinning pressure is usually unknown due to the challengesof accurately measuring particle size and fraction. Therefore, a fittingprocedure is frequently applied to find Pz from the limiting grain size usingEq. 2.38. This procedure can produce a fitting curve on the experiments;however, the estimated value of Pz may not be accurate. For instance, thisprocedure is applied to the current simulation data pretending not to knowthe value of Pz. The rate equation 2.38 is integrated from the initial grainsize to obtain the average grain radius as a function of time. The result ofthis fitting is shown in Fig. 9.10. Based on Andersen?s model, the best fit onthe data is obtained with Pz = 0.0035 and ?1 = 0.21. With this procedure,a pinning pressure is obtained that is three times smaller than the actualvalue. For comparison, a curve with Pz = 0.010 that gives the same limitingsize is also plotted but it does not describe the kinetics correctly since ?1 isunrealistically high. Further, based on measurements by Patterson and Liu[72], Rios [73] suggested ?1 = 0.47. Using this value for ?1 and the actualvalue of Pz = 0.010 leads to a clear underestimation of the limiting grainsize.0 2000 4000 6000 8000 10000 12000 14000 16000202530354045RTimeSimulation with Pz=0.010Pz=0.010?1=0.60Pz=0.0035?1=0.21Pz=0.010?1=0.47Figure 9.10: Fitting Andersen?s equation, Eq. 2.38, on the simulationdata for a 3D system with Pz = 0.010, ? = 0.667 and M =3.The dashed line represents Eq. 9.6.125Individual Grain Growth RateThe pinning pressure reduces the absolute value of the grain growth (orshrinkage) rate. Fig. 9.11 shows the normalized growth rate of the grains asa function of their relative size over four timesteps, from the initial stagesof growth to the total stagnation of the structure. A Weibull distribution isfitted on the grain size data and the resulting grain size distribution curves,fR are superimposed on the graphs.Hillert?s equation (Eq. 2.35) gives the value of dR/dt for grain growthwith a pinning pressure. In this equation, ? = 1 and Rcr are determined byusing the fR curves and solving the volume conservation equation for Rcr,i.e. by inserting Eq. 2.35 into Eq. 2.6:? ?0R2fR[??(1Rcr?1R)? Pz]dR = 0. (9.7)For a comparison with ideal grain growth, Hillert?s relationship (Eq. 2.5)produced a good fit on the scattered simulation data (Fig. 6.6a), whereucr = Rcr/R? = 9/8. However, during grain growth with pinning, ucr changesdue to the lack of self-similarity in the structure. ucr is found to be less than9/8 in the presence of the pinning pressure. The simulation data followsthe trend suggested by Hillert?s model; however, grains in the intermediaterange have higher absolute growth rates compared with those of the Eq. 2.35,suggesting zero growth. As R? approaches R?lim, grain growth rates decreaseand the horizontal section of the curve with dR/dt = 0 extends to higherand lower R/R?.The grain size distribution evolves during grain growth with pinning. Itchanges from the self-similar distribution of ideal grain growth and becomesnarrower and less skewed. Narrowing of the grain size distribution was alsopredicted by Hunderi and Ryum [70], where, they numerically solved thecontinuity equation (Eq. 2.7) with Hillert?s kinetic relationship. However,the shape of the distribution curve here is very different from that of theirwork shown in Fig. 2.9. The truncated shape of the distribution is a resultof the abrupt zero-growth region in Hillert?s equation. The main sourceof discrepancy here is the mean field assumption of Hillert?s relationship,126Figure 9.11: The normalized growth rate of individual grains as afunction of their relative size for 3D simulations with Pz =0.010 at time steps (a) tn = 100 (b) tn = 3000 (c) tn = 36000and (d) tn = 136000. The solid red lines represent Hillert?sequation with ? = 1 and the given ucr. The green curvesrepresent grain size distributions.where it supposes that a grain is surrounded by grains with the radius Rcr.In contrast, in the numerical simulation, grains that have a radius close toRcr still evolve depending on their local topological environment.9.2.4 The Topological Aspects of Grain GrowthAs shown in Fig. 6.2, the 2D ideal system closely follows the von Neumann-Mullins relationship. For a system with pinning pressure, an equivalentrelationship can be derived by modifying Eq. 2.13. The pining pressure127opposes the driving pressure of a curved grain boundary. Therefore, basedon Eq. 2.12, the rate of change in the area of a 2D grain is:dAdt= ?M?s(?K ? Pz) ds = ?M?(2pi ?pi3N ?Pzs?), (9.8)Here, the ? sign depends on the value of N . For grains with concavesegments (N < 6), the minus sign is considered, and for grains with convexcurvature (N > 6), the plus sign is taken into account. s is the perimeterof the grain and can be written as a function of the grain radius. Byapproximating s = 2piR, Eq. 9.8 can be rearranged as:RdRdt=M?6(N ? 6?6Pz?R)(9.9)The plus sign is used when N < 6 ? 6PzR? and the negative sign isapplied when N > 6 + 6PzR? . For N values in this range, dR/dt = 0. Therelationship between the size of a grain and its number of faces depends onthe morphology of the grain structure. The following relationship is obtainedfrom the ideal grain structure by curve fitting:R = R?(0.20N ? 0.14). (9.10)As shown in Fig. 9.1, the pinning pressure has little influence on theoverall shape of the grains. Therefore, R on the right hand side of Eq. 9.9can be replaced by Eq. 9.10 to obtain a relationship between growth rateand the number of faces.Fig. 9.12 shows the normalized growth rate of individual grains as afunction of their topological class in 2D. In the presence of particle pinning,the absolute value of RdR/dt decreases for all topological classes. As graingrowth continues toward the frozen structure, the growth rates of grainsclose to the six-sided class are influenced more than that of the grains thatare topologically further away from the central class. The reason is that theaverage curvature of grain boundary segments is zero for six-sided grainsand is higher for grains with greater or lesser number of faces.128Figure 9.12: The normalized growth rate as a function of the numberof faces of a grain for a 3D system with Pz/? = 0.0075 duringgrain growth until stagnation. The square markers representthe average value of the growth rate for each topological classand the solid lines represent Eq. 9.9.The normalized growth rates according to Eq. 9.9 are plotted in Fig. 9.12for comparison. The modified von Neumann-Mullins equation follows thegeneral trend of the simulation data, however, there is a large deviation atearlier timesteps where a higher driving pressure is available for the motionof the grain boundaries.A key assumption in the derivation of the von Neumann-Mullinsrelationship and its modified version is that the grain boundaries have theequilibrium angle of 120? at the triple junctions. It should be noted since129the pinning pressure in the friction pressure model does not change theorder parameter profile (Appendix A) and the grain boundary energy, thetriple junction angle remains also 120? with the pinning pressure and this isconfirmed in the present phase field simulations.For 3D grain growth, a similar theory can be proposed by modifyingMullins (Eq. 2.17) or Hilgenfeldt et al. (Eq. 2.21) relationships. For asystem with pinning pressure, the growth rate of a grain is:dVdt= M?? (K ?Pz?)dS. (9.11)The two sides of this equation are multiplied by V ?13 and can becompared with the left hand side of Eq. 2.21, i.e. V ?1/3dV/dt = M?GH(N).The integral term can be replaced by the GH(N) function, resulting in:V ?13dVdt= M?(GH(N)??RPz?). (9.12)Here, ?R = S/V13 , is a metric measure and depends on the shape andsize of the grain. For instance, for a sphere, ?R = (4pi/(4/3pi)13 )R = 7.79R.The ? sign is negative when GH(N) ? ?RPz/? > 0 and is positive whenGH(N) + ?RPz/? < 0.Calculations using Eq. 9.12 are now compared with the simulation data.Fig. 9.13 shows the normalized growth rate of individual grains as a functionof their topological class for a system with Pz = 0.010. The value of ?R =5.0R is chosen as a best fit to the scatter data for all timesteps. Sincethe effect of pinning pressure on the shape of the grains is negligible, therelationship between R and N is obtained by a cubic polynomial fit on thesimulation data for ideal grain growth:?R = 5.0? R?(8.1? 10?5N3 ? 6.6? 10?3N2 + 0.22N ? 0.74). (9.13)Eq. 9.12 follows the simulation data reasonably well. Here, a fittingprocedure was used to obtain ?R = 5.0R and to relate R to N . However,130Figure 9.13: The normalized growth rate of individual grains as afunction of the number of faces for the grain growth with Pz =0.01. The solid lines represent Eq. 9.12. Eq. 2.21 for idealgrain growth is plotted in (a) for comparison.such relationships can also be derived by simplifying geometrical assump-tions. Nonetheless, Eq. 9.12 clearly establishes a topological relationship forgrain growth with the effect of pinning pressure and paves the way for afundamental analytical model to describe the kinetics of grain growth withpinning.9.3 Grain Growth with Solute DragThe effects of solute drag on grain growth in 2D is studied by using thefriction pressure phase field model (Eq. 8.12). Grain growth simulations are131performed with m = 2, ? = 4, L = 1, ?x = 1 and ?t = 0.05. The systemsize is 2000 ? 2000 grid points. Cahn?s expression for solute drag pressure,Eq. 2.43 is used. Snapshots of the structure during grain growth are shownin Fig.9.14 for a system with a = 1 and b = 4000. Since the drag pressureis velocity dependent and reaches zero when the velocity is zero, no limitinggrain size was observed, and the structural coarsening continues at a slowerrate.Figure 9.14: Snapshots of 2D grain growth with solute drag pressure(a = 1, b = 4000 in Eq. 2.43).To study the effects of the solute drag pressure on the kinetics of graingrowth, a series of simulations with a = 1 and b varying from 200 to 4? 105are performed. The parameter b determines the driving pressure where thetransition from low to high velocity drag occurs. With lower b values, higherdriving pressures are required for the transition (Fig. 9.15a). The square ofthe average grain radius is plotted against time in Fig. 9.15b. An idealgrain growth curve in a system with the same properties is also plotted forcomparison. All curves start at an initial grain size, R0, coinciding with theideal grain growth line. The curves deviate from the ideal grain growth lineand do not show any sign of saturation. This behaviour can be describedphenomenologically by introducing a time exponent, ? lower than 0.5, forthe grain growth kinetics [2], i.e.,R?1/? ? R?01/? = M?kst, ? ? 0.5, (9.14)132where ks is a fitting parameter. For each simulation, the value of thetime exponent ? and ks are obtained by a non-linear fitting of Eq. 9.14on the average equivalent radius versus time data. The straight solid linesin Fig. 9.15 represent parabolic grain growth with ? = 0.5.Figure 9.15: The effect of solute drag pressure on the kinetics ofgrain growth for simulations with different b values in Cahn?sequation. (a) The normalized grain boundary velocity as afunction of the normalized driving pressure (K = Pd/?). Thecurvature distribution for initial grain size is superimposed onthe Pd/? axis. (b) The kinetics of grain growth. The solidlines represent the parabolic kinetics.Lower ? values indicate deviation from parabolic grain growth. Forsimulations with high b (e.g. b = 2?104), the time exponent is close to thatof ideal grain growth. However, as b decreases, the deviation from parabolicgrain growth increases. For simulation with b = 4000, the maximumdeviation is observed with ? = 0.35?0.01. Lowering the b values further (e.g.b = 200) reverses the trend and the grain growth becomes parabolic againwith a very low kinetic coefficient, k. Fig. 9.16 illustrates the relationshipbetween the grain growth exponent, ? and the solute drag parameter, b.The minimum value of the exponent is positioned at intermediate values ofb.Grain boundary velocity curves in Fig. 9.15a can be compared with that133Figure 9.16: The time exponent of grain growth (? in Eq. 9.14) as afunction of the solute drag parameter, b in Eq. 2.43.of grain growth kinetics in Fig. 9.15b. The curvature distribution (K =Pd/?) of the initial grain structure with K? = 0.013 is also superimposedon Fig. 9.15a. Different points on the grain boundary network experience adifferent drag pressures depending on the local driving pressure. In systemswith intermediate values of b (e.g. b = 4000), the majority of the grainboundary network has a driving pressure that falls in the transition regionfrom the high to the low velocity drag regime during grain growth. However,for curves with b = 200 or b = 2 ? 104 the major part of the length of thegrain boundary network is in the low or high velocity regimes, respectively.As grain growth continues, the curvature distribution moves towards lowercurvature values and the local velocity of the grain boundary follows thecurves in Fig. 9.15a. It is evident that the transition of velocity fromhigh to low solute drag regimes translates into grain growth kinetics withnon-parabolic behaviour. However, for very weak or strong solute drag,the majority of the boundary segments remain in the high or low velocityregimes, respectively, such that the grain growth kinetics is parabolic.1349.4 ConclusionsIn this chapter, the friction pressure phase field method was used tosimulate grain growth with the presence of a uniform pinning pressure.This simulation technique allows the pinning pressure obtained from smallscale models to be incorporated into large scale grain growth simulations.Relationships for Pz were available from a finite element model by Couturieret al. [8] and the small-scale phase field model presented in Chapter 7.Combining the results from the two length scales results in a limiting grainsize that coincides with the experimental data.The curvature distribution of the grain boundary network is furtherinvestigated to obtain the driving pressure of grain growth. It was shownthat, at the limiting structure, the average driving pressure for grain growthis not equal to but smaller than the pinning pressure (P?d = ?Pz). Thescaling factor ? = 0.34 for 2D and ? = 0.85 for 3D systems is independentof the pinning pressure in the investigated range. The general relationshipfor the limiting grain size (Eq. 9.4) obtained from curvature analysis can beapplied to any system with isotropic grain boundary properties. The ratioof ?1,lim/? obtained from the curvature analysis in 3D agrees well with theresults deduced from the simulated microstructures and is thus consistentwith the wide range of experimental data listed in Manohar?s review.All the grain boundaries in this phase field model have equal grainboundary energy and mobility. This is not always valid as physical systemscan have anisotropic grain boundary properties. The experimental studiesreported in Manohar?s review belong to steels and aluminium alloys. Thesesystems usually exhibit low to moderate amount of anisotropy. Nonetheless,it is not clear what would be the exact influence of anisotropy on thelimiting grain size and the grain growth kinetics. The friction pressure modelis indeed a first approximation for modelling complex physical systems.The good agreement between the simulation results and the experimentssuggests, however, that the effect of anisotropy in these systems is secondaryto the important geometrical corrections that has been suggested byconsidering curvature distribution of an isotropic grain boundary network.135Nevertheless, grain growth simulations should be extended to systems withsignificant anisotropies to quantify the role of pinning in these systems.Another core result of this study is related to developing an accuratedescription of grain growth kinetics, while conventional models for particlepinning fail to capture it. By analyzing the curvature of grain boundarynetwork, it was shown that the rate of grain growth, dR?/dt, is proportionalto P?d??Pz, with the proportionality constant ?2 obtained from ideal graingrowth. These parameters can be applied to isotropic systems with uniformpinning pressure in which the initial structure is close to the ideal grainstructure. Since the parameters are functions of the grain size and curvaturedistribution, they may not be generalized to systems with large variationin grain size, bimodal grain size distribution or inhomogeneous particledispersion. For such systems, one can use the proposed phase field modelto simulate the structure starting from a desired grain size structure withspatial and temporal variation in the pinning pressure.A fundamental step towards developing the kinetic theory of graingrowth with pinning is undertaken by proposing a topological relationshipfor the individual growth rate of grains. The modified Hilgenfeldt equationgiven here (Eq. 9.12) produces comparable results with the numericalsimulation. Furthermore, this study demonstrates that the evolution ofgrain boundary networks under uniform pinning pressure is fundamentallydifferent in 2D or 3D due to the different curvature distributions. Therefore,one can not obtain quantitative information about 3D systems from 2Dsimulations.Application of the friction pressure phase field model is also shown forgrain growth with solute drag in 2D. A phenomenological power-law kineticrelationship was used to fit the simulation data and to obtain the graingrowth time exponent, ?. A time exponent of 0.5, associated with idealgrain growth, has rarely been observed in materials [2]. The effect of solutedrag pressure on the motion of grain boundaries can justify the reduction inthe time exponent as can be seen in Fig. 9.16. Experimental observations ofgrain growth kinetics in materials with a small amount of impurities showtime exponents of ? = 0.2 ? 0.4 [2], which is consistent with the trends136of the simulations. These results need to be further confirmed with three-dimensional grain growth simulations.137Chapter 10Overall Conclusions10.1 Achievements and ContributionsThis work demonstrates a multi-scale modelling scheme for the simulationof interface drag phenomena in polycrystalline materials. In nano-scalesimulations, the interaction of a grain boundary with a particle distributionwas analyzed. The resulting pinning pressures were combined with theresults of a meso-scale model to describe the limiting grain size in materialswith particle pinning. The devised meso-scale simulation technique, theso called ?friction pressure? phase field model, is capable of incorporatinga drag pressure and simulating the evolution of the grain structure.Application of this model was also shown for grain growth with solutedrag according to Cahn?s model for the interaction of solute atoms and agrain boundary. With this model, the influence of any constant or velocity-dependent drag or the combination of solute drag with particle pinning ona grain structure, can be simulated.Analyzing the distribution of curvature of a grain boundary network pro-vides a compelling insight into the evolution of structure. The distribution ofdriving pressure of grain growth was measured with this technique for bothideal grain growth and grain growth with pinning pressure. This informationrevealed that, at the limiting grain structure, the average driving pressurewas not equal to the pinning pressure, as conventionally assumed, but the138maximum of the driving pressure distribution equals the pinning pressure.The kinetic constants of grain growth that were obtained from the curvatureanalysis can be used in future experiments to back-calculate a more accuratevalue for the pinning pressure compared with that of the mean field models.The present friction pressure model implies a generalized view of theinterface drag phenomena. Particle pinning and solute drag are not disparateprocesses, but are caused by the interaction of a moving interface withobstacles. Such obstacles can be small and mobile such as solute atomsor large and immobile such as second phase particles. This spectrum ofobstacles also encompasses solute clusters and very fine mobile particles,that have not yet been studied extensively, but is believed to be importantsources of drag. One can use the present multi-scale approach to firstsimulate the interaction of clusters with a grain boundary to obtain thevelocity-dependent drag pressure and then simulate the evolution of themicrostructure with the friction pressure model.The main engineering contribution of this work is on the correlation ofthe process parameters to the microstructure of a polycrystalline material.In applications where an accurate control of grain size is critical, thegeometrical/statistical constants derived from this work can be used, in theproposed physical model without further simulations, to estimate processparameters such as annealing time and temperature. Furthermore, thismulti-scale simulation methodology can be utilized for materials where morecomplicated pinning scenarios, such as particle dissolution and combinationof solute drag and particle pinning are involved.10.2 Future WorkAs a theoretical study, this work has provides a basis for further investigationof the effect of drag phenomena in the microstructural evolution ofpolycrystalline materials. Future projects can be defined to extend theapplications of this work to systems with direct technological applications.On a theoretical front, many questions have emerged regarding a generalizedtheory of interface drag and abnormal grain growth.139Advanced 3D microstructural characterization techniques are becomingmore feasible in materials science with the development of such techniquesas tomography, serial sectioning and 3D electron back-scattered diffraction.These methods enhance our understanding of the microstructure by pro-viding a way to measure detailed microstructural attributes and to replacethe traditional mean field measures often obtained from 2D cross-sections.Microstructural modelling can lead this paradigm shift by demonstratingthe applicability of such characterizations. The analysis of the distributionof curvature of the grain boundary network is an example.A few suggestions for future projects that directly benefit from this workis given below:Experiments on the kinetics of grain growth:The consistency of the values for parameters ?1 and ? with experimentalresults on limiting grain size has been demonstrated. The proposed kineticparameters of grain growth with pinning, together with ?2, needs to betested on a particular alloy system with stable, non-soluble particles tocompare with the kinetics of grain growth. To obtain the kinetic curve, thelaser-ultrasonic grain size measurement technique is the method of choice tobe applied to a metallic system with dispersed inert particles.The combined effect of solute drag and particle pinning:In many alloys with partial immiscibility, solute atoms both form in-termetallic particles and segregate at the grain boundaries. Duringthermal treatment, particles dissolve or precipitate, creating a transitionfrom predominantly one form of drag to another. This process may beaccompanied by drastic microstructural changes such as abnormal graingrowth. The friction pressure model allows both sources of drag to beimplemented on a structure and to simulate grain growth structures duringheating or cooling.140Investigations into abnormal grain growth:The explanations for the formation of abnormally large grains in materialswith second phase particles has remained controversial. The reason isthat the nucleation and growth of an abnormal grain depends on the localmorphology of the grains and the topological aspects of grain growth. Thefriction pressure model provides a tool to simulate a grain microstructurewith spatially and temporally variable pinning pressures. The curvatureand topological analysis will be used to investigate, in detail, the conditionswhich trigger nucleation of an abnormal grain and those which sustain itsgrowth. In addition, the anisotropy of the grain boundary network andlow-angle grain boundaries can be included in the phase field model.141Bibliography[1] J. W. Martin, R. D. Doherty, and B. Cantor, Stability ofMicrostructure in Metallic Systems. Cambridge University Press,2 ed., 1997. ? pages 1, 19, 21, 33[2] F. J. Humphreys and M. Hatherly, Recrystallization and RelatedAnnealing Phenomena, Second Edition. Pergamon, 2 ed., 2004. ?pages 1, 6, 19, 20, 37, 84, 132, 136[3] T. Gladman, The physical metallurgy of microalloyed steels. Instituteof Materials, 1997. ? pages 1[4] L. M. Fu, H. R. Wang, W. Wang, and A. D. Shan, ?Austenite graingrowth prediction coupling with drag and pinning effects in lowcarbon Nb microalloyed steels,? Materials Science and Technology,vol. 27, no. 6, pp. 996?1001, 2011. ? pages 2[5] C. Zener, ?Private communication to C. S. Smith, Transactions ofAIME,? vol. 175, pp. 15?51, 1948. ? pages 2, 19, 118[6] M. Hillert, ?Inhibition of grain growth by second-phase particles,?Acta Metallurgica, vol. 36, no. 12, pp. 3177?3181, 1988. ? pages 2,22, 27, 98[7] C. Wo?rner and A. Cabo, ?On the grain growth inhibition by secondphases particles,? Acta Metallurgica, vol. 35, no. 11, pp. 2801?2804,1987. ? pages 2, 22[8] G. Couturier, R. Doherty, C. Maurice, and R. Fortunier, ?3D finiteelement simulation of the inhibition of normal grain growth byparticles,? Acta Materialia, vol. 53, no. 4, pp. 977?989, 2005. ?pages x, 2, 29, 30, 31, 90, 91, 98, 118, 135142[9] N. Moelans, B. Blanpain, and P. Wollants, ?A phase field model forthe simulation of grain growth in materials containing finelydispersed incoherent second-phase particles,? Acta Materialia,vol. 53, no. 6, pp. 1771?1781, 2005. ? pages x, 2, 28, 53, 54, 55[10] P. A. Manohar, M. Ferry, and T. Chandra, ?Five decades of theZener equation,? ISIJ International, vol. 38, no. 9, pp. 913?924,1998. ? pages 2, 21, 22, 23, 118[11] M. Miodownik, E. A. Holm, and G. N. Hassold, ?Highly parallelcomputer simulations of particle pinning: zener vindicated,? ScriptaMaterialia, vol. 42, no. 12, pp. 1173?1177, 2000. ? pages ix, 2, 27,28, 102, 118[12] Y. Suwa, Y. Saito, and H. Onodera, ?Phase field simulation of graingrowth in three dimensional system containing finely dispersedsecond-phase particles,? Scripta Materialia, vol. 55, no. 4,pp. 407?410, 2006. ? pages ix, 2, 28, 92, 102, 118[13] J. W. Cahn, ?The impurity-drag effect in grain boundary motion,?Acta Metallurgica, vol. 10, no. 9, pp. 789?798, 1962. ? pages 2, 34,35[14] M. Hillert and B. Sundman, ?A treatment of the solute drag onmoving grain boundaries and phase interfaces in binary alloys,? ActaMetallurgica, vol. 24, no. 8, p. 3436, 1976. ? pages 2, 36, 38[15] L. A. Barrales-Mora, 2D and 3D Grain Growth Modeling andSimulation. Cuvillier Verlag, 2008. ? pages 2[16] M. Elsey, S. Esedolu, and P. Smereka, ?Large-scale simulation ofnormal grain growth via diffusion-generated motion,? Proceedings ofthe Royal Society A: Mathematical, Physical and EngineeringScience, vol. 467, no. 2126, pp. 381 ?401, 2011. ? pages ix, 2, 17, 18,83[17] E. A. Lazar, J. K. Mason, R. D. MacPherson, and D. J. Srolovitz,?A more accurate three-dimensional grain growth algorithm,? ActaMaterialia, vol. 59, no. 17, pp. 6837?6847, 2011. ? pages 2, 16, 17,18, 19, 51, 53, 83[18] J. Kepler, Strena seu: De nive sexangula. Ostwalds Klassiker derexakten Wissenschaften (1987), Geest & Portig, 1. aufl ed., 1611. ?pages 5143[19] G. F. V. Voort, Metallography: past, present, and future: (75thanniversary volume). ASTM International, 1993. ? pages 5[20] J. Burke and D. Turnbull, ?Recrystallization and grain growth,?Progress in Metal Physics, vol. 3, no. 0, pp. 220?292, 1952. ? pages6, 7[21] P. Feltham, ?Grain growth in metals,? Acta Metallurgica, vol. 5,no. 2, pp. 97?105, 1957. ? pages 7[22] M. Hillert, ?On the theory of normal and abnormal grain growth,?Acta Metallurgica, vol. 13, no. 3, pp. 227?238, 1965. ? pages 7, 9,24, 84[23] I. Lifshitz and V. Slyozov, ?The kinetics of precipitation fromsupersaturated solid solutions,? Journal of Physics and Chemistry ofSolids, vol. 19, no. 1?2, pp. 35?50, 1961. ? pages 9[24] N. Louat, ?The resistance to normal grain growth from a dispersionof spherical particles,? Acta Metallurgica, vol. 30, no. 7,pp. 1291?1294, 1982. ? pages 9, 22[25] D. Aboav and T. Langdon, ?The distribution of grain diameters inpolycrystalline magnesium oxide,? Metallography, vol. 1, no. 3?4,pp. 333?340, 1969. ? pages 10[26] O. Hunderi and N. Ryum, ?The kinetics of normal grain growth,?Journal of Materials Science, vol. 15, no. 5, pp. 1104?1108, 1980. ?pages 10[27] P. Rios, ?Comparison between a grain size distribution obtained by aMonte Carlo Potts model and by an analytical mean field model,?Scripta Materialia, vol. 41, no. 12, pp. 1283?1287, 1999. ? pages 10[28] P. Rios, ?Comparison between a computer simulated and ananalytical grain size distribution,? Scripta Materialia, vol. 40, no. 6,pp. 665?668, 1999. ? pages 18[29] P. R. Rios and K. Lu?cke, ?Comparison of statistical analyticaltheories of grain growth,? Scripta Materialia, vol. 44, no. 10,pp. 2471?2475, 2001. ? pages 10[30] J. von Neumann in Metals Interfaces, (Cleveland), p. 108, AmericanSociety for Metals, 1952. ? pages 11144[31] W. W. Mullins, ?TwoDimensional Motion of Idealized GrainBoundaries,? Journal of Applied Physics, vol. 27, no. 8, pp. 900?904,1956. ? pages 11[32] W. W. Mullins, ?Grain growth of uniform boundaries with scaling,?Acta Materialia, vol. 46, no. 17, pp. 6219?6226, 1998. ? pages 12, 13[33] T. Saetre, O. Hunderi, and N. Ryum, ?Modelling grain growth intwo dimensions,? Acta Metallurgica, vol. 37, no. 5, pp. 1381?1387,1989. ? pages 13[34] T. O. Saetre, ?On the theory of normal grain growth in twodimensions,? Acta Materialia, vol. 50, no. 6, pp. 1539?1546, 2002. ?pages[35] T. O. Saetre, ?A mean field theory of grain growth: I,? Modellingand Simulation in Materials Science and Engineering, vol. 12, no. 6,pp. 1267?1277, 2004. ? pages[36] T. O. Saetre, ?A mean field theory of grain growth: II,? Modellingand Simulation in Materials Science and Engineering, vol. 12, no. 6,pp. 1279?1292, 2004. ? pages 13[37] W. Mullins, ?Estimation of the geometrical rate constant in idealizedthree dimensional grain growth,? Acta Metallurgica, vol. 37, no. 11,pp. 2979?2984, 1989. ? pages 13, 14, 51[38] S. Hilgenfeldt, A. M. Kraynik, S. A. Koehler, and H. A. Stone, ?AnAccurate von Neumann?s Law for Three-Dimensional Foams,?Physical Review Letters, vol. 86, no. 12, pp. 2685?2688, 2001. ?pages 14, 51[39] R. D. MacPherson and D. J. Srolovitz, ?The von Neumann relationgeneralized to coarsening of three-dimensional microstructures,?Nature, vol. 446, no. 7139, pp. 1053?1055, 2007. ? pages 14, 15[40] P. R. Rios, M. E. Glicksman, and D. Lewis, ?Advances in GrainGrowth Theory,? in Recrystallization and Grain Growth Iv (E. J.Palmiere and B. P. Wynne, eds.), vol. 715-716, pp. 211?218,Stafa-Zurich: Trans Tech Publications Ltd, 2012. ? pages 15[41] E. A. Lazar, J. K. Mason, R. D. MacPherson, and D. J. Srolovitz,?Complete Topology of Cells, Grains, and Bubbles in145Three-Dimensional Microstructures,? Physical Review Letters,vol. 109, no. 9, p. 095505, 2012. ? pages 15[42] E. B. Matzke, ?The Three-Dimensional Shape of Bubbles inFoam-An Analysis of the Role of Surface Forces inThree-Dimensional Cell Shape Determination,? American Journal ofBotany, vol. 33, no. 1, pp. 58?80, 1946. ? pages 15[43] L. Barrales Mora, G. Gottstein, and L. Shvindlerman,?Three-dimensional grain growth: Analytical approaches andcomputer simulations,? Acta Materialia, vol. 56, no. 20,pp. 5915?5926, 2008. ? pages 15[44] M. Anderson, D. Srolovitz, G. Grest, and P. Sahni, ?Computersimulation of grain growth?I. Kinetics,? Acta Metallurgica, vol. 32,no. 5, pp. 783?791, 1984. ? pages 16[45] P. S. Sahni, G. S. Grest, M. P. Anderson, and D. J. Srolovitz,?Kinetics of the Q-State Potts Model in Two Dimensions,? PhysicalReview Letters, vol. 50, no. 4, pp. 263?266, 1983. ? pages 16[46] P. S. Sahni, D. J. Srolovitz, G. S. Grest, M. P. Anderson, and S. A.Safran, ?Kinetics of ordering in two dimensions. II. Quenchedsystems,? Physical Review B, vol. 28, no. 5, pp. 2705?2716, 1983. ?pages 16[47] D. Srolovitz, M. Anderson, P. Sahni, and G. Grest, ?Computersimulation of grain growth?II. Grain size distribution, topology, andlocal dynamics,? Acta Metallurgica, vol. 32, no. 5, pp. 793?802, 1984.? pages 16, 27[48] M. P. Anderson, G. S. Grest, and D. J. Srolovitz, ?Computersimulation of normal grain growth in three dimensions,?Philosophical Magazine Part B, vol. 59, no. 3, pp. 293?329, 1989. ?pages 16, 19[49] K. Fuchizaki, T. Kusaba, and K. Kawasaki, ?Computer modelling ofthree-dimensional cellular pattern growth,? Philosophical magazine.B. Physics of condensed matter. Structural, electronic, optical andmagnetic properties, vol. 71, no. 3, pp. 333?357, 1995. ? pages 16[50] D. Weygand, Y. Bre?chet, and J. Le?pinoux, ?A vertex dynamicssimulation of grain growth in two dimensions,? PhilosophicalMagazine Part B, vol. 78, no. 4, pp. 329?352, 1998. ? pages 17, 82146[51] D. Weygand, Y. Bre?chet, and J. Le?pinoux, ?A Vertex Simulation ofGrain Growth in 2D and 3D,? Advanced Engineering Materials,vol. 3, no. 1-2, pp. 67?71, 2001. ? pages 16[52] K. A. Brakke, ?The surface evolver,? Experimental Mathematics,vol. 1, no. 2, pp. 141?165, 1992. ? pages 16[53] F. Wakai, N. Enomoto, and H. Ogawa, ?Three-dimensionalmicrostructural evolution in ideal grain growth?general statistics,?Acta Materialia, vol. 48, no. 6, pp. 1297?1311, 2000. ? pages 16, 17,19[54] E. A. Lazar, R. D. MacPherson, and D. J. Srolovitz, ?A moreaccurate two-dimensional grain growth algorithm,? Acta Materialia,vol. 58, no. 2, pp. 364?372, 2010. ? pages 16, 18, 82[55] B. Merriman, J. K. Bence, and S. J. Osher, ?Motion of MultipleJunctions: A Level Set Approach,? Journal of ComputationalPhysics, vol. 112, no. 2, pp. 334?363, 1994. ? pages 17[56] M. Elsey, S. Esedoglu, and P. Smereka, ?Diffusion generated motionfor grain growth in two and three dimensions,? Journal ofComputational Physics, vol. 228, no. 21, pp. 8015?8033, 2009. ?pages 17[57] N. Moelans, F. Wendler, and B. Nestler, ?Comparative study of twophase-field models for grain growth,? Computational MaterialsScience, vol. 46, no. 2, pp. 479?490, 2009. ? pages 17, 51, 82[58] D. Zo?llner, ?Grain microstructure evolution in two-dimensionalpolycrystals under limited junction mobility,? Scripta Materialia,vol. 67, no. 1, pp. 41?44, 2012. ? pages 18, 82[59] R. Darvishi Kamachali and I. Steinbach, ?3-D phase-field simulationof grain growth: Topological analysis versus mean-fieldapproximations,? Acta Materialia, vol. 60, no. 6?7, pp. 2719?2728,2012. ? pages 18, 51, 53, 82[60] S. G. Kim, D. I. Kim, W. T. Kim, and Y. B. Park, ?Computersimulations of two-dimensional and three-dimensional ideal graingrowth,? Physical Review E (Statistical, Nonlinear, and Soft MatterPhysics), vol. 74, no. 6, pp. 061605?14, 2006. ? pages 18, 51, 82147[61] D. Zo?llner and P. Streitenberger, ?Three-dimensional normal graingrowth: Monte Carlo Potts model simulation and analytical meanfield theory,? Scripta Materialia, vol. 54, no. 9, pp. 1697?1702, 2006.? pages 18, 82[62] C. E. Krill and L. Q. Chen, ?Computer simulation of 3-D graingrowth using a phase-field model,? Acta Materialia, vol. 50, no. 12,pp. 3059?3075, 2002. ? pages 19, 41, 50, 52, 61[63] Y. Saito, ?Monte Carlo Simulation of Grain Growth inThree-dimensions,? ISIJ International, vol. 38, no. 6, pp. 559?566,1998. ? pages 19[64] D. Weygand, Y. Bre?chet, and J. Lepinoux, ?Reduced mobility oftriple nodes and lines on grain growth in two and three dimensions,?Interface Science, vol. 7, no. 3, pp. 285?295, 1999. ? pages 19, 23, 26[65] T. Senuma, ?Present status of and future prospects for precipitationresearch in the steel industry,? ISIJ International, vol. 42, no. 1,pp. 1?12, 2002. ? pages 19[66] M. F. Ashby, J. Harper, and J. Lewis, ?Interaction of CrystalBoundaries with Second-Phase Particles,? Transactions of theMetallurgical Society of Aime, vol. 245, no. 2, pp. 413?&, 1969. ?pages 20, 86[67] P. Hellman and M. Hillert, ?On the effect of second phase particleson grain growth,? Scandinavian Journal of Metallurgy, vol. 4, no. 5,pp. 211?219, 1975. ? pages 21, 91[68] C. Wo?rner, ?Some remarks on the zener drag,? Scripta Metallurgica,vol. 23, no. 11, pp. 1909?1912, 1989. ? pages 22[69] P. Rios, ?Overview no. 62: A theory for grain boundary pinning byparticles,? Acta Metallurgica, vol. 35, no. 12, pp. 2805?2814, 1987. ?pages 23[70] O. Hunderi and N. Ryum, ?On the stagnation of grain growth,? ActaMetallurgica, vol. 30, no. 3, pp. 739?742, 1982. ? pages ix, 24, 25,124, 126[71] I. Andersen and ?. Grong, ?Analytical modelling of grain growth inmetals and alloys in the presence of growing and dissolving148precipitates?I. Normal grain growth,? Acta Metallurgica etMaterialia, vol. 43, no. 7, pp. 2673?2688, 1995. ? pages 25, 26[72] B. Patterson and Y. Liu, ?Relationship between grain BoundaryCurvature and Grain Size,? Metallurgical and Materials TransactionsA, vol. 23, no. 9, pp. 2481?2482, 1992. ? pages 25, 71, 125[73] P. R. Rios, ?On the relationship between pinning force and limitinggrain radius,? Scripta Materialia, vol. 34, no. 8, pp. 1185?1188, 1996.? pages 25, 121, 125[74] M. Maalekian, R. Radis, M. Militzer, A. Moreau, and W. Poole, ?Insitu measurement and modelling of austenite grain growth in aTi/Nb microalloyed steel,? Acta Materialia, vol. 60, no. 3,pp. 1015?1026, 2012. ? pages 26[75] D. Srolovitz, M. Anderson, G. Grest, and P. Sahni, ?Computersimulation of grain growth-III. Influence of a particle dispersion,?Acta Metallurgica, vol. 32, no. 9, pp. 1429?1438, 1984. ? pages 26[76] J. Gao, R. G. Thompson, and B. R. Patterson, ?Computersimulation of grain growth with second phase particle pinning,? ActaMaterialia, vol. 45, no. 9, pp. 3653?3658, 1997. ? pages 26[77] N. Moelans, B. Blanpain, and P. Wollants, ?Phase field simulationsof grain growth in two-dimensional systems containing finelydispersed second-phase particles,? Acta Materialia, vol. 54, no. 4,pp. 1175?1184, 2006. ? pages 26, 27[78] G. Couturier, C. Maurice, and R. Fortunier, ?Three-dimensionalfinite-element simulation of Zener pinning dynamics,? PhilosophicalMagazine, vol. 83, no. 30, pp. 3387?3405, 2003. ? pages 29[79] M. Apel, B. Boettger, J. Rudnizki, P. Schaffnit, and I. Steinbach,?Grain Growth Simulations Including Particle Pinning Using theMultiphase-field Concept,? ISIJ International, vol. 49, no. 7, Sp. Iss.SI, pp. 1024?1029, 2009. ? pages x, 31, 32[80] I. Toda-Caraballo, C. Capdevila, G. Pimentel, and C. De Andre?s,?Drag effects on grain growth dynamics,? Computational MaterialsScience, vol. 68, pp. 95?106, 2013. ? pages 32149[81] K. Lu?cke and K. Detert, ?A quantitative theory of grain-boundarymotion and recrystallization in metals in the presence of impurities,?Acta Metallurgica, vol. 5, no. 11, pp. 628?637, 1957. ? pages 33[82] M. Hillert, ?Solute drag, solute trapping and diffusional dissipationof Gibbs energy,? Acta Materialia, vol. 47, no. 18, pp. 4481?4505,1999. ? pages 36[83] M. I. Mendelev and D. J. Srolovitz, ?A regular solution model forimpurity drag on a migrating grain boundary,? Acta Materialia,vol. 49, no. 4, pp. 589?597, 2001. ? pages 36[84] Y. J. M. Bre?chet and G. R. Purdy, ?Solute drag in ternary solidsolutions,? Acta Materialia, vol. 51, no. 18, pp. 5587?5592, 2003. ?pages 36[85] M. I. Mendelev and D. J. Srolovitz, ?Grain-boundary migration inthe presence of diffusing impurities: simulations and analyticalmodels,? Philosophical Magazine A, vol. 81, no. 9, pp. 2243?2269,2001. ? pages 37[86] M. Greenwood, C. Sinclair, and M. Militzer, ?Phase field crystalmodel of solute drag,? Acta Materialia, vol. 60, no. 16,pp. 5752?5761, 2012. ? pages 37[87] K. Gro?nhagen and J. A?gren, ?Grain-boundary segregation anddynamic solute drag theory?A phase-field approach,? ActaMaterialia, vol. 55, no. 3, pp. 955?960, 2007. ? pages xi, 37, 57, 58[88] S. G. Kim and Y. B. Park, ?Grain boundary segregation, solute dragand abnormal grain growth,? Acta Materialia, vol. 56, no. 15,pp. 3739?3753, 2008. ? pages 37[89] J. Li, J. Wang, and G. Yang, ?Phase field modeling of grainboundary migration with solute drag,? Acta Materialia, vol. 57,no. 7, pp. 2108?2120, 2009. ? pages 37, 58[90] H. Strandlund, J. Odqvist, and J. A?gren, ?An effective mobilityapproach to solute drag in computer simulations of migrating grainboundaries,? Computational Materials Science, vol. 44, no. 2,pp. 265?273, 2008. ? pages 37, 56150[91] J. Odqvist, B. Sundman, and J. A?gren, ?A general method forcalculating deviation from local equilibrium at phase interfaces,?Acta Materialia, vol. 51, no. 4, pp. 1035?1043, 2003. ? pages 38[92] G. Lame? and B. Clapeyron, ?Me?moire sur la solidification parrefroidissement d?un globe solid,? Annales de Chimie et de Physique,vol. 47, pp. 250?256, 1831. ? pages 40[93] J. Crank, The Mathematics of Diffusion. Oxford University Press,second edition ed., 1990. ? pages 40[94] F. Vasil?ev, Encyclopedia of Mathematics. Springer, 2001. ? pages40[95] J. Rowlinson, ?Translation of J.D. van der Waals? ?Thethermodynamic theory of capillarity under the hypothesis of acontinuous variation of density?,? J. Statist. Phys., vol. 20,pp. 197?245, 1979. ? pages 40[96] J. Gibbs, Collected works. Yale University Press, 1957. ? pages 40[97] R. Becker, ?Die Keimbildung bei der Ausscheidung in metallischenMischkristallen,? Annalen der Physik, vol. 424, no. 1-2, pp. 128?140,1938. ? pages 40[98] E. A. Guggenheim, ?Statistical thermodynamics of the surface of aregular solution,? Transactions of the Faraday Society, vol. 41, no. 0,pp. 150?156, 1945. ? pages[99] O. K. Rice, ?Shape of the Coexistence Curve near the CriticalTemperature,? The Journal of Chemical Physics, vol. 23, no. 1,pp. 164?168, 1955. ? pages 40[100] J. W. Cahn and J. E. Hilliard, ?Free Energy of a NonuniformSystem. I. Interfacial Free Energy,? The Journal of ChemicalPhysics, vol. 28, no. 2, pp. 258?267, 1958. ? pages 40, 46[101] J. W. Cahn, ?Free Energy of a Nonuniform System. II.Thermodynamic Basis,? The Journal of Chemical Physics, vol. 30,no. 5, pp. 1121?1124, 1959. ? pages 40[102] V. Ginzburg and L. Landau, ?On the theory of superconductivity,?Zh. Eksp. Teor. Fiz., vol. 20, p. 1064, 1950. English translation in:L. D. Landau, Collected papers (Oxford: Pergamon Press, 1965) p.54. ? pages 40151[103] S. M. Allen and J. W. Cahn, ?A microscopic theory for antiphaseboundary motion and its application to antiphase domaincoarsening,? Acta Metallurgica, vol. 27, no. 6, pp. 1085?1095, 1979.? pages 41, 45, 49, 71[104] A. Khachaturyan, Theory of Structural Transformations in Solids.New York, Wiley, 1987. ? pages 41[105] P. C. Hohenberg and B. I. Halperin, ?Theory of dynamic criticalphenomena,? Reviews of Modern Physics, vol. 49, no. 3, pp. 435?479,1977. ? pages 41[106] A. A. Wheeler, W. J. Boettinger, and G. B. McFadden, ?Phase-fieldmodel for isothermal phase transitions in binary alloys,? PhysicalReview A, vol. 45, no. 10, p. 7424, 1992. ? pages 41, 56[107] A. Karma, ?Phase-Field Formulation for Quantitative Modeling ofAlloy Solidification,? Physical Review Letters, vol. 87, no. 11,p. 115701, 2001. ? pages 41[108] L. Q. Chen and W. Yang, ?Computer simulation of the domaindynamics of a quenched system with a large number of nonconservedorder parameters: The grain-growth kinetics,? Physical Review B,vol. 50, no. 21, p. 15752, 1994. ? pages 41, 47, 51, 57[109] I. Steinbach, F. Pezzolla, B. Nestler, M. Seeelberg, R. Prieler, G. J.Schmitz, and J. L. L. Rezende, ?A phase field concept for multiphasesystems,? Physica D: Nonlinear Phenomena, vol. 94, no. 3,pp. 135?147, 1996. ? pages 41, 49[110] W. J. Boettinger, J. A. Warren, C. Beckermann, and A. Karma,?Phase-Field Simulation of Solidification,? Annual Review ofMaterials Research, vol. 32, no. 1, pp. 163?194, 2002. ? pages 41[111] N. Provatas, M. Greenwood, B. Athreya, N. Goldenfeld, andJ. Dantzig, ?Multiscale modeling of solidification: Phase-fieldmethods to adaptive mesh refinement,? Internatioanl Jounral ofModern Physics B, vol. 19, no. 31, pp. 4525?4565, 2005. ? pages 41,66[112] J. A. Warren, R. Kobayashi, A. E. Lobkovsky, and W. Craig Carter,?Extending phase field models of solidification to polycrystallinematerials,? Acta materialia, vol. 51, no. 20, pp. 6035?6058, 2003. ?pages 41152[113] M. Gururajan and T. Abinandanan, ?Phase field study of precipitaterafting under a uniaxial stress,? Acta Materialia, vol. 55, no. 15,pp. 5015?5026, 2007. ? pages 41[114] M. Greenwood, J. J. Hoyt, and N. Provatas, ?Competition betweensurface energy and elastic anisotropies in the growth of coherentsolid-state dendrites,? Acta Materialia, vol. 57, no. 9, pp. 2613?2623,2009. ? pages 41, 66[115] A. Artemev, Y. Jin, and A. Khachaturyan, ?Three-dimensionalphase field model of proper martensitic transformation,? ActaMaterialia, vol. 49, no. 7, pp. 1165?1177, 2001. ? pages 41[116] N. Moelans, B. Blanpain, and P. Wollants, ?An introduction tophase-field modeling of microstructure evolution,? Calphad, vol. 32,no. 2, pp. 268?294, 2008. ? pages 41[117] I. Steinbach, ?Phase-field models in materials science,? Modellingand Simulation in Materials Science and Engineering, vol. 17, no. 7,p. 073001, 2009. ? pages[118] H. Emmerich, The Diffuse Interface Approach in Materials Science:Thermodynamic Concepts and Applications of Phase-Field Models.Springer, 2003. ? pages[119] N. Provatas and K. Elder, Phase-Field Methods in Materials Scienceand Engineering. John Wiley & Sons, 2010. ? pages 41[120] J. W. Cahn, ?On spinodal decomposition,? Acta Metallurgica, vol. 9,no. 9, pp. 795?801, 1961. ? pages 44[121] N. Moelans, B. Blanpain, and P. Wollants, ?Quantitative Phase-FieldApproach for Simulating Grain Growth in Anisotropic Systems withArbitrary Inclination and Misorientation Dependence,? PhysicalReview Letters, vol. 101, no. 2, pp. 025502?4, 2008. ? pages 47[122] N. Moelans, B. Blanpain, and P. Wollants, ?Quantitative analysis ofgrain boundary properties in a generalized phase field model forgrain growth in anisotropic systems,? Physical Review B (CondensedMatter and Materials Physics), vol. 78, no. 2, pp. 024113?23, 2008.? pages 48, 49, 52, 64, 105, 107, 112, 158, 159153[123] D. Fan and L. Q. Chen, ?Diffuse-interface description of grainboundary motion,? Philosophical Magazine Letters, vol. 75, no. 4,pp. 187?196, 1997. ? pages 49, 64[124] ?Access e. V. , Aachen, Germany, www.micress.de.? ? pages 50[125] Y. Suwa, Y. Saito, and H. Onodera, ?Phase field simulation of storedenergy driven interface migration at a recrystallization front,?Materials Science and Engineering: A, vol. 457, no. 1-2, pp. 132?138,2007. ? pages 51[126] B. Zhu and M. Militzer, ?3D phase field modelling of recrystallizationin a low-carbon steel,? Modelling and Simulation in MaterialsScience and Engineering, vol. 20, no. 8, p. 085011, 2012. ? pages 51[127] Y. Suwa, Y. Saito, and H. Onodera, ?Phase-field simulation ofrecrystallization based on the unified subgrain growth theory,?Computational Materials Science, vol. 44, no. 2, pp. 286?295, 2008.? pages 51[128] T. Takaki, Y. Hisakuni, T. Hirouchi, A. Yamanaka, and Y. Tomita,?Multi-phase-field simulations for dynamic recrystallization,?Computational Materials Science, vol. 45, no. 4, pp. 881?888, 2009.? pages 52[129] Y. Suwa, Y. Saito, and H. Onodera, ?Three-dimensional phase fieldsimulation of the effect of anisotropy in grain-boundary mobility ongrowth kinetics and morphology of grain structure,? ComputationalMaterials Science, vol. 40, no. 1, pp. 40?50, 2007. ? pages 52, 61[130] I. McKenna, M. Gururajan, and P. Voorhees, ?Phase field modelingof grain growth: effect of boundary thickness, triple junctions,misorientation, and anisotropy,? Journal of Materials Science,vol. 44, no. 9, pp. 2206?2217, 2009. ? pages 52[131] N. Moelans, B. Blanpain, and P. Wollants, ?Pinning effect ofsecond-phase particles on grain growth in polycrystalline filmsstudied by 3-D phase field simulations,? Acta Materialia, vol. 55,no. 6, pp. 2173?2182, 2007. ? pages 53, 92[132] M. F. Ashby and R. M. Centamor, ?Dragging of Small OxideParticles by Migrating Grain Boundaries in Copper,? ActaMetallurgica, vol. 16, no. 9, p. 1081, 1968. ? pages 54, 56154[133] A. Harun, E. A. Holm, M. P. Clode, and M. A. Miodownik, ?Oncomputer simulation methods to model Zener pinning,? ActaMaterialia, vol. 54, no. 12, pp. 3261?3273, 2006. ? pages 55, 89, 91[134] Y. Suwa, Y. Saito, and H. Onodera, ?Phase-field simulation ofabnormal grain growth due to inverse pinning,? Acta Materialia,vol. 55, no. 20, pp. 6881?6894, 2007. ? pages 55[135] L. Vanherpe, N. Moelans, B. Blanpain, and S. Vandewalle, ?Pinningeffect of spheroid second-phase particles on grain growth studied bythree-dimensional phase-field simulations,? Computational MaterialsScience, vol. 49, no. 2, pp. 340?350, 2010. ? pages 55[136] Y. Suwa, ?Phase-field Simulations of Grain Growth inTwo-dimensional Systems Containing Finely Dispersed CoarseningParticles,? ISIJ International, vol. 52, no. 4, pp. 582?591, 2012. ?pages 56[137] S. Vedantam and A. Mallick, ?Phase-field theory of grain growth inthe presence of mobile second-phase particles,? Acta Materialia,vol. 58, no. 1, pp. 272?281, 2010. ? pages 56[138] P. R. Cha, S. G. Kim, D. H. Yeon, and J. K. Yoon, ?A phase fieldmodel for the solute drag on moving grain boundaries,? ActaMaterialia, vol. 50, no. 15, pp. 3817?3829, 2002. ? pages 56[139] A. A. Wheeler, W. J. Boettinger, and G. B. McFadden, ?Phase-fieldmodel of solute trapping during solidification,? Physical Review E,vol. 47, no. 3, p. 1893, 1993. ? pages 56[140] N. A. Ahmad, A. A. Wheeler, W. J. Boettinger, and G. B.McFadden, ?Solute trapping and solute drag in a phase-field modelof rapid solidification,? Physical Review E, vol. 58, no. 3, p. 3436,1998. ? pages 56[141] D. Fan, S. P. Chen, and L. Q. Chen, ?Computer simulation of graingrowth kinetics with solute drag,? J. Mater. Res, vol. 14, no. 3, 1999.? pages 57[142] D. Fan and L. Q. Chen, ?Topological evolution during coupled graingrowth and Ostwald ripening in volume-conserved 2-D two-phasepolycrystals,? Acta Materialia, vol. 45, no. 10, pp. 4145?4154, 1997.? pages 57155[143] D. Fan and L. Q. Chen, ?Computer simulation of grain growth usinga continuum field model,? Acta Materialia, vol. 45, no. 2,pp. 611?622, 1997. ? pages 61[144] P. M. Dupuy, H. Svendsen, L. E. Patruno, and H. A. Jakobsen,?Finite-difference stable stencils based on least-square quadricfitting,? Computer Physics Communications, vol. 181, no. 10,pp. 1671?1681, 2010. ? pages 62[145] N. Provatas, N. Goldenfeld, and J. Dantzig, ?Adaptive MeshRefinement Computation of Solidification Microstructures UsingDynamic Data Structures,? Journal of Computational Physics,vol. 148, no. 1, pp. 265?290, 1999. ? pages 66[146] C. W. Lan, C. M. Hsu, and C. C. Liu, ?Efficient adaptive phase fieldsimulation of dendritic growth in a forced flow at low supercooling,?Journal of Crystal Growth, vol. 241, no. 3, pp. 379?386, 2002. ?pages 66[147] P. Rios and G. Fonseca, ?Grain boundary pinning by Al6Mnprecipitates in an Al?1wt%Mn alloy,? Scripta Materialia, vol. 50,no. 1, pp. 71?75, 2004. ? pages 71[148] K. Song and M. Aindow, ?Grain Boundary Curvature in a ModelNi-Based Superalloy,? Metallurgical and Materials Transactions A,vol. 38, no. 1, pp. 1?6, 2007. ? pages 71[149] D. Kammer and P. Voorhees, ?The morphological evolution ofdendritic microstructures during coarsening,? Acta Materialia,vol. 54, no. 6, pp. 1549?1558, 2006. ? pages 71[150] W. Go`zacutedzacute and R. Holstrokyst, ?High Genus PeriodicGyroid Surfaces of Nonpositive Gaussian Curvature,? PhysicalReview Letters, vol. 76, no. 15, p. 2726, 1996. ? pages 71[151] M. E. Glicksman, ?Analysis of 3-D network structures,?Philosophical Magazine, vol. 85, no. 1, pp. 3?31, 2005. ? pages 83[152] Y. Bre?chet and M. Militzer, ?A note on grain size dependentpinning,? Scripta Materialia, vol. 52, no. 12, pp. 1299?1303, 2005. ?pages 99[153] G. B. Arfken and H. Weber, Mathematical methods for physicists.Amsterdam: Elsevier Academic Press, 6th ed. ed., 2005. ? pages 158156Appendix AThe Interpolation Function for Driving pressureHere, we demonstrate that the driving pressure function, Eq. 8.3 and itsresulting term (Eq. 8.6) for the phase field PDE satisfies the fundamentalsharp interface relationship between interface velocity and the applieddriving pressure, i.e. v = M?G. We consider an interface between twoorder parameters with an energy difference of ?Gij . An unknown termcan be assumed for the function g?(?i) that interpolates ?Gij across theinterface. Eq. 8.7 for two order parameters can be written as,??i?t= ?L[?f0??i? ??2?i + g?(?i)?Gij]. (A.1)For the order parameter ?i, the curvilinear coordinate system is definedwith a vector ni that is normal to the iso-?i profile. In the curvilinearsystem, the Laplacian term can be rearranged using ??i = (??i/?ni)ni and?2?i = ? ???i = (?2?i/?n2i +? ?ni ??i/?ni). Therefore, Eq. A.1 is writtenas:??i?t= ?L[?f0??i? ?(?2?i?n2i+ (? ? ni)??i?ni)+ g?(?i)?Gij](A.2)The weighting function g?(?) is chosen to satisfy two conditions: (i) thatthe application of the driving pressure (?Gij) does not change the shapeof the interface profile, and (ii) that the relationship between the interface157velocity and the driving pressure is linear (v = M?Gij). If the curvatureradius of the interface is much larger than the interface thickness, the curvedinterface has the same order parameter profile in the curvilinear coordinatesystem as a flat interface. For a flat interface ? ? ni = 0 and without adriving pressure (?Gij = 0), it remains stationary i.e. ??i/?t = 0. Thus,the following equality can be obtained from Eq. A.2:??i?t= 0??f0??i= ??2?i?n2i(A.3)Using Eq. A.3 and ??i/?t = ??i/?ni ? ?ni/?t simplifies Eq. A.2 to thefollowing equation:v?i =?ni?t= ?L[??? ? ni + g?(?i)?Gij?ni??i](A.4)where v?i is the velocity field (Eq. 8.9) associated with ?i and is equalto the interface velocity, v, between grain i and j. The gradient of a normalvector, ? ?ni, is equal to the curvature of the interface [153]. Moelans et al.[122] demonstrated that L? = M?. In the sharp interface limit, we have:v = M(?K + ?Gij), (A.5)By comparing Eqs. A.4 and A.5, g? can be obtained as:g?(?i) = ?M/L??i?ni(A.6)If the interface profile does not change when an additional drivingpressure is applied then we can use the profile obtained for a stationaryflat interface. The interface profile can be obtained by solving Eq. A.3 fortwo order parameters [122]:?i =12[1? tanh(?m2?ni)](A.7)The derivative, ??i/?ni, as a function of ?i is:158??i?ni= ?2?m2?(?i ? ?2i). (A.8)And according to Eq. 3.13:L =23?m2?M. (A.9)Consequently, by substituting Eqs. A.9 and A.8 into Eq. A.6 theweighting function can be obtained, i.e.,g?(?i) = 3(?i ? ?2i ) = 3?i(1? ?i). (A.10)During the derivation of g?, it is assumed that the applied drivingpressure does not change the order parameter profile. Therefore, the profileof a moving interface must be a solution of phase field PDE (Eq. 8.7).Transition of the coordinate system to a co-moving frame of reference allowsus to write the following equation for the moving ?i:?i =12[1? tanh(?m2?ni ? vt)](A.11)Rewriting the velocity of the diffuse interface in the curvilinear coordi-nate system (Eq. A.4) results in:v?i = ?L[??K + 3?i(1? ?i)?Gij?ni??i]. (A.12)The sharp interface velocity (Eq. A.5), together with Eqs. 3.12 and A.9,are substituted into Eq. A.12 to obtain the equality:??2m??i(1? ?i)?ni??i= 1. (A.13)The profile of a moving interface given in Eq. A.11 is indeed a solution ofEq. A.13. Therefore, this satisfies the phase field equation and demonstratesthat applying driving pressure with the given g? does not change the shapeof the order parameter profile. g? can be simplified further using ?i = 1? ?j[122]. Consequently, the driving pressure term in the phase field equation159(Eq. 8.6) is written as:?fd??i= g?(?i)?Gij = 3?i?j?Gij . (A.14)The contribution of the driving pressure to the energy functional, fd(Eq. 8.3), is obtained by the integration of Eq. A.14 with the condition thatfd(?i = 1) = Gi.160
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Phase field modelling of grain growth with particle...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Phase field modelling of grain growth with particle pinning and solute drag Shahandeh, Sina 2013
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Phase field modelling of grain growth with particle pinning and solute drag |
Creator |
Shahandeh, Sina |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | Second phase particles and solute atoms have been used as an important constituent in the design of materials and processes due to their ability to restrain the motion of grain boundaries. The drag effect occurs on a scale comparable to the particle diameter and interface thickness. However, to simulate grain growth with numerical efficiency one requires a model that captures the drag pressure on the interfaces but does not resolve the particles or solute segregation spike. In this work, a multi-scale modelling scheme is proposed to simulate grain growth with particle pinning and solute drag. The interaction of a grain boundary with an ensemble of particles is simulated to obtain the pinning pressure. A phase field model is then developed that incorporates the drag pressure in the meso-scale and simulates grain growth. The accuracy of the model is confirmed in comparison with analytical expressions. The application of the model is then presented for grain growth in two- and three-dimensional systems under the influence of particle pinning. Measuring the curvature of the grain boundary network reveals that in the completely pinned structure, the average driving pressure is not equal to but lower than the pinning pressure. The results of the nano-scale simulations for pinning pressure is combined with the results from the meso-scale to produce a limiting grain size that coincides with the experiments. This curvature analysis provides a kinetic model that describes the evolution of the structure more accurately than that of the mean field theories. The proposed phase field formulation is also applied to simulate grain growth in the presence of solute drag. The grain growth kinetics follows a phenomenological relationship that is described using a power law, with a time exponent in the range of 0.35 to 0.50. The deviation from ideal grain growth, associated with a time exponent lower than 0.50, and its correlation with the solute drag parameters is investigated. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-09-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-ShareAlike 2.5 Canada |
DOI | 10.14288/1.0103348 |
URI | http://hdl.handle.net/2429/45137 |
Degree |
Doctor of Philosophy - PhD |
Program |
Materials Engineering |
Affiliation |
Applied Science, Faculty of Materials Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-sa/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2013_fall_shahandeh_sina.pdf [ 15.66MB ]
- Metadata
- JSON: 24-1.0103348.json
- JSON-LD: 24-1.0103348-ld.json
- RDF/XML (Pretty): 24-1.0103348-rdf.xml
- RDF/JSON: 24-1.0103348-rdf.json
- Turtle: 24-1.0103348-turtle.txt
- N-Triples: 24-1.0103348-rdf-ntriples.txt
- Original Record: 24-1.0103348-source.json
- Full Text
- 24-1.0103348-fulltext.txt
- Citation
- 24-1.0103348.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0103348/manifest