ANALYSIS OF A MATHEMATICAL MODEL RELATED TO CZOCHRALSKI CRYSTAL GROWTH

This paper is devoted to the study of a stationary problem consisting of the Boussinesq approximation of the Navier–Stokes equations and two convection–diffusion equations for the temperature and concentration, respectively. The equations are considered in 3D and a velocity–pressure formulation of the Navier–Stokes equations is used. The problem is complicated by nonstandard boundary conditions for velocity on the liquid–gas interface where tangential surface forces proportional to surface gradients of temperature and concentration (Marangoni effect) and zero normal component of the velocity are assumed. The velocity field is coupled through this boundary condition and through the buoyancy term in the Navier– Stokes equations with both the temperature and concentration fields. In this paper a weak formulation of the problem is stated and the existence of a weak solution is proved. For small data, the uniqueness of the solution is established.


Introduction
In this paper we investigate the solvability of a stationary three-dimensional mathematical model describing the processes in the melt during a silicon single crystal pulling by the Czochralski method.The main feature of the Czochralski method (cf.e.g.[3,13,15]) is that the grown single crystal is pulled from its melt which is situated in a crucible (cf.Fig. 1).The device is rotationally symmetrical and, during the pulling, both the crucible and the crystal perform rotational motions around the symmetry axis and, at the same time, motions in the vertical direction which correspond to the crystal growth velocity and maintain the melt free surface in a constant position.The crystal growth velocity is usually considerably smaller than the typical Cross-section through the melt and crystal in a Czochralski apparatus.
velocities in the melt.The crucible is made of vitreous silica which leads to a contamination of the melt by oxygen.
The region occupied by the melt is assumed to be known a priori and will be denoted by Ω in the following.The boundary ∂Ω of Ω comprises the crucible wall Γ W , the melt free surface Γ LG and the interface Γ LS between the melt and the growing crystal (cf.Fig. 1).The mathematical model then consists of the following partial differential equations defined in Ω and boundary conditions prescribed on ∂Ω: We use the following notations: v is the velocity, p is the pressure, θ is the temperature, c is the oxygen concentration, I is the identity tensor, n is the unit outward normal vector to ∂Ω, and ∇ s denotes the surface gradient (in this paper, ∇ s • = (I − n ⊗ n) (∇•)| ∂Ω ).The constants α 1 , . . ., α 7 and the functions f 1 , f 2 , ϕ 1 , ϕ 2 , v b , θ b , c b will be specified in the next section.
The equations (1.1)-(1.4)can be derived from general balance laws for linear momentum, mass, energy, and mass of a constituent, respectively.The equations (1.1) and (1.2) represent the Boussinesq approximation of the Navier-Stokes equations.Among the boundary conditions, the most interesting one is the boundary condition (1.6) describing the so-called Marangoni effect, i.e. the fact that surface tension variations due to temperature and concentration gradients induce tangential surface forces on the melt free surface.A detailed derivation and explanation of the model (1.1)-(1.8)can be found in [8].The solvability of a rotationally symmetrical case of (1.1)-(1.8)defined in a simplified geometry and taking into account magnetic forces was investigated in [4].The formulation of (1.1)- (1.8) used in [4] requires to apply weighted Sobolev spaces, but on the other hand, it considerably simplifies the treatment of (1.2) and (1.6).
Remark 1.1.The constants α 1 , . . ., α 7 have the following physical meanings: Re Pr , where Re is the Reynolds number, Pr is the Prandtl number, Sc is the Schmidt number, Gr is the Grashof number (Gr θ for temperature, Gr c for concentration), and Ma is the Marangoni number (Ma θ for temperature, Ma c for concentration).The usual definitions of the functions f 1 , f 2 , ϕ 1 , ϕ 2 are for θ and c belonging to some bounded intervals (θ 1 , θ 2 ) and (c 1 , c 2 ), respectively.Outside these bounded intervals, the functions f 1 , f 2 , ϕ 1 , ϕ 2 have no physical meaning and can be defined arbitrarily.In the above relations, e g is a unit vector in the direction of the gravity, Bi is the Biot number, Rd is the radiation number, and k 0 is the segregation coefficient.
Let us mention the most important difficulties arising when investigating the weak solvability of (1.1)- (1.8).The first difficulty comes from the fact that full Dirichlet boundary conditions are prescribed only on a part of the boundary whereas, on the remaining part, only the normal component of the velocity is known a priori.That also causes difficulties in case of the Navier-Stokes equations since all methods used in the literature for proving their solvability are based on the assumption that there exists a divergencefree function v b satisfying the Dirichlet boundary conditions for the velocity such that the integral Ω v b •(∇ v) v dx is, in some sense, small (cf.e.g.(4.29), (4.35)).It is known that such assumption is fulfilled if full Dirichlet boundary conditions for the velocity are prescribed on the whole boundary and the flux through each connected component of the boundary vanishes (cf.[7,p. 287,Lemma 2.3]).However, in the case considered here, the existence of a suitable function v b is an open problem.That means that, in case of nonhomogenous mixed boundary conditions of the mentioned type, which often occur in various applications, it is not known how to prove the solvability of the stationary incompressible Navier-Stokes equations in general and it is rather surprising that such a fundamental problem still remains unsolved.
Another difficulty consists in the nonstandard boundary condition (1.6) for the velocity on Γ LG .Here, a suitable generalization in form of a linear functional defined on a subspace of H 1 2 (∂Ω) 3 has to be constructed for the occurring surface gradients of temperature and concentration.This generalization should be also appropriate for a numerical solution of (1.1)-(1.8)by means of the finite element method.
Finally, the investigations of (1.1)-(1.8)are complicated by the coupling between the equation (1.1) and the equations (1.3) and (1.4) which is realized through the buoyancy terms in (1.1) and through the boundary condition (1.6).In Theorem 4.5, we shall prove that the solvability does not depend on the magnitude of the constants in the coupling terms.Therefore, it suffices to prove the solvability only for those cases when these constants are small, which is, of course, much easier.
The plan of the paper is as follows.In Section 2, we formulate assumptions on the problem (1.1)-(1.8)and introduce some notations.In Section 3, we construct the mentioned generalization of the surface gradient, derive a weak formulation and show the equivalence between classical solutions and smooth weak solutions.In Section 4, we establish an equivalent operator formulation and prove the solvability of the weak formulation by applying the Leray-Schauder principle.The Leray-Schauder principle allows us to make a weaker assumption on the Dirichlet boundary condition v b for the velocity than usually made in the literature.Further, we prove that, for small data, the weak solution is unique.
Remark 2.1.The assumption on the existence of the extension m is made because of the treatment of the surface gradients in the boundary condition (1.6).It can be shown that this assumption is satisfied if Γ LG is a C 1,1 surface and if there exists a finite number of local Cartesian coordinate systems providing a description of ∂Ω (cf.[10, p. 14]) with the property that, in each of these coordinate systems, the projection of the respective part of Γ LG into the (x 1 , x 2 )-plane is a set with a Lipschitz-continuous boundary.In this case, we even have m ∈ W 1,∞ (IR 3 ) 3 .Another sufficient condition for the existence of m ∈ W 1,∞ (IR 3 ) 3 is the existence of a domain Ω with a C 1,1 boundary satisfying Γ LG ⊂ ∂ Ω.
We make the following assumptions on the data of the problem (1.1)-(1.8): ) ) ) Generally, the condition v b • n ≥ 0 on Γ LS is a technical assumption needed for proving both the existence and the uniqueness of the weak solution.In the case of the Czochralski method, however, this assumption is a natural condition which expresses that the crystal is really growing and not melting.

