Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://mgumus.chem.msu.ru/seminar/mat/Paper-Crystallization-Eng.pdf Äàòà èçìåíåíèÿ: Tue Apr 9 14:57:10 2013 Äàòà èíäåêñèðîâàíèÿ: Thu Feb 27 21:09:06 2014 Êîäèðîâêà:

Applied Mathematics, 2010, 1, 159-178
doi:10.4236/am.2010.13021 Published Online September 2010 (http://www.SciRP.org/journal/am)

Solidification and Structuresation of Instability Zones
1

Evgeniy Alexseevich Lukashov1, Evgeniy Vladimirovich Radkevich2
Soyuz Aircraft-Engine Scientific and Technical Complex (AMNTK SOUYUZ) Moscow 119270, Russia 2 Moscow State University Moscow 119991, Russia E-mail: elukashov@yandex.ru, evrad07@gmail.ru Received May 1, 2010; revised July 15, 2010; accepted July 17, 2010

Abstract
Two mathematical crystallization models describing structure formations in instability zones are proposed and justified. The first model, based on a phase field system, describes crystallization processes in binary alloys. The second model, based on a modified Biot model of a porous medium and the convective Cahn­ Hilliard model, governs oriented crystallization. Physical interpretation and numerical analysis are discussed. Keywords: Solidification, Structuresation of Instability Zones, Phase Field Model

1. Introduction
Unlike the main properties of oriented crystallization, properties responsible for the alloy structure have not yet been studied well. At the same time, owing to recent experimental results, many details of crystallization become known. In this paper, we propose the so-called "reconstruction" of oriented crystallization processes, i.e., a detailed theoretical description based on the known main properties. To reconstruct a process of binary alloy crystallization, one should begin with the question why the process "can live" in the stochastic instability. Perhaps, like in the case of complicated systems [1], the crystallization process can exist for a long time only due to solid structure formations in instability zones. Moreover, taking into account such structure formations, we are able to explain the solid phase growth ­ the crystallization mechanism. It is known [2] that the structure formation in an alloy obtained by the oriented crystallization method is characterized by the following properties. 1) The process proceeds in a solid­liquid domain ­ a dynamic porous medium­where the solid phase is represented by growing dendrites, whereas the liquid phase occupies the space between these dendrites. According to experimental results, the solid phase growth is of order O t , where t is time. 2) In the case of overlapping dendrites (in particular, their secondary branches), the melt solidification can lead to the contraction of melt and formation of internal stresses and micropores. 3) In turn, a solid-liquid crystallization zone appears



because of the instability of the crystallization front which can be caused by the following reasons: concentration overcooling, segregation of the melt components in view of the spinodal decomposition (i.e., phase transition with instable states) when the melt deeply penetrates into the metastable (or even labile) domain under high-speed (high-gradient) cooling in the interphase zone. 4) Properties of a new alloy are encoded in a seed crystal (a small piece of the solid phase) which, like the genetic code, determines the required properties of the crystallized part. The experimental results concerning the distribution of crystallization centers over the blank surface are represented in Figure 1, where it is seen that crystallization centers are concentrated on convex parts of the surface, but not on its concave parts. In both cases, one of the phases grows in time, whereas the other decreases. We also note that the picture demonstrates the structure ordering. The goal of this paper is to construct mathematical models reflecting Properties 1­4 and simulating the structure formation in alloys and, first of all, in instability zones. We propose two models (cf. Sections 3 and 4) with banding structure in the zone of instability. But, first, we emphasize that, within the frameworks of models where structure formations in the instability zone are not taken into account, it is impossible to obtain the experimental order O t of the solid phase growth (cf. Property 1). We illustrate this fact by considering the well-known statistical Kolmogorov model [3] describing



Copyright © 2010 SciRes.

AM


160

E. A. LUKASHOV ET AL.

Figure 1. The experimental results are represented the distribution of crystallization centers over the blank surface.

the process of metal crystallization (cf. Section 2).

2. Kolmogorov's Model of Metal Crystallization
2.1. Physical Interpretation
In metallurgy, it is important to know the crystal growth velocity under a random formation of crystallization centers. Under rather general assumptions, Kolmogorov [3] derived an expression for the probability p(t) that a randomly taken point P gets into the crystallized mass during the crystallization time-interval. With rather good approximation, we can assume that the mass crystallized in time t is equal to p(t). Then it is possible to find the number of crystallization centers formed during the whole process of crystallizaton.

whereas the probability of appearing more than one crystallization centers is of order o( t ), where o( t ) is infinitesimal in comparison with t . These probabilities are independent of the distribution of crystallization centers that are formed before time t (the process is Markovian) if only the freedom of V' from the crystallized mass at time t is not guaranteed. 2) Around the new-formed crystallization centers and around the entire crystallized mass, the mass grows with linear velocity c t, n k t c n depending on time t and direction n, n = 1. It is assumed that the endpoints of vectors c(n)n started at the origin and directed towards n form a convex surface. Note that the homogeneous dependence of the linear velocity c(t, n) on the direction n at all points is an essential restriction. In other words, we obtain formulas that are valid either in the case where the growth is uniform along all directions, or in the case of crystals of arbitrary shape but with the same spatial orientation. We also mention the case where all crystallization centers are formed at initial times, in mean, per volume unit. We obtain the corresponding formulas by taking into account that, in this case, t is the Dirac function 0 concentrated at the origin.

2.2. Mathematical Statement
Consider a domain V d , d = 2,3. Assume that at the initial time t = 0, the domain V is occupied by the so-called mother phase. At time t, some part V1(t) of V is occupied by a crystallized matter. Moreover, V1(t) enlarges in t as follows. 1) In a free part V/V1 of V, new crystallization centers appear, so that for any domain V V / V1 the probability of appearing a single crystallization center in V' during time t is equal to

2.3. Formulas
We introduce the mean (over all directions) velocity of
AM

t V t o t ,
Copyright © 2010 SciRes.


E. A. LUKASHOV ET AL.

161

the growth of crystallized mass c by the formula
cd = 1d c (n)ds, d = 2, 3, |S | s

where the integral is taken over the surface of unit sphere S in d with center at the origin, S 4 if d = 3 and S 2 if d = 2. Then the following assertions hold. 1) For a sufficiently large (in comparison with the size of a crystallization center) domain V the domain V1(t) occupied by the crystallized phase takes the form
V1 t V 1 e

structure formation in the zone of instability. A crystallization model based on the phase field conception was constructed in [6], where, in particular, a sawtooth solution to the temperature distribution problem in the phase transition domain was obtained. This result agrees with the qualitative description of autocrystallization phenomena in [7,8]. The goal of this section is to obtain a sawtooth solution to the temperature distribution problem for the following phase field system:
, x, t Q , t t



d Ad cd d



(1)

(7) (8) , (9) (10)

If (t) and c(t, n) are time-independent, we can set (t) = , k(t) = 1. In this case,
d




2

t

d 1

= 2 3 , t
0

d 1

,

(2)

t 0

x, ,

t 0



0

x,

which implies

| = 1, | = b ,
(3)

V1 t = V 1 e

A d c d t d 1 d 1 d

2) If all crystallization centers are formed at initial times, then
t t d t k d dt k d 0 t 0
t d d

(4)

If, in addition, k = 1, i.e., c (t, n) is independent of t, then
d t d ,

where is the temperature; is the specific concentration of the order function, equal to 1 in the liquid phase and to 1 in the solid phase; = const ; Q = (0, T) â , where n is a bounded domain with C - boundary, n 3; 0, T ; the functions 0 and 0 are sufficiently smooth for const > 0, and the function b is also sufficiently smooth. The system (7)­(10) describes slow crystallization processes [9] with an instable domain of intermediate aggregate state, where a structure formation appears.

(5)
d Ad cd t d

which implies

3.1. Wave Train Type Solutions and Singular Limit Problem
Here, we consider a more general case where 0 BV , but 0 B VC 1 (cf. [6]). In the case of diffusion, we say that is a domain of intermediate aggregate state (an IAS-domain) if 0 ( x, ) 0 weakly as 0 in some subdomain 0 of nonzero measure. cr In accordance with [6], an IAS-domain is formed by a large number M of domains of pure (solid and liquid) phases of small volume of order (i.e., M M and 0 as 0 ). The macroscopic description of an IAS-domain can be obtained by computing the weak limit of wave train type solutions as 0 . We formulate conditions imposed on IAS-domains. 1) The weak limit of the order functions ( x, t , ) as 0 vanishes identically in the transition zone * . t 2) In the domain *, corresponding to the regulat rization of the IAS-domain, the solution to the phase field system can be described in terms of the wave train
1

