N.ArmourS.Dost

1　Introduction

Crystal grow th is a solidification process that produces solid materials of single crystalline structure.Most bulk single crystals of semiconductors are grown from the liquid phase(melt and solution)[Müller and Ostrogorsky(1994);Dost and Lent(2007);Mechighel(2013)].

The solution grow th techniques such as Liquid Phase Electroepitaxy(LPEE),Liquid Phase Diffusion(LPD),and the Traveling Heater Method(THM)and the melt grow th techniques such as Bridgman,Czochralski(CZ),Vertical Gradient Freezing(VGF),and Float Zone(FZ)are used to grow single crystals of semiconductor materials that are extremely valuable in the development of electronic and optoelectronic devices[Dost and Lent(2007);Lappa(2005);Gulluoglu and Tsai(2000)].Wafers cut from such crystals require uniform electrical properties to ensure a higher reproducibility and yield of such devices[Deal(2004)].Consequently,there is enormous economic motivation to produce homogeneous crystals.Hence,developments in the science and technology of crystal grow th are continually sought to meet increasingly precise wafer specifications required by chip manufacturers.

The literature on crystal grow th processing,both experimental and/or numerical,is very rich.These studies have addressed various issues involving the application of magnetic fields,convective flow structures, flow stability, flow suppression,interface stability,and grow th rates.For instance,convective flow structures,the effect of applied magnetic fields on flow structures,grow th interface shapes,and grow th rates in solution crystal grow th techniques such as LPD may be found in[Yildiz,Dost and Lent(2005);Yildiz,Dost and Lent(2006);Yildiz,Dost and Yildiz(2006);Armour,Dost and Lent(2007);Yildiz and Dost(2007a);Yildiz and Dost(2007b);Armour,Yildiz,Yildiz and Dost(2008)],LPEE in[Armour,Sheibani,and Dost(2006);Liu,Okano,and Dost(2002)],and THM in[Okano,Nishino,Ohkubo and Dost(2002);Liu,Dost,Lent and Redden(2003);Kumar,Dost and Durst(2007)].Although theoretical,numerical,and/or experimental investigations of morphological instability in melt/solution crystal grow th over the past half-century have provided a great deal of insight into the physics of the phenomenon,accurate prediction of instabilities in semiconductors has proven,however,to be difficult.This is mainly attributable to the fact that semiconductors are facet-forming materials with highly anisotropic properties.Interfacial kinetics play a major role in the solidification of these materials,but the underlying kinetics mechanisms[Deal(2004);Sampath and Zabaras(2001);Mechighel(2013)]influencing morphological stability are still largely undiscovered.This is due in part to experimental difficulties associated with fine control of the process parameters and an inability to monitor processes occurring at the solid/liquid interface.

Recently,a series of dissolution experiments were performed in a crucible similar to that used in the LPD and Melt-Replenishment Czochralski(Cz)growth systems[Armour and Dost(2009)](schematically shown in Fig.1)with and without the application of a static vertical magnetic field.The measured concentration profiles from the samples processed with and without the application of magnetic field showed that the amount of silicon transport into the melt was slightly higher in the samples processed under magnetic field,and there was a substantial difference in dissolution interface shape indicating a change in flow structure in the melt.Without an applied field,a flat and stable interface was observed.In the presence of an applied field,however,the dissolution interface remained flat in the center but dramatically curved back into the source material near the crucible wall.This indicated a far higher dissolution rate at the edge of the silicon source.

In order to shed light on these experimental observations,numerical simulations using both axisymmetric and 3-D simulation domains were carried out for this dissolution process.Various levels of applied magnetic field were considered.The objectives were i)to predict the experimental observations on dissolution rates and enhanced mass transport near the wall in the melt,ii)to predict the flow patterns and concentration distributions in the melt,iii)to have a better understanding for the diffusion process in the melt which w ill be very beneficial in grow th of SiGe single crystals from the germanium side,iv)to exam ine the in fluence of convective flow on dissolution of silicon,and v)to determine the limit of applicability of the axisymmetric model in predicting these effects.

The present simulation results showed that both the axisymmetric and 3-D models could predict the experimentally observed enhanced dissolution rates under magnetic field and also the associated changes in flow structure in the melt.The axisymmetric model was also capable of predicting the dissolution interface shape under magnetic field.However,the enhanced silicon dissolution near the crucible wall under magnetic field was only predicted by the 3-D model;showing the threedimensionality of the flow field in the melt.

2　General Considerations and Mathematical Background

The setup used in the simulation is the arrangement where the silicon seed was floating on top of the germanium melt.In this case,the silicon seed covers the melt free surface.This arrangement is similar to the crucible stacking used in the LPD grow th system for SiGe material[Yildiz,Dost and Lent(2005)].A schematic for the material con figuration used in this work is shown in Fig.1.The grow th crucible is of cylindrical shape and thus the samples are of cylindrical shapes with 24 mm diameter as clearly indicated in Fig.1 in the model domain.Initially,the domain is equalized at 800?C and it is then suddenly immersed into an 1100?C isothermal region.

