Diffusional evolution of precipitates in elastic media

Diffusional evolution of precipitates in elastic media using theextended finite element and the level set methodsRavindra Duddu a ,⇑,David L.Chopp b ,Peter Voorhees c ,Brian Moran a ,1a Department of Civil and Environmental Engineering,Northwestern University,2145Sheridan Rd,Evanston,IL 60208,USAbDepartment of Engineering Sciences and Applied Mathematics,Northwestern University,2145Sheridan Rd,Evanston,IL 60208,USAc Department of Material Science and Engineering,Northwestern University,2145Sheridan Rd,Evanston,IL 60208,USA a r t i c le i nf o Article history:Received 9January 2010Received in revised form 22September 2010Accepted 1November 2010Available online 12November 2010Keywords:Phase transformation XFEM Level sets Dendritic evolution Anisotropy Grid rotationa b s t r a c tA sharp-interface numerical formulation using an Eulerian description aimed at modelingdiffusional evolution of precipitates produced by phase transformations in elastic media,ispresented.The extended finite element method (XFEM)is used to solve the field equationsand the level set method is used to evolve the precipitate–matrix interface.This new for-mulation is capable of handling microstructures with arbitrarily shaped particles and cap-turing their topological transitions without needing the mesh to conform with theprecipitate–matrix interface.The XFEM makes it possible to model the precipitate andthe matrix to be both elastically anisotropic and inhomogeneous with ease.The interfaceevolution velocity is evaluated using a domain integral scheme [1]that is consistent withthe sharp interface.Numerical examples modeling two distinct phases of particle evolu-tion,growth (dendritic evolution)and equilibration (Ostwald ripening)are presented.Toovercome the issue of grid anisotropy in growth simulations,a random grid rotationscheme is implemented in conjunction with a bicubic spline interpolation scheme.Grow-ing shapes are dendritic while equilibrium shapes are squarish and in this respect our sim-ulation results are in agreement with those presented in the literature [2–4].Ó2010Elsevier Inc.All rights reserved.1.Introduction‘‘Whenever a single phase system is put into a two-phase metastable state,a second phase nucleates,grows and coars-ens’’[5].The growth of second phase domains from a supersaturated matrix,and coarsening or Ostwald ripening of two phase binary alloys are examples of diffusional phase transformation in crystalline solids.Numerical modeling of micro-structure evolution in a binary alloy involves tracking the morphology of arbitrarily shaped precipitates in an elastically stressed matrix.Particle evolution has two distinct stages:(i)growth of the particle that satisfies a local mass balance rela-tion at each precipitate interface;(ii)coarsening of the particle (equilibration)that involves dynamic rearrangement of mass in the system so as to minimize the system energy.The particle evolution problem (growth and equilibration)falls into the class of free boundary problems and the current work is our first step towards developing a general numerical formulation to simulate microstructural evolution.The precipitate–matrix interface evolves with time and separates two regions where the governing equations differ significantly,leading to a solution with discontinuous derivatives through the interface.The computational challenges of 0021-9991/$-see front matter Ó2010Elsevier Inc.All rights reserved.doi:10.1016/j.jcp.2010.11.002⇑Corresponding author.Address:Department of Civil Engineering and Engineering Mechanics,Columbia University,610Seeley W.Mudd Building,500West 120th Street,Mail Code 4709,New York,NY 10027,USA.Tel.:+15129340967.E-mail addresses:rduddu@ ,r-duddu@ (R.Duddu),chopp@ (D.L.Chopp).1Present address:Division of Physical and Environmental Sciences and Engineering,King Abdullah University of Science and Technology,Thuwal,Saudi Arabia.1250R.Duddu et al./Journal of Computational Physics230(2011)1249–1264precipitate evolution models are significant due to the presence of:discontinuous strains across the interface;cubic elastic anisotropy;elastic inhomogeneity;discontinuous substrate gradients across the interface;moving interfaces;inter-particle interactions.Also,computing the interfacial curvature and the evolution velocity present a significant challenge.Therefore, specialized methods for handling moving boundaries and for solving diffusion and elasticity equations on irregularly shaped domains are needed for producing accurate simulations.The different approaches for simulating microstructure evolution can be broadly classified on the basis of how the interface is resolved(sharp vs diffused interface)and how it is represented (Lagrangian vs Eulerian).Previous work on simulating the evolution of two phase microstructures can be classified as:sharp interface methods in Lagrangian framework[2–4,6–11];diffused(or smooth)interface methods in Eulerian framework[12–18];sharp interface methods in Eulerian framework[19–22].The smooth interface approach is suitable for simulating microstructural evolution with large numbers of particles and for capturing topological transitions,such as particle merging,splitting,and vanishing.The phasefield method[23–27]falls under the class of smooth or diffused interface methods.The drawback of the phasefield method is the computational com-plexity involved with solving the evolution equations with a large number of parameters[19].The sharp interface approach based on the boundary integral formulation presented in[2–4,6,7,11]resolves the interface byfitting cubic splines through marker particles.A major advantage of the boundary integral formulation(BIF)is that it reduces the dimensionality of the problem.We need N marker particles only on the interface whereas smooth interface methods require an NÂN grid.How-ever,handling topological transitions with the interface tracking scheme is not straightforward.Also in practice,the numer-ical implementation of the BIF is made difficult by the complicated structure of the equations for the elasticfields and the presence of high order time step stability constraints(stiffness)that arise when interfaces with surface tension are tracked [2].This restricted the early work to systems with fewer particles[2,4,6].Later[28],the boundary integral formulation was extended to handle a sufficiently large number of particles.Generally,to reduce the complexity of the governing equations, the BIF generally makes certain assumptions of:quasi-steady state diffusion[2–4,6,7,11];elastic homogeneity and/or isot-ropy[2–4,6].The implementation of both elastic inhomogeneity and anisotropy using the BIF is more difficult,although it is noted in[2]that it should be possible in their general framework.To our knowledge,there are no existing numerical formu-lations that are general enough to model the evolution of particles considering non steady-state diffusion in an elastically inhomogeneous and anisotropic media.The proposed Eulerian formulation resolves the interface sharply and is also able to handle large numbers of particles and topological transitions without difficulty.The current numerical scheme utilizes the XFEM[29–31]to solve thefield equations and the level set method[32–34]to evolve the precipitate–matrix interface.The level set method has been pre-viously used to simulate epitaxial phenomena[22]and microstructural evolution during thin-film deposition[21].A major advantage of the proposed numerical method is that the precipitate–matrix interface and its evolution can be captured without needing the mesh to conform with the interface.A similar computational method has been used to study solidifi-cation problems,where the temperaturefield evolves according to classical time dependent heat conduction[19,20].Here, we study the effects of elastic stress and interfacial curvature on the morphology of the precipitate at a constant temper-ature.Although we solve the steady state substrate diffusion equation,the current numerical technique can be easily ex-tended to solve the time dependent diffusion equation[19],while considering an elastically inhomogeneous and anisotropic medium.Thus,the current numerical formulation can be extremely useful to study the complex physics behind evolving microstructures.An important numerical issue while simulating dendritic evolution of microstructures using mesh/grid based schemes is that of grid anisotropy.The evolution of a spherical particle whose radius is greater than a certain critical radius,undergoing diffusion-controlled growth in a supersaturated binary matrix,is known to be unstable[35].The evolution of the precipitate is affected by the anisotropy or noise arising from the regular computational grid.Although this anisotropy is very small,it might give rise to unphysical shapes as the numerical errors accumulate with time.Early lattice models presented in[36–38] have the problem of ter,concepts like random self-adaptive grid of anchor points for the diffusionfield[39] and grid generation methods popular infinite element calculations[40]were introduced to tackle the issue of grid anisot-ropy,but they turned out to be computationally expensive.A rotated lattice sandwich construction was suggested in[41]to reduce the anisotropy.The basic idea of this method is to implement the simulations on a stack of four lattices that are ro-tated against each other by random irrational angles.The diffusionfield is computed and the interface is locally advanced on each of these lattices at each time step.The new interface location is taken to be the average of the interfaces over the four lattices.Although this increases the computational cost by a factor of4,the method is simple and has produced good results.A simple yet effective idea for reducing the undesirable grid anisotropy is to rotate the grid by a random angle after every few time steps.Such a random grid rotation method using a twelve-point interpolation scheme was used in[42]to remove the anisotropy.A rotation scheme involving a45°rotation in a random direction(clockwise or anticlockwise)with respect to the original orientation was employed in[43].Another method of eliminating the grid anisotropy by quantitatively determining the numerical noise contributions was suggested in[44].However,determining the contribution of the grid anisotropy by matching phase-field calculations of the equilibrium shape(that supposedly could have some grid anisotropy)with analyt-ical Gibbs-Thompson calculations of that same shape(that does not include grid anisotropy)can be tedious.The random grid rotation scheme proposed in this article is somewhat similar to that described in[42].The remainder of this article is orga-nized as follows:Section2briefly presents the governing equations of the continuum model;Section3gives the detailed algorithm;Section4describes briefly the numerical methods employed and the implementation;Section5presents2D numerical results of equilibrium morphology and dendritic growth.Some concluding remarks are made in Section6.2.Formulation2.1.Model assumptionsIn this section we present the key assumptions which help simplify the two-dimensional precipitate evolution model pre-viously presented in [4,6].1.The precipitate-matrix interfaces are assumed to be coherent.The shape of the particles (Ni 3Al)is evolved in a mannerconsistent with the flow of mass in the system.2.The misfit or the transformation strain between particle and matrix is taken to be purely dilatational.3.The precipitate is modeled as an elastic solid particle in an elastic binary alloy matrix.This two-phase binary system(matrix and precipitate)is assumed to be elastically homogenous and anisotropic with cubic symmetry for benchmark calculations.The lattice parameters of both phases and the elastic constants are assumed to be independent of composition.4.The system is considered to be in a state of 2D plane strain and so there are no strains out of the plane of the particle.5.No external traction is applied to the system at infinity.Because the computational domain cannot be infinitely large,weconsider a domain sufficiently large to impose traction-free boundary conditions on the external boundary.We assume zero mass flux through the external domain boundary for equilibration phase of evolution.Hence,the total volume of the particle phase is a constant during the equilibration process.6.Precipitate growth is driven by the establishment of concentration gradients of substrate.Here,the substrate is Al and thediffusion of Al atoms in the Ni–Al matrix establishes the concentration gradient around the particle interface.7.The processes of diffusion of Al atoms and precipitation (which involves attachment of Al atoms at particle–matrix inter-face)occur at different time scales.Typically,the Ni 3Al particles are assumed to evolve via the isothermal diffusion of mass at a sufficiently slow rate,so that the quasi-stationary approximation for the diffusion field in the matrix holds[4,6].In the current article even though we use the quasi-stationary approximation to reproduce the results in literature[4,6],our numerical technique can be used to study the more general unsteady diffusion problem.8.The interfacial energy is assumed to be isotropic and independent of the interfacial strain.2.2.The free-boundary problemWe will now review the domains and the different boundaries for the two-dimensional free boundary problem.Let us consider the evolution of coherent particles (b )in an infinite matrix (a ).Let X be a regular region bounded by the curve @X =C .The domains containing the matrix and the precipitate phases are denoted by X a and X b respectively.The precip-itate domain is bounded by the precipitate–matrix interface ,C int (t ),at all times.The precipitate domain,X b (t ),is a subset of X and X =X b (t )S X a (t )at all times (see Fig.1).We use a Cartesian coordinate system and the spatial (Eulerian)coordinates are denoted by x .Let U (x ,t )be the displacement vector due to the misfit strain–induced deformation of the domain at time,t .Let T (x ,t )denote the stress tensor and E (x ,t )denote the total strain tensor measured with respect to the stress-free a .The stress field,T ,in the two phase binary system is determined using the equilibrium equations,the strain–displacement rela-tions and the constitutive equations (stress–strain relations).The equations of equilibrium in both phases (matrix and pre-cipitate)are,r ÁT ¼0:ð1ÞThe above equations are solved in the domain with boundaryconditions,R.Duddu et al./Journal of Computational Physics 230(2011)1249–12641251T Án ¼0on C ;s T t n ¼0on C int ;s U t ¼0on C int ;ð2Þwhere n is the unit outward normal to the interface.The terms s T t 2and s U t represent the jump in traction and displacements respectively across the interface.The first equation in (2)is the traction-free boundary condition imposed on the domain bound-ary.The last two equations in (2)impose the continuity of tractions and displacements at the interface.The total strain tensor,E ,measured with respect to the stress-free a ,can be written in terms of the displacement vector,U as:E ¼12ðr U þðr U ÞT Þ;ð3Þwhere superscript T denotes the transposition and r is the spatial gradient operator.Let us denote the magnitude of the purely dilatational misfit strain by and the elastic strain tensor by e E.The constitutive law for both the matrix and the pre-cipitate is given by,T ¼K e E;ð4Þwhere K is the tensor of elastic moduli.The elastic strains in the particle and matrix are defined by,e Eb ¼E b À 1in X b ;e E a ¼E a in X a ;)ð5Þwhere 1is the second-order unit tensor.Eqs.(1)–(5)define the two-dimensional plane strain elasticity problem to solve for the stress state in the phases.Let C be the mole fraction of component 2of the binary alloy (mole fraction of Al in the case of Ni–Al).The scalar substrate concentration field,C (x ,t ),in the matrix phase at steady state is given by a solution tor 2C a ¼0in X a :ð6ÞThe above equation is subject to the interfacial conditions given by the stress modified Gibbs–Thomson equation.By assum-ing local equilibrium,or negligible interfacial kinetics,the interfacial compositions on the matrix side of the interface are given by [28],C a ¼C a i þB a 12s T Áe E t ÀT a Ás E t þr K on C aint ;ð7Þwhere C a is the interfacial concentration in the matrix (a )phase,C ai is the equilibrium concentration in the matrix at a flatinterface in the incoherent or stress-free state,r is the isotropic interfacial energy,K is the interfacial curvature reckoned positive for a spherical b particle,B is a constant in the designated phase and is proportional to the second derivative of the molar free energy.The interfacial composition on the particle (Ni 3Al)side of the interface is equal to the composition in the particle domain,which is a constant given by,C b ¼mole fraction of Al mole fraction of Al þmole fraction of Ni ¼11þ3¼0:25:in X b ;C bi ¼0:25on C bint ;ð8Þwhere C b is the concentration in the precipitate (b )phase,C bi is the concentration in the matrix at the interface.It is apparentthat there is a jump in the value of substrate concentration,s C t =(C b ÀC a ),at the interface.The final boundary condition is the concentration at infinity,which for all computational purposes,is the domain boundary C .For equilibrium shape evo-lution of particles,we fix the concentration at infinity by requiring that the total volume of the particle phase remains con-stant (flux-free boundary),XN p i Z A i r C a Án dA ¼0;ð9Þwhere A i is the area of the i th particle–matrix interface and N p is the number of particles considered.The above equation states that for a multiple particle system,the sum of the fluxes of substrate at the particle–matrix interfaces should be equal to zero when there is no flux at exterior boundary of the domain.Let V be the normal interface speed (a scalar)with which the interface,C int ,evolves in time.From the mass balance of substrate at the interface we get the following expression for V (x C ,t )at a point x C on the interface,s C t V ðx C ;t Þ¼D r C j a Án ðx C ;t Þ;ð10Þ2The term s w t =w b Àw a denotes the jump in the value of w at the interface.1252R.Duddu et al./Journal of Computational Physics 230(2011)1249–1264where r C j a is the substrate concentration gradient evaluated on the matrix side of the interface at x C.Eqs.(6)–(10)specify an exterior diffusion problem with Dirichlet boundary conditions on the interface in a multiply-connected domain.All the Eqs.(1)–(10),put together define the two-dimensional microstructural evolution problem,which is a coupled system of partial differential equations in the unknowns,U and C.3.Solution strategyThe assumption of linear elasticity and the quasi-stationary approximation of the diffusionfield(see Section2.1)elimi-nate the need to track the stress and substrate concentrationfield history.This allows us to uncouple the system of partial differential equations that define the evolution problem.Thus,we can solve for the stressfield,T,the concentrationfield,C, given just the level set function,/,and then update/using the computed concentrationfield,C,in a sequential manner. Also,the above assumptions eliminate the need for rotation and interpolation of thefield variables for dendritic growth problems when using the random grid rotation scheme.However,if the assumptions are not made then we need to rotate and interpolate thefield variables when the grid is rotated in order to keep track of the history.The proposed numerical scheme is explicit in time and the algorithm used for solving the precipitate growth problem at time,t,is given below: pute the elastic stress tensor,T,in the particle and matrix due to deformation by solving the Eqs.(1)–(5)using XFEMdiscretization;2.Establish the interfacial concentrations,C a j Cint ,using the stress-modified Gibbs-Thompson Eq.(7);3.Solve the Laplace equation for diffusion(6)and determine the concentrationfield,C a in the matrix(a)phase using XFEMdiscretization;pute the normal interface speed,V,by evaluating the concentration gradients at the interface using the domain inte-gral scheme given in[1];5.Evolve the interface by an explicit time update of the level set function,/(11);6.Rotate the grid by a random angle and interpolate/onto a new regular unrotated grid using bicubic spline interpolation.Note that this step is required only for dendritic evolution problems.4.Numerical implementationTo simulate the evolution of microstructures we extend the numerical method combining the extendedfinite element method(XFEM)with the level set method presented in[45].This combined method eliminates the need for the mesh to coincide with the interface and subsequent remeshing as the interface evolves.The convergence of this combined numerical scheme was illustrated in[45]and its applicability to laminarflow driven biofilm growth is presented in[46].In this section we shall briefly describe the details of the numerical methods and their implementation relevant to the problem at hand.For more details on the numerical implementation we encourage the reader to look into the references provided herein.4.1.The level set methodThe level set and fast marching methods werefirst introduced in[47,48].The level set method is based upon representing an interface as the zero level set of some higher dimensional function,/.Here,the particle–matrix interface is maintained by the level set function,/(x,t),such that/<0in the precipitate(b)domain,and/>0in the matrix(a)domain.The interface is updated by solving the level set evolution equation,@/@tþV ext k r/k¼0;ð11Þwhere V ext is obtained from fast construction of extension velocities.The extended speed function,V ext,is evaluated using the normal interface speed function,V(for more details see[45]).The speed function,V,given by(10)is evaluated using the domain integral method presented in[1].The level set evolution Eq.(11)is solved in the whole computational domain using an upwindfinite difference scheme and fast marching methods as described in[32–34]forfinite difference based level set methods.4.2.The extendedfinite element approximationThe XFEM is used to compute the elastic stressfield,T,and to establish the substrate concentrationfield,C,in the domain. The strainfield,E,is discontinuous across the precipitate–matrix interface due to the misfit strain arising from phase trans-formation(see Eq.(5)).Therefore,the displacementfield,U,is continuous across the interface but its derivatives are discon-tinuous across the interface.Also,the substrate concentrationfield,C,is discontinuous across the interface(see Eqs.(7)and (8)).The‘‘step enrichment’’of the standardfinite element approximation permits the introduction of the discontinuity in the value of the unknownfields at the interface and also eliminates the need for thefinite element mesh to coincide with the precipitate–matrix interface.R.Duddu et al./Journal of Computational Physics230(2011)1249–12641253In the XFEM,the physical domain,X,is divided into elements where an unknownfield,generically named u in the fol-lowing equation,is given by,uðx;tÞ¼X nI¼1N IðxÞu IðtÞþX n eðtÞJ¼1S Jðx;tÞa JðtÞ;ð12Þwhere N I(x)are the standardfinite element shape functions,u I are nodal degrees of freedom,n is the total number of nodes in the mesh,S denotes the step enrichment function,a J are the nodal enrichment degrees of freedom(which are additional degrees of freedom at the enriched nodes),and n e(t)is the total number of enriched nodes in the mesh at time step t. The step enrichment function,S,described below is used to enrich thefinite element approximations of both the unknown fields U and C.4.3.Step enrichment functionA piecewise continuous C0enrichment function,S,can yield a continuous or discontinuous solution across the interface but requires Lagrange multipliers to apply the Dirichlet jump condition.This enrichment function is defined by: S Jðx;tÞ¼N JðxÞðHð/ðx;tÞÞÀHðÀ/ðx;tÞÞÞ;ð13Þwhere H(/)is the Heaviside step function given by,Hð/Þ¼1;/>0;0;/<0:ð14ÞThis enrichment does not enforce the value of the jump in the solution or the normal derivative across the interface,so two separate conditions are required to maintain well-posedness.Therefore,Lagrange multipliers are used to enforce the jump condition across the interface in the solution and the normal derivative condition across the interface is weakly enforced.The numerical scheme proposed herein uses Lagrange multipliers to impose zero jump in displacements(continuity of displace-ments)and to impose the values of C a and C b(Dirichlet boundary conditions on substrate)at the interface.A more detailed description of the implementation can be found in[49,50].4.4.Curvature computationEstablishing the interfacial concentrations,C a j Cint ,using(7)requires the explicit computation of the interfacial curvature,K,which can be evaluated using the level set function,/[32].Usingfinite differences we can approximate the local inter-facial curvature,K,using the formula given in[32,51].The expression for the local curvature at the grid points in terms of/is given by,K¼rÁr/k r/k¼@@x;@@yÁ/x/2xþ/2y1=2;/y/2xþ/2y1=2B@1C A¼/xx/2yþ/yy/2xÀ2/xy/x/yð/2xþ/2yÞ3=2;ð15Þwhere/x and/y are derivatives of/along x and y directions evaluated using central differences.Thus,computing the inter-facial curvature involves the evaluation of the second-order derivatives of/.For two phaseflows a curvature term appears in the surface tension integral that defines the Neumann jump conditions at the interface.In that case a weak form implemen-tation of the surface tension integral using the Laplace–Beltrami operator can be advantageous because it uses only thefirst-order derivatives of/[52].However,the current problem requires us to compute the local curvature explicitly in order to set the Dirichlet boundary conditions for the substrate concentration at the precipitate–matrix interface.Therefore,we have implemented the bicubic interpolation scheme given in[53]to evaluate the curvature.The curvature at a grid point is eval-uated usingfinite differences and the local curvature at a point on the interface is evaluated using element-wise bicubic interpolation.The bicubic interpolation method uses more data(4Â4grid=16grid points)to construct the interpolating polynomial and ensures the consistency of thefirst derivatives across voxel boundaries.The details of bicubic interpolation and the algorithm can be found in[53,54].4.5.Method of random grid rotationSince our combined numerical scheme based on the XFEM and the level set method is implemented on a structured mesh (that is generated only once),numerical anisotropy can accrue over time leading to unphysical shapes.Even though an unstructured mesh should not have a particular direction of anisotropy,the numerical noise could give rise to different mor-phology for different mesh instances.Also,the implementation of the level set method on an unstructured grid is very com-plicated and therefore,unattractive.The numerical anisotropy of structured meshes definitely affects the unstable evolution of dendrites in a supersaturated solution.However,the stable equilibration of particle shapes is unaffected by numerical 1254R.Duddu et al./Journal of Computational Physics230(2011)1249–1264。