V1 t V 1 e





(6)

We see that the mass growth is of power-like order O t , = 1, 2, 3, d = 1, 2, 3.



2.4. Conclusions
The Kolmogorov model is not suitable for describing crystallization of twocomponent mixtures. Indeed, within the frameworks of the Kolmogorov model, the fact that the mass growth is of power order implies that the velocity is finite at t = 0, which contradicts the initial stage of the spinodal decomposition generating an initial distribution of crystallization centers.

3. Model of Binary Alloy Crystallization
Based on the phase field system proposed in [4] and [5], we construct a model of binary alloy crystallization with
Copyright © 2010 SciRes.

A function

0

belongs to the class BVC if
0 0

0

is a function of

bounded variation ( BV ) and

= 1.

AM


162

E. A. LUKASHOV ET AL.

structure. In this case, the domain is divided into a large number of domains of "small" volume occupied by "pure phases" and transition zones between them. Remark 3.1 Condition (1) means, in particular, that for almost all t the limit order function belongs to BV () , but BVC () . Condition (2) is based on the conception proposed in [6]. According to this conception, the wave train structure is described by a chain of modified Stefan problems in domains of "small" volume occupied by "pure phases" and can be used for approximating the temperature in an IAS-domain. Such a structure is called the diffusion of the IAS-domain. Remark 3.2 A situation where the limit order function vanishes on a set of nonzero measure is not good since this case corresponds to instable solutions to the isothermal diffusion equation. It is clear that such solutions can exist only under rather special conditions. Therefore, we need to impose rather restrictive conditions on the geometry of domains , * , as well as t on the initial and boundary conditions. Remark 3.3 From the point of view of the theory of distributions, free boundary problems are problems about singularity propagation. Indeed, in the rigid-front situation, the limit order function is a Heaviside type function ( = 1 on t and = 1 on t ) and the limit temperature remains continuous, but with weak discontinuity on the free boundary t = t t . To formulate the singular limit problem, we suppose that 0 is a smooth surface of codimension 1 , 0 = , dividing into two parts Let
0

0

so that
0

= 0 0

boundary t ; is the outward (relative to t ) normal to t , V is the normal velocity of the front t , t = div ( ) |t is the mean curvature of the surface t , and 1 = 3 / 2 . We assume that t = t0, i.e., the front does not intersect the fixed boundary . Remark 3.4 The boundary conditions (13) and (14) can be interpreted as the Hugoniot type conditions corresponding to the problem of propagation of strong discontinuities of the limit order function and the problem of propagation of weak discontinuities of the limit temperature . This interpretation can be justified as follows. As is known, the necessary conditions for the existence of a shock wave type solution to a quasilinear hyperbolic equation generate an instable chain of Hugoniot type conditions. The same instability conditions (cf. [6]) are obtained for the boundary conditions on the free boundary if we use the classical definition (in ) of a weak solution to the phase field system. The boundary conditions in the interpretation of an IAS-domain as the limit of wave train type solutions are referred to as Hugoniot type conditions. Let us describe the geometric structure. Assume that, at t = 0 , the domain contains domains of pure (liquid or solid) phase 0, and also the melt domain * 0, occupied by a large number of pure phase domains of small volume i0, , i = 1,2, , M , where M is even. For the sake of simplicity, we consider the case of i quasispherical symmetry. Let 0, , i = 1, , M 1 , be i interfaces of domains 0, so that

0 0, i 0,

and

0

be the initial data such that

=

i 1 0,

i0, , =
M 0,

0 = 1 ( )
outside an -neighborhood of the surface 0 and 0 C () (cf. details in [6,9]). The singular limit problem is written as
= , t x , t > 0,
0 t

= 0, ,

0,

.


i We denote by D0, and assume that



i the domains bounded by 0,

(11) (12) (13) (14)

i i D0, D0,1 , i = 0, , M ,

where
0 M D0, = 0, , D0,1 . i Assume that 0, 1 such that

|


t =0

= ( x), x , | = b ,
0

| = 0, t

|t = 2V ,



are smooth surfaces of codimension

1 |t = t V .

c1 dist c1
0,



k 1 0,

,

k 0,

This problem is the well-known modified Stefan problem with the Gibbs-Thomson condition (14) on the free boundary. Here,

c2 , dist



c2 , , c3 ,

M 0,



(15)


where [ f ] |
t

0

x = 0 x

, x 0 ,

denotes the jump in f across the free

where k = 1, , M , (0,1) and the constants cl , c j > 0 are independent of . i Assume that 0, , i = 0, , M , satisfy the following geometric condition. i 0, C 3 uniformly in [0, 0 ] , M and
AM

Copyright © 2010 SciRes.


E. A. LUKASHOV ET AL.

163

M L = const as 0 ; moreover, the surfaces obtained after the limit passage occupy the mixture domain * bounded by C 3 -surfaces 0 and 0 . 0 Remark 3.5 If Condition (A) holds, then there exists a i function s 0 ( x, ) C 3 () such that any 0, is a level surface of this function. Formula (15) shows that there are no interactions (up to ( ) ) between neighboring waves such that the distance between them is not less than ( 1 ) with any constant > 0 . Thus, for sufficiently small t an asymptotic solution is expressed as the superposition of local solutions to the rigid-front problems (with one front i 0, ) (cf. [6,8]) (as shown in Formula (16)). As in the case of rigid-front solution, i ,c is a smooth extension of the auxiliary function i = i ( x, t , h) ,