In the present model,the liquid phase is taken as the Si-Ge solution domain(?)in the grow th crucible(Fig.1),and the solid phases represent the silicon polycrystalline feed(?S),and the ampoule-crucible wall(?Quatrz).The liquid phase is considered as a mixture of two viscous and heat conducting incompressible Newtonian fluids,and is assumed to be a dilute binary solution(Ge-rich)of the solute(Si)and the solvent(Ge).

Figure 1:Schematic view of the experimental setup(upper left),the simulation axisymmetric domain(upper right),the simulation 3-D domain(lower right),and a picture of the LPD grow th furnace(lower left).

Furthermore,in the present model the follow ing assumptions regarding the Si-Ge solution are considered:(i)thermo-physical properties,such as,the thermal con-ductivitykL,the dynamic viscosityμL,and the thermal diffusivityand the diffusion coefficientDLof the melt are taken as constants,(ii)only laminar flow regime is considered,(iii)no viscous dissipation,and(iv)the Boussinesq approximation holds:In the SixGe1?xliquid solution of the present setup,the density of the mixture is expressed as:whereρ0is the density of the reference(starting)liquid solution which is the molten germanium.

For the assumptions regarding the silicon source,thermo-physical properties,such as,the thermal conductivitykS,and the thermal diffusivity■ are taken as constants,moreover,no species diffusion in the silicon solid is considered(DS≈0).Similar thermo-physical properties for the quartz solids(ampoulecrucible wall)are taken constants.The materials(silicon solid,Si-Ge solution,and quartz)thermo-physical data were compiled from[Yildiz(2005);Yildiz and Dost(2005)].

Figure 1 shows schematically the vertical cross section of the present grow th system.Under the above assumptions,the three-dimensional time-dependent governing equations describing the melt( fluid) flow,heat and solute transport for the liquid phase are written in cylindrical coordinatesx(r,?,z).It is easy to write two-dimensional(axisymmetric)equations from three-dimensional equations by simply dropping the dependency on the azimuthal angle(?)in the field dependent variables.

2.1　The Governing Equations for the liquid phase

2.1.1The Navier-Stokes Equations

Under the above assumptions,the principles of conservation of mass and balance of momentum yield the following Navier-Stokes equations,in which the velocity pressure formulation is adopted:

whereρL,u(x,t),withx(r,?,z)are the mass density of the m ixture(satisfying the Boussinesq approximation),mass average velocity vector,which has a radial,an axial and an azimuthal(tangential)components denoted by(u,vandw)and the position vector(note that an axisymmetric flow field expressed in terms of the cylindrical coordinate systemx(r,?,z)where all flow variables are independent of the azimuthal angle?).For simplicity since the velocity in the solid phases is zero,as w ill be shown later,in the previous equations the flow velocity in the liquid phase is designed byu(x,t)insteaduL(x,t).

For the flow velocity field,the homogeneous Dirichlet condition(no-slip)was assumed on the entire crucible(excluding the interface(Γint)t,which is illustrated latter),thus:

As initial condition(at t=0),since the velocity field was speci fiedu0(x)=0(divergence–free velocity field)over the liquid domain ?tthus

2.1.2Energy transport

The balance of energy over the liquid domain ?twith the boundary Γtleads to the following time-dependent convection-conduction equation,where neither heat generation(Joule heating)nor absorption are considered,furthermore,no heat transfer by radiation is assumed:

whereTLrepresents the temperature in the liquid phase(?t).

The temperature boundary conditions on all the internal boundaries of the domain(Fig.1)(excluding the interface(Γint)t,which is illustrated latter)are assumed as continuity,which reads:

where(where,irefers tomaterial1 andmaterial2,i.e.solution/quartz).The conditionTini=800°Cwas speci fied as the initial condition,which reads:

2.1.3Solute transport

Since there is no mass transport through the crucible wall(impermeable),and by assum ing no solute(silicon)diffusion in the solid(Ds＜＜DL),the solute transport may then be given by the follow ing time-dependent mass transport equation:

whereCLandDLare the concentration and the coef ficient of the diffusion of the liquid solute in the liquid phase(Si-Ge melt),respectively.

Since it was assumed that the crucible boundaries are not permeable for species transport,thus,the solutal boundary conditions(excluding the interface(Γint)t,which is illustrated latter)associated with Eq.(8)are:

wherehCrepresents a specified species flux(Neumann condition).Furthermore,the following condition is specified as the initial condition:

2.1.4Electric charge balance equation and the magnetic body force

The crucible is subjected to an applied axial static magnetic fieldB=B0ezwith the uniform field intensity ofB0as shown in Fig.1.For metallic liquids,the magnetic body forceFEacting on the points of the liquid(domain ?t)may be taken simply as

where the contribution of electric charge is neglected,and the current densityJis given by the Ohm’s law by

The current densityJis governed by the conservation of current

with the assumption that the induced fields due to the applied magnetic field are negligible[Mechighel(2013)].This is a good approximation for metallic liquids since the magnetic Reynolds numbers are suf ficiently small[Giessler,Sievert,Krieger,Halbedel,Huelsenberg,Luedke and Thess(2005);Kakimoto and Liu(2006)].In the above equationsφis the electric potential andσethe electric conductivity of the liquid phase.In this case the magnetic body force becomes

