ALevelSetApproachfortheNumericalSimulationofDendriticGrowth
FrédéricGibou,1RonaldFedkiw,2RusselCaflisch,3andStanleyOsher3
ReceivedApril24,2002;accepted(inrevisedform)October17,2002
Inthispaper,wepresentalevelsetapproachforthemodelingofdendriticsolidifica-tion.ThesesimulationsexploitarecentlydevelopedsecondorderaccuratesymmetricdiscretizationofthePoissonequation,see[12].Numericalresultsindicatethatthismethodcanbeusedsuccessfullyoncomplexinterfacialshapesandcansimulatemanyofthephysicalfeaturesofdendriticsolidification.Weapplythisalgorithmtothesimulationofthedendriticcrystallizationofapuremeltandfindthatthedendritetipvelocityandtipshapesareinexcellentagreementwithsolvabilitytheory.Numericalresultsarepresentedinbothtwoandthreespatialdimensions.
KEYWORDS:Interfaces;levelsetmethod;ghostfluidmethod;Poissonequation;dendriticgrowth.
1.INTRODUCTION
Variousnumericalmethodshavebeendevelopedtosolvethedifficultproblemsassociatedwithdendriticcrystallization.Broadlyspeaking,therearetwoissuesthatasuccessfulnumericaltechniquemustaddress.First,itneedstotrackatopologi-callycomplex,movingsolid-liquidinterfaceinbothtwoandthreespatialdimen-sions.Second,itmustbecomputationallyefficientastheseproblemsareusuallyparabolicinnaturewithstringentrestrictionsonthetimestepandsmallspatialscalesthatrequiresufficientgridresolution.Thus,adesirableschemeshoulduseimplicittimesteppingwithasymmetricinversionmatrixaswellashighorderaccuratespatialdiscretizationsofboththeparabolicpartialdifferentialequationandtheinterfaceitself.
Theinterfacethatseparatesthetwophasescanbetrackedeitherexplicitlyorimplicitly.Themaindisadvantageofanexplicitapproach,e.g.,fronttracking(see,e.g.,[17]),isthatspecialcareisneededfortopologicalchangessuchasmergingorbreaking.Whilethisiseasilyovercomeintwospatialdimensions,explicitlytreating
1
MathematicsDepartment&ComputerScienceDepartment,StanfordUniversity,Stanford,California94305.E-mail:fgibou@math.Stanford.edu2
ComputerScienceDepartment,StanfordUniversity,Stanford,California94305.3
MathematicsDepartment,UniversityofCalifornia,LosAngeles,California90095-1555.
183
0885-7474/03/1200-0183/0©2003PlenumPublishingCorporation
184Gibou,Fedkiw,Caflisch,andOsher
connectivityinthreespatialdimensionscanbedaunting.Implicitrepresentationssuchaslevelset[26]orphase-field[18]methodsrepresentthefrontasalevelsetofacontinuousfunction.Topologicalchangesareconsequentlyhandledinastraightforwardfashion,andthusthemethodsarereadilyimplementedinbothtwoandthreespatialdimensions.Moreover,onecaneasilymodeladditionalphysics,e.g.,materialstrainorflowpastdendrites.SometimesEulerianmethods,suchasthelevelsetmethod,arecriticizedfornotaccuratelypreservingthemassofamaterial.However,thisartifacthasrecentlybeenremovedforlevelsetmethodswiththeaidofmasslessmarkerparticlesthatobtaintheaccuracybenefitsofafronttrackingmethodwithouttheaddedhindranceofaddressingconnectivity,see[8].Moreover,in[9],theparticlelevelsetmethoddevelopedin[8]wasusedtotracktopologicallycomplexair/waterinterfacessubjecttoavarietyofpinchingandmerging.Theseaccuracylimitationshavenotyetbeenaddressedforphase-fieldmethods,butweareoptimisticthatthenonphysicalmasslosspresentinphase-fieldmethodscanbealleviatedtoalargedegreeusingamethodsimilartothatproposedin[8].
AsimplelevelsetapproachtosolvingthesharpinterfaceproblemdescribedinSec.2.2laterwasfirstproposedin[5].Theyusedthelevelsetmethodtokeeptrackofthefrontandsolvedforthediffusionfieldusinganimplicittimediscre-tizationmethod.Inordertoapplythisimplicittimediscretizationaconstantcoef-ficientmatrixneedstobeinvertedateverytimestep.Theirmatrixwasnonsymme-tricandtheyusedaratherslowGauss–Seideliterativeschemetoinvertlimitingboththegridresolutionandthenumberofspatialdimensions,i.e.,theywerenotabletoaddressStefanproblemsinthreespatialdimensions.In[21],theauthorsimproveduponthealgorithmpresentedin[5],forexamplecomputingthevelocityinamoreaccuratemanner.Theynumericallysimulatedthestandardfour-foldani-sotropytestcase,andobtainedresultsinexcellentagreementwiththepredictionsofmicroscopicsolvabilitytheory.
Inboth[5]and[21]thediscretizationofthetemperatureneartheinterfaceproducesanon-symmetricmatrixthatneedstobeinvertedforimplicittimestep-ping.Thelackofsymmetrymakesthisapproachcomputationallyexpensive,althoughmethodslikeGMRES[29]andBICGSTAB[29]mighthelptoalleviatethesituation.In[10],asymmetricsecondorderaccuratediscretizationforthePoissonequationwasoriginallydevelopedandlaterpresentedandshowntobesecondorderaccuratein[12].Thisalgorithmwasinspiredbytheghostfluidmethod[11]andhasbeensuccessfullyusedbyavarietyofauthors(e.g.,[6]).ApplyingthisdiscretizationtechniquetothetemperaturefieldneartheinterfaceallowsonetousearobustandefficientPreconditionedConjugateGradient(PCG)[13]methodtoinverttheconstantcoefficientmatrixresultingfromtheimplicitdiscretizationintime(thisalgorithmistobecontrastedwith[22]wheretheauthorsobtainedonlyfirstorderaccuracyinthepresenceofajumpcondition.HeresecondorderaccuracyisobtainedfortheDirichletboundarycondition).In[12],numericalresultsshowedthatthisschemeissecondorderaccurateforthevariablecoefficientandconstantcoefficientPoissonequationandtheheatequa-tion.Inparticular,weshowedthatthisnewalgorithmconvergestosomeknownexactsolutions,e.g.,theFrank-Sphere.Inthispaper,weapplythisalgorithmtothe
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth185
modifiedStefanproblemtakingintoaccountcrystallineanisotropy,surfacetensionandmolecularkinetics.
Themaindifferencebetweenthephase-fieldandlevelsetapproachisthatthelevelsetmethodcanbeusedtoexactlylocatetheinterfaceinordertoapplydiscre-tizationsthatdependontheexactinterfacelocation.Incontrast,thephase-fieldmethodonlyhasanapproximaterepresentationofthefrontlocationandthusthediscretizationofthediffusionfieldislessaccuratenearthefrontresemblinganenthalpymethod[7].Formulatingaphase-fieldmodelrequiresanasymptoticexpansionanalysisbeperformedwithasmallparameterproportionaltotheinter-facewidth,W.ItisimportanttonotethatthegridsizeisproportionaltoWandonlyinthelimitasWQ0doesthephase-fieldmethodconvergetothesharpinter-facemodel.Inthatsense,thephase-fieldmethodisonlyafirstorderaccurateapproximationtothetruemacroscopicsharpinterfacemodel.Thatis,evenifthenumericsaresecondorderaccurateforagivenvalueofW,themodelisinerrorbyW’O(gx)sothatthemethodcanbenobetterthanfirstorderaccurateoverall.Infactin[20]itwasshownrigorouslythatifthegridsizeisnotproportionaltoW,thenumericalresultsaregenerallyincorrect.Thelevelsetmethoddoesnotneedthisextralevelofadaptivity.Theinterestedreaderisreferredto[18]andthereferencesthereinformoredetailsonthephase-fieldmethod.2.EQUATIONSANDNUMERICALMETHOD2.1.LevelSetEquationandNumerics
Thelevelsetequation
F·Nf=0,ft+W
(1)
Fisthevelocityfield,isusedtokeeptrackofwherefisthelevelsetfunctionandW
theinterfacelocationasthesetofpointswheref=0.Theunreactedandreactedmaterialsarethendesignatedbythepointswheref>0andf[0respectively.Tokeepthevaluesoffclosetothoseofasigneddistancefunction,i.e.,|Nf|=1,thereinitializationequationintroducedin[30]
fy+S(fo)(|Nf|−1)=0
(2)
isiteratedforafewstepsinfictitioustime,y.HereS(fo)isasmoothedoutsignfunction.ThelevelsetfunctionisusedtocomputethenormalnF=Nf/|Nf|andthemeancurvatureo=N·nFinastandardfashion.ThelevelsetadvectionequationandthereinitializationequationarediscretizedbyusingtheHJ-WENOtypeschemes[15],seealso[23,16].Formoredetailsonthelevelsetmethodsee,e.g.,[25,24].2.2.Sharp-InterfaceModel
Dendriticsolidificationthatincludeseffectsofundercooling,surfacetension,crystallineanisotropyandmolecularkineticscanbedescribedbythesharp-interfacemodel.ConsideraCartesiancomputationaldomain,W,withexteriorboundary,“W,
186Gibou,Fedkiw,Caflisch,andOsher
andalowerdimensionalinterface,C,thatdividesthecomputationaldomainintodisjointpieces,W−andW+.Thesharp-interfacemodelisgivenby
“T
=N·(nNT)“t
xF¥W,
xF¥C,
(3)(4)
Vn=[nNT·nF]=(nNT·nF)r−(nNT·nF)u
F·nwhereTdenotesthetemperature,Vn=VFthenormalvelocityattheinterface,and
thesubscriptsuandrdefinetheunreactedandreactedmaterialsrespectively(forexample,inthecaseofasolidificationprocess,thereactedmaterialwouldbethesolidandtheunreactedonewouldbethemeltbath).Thethermalconductivityn(xF)
−+
isassumedcontinuousoneachdisjointsubdomain,WandW,butmaybedis-continuousacrosstheinterfaceC.Furthermore,n(xF)isassumedtobepositiveandboundedbelowbysomeE>0.Ontheboundaryofthecomputationaldomain,“W,weconsidereitherDirichletboundaryconditionsofT(xF)=g(xF)orNeumannboundaryconditionsofNT·nF(xF)=h(xF),althoughonecouldalsoconsiderperiodicboundaryconditionsinastraightforwardway.Attheinterface,weimposetheGibbs–Thomsonrelation
TI=−Eco−EvVn,
(5)
whereoisthecurvatureofthefront,Ecthesurfacetensioncoefficient,EvthemolecularkineticcoefficientandVntheinterfacevelocity.ThisinterfaceconditionaccountsforthedeviationoftheinterfacetemperatureTIfromequilibrium.Inthecasewherethesurfacetensionattheinterfaceisanisotropic,onecantakeforexampleintwospatialdimensions
Ec=d0(1−15Ecos(4a)),
(6)
whereEistheanisotropystrength,andaistheanglebetweenthenormalattheinterfaceandthex-axis.Thisformularepresentsthestandardfour-foldanisotropy.TheGibbs–Thomsonrelationiscomputedateverygridpointneighboringtheinterfaceandthenlinearlyinterpolatedtothefrontusinganeasytoimplementlevelsetprocedure,see[22].WhencomputingtheGibbs–Thomsonrelation,weusethevalueofthenormalvelocityVnattimetn.
Inonespatialdimension,thediscretizationofthetemperatureisasfollows.Forgridpointsmorethanonegridcellawayfromtheinterface,weuseastandardbackwardEulerimplicittimediscretization
ni+1
2−TnTn+1ii
=gt
1T
n+1i+1
−Tn+1Tn+1−Tn+1iii−1
−ni−1
2gxgxgx
212
(7)
whichissecondorderaccurateinspacewheretheresolutionoftheinterfaceiscrucial,butonlyfirstorderaccurateintime.Moreover,secondorderaccuratespatialresolutionisdesirablesinceEq.(4)computestheinterfacevelocityusingthefirstderivativesofthetemperature.
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth187
Thisdiscretizationisnotvalidiftheinterface,C,cutsthroughthestencilofpoints:xi−1,xiandxi+1.Forexample,supposethattheinterfacelocation,xI,fallsinbetweenxiandxi+1.Thenwhendiscretizingatxi,wedonothaveavalidvalueforTi+1atxi+1sinceTwillgenerallyhaveakinkacrosstheinterface(butmayalsobediscontinuousinsomemodels).WecircumventthisdifficultybydefiningaghostvalueofTGi+1atxi+1rewritingEq.(7)as
ni+1
2−TnTn+1ii
=gt
1T
Gi+1
−Tn+1Tn+1−Tn+1iii−1
−ni−1
2gxgx
.
gx
212
(8)
WithDirichletboundaryconditionsontheinterfaceC,wecancomputetheinterfacetemperatureasTIanduseittodefinetheghostvalueTGi+1.Then,possible
G
candidatesforTi+1include
TGI,i+1=T
TI+(h−1)TiTGi+1=h
and
T
G
i+1
(9)(10)
2TI+(2h2−2)Ti+(−h2+1)Ti−1=
h2+h
(11)
correspondingtotheuseofconstant,linearandquadraticextrapolationacrosstheinterface,respectively.Hereh=(xI−xi)/gx.PluggingEq.(11)intoEq.(8)givesthestandardsecondorderaccuratenonsymmetricdiscretizationusedbymanyauthors,e.g.,[5]and[21].Aspointedoutin[12],thislocalquadraticextrapola-tionisunnecessary,howevermanyauthorsstilluseitbeingmisledbylocalTaylorexpansionanalysis.[12]showedthatlocallinearextrapolationusingEq.(10)isenoughtoobtainsecondorderspatialaccuracyinthecasewheretheexactinterfacelocationisknown.Moreover,theresultingdiscretizationissymmetricenablingtheuseoffastlinearsolverssuchasPCGwhenemployingimplicitdiscretizationintime.PluggingEq.(10)intoEq.(8)givesasecondorderaccuratesymmetricdiscre-tizationof
ni+1
2T−T=gt
n+1i
ni
−T2T−T21T−n1hgxgx
I
i
i−1
2i
i−1
gx
.(12)
ThisistheschemethatwewilluseintheexamplesofSec.3.NotethatsincethetemperatureiscomputedtosecondorderaccuracyandVn=[NT·nF],theinterfacelocationisonlyfirstorderaccurate.Thereforeinthecaseofamovingfronttheaccuracyofthemethodfirstorder.TheCrank–Nicholsonschemewouldlowerthetruncationerrorbuttheschemewouldstillbefirstorderaccurateoverallduetothecalculationofthevelocityfield.
188Gibou,Fedkiw,Caflisch,andOsher
Incertainsituations,nmayonlybeknownatthegridnodesandtheinterfaceinwhichcaseni+1inEq.(8)isdeterminedfromaghostvalue,nGi+1,andtheusual2averaging
n
i+1
2
ni+nGi+1=2
(13)
notingthattheghostvalueiseasilydefinedusinglinearextrapolation
nI+(h−1)ni
nGi+1=h
(14)
accordingtoEq.(10).
Inmultiplespatialdimensions,theequationsarediscretizedinadimensionbydimensionmannerusingtheonedimensionaldiscretizationoutlinedabove.Theinterestedreaderisreferredto[12]forthedetailsofthismethod.
Wefindtheinterfacenormalvelocityusingthejumpinthetemperaturegra-dientinthenormaldirectionacrosstheinterface.Hereagain,weusethevalueofthetemperatureattheinterface,TI,andassumethatthetemperatureprofileislocallylineartocomputethetemperaturegradientoneachsideoftheinterface.Moreprecisely,wecomputethederivativesTxandTyofthetemperatureinthexandydirectionandthencomputethetemperaturegradientinthenormaldirectionTn=NT·nF.OnceTnisdefinedatgridpointsadjacenttotheinterface,weextrapo-latethevaluesof(Tn)rfromthereactedsideoftheinterfacetotheunreactedsideandextrapolatethevaluesof(Tn)ufromtheunreactedsidetothereactedsidesothatboth(Tn)rand(Tn)uaredefinedateverygridpointinabandabouttheinter-face.Thisisaccomplishedwithconstantextrapolationinthenormaldirectiontotheinterfaceaccordingto
Iy±n·NI=0
(15)
whereIisthevariabletobeextrapolatedandtheequationissolvedinfictitioustimeytosteadystate.Thiswasfirstimplementedin[5]usinganequationthatappearsin[31].ThisisdoneseparatelytoadvectI=(Tn)rinonedirectionandtoadvectI=(Tn)uintheotherdirection.NotethatwedonotsolveEq.(15)directly,butinsteaduseanoptimalspatialmarchingprocedure[1]inordertoobtainIlocaltotheinterface.See[12]formoredetails.Onceboth(Tn)rand(Tn)uarebothdefinedateachgridnodeneartheinterface,wecomputethejumpinanodebynodefashionusingthenodalvaluesof(Tn)rand(Tn)u,i.e.,(Vn)i=(nTn)r,i−(nTn)u,i.
Thetimesteprestrictionthroughoutthispaperisgivenby
gt=min(gtH,gtL),
with
gtH=0.5min(gx,gy,gz)
(17)(16)
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth189
chosenforaccuracyconsiderationsintheheatequation,
gtL
ww2w1g++[0.5,xgygz
1
2
3
(18)
F=(w1,w2,w3)=Vn·nwithWFchosenforstabilityofthelevelsetequation.3.NUMERICALRESULTS
In[12]weaddressedthesimulationsofameltwithconstanttemperatureattheinterfaceandcomparedoursolutionstotheFranksphereexactsolution.Here,weshowthatourmethodcantakeintoaccountthedifferentphysicalparametersneededtofaithfullymodelcrystalgrowth.3.1.TopologicalChanges—MergingandBreaking
Thesetwoexamplesillustrateoneofthemainadvantagesofthelevel-setmethod,whichisthatithandlestopologicalchangesandcomplexinterfacialshapesinanaturalway.Figure1showssevenseedsgrowingunderaconstantnormal
Fig.1.Severalseedsofarbitraryshapesgrowunderthenormalvelocityvn=1andmergetogether.
Thisexamplewasranwith100gridpointsineachdirection.
190Gibou,Fedkiw,Caflisch,andOsher
Fig.2.Abulkofmaterialisevolvedwithaninterfacialnormalvelocityvn=−1andeventually
breaksapart.Thisexamplewasranwith100gridpointsineachdirection.
velocityandmerging.ThischaracteristichasbeenexploitedintheIslands-Dynam-icsmodelforepitaxialgrowthtomodelthemergingofislands[4,6,28].Asanote,aconstantvelocityattheinterfaceisphysicallyrelevantsince,forexample,itcanmodelinfiniteedgediffusionattheinterfaceofanisland[4,6,28].Figure2showsthebreakingofcrystals.Inthisexample,asinglecrystalisshrunkwithanegativeconstantnormalvelocity.Thisfeaturehasbeenappliedtothemodelingofepitaxialgrowthtoproperlyaccountforthethermaldetachmentofatomsfromislandedges,whereislanddissociationcanoccur[27].Themergingandthebreakingofcrystalsaredoneautomaticallywithoutanyspecialtreatmentsasopposedtowhathastobedonewithexplicitinterfacerepresentations.3.2.EffectofVaryingIsotropicSurfaceTension
Figure3demonstratestheeffectofaddingisotropicsurfacetensionattheinterfacebyimposingT=−Eco,andthecorrespondingsmoothingeffectonthecrystal.Figure3adepictstheevolutionofapuremeltwithEc=0.Sincethiscaseis
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth191
Fig.3.Effectofvaryingisotropicsurfacetension.WeimposetheGibbs–ThomsonrelationattheinterfacewithEv=0and(a)Ec=0,(b)Ec=0.0005,(c)Ec=0.001.Gridsizesusedare300×300.
mathematicallyunstable,theregularizationbuilt-intothelevelsetmethodallowsitscalculation,see[14].Ourmethodseemstointroducelessnumericaldiffusionthanpreviousalgorithmsusingthelevelsetmethod,see[5]foracomparison.Thiscomesinpartfromahigherorderaccuracyinthelevelsetevolution.Thistranslatesintothemoredetaileddendritestructureobtainedwithourmethod.Figures3bandcshowthesmoothingeffectofincreasingthesurfacetensionattheinterface.3.3.GridRefinement
Thisexample,takenfrom[5],teststheconvergenceofourmethodundergridrefinement.ConsiderasmallfrozenseedofmaterialplacedinasurroundingregionofundercooledliquidwithtemperatureT.=−0.5.Thecomputationaldomainis[−2,2]×[−2,2]andtheinitialshapeisgivenintermsofthefollowingparametricequations,
x(s)=(R+Pcos(8ps))cos(2ps)y(s)=(R+Pcos(8ps))sin(2ps)
whereR=0.1,P=0.02ands¥[0,2p].WeapplytheGibbs–Thomsonrelation,Eq.(5),attheinterfacewithEc=0.002andEv=0.002.Forcoarsergrids,Fig.4a,thenumericaldiffusionovercomesthephysicalregularization.However,aswerefinethegrid,theartificialnumericaldissipationbecomesnegligiblecomparedtothephysicalsurfacetensionandthezerolevelsetoffconvergestoasimilarshape.Theconvergenceoftheseplotsundergridrefinementarecomparabletothoseobtainedby[5],butweachieveconvergencewithfewergridpointsasourmethodhaslessartificialnumericaldissipation,duetoouruseofmoreaccuratenumericalmethodsinavarietyofplaces.
3.4.EffectofAnisotropicSurfaceTension
Physicalanisotropyforcesacrystaltogrowalongpreferreddirections.InthisexampleweimposetheconditionT=−Ecoattheinterfacewith
Ec=0.001(8/3sin4(2a−p/2)),
(19)
192Gibou,Fedkiw,Caflisch,andOsher
Fig.4.Growthhistoriesforfourgridresolutions:fromlefttorightandtoptobottom,120points,160points,200points,240pointsineachspatialdimension.WeapplytheGibbs–ThomsonrelationattheinterfacewithEc=0.002andEv=0.002.Aswerefinethegridtheartificialnumericaldissipa-tionbecomesnegligiblecomparedtothephysicalmechanismsandthezerolevelsetoffconvergestoasimilarshape.
whereaistheanglebetweenthenormaltotheinterfaceandthex-axis.Thisfour-foldanisotropywillcausethecrystaltogrowalongthefourdiagonalaxesasdetailedin[2].ConsidertheinitialseedtobeanirregularpentagononadomainW=[−1.5,1.5]×[−1.5,1.5]onagridwith200pointsineachspatialdimension.Weletthecrystalgrowtoafinaltimetfinal=0.4.Figure5demonstratesthatthecrystalindeedgrowsalongthepreferreddiagonaldirections.Notethatthedendritegrowingfromtheinitialsingularityatthetopcornerofthepentagonisproperlylimitedbytheanisotropicsurfacetensionthatattainsamaximumalongthey-axis.3.5.GridOrientationEffectswithAnisotropicSurfaceTension
Thistestdemonstratesthattheartificialgridanisotropyisnegligible.ConsideradomainW=[−1,1]×[−1,1]with100gridpointsineachspatialdimensionon
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth
1.5193
10.50–0.5–1–1.5–1.5–1–0.500.511.5Fig.5.Effectofanisotropicsurfacetension.TheinterfaceconditionisthefourfoldanisotropyboundaryconditionT=−EcowithEc=0.001(8/3sin4(2a−p/2)).Thecrystalgrowsalongthepre-ferreddiagonaldirectionsandthedendritegrowingfromtheinitialsingularityatthetopcornerofthepentagonisproperlylimitedbytheanisotropicsurfacetensionthatattainsamaximumalongthey-axis.
whichweseedacrystalwithinitialshapedescribedbythefollowingparametricequation:
x(s)=0.4cos(4s)cos(s)y(s)=0.4cos(4s)sin(s)
wheres¥[0,2p].TheinterfaceconditionisthefourfoldanisotropyboundaryconditionT=−EcowithEc=0.001(8/3sin4(2(a−a0))),whereaisdefinedasinSec.3.4anda0isthepreferreddirectionofgrowth.Figure6adepictstheevolutionofthelevelsetfunctionuptoatimet=0.04whena0=0.Thepreferreddirectionsarethexandyaxeswherethedendritetipssharpen.Onthediagonalaxes,thesurfacetensionattheinterfaceattainsamaximumvaluethatforcesthedendritetipstowiden.Figure6billustratestheevolutionofthesameinitialdatawhena0=p/4.Thepreferreddirectionsarealongthediagonaldirectionsasexpected.Moreover,theshapeisthatofFig.6arotatedbyp/4demonstratingthatthearti-ficialgridanisotropyisnegligible.Inthisexample,wecomputethenormalvelocitycomponentinfourdifferentcoordinatesdirectionsasdetailin[5].3.6.EffectofVaryingtheThermalConductivity
Thisexampleillustratestheimpactofvaryingthethermalconductivity.The
F],sothestrongerthethermalvelocityattheinterfaceisgivenbyVn=[nNT·n
conductivityn,thefasterthegrowth.HereweconsideradomainW=[−0.5,0.5]×[−0.5,0.5]with100gridpointsineachdirection.Theinitialdataisgivenbythefollowingparametricequations:
194
11Gibou,Fedkiw,Caflisch,andOsher
0.80.80.60.60.40.40.20.200–0.2–0.2–0.4–0.4–0.6–0.6–0.8–0.8–1–1–0.8–0.6–0.4–0.200.20.40.60.81–1–1–0.8–0.6–0.4–0.200.20.40.60.81Fig.6.Gridorientationeffectswithanisotropicsurfacetensiononagridwith100pointsineachspatialdimension.TheinterfaceconditionisthefourfoldanisotropyboundaryconditionT=−0.001(8/3sin4(2(a−a0)))owith(left)a0=0and(right)a0=p/4.Theshapeofthecrystalintherightfigureisthatofthecrystalintheleftfigurerotatedbyp/4demonstratingthattheartificialgridanisotropyisnegligible.
x(s)=0.2(0.5+0.2sin(6s))cos(s)y(s)=0.2(0.5+0.2sin(6s))sin(s)
wheres¥[0,2p]andweletthecrystalgrowtoafinaltimeoft=0.025.InthisexampleweimposetheconditionT=−EcoattheinterfacewithEc=0.001.Figure7depictstheresultswithtwodifferentvaluesofthethermalconductivityn.
0.50.50.40.40.30.30.20.20.10.100–0.1–0.1–0.2–0.2–0.3–0.3–0.4–0.4–0.5–0.5–0.4–0.3–0.2–0.100.10.20.30.40.5–0.5–0.5–0.4–0.3–0.2–0.100.10.20.30.40.5Fig.7.Effectofvaryingthethermalconductivityonagridwith100gridpointsineachspatial
F]islargerforhighervalueofdimension.ThisdemonstratesthattheinterfacevelocityVn=[nNT·n
thethermalconductivity.n=0.5fortheleftfigureandn=1fortherightfigure.
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth195
3.7.DifferentDiffusionConstants
Oneadvantageofusingalevelsetformulationistheabilitytosimulatethegrowthofacrystalwithadifferentthermalconductivitythanthemediumthatitisin.Thisillustratestheversatilityofourlevelsetapproach.Thephase-fieldfor-mulation,ontheotherhand,requiresonetoperformasymptoticexpansiontoincludedifferentthermalconductivityoneachsideofthefront.Consequently,thishasbeenachallengeuptorecently[19].
WeconsiderthesameproblemasinSec.3.6butwithdifferentthermalcon-ductivitiesoneachsideoftheinterface.Figure8depictstheresultswhenthethermalconductivityinsidenr=1andthethermalconductivityoutsidenu=1.2.Notethatinourapproachthetwo-phaseproblemisseparatedintotwodistinctonesinceweimposeaDirichletconditionattheinterfaceofeachsub-domain.There-fore,onecouldtaketherationr/nuaslargeaswantedandthealgorithmwouldperformaswellincontrastwithphase-field.Wehavechosenamoresubtlediffer-enceintheratiotofittheresultsinthegraph.In[21],theauthorsalsostudieddiscontinuousthermalconductivitieswiththeirnonsymmetriclevelsetapproachandshowedthattheirnumericalresultsconvergedtothelinearizedsolvabilitytheoryofBarbieriandLanger[3].3.8.ComparisontoSolvabilityPredictions
Anoftenusedexampletotestthevalidityofthephase-fieldapproachistocomputethedendritetipvelocityandtipshapeofthestandardfourfoldanisotropycrystallinegrowth.ConsideradomainW=[−8,8]×[−8,8]withagridspacingDx=0.01.Weseedadiskofradius0.15inanundercooledbathwithtemperatureTundercool=−0.65.WeimposetheGibbs–Thomsonrelation5attheinterfacewith
2tip=Vtipd0/DreachesE=0.05,d0=0.01andEv=0.ThedimensionlesstipvelocityV
0.50.50.40.40.30.30.20.20.10.100–0.1–0.1–0.2–0.2–0.3–0.3–0.4–0.4–0.5–0.5–0.4–0.3–0.2–0.100.10.20.30.40.5–0.5–0.5–0.4–0.3–0.2–0.100.10.20.30.40.5Fig.8.Effectofhavingdifferentthermalconductivityoneachsideoftheinterfaceonagridwith100gridpointsineachspatialdimension.Thefinaltimeis0.025.Fromlefttoright:nr=1andnu=1,nr=1,andnu=1.2.
196
8Gibou,Fedkiw,Caflisch,andOsher
6420–2–4–6–8–8–6–4–202468Fig.9.Sequenceofdendriticshapesinthecaseofstandardfourfoldanisotropy.Thesolutioniscomputedona[−8,8]×[−8,8]domainwith400pointsineachdirections,theundercoolingis−0.65andtheanisotropystrengthisgivenviatheGibbs–ThomsonrelationTI=0.01(1−0.75cos4a),whereaistheanglebetweenthenormaltotheinterfaceandthex-axes.Thefinaltimeist=1.42.
0.1Level Set Solvability0.090.08Vtipd0/D0.070.060.050.0405000tD/d2
0
1000015000Fig.10.Timeevolutionofthedimensionlesstipvelocityinthecaseofstandardfourfoldaniso-tropy.Thesolutioniscomputedona[−8,8]×[−8,8]domainwith400pointsineachdirections,theundercoolingis−0.65andtheanisotropystrengthisgivenviatheGibbs–ThomsonrelationTI=0.01(1−0.75cos4a),whereaistheanglebetweenthenormaltotheinterfaceandthex-axes.
2tip=Vtipd0/DreachesasteadystatevalueThefinaltimeist=1.42.ThedimensionlesstipvelocityV
of0.047consistentwithsolvabilitytheory.
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth197
Fig.11.3DStefanproblem.Thecontoursshowtheevolutionoftheinterfacelocationfromtimet=0s(topleft)tot=0.14s(bottomright).WeapplytheGibbs–ThomsonrelationattheinterfacewithEc=0.002andEv=0.0.002.Thiscomputationuses100gridpointsineachspatialdirection.
198Gibou,Fedkiw,Caflisch,andOsher
asteadystatevalueof0.047consistentwithsolvabilitytheory(seeFig.10)andFig.4illustratesthedendriticshapesofthestandardfourfoldanisotropy.3.9.DendriticGrowthin3D
Thefollowingexampleillustratesourmethod’spotentialformodelingcrystalgrowthinthreespatialdimensions.Weobtainafairamountofdetailinthedendritestructureswithlittlecomputationaleffort,sinceourmethodyieldsasymmetriclinearsystemthatcanbeinvertedefficiently(numericalexperimentswereperformedonaPentiumIIIlaptop).
LetW=[−1.5,1.5]×[−1.5,1.5]×[−1.5,1.5].Theinitialdataisgivenbyf=min(f1,min(f2,min(f3,min(f4,min(f5,f6))))),wherethefi’sarespheresofradius0.1andcenteredoffthex-axisory-axisorz-axisby±0.05.Initially,T=0insideW−,andT=−0.5outsideW−.DirichletboundaryconditionsofT=−0.5areenforcedon“W.WeapplytheGibbs–Thomsonrelation,Eq.(5),attheinterfacewithEc=0.002andEv=0.002.Theevolutionofthisinitially‘‘bumpy’’sphereispresentedinFig.11.Thesnapshotsaregivenevery0.014sfromtheinitialtimet=0stothefinaltimet=0.14.4.CONCLUSION
Inconclusion,thelevelsetmethodshouldbeconsideredasamethodofchoiceforthestudyofdendriticcrystallization.Wehaveshownthatouralgorithmcanproduceaccuratesolutionsthatcanbecomparedfavorablywithsolvabilitytheory.Inaddition,ourspatialdiscretizationofthetemperatureyieldssymmetricmatricesfortheimplicittimediscretizationthatcanbeinvertedwithfastlinearsolvers,makingitveryefficient.ACKNOWLEDGMENTS
F.G.wassupportedinpartbyanNSFpostdoctoralfellowship#DMS-0102029,aFocusedResearchGroupGrant#DMS0074152fromtheNSFandbyONRN00014-01-1-0620.R.F.wassupportedinpartbyanONRYIPandPECASEawardN00014-01-1-0620andNSF-0106694.R.C.wassupportedinpartbyaFocusedResearchGroupGrant#DMS0074152fromtheNSF.S.O.wassupportedinpartbyONRN00014-97-0027andNSFDMS-0074735.REFERENCES
1.Adalsteinsson,D.,andSethian,J.(1999).Thefastconstructionofextensionvelocitiesinlevelsetmethods.J.Comput.Phys.148,2–22.
2.Almgren,R.(1993).Variationalalgorithmsandpatternformationindendriticsolidification.J.Comput.Phys.106,337.
3.Barbieri,A.,andLanger,J.(1989).Predictionsofdendriticgrowthratesinthelinearizedsol-vabilitytheory.Phys.Rev.A.39,5314.
4.Caflisch,R.,Gyure,M.,Merriman,B.,Osher,S.,Ratsch,C.,Vvedensky,D.,andZinck,J.(1999).Islanddynamicsandthelevelsetmethodforepitaxialgrowth.Appl.Math.Lett.12,13.
ALevelSetApproachfortheNumericalSimulationofDendriticGrowth199
5.Chen,S.,Merriman,B.,Osher,S.,andSmereka,P.(1997).AsimplelevelsetmethodforsolvingStefanproblems.J.Comput.Phys.135,8.
6.Chen,S.,Merriman,B.,Kang,M.,Caflisch,R.,Ratsch,C.,Cheng,L.-T.,Gyure,M.,FedkiwR.,Anderson,C.,andOsher,S.(2001).Levelsetmethodforthinfilmepitaxialgrowth.J.Comput.Phys.167,475–500.
7.Chorin,A.(1985).Curvatureandsolidification.J.Comput.Phys.58,472.
8.Enright,D.,Fedkiw,R.,Ferziger,J.,andMitchell,I.(inpress).Ahybridparticlelevelsetmethodforimprovedinterfacecapturing.J.Comput.Phys.
9.Enright,D.,Marschner,S.,andFedkiw,R.(2002).Animationandrenderingofcomplexwatersurfaces,SIGGRAPH2002.ACMTrans.onGraphics21,736–744.
10.Fedkiw.R.ASymmetricSpatialDiscretizationforImplicitTimeDiscretizationofStefanType
Problems,unpublished,June1998.
11.Fedkiw,R.,Aslam,T.,Merriman,B.,andOsher,S.(1999).Anon-oscillatoryEulerian
approachtointerfacesinmultimaterialflows(theghostfluidmethod).J.Comput.Phys.152,457.
12.Gibou,F.,Fedkiw,R.,Cheng,L.-T.,andKang,M.(2002).Asecondorderaccuratesymmetric
discretizationofthepoissonequationonirregulardomains.J.Comput.Phys.176,205.
13.Golub,G.,andVanLoan,C.(1989).MatrixComputations,TheJohnsHopkinsUniversity
Press,Baltimore.
14.Harabetian,E.,andOsher,S.(1998).Regularizationofill-posedproblemsviathelevelset
approach.SIAMJ.Appl.Math.58,1689.
15.Jiang,G.-S.,andPeng,D.(2000).WeightedENOschemesforHamilton–Jacobiequations.
SIAMJ.Sci.Comput.21,2126–2143.
16.Jiang,G.-S.,andShu,C.-W.(1996).EfficientimplementationofweightedENOschemes.
J.Comput.Phys.126,202–228.
17.Juric,D.,andTryggvason,G.(1996).Afronttrackingmethodfordendriticsolidification.
J.Comput.Phys.123,127.
18.Karma,A.,andRappel,W.-J.(1997).Quantitativephase-fieldmodelingofdendriticgrowthin
twoandthreedimensions.Phys.Rev.E.57,4323.
19.Karma,A.(2001).Phase-fieldformulationforquantitativemodelingofalloysolidification.
Phys.Rev.Lett.87,11.
20.Merriman,B.,Bence,J.,andOsherS.(1994).Motionofmultiplejunctions,alevelset
approach.J.Comput.Phys.112,334.
21.Kim,Y.-T.,Goldenfeld,N.,andDantzig,J.(2000).Computationofdendriticmicrostructures
usingalevelsetmethod.Phys.Rev.E.62,2471.
22.Liu,X-D.,Fedkiw,R.,andKang,M.-J.(2000).Aboundaryconditioncapturingmethodfor
Poisson’sequationonirregulardomain.J.Comput.Phys.160,151–178.
23.Liu,X-D.,Osher,S.,andChan,T.(1996).Weightedessentiallynon-oscillatoryschemes.
J.Comput.Phys.126,200–212.
24.Osher,S.,andFedkiwR.(2001).Levelsetmethods:Anoverviewandsomerecentresults.
J.Comput.Phys.169,475.
25.Osher,S.,andFedkiw,R.(2002).LevelSetMethodsandDynamicImplicitSurfaces,Springer-Verlag,NewYork.
26.Osher,S.,andSethian,J.(1988).Frontpropagatingwithcurvature-dependentspeed:Algo-rithmsbasedonHamilton–Jacobiformulation.J.Comput.Phys.79,12.
27.Petersen,M.,Ratsch,C.,Caflisch,R.,andZangwill,A.(2001).Levelsetapproachtoreversible
epitaxialgrowth.Phys.Rev.E.64,061602.
28.Ratsch,C.,Gyure,M.,Caflisch,R.,Gibou,F.,Petersen,M.,Kang,M.,Garcia,J.,and
Vvedensky,D.(2002).Level-setmethodforislanddynamicsinepitaxialgrowth.Phys.Rev.B.65,195403.
29.Saad,Y.(1996).IterativeMethodsforSparseLinearSystems,PWSpublishingcompany,
Boston.
30.Sussman,M.,Smereka,P.,andOsher,S.(1994).Alevelsetapproachforcomputingsolutions
toincompressibletwo-phaseflow.J.Comput.Phys.114,146.
31.Zhao,H.-K.,Chan,T.,Merriman,B.,andOsher,S.(1996).Avariationallevelsetapproachto
multiphasemotion.J.ComputPhys.127,179.
因篇幅问题不能全部显示，请点此查看更多更全内容