t0, = t, denotes the domain bounded by t0 and , tM1 = t, denotes the domain bounded by tM and , , . The small corrections c1( j ) ( x, t , h) are simultaneously corrections of order ( ) for temperature which can be computed as solutions to the linearized chain of modified Stefan problems with the GibbsThomson condition (cf. [6]). For the sake of convenience, we impose the following condition (cf. [6]). (A') There exist functions s (1) ( x, t , ) and (2) s ( x, t , ) that describe respectively the surfaces ti, with even and odd superscripts for t0 . We denote by it , the domain bounded by the surfaces ti,1 and ti, , i = 1, , M , and introduce the notation
* t ,

i = s



ni



x , t , ih



s0

ni

, h = , 0, 0 .




=


i =1

M

it , .

( We recall that the family of functions { i } and s0 j ) , j = 1,2 , is defined as a solution to the chain of modified Stefan problems with the Gibbs-Thomson condition

Constructing formal asymptotic solutions, we find
( s ( j ) ( x, t , ) = s0 j ) ( x, t , h) c1( j ) ( x, t , h),

i = i , x it , , t > 0, t

h = , 0, 0 , j = 1, 2,

(17) (18) (19)
,



i 1 i 1 0 t , i t , 0

|

= i |

i t 1 0 ,

, ,

so that | x s ( j ) |> 0 uniformly with respect to x for any h [0, h0 = 0 ] and

i t ,

* t ,

n = x, s0 i x, t , h = ih , ni = 1, i = 2k ,

i |

i 1

=

i 1 i 0 t ,

|





ni = 2, i = 2k 1, 0 i M .

i 1 0 i 1 t ,

|



i | i 1

i 1 0 t ,

= (1)i 1 2V

i 1

(20) (21) (22) (23)

It is obvious that
s
(1)

|t =0 = s

(2)

|t =0 = s 0 ( x, )
( ni ) 0

i | i

i 0 t ,

i 1 | i
ti,1 0

i 0 t ,

= (1)i 2V i ,
i 1 t

and, with accuracy ( ) ,

i = s

( ni ) 0

/ | s

||

(1) 1 i |
i 1

=

V

i 1

,

ti,

(1) i 1 i |

ti, 0

= it V i ,

with the initial and boundary (on ) conditions. Here, i = 0, , M 1 . We set

1 t ,

=

M 1 t ,

= ,

so that the condition (18) (the condition (21)) vanishes for i = 0 ( i = M 1 ). Furthermore,
V i = (| s
( ni ) 0

|)

1

s0 i | t

(n )

ti,

, it = div i |

ti,

are outward normals to Dti, . For fixed > 0 and sufficiently small t > 0 the classical solvability of the chain of modified Stefan problems with the Gibbs­Thomson condition is established in the same way as in [10]. At the same time, based only on the limit problem below, it is impossible to formulate the initial conditions for the temperature 0 ( x, ) in such a way that these conditions have sense as 0 because the classical sol v ability of the modified Stef an p r o b lem with th e Gib b s- Th omson condition assumes conjugate conditions on the initial i surface 0, for any M . However, we can over-

1as ( x, t , ) = (1) i ( i ) { 0as ( x, t , h) i ( i , x)},
1 ( x, t , ) = ( 2
as 0 i =0 i 1, c

M

1 i ,c ) ( 2

2

M

i =0

i 1, c

i ,c ) ( i ), x it ,1 it , .

(16)

Copyright © 2010 SciRes.

AM


164

E. A. LUKASHOV ET AL.
s0
j

come these difficulties if find a model problem for a weak limit of temperature as 0 . Thus, we choose the initial data





x, t , h s0 x, t hs0
j

j





x, t , h , j 1, 2

(29)

|
|

t =0

= 1as ( x,0, ) ( 2 ),
= 0as ( x,0, ) ( ),

(24) (25) (26)

where the functions s0 , s0 and their third order derivatives are uniformly bounded for h [0, h0 ] . As a consequence, we find

t =0

V i = it (h).
By (30) and (22), (23), we have
(1) i 1 1 i |
ti,1 0

(30)

s

( j)

|t =0 = s ( x, ),
0

where ( x,0, ) , ions s 0 ( x, ) are satisfied for fixed these conditions by
as 0

( x,0, ) , and the smooth functsuch that conjugate conditions are > 0 . We will be able to specify obtaining the limit problem.
as 1

=

i 1 t

V

i 1

,

(1) i 1 i |

ti, 0

= it V i ,

3.2. Limit Problem
The evolution of solutions can proceed in two different ways depending on the initial data:
t s
(1) (1)

which implies

i |

ti,

= (h).

From Condition (B) it follows that
x, t , h 1 , 1 1 , x , , t 0, T , (31) t where 1 i1 for x it , . By the Gibbs-Thomson law,
ti ti 1 Vvi Vvi1 i 1 1 i1 | h h
i t , 1
i t , 0

|t

=0 =0

t s

(2) (2)

|t =0 < 0,

(27) (28)

t s
( j)

|t

t s

|t =0 > 0,

where s , j =1, 2, are the functions in Condition (A'). In the case (27), the boundaries move in the opposite directions. Consequently, the wave train type structure exists only during a small time interval since the domain t2,k or t2,k 1 vanishes for t . A similar situation for the classical Stefan problem was treated in [6]. In teh case of the phase field system, from (27) it follows that an "overheated" or "overcooled" domain appears in * . t To find conditions for the existence of wave train type solutions in some finite time interval independent of , we consider the case (28), where the boundaries move in the same direction. Assume that the following condition holds. (B) There exists T >0 such that for any 0 t T there exist functions i ( x, t , h) , i = 0, , M 1 , such that the function x, t , (defined by i for i x t , ) is continuous and is uniformly bounded for [0, 0 ] . Fur t her more, i C1 (Qi ) un ifor ml y for

i1 |

i 1 t , 0

.

Since C uniformly with respect to h , we find i1 C Qi . We need the following assertion. Lemma 3.1 1) Let i be partition points in the interval [0, L] , 0 < 1 < < M , and let h = max i i i 1 . Suppose that M is even, F C 0, L , and F C1 i 1 , i for any i =1, , M . Then
3



i 1 F i i =0

M

const uniformly for M 2. F ( )

2) Assume that F ( ) C ([0, L]) and C 2 ([ i 1 , i ]) for any i = 1, , M . Then
1 F 0 F 2 i 0 uniformly for even M 2.
i 1 F i M



M



h



[0, 0 ] , where Qi =

t[0,T ]





i t ,

, and

i t ,

C3 .

We list some consequences of Condition (B). Since the functions i are smooth, it is obvious that

i |

ti, 0

i |

ti,1 0

= (h).

Therefore, taking into account the Gibbs--Thomson law (22), (23), we find
it V i
i 1 t

To prove the lemma, it suffices to group the terms in F ( i ) F ( i 1 ) in such a way that to represent them as differences of derivatives. Note that for passing to the limit in the wave train as 0 , we need a suitable well-defined notion of a weak solution. We give such a definition in accordance with [6]. Definition 3.1 A pair of functions

V

i 1

= (h).

L2 0, T ;W21 L 0, T ; L2 ,
W22,1 Q L 0, T ;W21 L4

Since the surfaces we have

i t ,

are smooth and V i 1V i > 0 ,




AM

Copyright © 2010 SciRes.


E. A. LUKASHOV ET AL.

165

is called a weak solution to the problem (49) if for any test functions ( x, t ) , g ( x, t ) = ( g1 ( x, t ), , g n ( x, t )) satisfying (38) the functions and satisfy the equation
I =

sum and using (31) together with Condition (B), we find
J = ( g , s0 ( it V i ) A 2 , ( i )) (h 1 h) = 0.
i=0 M

Q


, t dxdt
0

(35) We again obtain the relation (30) since the first sum in (35) has order (h 1 ) . Taking into account (29) and passing to the limit as 0 , we see that (35) implies (30) in the entire domain
* = lim*, . t t
0



0 x,0 )dx = 0



(32)

and the integral identity (33) where
e = 1 1 2 W , W = 2 1 2





2

4,

Consequently,
s0
1

and g x is the matrix with entries ( g x )ik = gi / xk . We set
i = ti, , i = 0, , M ,
t[0,T ]

s s0 = div 0 , s t 0

x * , t > 0. t

(36)

and



M 1



M 1

= [0, T ].

We consider the integral identity (32). We first compute the weak limit of wave train in the derivative t in the heat equation. Lemma 3.3 Let

Then we substitute (34) in the integral identity (33). We need the following assertion.
Lemma 3.2 Suppose that ( , x) S, ( x) C 2 () , | |= 0 , and
dist t , const > 0.

( x, t , ) = 1as ( x, t , ) ( 2 ),
where 1as is defined by formula (34), and
c1hdist ( i ,

i 1

)c2 h, i = 0, , M 1,

Then for any function g C Q
1


x, dx,

where the constants c1 and c2 are independent of . Then
(
M , ) = 2({ (1) j =0 t j 1

lim
0

1

, x g x, t dxdt Q
1

s



V j ( j )}, ) C1 (h 1 h)

=

T



A x

x g

for any functions C Q
1



(37) such that (38)

where s = (t ) s1 , =| |1 ,
A = , x d ,


, g C1 Q , | = g | = 0, |t =T = g |t =T = 0
Here,



1 C1 s0 s0
0

and T is the domain bounded by By Lemma 3.2,
M i=0

and T .



2




h 0

is a possible contribution of the terms depending on the first corrections to the phase s0 relative to h . We set F
k

J = ( g , s0 ( it V i ) A 2 , ( i ))

(1) i ( g , s0 A , ( i )) (h 1 h) = 0.
i =0

M



k

Vvk d k .

Applying assertion (a) of Lemma 3.1 to the second
J t g , dxdt e divgdxdt
Q Q

Applying assertion (b) of Lemma 3.1 to (37), we find

Q



, g x div g dxdt 0 ,
as 0

(33)

1as x, t ,

i 1 i i 0

M

2



x, t , h

i i , x
i 0

M

.

(34)

Copyright © 2010 SciRes.

AM


166

E. A. LUKASHOV ET AL.

(

, ) = V 0 d 0 V M d 0 M t

M

C1 (h 1 h).

= 0, x * , t 0, t
s s0 = s0 div 0 , s t 0 x * , t > 0, t

(44) (45)

(39) We recall that, by Condition (B), x, t , is bounded in L (0, T ;W21 (* t with respect to and, consequently, * verges in L (0, T ;W21 (* )) ; moreover, by t 0 as 0 for x * in the t L2 ((0, T ) * ) -convergence. Thus, t





the family )) , uniformly -weakly con(31), we have sense of the

|* = 0,
t

| * = Vn , t 0, n t

(46) (47)

|t =0 = 0 ( x), x \ * , 0
s0 |t = 0 = s 0 ( x), x * , | = b , 0 where
* = t t , t t = x , s0 x, t = 0 , t = x , s0 x, t = L ,