Combining(11)and(12)leads to the follow ing Poisson’s equation

that governs the electric potential distribution.

2.2　Solid phases

The problem can be greatly simpli fied when assum ing that there is no species diffusion in the solid,so thatus=0(whereusis the velocity vector in the solid phases).The only governing equation needed is the energy balance:

whereρs,cps,andksare the density,specific heat and thermal conductivity of the solid phases(Fig.1),respectively.Recall that the solid phases in the model are the polycrystalline Si source;and the quartz crucible-ampoule which are modeled as a rigid,heat conducting solid,as shown in Fig.1.

In the silicon solid phase,?solidrepresents ?S,(Fig.1),and the temperature is denoted byTS,the conditions on Γint(illustrated later),continuity conditions on silicon source/crucible boundary and top silicon/vacuum boundary together provide the required boundary conditions for resolution of the heat transport problem.(The vacuum above the surface is assumed to be a gas of negligible viscosity and conductivity).

For the quartz crucible,?solidrepresents ?Quartz,(Fig.1),and the temperature is denoted byTQuartz,the thermal boundary conditions for the outer vertical wall,and bottom surface of the quartz crucible are expressed as a speci fied temperature(1100?C).At the top surface of the domain the flux was set to zero.Furthermore,perfect thermal contacts and continuous heat flux at the silicon seed/crucible,solution/crucible and inner crucible boundaries were assumed.

2.3　Conditions at silicon source/melt(dissolution)interface

The interfacial system of the present study consists of two non-reacting components,Ge and Si,and two phases;the liquid phase(L)(Si-Ge melt assumed as a dilute binary liquid m ixture(Ge-rich)of the solute(Si)and the solvent(Ge))and the solid phase(S)(silicon source).

The solid/liquid(source/melt or dissolution)interface is a complicated boundary where the imposed conditions must account for interfacial transport of energy and mass.Indeed,the dissolution interface is not only a source and sink for thermal energy and components,but it moves as phase change takes place.The conditions required to establish such an interface can be determined using thermodynamics,(i.e.assuming the interface is an equilibrium state of the system).Furthermore,the driving force for phase transformation(i.e.dissolution in the present study)can be described by a deviation from this equilibrium state[Mechighel(2013)].

General boundary conditions at a solid/liquid interface may be written as:

whereis the melting temperature of pure Ge,mLis the liquidus slope given by the two-component phase diagram,Υ is the local mean curvature of the interface(positive for convex surfaces viewed from the liquid),ΔS0is the entropy of fusion per unit volume for pure(Ge),γSLis the interfacial energy per unit area between the solid and liquid phases,nis the unit normal vector to a point solid/liquid interface pointing toward the liquid;uintis the velocity of the solid/liquid interface in the reference( fixed)frame at that point,andLfV＜0 is the latent heat of fusion(enthalpy of fusion)for the solid(Si source)per unit volume,assumed identical for the solute andis the equilibrium distribution coef ficient of the solute.Here all quantities are associated with the solid/liquid(S/L)interface boundary.

Assumptions at the dissolution interface

-The inclusion of the liquid velocityuat the dissolution interface in previous equations is necessary if the density of the solid and liquid are not equal[Deal(2004)].The resulting shrinkage or expansion upon dissolution w ill alter the velocity of the interface and affect the transport of components and heat,anduaccounts for this.However,since the experimentally observed dissolution rate of silicon source into the solution is very small,so the velocity components(u,vandw)of fluid particles at the dissolution interface can be neglected in the present system(thusu·n=0)[Armour and Dost(2009)].

-As mentioned earlier,assuming the interface is at an equilibrium state of the system,means that the local thermodynamic equilibrium(TS=TL=Teq)is reached at the interface,the solute concentration at the dissolution interface is determined from the Si-Ge binary phase diagram,namelyCeq=fn(Teq),which is written for the liquidand solidequilibrium compositions for the Ge-rich side of the phase diagram.On the dissolution interface we adopt the following saturation concentration(in molar concentration of silicon)predicted from the Si-Ge phase diagram:

and

where

whereandTdisare respectively the equilibrium mass fractions of the solute and the temperature at the dissolution interface

2.4　Symmetry conditions(used only for the axisymmetric model):

At the axis of symmetry,we use symmetry conditions.Also for physical( finite)results,it is required that the radial velocity component be zero[Yildiz and Dost(2005)].Thus,

3　Numerical resolution

3.1　The ST-GLS formulation for unsteady incompressible fluid flow

