- Email: [email protected]
Computational modeling of turbulent melt flow in CdZnTe crystal growth
Computational Materials Science 25 (2002) 316–328 www.elsevier.com/locate/commatsci
Computational modeling of turbulent melt flow in CdZnTe crystal growth ern R. C y a, P. Jelınek b, P. Prikryl
c,*
a
b
Department of Structural Mechanics, Faculty of Civil Engineering, Czech Technical University, Th akurova 7, 166 29 Prague 6, Czech Republic Institute of Physics, Academy of Sciences of the Czech Republic, Cukrovarnick a 10, 162 53 Prague 6, Czech Republic c Mathematical Institute, Academy of Sciences of the Czech Republic, Zitn a 25, 115 67 Prague 1, Czech Republic Received 14 December 2001; accepted 17 January 2002
Abstract A computational model of CdZnTe crystal growth including turbulent phenomena in the melt is presented as a erny, Kalbac and Prikryl [Comp. Mater. Sci. 17 (2000) 34]. The ternary CdZnTe continuation of the previous work by C system is treated as pseudobinary in the temperature range under consideration. The phase interface between the solid and liquid is modeled as a discontinuity surface. The k–e model for low turbulent Reynolds numbers is used for the description of the turbulent processes. The numerical solution of the model is done using the Galerkin finite element method in two-dimensional approximation with cylindrical symmetry. The problem of moving boundary is solved by a front-fixing method. The results obtained with the turbulent model are compared to those of the previous laminar model and the influence of including the turbulence into the model is discussed. It is shown that the laminar and the turbulent models exhibit significant differences in the structure of melt flow, concentration profiles, and also in the movement of the phase interface for larger time periods. Ó 2002 Elsevier Science B.V. All rights reserved. PACS: 02.70.Dh; 47.27.Eq; 47.27.Te; 73.61.Ga; 81.10.Fq Keywords: CdZnTe crystal growth; Computational modeling; Vertical Bridgman method; Turbulent melt flow; k–e turbulence model; Low-Reynolds-number model
1. Introduction Computational modeling of the process of crystal growth is an important aid in optimizing the real experiment to produce high-quality single *
Corresponding author. Tel.: +420-22-22090-728; fax: +42022-2221-1638. E-mail address: [email protected] (P. Prikryl).
crystals of semiconducting materials. The computational simulations can save significant amount of time and money in both the design and the optimization of crystal growth processes, and can improve the understanding of macroscopic transport effects in many systems as well. Today, still more complex models are developed describing more complicated systems close to the experimental reality.
0927-0256/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 2 ) 0 0 3 4 3 - 9
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
The natural convection driven by the gradient of temperature and concentration belongs to the factors which can affect the shape and stability of the solid/liquid interface in all crystal growth techniques. Generally, larger driving forces can result in the appearance of turbulent flow patterns, and consequently in the formation of very irregular structures instead of the desired single crystals. In addition, the turbulence itself can enhance the transport processes in the system and lead to the production of heat due to the dissipation of kinetic energy, which can have an impact on the process of the crystal growth from the melt. Therefore, it seems highly desirable to analyze the possibility of occurrence of such turbulent processes in the case of crystallization. In this paper we introduce a new model of crystal growth from the melt based on the laminar model presented in [1] previously that takes also the turbulent processes into account. Two different k–e models for low turbulent Reynolds numbers are employed and the results obtained for the CdZnTe crystal growth by the Bridgman method with these two models and the previous laminar model are compared. The computational results indicate that the turbulent model provides a more adequate description of the processes under consideration.
2. Mathematical model We will consider the following model situation as sketched in Fig. 1. A cylindrical ampoule of length L and radius R containing a binary alloy in the liquid state is encircled by two furnaces, which is implemented in the model by imposing the effective boundary conditions for the temperature at the lateral area of the ampoule. Both the base sides are thermally insulated. The upper furnace is kept at a fixed temperature well above the melting point of the material, whereas the lower furnace is adjusted to a temperature below the melting point. Therefore, a temperature gradient is realized in the area of the adiabatic zone between the furnaces. Typically, the solidification process begins at the bottom of the ampoule. However, both the ampoule and the facilities setting up the boundary
317
Fig. 1. Schematic figure of the experimental situation.
conditions for temperature can move up and down to generate the required temperature changes in the ampoule, and thus the solidification process can also be initiated from any intermediate state when the solid/liquid interface is somewhere between the top and the bottom of the ampoule. It is believed that all incompressible fluid flows, whether they possess laminar or turbulent character, can be described mathematically by the three-dimensional Navier–Stokes equations. Hence, the most natural approach for modeling turbulence would be to discretize and solve the Navier–Stokes equations on a very fine grid so that we would be able to represent even the tiniest eddies occurring in the turbulent flow. This is what is called direct numerical simulation (DNS) of turbulence. However, the requirements imposed by this approach on the memory capacity and the speed of the computer used are so vast that the modeling approach is changed mostly in the way we describe shortly in what follows (see also [2,3] and Chapter 10 of [4] in this respect). In most practical studies of turbulent flow, the main interest is in the macroscopically observed
318
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
mean values rather than in every tiny detail of the flow field. In addition, a certain decomposition into different scales is usually natural and inherent in the flow itself. Therefore, we focus our attention on computing the mean flow characteristics using a feasible number of mesh points rather than performing a complete DNS and trying to obtain every microscopic detail. In fact, this idea––that goes back to Reynolds––is essential for all common models of turbulence employed in the everyday practice today. Proceeding in this way, we decompose all the quantities occurring in the Navier–Stokes equations into a mean part and the small variations known as fluctuations and try to arrive to a system of equations describing the behavior of the mean parts. In the most general formulation, the mean part is formed by applying a filter possessing certain required properties (see [4] for a detailed description) but for our purposes this filter may be just a componentwise temporal average. If we apply this procedure to the governing equations of our previous laminar model [1] we obtain the governing equations for the time averaged variables of the turbulent model. However, these equations will contain new unknown terms (Reynolds stresses) which represent the turbulent strain arising from the mutual interaction of the particular components and depend on the velocity fluctuations. Hence the resulting equations do not constitute a closed system any more (we now have more unknowns than equations) and we thus need to postulate further equations describing the relations between the Reynolds stresses and the mean quantities we would like to determine. Solving this closure problem in turbulence modeling we are forced to employ hypotheses and approximations which are only of heuristic nature and can be obtained from empirical information and experimental data. Various approaches are used in this respect, all of them being based on the assumption that the approximation of the Reynolds stress tensor used in the model depends on the filtered velocities exclusively. Excellent reviews of such models allowing the closure of the set of governing equations can be found, for example, in [5] or [6]. The most widely used model today probably is the two-equation k–e model of Laun-
der and Spalding [7], which closes the set of the filtered balance equations by introducing two additional variables, the turbulent kinetic energy k and the dissipation rate e, and adding two supplementary equations for the determination of them. This model has been extended later to apply to more complex situations and an up-to-date study of it is the subject of the book by Mohammadi and Pironneau [8], where a thorough derivation and analysis of the k–e turbulence model may be found. In our case the natural convection generated in the liquid by the temperature and concentration gradients must take the buoyancy effects into account correctly. Near the heating or cooling walls, the flow can be strongly affected by the viscous forces as well. One therefore attempts to modify the original k–e model in such a way that damping effects near the wall are accounted for or a boundary layer near the wall is described in an appropriate way [9]. In this paper we employ the strategy followed in the low-Reynolds-number models where the relevant constants in the k–e model are multiplied by appropriate modifying functions that take the distance to the nearest wall into account. In particular, our simulations have been done with the low-Reynolds-number model of Lam and Bremhorst [10] and with that of Davidson [11]. Owing to the above described modeled situation we can reduce our problem to a cylindrical-symmetry case. Therefore, the set of the governing equations for the filtered mean quantities can be written in the liquid part as (see [1,4,8]) ovlr vlr ovlz þ þ or r oz blT l 1 blT ðT l T0l Þ þ blC ðCAl CA;0 Þ l oT oT l oT l þ vlz þ vlr ot or oz
blC l 1 blT ðT l T0l Þ þ blC ðCAl CA;0 Þ l oCA oC l oC l þ vlr A þ vlz A ¼ 0; ot or oz þ
ð1Þ
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
ql
ql
ovlr ovl ovl þ vlr r þ vlz r ot or oz 2 l l op o vr 1 ovlr vlr o2 vlr þ ðgl þ glt Þ ¼ þ þ or or2 r or r2 oz2 2 l 1 o vr 1 ovlr vlr o2 vlz þ ðgl þ glt Þ þ þ 3 or2 r or r2 oroz oðgl þ glt Þ ovlr vlr 1 ovlr vlr ovlz þ2 þ þ þ ; or or r 3 or r oz ð2Þ ovlz ovl ovl þ vlr z þ vlz z ot or oz 2 l opl o vz 1 ovlz o2 vlz ¼ þ ðgl þ glt Þ þ þ oz or2 r or oz2 2 l 1 o vr 1 ovlr o2 vlz þ 2 þ ðgl þ glt Þ þ 3 oroz r oz oz oðgl þ glt Þ ovlz ovlr oðgl þ glt Þ þ þ þ2 or oz or oz l ovz 1 ovlr vlr ovlz þ þ oz 3 or r oz l Þ; þ ql0 blt gðT l T0l Þ ql0 blC gðCAl CA;0
ql
ql
319
oCAl oC l oC l þ vlr A þ vlz A ot or oz 1 o gl oCAl oT l rql Dl þ t þ rql dl ¼ r or rc or or l l l o l gt oCA l l l oT q D þ þqd þ ; oz rc oz oz
ð5Þ
ok l ok l ok l þ vlz þ vlr ot or oz l 1 o gt ok l o glt ok l l l r g þ g þ ¼ þ r or oz rk or rk oz þ P þ G ql el ;
ð6Þ
l l oel l oe l oe q þ vz þ vr ot or oz 1 o glt oel o glt oel l l r g þ g þ ¼ þ r or oz re or re oz l
2
þ c1e f1 ðP þ GÞ ql c2e f2 ð3Þ
l oT oT l oT l ðql clp blT pl Þ þ vlz þ vlr ot r oz l l 1 o g oT oC l r Kl þ t þ rBl A ¼ r or rT or or l l o g oT oC l Kl þ t þ Bl A þ oz rT oz oz l l oCA oC oC l blC pl þ vlr A þ vlz A ot or oz l l l oH g oCA oT l ql þ dl þ Dl þ t or rc or or l l l oH g oCA oT l ql þ dl þ Dl þ t oz rc oz oz l l l l ovr 1 ovr vr ovlz ovr l l þ 2ðg þ gt Þ þ þ or 3 or r oz or l l l l ovr ovz ovr ovz þ ðgl þ glt Þ þ þ oz or oz or l l l l ovz 1 ovr vr ovlz ovz l l ; þ 2ðg þ gt Þ þ þ oz 3 or r oz oz ð4Þ
el : kl
ð7Þ
The meaning of the symbols used in the above equations is as follows. The terms P and G represent the production of turbulent kinetic energy and we have (see [3]) l l ovr 1 ovlr vlr ovlz ovr P ¼ 2ðgl þ glt Þ þ þ or 3 or r oz or l l l l ovr ovz ovr ovz þ ðgl þ glt Þ þ þ oz or oz or l l l ovz 1 ovr vlr ovlz ovz ; þ 2ðgl þ glt Þ þ þ oz 3 or r oz oz ð8Þ oT l oC l þ ql0 blC g A : ð9Þ oz oz Further, v is the velocity, bT and bC the thermal and the solutal expansion coefficient, respectively, T is the temperature, CA the concentration of the component A of the mixture, q is the density, p the pressure, g the dynamic viscosity and gt the turbulent dynamic viscosity, g ¼ j~ gj is the gravity, K is the thermal conductivity, cp the specific heat, B the Dufour coefficient, H the difference in partial G ¼ ql0 blT g
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
320
specific enthalpies of the components of the mixture, D is the diffusion coefficient, d is the Soret coefficient, rc is the turbulent Schmidt number, rT the turbulent Prandtl number, k the turbulent kinetic energy, and e means the dissipation rate of k. The superscript l denotes the liquid phase, q0 ; T0 ; CA;0 are values corresponding to a reference state. The turbulent dynamic (eddy) viscosity glt is defined as 2
glt ¼ cl fl
kl el
ð10Þ
in the framework of the low-Reynolds-number k–e models, where cl is an empirical constant and fl is a function depending on the chosen model. The functions f1 , f2 and fl in the models under consideration are given as follows. Lam–BremhorstÕs model [10]: 3 0:05 f1 ¼ 1:0 þ ; f2 ¼ 1:0 expðRe2s Þ; fl 20:5 2 fl ¼ ð1:0 expð 0:0165Red ÞÞ 1:0 þ : Res DavidsonÕs model [11]: 3 0:14 f1 ¼ 1:0 þ ; fl f2 ¼ ½1 0:27 expððRes Þ2 Þ ½1 expðRed Þ ; " # 3:4 fl ¼ exp : 2 ð1 þ 0:02Res Þ
re ¼ 1:3;
c1e ¼ 1:44; rc ¼ 0:9;
c2e ¼ 1:92;
rT ¼ 0:9;
rk ¼ 1:0; pffiffiffiffi d Red ¼ k l l ; g
2
kl Res ¼ l l ; ge
where vamp is the velocity of motion of the ampoule, taken as positive for the downward direction, the other quantities being denoted in the same way as in the liquid, and the superscript s denotes the solid phase. At the solid/liquid interface assumed to have the form z ¼ Zðr; tÞ, the following balance equations can be written: ð16Þ
ql vlr viN ¼ ql vlr ðnr vlr þ nz vlz Þ þ nr ðpl p0 Þ
l ovlr ovr ovlz nz ðgl þ glt Þ þ or oz or l l l 2 ovr vr ovz þ ðgl þ glt Þnr þ þ ; ð17Þ 3 or r oz
2nr ðgl þ glt Þ
ql vlz viN þ qs vamp viN
where d is the distance to the nearest wall. This completes the description of the model equations in the liquid phase. In the solid we have [1] vsr ¼ 0;
ð12Þ
ð13Þ ps ¼ p0 ¼ const:; s oT oT s qs csp vamp ot oz 1 o oT s oC s o oT s rK s Ks ¼ þ rBs A þ r or oz or or oz s s oC oCA oC s þ B s A qs H s vamp A oz ot oz s s s oH oC oT qs Ds A þ ds þ or or or s s s oH s oT s s oCA q D þ þd ; ð14Þ oz oz oz s oCA oC s qs vamp A ot oz 1 o oC s oT s rqs Ds A þ rqs ds ¼ r or or or s o oC oT s qs Ds A þ qs ds þ ; ð15Þ oz oz oz
ðql qs ÞviN ¼ ql ðnr vlr þ nz vlz Þ þ nz qs vamp ;
Some of the constants and functions are common for both the models: cl ¼ 0:09;
vsz ¼ vamp ¼ const:;
ð11Þ
¼ ql vlz ðnr vlr þ nz vlz Þ qs nz v2amp þ nz ðpl p0 Þ l l ovr ovlz l l ovz l l nr ðg þ gt Þ 2nz ðg þ gt Þ þ oz oz or l l l 2 ovr vr ovz þ ðgl þ glt Þnz þ þ ; ð18Þ 3 or r oz
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
the solid into the liquid, LðT0 Þ is the latent heat of melting at the reference temperature T0 , and p0 is a reference pressure. The above equations can be developed analogously as with the previous laminar model [1], however, the turbulent processes have been taken into account here. Besides the balance conditions, two additional conditions at the solid/liquid interface are formulated. They are given by the liquidus and solidus curves of the phase diagram,
½ql LðT0 Þ þ ðql clp qs csp ÞðT T0 Þ þ ðql qs Þcsp T0 pl þ p0 viN ¼ ql LðT0 Þðnr vlr þ nz vlz Þ þ ½ql ðnr vlr þ nz vlz Þ þ qs nz vamp csp T0 þ ½ql clp ðnr vlr þ nz vlz Þ þ qs csp nz vamp ðT T0 Þ pl ðnr vlr þ nz vlz Þ p0 nz vamp l glt oT l l l oCA þB nr K þ rT or or gl oT l oC l þ Bl A nz K l þ t rT oz oz oT s oC s þ nr K s þ Bs A or or s s s oT s oCA þ nz K þB oz oz l l gt oCAl l oT l l l þd D þ nr q H rc or or l l l gt oCA l oT l l l þd D þ nz q H rc oz oz s s oC oT þ nr q s H s D s A þ ds or or oC s oT s þ nz qs H s Ds A þ ds ; oz oz
321
T l ¼ Tm ðCAl Þ;
ð21Þ
T s ¼ Tc ðCAs Þ;
ð22Þ
where Tm ðCAl Þ is the liquidus curve and Tc ðCAs Þ is the solidus curve. Finally, we formulate the local thermodynamic equilibrium condition in the form Tl ¼ Ts
ð19Þ
ðql CAl qs CAs ÞviN ¼ ql CAl ðnr vlr þ nz vlz Þ þ qs CAs nz vamp gl oCAl oT l nr ql Dl þ t þ dl rc or or l l g oCA oT l þ dl nz ql Dl þ t rc oz oz s s oC oT þ nr qs Ds A þ ds or or oC s oT s ; þ nz qs Ds A þ ds oz oz ð20Þ
where viN is the velocity of the interface in the normal direction, nr , nz are the components of the unit normal vector to the interface, directed from
ð23Þ
at the phase interface. Appropriate initial and boundary conditions corresponding to the particular crystal growth techniques are formulated in a common way corresponding to the particular experimental setup. The conditions for T, C, and v are analogical to those used in the previous laminar model [1]. We employ constant initial conditions for k and e and Dirichlet boundary conditions for these variables (these will be discussed later, see also [8]). The boundary conditions at the axis of symmetry are treated in the standard way. If we consider a system cooled from bottom, then the boundary and initial conditions can be written in the following manner: vlr ðr; L; tÞ ¼ vlz ðr; L; tÞ ¼ vamp ;
r 2 ½0; R ;
vlr ðR; z; tÞ ¼ vlz ðR; z; tÞ ¼ vamp ;
z 2 ½ZðR; tÞ; L ;
ð24Þ
ð25Þ vlr ð0; z; tÞ ¼ vamp ; ovlz ð0; z; tÞ ¼ 0; or
z 2 ½Zð0; tÞ; L ;
z 2 ½Zð0; tÞ; L ;
T l ðR; z; tÞ ¼ T2 ðz; tÞ > Tm ðCAl Þ;
ð26Þ ð27Þ
z 2 ½z2 ðtÞ; L ; ð28Þ
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
322
ok l ð0; z; tÞ ¼ 0; or
z 2 ½Zð0; tÞ; L ;
ð29Þ
oel ð0; z; tÞ ¼ 0; or
z 2 ½Zð0; tÞ; L ;
ð30Þ
k l ðR; z; tÞ ¼ k1 ¼ const:;
z 2 ½ZðR; tÞ; L ;
ð31Þ
e ðR; z; tÞ ¼ e1 ¼ const:;
z 2 ½ZðR; tÞ; L ;
ð32Þ
k l ðr; Z; tÞ ¼ k1 ¼ const:;
r 2 ½0; R ;
ð33Þ
el ðr; Z; tÞ ¼ e1 ¼ const:;
r 2 ½0; R ;
ð34Þ
l
l
k ðr; Zðr; tÞ; tÞ ¼ k1 ¼ const:; l
e ðr; Zðr; tÞ; tÞ ¼ e1 ¼ const:;
r 2 ½0; R ;
ð35Þ
r 2 ½0; R ;
ð36Þ
z2;0 z þ T2 ðz; 0Þ z2;0 z1;0 z z1;0 ; r 2 ½0; R ; z 2 ½z1;0 ; z2;0 ; z2;0 z1;0
T s;l ðr; z; 0Þ ¼ T1 ðz; 0Þ
ð49Þ T s ðr; z; 0Þ ¼ T1 ðz; 0Þ;
r 2 ½0; R ; z 2 ½0; z1;0 ; ð50Þ
l CAl ðr; z; 0Þ ¼ C1;0 ðz; 0Þ;
r 2 ½0; R ; z 2 ½Z0 ðrÞ; L ; ð51Þ
s ðz; 0Þ; CAs ðr; z; 0Þ ¼ C1;0
r 2 ½0; R ; z 2 ½0; Z0 ðrÞ ; ð52Þ
k l ðr; z; 0Þ ¼ const:;
z 2 ½z2 ðtÞ; L ; r 2 ½0; R ;
s;l
oT ðR; z; tÞ ¼ 0; or
z 2 ½z1 ðtÞ; z2 ðtÞ ;
T s ðR; z; tÞ ¼ T1 ðz; tÞ < Tc ðCAs Þ;
z 2 ½0; z1 ðtÞ ;
ð38Þ
oT l ðr; L; tÞ ¼ 0; oz
r 2 ½0; R ;
ð39Þ
oT s ðr; 0; tÞ ¼ 0; oz
r 2 ½0; R ;
ð40Þ
oT s;l ð0; z; tÞ ¼ 0; or
el ðr; z; 0Þ ¼ const:;
z 2 ½z2 ðtÞ; L ; r 2 ½0; R ;
ð54Þ
where z1 ðtÞ ¼ z1;0 þ vT t, z2 ðtÞ ¼ z2;0 þ vT t, z1;0 6 z2;0 , and vT is the velocity of motion of the facilities asserting the boundary conditions for temperature.
3. Computer implementation z 2 ½0; L ;
CAl ðr; L; tÞ ¼ C1 ¼ const:;
r 2 ½0; R ;
ð41Þ ð42Þ
oCAs ðr; 0; tÞ ¼ 0; oz
r 2 ½0; R ;
ð43Þ
oCAs;l ð0; z; tÞ ¼ 0; or
z 2 ½0; L ;
ð44Þ
oCAs;l ðR; z; tÞ ¼ 0; or
z 2 ½0; L ;
ð45Þ
Zðr; 0Þ ¼ Z0 ðrÞ;
r 2 ½0; R ;
ð46Þ
vlr ðr; z; 0Þ ¼ vlz ðr; z; 0Þ ¼ vamp ; r 2 ½0; R ; z 2 ½0; L ; ð47Þ l
ð53Þ
ð37Þ
T ðr; z; 0Þ ¼ T2 ðz; 0Þ;
r 2 ½0; R ; z 2 ½z2;0 ; L ; ð48Þ
Our approach here is analogous to that employed with the previous laminar model, the principal difference being the necessity to process the equations for k and e in addition. Therefore, we are rather brief in this Section referring the erny et al. [1] interested reader to the paper by C where all the delicate details of the discretization procedure have been discussed and necessary references quoted. To begin the survey of our computational procedure we note that we use the penalty method for the continuity Eq. (1), which means that a slightly modified form of this equation is used in the computer implementation, where an artificial term ep pl with a small ep is added (see [12] in this respect). This makes it possible to eliminate the pressure from the system of equations to be solved in the liquid and to obtain a system for velocities, k and e only. The motivation here is the same as in
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
[1]: the elimination of pressure from the unknown variables not only makes an appropriate discretization feasible but also reduces the computing time (we are not interested in the pressure field in our practical simulations). For the discussion of the application of the penalty method the interested reader is referred, e.g., to [13,14]. The system of equations obtained thus in the liquid and that in the solid are discretized by using the Galerkin finite element method with bilinear 4 node quadrilateral elements in space, linear basis functions in time. The space discretization leads to two systems of first-order differential equations, one of them corresponding to the solid, the other one to the liquid phase. These are then subject to the time discretization. As the overall procedure has been described and critically discussed in [1] we refer the reader there for more details again. However, we feel necessary to repeat that the practical situations for which our model is intended are strongly diffusion dominated and that the Galerkin FEM does not suffer from spurious node-to-node oscillations here. Finally, the matrices and vectors of the solid and liquid phases are coupled together to obtain the global system of algebraic equations that represents the whole space region ½0; R ½0; L and describes the change of the system variables from time t to t þ Dt in our discrete model. In this way we get (similarly as in [1]) ½HP fwgtþDt ¼ ½PF fwgt þ fFg;
ð55Þ
where D s s s ; . . . ; Tm1 ; CA;m1 ; fwg ¼ T1s ; CA;1 l vlr;1 ; vlz;1 ; T1l ; CA;1 ; k1l ; el1 ; . . . ; ET l vlr;n ; vlz;n ; Tnl ; CA;n ; knl ; eln :
ð56Þ
Here, the subscripts 1; . . . ; m 1 and 1; . . . ; n are the indices of the nodal points, and [HP], [PF], {F} are the resulting global matrices and the global vector which correspond to both the solid and liquid phases and are localized in accordance with the ordering of the global vector fwg. In the combination of matrices and vectors we employ the local thermodynamic equilibrium condition (23) and the FEM analogs of the interface bal-
323
ances of internal energy and mass of the component A so that we obtain nonzero components of the right-hand side vector {F} at the positions corresponding to these equations. The system (55) is doubly nonlinear as its terms [HP], [PF], and [F] are not only functions of fwg but also of the a priori unknown position Z of the phase interface, which is also to be determined in the course of the solution. The moving boundary problem (55), (21) and (22) is solved now by the successive approximation method and a front-fixing method that transforms the two-dimensional variable space regions ½0; R ½0; Zðr; tÞ , ½0; R ½Zðr; tÞ; L corresponding to the solid and liquid phases, respectively, to the fixed space domains ½0; 1 ½0; 1 . Due to the character of the modeled crystal growth techniques and the assumed form of the phase interface it is possible to use analogous transformations here as are those for the one-dimensional problems introduced by Landau [15]. Hence, we put f¼
r ; R
g¼
z ; Zðr; tÞ
ð57Þ
in the solid phase and f¼
r ; R
g¼
z Zðr; tÞ ; L Zðr; tÞ
ð58Þ
in the liquid phase, transform the problem and work with the independent variables f; g instead of r; z. The transformations are applied to all the element matrices and vectors defined before, however, the system of equations (55) remains formally the same and has similar nonlinear properties after the transformation. The resulting system of nonlinear algebraic equations is then solved at each time step by two nested iteration procedures using the successive approximation method with underrelaxation. The inner iterations are due to the nonlinearity of the field equations, the outer iterations are performed with the aim to satisfy all the conditions prescribed at the a priori unknown phase interface for the particular time instants. Our numerical tests have shown that only one possible choice of the outer iterative scheme leads to a stable computational algorithm. This scheme consists in applying the phase-diagram
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
324
and the local equilibrium conditions (21)–(23) as Dirichlet conditions at the interface and using the interface internal energy balance equation as the convergence criterion. The situation here is similar to that in our previous laminar model. Therefore, also the overall scheme of one time step of our algorithm with two nested iteration procedures is analogous but we feel it reasonable to present it here for the convenience of the reader: 1. Put Z ð1Þ ðt þ DtÞ ¼ ZðtÞ, Z_ ð1Þ ðt þ DtÞ ¼ Z_ ðtÞ, ð1Þ fwgtþDt ¼ fwgt . ð2Þ 2. Compute fwgtþDt from (55) using Z ð1Þ ðt þ DtÞ, ð1Þ Z_ ð1Þ ðt þ DtÞ and fwgtþDt in the matrices [HP], [PF]. 3. If vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPn ð2Þ 2 2 u i¼1 ðwi Þ ðwð1Þ Þ i t > di ; Pn ð2Þ 2 i¼1 ðwi Þ where di is a user-defined error tolerance for the internal iteration cycle, put ð1Þ
ð2Þ
fwgtþDt ¼ fwgtþDt ; and go back to Step 2. Otherwise, put ð1Þ
ð2Þ
fwgtþDt ¼ fwgtþDt and continue with Step 4. 4. Compute Z_ ð2Þ ðt þ DtÞ using the finite element analog of the interface energy balance equation ð1Þ (19) and fwgtþDt . 5. If Z_ ð2Þ ðt þ DtÞ Z_ ð1Þ ðt þ DtÞ > de ; Z_ ð1Þ ðt þ DtÞ where de is a user-defined error tolerance for the external iteration cycle, then choose a 2 ½0; 1 , put Z_ ð1Þ ðt þ DtÞ ¼ aZ_ ð2Þ ðt þ DtÞ þ ð1 aÞZ_ ð1Þ ðt þ DtÞ; Z ð1Þ ðt þ DtÞ ¼ ZðtÞ þ DtZ_ ð1Þ ðt þ DtÞ; and go back to Step 2. Otherwise, put ð1Þ
fwgtþDt ¼ fwgtþDt ;
Z_ ðt þ DtÞ ¼ Z_ ð2Þ ðt þ DtÞ; 1 Zðt þ DtÞ ¼ ZðtÞ þ Dt½Z_ ð2Þ ðt þ DtÞ þ Z_ ðtÞ 2 and go to the next time step. The computer implementation of the numerical algorithm was performed in Fortran 77 and a computer code has been developed. The numerical experiments discussed in the next Section were done on on a PC with two processors PII 450 MHz under OS Linux and on a workstation HP C180 under OS HP–UX.
4. Numerical simulations The main object of the computer simulations we performed has been to compare the results produced by the model presented in this paper with those obtained using our previous model of CdZnTe crystal growth with laminar melt flow from [1]. As before in [1], a set of numerical calculations was performed first in order to optimize the numerical parameters of the model. We focused namely on the type of the iteration procedure, the magnitude of the pressure penalty ep , the initial values of k and e, the time step of simulation and the optimal values of the relaxation coefficients for particular variables, which strongly affect the rate of convergence. From these extensive numerical tests the following optimum parameters were obtained: ep in the range ½1010 ; 108 kg1 m2 s2 , the initial values k 102 m2 s2 and e in the range ½104 ; 103 m2 s3 , and the relaxation coefficients for velocities about 0.6, for temperature and concentration about 0.8 and for turbulent kinetic energy and its dissipation rate about 0.3. We observed that particularly the initial values of k and e had a significant influence on the stability of the algorithm. Also for the values of the pressure penalty which laid outside the above mentioned optimum interval a stable solution was not achieved. A typical number of inner iteration steps for achieving the convergence of all variables with di ¼ de ¼ 103 was of the order of 102 and the
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
number of outer iterations of the order of 101 . The optimum finite element mesh in both the liquid and the solid that guaranteed sufficient precision in a still acceptable CPU time was 20 20. The time step was controlled by limiting its maximum value by Dtmax , the optimum value being Dtmax 1 s. After finding the optimal numerical parameters we performed a set of numerical simulations with both the laminar and the turbulent models corresponding to a real experimental situation. All our numerical computations modeled the CdZnTe crystal growth using the vertical Bridgman method. The common basic parameters of the experimental setup were as follows: the radius of the ampoule R ¼ 37:5 mm, the length of the ampoule L ¼ 100 mm, the temperatures of the Bridgman furnaces 1125 and 975 °C, the distance between the furnaces 10 mm, the velocity of the ampoule vamp ¼ 106 m s1 or vamp ¼ 107 m s1 . The thermophysical properties of CdZnTe used in our simulations can be found in Table 1 of [1]. Fig. 2 shows the development of the shape of the solid/liquid interface during the CdZnTe crystal growth in the turbulent model. We can observe that an S-shape is characteristic for the initial time period studied, whereas a C-shape is established after achieving a more developed state. A comparison with the results obtained with the laminar model from [1] shows significant differ-
Fig. 2. Development of the phase interface in time for vamp ¼ 106 m s1 .
325
Fig. 3. Comparison of the position of the phase interface between the turbulent and the laminar model for vamp ¼ 106 m s1 .
ences in the case of vamp ¼ 106 m s1 (see Fig. 3). For the case of vamp ¼ 107 m s1 no differences were found in the whole range of the simulated time interval (10,000 s). Reverse movement (remelting) of the phase interface near the wall was observed with the turbulent model that can be seen in Fig. 2. Fig. 4 shows the evolution of the velocity fields characteristic for our turbulent model. The typical time stable flow pattern for the laminar model is shown in Fig. 5. While the laminar flow model always leads to a very simple flow structure with only one recirculating region, the structure of the flow calculated by the turbulent model is more complicated. Isolated recirculating flow regions arise in the domain in the vicinity of the phase interface, while on the opposite side practically no flow is observed. Especially in the vicinity of the wall and the phase interface, instabilities were observed. Also, the magnitude of the velocities was about one order lower for the turbulent flow model compared to both the velocity of the interface and the velocities in the laminar model. The computed concentration profiles are different for both the turbulent and the laminar model (see Fig. 6). The turbulent profiles are affected by the more intensive transport processes.
326
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
Fig. 6. Comparison of the concentration profiles at the axis of the ampoule for vamp ¼ 106 m s1 at time t ¼ 10 000 s.
Fig. 4. Time series of the turbulent flow pattern for vamp ¼ 106 m s1 .
Moreover, some concentration peaks were found in the produced crystal, in particular near the wall. These peaks were more pronounced in the case of the turbulent model, where also variable concentration along the phase interface was found. Finally, we note that we also compared the two above mentioned turbulent models, namely that of Davidson and that of Lam and Bremhorst. No significant differences in the structure of the velocity field or of the shape of the solid/liquid interface were found. Sometimes we had problems with obtaining negative values of the turbulent kinetic energy and its dissipation rate using the Lam–Bremhorst model. Therefore, DavidsonÕs model has been used for our practical simulations.
5. Discussion
Fig. 5. Stable structure of the laminar flow pattern for vamp ¼ 106 m s1 .
In our numerical simulations, we tested the influence of turbulent natural convection due to the thermal and solutal expansion in the liquid on the temperature, concentration and velocity fields, and on the position and velocity of the moving solid/ liquid interface. We believe that the results obtained have shown that the proposed computational model involving
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
turbulent motion in the melt gives more realistic simulation than the previous laminar model. The turbulent model simulates the whole cascade of kinetic energy, which is successively transferred from the largest eddies that are on the order of the flow domain to the smallest eddies (measuring a few millimeters or less), where the turbulent kinetic energy is converted into the internal heat of the fluid by viscous dissipation. The laminar model is not able to describe such processes in a proper way due to its disability to simulate processes running in smaller scales than those corresponding to the size of the computational mesh used. However, the small dissipative eddies play a significant role in the transport processes occurring in the system and hence they have an important impact on the whole process of crystallization. We have found distinct differences in the velocity patterns and the profiles of concentration between the two models. The turbulent melt flow was characterized by nonstationary oscillations contrary to the stable laminar pattern with two inversely oriented vortices in the ampoule. Instabilities of the turbulent flow pattern were observed close to the wall. This fact results in re-melting of the created solid phase, the oscillation of the phase interface and consequential appearance of defects in the crystals grown. This result is in agreement with the experimental observations of enhanced occurrence of defects near the wall. This phenomenon, which can explain some perturbations during the crystallization, was not (and could not be) observed with the laminar model. The results obtained with the turbulent model indicate that the flow of the melt is mainly affected by the movement of the phase interface and that the buoyancy effects play a middling role only. The buoyancy effects tend to establish a stable structure of the flow pattern, usually two vortices in the domain as we observed with our laminar model, but in the case of the turbulent fluid flow this tendency is dumped by additional viscous dissipation. Therefore, no stable structures have been observed in the turbulent flow field. The turbulent concentration profiles showed marked changes in both time and space. Peaks in the concentration field were observed in the produced solid phase. The observed instabilities near
327
the wall are in agreement with the experimental experience. From this point of view it seems to be clear that the turbulent model should be used to obtain a realistic representation of the events being in progress during the crystal growth. In the future work, we intend to perform a critical comparison with the real experimental data with the aim to assess the importance of particular factors influencing the process of crystallization and to provide more realistic values of important constants (rT ; rc ) used in simulating the CdZnTe material.
6. Conclusions The proposed turbulent model has shown a new possibility how to describe the more complicated transport processes that play a significant role in the production of CdZnTe single crystals. The turbulent model has revealed new facts about the mutual relations between the particular factors involved in the process of crystal growth, such as between the buoyancy effects and the interface kinetics, or about the diffusion of a compound of the mixture that have a pronounced influence on the final result of the production process. We conclude that the turbulent model is more appropriate to describe the transport processes in the melt during crystallization than the laminar model.
Acknowledgements This research has been supported by the Grant Agency of the Czech Republic, under grants # 202/ 99/1646 and # 201/01/1200.
References erny, A. Kalbac, P. Prikryl, Comp. Mater. Sci. 17 [1] R. C (2000) 34. [2] J.O. Hinze, Turbulence: An Introduction to Its Mechanism and Theory, McGraw-Hill, Delft, 1959. [3] W. Frost, T.H. Moulden, Handbook of Turbulence, Plenum Press, New York, 1992.
328
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
[4] M. Griebel, T. Dornseifer, T. Neunhoeffer, Numerical Simulation in Fluid Dynamics: A Practical Introduction, SIAM, Philadelphia, 1998. [5] W. Rodi, Turbulence Models and Their Application for Hydraulics, International Association for Hydraulics Research, Delft, 1980. [6] M. Lesieur, Turbulence in Fluids, Kluwer, Dordrecht, 1997. [7] B. Launder, D. Spalding, Mathematical Models of Turbulence, Academic Press, New York, 1972. [8] B. Mohammadi, O. Pironneau, Analysis of the K-Epsilon Turbulence Model, Wiley, Chichester, 1993.
[9] W.C. Patel, W. Rodi, G. Scheuerer, AIAA J. 23 (1985) 1308. [10] C.K.G. Lam, K.A. Bremhorst, ASME J. Fluids Eng. 103 (1981) 456. [11] L. Davidson, Numer. Heat Transfer 18 (1990) 129. [12] F. Thomasset, Implementation of Finite Element Methods for Navier–Stokes Equations, Springer, New York, 1981. [13] T.J.R. Hughes, W.T. Liu, A.J. Brooks, J. Comput. Phys. 30 (1979) 1. [14] M. Bercovier, M. Engelman, J. Comput. Phys. 30 (1979) 181. [15] H.G. Landau, Quart. Appl. Math. 8 (1950) 81.
Computational modeling of turbulent melt flow in CdZnTe crystal growth ern R. C y a, P. Jelınek b, P. Prikryl
c,*
a
b
Department of Structural Mechanics, Faculty of Civil Engineering, Czech Technical University, Th akurova 7, 166 29 Prague 6, Czech Republic Institute of Physics, Academy of Sciences of the Czech Republic, Cukrovarnick a 10, 162 53 Prague 6, Czech Republic c Mathematical Institute, Academy of Sciences of the Czech Republic, Zitn a 25, 115 67 Prague 1, Czech Republic Received 14 December 2001; accepted 17 January 2002
Abstract A computational model of CdZnTe crystal growth including turbulent phenomena in the melt is presented as a erny, Kalbac and Prikryl [Comp. Mater. Sci. 17 (2000) 34]. The ternary CdZnTe continuation of the previous work by C system is treated as pseudobinary in the temperature range under consideration. The phase interface between the solid and liquid is modeled as a discontinuity surface. The k–e model for low turbulent Reynolds numbers is used for the description of the turbulent processes. The numerical solution of the model is done using the Galerkin finite element method in two-dimensional approximation with cylindrical symmetry. The problem of moving boundary is solved by a front-fixing method. The results obtained with the turbulent model are compared to those of the previous laminar model and the influence of including the turbulence into the model is discussed. It is shown that the laminar and the turbulent models exhibit significant differences in the structure of melt flow, concentration profiles, and also in the movement of the phase interface for larger time periods. Ó 2002 Elsevier Science B.V. All rights reserved. PACS: 02.70.Dh; 47.27.Eq; 47.27.Te; 73.61.Ga; 81.10.Fq Keywords: CdZnTe crystal growth; Computational modeling; Vertical Bridgman method; Turbulent melt flow; k–e turbulence model; Low-Reynolds-number model
1. Introduction Computational modeling of the process of crystal growth is an important aid in optimizing the real experiment to produce high-quality single *
Corresponding author. Tel.: +420-22-22090-728; fax: +42022-2221-1638. E-mail address: [email protected] (P. Prikryl).
crystals of semiconducting materials. The computational simulations can save significant amount of time and money in both the design and the optimization of crystal growth processes, and can improve the understanding of macroscopic transport effects in many systems as well. Today, still more complex models are developed describing more complicated systems close to the experimental reality.
0927-0256/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 2 ) 0 0 3 4 3 - 9
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
The natural convection driven by the gradient of temperature and concentration belongs to the factors which can affect the shape and stability of the solid/liquid interface in all crystal growth techniques. Generally, larger driving forces can result in the appearance of turbulent flow patterns, and consequently in the formation of very irregular structures instead of the desired single crystals. In addition, the turbulence itself can enhance the transport processes in the system and lead to the production of heat due to the dissipation of kinetic energy, which can have an impact on the process of the crystal growth from the melt. Therefore, it seems highly desirable to analyze the possibility of occurrence of such turbulent processes in the case of crystallization. In this paper we introduce a new model of crystal growth from the melt based on the laminar model presented in [1] previously that takes also the turbulent processes into account. Two different k–e models for low turbulent Reynolds numbers are employed and the results obtained for the CdZnTe crystal growth by the Bridgman method with these two models and the previous laminar model are compared. The computational results indicate that the turbulent model provides a more adequate description of the processes under consideration.
2. Mathematical model We will consider the following model situation as sketched in Fig. 1. A cylindrical ampoule of length L and radius R containing a binary alloy in the liquid state is encircled by two furnaces, which is implemented in the model by imposing the effective boundary conditions for the temperature at the lateral area of the ampoule. Both the base sides are thermally insulated. The upper furnace is kept at a fixed temperature well above the melting point of the material, whereas the lower furnace is adjusted to a temperature below the melting point. Therefore, a temperature gradient is realized in the area of the adiabatic zone between the furnaces. Typically, the solidification process begins at the bottom of the ampoule. However, both the ampoule and the facilities setting up the boundary
317
Fig. 1. Schematic figure of the experimental situation.
conditions for temperature can move up and down to generate the required temperature changes in the ampoule, and thus the solidification process can also be initiated from any intermediate state when the solid/liquid interface is somewhere between the top and the bottom of the ampoule. It is believed that all incompressible fluid flows, whether they possess laminar or turbulent character, can be described mathematically by the three-dimensional Navier–Stokes equations. Hence, the most natural approach for modeling turbulence would be to discretize and solve the Navier–Stokes equations on a very fine grid so that we would be able to represent even the tiniest eddies occurring in the turbulent flow. This is what is called direct numerical simulation (DNS) of turbulence. However, the requirements imposed by this approach on the memory capacity and the speed of the computer used are so vast that the modeling approach is changed mostly in the way we describe shortly in what follows (see also [2,3] and Chapter 10 of [4] in this respect). In most practical studies of turbulent flow, the main interest is in the macroscopically observed
318
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
mean values rather than in every tiny detail of the flow field. In addition, a certain decomposition into different scales is usually natural and inherent in the flow itself. Therefore, we focus our attention on computing the mean flow characteristics using a feasible number of mesh points rather than performing a complete DNS and trying to obtain every microscopic detail. In fact, this idea––that goes back to Reynolds––is essential for all common models of turbulence employed in the everyday practice today. Proceeding in this way, we decompose all the quantities occurring in the Navier–Stokes equations into a mean part and the small variations known as fluctuations and try to arrive to a system of equations describing the behavior of the mean parts. In the most general formulation, the mean part is formed by applying a filter possessing certain required properties (see [4] for a detailed description) but for our purposes this filter may be just a componentwise temporal average. If we apply this procedure to the governing equations of our previous laminar model [1] we obtain the governing equations for the time averaged variables of the turbulent model. However, these equations will contain new unknown terms (Reynolds stresses) which represent the turbulent strain arising from the mutual interaction of the particular components and depend on the velocity fluctuations. Hence the resulting equations do not constitute a closed system any more (we now have more unknowns than equations) and we thus need to postulate further equations describing the relations between the Reynolds stresses and the mean quantities we would like to determine. Solving this closure problem in turbulence modeling we are forced to employ hypotheses and approximations which are only of heuristic nature and can be obtained from empirical information and experimental data. Various approaches are used in this respect, all of them being based on the assumption that the approximation of the Reynolds stress tensor used in the model depends on the filtered velocities exclusively. Excellent reviews of such models allowing the closure of the set of governing equations can be found, for example, in [5] or [6]. The most widely used model today probably is the two-equation k–e model of Laun-
der and Spalding [7], which closes the set of the filtered balance equations by introducing two additional variables, the turbulent kinetic energy k and the dissipation rate e, and adding two supplementary equations for the determination of them. This model has been extended later to apply to more complex situations and an up-to-date study of it is the subject of the book by Mohammadi and Pironneau [8], where a thorough derivation and analysis of the k–e turbulence model may be found. In our case the natural convection generated in the liquid by the temperature and concentration gradients must take the buoyancy effects into account correctly. Near the heating or cooling walls, the flow can be strongly affected by the viscous forces as well. One therefore attempts to modify the original k–e model in such a way that damping effects near the wall are accounted for or a boundary layer near the wall is described in an appropriate way [9]. In this paper we employ the strategy followed in the low-Reynolds-number models where the relevant constants in the k–e model are multiplied by appropriate modifying functions that take the distance to the nearest wall into account. In particular, our simulations have been done with the low-Reynolds-number model of Lam and Bremhorst [10] and with that of Davidson [11]. Owing to the above described modeled situation we can reduce our problem to a cylindrical-symmetry case. Therefore, the set of the governing equations for the filtered mean quantities can be written in the liquid part as (see [1,4,8]) ovlr vlr ovlz þ þ or r oz blT l 1 blT ðT l T0l Þ þ blC ðCAl CA;0 Þ l oT oT l oT l þ vlz þ vlr ot or oz
blC l 1 blT ðT l T0l Þ þ blC ðCAl CA;0 Þ l oCA oC l oC l þ vlr A þ vlz A ¼ 0; ot or oz þ
ð1Þ
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
ql
ql
ovlr ovl ovl þ vlr r þ vlz r ot or oz 2 l l op o vr 1 ovlr vlr o2 vlr þ ðgl þ glt Þ ¼ þ þ or or2 r or r2 oz2 2 l 1 o vr 1 ovlr vlr o2 vlz þ ðgl þ glt Þ þ þ 3 or2 r or r2 oroz oðgl þ glt Þ ovlr vlr 1 ovlr vlr ovlz þ2 þ þ þ ; or or r 3 or r oz ð2Þ ovlz ovl ovl þ vlr z þ vlz z ot or oz 2 l opl o vz 1 ovlz o2 vlz ¼ þ ðgl þ glt Þ þ þ oz or2 r or oz2 2 l 1 o vr 1 ovlr o2 vlz þ 2 þ ðgl þ glt Þ þ 3 oroz r oz oz oðgl þ glt Þ ovlz ovlr oðgl þ glt Þ þ þ þ2 or oz or oz l ovz 1 ovlr vlr ovlz þ þ oz 3 or r oz l Þ; þ ql0 blt gðT l T0l Þ ql0 blC gðCAl CA;0
ql
ql
319
oCAl oC l oC l þ vlr A þ vlz A ot or oz 1 o gl oCAl oT l rql Dl þ t þ rql dl ¼ r or rc or or l l l o l gt oCA l l l oT q D þ þqd þ ; oz rc oz oz
ð5Þ
ok l ok l ok l þ vlz þ vlr ot or oz l 1 o gt ok l o glt ok l l l r g þ g þ ¼ þ r or oz rk or rk oz þ P þ G ql el ;
ð6Þ
l l oel l oe l oe q þ vz þ vr ot or oz 1 o glt oel o glt oel l l r g þ g þ ¼ þ r or oz re or re oz l
2
þ c1e f1 ðP þ GÞ ql c2e f2 ð3Þ
l oT oT l oT l ðql clp blT pl Þ þ vlz þ vlr ot r oz l l 1 o g oT oC l r Kl þ t þ rBl A ¼ r or rT or or l l o g oT oC l Kl þ t þ Bl A þ oz rT oz oz l l oCA oC oC l blC pl þ vlr A þ vlz A ot or oz l l l oH g oCA oT l ql þ dl þ Dl þ t or rc or or l l l oH g oCA oT l ql þ dl þ Dl þ t oz rc oz oz l l l l ovr 1 ovr vr ovlz ovr l l þ 2ðg þ gt Þ þ þ or 3 or r oz or l l l l ovr ovz ovr ovz þ ðgl þ glt Þ þ þ oz or oz or l l l l ovz 1 ovr vr ovlz ovz l l ; þ 2ðg þ gt Þ þ þ oz 3 or r oz oz ð4Þ
el : kl
ð7Þ
The meaning of the symbols used in the above equations is as follows. The terms P and G represent the production of turbulent kinetic energy and we have (see [3]) l l ovr 1 ovlr vlr ovlz ovr P ¼ 2ðgl þ glt Þ þ þ or 3 or r oz or l l l l ovr ovz ovr ovz þ ðgl þ glt Þ þ þ oz or oz or l l l ovz 1 ovr vlr ovlz ovz ; þ 2ðgl þ glt Þ þ þ oz 3 or r oz oz ð8Þ oT l oC l þ ql0 blC g A : ð9Þ oz oz Further, v is the velocity, bT and bC the thermal and the solutal expansion coefficient, respectively, T is the temperature, CA the concentration of the component A of the mixture, q is the density, p the pressure, g the dynamic viscosity and gt the turbulent dynamic viscosity, g ¼ j~ gj is the gravity, K is the thermal conductivity, cp the specific heat, B the Dufour coefficient, H the difference in partial G ¼ ql0 blT g
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
320
specific enthalpies of the components of the mixture, D is the diffusion coefficient, d is the Soret coefficient, rc is the turbulent Schmidt number, rT the turbulent Prandtl number, k the turbulent kinetic energy, and e means the dissipation rate of k. The superscript l denotes the liquid phase, q0 ; T0 ; CA;0 are values corresponding to a reference state. The turbulent dynamic (eddy) viscosity glt is defined as 2
glt ¼ cl fl
kl el
ð10Þ
in the framework of the low-Reynolds-number k–e models, where cl is an empirical constant and fl is a function depending on the chosen model. The functions f1 , f2 and fl in the models under consideration are given as follows. Lam–BremhorstÕs model [10]: 3 0:05 f1 ¼ 1:0 þ ; f2 ¼ 1:0 expðRe2s Þ; fl 20:5 2 fl ¼ ð1:0 expð 0:0165Red ÞÞ 1:0 þ : Res DavidsonÕs model [11]: 3 0:14 f1 ¼ 1:0 þ ; fl f2 ¼ ½1 0:27 expððRes Þ2 Þ ½1 expðRed Þ ; " # 3:4 fl ¼ exp : 2 ð1 þ 0:02Res Þ
re ¼ 1:3;
c1e ¼ 1:44; rc ¼ 0:9;
c2e ¼ 1:92;
rT ¼ 0:9;
rk ¼ 1:0; pffiffiffiffi d Red ¼ k l l ; g
2
kl Res ¼ l l ; ge
where vamp is the velocity of motion of the ampoule, taken as positive for the downward direction, the other quantities being denoted in the same way as in the liquid, and the superscript s denotes the solid phase. At the solid/liquid interface assumed to have the form z ¼ Zðr; tÞ, the following balance equations can be written: ð16Þ
ql vlr viN ¼ ql vlr ðnr vlr þ nz vlz Þ þ nr ðpl p0 Þ
l ovlr ovr ovlz nz ðgl þ glt Þ þ or oz or l l l 2 ovr vr ovz þ ðgl þ glt Þnr þ þ ; ð17Þ 3 or r oz
2nr ðgl þ glt Þ
ql vlz viN þ qs vamp viN
where d is the distance to the nearest wall. This completes the description of the model equations in the liquid phase. In the solid we have [1] vsr ¼ 0;
ð12Þ
ð13Þ ps ¼ p0 ¼ const:; s oT oT s qs csp vamp ot oz 1 o oT s oC s o oT s rK s Ks ¼ þ rBs A þ r or oz or or oz s s oC oCA oC s þ B s A qs H s vamp A oz ot oz s s s oH oC oT qs Ds A þ ds þ or or or s s s oH s oT s s oCA q D þ þd ; ð14Þ oz oz oz s oCA oC s qs vamp A ot oz 1 o oC s oT s rqs Ds A þ rqs ds ¼ r or or or s o oC oT s qs Ds A þ qs ds þ ; ð15Þ oz oz oz
ðql qs ÞviN ¼ ql ðnr vlr þ nz vlz Þ þ nz qs vamp ;
Some of the constants and functions are common for both the models: cl ¼ 0:09;
vsz ¼ vamp ¼ const:;
ð11Þ
¼ ql vlz ðnr vlr þ nz vlz Þ qs nz v2amp þ nz ðpl p0 Þ l l ovr ovlz l l ovz l l nr ðg þ gt Þ 2nz ðg þ gt Þ þ oz oz or l l l 2 ovr vr ovz þ ðgl þ glt Þnz þ þ ; ð18Þ 3 or r oz
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
the solid into the liquid, LðT0 Þ is the latent heat of melting at the reference temperature T0 , and p0 is a reference pressure. The above equations can be developed analogously as with the previous laminar model [1], however, the turbulent processes have been taken into account here. Besides the balance conditions, two additional conditions at the solid/liquid interface are formulated. They are given by the liquidus and solidus curves of the phase diagram,
½ql LðT0 Þ þ ðql clp qs csp ÞðT T0 Þ þ ðql qs Þcsp T0 pl þ p0 viN ¼ ql LðT0 Þðnr vlr þ nz vlz Þ þ ½ql ðnr vlr þ nz vlz Þ þ qs nz vamp csp T0 þ ½ql clp ðnr vlr þ nz vlz Þ þ qs csp nz vamp ðT T0 Þ pl ðnr vlr þ nz vlz Þ p0 nz vamp l glt oT l l l oCA þB nr K þ rT or or gl oT l oC l þ Bl A nz K l þ t rT oz oz oT s oC s þ nr K s þ Bs A or or s s s oT s oCA þ nz K þB oz oz l l gt oCAl l oT l l l þd D þ nr q H rc or or l l l gt oCA l oT l l l þd D þ nz q H rc oz oz s s oC oT þ nr q s H s D s A þ ds or or oC s oT s þ nz qs H s Ds A þ ds ; oz oz
321
T l ¼ Tm ðCAl Þ;
ð21Þ
T s ¼ Tc ðCAs Þ;
ð22Þ
where Tm ðCAl Þ is the liquidus curve and Tc ðCAs Þ is the solidus curve. Finally, we formulate the local thermodynamic equilibrium condition in the form Tl ¼ Ts
ð19Þ
ðql CAl qs CAs ÞviN ¼ ql CAl ðnr vlr þ nz vlz Þ þ qs CAs nz vamp gl oCAl oT l nr ql Dl þ t þ dl rc or or l l g oCA oT l þ dl nz ql Dl þ t rc oz oz s s oC oT þ nr qs Ds A þ ds or or oC s oT s ; þ nz qs Ds A þ ds oz oz ð20Þ
where viN is the velocity of the interface in the normal direction, nr , nz are the components of the unit normal vector to the interface, directed from
ð23Þ
at the phase interface. Appropriate initial and boundary conditions corresponding to the particular crystal growth techniques are formulated in a common way corresponding to the particular experimental setup. The conditions for T, C, and v are analogical to those used in the previous laminar model [1]. We employ constant initial conditions for k and e and Dirichlet boundary conditions for these variables (these will be discussed later, see also [8]). The boundary conditions at the axis of symmetry are treated in the standard way. If we consider a system cooled from bottom, then the boundary and initial conditions can be written in the following manner: vlr ðr; L; tÞ ¼ vlz ðr; L; tÞ ¼ vamp ;
r 2 ½0; R ;
vlr ðR; z; tÞ ¼ vlz ðR; z; tÞ ¼ vamp ;
z 2 ½ZðR; tÞ; L ;
ð24Þ
ð25Þ vlr ð0; z; tÞ ¼ vamp ; ovlz ð0; z; tÞ ¼ 0; or
z 2 ½Zð0; tÞ; L ;
z 2 ½Zð0; tÞ; L ;
T l ðR; z; tÞ ¼ T2 ðz; tÞ > Tm ðCAl Þ;
ð26Þ ð27Þ
z 2 ½z2 ðtÞ; L ; ð28Þ
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
322
ok l ð0; z; tÞ ¼ 0; or
z 2 ½Zð0; tÞ; L ;
ð29Þ
oel ð0; z; tÞ ¼ 0; or
z 2 ½Zð0; tÞ; L ;
ð30Þ
k l ðR; z; tÞ ¼ k1 ¼ const:;
z 2 ½ZðR; tÞ; L ;
ð31Þ
e ðR; z; tÞ ¼ e1 ¼ const:;
z 2 ½ZðR; tÞ; L ;
ð32Þ
k l ðr; Z; tÞ ¼ k1 ¼ const:;
r 2 ½0; R ;
ð33Þ
el ðr; Z; tÞ ¼ e1 ¼ const:;
r 2 ½0; R ;
ð34Þ
l
l
k ðr; Zðr; tÞ; tÞ ¼ k1 ¼ const:; l
e ðr; Zðr; tÞ; tÞ ¼ e1 ¼ const:;
r 2 ½0; R ;
ð35Þ
r 2 ½0; R ;
ð36Þ
z2;0 z þ T2 ðz; 0Þ z2;0 z1;0 z z1;0 ; r 2 ½0; R ; z 2 ½z1;0 ; z2;0 ; z2;0 z1;0
T s;l ðr; z; 0Þ ¼ T1 ðz; 0Þ
ð49Þ T s ðr; z; 0Þ ¼ T1 ðz; 0Þ;
r 2 ½0; R ; z 2 ½0; z1;0 ; ð50Þ
l CAl ðr; z; 0Þ ¼ C1;0 ðz; 0Þ;
r 2 ½0; R ; z 2 ½Z0 ðrÞ; L ; ð51Þ
s ðz; 0Þ; CAs ðr; z; 0Þ ¼ C1;0
r 2 ½0; R ; z 2 ½0; Z0 ðrÞ ; ð52Þ
k l ðr; z; 0Þ ¼ const:;
z 2 ½z2 ðtÞ; L ; r 2 ½0; R ;
s;l
oT ðR; z; tÞ ¼ 0; or
z 2 ½z1 ðtÞ; z2 ðtÞ ;
T s ðR; z; tÞ ¼ T1 ðz; tÞ < Tc ðCAs Þ;
z 2 ½0; z1 ðtÞ ;
ð38Þ
oT l ðr; L; tÞ ¼ 0; oz
r 2 ½0; R ;
ð39Þ
oT s ðr; 0; tÞ ¼ 0; oz
r 2 ½0; R ;
ð40Þ
oT s;l ð0; z; tÞ ¼ 0; or
el ðr; z; 0Þ ¼ const:;
z 2 ½z2 ðtÞ; L ; r 2 ½0; R ;
ð54Þ
where z1 ðtÞ ¼ z1;0 þ vT t, z2 ðtÞ ¼ z2;0 þ vT t, z1;0 6 z2;0 , and vT is the velocity of motion of the facilities asserting the boundary conditions for temperature.
3. Computer implementation z 2 ½0; L ;
CAl ðr; L; tÞ ¼ C1 ¼ const:;
r 2 ½0; R ;
ð41Þ ð42Þ
oCAs ðr; 0; tÞ ¼ 0; oz
r 2 ½0; R ;
ð43Þ
oCAs;l ð0; z; tÞ ¼ 0; or
z 2 ½0; L ;
ð44Þ
oCAs;l ðR; z; tÞ ¼ 0; or
z 2 ½0; L ;
ð45Þ
Zðr; 0Þ ¼ Z0 ðrÞ;
r 2 ½0; R ;
ð46Þ
vlr ðr; z; 0Þ ¼ vlz ðr; z; 0Þ ¼ vamp ; r 2 ½0; R ; z 2 ½0; L ; ð47Þ l
ð53Þ
ð37Þ
T ðr; z; 0Þ ¼ T2 ðz; 0Þ;
r 2 ½0; R ; z 2 ½z2;0 ; L ; ð48Þ
Our approach here is analogous to that employed with the previous laminar model, the principal difference being the necessity to process the equations for k and e in addition. Therefore, we are rather brief in this Section referring the erny et al. [1] interested reader to the paper by C where all the delicate details of the discretization procedure have been discussed and necessary references quoted. To begin the survey of our computational procedure we note that we use the penalty method for the continuity Eq. (1), which means that a slightly modified form of this equation is used in the computer implementation, where an artificial term ep pl with a small ep is added (see [12] in this respect). This makes it possible to eliminate the pressure from the system of equations to be solved in the liquid and to obtain a system for velocities, k and e only. The motivation here is the same as in
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
[1]: the elimination of pressure from the unknown variables not only makes an appropriate discretization feasible but also reduces the computing time (we are not interested in the pressure field in our practical simulations). For the discussion of the application of the penalty method the interested reader is referred, e.g., to [13,14]. The system of equations obtained thus in the liquid and that in the solid are discretized by using the Galerkin finite element method with bilinear 4 node quadrilateral elements in space, linear basis functions in time. The space discretization leads to two systems of first-order differential equations, one of them corresponding to the solid, the other one to the liquid phase. These are then subject to the time discretization. As the overall procedure has been described and critically discussed in [1] we refer the reader there for more details again. However, we feel necessary to repeat that the practical situations for which our model is intended are strongly diffusion dominated and that the Galerkin FEM does not suffer from spurious node-to-node oscillations here. Finally, the matrices and vectors of the solid and liquid phases are coupled together to obtain the global system of algebraic equations that represents the whole space region ½0; R ½0; L and describes the change of the system variables from time t to t þ Dt in our discrete model. In this way we get (similarly as in [1]) ½HP fwgtþDt ¼ ½PF fwgt þ fFg;
ð55Þ
where D s s s ; . . . ; Tm1 ; CA;m1 ; fwg ¼ T1s ; CA;1 l vlr;1 ; vlz;1 ; T1l ; CA;1 ; k1l ; el1 ; . . . ; ET l vlr;n ; vlz;n ; Tnl ; CA;n ; knl ; eln :
ð56Þ
Here, the subscripts 1; . . . ; m 1 and 1; . . . ; n are the indices of the nodal points, and [HP], [PF], {F} are the resulting global matrices and the global vector which correspond to both the solid and liquid phases and are localized in accordance with the ordering of the global vector fwg. In the combination of matrices and vectors we employ the local thermodynamic equilibrium condition (23) and the FEM analogs of the interface bal-
323
ances of internal energy and mass of the component A so that we obtain nonzero components of the right-hand side vector {F} at the positions corresponding to these equations. The system (55) is doubly nonlinear as its terms [HP], [PF], and [F] are not only functions of fwg but also of the a priori unknown position Z of the phase interface, which is also to be determined in the course of the solution. The moving boundary problem (55), (21) and (22) is solved now by the successive approximation method and a front-fixing method that transforms the two-dimensional variable space regions ½0; R ½0; Zðr; tÞ , ½0; R ½Zðr; tÞ; L corresponding to the solid and liquid phases, respectively, to the fixed space domains ½0; 1 ½0; 1 . Due to the character of the modeled crystal growth techniques and the assumed form of the phase interface it is possible to use analogous transformations here as are those for the one-dimensional problems introduced by Landau [15]. Hence, we put f¼
r ; R
g¼
z ; Zðr; tÞ
ð57Þ
in the solid phase and f¼
r ; R
g¼
z Zðr; tÞ ; L Zðr; tÞ
ð58Þ
in the liquid phase, transform the problem and work with the independent variables f; g instead of r; z. The transformations are applied to all the element matrices and vectors defined before, however, the system of equations (55) remains formally the same and has similar nonlinear properties after the transformation. The resulting system of nonlinear algebraic equations is then solved at each time step by two nested iteration procedures using the successive approximation method with underrelaxation. The inner iterations are due to the nonlinearity of the field equations, the outer iterations are performed with the aim to satisfy all the conditions prescribed at the a priori unknown phase interface for the particular time instants. Our numerical tests have shown that only one possible choice of the outer iterative scheme leads to a stable computational algorithm. This scheme consists in applying the phase-diagram
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
324
and the local equilibrium conditions (21)–(23) as Dirichlet conditions at the interface and using the interface internal energy balance equation as the convergence criterion. The situation here is similar to that in our previous laminar model. Therefore, also the overall scheme of one time step of our algorithm with two nested iteration procedures is analogous but we feel it reasonable to present it here for the convenience of the reader: 1. Put Z ð1Þ ðt þ DtÞ ¼ ZðtÞ, Z_ ð1Þ ðt þ DtÞ ¼ Z_ ðtÞ, ð1Þ fwgtþDt ¼ fwgt . ð2Þ 2. Compute fwgtþDt from (55) using Z ð1Þ ðt þ DtÞ, ð1Þ Z_ ð1Þ ðt þ DtÞ and fwgtþDt in the matrices [HP], [PF]. 3. If vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPn ð2Þ 2 2 u i¼1 ðwi Þ ðwð1Þ Þ i t > di ; Pn ð2Þ 2 i¼1 ðwi Þ where di is a user-defined error tolerance for the internal iteration cycle, put ð1Þ
ð2Þ
fwgtþDt ¼ fwgtþDt ; and go back to Step 2. Otherwise, put ð1Þ
ð2Þ
fwgtþDt ¼ fwgtþDt and continue with Step 4. 4. Compute Z_ ð2Þ ðt þ DtÞ using the finite element analog of the interface energy balance equation ð1Þ (19) and fwgtþDt . 5. If Z_ ð2Þ ðt þ DtÞ Z_ ð1Þ ðt þ DtÞ > de ; Z_ ð1Þ ðt þ DtÞ where de is a user-defined error tolerance for the external iteration cycle, then choose a 2 ½0; 1 , put Z_ ð1Þ ðt þ DtÞ ¼ aZ_ ð2Þ ðt þ DtÞ þ ð1 aÞZ_ ð1Þ ðt þ DtÞ; Z ð1Þ ðt þ DtÞ ¼ ZðtÞ þ DtZ_ ð1Þ ðt þ DtÞ; and go back to Step 2. Otherwise, put ð1Þ
fwgtþDt ¼ fwgtþDt ;
Z_ ðt þ DtÞ ¼ Z_ ð2Þ ðt þ DtÞ; 1 Zðt þ DtÞ ¼ ZðtÞ þ Dt½Z_ ð2Þ ðt þ DtÞ þ Z_ ðtÞ 2 and go to the next time step. The computer implementation of the numerical algorithm was performed in Fortran 77 and a computer code has been developed. The numerical experiments discussed in the next Section were done on on a PC with two processors PII 450 MHz under OS Linux and on a workstation HP C180 under OS HP–UX.
4. Numerical simulations The main object of the computer simulations we performed has been to compare the results produced by the model presented in this paper with those obtained using our previous model of CdZnTe crystal growth with laminar melt flow from [1]. As before in [1], a set of numerical calculations was performed first in order to optimize the numerical parameters of the model. We focused namely on the type of the iteration procedure, the magnitude of the pressure penalty ep , the initial values of k and e, the time step of simulation and the optimal values of the relaxation coefficients for particular variables, which strongly affect the rate of convergence. From these extensive numerical tests the following optimum parameters were obtained: ep in the range ½1010 ; 108 kg1 m2 s2 , the initial values k 102 m2 s2 and e in the range ½104 ; 103 m2 s3 , and the relaxation coefficients for velocities about 0.6, for temperature and concentration about 0.8 and for turbulent kinetic energy and its dissipation rate about 0.3. We observed that particularly the initial values of k and e had a significant influence on the stability of the algorithm. Also for the values of the pressure penalty which laid outside the above mentioned optimum interval a stable solution was not achieved. A typical number of inner iteration steps for achieving the convergence of all variables with di ¼ de ¼ 103 was of the order of 102 and the
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
number of outer iterations of the order of 101 . The optimum finite element mesh in both the liquid and the solid that guaranteed sufficient precision in a still acceptable CPU time was 20 20. The time step was controlled by limiting its maximum value by Dtmax , the optimum value being Dtmax 1 s. After finding the optimal numerical parameters we performed a set of numerical simulations with both the laminar and the turbulent models corresponding to a real experimental situation. All our numerical computations modeled the CdZnTe crystal growth using the vertical Bridgman method. The common basic parameters of the experimental setup were as follows: the radius of the ampoule R ¼ 37:5 mm, the length of the ampoule L ¼ 100 mm, the temperatures of the Bridgman furnaces 1125 and 975 °C, the distance between the furnaces 10 mm, the velocity of the ampoule vamp ¼ 106 m s1 or vamp ¼ 107 m s1 . The thermophysical properties of CdZnTe used in our simulations can be found in Table 1 of [1]. Fig. 2 shows the development of the shape of the solid/liquid interface during the CdZnTe crystal growth in the turbulent model. We can observe that an S-shape is characteristic for the initial time period studied, whereas a C-shape is established after achieving a more developed state. A comparison with the results obtained with the laminar model from [1] shows significant differ-
Fig. 2. Development of the phase interface in time for vamp ¼ 106 m s1 .
325
Fig. 3. Comparison of the position of the phase interface between the turbulent and the laminar model for vamp ¼ 106 m s1 .
ences in the case of vamp ¼ 106 m s1 (see Fig. 3). For the case of vamp ¼ 107 m s1 no differences were found in the whole range of the simulated time interval (10,000 s). Reverse movement (remelting) of the phase interface near the wall was observed with the turbulent model that can be seen in Fig. 2. Fig. 4 shows the evolution of the velocity fields characteristic for our turbulent model. The typical time stable flow pattern for the laminar model is shown in Fig. 5. While the laminar flow model always leads to a very simple flow structure with only one recirculating region, the structure of the flow calculated by the turbulent model is more complicated. Isolated recirculating flow regions arise in the domain in the vicinity of the phase interface, while on the opposite side practically no flow is observed. Especially in the vicinity of the wall and the phase interface, instabilities were observed. Also, the magnitude of the velocities was about one order lower for the turbulent flow model compared to both the velocity of the interface and the velocities in the laminar model. The computed concentration profiles are different for both the turbulent and the laminar model (see Fig. 6). The turbulent profiles are affected by the more intensive transport processes.
326
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
Fig. 6. Comparison of the concentration profiles at the axis of the ampoule for vamp ¼ 106 m s1 at time t ¼ 10 000 s.
Fig. 4. Time series of the turbulent flow pattern for vamp ¼ 106 m s1 .
Moreover, some concentration peaks were found in the produced crystal, in particular near the wall. These peaks were more pronounced in the case of the turbulent model, where also variable concentration along the phase interface was found. Finally, we note that we also compared the two above mentioned turbulent models, namely that of Davidson and that of Lam and Bremhorst. No significant differences in the structure of the velocity field or of the shape of the solid/liquid interface were found. Sometimes we had problems with obtaining negative values of the turbulent kinetic energy and its dissipation rate using the Lam–Bremhorst model. Therefore, DavidsonÕs model has been used for our practical simulations.
5. Discussion
Fig. 5. Stable structure of the laminar flow pattern for vamp ¼ 106 m s1 .
In our numerical simulations, we tested the influence of turbulent natural convection due to the thermal and solutal expansion in the liquid on the temperature, concentration and velocity fields, and on the position and velocity of the moving solid/ liquid interface. We believe that the results obtained have shown that the proposed computational model involving
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
turbulent motion in the melt gives more realistic simulation than the previous laminar model. The turbulent model simulates the whole cascade of kinetic energy, which is successively transferred from the largest eddies that are on the order of the flow domain to the smallest eddies (measuring a few millimeters or less), where the turbulent kinetic energy is converted into the internal heat of the fluid by viscous dissipation. The laminar model is not able to describe such processes in a proper way due to its disability to simulate processes running in smaller scales than those corresponding to the size of the computational mesh used. However, the small dissipative eddies play a significant role in the transport processes occurring in the system and hence they have an important impact on the whole process of crystallization. We have found distinct differences in the velocity patterns and the profiles of concentration between the two models. The turbulent melt flow was characterized by nonstationary oscillations contrary to the stable laminar pattern with two inversely oriented vortices in the ampoule. Instabilities of the turbulent flow pattern were observed close to the wall. This fact results in re-melting of the created solid phase, the oscillation of the phase interface and consequential appearance of defects in the crystals grown. This result is in agreement with the experimental observations of enhanced occurrence of defects near the wall. This phenomenon, which can explain some perturbations during the crystallization, was not (and could not be) observed with the laminar model. The results obtained with the turbulent model indicate that the flow of the melt is mainly affected by the movement of the phase interface and that the buoyancy effects play a middling role only. The buoyancy effects tend to establish a stable structure of the flow pattern, usually two vortices in the domain as we observed with our laminar model, but in the case of the turbulent fluid flow this tendency is dumped by additional viscous dissipation. Therefore, no stable structures have been observed in the turbulent flow field. The turbulent concentration profiles showed marked changes in both time and space. Peaks in the concentration field were observed in the produced solid phase. The observed instabilities near
327
the wall are in agreement with the experimental experience. From this point of view it seems to be clear that the turbulent model should be used to obtain a realistic representation of the events being in progress during the crystal growth. In the future work, we intend to perform a critical comparison with the real experimental data with the aim to assess the importance of particular factors influencing the process of crystallization and to provide more realistic values of important constants (rT ; rc ) used in simulating the CdZnTe material.
6. Conclusions The proposed turbulent model has shown a new possibility how to describe the more complicated transport processes that play a significant role in the production of CdZnTe single crystals. The turbulent model has revealed new facts about the mutual relations between the particular factors involved in the process of crystal growth, such as between the buoyancy effects and the interface kinetics, or about the diffusion of a compound of the mixture that have a pronounced influence on the final result of the production process. We conclude that the turbulent model is more appropriate to describe the transport processes in the melt during crystallization than the laminar model.
Acknowledgements This research has been supported by the Grant Agency of the Czech Republic, under grants # 202/ 99/1646 and # 201/01/1200.
References erny, A. Kalbac, P. Prikryl, Comp. Mater. Sci. 17 [1] R. C (2000) 34. [2] J.O. Hinze, Turbulence: An Introduction to Its Mechanism and Theory, McGraw-Hill, Delft, 1959. [3] W. Frost, T.H. Moulden, Handbook of Turbulence, Plenum Press, New York, 1992.
328
erny et al. / Computational Materials Science 25 (2002) 316–328 R. C
[4] M. Griebel, T. Dornseifer, T. Neunhoeffer, Numerical Simulation in Fluid Dynamics: A Practical Introduction, SIAM, Philadelphia, 1998. [5] W. Rodi, Turbulence Models and Their Application for Hydraulics, International Association for Hydraulics Research, Delft, 1980. [6] M. Lesieur, Turbulence in Fluids, Kluwer, Dordrecht, 1997. [7] B. Launder, D. Spalding, Mathematical Models of Turbulence, Academic Press, New York, 1972. [8] B. Mohammadi, O. Pironneau, Analysis of the K-Epsilon Turbulence Model, Wiley, Chichester, 1993.
[9] W.C. Patel, W. Rodi, G. Scheuerer, AIAA J. 23 (1985) 1308. [10] C.K.G. Lam, K.A. Bremhorst, ASME J. Fluids Eng. 103 (1981) 456. [11] L. Davidson, Numer. Heat Transfer 18 (1990) 129. [12] F. Thomasset, Implementation of Finite Element Methods for Navier–Stokes Equations, Springer, New York, 1981. [13] T.J.R. Hughes, W.T. Liu, A.J. Brooks, J. Comput. Phys. 30 (1979) 1. [14] M. Bercovier, M. Engelman, J. Comput. Phys. 30 (1979) 181. [15] H.G. Landau, Quart. Appl. Math. 8 (1950) 81.