Phase-field Simulation of Dendritic Solidification

Phase-field Simulation of Dendritic Solidification

Christer Andersson

Stockholm 2002

Doctoral Dissertation

Royal Institute of Technology

Akademisk avhandling som med tillstånd av Kungl Tekniska Högskolan framlägges till offentlig granskning för

avläggande av teknologie doktorsexamen fredagen den 30 augusti 2002 kl 10.00 i seminarierum D3, Kungl

Tekniska Högskolan, Lindstedtsvägen 3, Stockholm.

ISBN 91-7283-343-2

TRITA-NA-0219

ISSN 0348-2952

ISRN KTH/NA/R--02/19

© Christer Andersson, juli 2002

Weconsiderthesolidicationofapurematerialfromitsundercooledmelt, tak-ing into accounteects of surfacetension, kineticundercoolingand crystalline anisotropy. The classicalmathematical formulationconsists of two heat equa-tions coupled with boundary conditions prescribed at the moving solid-liquid interface. Analyticsolutionsare known onlyfor afew special cases. Sincethe problem has important industrial applications one is interested in numerical methodsforaccuratesimulation.

Theclassicalformulationisnotwellsuitedforcomputationsinceitrequires explicittrackingofgeometricallycomplexfreeboundaries. Phase-eldmethods haveemergedasapopularalternative. Theyarebasedonthenotionthatthere existsasmoothapproximationofanindicatorfunction,orphase-eld,uniquely identifying the phaseat a given locationand time. Phase-eldmodelsare de-rivedfrom thermodynamicsand are formulated asreaction-diusionequations validin theentire domain. Freeboundariesarereplacedbydiuse interfacesof widthcorrespondingtotransitionsinthephase-eldvariable. Thereisnoneed to explicitly track interfacessincethey areimplicitly dened by thecomputed solution.

Themain purpose ofthis work isto nd computationallyeÆcientmethods forphase-eldsimulationofdendriticsolidication. Tothisendnitedierence approximations and their implementations are considered. A comparison be-tweenarstorderaccurate semi-explicittime-steppingscheme andanimplicit schemebasedonBDF-discretizationispresented. Numericalexperimentsintwo spacedimensionsshowthatforsuÆcientlylowerrortolerancestheBDF-scheme givesasolutioninlesssimulationtime.

Tofurtherreducesimulationtimesfortheimplicitschemeanimplementation has been written for parallel distributed memory architectures. Performance measurementsonanIBMSP2showthataparalleleÆciencyexceeding70%can beobtainedforgridsassmallas200200pointson16processors.

Agreement between phase-eld model and classical formulation is investi-gated in thelimit ofvanishing diuse interfacethickness. It isshownhowthe accuracyin theinterfaceconditionscanbeincreasedfrom rstto secondorder in bymodifyingasinglescalarparameter. Simulationcanbeperformedwith larger, andhencesmallercomputationalgrids,withoutlossofaccuracy.

Analternativelimitwheretheundercoolingtends tozeromuch fasterthan the diuse interface thickness is considered. The analysis shows that the in-terfacialtemperaturebecomesindependentof to leadingorderonlywhenthe modication is applied. At low undercoolingsthe smallest geometric scale in the dendrite increases, and solidication proceeds slower. The correction pre-dicted by analysis allows simulation for smaller undercoolings with the same computationalresources.

Firstand foremost I would liketo thank myadvisor Gunilla Kreiss for invalu-ablesupport, guidanceandencouragementthroughoutmyresearch. Thisthesis wouldnothavebeenpossiblewithouther.

Iamalsogratefulto GustavAmberg, IrinaLoginova,Robert Tonhardtand John



Agrenformanystimulatingandinformativediscussionsonphase-eld mod-elsandtheunderlyingphysicalprocessesofsolidication.

Last,butnotleast,Iwishtothankmyfamilyandfriendsforputtingupwith meoverthepastcoupleofyears.

FinancialsupportfromtheNationalGraduateSchoolinScienticComputing (NGSSC),theSwedishResearchCouncil(TFR)andtheEuropeanCommunity's TMR Programmeisgratefully acknowledged.

1 Introduction 1

2 Solidication of a pure material 3

2.1 TheStefanproblem . . . 3

2.2 Relatedphysicalproblems . . . 5

2.3 Directnumericalsimulation . . . 5

3 The phase-eldapproach 7 3.1 Derivationofaphase-eldmodel . . . 7

3.2 Non-dimensionalization . . . 10

3.3 Model properties . . . 11

3.4 Inclusionof anisotropy . . . 12

4 Numericalsimulation 15 4.1 Finitedierencemethods . . . 15

4.2 AsymptoticsandimplicationsoncomputationaleÆciency . . . . 16

4.3 Parallelcomputing . . . 17

4.3.1 Scalabilitymetrics . . . 18

5 Computeddendrites 21 6 Summaryof papers 25 6.1 Paper1: Time-steppingschemesforphase-eldsimulationof den-driticsolidication . . . 25

6.2 Paper2: Computationof dendritesonparallel distributed mem-oryarchitectures . . . 26

6.3 Paper3: Thirdorderasymptoticsofaphase-eldmodel . . . 26

6.4 Paper4: Analysisof aphase-eld model forsolidication in the limitofvanishingundercooling . . . 27

Paper1 C. Andersson, Time-stepping schemes for phase-eld simu-lation of dendritic solidication, Tech. Rep. TRITA-NA-0216, Department of Numerical Analysis and Computer Science, RoyalInstituteofTechnology,2002.