Since the dissolution interface velocity is assumed to be negligible in the present study so the space-time Galerkin/least-square(ST-GLS) finite element formulation for incompressible fluid flows is used here.This formulation is used in order to prevent numerical instabilities encountered when using the standard Galerkin finite element method(details on this technique can be found in works by Tezduyar and Coworkers[Tezduyar(1992);Tezduyar,Mittal,Ray and Shih(1992)].In the present work,this stabilized finite element formulation was utilized usingequalorderinterpolation velocity-pressure elements as proposed by[Tezduyar,Mittal,Ray and Shih(1992)].

Firstly,the formulation assumes that the spatial domain is fixed in time;this implies that the subscripttis dropped from the symbols ?tand Γt[Tezduyar(1992)].In the space-time finite element formulation,the time interval(0,tmax)is partitioned into subintervalsIn=(tn,tn+1),wheretnandtn+1belong to an ordered series of time levelsThe space-time slabQnis de fined as the space-time domain ?×In.The lateral surface ofQnis denoted byPn;this is the surface described by the boundary Γ,asttraversesIn.Similarly,Pnis decomposed into(Pn)Dand(Pn)Nwith respect to the type of boundary condition being imposed(as shown in Eq.3,since we have only the Dirichlet part,including the dissolution interface(Eq.16),thus

Finite element discretization of a space-time slabQnis achieved by dividing it into elements,where(nel)nis the number of elements in the space-time slabQn.Associated with this discretization,by adopting the standard notation of[Tezduyar(1992);Tezduyar,Mittal,Ray and Shih(1992)])we define the following finite element interpolation and variational function spaces for velocity and pressure:

For various ways of calculating matrixsee for instance[Tezduyar(1992);Shakib(1988);Hauke and Hughes(1998);Polner(2005);F?rster(2007)].In the present study,this was implemented in COMSOL package(see[COMSOL Multiphysics Modeling Guide(2008)]).

The finite element interpolation functions are discontinuous in time.The fully discrete equations can be solved one space-time slab at a time(a fractional-step procedure).The memory needed for the global matrices involved in this method is quite extensive.Iteration methods can be employed to substantially reduce the cost involved in solving the linear equation systems arising from the space-time finite element discretization.Here,resolution was achieved using the generalized minimal residual(GMRES)method[Saad and Schultz(1986);COMSOL Multiphysics Modeling Guide(2008)].

3.2　Stream line Upwind Petrov-Galerkin(SUPG)formulation for the heat transport

For the heat transport equation,we adopt the SUPG stabilization given by[Ed Akin and Tezduyar(2004)].Consider the initial boundary value problem given by the unsteady convection-conduction equation,Eq.5,and its associated boundary and initial conditions,Eqs.6 and 7.The principle of the SUPG stabilization technique is de fined by taking a perturbation as the following:This w ill be introduced in the weak form of the problem,with being the weighting function.Thus

whereR(w)is the residual of the equation of energy,Eq.5.

Using a suitably-de fined finite-dimensional trial solution and interpolation function spaces(further details,see[Ed Akin and Tezduyar(2004)]),the stabilized SUPG finite element formulation of the previously written energy equation with boundary and initial conditions can be written as follows.

Findsuch that

Herenelis the number of elements and ?eis the element domain corresponding to elementis the SUPG stabilization parameter.For various ways of calculatingsee for instance[Tezduyar and Osawa(2000);Ed Akin and Tezduyar(2004)].In the present study,this was implemented in COMSOL package.

Finally,the finite element discretization of this weak form yields a system of semidiscrete equations fort∈(0,tmax).In order to trace the transient response,this system of sem i-discrete equations can be advanced in time by suitable finite difference schemes such as theθfamily methods.A fully implicit method known as"Backward Differentiation Formulas"(BDF)is used[Donea and Huerta(2003);COMSOL Multiphysics Modeling Guide(2008)].

3.3　SUPG FEM for solute transport equation

The same SUPG stabilization technique is used for the solute transport equation.Using a suitably-de fined finite-dimensional trial solution and interpolation function spaces,the stabilized finite element formulation of the previously written species equation(Eq.8)with boundary and initial conditions(Eqs.9 and 10)can be written as follows:

Finduch that

Hereis the SUPG stabilization parameter(for solute transport equation).For various ways of calculatingsee for instance[Franca;Frey and Hughes(1992);Tezduyar and Osawa(2000)].In the present study,this was implemented in COMSOL package.

3.4　Heat transport by conduction in the solid phases

The classical Standard Galerkin formulation[Donea and Huerta(2003);COMSOL package]is used for Eq.14,with the corresponding boundary and initial conditions.Further details about the mathematical modeling and solution methodology are illustrated in the article by[Mechighel,Armour,Dost and Kadja,(2011)].

4　Results and discussion

The 3-D simulation conducted in the absence of magnetic field exhibited a decaying dissolution phenomenon,and showed an expected diffusion-dominated behavior in the region close to the dissolution interface,in agreement with the experimental observation of[Armour and Dost(2009)](Fig.2).Transport into the Ge-rich melt is relatively slow and continues to slow down as the concentration gradient flattens.Mathematically,this means that in Eq.8 the diffusive term,is dominant compared with the convective term,(u·?CL),in this region.In fact,this is due to the physics of the experimental con figuration.In this system,the silicon source placed at the top of the melt dissolves,under the applied temperature pro file,into the Ge-rich melt,and then diffuses downward in the opposite direction to the gravity-induced buoyancy force.Due to the large density difference between the Si solute and the Ge-rich melt,the lighter silicon solute is buoyant,and the diffusion of silicon acts to stabilize the melt against natural(thermosolutal)convection.This makes the silicon transport in this system diffusion dominated,and naturally leads to slower dissolution rates.At the center of the melt and near the crucible wall,however,convective flow is present as illustrated in the follow ing manner.

The nature of the flow in the melt can be quantified by determining the relative contributions of convection and diffusion by calculating the ratio of the buoyancy force to the viscous force,known as the thermal Grash of number,GrT=■ForGrTvalues less than 104,viscosity dominates and convection is minimal(Stokes flow).However,ifGrTis greater than 104,convection is present,and the diffusion boundary layer is truncated b for the presents tu dy(with the melt radiusRof12mm and a temperature gradientGTof 5 K/cm)we calculateThis means that a relatively weak convective flow(thermal convection)was present in the melt.In addition,an upper bound for the size of the diffusion boundary layer is estimated by■using the smallest dissolution rate.ForRdis=2.5 mm/hr,this is 36 mm.Therefore,since convection was present,the diffusion boundary layer in the system was less than 36 mm.

Figure 2:Average height of silicon dissolved into Ge-rich melt:(a)experiment conducted with no field and under a 0.8 Tesla magnetic field level[Armour and Dost(2009)];(b)comparison between the numerically predicted and measured dissolution heights.Prediction is in good agreement with experiment.Indeed,the maximum relative error between the simulation and experimental results is less than 6%.

To provide further insight into the nature of the melt flow and mixing in the present crucible con figuration,numerical simulations using both axisymmetric and 3-D models were carried out without and with the application of magnetic field.The numerically predicted concentration and melt flow velocity fields are presented in Figs.3,4,5.The components of melt flow velocity(u,vandw)are also graphically shown in Figs.6,7,8.

Computed flow patterns in the melt after 20 minutes of dissolution are graphically represented by flow velocity vectors in Fig.3.The rotation of the thermosolutal convection cell is counterclockwise(Fig.3a).In 3-D simulation,this is an annular roll(Fig.3band 3d).Near of the sample) fluid moves up toward the dissolution interface,and then it travels along the dissolution interface toward the center, finally turning downwards in the direction of the crucible bottom surface.Along the dissolution interface near the center of the sample, flow is very weak,and consequently,the solute transport is diffusion dominated in this region.Slow dissolution(lower dissolution rate)in this region results in a semi-mixed melt.