Weak Formulation
To derive a weak formulation of (1.1)-(1.8),we first assume that the functions v, p, θ and c are a classical solution of our problem.We multiply the equations (1.1)-(1.4)by arbitrary functions , respectively, integrate them over Ω, use the identities apply the Gauss integral theorem (cf.[5, p. 33]) and substitute the Neumann boundary conditions and the condition (1.2).Then we obtain ) ).However, this new formulation makes it possible to introduce more general solutions in a natural way.In fact, all terms with the exception of the two last ones in (3.5) are well defined for functions v, w, p, λ, θ, η, c, q belonging to the Sobolev spaces H 1 (Ω) 3 , L 2 (Ω) and H 1 (Ω), respectively.Therefore, it remains to generalize the surface integrals of the type [10, p. 72, Theorem 3.8] and [10, p. 69, Theorem 3.4]), the terms in the square brackets in (3.9) are elements of L 2 (Ω) 3 and hence the mapping d is well defined.Using the Hölder inequality, we obtain and hence it follows from the Sobolev imbedding theorems that which implies the continuity of d.
Let us assume for a moment that m ∈ C ∞ (Ω) 3 .Then we have for any According to the trace theorems (cf.[10, p. 84, Theorem 4.2]), the righthand side of (3.13) is defined for any and H 1 (Ω), we obtain the property (3.10).Now, due to (3.10), we have The constant Thus, we have which means that the mapping d s is continuous.By (3.10), we obtain Using the density of H 2 (Ω) into H 1 (Ω), we observe that, for any ζ ∈ H 1 (Ω) and w ∈ W, the value <d s (ζ), γ(w)> is independent of the choice of m.Since, for any w ∈ W, we have n According to the Friedrichs inequality (cf.[10, p. 20, Theorem 1.9]), | • | 1,Ω is a norm on W, equivalent to • 1,Ω , and hence the constant C s is finite by (3.12).
Thus, we see that the functional d s (ζ) is a reasonable generalization of the surface gradient on Γ LG for functions ζ ∈ H 1 (Ω) and hence, replacing the last two integrals in (3.5) by −α 3 <d s (θ), γ(w)> − α 4 <d s (c), γ(w)>, we can define a weak formulation of (1.1)-(1.8).First, however, let us introduce the following notations: Then the functions v, p, θ, c are a weak solution of the problem (1.1) and The reason is that, for a plane Γ LG , any functions v, w ∈ C 2 (Ω) 3 with v • n = w • n = 0 on Γ LG satisfy w • (∇v) T n = 0 on Γ LG , which makes it possible to use the identity instead of the identity (3.1).Let us mention that, for a plane Γ LG with a normal vector n, the formula (3.9) can be simplified to The following two theorems show that the weak solution really is a meaningful generalization of the classical solution.
Since C ∞ (Ω) 3 is a dense subspace of H 1 (Ω) 3 , we deduce that (3.24) holds for any w ∈ H 1 (Ω) 3  2 Particularly, (3.25) holds for any w ∈ C ∞ 0 (Ω) 3 with a vanishing right-hand side and since the terms in the square brackets are continuous, we infer that the functions v, p, θ, c fulfil the differential equation (1.1) in the classical sense.Then it follows from (3.25) that