Paper2 C.Andersson,Computationofdendritesonparalleldistributed memory architectures, in Simulation and Visualization on the Grid, B. Engquist, L. Johnsson, M. Hammill, and F. Short, eds., vol. 13 of Lecture Notes in Computational Science and Engineering,SpringerVerlag,June2000,pp.195{208.

Paper3 C. Andersson, Third order asymptotics of a phase-eld model, Tech. Rep. TRITA-NA-0217,Department ofNumerical Analysisand ComputerScience, RoyalInstituteof Technology, 2002.

Paper4 C.Andersson,Analysisofaphase-eldmodelforsolidication in the limit of vanishing undercooling, Tech. Rep. TRIT A-NA-0218,DepartmentofNumericalAnalysisandComputerScience, RoyalInstituteofTechnology,2002.

Introduction

Solidication anddendritic growth isveryimportant from apracticalpoint of view. Some properties of solids, e.g. ductility, electrical conductivity and me-chanicalstrength,aredeterminedbythemicroscopicstructuresproducedupon solidication. Onewouldliketogaincontrolofthestructureformationtoobtain thedesiredpropertiesinmanufacture. Applicationscanbefoundinforexample castingandsemiconductorproduction.

Thepattern formationprocessisnotyetfullyunderstood. It isforexample known that the dendrite tip radius multiplied by the tip velocity is constant, butitisnotclearhowthisconstantisselectedforagivenundercooling. Froma scienticpointofviewonewouldliketounderstandtheunderlyingmechanisms inorder tocreatebettermodels.

Hereweespeciallyconsiderthesolidicationofapurematerialfromits under-cooledmelt,taking intoaccounteects ofsurfacetension,kineticundercooling andcrystallineanisotropy. Theclassicalmathematicalformulationhastheform ofaStefanproblem. Twodiusionequationsgovernthetransportofheatinthe solid and liquidphases, and are coupledby boundary conditionsprescribed at thesolid-liquidinterface. Duringsolidicationtheinterfaceswillmove,resulting inafreeboundaryproblem.

In most situtations it is impossible to nd analytic solutions to the Stefan problemofpracticaluse. Forthisreasonaconsiderableeorthasbeenputinto developingnumericalalgorithmsoverthelasttwentyyears. Inparticular, phase-eldmethodshavebecomeincreasinglypopularsincetheyallowsimulation with-out explicit treatment of free boundaries. Chapter 2 givesa brief overview of theclassicalmathematical descriptionanddiscussessomeof thediÆculties as-sociatedwithdirectnumericalsimulation. Thederivationofaphase-eldmodel fromthermodynamicsisoutlinedinChapter3. Thederivationistolargeextent basedontheworkofWanget al. [40].

The main purpose ofthis work isto nd computationallyeÆcientmethods fornumericalsimulationofdendriticsolidication.Itcanberoughlydividedinto

two parts. The rstparttreats thenumericalalgorithms and their implemen-tation on parallel distributed memory computers. Section 4.1givessomebrief preliminaries onapproximationby nitedierencemethods. Paper1compares twodierentapproachestotime-steppingwithrespecttoachievedaccuracyand total executiontime. All aspects ofthenumericalschemesarecoveredin some detail,withtheemphasisonlinearandnon-linearsolvers.

Themoreaccurateofthetwotime-steppingschemeisexpensivein termsof oating point operations and also requires morestorage. Paper2 describes a parallel implementationaimedatreducingthesimulationtimeand/orallowing largerproblemsto besolved. TheeÆciencyoftheparallel codeismeasuredby themetricsintroducedin Section4.3.

Thesecond partof thethesis dealswiththe analysis ofphase-eld models. Phase-eldmodelsarederivedfromthermodynamics,anditisnecessarytoshow thattheresultingequationsyieldsolutionssatisfyingthefreeboundaryproblem onewishestosolve. KarmaandRappelnotedthatamoreaccurateanalysiscan beexploitedtoreducethesizesofcomputationalgridswithoutlosingagreement withthefreeboundaryformulation[21]. Themodicationrequiresonlyasmall correctionto asinglephysicalparameter,asoutlinedin Section4.2.

Inpaper3asymptoticapproximationstosolutionsofaphase-eldmodelare constructed. It is shown that the correction suggested by Karma and Rappel [21] increases accuracy in the interface conditions from rst to second order. This propertydoesnotfollowautomatically from moreaccurateanalysis. The analysisshowsthatthefthdegreepolynomialmodelingthereleaseoflatentheat originally proposedin[40]can neverleadtothirdorderaccuracy. Furthermore, iftherearepolynomialsinthephase-eldmodelleadingtothirdorderaccuracy theyareofdegreetenorhigher. Polynomialsofsuchhighdegreearenotdesirable forcomputation beacuseofsensitivitytoround-oerrors.

Paper 4 presents an alternative asymptotic analysis where the limits are chosentobemoreappropriatefornumericalsimulationatlowundercoolings. It isshownthatthecorrectiondeterminedinpaper3increasestheorderofaccuracy byalmostoneinthislimit. Firstorderconvergenceisobtainedinsteadof\almost zero". Thelimititself isnotofpracticalinterestbutallowsinvestigationofthe phase-eldmodelatlowundercoolings.

Solidication of a pure

material

2.1 The Stefan problem

Let us consider the simple exampleof asingle growing crystal. We can think of the computational domain asconsisting of the regions (t) and

+ (t) correspondingto thesolidand liquidphasesatagiventime t. Thetwophases areseparatedbyafreeboundary (t)representingthesolid-liquidinterface,see Figure2.1. isheredenedastheanglebetweentheoutwardnormalandsome referencedirection. Itisneededtoaccountforpreferredgrowthdirectionsinthe crystal. Inthreedimensionsshouldbeavectoroftwocomponentstoproperly accountforanisotropy.

Givenaninitialinterfacelocation,temperaturedistribution,andouter bound-aryconditions,wewishtodeterminetheinterfacelocation (t)andtemperature T(x;t). Taking intoaccounteects ofsurfacetension,kineticundercoolingand anisotropyone canformulateaStefanproblemin physicalvariables,

8 < : @ t T = Dr 2 T; x2  (t); cD@ n T   + = Lv n ; x2 (t); T = T M e ()T M K =L v n =(); x2 (t) (2.1) where T M

is the melting temperature, D the thermal diusivity, c the heat capacity,L thelatentheat offusion,v

n

thenormalvelocityoftheinterface,K theinterfacialcurvaturetakentobepositiveifthecenterofcurvatureis inside the solid, () the interfacial mobility and e() a function depending on the surfacetension().

In (2.1)it has been assumedthat all physical parametersare constantand equalinthebulksolidandliquidphases. Thisisavalidassumptionforcertain materials, such assuccinonitrile(SCN)and will reducetheamountofnotation needed. Theassumptionis notvitalforthephaseeldmodeling[3].

Ω

Γ

Ω

(t)

+

-n

^

θ

(t )

( )

t

Figure2.1. Schematicoverviewoftheconsideredgeometry

The rst equation of (2.1) describes the diusion of heat inside the bulk solid and liquidphases. Thesecond equation givesanexpression forthe inter-facialvelocity, proportionalto thediscontinuityin thenormalderivativeof the temperatureacrosstheinterface. Itisneededforconservationofenergyand cor-respondstoreleaseand/orabsorptionoflatentheat. Thethirdequation,orthe Gibbs-Thomsoncondition,modelsthechangeofmeltingtemperatureaccording to kineticandcapillaryeects.

TheGibbs-Thomsonconditionof (2.1)isnotpresentin theclassicalStefan problem often seen in mathematical textbooks. Instead the third equation is replaced by T = T

M

, i.e. it is assumed that the phase transition takesplace at equilibrium. Themorecomplicatedboundaryconditionisneededinorderto obtaindendriticstructures. Inan attempttomotivatewhytheclassicalStefan problem is not suÆcientin thecurrentcontext, and explain why the dendritic structures occur, theunderlying physicsof the\pattern selectionprocess"will bebrie y overviewed. Moredetailed descriptionsof theunderlyingphysicscan befoundin[27,28].

A pure liquidis held at a temperaturebelow meltingtemperature. This is ametastable stateand solidication can beinitiated by disturbingthesystem, e.g. by adding a crystal seed. As solidication proceeds latent heat will be released at the solid-liquid interface, and the growth rate will depend on how fast this heat can be conducted away. From this point of view it is favorable to havealargeinterfacialareasince itimpliesfaster dissipation. Ontheother hand, the interface is associated with a surface energy which would in turn imply that asmall interfacialareaisto prefer. Theactualshapeof thecrystal is determined by a compromise between these two competing eects. Surface energyhasastabilizingeectonthegrowthsinceitdeterminesasmallestscale in the geometry. The processis unstablein thesense that small perturbations in initialdatacan signicantlyalterthetime-dependentshapeofthecrystal.

2.2 Related physical problems

TheStefan problem is signicantfor applications in semiconductorproduction where one is interested in the growth of pure crystals. In order to simulate realistic casting(2.1) needsto becomplementedwith equationsto accountfor thediusionofalloyingelements. Foreach additionalcomponentinthealloya diusion equationforaconcentrationisadded tothe Stefanproblem, andalso an interface condition prescribing a discontinuity in the concentration across the interface. The phase eld approach has been successfully used to model solidication of binary alloys, both in isothermal [42, 43] and non-isothermal settings[10,30].

For casting applications it may be unrealistic to assume that the melt is stationary during the solidication process. Phase-eld models have been for-mulatedwhich takeintoaccount uid owby accountingfor convectiveeects inthederivationofthemodel[4,37, 38].

2.3 Direct numerical simulation

Despite the apparentsimplicityof equation (2.1) analytic solutionsare known only for afew very special cases. For example aplanarinterface movingwith constant velocity for undercooling
= c(T

M T 1 )=L > 1, where T 1 is the uniform temperature of the liquid far away from the interface. The Ivantsov parabolasare amoreusefulfamilyof analyticsolutionsin twoand threespace dimensions. Intheabsenceofsurfacetensiontheygiveaparaboloidaldendrite tipandtheassociatedtemperatureforagivenundercooling. TheIvantsov solu-tionsareformulatedin termsoftheproductbetweentipvelocityandtipradius (the Peclet number) [28]. They are however not capable of determining the velocityor curvatureseparately.

Numericalcomputationsisthereforeanecessarytoolformodeling solidica-tionin realisticsettings. The majordiÆculties of direct simulation are associ-atedwithkeepingtrackofthemovingsolid-liquidinterfacessothatthecorrect boundaryconditionsareprescribed.

Front-trackingmethodstnaturallyintheframeworkoffreeboundary prob-lems, see e.g. [2,20]. The interface is representedby adiscrete set of marker points which are connected by interpolants to obtain the location of the free boundary. Thetime-dependentproblemisthensolvedbyadvectingthemarker points accordingto thevelocitygiven by(2.1), and solvingtheheat equations insolid andliquidgiventhelocationoftheboundary.

In[20]anapproachbasedontheimmersed boundarymethod[35]isusedto transfer data between a xed temperature grid and interfaces. Time-stepping essentially reduces to solving a system of non-linear equations for the normal velocities ofthemarkerpointsguaranteeingthat theGibbs-Thomson condition

holds. Almgrenconsidersadierentalgorithmin[2]wheremarkerpointsare ad-vectedbyminimizinganenergyfunctionalchosensuchthattheGibbs-Thomson conditionissatised.

Dueto thecomplicatedinterfacialshapeswhichariseitwill benecessaryto add anddelete pointsduring simulation. Merginginterfacespresent aproblem sinceitbecomesnecessaryto redenetheconnectivityofmarkerpointsto rep-resentthe newgeometry. This becomes increasingly diÆcult withthe number of spacedimensions.

Overthepastdecade phase-eldmethods haveemergedasapopular alter-nativeforcomputation. Usingideasfromcontinuummodelsofphasetransitions asetofreaction-diusionequationsvalidintheentire computationaldomain is derivedfromthermodynamics. Thesharpinterfacesof(2.1)arereplacedbythin diuse interfacesof width  implicitlygiven bythe solution. Thusthe needto explicitlyrepresentmovingboundariesiseliminated.

Level set methods can be used as an alternative to phase-eld modeling [11, 33]. Front-trackingisavoidedbyintroducingasigneddistance function,or levelsetfunction,whichforeverypointinthedomaingivesthedistance tothe closestsolid-liquidinterface. Incomputation thelevelsetisadvectedaccording to anextensionoftheinterfaceconditions. Thetemperatureisupdatedsothat the Stefan problem is satised by requiring that the temperature interpolant satises theGibbs-Thomsonconditionatsolid-liquidinterfaces[11].

Comparedto phase-eld methods thelevel set approach hasthe advantage of being independent of the articial diuse interface width . Simulation is performed using the sharp interface formulation so there is no need to derive a new set of equations. There are however some additional technicalities in the formulation of numerical schemes. Interpolation is needed to ensure that the interfaceconditions hold. It is also necessaryto \reinitialize" thelevelset function sothat itremainsasigneddistancefunction[11].

Todatephase-eldmethodshavebeenusedmuchmoreextensivelythanlevel setmethods. Thephase-eldmodelconsideredinthisworkisbrie yoverviewed in Chapter3.

The phase-eld approach

3.1 Derivation of a phase-eld model

There are several dierent phase-eld models describing the solidication of a pure material from its undercooled melt. They are all based on the key as-sumptionthat thereis asmoothapproximationof anindicatorfunction taking on distinct valuesin the bulk solidand liquid phases, the socalled phase-eld variable.

Thederivationproceedsbyformulatingthefreeenergy/entropyofthesystem and requiring that the temperature and phase-eld evolve such that the free energy decreases [29] or the entropy increases [34, 40]. The latter approach resultsinthermodynamicallyconsistentphase-eldmodelsaprioriguaranteeing that the second law of thermodynamics holds. It is notclear from simulation whether the thermodynamically consistentmodels give betterresults thanthe ones basedonderivationfrom freeenergy[24]. Sincenon-isothermalsituations areconsideredtheentropyishowevermoreappropriatefromaphysicalpointof view[34].

For both approaches a set of reaction-diusion equations is obtained but with dierent non-linear reactionterms. The dierence is also aected bythe choice of constant valuestaken by thephase-eld, which has consequences for theinterpolantsG()andp()discussedbelow.

A brief outline of the derivation of a thermodynamically consistent phase-eld model is now given. The presentation is based on the work of Wang et al. [40]wherethederivationiscarriedoutinfulldetail,andto lesserextenton themoregeneralworkofPenroseandFife[34]. Toavoidsometechnicalitieswe considertheisotropiccasewheretheangulardependenceofe()and()isnot takenintoaccount. Thesameapproachworksalsofortheanisotropiccase[32], asdiscussed inSection3.4.

Webeginbypostulatingtheexistenceof aphase-eldvariable, whichis an approximationofanon-conservedindicatorfunction suchthat

(x;t)=  0 x2 (t) 1 x2 + (t) ; (3.1)

c.f. Figure2.1. Furthermore,itisassumedthat(x;t) variessmoothlybetween 0and 1onlyin thin regions



(t)of width O()correspondingto diuse solid-liquidinterfaces. Thesmallpositiveparameterisrelatedtophysicalquantities below.

To derive evolution equations for phase-eld and temperature, we assume that the entropy and internal energy of any given subvolume V  can be written as S = Z V " s(e;) Æ 2 2 r
2 # dV; (3.2) E = Z V e(T;)dV; (3.3)

respectively. Here s is an entropy density and e an internal energy density. The second term in the integrand of (3.2) is a gradient entropy term needed to account for contributions from the solid-liquid interface, analogous to the Landau-Ginzburg energy term in the free energy [34]. According to the rst and second laws of thermodynamics, we assume that T and  evolve in such a way that the energy is conserved and the local entropy production in V is non-negative. Thesecondrequirementismetif

@ t = ÆS Æ ; (3.4)

where >0canbeinterpretedastherelaxationtimeofthephase-eld[34,40]. Toformulateasetofpartialdierentialequationstheinternalenergydensity is written e(T)=e (T)+p()L(T)=e + (T)+ p() 1)L(T); (3.5) where e 