x, t lim x, t , 0, x . t
0

def

It is obvious that (31) does not contradict sign of the leading term of corrections j s0 ) of velocities is independent of C1 =0 . On the other hand, in the domain ing heat equation has the free term C1 . fact, one should prove that
s0 s0 h .
1 2

(39) if only the (depending on j and then * , the limitt To verify this

n denotes the outward normal to | s0 |1 s0 / t | * , and s 0 ( x) s 0 ( x, 0) .
t

* , t

Vn =

The proof is given below in the spherically symmetric case. Now, we continue computations in the integral identity (32). Integrating by parts
I
def

Q



, t



dxdt x, 0 0 dx,


(40)

we find (41) where



(i )

= |

i Q

.

By (31) the integrals over it , and i , i =1, , M , converge to zero as 0 . We recall that, by Definition 3.1,
I I t dxdt 0 x,0 dx 0.
Q


(42)

Taking into account (30), (39), and (41), we arrive at the required result as 0 :
= , t x \ * , t > 0, t

(43)

Thus, the problem (43)-(46) can be interpreted as two classical one-phase Stefan problems joined by Equation (45). Such an interpretation leads to the problem about the mixture domain for processes with surface tension (cf. [6,8]). The conditions (45), (46) and = 0 on * are t conditions of Hugoniot type since they should be satisfied for the existence of the solution under consideration. The operator on the right-hand side of (45) degenerates along the direction s0 , i.e., along y1 if we introduce the new coordinates y1 = s0 , y2 , , yn , where y2 , , yn are the coordinates on the surface s0 = const . Equation (45) is ultraparbolic. As is known [11], a homogeneous ultraparabolic equation has no real analytic solutions with respect to t and y1 , except for the case where the solution is independent of the tangent variables. Further, we need to solve the Cauchy problem (46) for the heat Equation (43) relative to y1 with the initial conditions on the surface * . For sufficiently small y1 and t t this ill-posed problem has a solution only for real analytic surfaces and initial data [11]; moreover, in this case, the values of on the external boundary and at the initial time are uniquely determined by the values on
dx dt

T I 0

0 t ,







0

t
0



0

dx

tM 1 ,








M 1

t
MT





M 1



0 t ,







v0

d 0

M









M 1

vM

d

M

i 1

0 i , t

i dxdt t

(41)


i

i

d i vi

i 1





i

d vi 1

i 1



* 0,



x,0 0 dx,

Copyright © 2010 SciRes.

AM


E. A. LUKASHOV ET AL.

167

3.3. Example of a Structured Domain
Assume that n = 3, ={x, R < r < R } , where r =| x | , R > 0 , and * ={r (0) < r < r (0)} . Then Equation 0 (45) becomes the first order equation s0 2 s0 = , r * = r (t ) < r < r (t ) , t > 0. t r r t (48)

It is easy to solve the problem (48) with the initial condition
s0 |t = 0 = s 0 (r ).

0 where 0 < R < r00 < r10 < < rM < R . Assume that 0 0 0 rj 1 rj = h and the differences r00 R and R rM 0 are sufficiently small; moreover, s (r ) is real analytic, s 0 / r > 0 , 0 ( x ) and b are special data corresponding to the solution to the Cauchy problem for the heat Equations (43), (46). We show that Condition (C) implies Condition (B) and the equality C1 =0 in (39). For this purpose, we return to the main problem (cf. (17)-(23))

i = i , x it , , t > 0, t

(50) (51) (52)
i 1

Namely,

s0 r , t = s 0 r
along the characteristics
r r ,t =


0
0



i 1 i 1 0 t , i t , 0

|

= i |

i t 1 0 ,

, , ,

i |

=

i 1 i 0 t ,

|



0




r
0

2

4t , r 0 r r 0





i 1

i 1 0 i 1 t ,

|



i | i 1

i 1 0 t ,

= (1)i 1 2V

(53) (54) (55) (56)

for any smooth function s 0 (r ) such that sr0 > 0 . Now, (43), (46) with
Vn = 2 / r |
* t

i | i

i 0 t ,

i 1 | i
ti,1 0

i 0 t ,

= (1)i 2V i ,
i 1 t

(1) i 1 1 i |

=

V

i 1

,

is the Cauchy problem (with respect to r) in two domains

(1) i 1 i |

Q1 = R < r < r t , t > 0 , Q2 = r t < r < R , t > 0.
To formulate the solvability conditions for this illposed problem, we recall the well-known fact (cf., for example, [11]): for the local existence of a solution to (43), (46) it is sufficient that the curves r (t ) be real analytic functions with respect to t , i.e., r (0) > 0 and t < r2 (0) / 4. Consequently, for sufficiently small 0 > 0 and T0 = T0 ( 0 ) , in the domains
Q1* = r 0 0 < r < r t ,t < T0 ,
* Q2 = r t < r < r 0 0 ,t < T0

ti, 0

= it V i .

Let i = i (t , h) be functions such that it , = {r , r = i (t , h)} . In the spherically symmetric case, we have
it = 2/ i .

Therefore, taking into account (30) and choosing i directed in the opposite direction relative to the normals (with respect to Dti, ), we find
V i = 2/ i (h).



We make the change of variables i = wi / r . Then the equality
t = , x it , , t > 0

there exists a real analytic solution to the corresponding Cauchy problem. Thus, in order to solve the limit problem (43)-(46), we need to impose the following condition. (C) Suppose that is a spherically symmetric layer in 3 , the initial and boundary data of the problem
t t = ,

takes the form
wi 2 wi = , r t r 2
i 1

t , i t
def

, t > 0.

(57)

Since
V i = 2 i1 (1 hv i ), v i = i ( it V i )/2h,

2 t = 2 3 |t = 0 = 0 ( x, ), |t = 0 = 0 ( x, ), | = 1, | = b
are spherically symmetric, and
i 0,

(49)

the conditions (51), (54) can be written as follows:
wi 1 | r wi | r
i 1 0 t ,



wi | r

i 1 0 t ,

= 1 4 1 hv
i



i 1



,

(58) (59)

= x , | x |= ri0 ,





i 0 t ,



wi 1 | r

i 0 t ,

= 1

i 1

4 1 hv i ,



Copyright © 2010 SciRes.

AM


168

E. A. LUKASHOV ET AL.

wi 1 | wi |

i t ,1 0

= wi |

i t 1 0 ,

, .

(60) (61)

tion. Denote by R( z , t , h) a solution to the equation
s0 R, t hs1 R, t = z.

i t , 0

= wi 1 |

i t , 0

We show that the problem (57), (58) has a solution satisfying the following properties: 1) wi = (h) uniformly with respect to i , i ^ 2) for any t the values wi 1 wi r pi are deter2 mined, with accuracy (h ) , by the values of some ^ function wi C1 0 , M on the grid {0 , , M } . We note that the first property is related to (56) and (30). We look for a solution wi to the problem (57), (58) in the form wi = ai r i 1 bi t ui t , r , h , (62) where the first tion (58) and problems: ui t two terms correspond to the Stefan condiui is a solution to the following chain of
2 ui = ai r 2
bi , i = 1, , M ,

By construction, i = R (ih, t , h) and, uniformly with respect to i up to order (h) , the functions v i are traces of some C1 -function v on the surfaces r = i . We note that R / z > 0 . Furthermore,

b j = 2h (1)
k =2

j

k

R | z

z = h ( k 2)

(h 2 ) = (h).

By Lemma 3.1, the last estimate is uniform with respect to j . Further,
b
j2

b j = 2(1)

j 1

(

j 1

2 j

j 1

) (h 3 ) = (h 2 ),

(66) and this estimate is also uniform with respect to j . Now, we see (67) Furthermore, from (66) and Lemma 3.1 it follows that
b
j 2l

i 1

(63)

b j = (h 2 )



uj u

j 1



|r

= j

u j u j 1 = 0, |r r r

= j

= 0,

(64)

j = 0, , M .

uniformly with respect to j and l . In particular, from this estimate, the equality (67), and the condition b1 = 0 we find R b2l = 2h | z = (2l 1) h (h 2 ), b2l 1 = (h 2 ). z We consider a broken line such that its linear parts are defined as ai (r i 1 ) bi on the segments [ i 1 , i ] . It is obvious that bi are the values of at the points r = i 1 . Consequently, is not symmetric with respect to the zero line (it is directed toward to the domain of positive values). However, the broken line can R |z = h (i 1) in be centered by decreasing its values h z each segment [ i 1 , i ] . It is obvious that this is equivalent to the existence of functions R | z = z ( r ,t , h ) m=h z in . Here, z = z (r , t , h) satisfies the equation R ( z , t , h) = r . We set
1 = m, U i = ui m.