Existence and Uniqueness of the Weak Solution
In this section, we investigate the existence and uniqueness of the weak solutions of the problem (1.1)-(1.8).First, in Theorem 4.1, we show that the pressure can be eliminated from the weak formulation (3.19)-(3.23)and we can confine ourselves to investigations for the functions v, θ and c.Then we construct an operator formulation which enables to perform a proof of the weak solvability for small values of the constants α 1 , . . ., α 4 applying the Leray-Schauder principle.Using a simple scaling argument, we extend this existence result to arbitrarily large constants α 1 , . . ., α 4 .Finally, we show that the weak solution is unique for small data. and Then there exists a unique function p ∈ L 2 0 (Ω) such that the functions v, p, θ, c are a weak solution of the problem (1.1)- (1.8).
Lemma 4.1.The spaces V, Θ and C are separable Hilbert spaces for the scalar products a 1 (•, •) and a 2 (•, •), respectively, and the space X is a separable Hilbert space for the scalar product where U = (v, θ, c), W = (w, η, q) and U, W ∈ X.The norm U X ≡ (U, U) X is equivalent to the norm v 1,Ω + θ 1,Ω + c 1,Ω , i.e., there exist constants C 1 and C 2 depending only on Γ W , Γ LG and Ω such that For any sequence Proof.Using the Friedrichs and Korn inequalities (cf.[10, p. 20, Theorem 1.9] and [11, p. 97, Lemma 3.1]), we deduce that the bilinear forms a 1 and a 2 are scalar products on the respective spaces and that the norms induced by these scalar products are equivalent to the norm • 1,Ω .Thus, for proving that the spaces V, Θ and C are Hilbert spaces for the mentioned scalar products, it suffices to show that these spaces are closed subspaces of H 1 (Ω) 3 and H 1 (Ω), respectively, with respect to the norm • 1,Ω .Consider any u ∈ H 1 (Ω) 3 and assume that there exists a sequence which implies that div u = 0 (cf.[10, p. 56, Proposition 1.1]).The continuity of the trace operator gives γ(u n ) → γ(u) in L 2 (∂Ω) 3 and also γ( Since H 1 (Ω) 3 ⊂ V and H 1 (Ω) ⊂ Θ , C and since, by the Hahn-Banach theorem (cf.[12, p. 186, Theorem 4.3-A]), any functional belonging to V , Θ or C can be extended to a functional belonging to H 1 (Ω) 3 or H 1 (Ω) , respectively, we obtain the equivalence (4.8).
In the following lemma, we shall use the fact that, according to the Sobolev imbedding theorems and the trace theorems (cf.[10, pp.69 and 84]), there exist finite constants which depend only on Ω.
Lemma 4.2.The trilinear form b 1 is continuous on H 1 (Ω) 3 3 and satisfies the inequalities The trilinear form b 2 is continuous on H 1 (Ω) 3 ×H 1 (Ω)×H 1 (Ω) and satisfies the inequalities Proof.Using the imbedding H 1 (Ω) ⊂ L 4 (Ω) (cf.[10, p. 69]) and applying the Hölder inequality, we obtain (4.9).For proving the inequalities (4.10) and (4.11), we use the identities and the Gauss integral theorem.In case of the trilinear form b 2 , we can proceed in the same fashion.Lemma 4.3.Choose any U, U, W ∈ X, U = (v, θ, c), U = ( v, θ, c), W = (w, η, q), and denote Then there exists a constant C independent of U, U and W such that where Proof.According to the Taylor formula and (2.2), we have and hence it follows from the Hölder inequality that the mappings F 1 and F 2 are well defined and that Further, using (2.4) and (2.5), we infer that the mappings Φ 1 and Φ 2 are well defined too and that Thus, applying the trace theorems and Lemma 4.2, we obtain (4.15).Using (2.2)-(2.5), the Taylor formula and the Hölder inequality, we deduce that The operator A is compact, strongly continuous and satisfies for any U ≡ ( v, θ, c) ∈ X and any W ≡ (w, η, q) ∈ X where C is independent of U and W. The operator L is linear and continuous and satisfies for any U ∈ X Then, according to Lemmas 4.3, 3.1 and 4.1, we have a(U, •), l(U, •) ∈ X for any U ∈ X and hence it follows from the Riesz representation theorem that there exist uniquely determined operators A, L : X → X satisfying Clearly, the operator L is linear and the equivalence (4.24) holds.Consider any sequence Then, according to (4.16) and (4.7), [10, p. 106, Theorem 6.1]), the compactness and linearity of the trace operator γ : H 1 (Ω) → L 2 (∂Ω) (cf.[10, p. 107, Theorem 6.2]) and the fact that any weakly convergent sequence is bounded, we deduce that A(U n ) → A(U) in X.Thus, the operator A is strongly continuous and since the space X is (as a Hilbert space) reflexive, the operator A is compact.
The inequalities (4.27) and (4.28) follow using Theorem 3.1 and the inequality (4.7).Since L is linear, its continuity is a consequence of (4.27).
Since the solutions of the problem (4.2)-(4.5)do not depend on the particular choice of the functions v b , θ b and c b , we can use the following lemma to reduce the influence of the second term on the right-hand side of (4.26).For proving the solvability of the problem (4.2)-(4.5) in case of small constants α 1 , . . ., α 4 , we shall apply the Leray-Schauder principle which is formulated in the following theorem.
Then there exists a positive constant M such that the problem (4.2)-(4.5) is solvable for any constants α 1 , α 2 , α 3 , α 4 ∈ (0, M) and α 6 , α 7 ∈ IR + and for any functions Proof.Consider any constants α 5 , α 6 , α 7 ∈ IR + and any functions ) and (4.29).Let ε be the constant from the assumption (4.29) and choose any ϑ ∈ (ε, 1).Set and consider any constants α 1 , α 2 , α 3 , α 4 ∈ (0, M).According to [7, p.  Further, let A, L : X → X be the operators from Theorem 4.2.According to (4.27), we have (2 + ϑ) U X ≤ 3 L U X for any U ∈ X and since L is linear, there exists a continuous linear inverse operator L −1 .Let us assume that the problem (4.2)-(4.5) is not solvable.Then, due to the fact that the operator L −1 A is compact, it follows from Theorems 4.2 and 4.3 that the solutions of the equations U = −λ L −1 A(U), λ ∈ (0, 1), form an unbounded set.Thus, there exist sequences {λ n } ∞ n=1 ⊂ (0, 1) and and U n X > n.Denoting U n = ( v n , θ n , c n ), we obtain by (4.28) and (4.26) Let us set U n = U n / U n X and denote U n = ( v n , θ n , c n ).Since any Hilbert space is reflexive and the sequences {λ n } ∞ n=1 and { U n } ∞ n=1 are bounded, we can assume without loss of generality that there exist λ 0 ∈ <0, 1> and U ≡ ( v, θ, c) ∈ X such that λ n → λ 0 and U n U in X.Clearly, U X ≤ 1 (cf.[12, p. 209]).By (4.31), we have and hence, applying (4.8), (4.10) and the arguments used for proving the strong continuity of A in the proof of Theorem 4.2, we obtain as a limit of (4.32) Consider any w ∈ V and set W = (w, 0, 0).Then W ∈ X and it follows from (4.30) that, for any Since the first term on the left-hand side converges to zero, we deduce from (4.25) that λ n b 1 ( v n , v n , w) → 0. Thus, using the same arguments as above, we obtain which means that the inequality (4.33) holds for any v b satisfying the conditions (4.1).This is, however, not possible since, according to the assumption (4.29), there exists v b satisfying (4.1) and the inequality Therefore, the problem (4.2)-(4.5) is solvable.
Remark 4.2.The Leray-Schauder principle was already used by Ladyženskaja [9] for proving the weak solvability of the stationary incompressible Navier-Stokes equations.However, the proof given in [9] is not correct.The substantial shortcoming is an incorrect argumentation why the inequality Then, according to (4.26) and (4.28), the operator A ≡ A + L satisfies

the set of all the solutions of the problem (4.2)-(4.5). Then we have for any positive real numbers
The functions f 1 , f 2 , ϕ 1 , ϕ 2 satisfy the relations (2.2)-(2.5).
Proof.The statement immediately follows by writing θ and c as θ (θ/ θ) and ĉ (c/ĉ), respectively, in (4.2)-(4.5).Now, as a consequence of the two above theorems, we obtain the main result of this paper.Defining the functions f 1 , f 2 , ϕ 1 , ϕ 2 as in Theorem 4.5, it follows from Theorem 4.4 that the set is nonempty.Thus, according to Theorem 4.5, the set is nonempty as well and hence the problem The remainder of the paper is devoted to investigations of the uniqueness of the weak solution.Similarly as for the Navier-Stokes equations, we shall be able to prove the uniqueness for small data.First, let us introduce some notations.In consequence of the trace theorems and the Friedrichs inequality (cf.[10, pp. 84 and 20]) there exist finite constants depending only on Γ W , Γ LG and Ω.Further, we denote For proving the uniqueness of a weak solution, we shall need certain a priori estimates which we establish in the following lemma.

3
be an arbitrary function satisfying γ(w) = z and let us set <d s (ζ), z> = d(ζ, w) .(3.14) Then (3.14) defines a continuous linear mapping d s : H 1 (Ω) → γ( W) which does not depend on the choice of the extension m in the definition of the mapping d.Moreover, for ζ ∈ H 2 (Ω), we have

Remark 3 . 1 .Remark 3 . 2 .Remark 3 . 4 .
.23) The weak solution does not depend on the particular choice of the functions v b , θ b and c b satisfying (3.17).In view of (3.21) and (2.7), we have b(v, λ) = 0 ∀ λ ∈ L 2 (Ω) and hence any weak solution satisfies the condition div v = 0. Remark 3.3.Since the pressure p is determined by (1.1), resp.by (3.20), up to an arbitrary additive constant, we consider only pressures with zero mean value.If Γ LG is plane, then, in (3.20), the bilinear form a 1 can be replaced by Ω ∇v • ∇w dx .