However,a stronger flow along the vertical wall(sides)of the domain leads to a better mixing in the melt and gives rise to relatively faster dissolution rates of the source at the edge.Indeed,the concentration pro file around the dissolution interface shows that silicon is being mixed away from the region near the crucible wall towards the center of the melt(Fig.3b and 3c).

Figure 3:Results after 20 m in of dissolution with no magnetic field.Arrows indicate the flow structure,and isolines/isosurfaces illustrate the solute concentration distribution.In the 3-D simulation the dissolution pattern clearly indicates that it is more concentrated in the center.

The stable flow structure(caused by silicon buoyancy in the melt)results in a flat dissolution interface.The most significant difference with the application of the magnetic field is the shape of the dissolution interface as seen in Fig.9(taken from reference[Armour and Dost 2009]).In Fig.9,we see the photos of the polished surfaces of the cross-sections of two processed samples,which are obtained by cutting through the center of the cylindrical processed samples at the end of the dissolution experiment.In this figure,the photo of the cross-section of the processed sample on the right(under magnetic field)clearly shows the enhanced dissolution(more silicon dissolution)near the crucible wall compared with that in the photo on the left(with no magnetic field).On the left,we see that the silicon dissolution interface remains almost flat,but on the right,more silicon dissolved near the wall leading to a curved interface that is visible in the photo.

Figure 4:Results after 20 m in of dissolution under a 0.2 Tesla field.Arrows indicate the flow structure and isolines/isosurfaces illustrate the solute concentration distribution.In the 3-D simulation the dissolution pattern clearly indicates that it is more concentrated in the center(but is less pronounced than the simulation with no applied field).

As mentioned earlier,the change observed in the interface shape indicates the effect of the applied magnetic field on the melt flow.This expected change in the melt flow is well predicted by the present 3-D numerical simulation(as shown in Figs.3 and 5).The rotation of the thermosolutal convection roll isclockwiseunder magnetic field(Figs.5b and 5d).One may explain the numerically predicted flow and dissolution patterns as follows.The downward strong warm flow,due to the combined effect of gravitational and magnetic body forces,is observed in the regions near to the crucible wall(Fig.5b and 5d),due to solute rejection in the melt.The flow starts from the region near the interface edge(i.e.edge of the silicon source)and moves towards the crucible bottom surface,supplying and enriching the melt with the dissolved solute.The solute transport by convection in this region leads to further dissolution from the source(faster dissolution at the edge of the silicon source).Near the crucible bottom surface,a warmer downward flow turns towards the center of the melt and brings the dissolved solute to the center(Fig.5b and 5d)which results in well-mixed regions along the sides and the center of the melt.The observed stronger downward flow(driven by the gravitational body force)at the center was reversed in direction with the applied magnetic field(Fig.3b and 5b).In the region close to the center of the interface,the melt flow is almost absent and the solute transport is diffusion-dominated.This naturally leads to far slower dissolution rates at the center and results in semi-mixed region at the center of the melt near the dissolution interface.It is important to note that the axisymmetric model does not predict these trends(Fig.3a and 5a).This indicates the significant impact of the tangential flow component on the interface shape.