We note that this chain is similar to that considered in [12] and differs by only the dependence of fi = ai i 1 bi in (63) on t . However, because of this dependence, it is obvious that the contribution of this chain to the solution is of order (h 2 ) . To solve Equation (63), we first compute the coefficients ai and bi . From (58) and (63) it follows that
ai = 2 1 bj = 2
i 1

1

hv i , b1 = 0,



k =2

1

j

k

1 hv



k 1





k 1

1 hv



k 2





k 2

,

j = 2, , M . Assume that s0
j





x, t , h s1 x, t h , j 1, 2,

(65)

j where the functions s0 are defined in (29). At the first glance, this assumption can lead to a contradiction in the equation for velocity correction (the linearized Gibbs j Thomson equation for s0 ) if the functions i computed under this assumption do not satisfy Conditions 1) and 2). However, it turns out that there is no contradict-

Then for U i we have the problem of the form (63) with the right-hand sides
Gi = ai
i 1 2 m m , r , . bi i 1 i t r 2

(68)

b

j 1

b j = 2(1)

j 1

(h

R h 2 2 R h 2 ( Rv)) | 2 2 z z z

z = ( j 1) h

(h 3 ).

(67)

Copyright © 2010 SciRes.

AM


E. A. LUKASHOV ET AL.

169

To construct an asymptotic expansion of U i , we solve a chain of problems. We look for a solution in the form
U i = ci r

conditions
s n j 2 s ni 0 0 t r r




i

r

i r ci1 c1i r i 1 i
i 1

r r

i 1



2



r i

1
i

i

1 s0 h r

r i

h .



2

,

(70) Our analysis shows that, with accuracy (h) , the right-hand side of (70) is the trace of a function of class C 1 . Therefore, from (69) and the conditions 1 s0 we find
s1 2 s1 t r r
r i t 0

where dots denote polynomials of higher degree. We note that polynomials of degree higher than 2 admit the estimate O(h3 ) and the coefficients ci are determined by the relations
ci = 2(1)
i 1

s0

2



t 0

0
s0 r



i 1

(h), i = 1, , M .

The contribution to the solution U i of terms of order O(h) in Gi is estimated by (h 3 ) . Consequently,
Ui U i h



1 R h z

z i 1 h

r i

h . (71)


3

Let
i t , h ri t hri t , h ,

and the function

U i ci r

i 1



i

r



so that
ri t , h ri t

is defined by a sequence that is symmetric with respect to the zeros of parabolas of order mod (h 3 ) because
ai
i 1





1

ai 1 i = (h).

Hence U i = (h 2 ) for r ( i 1 , i ) and the values of 1 at the points j are given by the relation

uniformly with respect to i = 0, , M . Taking into account Equation (48), we obtain
ri = g
2

ih

4t ,

1 |

r= j

= (1) j h

R | z

z = ( j 1) h

(h 2 ),

j = 1, , M .

(69) Thus, the problem (57), (58) has a solution with properties 1) and 2). It remains to construct in the domains R r 0 (t ) and M (t ) r R . We note that constructing 1 , we defined mod (h) the values of and r at the points r = 0 (t ) and r = M (t ) . As in the case (43)-(46), this fact completes the construction of . Now, it is again required to solve the Cauchy problem with respect to r for the heat equation. Nevertheless, by Condition (C), the analyticity condition (necessary for solvability) is already valid. Thus, by (69), the functions
i t 1
i 2 r i

where g ( z ) is the inverse of s 0 , i.e., s 0 ( g ( z )) = z . Thus, ignoring terms of order (h) , we can transform (71) as follows: s1 2 s1 1 = , s1 | r r r t
t =0

= 0,

i.e., our assumption about s1 (r , t ) is valid. We note that, in view of (65), the value C1 in (37), (39) is equal to zero and consequently, right-hand side of the heat equation in * vanishes. t Thus, Condition (C) implies the validity of Condition (B). As a result, we find (43)-(46) as the limit of the chain of Stefan problems with the Gibbs-Thomson condition. We formulate the initial conditions. We assume that Conditions (A) and (C) are satisfied. Let

,
i t ,

|t =0 = 1as x, 0, O
0 0

2



, S j |t = 0 = s

0

x,

,

with accuracy (h ) , are traces on the surfaces of some function

where s ( x, ) = s (r ) . Let |t = 0 in the domains i0, = {ri0 1 < r < ri0 } , i = 1, , M , is defined by

x, t , h h
1



(i ) |

t =0

= (1)

i 1

2 h ( (r ri0 1 ) (h 2 )), r 2( s 0 )r

of class C . Owing to this fact, we can compute the first correction for the phase s0 (r , t ) . Indeed, substituting (29) into (56), we obtain the linearized Gibbs-Thomson
Copyright © 2010 SciRes.

0 and, in the domains R < r < r00 and rM < r < R , we set

|t =0 = |t =0 ,
AM


170

E. A. LUKASHOV ET AL.

where is a solution to a special Cauchy problem (relative to r ) for the heat Equation (43). Theorem 3.1 Under the above assumptions, there exists an asymptotic solution to the phase field system satisfying Condition (B), and it is possible to pass to the limit in (49) as 0 in the sense of Definition 3.1. The limit problem
= , x \ * , t > 0, t t

(72) (73) (74)

= 0, x * , t 0, t
s s0 = s0 div 0 , s t 0 x * , t > 0, t

|* = 0,
t

| * = Vn , t 0, n t

(75)

|t =0 =
where

0 0

s0 |t = 0 = s

x, x,

x \ * , 0 x * , | = b , 0

(76)

* = t t , t

t = x , s0 x, t = 0 , t = x , s0 x, t = L ,

n denotes the outward normal to * , Vn = t | s0 |1 s0 / t | * and s 0 ( x) s 0 ( x, 0) , possesses a t solution, at least for sufficiently small (but independent of ) time. The above case can be explained by the fact that, outside the layer r0 r rM , the order function of the original problem takes different values: 1 for r < r0 and 1 for r > rM . It is obvious that all the arguments remain valid in the case where takes the same values ( 1 or 1 ) for r [r0 , rM ] . This means that M is odd. Then we again obtain a limit problem of the form (72)-(75). The limit passage can be justified in the same way as above, by solving a chain of problems which can be reduced to the chain of problems (63). In both cases ( M is even or odd), the problems are ill-posed. However, as was noted in [6], such a wave train type structure appears in numerical experiments as solutions to the phase field system with the initial data 0 = 0 for R r R and
M j 1 1 th j =0 0 = M j 1 th j =0
Copyright © 2010 SciRes.

Figure 2(b) presents the graphs of solutions to the phase field system with spherically symmetric initial data for M =19 and =102 at different times. One can see that the temperature in the mixture domain is of sawtooth form. Such a function is the leading part of the asymptotic expansion (62) of the solution to the chain of modified Stefan problems with the Gibbs-Thomson condition. In the numerical analysis performed by O. A. Vasil'eva, = 0 on the external boundaries. This leads to an effect presented in the figure for time t = 0.02 : the sawtooth structure begins to break down under the influence of boundary data. However, the order function is more stable and preserves its shape. Figure 2(a) presents the graphs of solutions to the phase field system with spherically symmetric initial data for M = 7 and =102 at different times. The temperature has sawtooth shape in the IAS-domain, whereas it is periodic with amplitude l = 1 at center. Such a function is the leading part of the asymptotic expansion (62) of the solution to the chain of modified Stefan problems with the Gibbs-Thomson condition. The sawtooth structure "moves" to the center and begins to break down under the influence of nonspecial boundary data. The order function preserves its shape in this case. Figure 3(a) presents the graphs of solutions to the phase field system with spherically symmetric initial data for M = 7 and =102 at different times. Figure 3(b) presents the graphs of solutions to the phase field system with spherically symmetric initial data for M =19 and =102 at different times.

