Numerical Simulation of Dendrite Evolution during Solidification Process

Vol.12 No. 1J. Iron & Steel Res. , Int.Jan, 2005Numerical Simulation of Dendrite Evolution during Solidification ProcessLI Qiang'，GUO Qiao-yi?，REN Chuan-fu3(1. Shenyang University of Technology, Shenyang 110023, China;2. Shenyang Ligong University,Shenyang 110168, China; 3. Shenyang University, Shenyang 110044, China)Abstract: In order to precisely describe the dendrite evolution during solidification process, especially in micro-scale, a continuous method is presented to deal with the discontinuous physical properties beside the solid/liquidinterface. In this method, the physical properties are used as averaging physical properties of solid phase and liquidphase in the interface zone, which can smooth the property gap between solid and liquid phases, and make theproperties from liquid phase to solid phase. The simulated results show that the method can represent the side-branches and the solute micro-segregation well.Key words: numerical simulation; continuous method; dendrite evolution; micro segregationSolidification is a very complicated process indifferent values. So some researchers used othercluding heat and solute transfer， phase transitionmethods to deal with the discontinuous properties,and interface migration processes which directly in-Damir JuricCL3] adopted the Front-Tracking methodfluence solidification microstructure. Solidification(FT) using complex expression to define the posi-microstructure is an important index to evaluate me-tion of solid/ liquid interface, which limits the calcu-chanical and other physical properties of the finallating scale. The sharp -interface method (SI), usedproduct. In the past, most of researchers used ex-by Udaykumar H S14, has too much calculation ofperimental and theoretical methods to study the so-the interface temperature by iterating the solute con-lidification theory[1,2].However, experimentalcentration. The phase field (PF) method requiresmethods are expensive and cannot avoid the limita-the mesh size to be smaller than the interface thick-tion caused by the conditions, while the theoreticalnessl1s,16], which also limits the calculation scale.methods cannot describe the solidification micro-As stated by Christof Eck et al[1]，these models re-structure preciselyt3. The numerical simulationquire much more details of dendrite in order to ob-method can overcome the problems and predict the tain reasonable results, which are only suitable tosolidification microstructure accurately, so there ismicro-scale. In this paper, Cellular Automatonmuch interest on this method in recent years4- 17] .(CA) method is adopted because this method is notDendrite is one of the most common morpholo-limited by simulated size and CA rules are simple togies of solidification microstructure and has a directunderstand.impact on the final casting mechanical propertiesConsidering the discontinuous properties in in-through the local cooling conditions. Therefore, terface zone, a continuous method is developed, assimulation of the dendrite formation and evolutionillustrated in Fig. 1. In this method, the propertieshas attracted much attention[4- 17].of interface can be treated as a continuous functionWhen liquid transforms into solid phase, heat of solid fraction to smooth property gap between theconductivity, specific heat and density of the liquid solid and liquid phases, so that the unstable heat andmaterial change greatly. In order. to solve these dis-solute transfers can be reflected to precisely describecontinuous problems, some researchers took the sol-the solid fraction, interface velocity and solute con-id and liquid physical properties as a constant and re-centration redistribution in interface zone. Mean-garded their values as equalC4- 12]. However, this apwhile, the averaging solute concentration is adoptedproximate simulation make the results deviate from to keep the total solute concentration conservation.real dendrite morphology because the heat transferand solute diffusing rates in different phases have of th中国煤化工。-only requires its ve-Bography: LI Qiang(1973-), Male Doctor, Letureship; E mail: qli@imr. acTHCNM H G,No.1Numerical Simulation of Dendrite Evolution during Solidification ProcessIn the interface cell, solute conservation equa-Volume elementLiquid.tion is given[14-6] as:physical,2(0propertyLiquidphysica(1-k,)C*Vx=D.(N).-D(),(3)l propertyIn addition, the chemical potential is assumedSold Liquidas equilibrium beside the solid/liquid interface,x_(0x(0)which causes the redistribution of solute concentra-tion, as follows:hak.f+kxhC,=kpc*(4)f+f=1Boundary conditions:SoldpaC=0aT_(5)an)第=q q=h(T" - Tm)Liquid : Iteiface i0 f(x)→5,T(1=0,x,y)= To(x,y), C(=0,x,y)=Co(x,y)Initial conditions: .(6)Fig. 1 Continuous method to deal withdiscontinuous propertieswhere ρ is density; Cp is specific heat; k is thermalconductivity; Ln is latent heat released by phaselocity and which cell it' s located. Therefore, there istransition from liquid to solid phase; f, is solid frac-no much calculation for tracking and can be extendedtion; T is temperature, t is time; C is solute concen-to macro-scale simulation. Moreover, the stochastictration; D is solute diffusion coefficient, and in solidmethod is used to show the nucleation and growth,and liquid cell, D is equal to D and D, respective-thus the evolution can reflect the random feature ofy; C* is solute concentration in interface cell; Vx issolidification microstructure.normal velocity of solid/ liquid interface; kp is soluteContinuous Method of Heat and Speciespartition coefficient in interface cell, which is as-sumed as a constant and reflects solutal concentra-Transfertion gap between solid and liquid phase beside inter-Dendrite evolution process involves two primaryface; n is normal direction of solid/ liquid interface;aspects: heat and solute transfer model, and nuclea-Tm is temperature of boundary cell; h is heat trans-tion and interface migrating process in undercoolingfer rate on boundary; Tsur is surrounding tempera-liquid, and these two processes always interact andture, assumed as a constant Tsr = 298 K; subscriptcorrelate with each other. To calculate precisely,s and 1 stand for the solid and liquid phase respec-the model must consider the interactional relationtively.between heat and solute transfer, and microstruc-Eqn. (3) can be discretized, and Eqn. (7) andture.Eqn. (8) are obtained[4-6]:The model is based on the following assump-__D1-C(i-1.j).)f(i-1,j)+tions:①The density of material keeps constantV.- a(1-k,t(-C(i,j)during solidification; ②There are only three possi .1- CL(i+1,j)、f(i+1,)]+-kpD.ble phase states, i. e. liquid phase, solid phase andC(i,) )a(1-kp)interface，which correspond to the liquid cell, solidC,(i-1,j))f,(i-1,j)+cell and interface cell respectively; ③Solute cannot[(1-kp.C*(i,j)，transfer across the solid/liquid interface. At the be-C,(i+1,j)ginning of solidification, the initial temperature fieldk。.C(i,j))f。(i+1,i)](7)and solute concentration of liquid is assumed to be u-v,=-D1-[(1-C(i,-)f(i,j-1)+niform.a(1-kpC'(i,j)Energy conservation equation is expressed as:1- C(i,j+1))(i,;+1)]+a&-D,djC"(i,j)pC, i=v●(kVT)+pLo(1)dt中国煤化工1>+Solute concentration conservation equation isYHCNMH Ggiven as:k。●C"(i,j)，f,(i,j+1)]<8)aC3t= V (DVC)(2)where fi and f, are liquid fraction and solid fraction30●J. Iron & Steel Res.，Int.Vol.12respectively, and fi+f.=1.+η* sin0，and η is amplitude of the noise and equalsBased on KGT model, the solute concentrationto0.1, as in Ref.[4] and Ref. [5]. .in interface cell is given as followsI1,3.6-6;The anisotropy of solid/liquid interface can beT- T+rkf(0,p)given as follows'l] :C°=C+m(9)f(0,p)=1+δ. cos[4(φ- -0)](14)where m is liquidus slope; k is curvature of solid/liq-To keep the above equations convergent, timeuid interface; r is Gibbs Thomson coefficient; Jstep is very important to save the calculating time.(q,日) is anisotropic function related to preferentialIn this model, the time step is selected as a variablegrowth angle (θ) and normal growth angle (p); Coautomatically, which depends on the meshed sizeis initial liquid solute concentration; To is meltingand physical properties. In solid/liquid interface intemperature; and T can be obtained from Eqn. (1).one time step, the interface can' t move over to aThe solid/liquid interface curvature is used tomesh cell., For the same reason, the time step mustcalculate each interface cell curvature, as fol-make Eqn. (1)- (3)，(5)-(7) stable. In this pa-lows'4-8]:per, the time step is given as follows-4-8] :f.+ gf。(k)pdi= min[(15)的部[1-2:2- N+I 」(10)言(0.号V2)]where Vmax is maximum velocity at that step; a iswhere N is total number of neighbor cells, and threemesh size in simulated zone; min is a function to se-nearest layer of cell are used to calculate the averagecurvature; f,(k) is solid fraction of neighbor celllect a minimum value from .，.. and . Thek"; a is length of a cell. .coefficient 1/5 means that the interface cannot moveIn both sides of the solid/liquid interface cell,over to 1/5 length of the CA cell. .solute concentration is higher than the neighbor cell,2 Dendrite Growth Modelso solute will diffuse into solid phase and liquidphase. In each finite difference element, the soluteTo simplify the process, a solid seed is input in-concentration distribution is not continuous but dis-to the center of the ingot. At a certain cooling rate,crete. Therefore, in these zones, the average meth-the solid seed will randomly capture a neighbor liq-od is adopted to deal with the solute concentration.uid cell and become the interface cell. A neighborAssuming that there are three kinds of solute con-cell may be captured when it is in liquid phase andcentration: solid, liquid and averaging solute con-its temperature is below liquidus. The capturingcentrations, the averaging solute concentration C isprocess is only in the four nearest neighbors. Thegiven as follows:capturing possibility is given as follows:C=f.C.+f{CPc= cosq(16)where C. is solute concentration in solid-side of inwhere Pc is capturing possibility.terface cell; and C is solute concentration in liquid-For every liquid neighbor cells, when the tem-side of interface cell, and they can be given by:perature is below the liquidus, a random number R。C,=. k,C"(0≤R,≤1) will be generated. If R,> Pc, then theC=(C- f,C.)/filiquid cell will be captured, or it will keep the origi-where C" can be obtained from Eqn. (8).nal state. When captured, it will become an inter-The solute concentration can be calculated byface cell and its growth direction will be consistentconsidering the solute concentration flux in eachwith the solid phase.time step-1When liquid temperature is above the liquidus,In the interface cell，the solid fraction can bethere is only the liquid phase in whole simulatedcalculated4.5]:zone and the liquid phase state is equal to zero. Asthe liquid cell temperature is cooled below the liqui-8f.=(V.+V,-V.v,当)台f:=f:'+β.8f。dus,中国煤化工quid to be captured.(13)Once.YHCNMHG the liquid cell willwhere δf, is increase of solid fraction in interfacebe translormed Into intertace cell and its phase statecell; fi: is solid fraction at n"h time step; f-1 is solidwill be equal to phase state of its parent cellfraction at (n- 1) time step; dt is time step; β=1Once the liquid phase is completely solidified, itNo.1Numerical Simulation of Dendrite Evolution during Solidification Process●31●will grow in a certain direction. The dendriteTable 1 Simulated parameter selection in modelgrowth direction reflects the dendrite growth anisot-Density/(kg. m-3)7 300ropy. In this paper, the growth angle is the normalThermal conductivity of solid phase/(W. m-?●K-1)33velocity angle of solid/liquid interface, and the nor-Thermal conductivity of liquid phase/(W.m-2●K-1)8mal velocity angle is:Specic heat of solid phase/(kJ.kg-' .K~1)650Specific heat of liquid phase/(kJ.kg-1.K-1)80φ= arctan((17)Solute dffusiviy in solid phase/(m2●s-1)5X 10-10(削)Solute dffusivity in liquid phase/(m2●s 1)2X10-9where the preferential growth direction is a randomLatent heat/(kJ. kg-1)284function whose value is between zero and one. V,Liquid slope- 82and V, can be obtained from Eqn. (7). At the begin-Thomson-Gibbs cefficient1.9X10-7ning, Vx and V, is equal to zero, here φ=0. .Partition cofficient0.6The calculation step is as follows:①Based onMelting temperature/K1763the initial condition and boundary condition to deter-mine the time step by Eqn. (14);②Calculating thetemperature field according to Eqn. (1);③Determi-ning the state of each cell based on the local condi-tion by CA transition rules; ④Calculating the ki-0.2mmnetics of solid/liquid interface based on Eqn. (3),Eqn. (4) and Eqn. (7)- Eqn. (9), including the sol-(White color stands for non-liquid zone,id/liquid interface velocity， solute concentration,and black color stands for liquid zone)solid fraction and curvature;⑤Calculating the sol-Fig.2 Phase state evolution during solidificationute concentration by Eqn. (2);⑥Selecting the ini-tial condition as the final data;⑦Repeating the a- then into complete solid cell slowly. In Fig.2 (c),bove steps until the total calculating steps reach the the dendritic contour continues growing along x andpreset number of steps.y axes since the temperature decreases continuouslyat 8 000 steps. Fig.2 (d) represents that the liquid3 Result and Discussionphase only exists in the four corners of simulatedIn this paper, a 2D square is used to simulatezone at 12 000 steps, and the velocity of interfacethe dendrite evolution. The simulated zone is select-and the solid capturing possibility are nearly theed as square shape with a size of 1 mmX1 mm, insame between the two axes. From Fig.2, it can bewhich, there are 10 000 cells with size of 0. 01 mmXseen that the solid/ liquid interface velocity and cap-0.01 mm. The cooling rate is 5000 C/s. The initialturing possibility in various directions are different.liquid temperatureis 1 490 C, and initial solute con-Fig.3 is the averaging solute concentration discentration is 0. 6%. The simulated material is a Fe tribution in the whole simulated zone, reflecting theC binary alloy and its physical properties are listed inmicro- segregation forming process, and the color barTable 1. The total calculating steps are 20 000.stands for the average solute concentration of eachFig. 2 shows phase states at 2000, 6 000，8 000cell at different steps.and 12 000 steps respectively. Fig. 2 (a) shows thatThat distribution can be divided into threethe non-liquid zone is nearly the same in different di-zones: solute concentration near to 0.4% in dendriterections at 2 000 steps, indicating that dendrite sizecorresponding to solid phase; solute concentrationis very small at the beginning and the capturing pos-near to 0. 6% corresponding to the liquid phase; andsibility is nearly the same in each direction. The ani-solute concentration over 0. 7% corresponding to in-sotropic growth appears with increasing calculatingterface cell in solid phase outside.step, as shown in Fig.2 (b), and the growth direc-Fig.3 (a) and (b) indicate that the solute con-tion of dendrite paralleling to x and y axes growscentration in solid phase is nearly the same at the be-more rapidly than that in other directions at 6 000ginnin中国煤化工while，there is nosteps, which means that the x and y axes are prefer-sidebrYHCNMH(Gpeed of dendrite a-ential directions, and between the two axes, the ve-long xt aund y axes IS iguer tiall that in other direc-locity at solid/ liquid interface is so small that it leadstions. Before 8 000 steps, dendrite becomes coarserto the transformation of liquid into interface cell and and coarser in its primary dendritic arms. The seco-●32●J. Iron & Steel Res.，Int.Vol.120.700 0001 4661 464间b)(a)(b)0.400 000,14651 460(0)00000 a(心)1455+十I 0.400 000问)00000中三喜0.700 00001 4551450(e)(I)0.2mm0.000000.2 nm1450.1445(a) 2 000;(b) 8000; (c) 10 000;(d) 12 000;(e) 14 000; (f) 16 000Fig. 3 Average solute concentration of each element(b) 8 000;(c) 10 000;during solidification at different step(d) 12 00;(e) 14 000;(f) 16 000Fig. 4 Temperature field evolution during solidificationndary dendritic arms appear gradually in both sidesat different stepof the primary dendritic arm with solidifying time in-creasing, as shown in Fig.3 (c) to (f). On the otherfield at different calculating steps, and the color barhand, it can be seen that the solute concentration a-stands for the temperature of each cell. From Eqn.head of dendrite is higher than that far away from(1)，the following relation can be deduced easily:aTthe interface and the solute concentration in front isatV●(aVT)higher than that in the root. The solute concentra-tion depends on the solute diffusivity in liquid and=V●(a.VT)(18)solid phase. With the dendrite advancing，the diam-eters of primary dendritic arms no longer grow, and=V●(aVT)+Q .athe concave parts of dendrite become unstable.af.Firstly, the concave part of dendrite leads to the ap-where Q=pLi ; a is heat transfer rate:)t'pearance of sidebranches. Fig.3 (e) and (f) showk\k。f。●k.+f●kythat the solute concentration around primary den-a1-pCpl，as=pips, a1=p(f。●Cp+f. Cn)一(19)dritic and secondary dendritic arms is much higherIf cooling rate in solid phase is assumed as athan other parts because the liquid solute concentra-constant, the cooling rate in liquid phase and inter-tion cannot diffuse into the neighboring solid phase.face can be calculated as follows:Therefore, solidification rate of the solute -rich cellas = const, al = const●is slower than that of solute -poor cell, which causesmicro-segregation. In addition, Fig.3 (e) and (f)Cp(f. k.+f。. k.)a1= const●(20)show that the sidebranch size in the tip side of den-(fi.C+f,. C%)●k,dritic arm is bigger than that at root because soluteFrom Eqn. (19), it can be obtained that a,>a1is accumulated in the root of secondary dendritic> an ，which means that the heat transfer rate in solidarms which decreases the undercooling and the di-phase is higher than that in others, and the heatameter of dendritic arm root.tran中国煤化工smallest, indicatingFig. 4 is the temperature field during solidifica-thatphase decreases moretion, showing the influences of different physicalCCHCNMHGe and interface zonesproperties of the liquid cell, the solid cell and inter-if the cooling conditions are the same.face cell on heat transfer. Fig. 4 is the temperatureFig. 4 (a) shows that the temperature of the in-No.1Numerical Simulation of Dendrite Evolution during Solidification Process●33●got is about 1 466 C except the center part whichd. Swizerland: Trans Tech Publications Ltd, 1998.has a temperature of about 1 465 C. Fig.4 (b)[2] Gungor M N. A Statistically Significant Experimental Tech-nique for Investigating Microsegregation in Cast Alloys [J].shows that the temperature difference in the wholeMetallurgical Transaction, 1989, 20A(11): 2529 2533.simulated zone enlarges for different physical prop-[3] Kurz W, Giovanola B, Trivedi. Theory of Microstructure Deerties, in which liquid temperature is about 1 464 Cvelopment during Rapid Solidification [UJ]. Acta Mater, 1986，34(5): 823 -830. .and in the center part of ingot, the temperature is a-[ 4] Dilthey U, Pavlik V. Proceedings of the Modeling of Casting,bout 1 460 C. The temperature differences becomeWelding and Advanced Solidifcation Proceses-VIlI [C].stabile as the solidifying time increasing, as shownThomas BG , Beckermann C, Eds. 1998: 589-596.in Fig.4 (c) to (f). Fig.4 (c) and (d) also reflect[5] Dithey U, Pavlik V, Reichel T. In Mathematical Modeling ofWeld Phenomena 3 [C]. Cerjak H, Bhadeshia H K D, Eds.that the heat transfer rate in solid cell is higher thanUK; The University of Cambridge. 1997. 85-105.in liquid cell and interface cell so that the tempera-[6] Nastac L. Numerical Modeling of Solidification Morphologiesture rises gradually from the center of simulatedand Segregation Patterns in Cast Dendritic Alloy [J]. Actazone to the edge of dendrite and the surrounding liq-Mater, 199 47(17): 4253-4262.uid. Moreover, the heat transfer rate in the primary[7] ZHU M F, Kim J M, HONG C P. Modeling of Globular andDendritie Strueture Evolution in Soldifcation of an Al-7mass%dendritic arm is higher than that in secondary den-Si Alloy [J]. ISU Inter, 2001, 41(9): 992-998.dritic arm or tertiary dendrite arm, as shown in[8] ZHU MF, HONG C. A Modified Cellular Automaton ModelFig.4 (e) and (f). In addition, the temperature infor the Simulation Dendrite Growth in Solidification Alloys [J].the root of dendrite decreases more quickly than thatISU Inter, 2001, 41(5): 436-445.in the tip. Fig.4 (e) and (f) show that the tempera-9] Sang Hyun Cho, Toshimitsu Okane, Takaeru Umeda. TheContribution of the Nucleation Process to Grain Formation inture is below 1450 C and 1 445 C in the center partCalculating Solidification Microstructure by CA-DFD[J]. Sci-of ingot, respectively, and the difference betweenence and Technology of Advanced Materials, 2001, 2(1); 241-the edge of dendrite and its neighbor liquid cell isnot very clear. The temperature surrounded by the[10] Jacot A, Rappaz M. A Pseudo Front Tracking Technique forthe Modeling of Solidification Microstructures in Multi-Com-second dendrite arm and tertiary dendrite arm isponent Alloys [J]. Acta Mater, 2002, 50(8): 1909-1926.higher than the dendrite temperature, which eflects[11]Gonzalez- Cinca R. Dendritic Shape at High Undercolingsthe unstable heat transfer during solidification.[J]. Physica A: Statistical Mechanics and Its Applications,2002, 314(1-4): 284 290.4 Conclusions[12] Kobabyshi R. Modeling and Numerical Simulation of DendriticCrystal Growth [J]. Physica D, 1993, 63(3-4): 410-423.(1) The continuous model can describe the[13] Damir Juric, Gretar Tryggvason. A Front Tracking Methoddendrite evolution during solidification, includingfor Dendrite Solidification [J]. J Comput Phys, 1996, 123the interface instability, secondary and tertiary den-(1): 127 148.dritic arm appearing.Udaykumar H s, Mao L. Shape lnterface Simulation of Den-(2) Micro-segregation is also simulated in localdritic Solidification [J]. Int J Heat Mass Transfer, 2002, 45(24): 4793-4808.zone surrounded by sidebranches and tertiary den-[15] Beckmann C, Li Q, Tong X. Microstructure Evolution in Edrite arms. The results show that heat is preferen-quiaxed Dendrite Growth [J]. Science and Technology of Ad-tially extracted along preferential grain growth di-vanced Materials, 2001, 2(1): 117-126.[16] Alain Karma, Wouter Jan Rappel. Phase Field of Dendriticrection.Sidebranching with Thermal Noise [J]. Physical Review E,(3) This method can deal with the discontinu-1999, 60(4): 3614-3625.ous physical properties during the phase transition.[17] Christof Eck, Peter Knabner, Segrey Korotov. A Two ScaleMoreover, the method is simple to understand andMethod for the Computation of Solid Liquid Phase Transitionsdon't need to track the position of solid/ liquid inter-with Dendritic Microstructure [J]. J Comput Phys, 2002, 178(1): 58-80.face during solidification.[18] LI Qiang, GUO Qiao-yi, LI Dian-zhong, et al. ContinuousReferences:Method to Describing Dendrite Evolution during Soldificaion[J]. Chinese Phys Lett, 2004, 21(1); 143-145 (in Chinese).[1] Kurz w, Fisher DJ. Fundamentals of Solidification [M]. 4th中国煤化工MYHCNMHG