Figure 5:Results after 20 m in dissolution with a 0.8 Tesla applied field.Arrows indicate the flow structure and isolines/isosurfaces illustrate the concentration pro file.In the 3-D simulation the dissolution pattern indicates that it is more concentrated near the edge.

Figure 6:Results after 20 m in of dissolution with no applied field.Slices indicate the radial and tangential velocity components plotted at different radial positions(z=2 mm,12 mm and 23 mm).

Figure 7:Results after 20 min of dissolution with 0.2 Tesla applied field.Slices indicate the radial and tangential velocity components plotted at different radial positions(z=2 mm,12 mm and 23 mm).

Figure 8:Results after 20 m in of dissolution with 0.8 Tesla applied field.Slices indicate the radial and tangential velocity components plotted at different radial positions(z=2 mm,12 mm and 23 mm).

The magnetic field appears to be acting to m ix silicon away from the region near the crucible wall into the center,a direction inverse to the case with no applied field(Fig.3d and 5d).This action creates a higher concentration gradient near the crucible wall at the edge of the silicon source,and increases dissolution(Fig.5c).Thus,in both cases(with and w ithout magnetic field)the melt flow gives rise to mixed and semi-mixed regions near the dissolution interface.Under the 0.8 Tesla magnetic field,much higher dissolution rates were realized at the edge of the source material,which is in agreement with the experimental observations(Fig.9).

In the 3-D simulations,magnetic field levels under(B0＜0.2 Tesla)appear not to have a significant impact on the flow structure.For magnetic field levels in the range of(0.2≤B0＜0.6 Tesla),the applied field does not have significant effects on the flow structure nor on the interface shape but causes a non-uniform concentration distribution in the melt(Fig.4b and 4c).Magnetic field levels in the range of(0.6≤B0≤0.8 Tesla),however,have a significant impact on the flow structure and lead to a more uniform concentration distribution in the melt(Figs.5b and 5c).

Figure 9:Polished cross-sections of the quenched samples from the experiments conducted with no field on the left and that with a 0.8 Tesla field on the right[Armour and Dost 2009].On the left,we see that the silicon dissolution interface remains almost flat,but on the right,more silicon dissolved near the edge of the dissolution interface(curved area).

Due to a small increase in the amount of dissolved silicon(Fig.2),it appears that the(0.2≤B0＜0.6 Tesla)range of applied field does not have a significant impact on the flow structure.Particularly the 3-D simulation results well predict this observation.In the center region and near the crucible bottom,both the radial and tangential flow velocity components appear to be smoothed out and slightly damped(Figs.6c,6d,7c and 7d).In the regions ahead of the dissolution interface near the crucible wall,however,these flow velocity components appear to be smoothed and slightly enhanced(Fig.6c,6d,7c and 7d).Also,except the regions near to the wall,where the axial flow appears to be enhanced,the sim ilar trend may be observed for the axial flow(Figs.6e and 7e).

The effect of the 0.8 Tesla level of applied field on flow structures is signi ficant as shown by the 3-D simulation.In the center region and near the crucible bottom,both the radial and tangential flow velocity components appear to be weakened and are inverted in direction(Figs.6c,6d,8c and 8d).In the regions near the interface and close to the crucible wall,however,these two flow velocity components appear to get slightly stronger and are inverted(Figs.6c,6d,8c and 8d).Furthermore,except the regions near to the crucible wall where the axial flow is stronger and inverted in direction(Figs.6e and 8e),the axial flow velocity appears to be weakened(for instance at the positions z=2 mm and 12 mm)and almost suppressed(for instance at z=23 mm)in the rest of the domain.These observations were noted by[Armour and Dost(2009)].

In most crystal grow th systems,a static magnetic field is frequently utilized to suppress thermosolutal convective fluid flow in the melt.In the present crucible con figuration,however,it appears that the melt flow is not suppressed.Instead,the already weak stable flow structure is strengthened under the effect of applied magnetic field due to the buoyancy of the silicon solute.Indeed,the applied field strengthens and inverts the upward flow near to the crucible wall and weakens and inverts the downward flow in the center region and near the crucible bottom(Figs.3 and 5).

In fact,an external magnetic field,aligned perfectly with the axis of the grow th cell(z-direction),gives rise to a magnetic body force in the radial plane that balances the vertical gravitational body force,and consequently strengths both the radial and the tangential flow velocity components,and suppresses the axial flow.This fluid flow behavior can easily be understood from the fact that maximum electromagnetic damping of the fluid flow occurs when the velocity fielduis oriented orthogonal to the magnetic field vectorbecause of the particular form of the electromagnetic damping forceIt should be noted that the axisymmetric model fails to predict these trends.