3.4. Comments and Conclusions
Based on the phase field system, it is possible to detect a banding structure formation in instability zones. However, to construct the mathematical model, we need to impose some restricted conditions. 1) The existence conditions are very restrictive, which can be explained by the geometry of domain and the initial and boundary conditions. Note that the initial and boundary data are determined by the solution to the limit problem. 2) A standard definition of a weak solution can turn out to be not suitable. However, we can avoid these difficulties by introducing a special definition of a weak solution, which is important for nonlinear problems. 3) As was shown, a wave train type solution exists only for special boundary and initial data providing the existence of an asymptotic solution to the chain of Stefan problems with the Gibbs-Thomson condition for sufficiently small (but independent of ) times. This fact allows us to pass to the limit of the chain of Stefan problems with the Gibbs-Thomson condition (in the sense of
AM

r rj0 , M is odd, r rj0 M is even. ,


E. A. LUKASHOV ET AL.

171


r r


r r


r r


r r





r

r

r

r

Figure 2. (a) Presents the graphs of solutions to the phase field system with spherically symmetric initial data for M =7 and =102 at different times; (b) Presents the graphs of solutions to the phase field system with spherically symmetric initial data for M = 19 and =102 at different times.

Definition 3.1) and derive the limit problem (43)-(46). 4) As we shown in the above examples, the temperature ( x, t , ) is small ( 0 as 0 ) and has special "periodic" structure in the stratified domain. 5) Even in the rigid-front case, the solid phase growth is of order ln (1 / t ) , which is lower than the order obtained in experimental way. Thus, a banding structure in the phase stratification domain of a binary alloy was constructed under extremely restrictive conditions on the geometry of domain and the initial and boundary conditions. Furthermore, the order O( t ) of the solid phase growth obtained in experiments is not achieved in this model. In view of these facts, it is necessary to look for other mathematical models describing qualitative experimental properties of crystallization. In the following section, for such a model we consider the convective Cahn-Hilliard equations in a
Copyright © 2010 SciRes.

porous medium of an overcooled melt.

4. Oriented Crystallization Model
There is a huge experimental literature on various structure formations in melt crystallization. Based on experimental results, one can conjecture that complex structure formations in crystallization are caused by the evolution of instabilities during phase transition processes which, in turn, is caused by different reasons and can be realized in different ways. We list some of such reasons. 1) concentration overcooling, 2) convective flows deforming the temperature field (gravity and thermocapillary convection), 3) phase stratification. In addition, elastic properties of the solid phase, thin phase boundary, and adsorption phenomena can also
AM


172

E. A. LUKASHOV ET AL.

contribute to this effect.

4.1. Modified Convective Cahn-Hilliard Model in a Porous Melt
To construct a mathematical model governing the reconstruction of oriented crystallization (cf. [7,8]), a modified Biot model of a porous medium [13] was used for describing a liquid-solid zone and the convective CahnHilliard model of spinodal decomposition [14,15] was used for describing segregation. In the model, we consider a binary eutectic alloy. For variables we take the concentration of the component A or the component B of the binary alloy, the temperature, the growth velocity of the solid phase, the contraction, the convection velocity of the liquid phase. The model includes the laws of conservation of mass and impulse for liquid and the law of conservation of total impulse for liquid and solid phases. In accordance with the physical interpretation, the model also includes a modified Cahn-Hilliard equation [14] and the heat equation [7], regarded as a generalization of the Stefan problem [9]. Using a nonisothermic modification of the Cahn-Hilliard model, proposed in [7], we can construct a model that take into account the following physical effects. Because of crystallization and melting, the temperature can vary. In turn, variations of temperature lead to variations of velocity and changes of the medium composition. An equilibrium phase transition is realized at the melting temperature, whereas a nonequilibrium phase transition can be realized at different temperatures depending on the depth of penetration into metastable or labile regions. This fact shows that the modified Cahn-Hilliard model should include temperature-dependent parameters. Then both heat-mass transfer equations will govern mutually dependent processes. The model reflects the structure of a liquid-to-solid transition zone of the crystallization front. It consists of an outer viscous layer (the hydrodynamic Prandtl layer) and a diffuse layer (the Nernst diffusion layer). In the case of a condensed system, the thickness of the Nernst layer is less than the thickness of the Prandtl layer by three orders and the heat-mass transfer laws can be assumed to be linear (the Fick and Fourier laws). On the boundary of the diffuse layer, near the solid phase, the volume strongly varies while a liquid-to-solid transition. Therefore, it is necessary to take into account elastic forces, which can be done within the framework of continuum mechanics.
Copyright © 2010 SciRes.

We introduce the following notation: c is the mole concentration of the component B in the binary alloy (In our case, the mole concentration of Sn in the liquid phase), z is the contraction, wl is the convection velocity of the liquid phase, u is the averaged displacement in the solid phase, ws is the mean growth velocity of the solid phase, v is the averaged fictitious displacement in the liquid relative to the solid phase, T is the temperature. Furthermore, we set
w = wl ws .

4.2. One-Dimensional Case
In this case, the model is represented (cf. [8]) by a system of differential equations which can be divided into the three subsystems:
s ut = w , vt = w, ws l w = 2 2 M u Mv g , t x x t x l s l w t add wt = Dw Mu x Mvx x g , (77)










l t
x

l w s w
x

s



x

= 0,

(78)

ct w(c ckr )c

2 = [ M D ( F (c , T ) u x c ( 2 )

2 ( F1 (c, T )c x ) x 4 ( F2 (c, T )c xx ) xx ) x ] x , (T c) t = D0 Txx ,

(79) This system describes processes in the diffuse and Prandtl layers in dimensionless variables c , T , wl , ws , u , v , z . The system (77) is a model of wave propagation in a porous skeleton filled with a liquid (a simplified version of the Biot model). The system (78) is the continuity equation and describes the evolution of contraction. The system (79) presented by the convective CahnHilliard model and the heat transfer equation describes the formation and growth of Gibbs grains. Note that we use equations of continuum mechanics to describe processes in the Prandtl layer, whereas for diffusion and heat processes we use the modified CahnHilliard model where hydrodynamic processes and elastic-plastic state of the solid phase are taken into account. Let's note that all constructions of the previous chapter were maded for this one-dimensional case but more technically.
AM


E. A. LUKASHOV ET AL.

173

The model contains a number of dimensionless parameters. The elasticity modulus of the solid phase = 2 is assumed to be a function of concentration and temperature: = (c, T ) . The parameter z is expressed by the formula
z= V V s V i , V

composition and temperature. The function F can be approximated by a cube polynomial in c at a fixed temperature:
c c c ccr c F c, T = 3 c ccr ,









, T T0 T T0 ,

where V is the total melt volume, V s is the current volume of the solid phase, and V l is the current volume of the liquid phase. The concentration c is expressed as

where T0 = 400 K and c , ccr are functions of temperature T which will be specified below,
Tmin T Tmax , Tmin = 233,15K , Tmax = 456,15 K .

c=
where M A A ( B ) and melt: M l = The extra expressed in

s mB M B , ml M l

( M B ) is the atomic mass of the component M l is the averaged atomic mass of the (1 c) M A cM B . l variable y (c ) of the form y = m B /m B is terms of concentration as follows:
y c = qc , 1 c Q q y c

We cmin cmid c0 We

define three concentration values: equals to c(Tmin ) , equals to ccr (Tmin ) , equals to cmid in our experiment. set
cmin = 0, 04, cmid = 0, 43.

We introduce c (T ) as the roots of the equation
T =
clust

c2

clust

c

clust

,

(80)

where

which implies



clust

=



Tmin T0 cmin c0



2

,

clust

= 2

clust 0

c,

l c, z =
2/3

1 z R 1 y c





R,



clust

= T0

2 clust 0

c

is a bound for the surface of the where (1 y (c))) solid phase in the liquid-solid region, 1 is a parameter,
R= mSn m M = 0, 74, q = Pl = 0, 25, Q = sV mSn M
Pl Sn

= 1, 75

are constants, and we adopt the normalization condition

s = B = 1.
Further, is the mean density. M = 6,2 is the mobility (fluidity) of the liquid. 2 is the inverse of the relaxation time of fluidity (estimated as 10 3 ), g = 1/30 is the acceleration of gravity, D is the interphase friction coefficient, estimated as
D = 1e
/ T

and define ccr (T ) as a linear function. The function F1 (c, z , T ) 1000 is interpreted as viscosity. At the first step, it is assumed to be constant. The structure of the interphase boundary at atomic level is characterized by the function F2 . We set F2 0 and = 10 4 in the numerical experiment. The model also contains some additional relations dictated by the physical interpretation of the problem. In particular, the model contains the "extra" density add such that