3 ,
where we used the fact that (I − m ⊗ m) w ∈ W ∀ w ∈ W and (I − n ⊗ n) ∇ s = ∇ s .Again, the terms in the square brackets are continuous and hence we deduce that the Neumann boundary condition (1.6) is also fulfilled in the classical sense.The validity of the equations (1.3) and (1.4) and of the Neumann boundary condition in (1.7) can be proven in the same fashion.For proving the Neumann boundary condition in (1.8), we have to apply Proposition 1.1 from [10, p. 56], since n| Γ LS is not continuous in general.The fulfilment of the Dirichlet boundary conditions immediately follows from (3.19).

( 4 .Remark 4 . 3 .
33) is valid for any extension v b with the same function v. Relations analogous to (4.33) and (4.34) can be also obtained when proving the weak solvability of the stationary incompressible Navier-Stokes equations using the contradiction argument developed by Leray (cf.[6, p. 55]).Leray himself considered only problems with Dirichlet boundary conditions v b prescribed on the whole boundary ∂Ω and therefore he could use the fact that v ∈ V. Assuming, in addition, that the flux of v b through each connected component of ∂Ω vanishes, he then showed that, in consequence of (4.34) (with λ 0 = 0 and V replaced by V), any divergence-free extension v b of v b satisfies b 1 ( v, v, v b ) = 0.In case of mixed boundary conditions for the velocity of the type we deal with, a proof of such a conclusion still remains an open problem.