(T) arethe classicalinternal energydensities in liquid and solid, and the latentheat L(T)=e (T) e

+

(T)istheir dierence. Thefunction p() is at this point unknownbut will bespeciedbelow. It isrequired to satisfythe normalizationp(0)=0andp(1)=1. Itwillbeassumedthattheinternalenergy densityin theliquidphaseislinearlydependentonthetemperature.

Toidentifythefunctionalderivativein (3.4)theinternalenergyisrelatedto theHelmholtzfreeenergythroughthermodynamicrelations[34],

f(T;)=T Z T TM e + ()  2 d p() 1
Q(T)+G() ! ; (3.6)

where Q(T)= Z T TM L()  2 d; (3.7)

and G()is anotherfunction tobespecied. Thefree energyisrequiredto be continuous with respect to temperatureat themelting point T

M

, andto have localminimainthesolidandliquidphases. Weobtaintheadditionalconstraints

G(0)=G(1);  G 0 () p 0 ()Q(T)  =0;1 =0;  G 00 () p 00 ()Q(T)  =0;1 >0: (3.8)

From (3.2){(3.6), and the rst and second law of thermodynamics the dy-namicequationscanbewritten

@ t =Q(T)p 0 () G 0 ()+Æ 2 r 2 ; @ t T+ 1 c (p() 1)L(T)
=Dr 2 T: (3.9) Expressionsforthefunctionsp()andG()areneededtofullyspecifythemodel. Following[40,ModelI]wechooseG()tobeasymmetricdoublewellpotential,

G()= 1 4a  2 (1 ) 2 = 1 4a g(); (3.10)

whereaistheinversedepthofthepotential,and

p()= R  0 g()d R 1 0 g()d = 3 (6 2 15+10): (3.11)

These functions satisfythe necessarynormalizations and(3.8) forall tempera-tures. Hencethefreeenergy(3.6)haslocalminimaat=0;1forallvaluesofT. Thiscorrespondstooneofthesolidorliquidphasesbeingstable,andtheother phasemetastable,regardlessofthesupercoolingoftheliquidorsuperheatingof thesolid. The choice ofp()hasthe advantagethat theconstant valuestaken in thebulkphasesare givenby (3.1)regardlessofT [34], but is notphysically realisticforT farfrom meltingtemperature.

Wenotethat byletting thelatentheatbeafunction ofT in (3.5)theheat capacities in the solid and liquid phases are allowed to depend dierently on temperature. Forsimplicityitwill nowbeassumed thatthelatentheat is con-stant,and thatL(T)=L(T

M )=L

0

,whichcorrespondsto equalconstantheat capacitiesc =c + =c. Equation(3.7)becomes Q(T)= L 0 T M " T T M T M +O  T T M T M  2 !# ; (3.12)

after TaylorexpansionaroundT =T M . IfjT T M jT M equations(3.9)and (3.12)resultinthelinearapproximation

@ t = L 0 T 2 M (T T M )p 0 () 1 4a g 0 ()+Æ 2 r 2 ; @ t T+ L 0 c p()
=Dr 2 T: (3.13)

SincethelinearizationisrelevantonlyforT closetoT M

itisnotalimitationto havelocalminimaofthefreeenergy(3.6)foralltemperaturesT.

The key assumptionsneeded for deriving thephase-eld model are the ex-istence of a phase-eld variable dened by (3.1) evolving according to (3.4), and asystem entropy on the form (3.2). All other assumptions can either be generalizedorfollowfrom standardphysicalmodeling[34].

Oneshould notethat a,Æand introducedduring thephase-eldmodeling are notphysically measurable quantities. Byconsideringthe phase-eldmodel under equilibrium conditions T =T

M

it is possibleto relatea and Æ to inter-face thicknessandsurfacetension. An asymptoticanalysis ofthe phase-eld modelinthelimitofvanishinginterfacethicknessthengivesanexpressionfor. Furthermore,theasymptotics showthat (2.1)and (3.13)are indeed equivalent in thelimitofvanishinginterfacethickness.

Detailsontheasymptoticanalysisofthis phase-eldmodelcanbefoundin [40]. In[32]anasymptoticanalysisofthemodelisgivenforanisotropicsurface energyandmobility. Caginalphasshownhowphase-eldmodelscangiveriseto other freeboundaryproblemsbyconsidering dierentasymptoticlimits[8]. In [9]convergenceto(2.1)isprovedrigorouslyintheone-dimensionalandradially symmetriccases.

Asymptoticsisanecessarytoolintheformulationofphase-eldmodelsalso for related problems, such asthesolidication of binary alloys [10, 43]. It has proved useful also for improving the eÆciency in numerical computation, as discussed inSection4.2.

3.2 Non-dimensionalization

Basedonresultsfromasymptoticanalysisthephase-eldmodelcanbewritten on dimensionless form to reduce the numberof parameters. Weintroduce the dimensionless temperature u= T T M T M T 1 ; where T 1

is somereferencetemperaturesuch astheinitial temperatureonthe boundaries. Space andtimearescaledby

x= x dim w ; t= t dim w 2 =D ;

wherethesubscriptdimdenotesadimensionalquantity,wisareferencelength scale,andw

2

=Disthethermaldiusiontime. Deningdimensionless variables

= p aÆ w =  dim w ; (3.15) m= T m DL 0 ; (3.16) = wL 2 0 6 p 2cT M ;
= c(T m T 0 ) L 0 ; (3.13)canberewrittenas  2 m @ t =
up 0 () 1 4 g 0 ()+ 2 r 2 ; (3.17) @ t  u+ 1
p()  =r 2 u: (3.18)

isherethedimensionlessundercooling(orinverseStefannumber)which typ-icallyliesbetween0and1. 
canbeinterpretedasascaledinversecapillary length, and m is the interfacial mobility. The diuse interface thickness is of orderas givenby(3.15).

We note that (2.1) can also be expressed using the dimensionless scalings introducedabove. TheStefanproblemtakestheform

8 > < > : @ t u= r 2 u x2  (t) @ n u   + = v n =
; x2 (t) u= 1 6 p 2
 K+ vn m  x2 (t) ; (3.19) wherev n

andKshould beinterpretedasdimensionless quantities.

3.3 Model properties

Thephase-eldmodel(3.17){(3.18)andtheequivalentdimensionalformulation (3.13)arevalidintheentirecomputationaldomain. Equation(3.18)corresponds to conservation of energy. Thenon-linear termsare neededto account forthe release/absorptionoflatentheatatthesolid-liquidinterface. Equation(3.17)is the evolutionequation for the phase-eld, and doesnot haveadirect physical interpretation.

Sincethephase-eldmodelisvalideverywherethereisnoneedtoexplicitly track geometrically complex solid-liquid interfaces. The location is implicitly givenbythedenitionofthephase-eldvariable. Onecouldforexampleletthe locationofthefreeboundarybedened by (t)=fy(t):(y(t);t)=1=2g,and thediuseinterfacebegivenby



neitheroftherepresentationsareneededduringcomputationtheyaremostlikely of interestduring post-processingto determinetheshapeofthecrystal.

Equivalencebetweenthephase-eldequationsandthesharpinterfacemodel isformallyshowninthelimit!0,thusresultinginasingularlyperturbedset of reaction-diusionequations. Forpurposesofcomputation itisimpossibleto let !0sincethegridresolution mustbeat leastof order inthe vicinity of interfacestoresolvethetransitionoffromzerotoone. EventhoughtheStefan problemisonlyanapproximationofaphysicalprocessitisimpossibletochoose a physically realisticvalueof in simulation. Interfaceshaveanite width on the scale of Angstromswhich can for examplebe comparedto theradius of a SCN dendritetipwhichisontheorderof0:1mm[28].

Thereisalsoanupperboundonthevalueofdeterminedbythedynamical properties of the phase-eld model. If a toolarge value is used in simulation unphysical oscillatoryintermediatestatesappearin thetransition regionofthe phase-eld variable, and themodel losesitsmeaning [1, 12, 14]. These restric-tionsoncanberelaxedbymoreaccurateasymptoticanalysis,asdiscussedin Section 4.2,andpapers3and 4.

Themaindrawbackofthephase-eldmethodsisthat theaccuracydepends on the parameter . Inpractice thevalue of  is determined as acompromise betweendesiredaccuracy andavailable computational resources. We shall dis-cuss this issuefurther in Section 4.2 where implicationsof asymptoticanalysis oncomputationaleÆciencyis considered.

Thereareanumberofother dieringlengthandtimescalesin theproblem which make the formulation of eÆcient numerical schemes challenging. The spatial discretizationmustbechosensuchthatthesmall-scalelocalgeometryis resolved while the gridshould coverasuÆciently large domain to capture the entirecrystal. Atlowundercoolingsthetemperatureeld extends muchfaster than thesolidicationfrontpropagates[6]. Toensurethatcomputed dendrites areunaectedbyboundaryconditionsimposedontemperatureitisnecessaryto discretize largedomains. Since thecrystalgrowthis slowlong-timeintegration becomesnecessary.

Thetimesteprestrictionsaregovernedbyconditionsinthevicinityof inter-faces. Changesawayfrom interfaces,where the phase-eldmodelreduces to a standardheatequation,occuronamuchslowertimescale. Choosingtimesteps becomes evenmorediÆcult foralloysolidication where thediusion times of heat andsolutecanvary greatly.

3.4 Inclusion of anisotropy

For manymaterials, includingmetals,the surfaceenergyof thesolid-liquid in-terfaceandthekineticcoeÆcientdepend onorientation. Theseeects werenot accountedfor in thederivation outlined in Section 3.1. Since anisotropyhasa crucialimpactontheshapeofthedendritesitisnecessarytomodifythe phase-eld modelforsimulationinrealisticphysicalsettings. In[32]McFaddenet. al.

derive a phase-eld model with anisotropy in twospace dimensions using the sameapproachasin Section3.1.

Letting the parameters and Æ from (3.9) depend on  (c.f. Figure 2.1) a modelcanbeformulatedwhichisequivalentto theStefanproblem(2.1)inthe limit!0. Themodiedversionof (3.9)becomes

()@ t =Q(T)p 0 () G 0 ()+  Æ 2 r 2 a (); @ t T+ 1 c (p() 1)L(T)
=Dr 2 T;

wheretheordinaryLaplacianoftherstequationhasbeenreplacedby

r 2 a ()= @ @x  () 0 () @ @y  + @ @y  () 0 () @ @x  +r  2 ()r
: (3.21)

Here()isafunctionchosentomodeltheangulardependenceofthecoeÆcient in front of the Landau-Ginzburg term in the entropy functional (3.2). It is typicallychosenas Æ()=  Æ()=  Æ[1+ cos(k)]; (3.22) in the notationof [44]. isthe relativestrengthof theanisotropy,the integer k amodenumberdeterminingthenumberofsymmetryaxes,and



Æaconstant governingthesizeofthegradientcontributions.

Oneofthemainreasonsforusingthephase-eldapproachistoavoidexplicit representationofinterfaces. Thereforetheangle mustbeformulatedinterms ofthephase-eldvariable. Wedothisbyidentifyingtheinterfacewithlevelsets ofthephase-eld. Theunit normalcanthenbewritten

b n= r krk 2 =cos(  )bx+sin(  )by; where 

istheanglebetweenbnandthex-axis. Iftheanglebetweenthereference directionandthex-axisis

0

,itfollowsfromsimpletrigonometrythat

=    0 =tan 1 @ y  @ x  !  0 : (3.24)

ThemodiedLaplacian(3.21) isclearlynon-linear in ,whichis signicantin theformulationofnumericalschemes.

In the dimensionless phase-eldmodelthe choice of () introduces an an-gular dependence also in the kinetic undercoolingterm of theGibbs-Thomson condition [6, 32]. Byletting () :=(), where () is given by (3.22), the contributionbecomesindependentoftheanglebetweenthesolid-liquidinterface and the reference direction [6]. If weinstead choose () = the only modi-cation of (3.17)is that theLaplacianshould bereplacedby (3.21). Equation (3.18)doesnotneedto bechangedtoaccountforanisotropy.

Numerical simulation

AfewdierentapproachestodirectsimulationoftheStefanproblemwerebrie y overviewedin Section 2.3. The phase-eld approach hasbeen gainingin pop-ularityand is emergingas the method of choice for simulationof dendritic so-lidication. Thereis anextensiveamountof literatureontheuseofphase-eld models in computation. Theyhavebeenused to model anumber of dierent physical phenomena, including solidication of pure materials [22, 25, 44], so-lidicationofalloys[5,10, 30,42],solidtosolidphasetransformations[41]and solidication under external ow [39, 38] just to name afew. Several dierent numericalapproximationschemeshavebeenused. Althoughuniformnite dier-encediscretizationsseemtobethemostcommonthereareexamplesofadaptive algorithmsbasedonbothnitedierences[6]andnite elements[36,39].

In this chapter some aspects of numerical approximation will be covered. Theaim isnotto fullycovernite dierencetechniquesbut to introducesome conceptswhich arediscussedingreaterdetailin thepapersbelow.

4.1 Finite dierence methods

There are anumberof generaltechniquesfor numericalapproximationof par-tial dierentialequations (PDE),including nite dierence, nite elementand spectralmethods. Theyallresultinasetofalgebraicequations,whosesolution constitutes anapproximationto the solution of thecontinuous problem. Here wehavechosentoworkwithnitedierencemethods,whichleadtostructured systemsapproximatingthePDE. Thus hierarchicalcomputermemory systems canbeeÆcientlyused. Identicationofdatadependenciesbecomestrivialwhich isusefulwhenformulatingparallel algorithms,asdiscussedin Section4.3.

Finitedierencetechniquesareconceptuallysimple. Algebraicequationsare obtained by replacing partial derivativesby dierence quotas. By solvingthe system anapproximate valueof the solutionis obtainedin everypoint onthe computational grid. For time-dependent PDE wethink of the computational

gridasadiscreterepresentationofspaceonly. Starting fromsomeinitial state, thesolutionisadvancedintimebysolvingasetofalgebraicequationspertime step.

Numericalschemescanbeconstructedviathemethodoflines. Byreplacing spatial derivativeswithdierencequotasasetofordinarydierentialequations (ODE)for thegridvaluesisobtained. Onecanthenproceedwiththetemporal discretization using standard numerical methods developed for ODE. We dis-tinguish between implicit and explicit time-stepping schemes. For the latter the algebraicsystemsare trivialandthesolutionat thenexttime levelcanbe explicitly formulatedusing the current approximation. Forimplicit schemes it is necessaryto actually solvea systemto update the solution. In thepresent application thesystemsarenon-linear.

For analysis of nite dierencemethods we areprimarily interestedin sta-bilityand accuracy. Stabilityimposes conditionson the relationbetweentime stepandgridspacingwhichmustholdtopreventuncontrolledgrowthoferrors. TheaccuracyisameasureofhowwellthediscretizationapproximatesthePDE and hencealsoatwhichratetheapproximationconvergesto theexactsolution asgridspacingandtimestepdecrease. Pleasereferto[19]forproperdenitions and adetailed expositionontheanalysis ofnitedierencemethods.

In computation onewould like to determine an approximate solution to a given tolerance as fast as possible. The stability and accuracy of numerical schemeshaveacrucialimpactsincetheygoverntherangeoftimestepsandgrid spacings resulting in thedesired accuracy. Forexplicit schemesthe simulation timeincreaseswiththenumberoftimestepstaken,andtheexecutiontimeper time stepgrowslinearlywiththenumberofgridpoints. Thisisnotnecessarily trueforimplicitschemes,whenasystemofequationsmustbesolvedeachtime step. Thestabilityconditionsforimplicitschemes arehoweveroftenmuch less restrictive.

Inpaper1twodierenttimesteppingschemesarecomparedwithrespectto achievedaccuracyandtotalexecutiontime.

4.2 Asymptotics and implications on computational eÆciency

Asymptoticanalysisofphase-eldmodelsisnecessarytoshowequivalencewith correspondingsharpinterfaceproblems but alsohassignicantimplications on computationaleÆciency. ThiswasrstobservedbyKarmaandRappel[21]and hassinceattractedconsiderableattention,seee.g. [3,13,22,31].

Theearlyanalysisveriedthat solutionstothephase-eldequationssatisfy the Stefan problem with errors of size O() [7, 8, 9, 40]. Thus convergence is guaranteed in the limit  ! 0. More accurate asymptotics have subsequently determinedtheinterfaceconditionstoO(

2

),andrequirementsonthephase-eld modelforsecondorderagreementwiththeStefanproblem[3,13,21]. Sincethe

gridresolutionisdeterminedbythisallowssimulationonsmallergridswithout increasingerrorsintheinterfaceconditions.

KarmaandRappelformulatedtheobservationsintermsofthesharp-interface andthin-interfacelimitswhichhavemoredirectphysicalinterpretation. Inthe rstorderasymptotics,orthesharp-interfacelimit,thetemperatureisconstant inthediuseinterface. InthesecondorderasymptoticstheO()-contributionto thetemperaturevariesacrosstheinterface. Takingthecontributionintoaccount canthereforebeseenasallowingasmall variationofu,or equivalently,having athininterfaceofwidthdierentfromzero.

Theiranalysisshowedthatthephase-eld modelgivesanerrorof sizeO() in the Gibbs-Thomson condition (3.19) because the mobility m is determined by rst order asymptotics. To obtain second order agreement it is therefore necessarytoaddacorrectionpriortosimulation. Forthemodelconsideredhere thecorrectiontakestheform

1 m = 1 m 0 + 209 35 ; (4.1) wherem 0

isthevaluepredictedby(3.16).

Itcanbeshownthatsimulationwithoutapplying(4.1)isequivalentto assum-ingadiuseinterfacethicknessmuchsmallerthanthecapillarylength[21,22]. Since the capillary lengthis a measure of the smallestgeometric scale, the ef-fectsofthethininterfacecorrection(4.1)aresignicant. Computationonsmaller gridsismadepossiblebymodifyingasinglescalarparameterpriortosimulation. Alargervalueofcanalsobeusedwithoutintroducingunphysicalintermediate states in the phase-eld variable. Furthermore, the correction allows simula-tion in the absence of kinetic undercooling (m

0

= 1), which is a reasonable assumptionatlowundercoolings[3].

Theexactformof (4.1)dependsonpropertiesofthefunctionsp()andg() inthephase-eldmodel. Almgren[3]hasshownthatthefunctions mustsatisfy certainintegralconstraintsforathin interfacecorrectionto exist. His analysis wasperformedforthemoregeneralproblemofbinaryalloysolidication.

4.3 Parallel computing

Parallel computing is a cost eÆcient way of overcoming limitations of current hardwaretechnologysuchasprocessorclockfrequencyandmemorybandwidth. By letting several computers cooperateto solvea single computational task a problemcanbesolvedfasterthanasingleconventionalcomputerwouldallow. It alsoenableslargerproblemsizes,whichcanbeexploitedtoincreasetheaccuracy orto simulatemorecomplexsystems.

There are several types of parallel architectures. The two most common ones are the sharedand distributed memorymachines which dier in the way computer memory is handled. On a shared memory computer all processors share thesameaddressspace and haveaccess to alargeglobal memory. On a

distributedmemorycomputereachprocessorhasaccesstoitsownlocalmemory only.

Hereweconsiderdistributedmemoryarchitectures. Sinceeachprocessorhas itsownlocalmemorythebandwidthscaleswiththenumberofcomputersused. Non-localdata can howevernotbeaccessed directly. Thereis aneed for com-munication-ormessagepassing-betweenprocessors. Oneofthecriticalissues whendesigningaprogramistokeeptheratioofcommunicationtocomputation small sincecommunicationisessentiallyanundesirableoverhead. Forthis rea-son itmaybenecessaryto choosedierent algorithmsthanonewould haveon conventionalcomputerarchitectures.

Wedistinguishbetweenlocalcommunicationinvolvingafewprocessors,and globalcommunicationinvolvingallprocessors. Onewouldliketoavoidthelatter if possiblesince itmeans that all computers must reach a certainstage in the computation before the execution in any processor can proceed. At these so called synchronization points all processors remain idle while waiting for the slowest one to complete its task. Thus it is important to have a good load balance so that an equal amount of work is performed on all processors. By choosingan appropriatedistributionofdata onto theparallel machine thiscan be achieved. The data distribution also aects the amount of communication needed.

4.3.1 Scalability metrics

The idealparallelprogramrunning onpprocessorsshould nishptimesfaster than onasingleprocessor. Due tocommunicationoverhead thisis hardlyever the case. To measure how eÆciently a parallel algorithm exploits the com-putational resources of a parallel machine one therefore investigates dierent scalabilitymetrics[26]. Themostwellknownisthespeedupdened by

speedup= T 1 (W) T(p;W) : HereT 1

(W)isthetimeneededtosolveaproblemofsizeW onasingleprocessor using the fastest known serial algorithm,and T(p;W)the time neededfor the parallel algorithmrunningonpprocessorstosolvethesameproblem

1 .

TheparalleleÆciencyisdenedasthespeedupoverthenumberofprocessors. A paralleleÆciency ofonecorrespondsto aperfectlyscalablealgorithm, which runs p times faster on p processors. For an algorithm which is not perfectly scalablethereisanasymptoticlimitonthespeedupwhichcanbeobtained. The famous Amdahl's lawstatesthat if afraction f ofthe code isinherentlyserial andcannotbeparallelizedthespeedupcanneverexceed1=f,evenonaninnite numberofprocessors[18].

1

Theproblemsizeismeasuredasthenumberof oatingpointoperationsrequiredbythe fastestknownserialalgorithmrunningonasingleprocessor.

Thedenition ofthespeedupis basedontheassumptionthat aproblemof agivensizeW istobesolved. Itdoesnotaccountfortheabilitytosolvelarger problems on a parallel machine. The scaled speedup is an alternativemetric, whichcompareshowlongagivenparallelprogram wouldhavetakento runon ahypotheticalserial processorwith thesame amount of RAM memoryasthe parallelmachine[18]. Theproblemsizeisincreasedlinearlywiththenumberof processors,andthescaledspeedupisdened by

scaledspeedup= T

1 (W
p) T(p;W
p)

Byincreasing the problem size memory and surface-to-volume eects are kept constant. The limitations of Amdahl's law no longer apply and the parallel eÆciencyapproaches1 f asthenumberofprocessorstendstoinnity[18].

Whatkindofscalabilitymetrictousedependsonwhatonewishestousethe increasedcomputationalresourcesfor. A morethorough coverageof scalability analysisandothermetricscanbefoundine.g. [16,26].

Computed dendrites

In this chapter three examplesof computed dendrites are given. They are in-tended to exemplify the type of dendritic structures observed in simulation. Hence only computational parameters and plots are given and no conclusions are drawn. Additional simulations can be found in papers 1, 2 and 4 where dierentaspectsofthephase-eldmodelareinvestigatednumerically.

All simulationsin this chapter usedimensionless parameters=400,
= 0:5 and m= 0:05,corresponding to pure nickel at acertain undercooling[44]. Thediuseinterfacethicknessistakento be=0:005. Anisotropyis modeled by(3.22)with =0:015,andk=4or6dependingonsymmetry. Fromphysics onewouldexpectthefourfoldsymmetrytobemoreappropriatesincethatbetter approximatestheBCC-latticeobservedin solidnickel.

Theplotsshowningures5.1{5.3havebeencomputedusingthesemi-explicit scheme described in paper1. The uniform discretization satises
x=
y= =0:005,andaconstanttimestep
t=10

4

hasbeenemployed. Thegridsize wasincreasedduringthecourseofthesimulations.Whenthetemperatureatan outerboundarydiersfromtheinitialvalue 1bymorethan10

3

anadditional 50gridpointsareaddedinthatdirectiontoensurethatboundaryconditionsdo notaect the growthof thedendrite. This adaptivityis ofcourse ad-hoc,and likelytoin uencethecomputedsolution. Sincetheonlypurposeistoexemplify computeddendritesithasneverthelessbeenusedtoreduceexecutiontime. Each simulationtook aboutvedaysonasingleworkstationwithalargeamountof RAMmemory.

In addition to the phase-eld equations described in Chapter 3 a random perturabtion of 1% amplitudewasadded to the non-linearterms in the right-hand side of equation. This serves to stimualte sidearm growth, and can be physicallyinterpretedasaphenomenologicalmodelformaterialimpurititesand thermalnoise(3.18)[23,44].

Thesimulationshavenotbeencompared toactual nickeldendritesbut the dendriticstructuresobservedexhibitthequalitativebehaviourobservedin phys-icalexperiments,seee.g. [15].

Figure5.1.Solidicationofpurenickelunderfourfoldanisotropystartingfrom circularinitialseedofradius0:05. Thesolidlinesshowthelocationofthe solid-liquidinterfaceattimelevelst=0;0:3;0:6;0:9;1:2.Thebackgroundcolor repre-sentstheheatdistributionatthenaltimet=1:2.

Figure5.2. Solidicationofpurenickelundersixfoldanisotropystartingfrom circularinitialseedofradius0:05. Thesolidlinesshowthelocationofthe solid-liquidinterfaceattimelevelst=0;0:3;0:6;0:9;1:2.Thebackgroundcolor repre-sentstheheatdistributionatt=1:2.

Figure5.3. Solidicationofpurenickelstartingfromtwocircularinitialseeds of radius0:05. Thesolidlinesshowthe locationofthesolid-liquidinterfaceat time levelst=0;0:25;0:5;0:75;1:0. Thebackground colorrepresentsthe heat distributionatt=1:0.

Summary of papers

6.1 Paper 1: Time-stepping schemes for phase-eld simulation of dendritic solidication

Inthispapertwodierentapproachesto time-steppingareinvestigatedforthe phase-eld model discretized by nite dierences on uniform Cartesian grids. Therst schemeis semi-explicit, computationallycheappertime stepandhas low storage requirements. Since the phase-eld variable is updated using an explicit method there are however step size restrictions. The second scheme is based onthe implicit BDF-methods. It leads to increased accuracy and re-laxedstepsizerestrictionsattheexpenseofincreasedstoragerequirementsand computationalcomplexity.

Forthesemi-explicitschemeone needstosolvealinearsystemofequations each time step, which can bedone eÆciently using a fastPoisson solver. The implicit scheme leadsto non-linearsystemssolvedbya Newton-GMRES algo-rithm. Solversare discussed in some detail. Extensions to adaptivegrids and parallelcomputersareconsidered,althoughneitherareusedinthis paper.

Numerical experiments showing that the semi-explicit scheme is only rst orderaccurateintimearepresented. Convergenceratesofrst,secondandthird order BDF-schemes are experimentallyvalidated. The importance of accurate timediscretizationisshownbyasimpleone-dimensionalexamplewherethe semi-explicitschemefailstocapturethecorrectpropagationspeedoftheinterface.

Byexperimentsusingrealisticparametersitisshownthat onebenetsfrom usingtheBDF-schemeintwodimensionsforsuÆcientlylowerrortolerances.

6.2 Paper 2: Computation of dendrites on parallel distributed memory architectures

In paper1it was shownby computationalexperiments that onebenetsfrom using the implicit scheme which is more computationally expensive per time step for suÆciently low error tolerances. In this paper the possibility to use parallel computers to obtain simulation times comparable to the semi-explicit scheme also for largererror tolerancesis investigated. The programdescribed is implemented in C and uses MPI [17] for message passing in a distributed memory setting. Theimplementationisveriedbyreproducing acomputation previouslypublishedin [6].

ThekeycomputationalkernelsoftheNewton-GMRESsolverareidentiedin termsoflinearalgebraoperations. Adatadistributionissuggestedwhichresults in alowcommunicationovercomputation ratio. Thescalability isinvestigated onamassivelyparallelIBMSP2,aswellasnetworksofSunUltra5workstations whichdonothaveaccesstoanoptimizedcommunicationsnetwork. Theparallel eÆciency exceeds 75% on 16 IBM SP2 processors even for grids as small as 200200 points. For the same grid size an eÆciency of more than 70% is achievedusinguptoeightSuncomputers. Investigationsofthescaledspeedup demonstrate the possibilityof using parallel computers foreÆcient simulation onlargergrids.

6.3 Paper 3: Third order asymptotics of a phase-eld model

AsrstobservedbyKarmaandRappelin[21]theagreementbetweenthe phase-eld andsharpinterfacemodelscan beincreasedfromrsttosecondorderin  bymodicationofthedimensionlessparametermonly. Themodicationallows simulationonsmallergridswithoutlossofaccuracy. Inthispaperweinvestigate whether it is possible to obtaina third order accurate phase-eld model by a similar modication. Thefocusisthereforenotonthenumericalalgorithmsor implementation,as inpapers1and 2,but onthephase-eldmodelitself.

To identify the interface conditionssatisedby solutions to the phase-eld equations a formal asymptotic analysis is performed in the limit of vanishing diuse interfacethickness. Theasymptoticsiscarriedoutusingthe dimension-less quantities from Section 3.2 ratherthan thedimensional phase-eld model. Fortheanalysisthefunctionp()modelingthereleaseoflatentheatis unspec-ied, but requiredto satisfy certainconstraintsimposed byassumptions made in derivingthemodel.

In the process of constructing asymptotic approximations additional con-straintsonp()appearwhichmust be satisedforthetemperature tobe con-tinuous across the interface. Conditions for the existence of a thin interface

correction are identied. An explicitexpression for the correction resulting in secondorderaccuracyisgiven.

Nothin interfacecorrectionleadingtothird orderaccuracyhasbeenfound. Due to diÆculties evaluating the interfacial temperature for general p() we have restricted the investigation to all possible polynomials of degree nine or less. That allows us to determine the interfacial conditions by using a C++ programforsemi-symboliccomputation.

Fromtheanalysis andthesemi-symboliccodeweconcludethat ifthereare polynomialsp()leadingtothirdorderagreementtheymustbeofdegreetenor higher. Forcomputationonewould prefernotto usepolynomialsof such high degreeduetoround-oerrors.

6.4 Paper 4: Analysis of a phase-eld model for solidication in the limit of vanishing

undercooling

Inthelast paperanalternativeasymptoticanalysis ispresentedin thelimitof vanishing undercooling
. The limitcan be motivated by an interest to per-formnumericalsimulationinthelowundercoolingregime. Astheundercooling decreasesthesmallestgeometric scaleincreases,andthesolidication proceeds slower. Thuslargerdomainsmustbediscretizedbeforedendriticstructuresare observed, andlongtime integrationis needed forthem to arise. It is desirable to haveadiuseinterfacethickness&
toreducecomputationalcomplexity. Insimulationitisoftennecessarytouseasmallertoavoidtheappearanceof unphysicaloscillationsinthephase-eldvariable.

Fortheanalysisweassumethat
approacheszeromuchfasterthan. The constructionofasymptoticapproximationsisbasedonthesametechniquesused inpaper3. Fortechnicalreasonsstrongerassumptionsontheinterfacialvelocity andcurvaturearerequired. Asaconsequencethelimit
!0correspondsto a stationary planarinterface. Thelimititself is notinteresting from a computa-tionalpointofview.Theasymptoticsismerelyusedtoinvestigatetheagreement betweenthesharpinterfaceandphase-eldmodels.

It is shown that the Gibbs-Thomson condition has an error of size . By applyingthesamethininterfacecorrectionderivedinpaper3thisleadingorder errorcanbeeliminated,resultingin anerrorofsize
. Sinceitisassumedthat  approacheszero muchslowerthan
theorder of accuracyis increased from \almost zero"to one. The analysis presentsa motivation forthe usefulnessof thethininterfacecorrectionatlowundercoolingswhichdoesnotdirectlyrelyon scaling argumentsor theapproximationm

0

=1. Furthermore,thecorrection allowslargertobeusedbefore unphysical intermediatestatesappear.

Numerical experiments supporting the conclusions of the analysis are pre-sented.

[1] A.S. Almgren andR. F. Almgren, Phaseeld instabilitiesand adap-tive mesh renement, in Modern methods for modeling microstructure in materials,SIAM,1996.

[2] R.F.Almgren,Variationalalgorithmsandpatternformationindendritic solidication, J.Comp.Phys.,106(1993),pp.337{354.

[3] , Second-order phase eld asymptotics for unequal conductivities, SIAMJ.Appl.Math, 59(1999),pp.2086{2107.

[4] C. Beckermann, H. J. Diepers, I. Steinbach, A. Karma, and X. Tong, Modeling melt convection in phase-eld simulations of solidi-cation, J.Comp.Phys.,154(1999),pp.468{496.

[5] W.J.Boettiner,A.A.Wheeler,B.T.Murray,andG.B.McF ad-den,Prediction ofsolutetrapping athigh solidication ratesusingadiuse interfacephase-eldtheoryofalloysolidication,Mat.Sci.andEng.,(1994), pp.217{223.

[6] R. J. Braun and B. T. Murray, Adaptive phase-eld computations of dendriticcrystalgrowth,JournalofCrystalGrowth,(1997),pp.41{53. [7] G.Caginalp,An analysis ofaphaseeldmodelofafree boundary,Arch.

Rat.Mech.Anal, 92(1985),pp.205{245.

[8] ,StefanandHele-Shawtypemodels asasymptoticlimits ofthe phase-eldequations,Phys.Rev.A,39(1989),pp.5887{5896.

[9] G. Caginalp and X. Chen, Phase eld equations in the singular limit of sharp interface problems, in On the evolution of phase boundaries, M.GurtinandG.McFadden, eds.,vol.43ofTheIMA Volumesin mathe-maticsanditsapplications,Springer-Verlag,1992,pp.1{27.

[10] G. Caginalp and W. Xie, Phase-eld and sharp-interface alloy models, Phys.Rev.E, 48(1993),pp.1897{1909.

[11] S.Chen, B.Merriman,S.Osher,andP.Smereka,Asimple levelset methodfor solvingStefan problems. preprint,1996.

[12] U.EbertandW.vanSaarloos,Breakdownofthestandardperturbation theory andmoving boundary approximation for "pulled" fronts,Tech.Rep. MAS-R0003, National Research Institute for Mathematics and Computer Science (CWI),2000.

[13] P. C. Fife and O. Penrose, Interfacial dynamics for thermodynam-ically consistent phase-eld models with nonconserved order paramter, Electronic Journal of Dierential Equations, 1995 (1995), pp. 1{49. http://ejde.math.swt.edu.

[14] K.GlasnerandR.Almgren,Dualfrontsinaphaseeldmodel,Physica D, 146(2000),pp.328{340.

[15] M.Glicksman,Freedendriticgrowth, Mater.Sci.Eng., 45(1984). [16] A. Grama, A. Gupta, and V. Kumar, IsoeÆciency function: A

scal-ability metric for parallel algorithms and architectures, IEEE Parallel and DistributedTechnology,1(1993),pp.12{21.

[17] W. Gropp, E. Lusk, and A. Skjellum, Using MPI: portable parallel programmingwiththe messagepassinginterface,TheMITPress,1994. [18] J. L. Gustafson,G.R. Montry, andR. E. Benner,Development of

parallel methods for a 1024-processor hypercube,Siam J. Sci.Stat. Comp., 9(1988),pp.609{638.

[19] B. Gustafsson, H.-O. Kreiss, and J. Oliger,Time Dependent Prob-lems andDierence Methods,JohnWileyandSons,1995.

[20] D. Juric and G. Tryggvason, A front-tracking method for dendritic solidication,J.Comp.Phys.,123(1996),pp.127{148.

[21] A. Karma and W.-J. Rappel, Phase-eld method for computationally eÆcient modeling of solidication with arbitrary interface kinetics, Phys. Rev.E,53(1996).

[22] ,Quantitativephase-eldmodelingofdendriticgrowthintwoandthree dimensions, Phys.Rev.E,57(1998),pp.4323{4349.

[23] ,Phase-eldmodelofdendriticsidebranchingwiththermalnoise,Phys. Rev.E,10(1999),pp.3614{3625.

[24] Y.-T. Kim,N.Provatas,N. Goldenfeld,andJ.Dantzig,Universal dynamics of phase-eldmodels for dendriticgrowth. preprint,1998. [25] R. Kobyashi,Anumerical approach tothree-dimensional dendritic

[26] V.KumarandA.Gupta,Analyzingscalabilityofparallel algorithmsand architectures, Journal of Parallel and Distributed Computing, 22 (1993), pp.379{391.

[27] W. Kurz and D. J. Fisher,Fundamentalsof Solidication, TransTech Publications,1986.

[28] J. S. Langer,Instabilities andpattern formation in crystal growth, Rev. Mod. Phys.,52(1980),pp.1{28.

[29] , Models of pattern formation in rst-order phase transitions, in Di-rectionsin CondensedMatterPhysics,G.GrinsteilandG. Mazenko,eds., WorldScientic,Singapore,1986,pp.164{.

[30] I.Loginova,G.Amberg,andJ. 

Agren,Phase-eldsimulationsof non-isothermalbinary alloy solidication, ActaMaterialia, 49(2001),pp.573{ 581.

[31] G.B.McFadden, A.A.Wheeler, andD.M.Anderson,Thin inter-faceasymptotics foran energy/entropy approachtophase-eldmodels with unequal conductivities,PhysicaD, 144(2000),pp.154{168.

[32] G.B.McFadden,A.A. Wheeler,R.J. Braun,andS. R.Coriell, Phase-eld models for anisotropic interfaces, Phys. Rev. E, 48 (1993), pp.2016{2024.

[33] S.OsherandJ.A.Sethian,Frontspropagatingwithcurvature-dependent speed: Algorithms basedon Hamilton-Jacobi formulations,J.Comp.Phys., 79(1988),pp.12{49.

[34] O. Penrose and P. C. Fife, Thermodynamically consistent models of phase-eldtype for the kinetics of phase transitions, PhysicaD, 43(1990), pp.44{62.

[35] C.S.Peskin,Numericalanalysisofblood owintheheart,J.Comp.Phys., 25(1977),pp.220{252.

[36] N. Provatas, N. Goldenfeld, and J. Dantzig, Adaptive mesh re-nement computation of soldication microstructures using dynamic data structures,J.Comp.Phys.,(1998).

[37] R. Th 

onhardt and G. Amberg, Phase eld simulation of dendritic growth inashear ow, J.CrystalGrowth,194(1998),pp.406{425. [38] , Dendritic growth of randomly oriented nuclei in a shear ow.

ac-ceptedforpublicationbyJ.CrystalGrowth,2000. [39] R.

T

onhardt, Convective Eects on Dendritic Solidication, Ph.D. dis-sertation,RoyalInstituteofTechnology,December1998.

[40] S.-L. Wang, R.F. Sekerka, A. A. Wheeler, B. T. Murray, S. R. Coriell, R. J. Braun, and G. B. McFadden, Thermodynamical ly-consistent phase-eld models for solidication, Physica D, 69 (1993), pp.189{200.

[41] Y. Wang, L.-Q. Chen, and A. Khachaturyan, Computer simulation ofmicrostructureevolutionincoherentsolids,inSolid!SolidPhaseT rans-formations, W. Johnson, J. Howe, D. Lauhlin, and W. Soa, eds., The Minerals, Metals&MaterilasSociety,1994.

[42] J. A. Warren and W. J. Boettinger, Prediction of dendritic growth andmicrosegregationpatternsinabinaryalloyusingthephase-eldmethod, Actametall.mater.,43(1995),pp.689{703.

[43] A. A. Wheeler, W. J. Boettinger, and G. B. McFadden, Phase-eldmodel for isothermal phasetransitions inbinaryalloys,Phys.Rev.A, (1992).

[44] A. A.Wheeler,B. T.Murray,andR.J. Schaefer,Computationof dendrites usingaphase eldmodel,PhysicaD, 66(1993),pp.243{262.