add

=

z

l

1

y c



1/ 3

,

(81)

where = 0,055 and, as a rule, ( z ) is small for small z .

1

e

1/ T

1 y c

2/3

,

4.3. Numerical Analysis of the Model Describing Oriented Crystallization. One-Dimensional Case
The numerical results obtained by Rykov and Zaitsev [16] are presented in Figures 3(a-c). Note, that the spatial x -axis is directed upward, whereas the t -axis is directed rightward along the horizontal line. The systems presented in Figures 3(a-c) differ by the value of the parameter . The numerical results show
AM

where = 1 , 1 = 5, M D is the diffuse mobility of the component B (estimated as M D = 1/T ), is the ratio of the melting enthalpy to the heat capacity of the solid phase at a constant pressure. The function F determines steady, metastable, and labile states of the system "melt-alloy" depending on the
Copyright © 2010 SciRes.


174

E. A. LUKASHOV ET AL.

(a)

(b)

(c)

Figure 3. Numerical simulation of the model describing oriented crystallization (one-dimensional case). The spatial x-axis is directed upward, whereas the t-axis is directed rightward along the horizontal line.

that the balance of convective and diffusive generates a modulated wave of formation of grains, which differentiate the spinodal decomp mechanism from the classical case where the Hilliard model possesses a periodic solution.

terms crystal osition Cahn-

4.4. Comments
1) In our model, for the sake of simplicity, we assume that porosity is constant, passing its functions to the contraction z. On this stap of the model construction we will elucidate the change range of the porosity when the modification of Biot model don't lose the hyperbolicity. It allow us on the next stap to pass to porosity as a problem variable, expressed the contraction as the function of porosity. 2) In the model (77)-(79), the convention is equal to zero at the initial time, w |t = 0 = 0 , and then it can be regarded as reaction to 1) the force of interphase friction between liquid and solid phases and 2) the gravity force. Thereby we specify the effective force in the convective Cahn-Hilliard model [14,15]. 3) The initial distribution of crystal grains (which, unlike [14], is not given here) depends on only contraction, whereas the further distribution is determined by the process. So, no restrictive conditions are imposed on the initial-boundary data, unlike the case of the phase field system and the one-dimensional convective Cahn-Hilliard model. 4) In the subsystem (79), we took into account the results of [17]. Note that the above constructions remain also valid for the modified model (77)-(79) obtained from the two-dimensional model (cf. below) in the radial-symmetric case.

c is the mole concentration of Sn in the liquid phase, z is contraction, w l is the liquid phase velocity u l is the averaged displacement in the solid phase, w s is the mean growth velocity of the solid phase (the averaged velocity of microfronts), w = w s wl is the averaged fictitious displacement in the liquid phase relative to the solid phase, T is the temperature. The system of two-dimensional equations can be written in the form

us1

t
t

= ws1 , = w1 ,
t



u

s2 t
t



= ws 2 , = w2 ,

(82) (83)

u1

u2

( ws1 )t l ( w1 )

=[( (c, T ) 2 (c, T ) 2 M (T ))(us1 ) ( (c, T ) 2 M (T ))(us 2 ) M (T )((u1 ) x (u2 ) y )]
x y

x

(84)

[ (c, T )((us1 ) y (us 2 ) x )] y ,

( ws 2 )t l ( w2 )t g
=[ (c, T )((us1 ) y (us 2 ) x ]
2 x x y

( (c, T ) M (T ))(us1 )
2

(85)

[( (c, T ) 2 (c, T ) M (T ))(us 2 ) M (T )((u1 ) x (u2 ) y )] y ,

4.5. Two-Dimensional Case
Introduce the notation:
Copyright © 2010 SciRes.

l ( ws1 )t add l ( w1 )t = D(c, T ) w1 [ M (T )( ((us1 ) x (us 2 ) y ) (u1 ) x (u2 ) y )))]x , l ( ws 2 )t add l ( w2 )t l g = D(c, T ) w2 [ M (T )( ((us1 ) x (us 2 ) y ) (u1 ) x (u2 ) y )))] y ,

(86) (87)

AM


E. A. LUKASHOV ET AL.

175
s2

( l )t ( l w1 ) x ( l w2 ) y s ( ws1 ) x s ( ws 2 ) y = 0 (88) ct w1 f (c, T ) x w2 f (c, T )
y x

u1 = w1 = us1 = ws1 = x u2 = x w2 = x u = x ws 2 = 0 for x 0 and x 1,

={M D (T )[ F (c, T ) 2 ( F1x (c, us , T )cx ) 2 ( F1 y (c, us , T )c y ) y
el

u2 = w2 = us 2 = ws 2 = y u1 = y w1 = y u = y ws1 = 0 for y 0 and y 2.

s1

Eel 6 (6) ]x }x (89) c {M D (T )[ F (c, T ) 2 ( F1x (c, us , T )cx ) x 2 ( F1 y (c, us , T )c y ) y
el

Eel 6 (6) ] y } c

2) The initial and boundary conditions on z have the form z 0, x, y = 0, n z | = 0. 3) The initial conditions on c are as follows:
c 0, x, y = c


y

(T c)t D0 (Txx Tyy )

(90)
2

where
Eel = ( (c, T ) 2 M (T )((us1 ) x (us 2 ) y ) M (T )[(u1 us1 ) x (u2 us 2 ) y ]((us1 ) x (us 2 ) y ) 1 2 (c, T )[((us1 ) x ) 2 ((us1 ) y (us 2 ) x ) 2 ((us 2 ) y ) 2 ], 2
(6)

for ( x, y ) ( xz1 , xz 2 ) (0, yz ) , xz1 =1 / 3 , xz 2 = 2 / 3 , y z =1 / 3 and
c o, x, y = ccr Ti
nt



,

= F2x c, us , T c



F
xx xx y 2

c, u s , T c


yy

yy

.

where Tint = 300 , in the remaining domain. The boundary conditions on C are written as n c | = 0, n * | = 0, where

The system is considered in the rectangle = [0,1] [0, 2] . The boundary conditions are specified by numerical experiments. Here, we write out general boundary conditions. 1) The vector-valued functions u and w satisfy the initial conditions
u s 0, x, y = u 0, x, y = 0,

* = F c, T

2

2

F
x 1 y

c, u s , T c
y


x

x

F
y 1

c, u s , T c



Eel , el c

which corresponds to
ws 0, x, y = w 0, x, y = 0.

Based on the one-dimensional model, we impose the boundary conditions
u s = ws = 0 for x = 0 and y = 0,

n u s = n ws = 0 for x = 1 and y = 2.

2 Eel = us1 x us 2 y c c 2 2 2 us1 x us1 y us 2 x 2 us 2 y c c, T = 20 c T c T , c c, T = c T c T , c i.e., for x = 0 and x = 1 the second condition the form x F c, T 2 vis cxxx 2 cxyy















2

,









takes

We assume that the displacements and velocities satisfy the following conditions on all four boundaries:
n u = n w = 0.

2 el x us1 x us 2 y c 2 2 us1 x us1 y u el c









s2 x





2

2 u



s2 y





2



At the same time, it is natural to impose the impermeability condition on all the boundaries. Therefore, the boundary conditions can be modified as follows:
u x = wx = u
sx

= 0, where
dc x F c, T = x c dT dc c c x c cr x dT xT c ccr T c c



= ws = ws



x

= xu y = x w





cc





y

= x usy = x wsy = 0 for x 0 and x 1, u y = wy = u
s











y



y

= y u x = y wx

= y usx = y wsx = 0 for y 0 and y 2