Remark 4 . 4 .
Let ϑ, α 1 , . . ., α 7 and f 1 , f 2 , ϕ 1 , ϕ 2 , v b , θ b , c b be defined as in the proof of Theorem 4.4 and let us assume that the assumption (4.29) holds for any v ∈ V with the same function v b , i.e., that
[8] Theorem 4.4 can be proved by means of Galerkin's method.The validity of (4.29) and (4.35) is known if Γ LG = Ø and the flux of v b through each connected component of ∂Ω vanishes (cf.[7, p. 287,  Lemma 2.3]).Unfortunately, if meas 2 (Γ LG ) = 0, such general results are not available.Of course, the assumptions (4.29) and (4.35) are satisfied if α 5 and v b are sufficiently small but, in practical cases, this requirement is usually too restrictive.Some sufficient conditions for the validity of (4.29) and (4.35) were derived in[8], however, their fulfilment is generally not easy to verify.Nevertheless, it can be shown that, for a rotationally symmetrical domain Ω corresponding to the Czochralski method (cf.Fig.1) and for v b representing rotational motions of the crucible and the crystal, assumptions (4.29) and (4.35) always hold.This result is very important in context of the Czochralski method since the non-rotational components in the Dirichlet boundary conditions for the velocity can be mostly neglected and it is sufficient to consider them only in the Neumann boundary condition(1.8).That leads to a problem which is solvable under the assumptions from Section 2.For proving the solvability of (4.2)-(4.5)for large values of α 1 , . . ., α 4 , the following property of the solution set is essential.For arbitrary positive real numbers α 1 , . . ., α 7 and arbitrary functions f 1