The present 3-D numerical simulation results con firm the experimental observations of[Armour and Dost(2009)].The observed dissolution structure may be valuable to the LPD crystal grow th of Si-Ge system,where the curvature of the grow th interface evolves with time.It may be possible to use this effect to better control interface shape.

5　Conclusions

The numerical simulations conducted for the dissolution process of silicon into germanium melt using both axisymmetric and 3-D simulation domains led to the following conclusions.

-Transport of silicon in the Ge-rich melt where silicon is being dissolved from the top exhibits a diffusion-dominated behavior(in the region near to the dissolution interface and the center).This behavior is due to the silicon solute buoyancy in the heavier Ge-rich melt,which has been previously explained experimentally by[Armour and Dost(2009)]and predicted numerically by[Yildiz,Dost and Lent(2005)].

-The present 3-D simulation predicts that the application of a static magnetic field strengthens the radial and tangential flow velocity components in the region near the crucible wall at the edge of the silicon source,and thus leads to a stronger convective flow pattern in the melt close to this region.Consequently,this change in the melt flow gives rise to a significant mixing of silicon away from the region near the crucible wall into the center of the melt.This phenomenon may have an application in controlling grow th interface geometry and maintaining a constant grow th interface curvature during the LPD grow th of SiGe single crystals.

-Axial mixing in the melt does not appear to significantly increase and the concentration pro file seems to evolve in a diffusion-dominated manner.

-The numerical results of the present simulations indicate that the silicon dissolution was slightly enhanced under an applied vertical magnetic field.This enhancement peaked for the field levels between 0.2 to 0.5 Tesla.The magnetic level of 0.8 Tesla appears to be an optimum value for obtaining uniform concentration distribution.This observation can be attributed to the altered flow structure in the melt due to the applied magnetic field.

-Simulation results also show that,as the magnitude of the applied magnetic field increases,the flow pattern becomesthree-dimensional.The radial and tangential flow velocity components increase,causing higher dissolution rates at the edge of the silicon source.These observations were well predicted by the 3-D model while the axisymmetric model failed to do so.

Acknowledgement:The financial support provided by the Natural Sciences and Engineering Research Council of Canada(NSERC)and the Canada Research Chairs(CRC)Program is gratefully acknowledged.

Armour,N.;Dost,S.(2009):Effect of an applied static magnetic field on silicon dissolution into a germanium melt.Journal of Crystal Growth,vol.311,pp.780–782.

Armour,N.;Dost,S.;Lent,B.(2007):Effect of free surface and gravity on silicon dissolution in germanium melt.Journal of Crystal Growth,vol.299,pp.227–233.

Armour,N.;Sheibani,H.;Dost,S.(2006):Grow th of cadmium zinc telluride by liquid phase electro epitaxy.Crystal Research and Technology,vol.41,no.10,pp.939–945.

Armour,N.;Yildiz,M.;Yildiz,E.;Dost,S.(2008):Liquid Phase Diffusion Grow th of SiGe Single Crystals under Magnetic Fields.ECS Transactions,vol.16,no.10,pp.135–146.

COMSOL Multiphysics Modeling Guide(2008):by COMSOL AB.Tegnérgatan 23 SE-111 40 Stockholm,Sweden.

Deal,A.(2004):Enhanced morphological stability in Sb-doped Ge single crystals,PhD thesis,Florida University,USA.

Donea,J.;Huerta A.(2003):Finite Element Methods for Flow Problems.John Wiley&Sons,Ltd.

Dost,S.;Lent,B.(2007):Single Crystal Growth of Semiconductors from Metallic Solutions.Elsevier,Amsterdam,the Netherlands.

Ed Akin,J.;Tezduyar,T.E.(2004):Calculation of the advective lim it of the SUPG stabilization parameter for linear and higher-order elements.Computer Methods in Applied Mechanics and Engineering,vol.193,pp.1909–1922.

F?rster C.(2007):Robust methods for fluid-structure interaction with stabilized finite elements,PhD thesis University of Stuttgart,Germany.

Franca,L.P.;Frey,S.L.;Hughes,T.J.R.(1992):Stabilized finite element methods:I.Application to the advective-diffusive model.Computer Methods in AppliedMechanics and Engineering,vol.95,pp.253–276.

Giessler C.,Sievert C.,K rieger U.,Halbedel B.,Huelsenberg D.,Luedke U.,Thess A.(2005):A Model for Electromagnetic Control of Buoyancy Driven Convection in Glass Melts.FDMP:Fluid Dynamics&Materials Processing,vol.1,no.3,pp.247-266.

Gulluoglu,A.N.;Tsai,C.T.(2000):Effects of grow th direction on tw in formation in GaAs crystal grown by vertical gradient freeze method.CMES:Computer Modeling in Engineering and Sciences,vol.1,no.1,pp.85-89.

Hauke,G.;Hughes,T.J.R.(1998):A comparative study of different sets of variables for solving compressible and incompressible flows.Computer Methods in Applied Mechanics and Engineering,vol.153,pp.1–44.