( c c







dc c ccr x c xT for T T0 dT

or, in the other notation of the axes,
Copyright © 2010 SciRes.

and
AM


176

E. A. LUKASHOV ET AL.

x F c, T = 3 c ccr



2

dccr x c dT xT for T T0 ;

dc =m dT

1
2 clust

4

clust

(

clust

T

,

dccr c0 = dT T0 Similarly, for y = 0 and ary condition on c takes the

= cmid Tmin y = 2 the second boundform

methods of [16,18], performed a numerical analysis of the model. The numerical results are reproduced here under their kind permission. Figures 4(a-d) and 5(a-d) illustrate the numerical results and show a complicated dynamics of the crystallization process. To test the crystallization model (82)-(89), the following dimensionless values of the main parameters were taken on the basis of their physical sense:

y * = 0.

= 6 , 7 3; R = 0, 7 4; q = 0, 2 5;
Q = 1, 75; M = 6, 2; = 1 103 ; g = 0; s = 1; D0 = 1; = 7,2; = 1 10 4 ;
F2 ( x, y ) 0; F1 ( x, y ) = 1 103.

4) For the temperature T we impose the initial conditions T 0, x, y = Tmin Tmax Tmin y yend , yend = 2 and the boundary conditions
T t , 0, y = T t ,1, y = Tmin Tmax Tmi nT |
y =0, y = yend n



t 1 ,

= 0,

where Tmin = 233.15K , Tmax = 456.15K , = 100 is a parameter.

4.6. Numerical Analysis of the Model of Oriented Crystallization, Two-Dimensional Case
N. A. Zaitsev, Yu. G. Rykov, and V. Lysov, based on the

Time-development of crystallization process in the case of isotropic surface tension of crystal grains. The banding structure is transformed to the equiaxial structure. (a) t = 4 , (b) t =8 , (c) t =18 , (d) t = 22 (Figure 4). Figures 4(a-d) presents the situation where the surface tension of crystal grains is isotropic. In this case, the chemical potential has the form

* = F c, T 2 c.

(91)

The banding structure is formed at initial times and

(a)

(b)

(c)

(d)

Figure 4. Time-development of crystallization process in the case of isotropic surface tension of crystal grains. The banding structure is transformed to the equiaxial structure: (a) t = 4; (b) t = 8; (c) t = 18; (d) t = 22.
Copyright © 2010 SciRes. AM


E. A. LUKASHOV ET AL.

177

(a)

(b)

(c)

(d)

Figure 5. Time-development of crystallization process in the case of anisotropic surface tension of crystal grains (formation of a dendrite liquid-solid damain): (a) t = 1; (b) t = 2; (c) t = 4; (d) t = 8.

then is developed to an equiaxial like structure. There is a certain analogy with the crystallization of eutectics when one of the phases splits into small cells [19] or the dendrite growth [20] when the distance between secondary branches of dendrites in a perfectly solidified mass is much larger than that at initial times. We also can imagine a similar situation where a jet breaks down into drops when the absolute value of the surface tension is rather large. Time-development of crystallization process in the case of anisotropic surface tension of crystal grains (formation of a dendrite liquid-solid domain): (a) t = 1 , (b) t = 2 , (c) t = 4 , (d) t = 8 (Figure 5). Figures 5(a-d) illustrates the numerical experiment in the case of an anisotropic surface tension, In this case, the following formula is used instead of (91):

developed to the dendrite structure.

4.7. Conclusions
The aforesaid shows that the proposed mathematical model of crystallization can be viewed as a mathematical reconstruction of various experiments. In particular, the following result of numerical analysis agrees with experimental observations: it is seen in Figures 4(a), (b) and Figures 5(a), (b) that the banding structure is the first structure formation to appear in the instability zone and the subsequent reformation of structure proceeds because of arising waves similar to the Marangoni instability wave at the boundaries of bands (cf. [2,4,5,21,22]). We also note that the numerical results concerning the influence of the crystallographic orientation of a growing crystal on the structure formation in alloys (cf., for example, [19]) can be regarded as confirmation of the verifiability of our mathematical crystallization model.

* = F c , T 2 x A A A
A =1, 5 E;







c c c , x c c ;

1 0 A = 0 14

5. Acknowledgements
The work was financially supported by the Russian Foundation for Basic Research (grants No. 09-01-12024
AM

where E is the identity matrix. In this case, the banding structure, deformed because of overcrystallization, is
Copyright © 2010 SciRes.


178

E. A. LUKASHOV ET AL. Stefan Problem," In: Free Boundary Problems Involving Solids, Longman Sci. Techn., Harlow, 1993, pp. 135-142. [11] O. A. Oleynik and E. V. Radkevich, "Second Order Equations with Nonnegative Characteristic Form," Am. Math. Soc., Providence, 1973. [12] A. A. Lacey and A. B. Tayler, "A Mushy Region in a Stefan Problem," IMA J. Appl. Math., Vol. 30, No. 3, 1983, pp. 303-313. [13] M. A. Biot, "Mechanics of Deformation and Acoustic Propagation in Porous Media," Journal of Applied Physics, Vol. 33, No. 4, 1962, pp. 1482-1498. [14] S. J. Watson, F. Otto, B. Y. Rubinstein and S. H. Davis, "Coarsening Dynamics for the Convective Cahn-Hilliard Equation," University of Bonn, Preprint, 2003. [15] A. A. Golovin, S. H. Davis and A. A. Nepomnyashchy, "A Convective Cahn-Hilliard Model for the Formation of Facets and Carners in Crystal Growth," Physical D, Vol. 118, 1998, pp. 202-230. [16] N. A. Zaitsev and Yu. G. Rykov, "Numerical Analysis of a Model Describing Metal Crystallization I. One-Dimensional Case," Preprint, Keldysh Institute of Appl. Math. RAS, No. 72, 2007. [17] W. Dreyer and B. Wagner, "Sharp-Interface Model for Eutectic Alloys. Part I. Concentration Dependent Surface Tension," Preprint, 2003. [18] N. A. Zaitsev and Y. G. Rykov, "Numerical Analysis to a New Model of the Matel Solidification, 1-D Case," Mathematical Simulation, 2010, to Appear. [19] N. P. Lyakishev and G. S. Burkhanov, "Metal Monocrystals [in Russian]," Eliz, Moscow, 2002. [20] P. E. Shalin et al., "Monocrystals of Heat-Resistance Nickel Alloys [in Russian]," Mashinostroenie, Moscow, 1997. [21] J. W. Matthews and A. E. Blacslee, "Defects in Epitaxial Multilayers. I. Misfit Dislocations," Journal of Crystal Growth, Vol. 27, No. 1, 1974, pp. 118-125. [22] V. N. Vigdorovich, A. E. Vol'pyan and G. M. Kurdyumov, "Oriented Crystallization and Physic-Chemical Analysis [in Russian]," Chemistry, Moscow, 1976.

and 09-01-00288).

6. References
[1] V. A. Avetisov, "p-Adic Description of Characteristic Relaxation in Complex System," Journal of Physics A: Mathematical and General, Vol. 36, No. 15, 2003, pp. 4239-4246. E. N. Kablov, "Cast Gas-Turbine Engine Blades (Alloys, Technology, Covering) [in Russian]," Moscow, 2001. A. N. Kolmogorov, "On the Statistical Theory of Crystallization of Metals [in Russian]," Izv. Akad. Nauk SSSR, Ser. Mat., No. 3, 1937, pp. 355-359. F. M. Shemyakin and P. F. Mikhalev, "Physic-Chemical Periodic Processes [in Russian]," Akad. Nauk SSSR, Moscow, 1938. I. Z. Bezbakh, B. G. Zakharov and I. A. Prokhorov, "Radiographical Characterization of Microsegregation in Crystals [in Russian]," In: Proceedings of the 6th International Conference "Growth of Monocrystals and Heat-Mass Transfer", Vol. 2, Obninsk, 2005, pp. 352361 V. G. Danilov, G. A. Omel'yanov and E. V. Radkevich, "Hugoniot-Type Conditions and Weak Solutions to the phase Field System," European Journal of Applied Mathematics, Vol. 10, 1999, pp. 55-77. E. V. Radkevich, "Mathematical Aspects of Nonequilibrium Processes [in Russian]," Tamara Rozhkovskaya Publisher, Novosibirsk, 2007. N. N. Yakovlev, E. A. Lukashev and E. V. Radkevich, "Problems of Reconstruction of the Process of Directional Solidification [in Russian]," Dokl. Akad. Nauk, Ross. Akad. Nauk, Vol. 421, No. 5, 2008, pp. 625-629; English translation: Dokl. Phys., Technical Physics, Vol. 53, No. 8, 2008, pp. 442-446. V. Visintin, "Models of Phase Transitions," Birkhauser, Boston, 1996.

[2] [3]

[4]

[5]

[6]

[7]

[8]

[9]

[10] E. V. Radkevich, "The Gibbs--Thomson Effect and Existence Conditions of Classical Solution for the Modified

Copyright © 2010 SciRes.

AM