Kakimoto,K.;Liu,L.(2006):Flow Instability of Silicon Melt in Magnetic Fields.FDMP:Fluid Dynamics&Materials Processing,vol.2,no.3,pp.167-173.

Kumar,V.;Dost,S.;Durst,F.(2007):Numerical modeling of crystal grow th under strong magnetic fields:An application to the travelling heater method.Applied Mathematical Modelling,vol.31,pp.589–605.

Lappa,M.,(2005):Review:Possible strategies for the control and stabilization of Marangoni flow in laterally heated floating zones.FDMP:Fluid Dynamics&Materials Processing,vol.1,no.2,pp.171??187.

Liu,Y.C.;Okano,Y.;Dost,S.(2002):The effect of applied magnetic field on flow structures in liquid phase electroepitaxy–a three dimensional simulation model.Journal of Crystal Growth,vol.244,pp.12–26.

Liu,Y.C.;Dost,S.;Lent,B.;Redden,R.F.(2003):A three-dimensional numerical simulation model for the grow th of CdTe single crystals by the traveling heater method under magnetic field.Journal of Crystal Growth,vol.254,pp.285–297.

Mechighel,F.(2013):Modélisation de la convection lors d’un changement dephase : Stabilisation par champ magnétique.(editor MorelP.)Presses Académiques Francophones‘PAF’,Saarbrücken,Germany.

Mechighel,F.;Armour,N.;Dost,S.;Kadja,M.(2011):Mathematical modeling of the dissolution process of silicon into germanium melt.TWMS Journal of Applied and Engineering Mathematics(Turkic Word Mathematical Society),vol.1,no.2,pp.200–222.

Müller,G.;Ostrogorsky,A.G.(1994):Convection in melt growth.Chapter 13,Handbook of Crystal Grow th.(editor Hurle D.T.J.),North-Holland/Elsevier,vol.2,pp.711??814.

Okano,Y.;Nishino,S-S.;Ohkubo,S-S.;Dost,S.(2002):Numerical study of transport phenomena in the THM grow th of compound sem iconductor crystal.Journal of Crystal Growth,vol.238–239;pp.1779–1784.

Polner,M.(2005):Galerkin Least-Squares Stabilization Operators for the Navier-Stokes Equations:A Unified Approach,PhD thesis University of Twente,AE Enschede,The Netherlands.

Saad,Y.;Schultz,M.H.(1986):GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systems.SIAM Journal of Scientific and Statistical Computing,vol.7,no.3,pp.856–869.

Sam path,R.;Zabaras,N.(2001):Numerical study of convection in the directional solidi fication of a binary alloy driven by the combined action of buoyancy,surface tension,and electromagnetic forces.Journal of Computational Physics,vol.168,pp.384–411

Shakib,F.(1988):Finite element analysis of the compressible Euler and Navier-Stokes equations,Ph.D.Thesis,Department of Mechanical Engineering,Stanford University,Stanford,California.

Tezduyar,T.E.(1992):Stabilized Finite Element Formulations for Incompressible Flow Computations.Advances in Applied Mechanics,vol.28,pp.1–44.

Tezduyar,T.E.;M ittal,S.;Ray,S.E.;Shih,R.(1992):Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity pressure elements’.Computer Methods in Applied Mechanics and Engineering,vol.95,pp.221–242.

Tezduyar,T.E.;Osawa,Y.(2000):Finite element stabilization parameters computed from element matrices and vectors.Computer Methods in Applied Mechanics and Engineering,vol.190,pp.411-430.

Yildiz,M.(2005):A Combined Experimental and Modeling Study for the Grow th of SixGe1?x Single Crystals by Liquid Phase Diffusion(LPD),PhD thesis,University of Victoria,Canada.

Yildiz,E.;Dost,S.(2007a):A numerical simulation study for the combined effect of static and rotating magnetic fields in liquid phase diffusion grow th of SiGe.Journal of Crystal Growth,vol.303,pp.279–283.

Yildiz,M.;Dost,S.(2005):A continuum model for the Liquid Phase Diffusion grow th of bulk SiGe single crystals.International Journal of Engineering Science,vol.43,pp.1059–1080.

Yildiz,M.;Dost,S.(2007b):Incorporation of surface tension to interface energy balance in crystal grow th.Crystal Research and Technology,vol.42,no.9,pp.914–919

Yildiz,M.;Dost,S.;Lent,B.(2005):Grow th of bulk SiGe single crystals by liquid phase diffusion.Journal of Crystal Growth,vol.280,pp.151–160

Yildiz,M.;Dost,S.;Lent,B.(2006):Evolution of the grow th interface in liquid phase diffusion grow th of bulk SiGe single crystals.Crystal Research and Technology,vol.41,no.3,pp.211–216.

Yildiz,E.;Dost,S.;Yildiz,M.(2006):A numerical simulation study for the effect of magnetic fields in liquid phase diffusion grow th of SiGe single crystals.Journal of Crystal Growth,vol.291,pp.497–511.