The steady state heat equation with mixed nonlinear boundary conditions—an example in crystallography

COMPUTER METHODS IN APPLIED h4ECHANICS NORTH-HOLLAND PUBLISHING COMPANY

THE STEADY BOUNDARY

AND ENGINEERING

STATE HEAT EQUATION CONDITIONS-AN

WITH MIXED

EXAMPLE

35 (1982) 47-66

NONLINEAR

IN CRYSTALLOGRAPHY

P. WITOMSKI Laboratoire

d’lnforrnatique

et de Mathkmatiques Appliquges de Grenoble, ZMAG, Grenoble Cedex, France Received Revised

manuscript

5 October

1981

received

4 January

Grenoble, F-38041,

1982

This article deals with the resolution of a nonlinear boundary value problem arising from a crystaliisation experiment. The model developed requires that the heat equation be solved in an open bounded set 0 of R3 with mixed non-linear boundary conditions. The existence of a unique solution is shown using monotone operators. This same result can be found by a constructive method using lower and upper solutions. The numerical analysis of the problem is presented. The axial symmetry of revolution of R is allowed for by the use of rotationally invariant finite elements which allow the problem to be reduced to one in two dimensions only.

1. The physical problem

The model treated arises from an experiment in crystal growth by the Czochralski method [2,9]. This is a method in which single crystals are pulled from a liquid melt. The material to be crystallised is heated above its melting point in a crucible and a monocrystalline germ is placed on the surface of the melt using a vertical holder. The crystallisation occurs around this germ. A pulling mechanism allows the holder to be raised as the crystal forms. The final product is a cylindrical monocrystal. The form of the liquid-solid interface is an important factor in the quality of the monocrystal. We have previously proposed a model which simulates the growth of the monocrystal and which allows the form of the interface to be determined [16]. We calculate here the heat distribution in the melt and in the solid crystal for a stationary state of the interface. 1.1.

The system of equations for thermal exchange

Fig. 1 shows a schematic diagram of the oven in which the monocrystal is formed. The experimental conditions allow the heat equation to be taken in a stationary form in 0, and 0, (see [16]): AT=0

in 0, and L&,

(1)

where T is the temperature. The boundary conditions arise from heat exchange by radiation (term in T4), free convection (term in (T - T,)l.=) and by the crucible temperature [15]. They fall in two classes: 00457825/82/0000-0000/$02.75

@

1982 North-Holland

48

P. Witomski,

Heat equation

with nonlinear boundary conditions

crysta rface

Fig. 1. Schematic

diagram of the crystal growth model.

(a) Dirichlet conditions:

T/n = Tc, the crucible temperature;

(2)

Tl,-, = T,,

(3)

and the melting temperature;

(b) Neumann conditions: -A

dT

‘dn

=

a(T

-

TJ'."

+

p,T4

-

TJ'=

+

/_L~T~-

-

(4)

~1 ;

r3

and -A

aT

’ an

=

a(T

7s

;

where cr, A,, A, are positive constants, T, is the ambient temperature positive bounded functions (see [l&16] for more details). The conditions at the liquid-solid interface TI are T = Tf ,

(5)

t-2

and p,, pl, r),, T,+are

(6) (7)

P. Witomski, Heat equation with nonlinear boundary conditions

I2

interface r;

m2

Fig. 2. The liquid-solid

interface.

where f represents the pulling speed and 8 the angle between the normal q1 to .R1 and the Oz.-axis (see Fig. 2). Also r/s = - nl. The system of equations (l)-(7), represents a free surface problem depending on the parameter f. A supplementary condition is imposed on the system: the pulling speed must be such that the interface T1 is always at m2 (see Fig. 2). This condition completes the formulation of the problem. The solution of this model is carried out by iteration on the interface. At each iteration step, for a fixed interface, two non-linear problems of the same type are solved: (PI) consisting of (1) (3) (5) in J&; (P2) consisting of (1) (2) (4) in &. We propose to solve P1 which, after change of variables can be written:

IAT=0

pI

Tb,=O,

I -A

and redesignated

2. Mathematical

inR,,

dT

I

dn

=p(T+

Tf)4+a(T+

Tf- T,)‘.25-q

I-2

Bi, as shown.

formulation

of the problem

We first recall some useful results (see [l,

101).

2.1. The trace operator and Green’s formula In what follows, 0 represents a regular bounded open set of W3 with boundary l-‘. Let yo: u + ujr be the trace operator on %‘(fi). y. extends into a bounded linear operator from

50

P. Witomski,

Heat equation

with nonlinear boundary

conditions

H’(R) onto H*‘*(r) (see [lo, 171 for the definition yl : u + (du/&z)l, linear operator duality between

of the spaces). For u E s’(o) the operator associates with u its normal derivative to IY y1 is extended into a bounded from H’(0, A) = {u E H’(R), Au E L*(R)} into H-l’*(r). Introducing the H-“*(r) and H”*(r), written ( , ),- Green’s formula becomes

tlu~H’(0)

I

n

VuEH’(R,A),

and

Av.udn+

Vu - Vv dR = (ylv,

you>r .

Now let rl be a regular open set of fl of non-zero measure. For h E L*(r), Z&h represents the function on L*(r) equal to h on T’l and zero elsewhere. Let i be the injection of H”*(r) into L*(r). Z,i takes values in H”*(r) and the projection operator ul so defined is continuous. We put a2 = 1 - pl and let a? represent the adjoint of gl (similarly a: is the adjoint of a*). 2.2. A new formulation for 95 We now regard the original physical problem so as to be able to use the theory of monotone

as belonging to a class of non-linear operators [3].

2.2.1. Modification of the non-linear condition on r2 If m represents a point on T’, we have seen that the condition C?r at the end of Section 1 can be written

-

E

(m) = + @(u(m) + Tf - u,)‘.*’ + p(m)(u(m)

where we are using (Y, p, A and q, introduced Forxru,-TfandmEr2weset f(x, m) = (l/A)(a(x

f(u,g(x, m) =

functions

T,, m),

f(x, m), i f(O, m),

+ T,)4 - q(m))

+ T,)4 - q(m)).

and Tr > u, > 0. The function

xlu,-

x + f(x, m) increases

for

Tf,

x E [u, - T,, 01 , xro.

g(x, m) is extended to include the zero on R x (r\I’*). x + g(x, m) is monotone, bounded on R. An operator cp from L*(r) to L’(r) is associated with g and defined by:

&)(m> = s(u(m>,4

of problem

in (5) above.

+ Tr - u,)~.~’+ p(m)(x

p(m) and v(m) are positive all fixed m E r,. Put

in the statement

problems

increasing

and

51

P. Witomski, Heat equation with nonlinear boundary conditions

where m E r2 a.e. PROPOSITION

2.1. cp is a monotone Lipschitzian operator from L2(r) into L*(r).

PROOF. Let u E L2(r); since 7) and /.L are in L”(T) it follows that p(u) belongs to L2(r).

lg(x,m)-g(y,m)l4clx-yl

onmr

SO

IId4 - (PbJ)lh)22cllu- 4L2u-)~ The monotony fixed m.

of cp arises simply because g(u, m) is monotone

increasing with respect to u for

2.3. A class of mixed non-linear problems The condition on r2 is expressed in terms of the operator

q by

--&l&z = q(u). The experimental conditions impose upon 40 the hypotheses: (a) rp(O)lO; (b) V(U, - T,) 5 0 . The problem PI can be put into the form of the following more general problem. Find u E H’(R, A) such that -Au=O, 2% . UlYOU = 0, aZ(r1u

+ VPYOU)) = 0,

We now give the results concerning the uniqueness and existence of a solution for this class of problem. Two methods are described. The first uses monotone operators; the second which is based on the techniques of lower and upper solutions (see [13] for example) leads quite naturally to a numerical method. Refs. [7,8] should be consulted for descriptions of the use of optimisation methods.

3. The existence and uniqueness of a solution to the problem CPM.Monotone method This section depends essentially on the following theorem for monotone

operators.

THEOREM 3.1 (Surjectivity theorem, see [3]). Let V be a reflexive Banach space. Let A be an operator from V to its dual space V* having the following properties:

52

l? Witomski, Heat equation with nonlinear boundary conditions

(a) A is monotone, bounded and hemicontinuous; (b) Lirnllv+ (Av, v)l~~v~~ = + ~0, where v E V, then A is surjective from V onto V”. The non-linearity of the problem pose the problem on the space

w = {g

E

H”*(r),

a1g =

PM resides

in its Neuman

at the boundary.

We

0).

W, with the induced norm, is a Hilbert space. Let g E H1’2(r). Introduce u, as the solution -Au,=O,

condition

of’

YOU,= g .

The operator R : g + u, is an isomorphism We now introduce the operator T : f~I”~(r)

Tg = y&

from H1’2(r) + H-1’2(lJ

to the space {u E H’(0), defined by

Au = 0).

+ dyo&)

and A : W --, W* with values Ag = (i*aZ Ti)g , where

i is the injection

THEOREM

from

3.2 (Uniqueness

W into JP2(r). theorem).

Note that (Agl, g&

= (Tg,, g2),

5% has one and only one solution.

PROOF. We show first that the function A satisfies the hypothesis of Theorem 3.1. A is bounded and hemicontinuous because cp is Lipschitzian and continuous [16]. We have

6% - Aa, where (, monotone

jw denotes property

gl -

g&v = (Tgl - Tg,, gl

the duality of cp we find

(Agl - Ag,, gl - g,),

=

between

W*

-

g&

2

and

W. Using

Green’s

formula

and

the

- g&y

gd>llh .

Friedrichs’ inequality [ll] and the fact that H1’2(r) onto {u E H’(0), Au = 0} give Ag2, gl -

gdr

(ylR(gl - g2) + cp(gl) - cp(gz), gl

2 IlWWl

(Agl -

-

R is by construction,

an isomorphism

from

const.- llgl- g&w) .

A is therefore monotone as well as coercive. Then by Theorem that A(t) = 0. As a direct consequence, R(t) is the only solution

3.1 above, of CPM.

t E W exists such

P. Witomski,Heat equation with nonlinear boundary conditions

We now set out the variational formulation

53

of PM.

PROPOSITION

3.3. Let H:.,(R) = {u E H’(R), (TVyou = 0). Then, u is the solution of ??)Mif and only if u E H;,(R) and satisfies, for all v in H;,(a),

I

R

Vu * Vv dR + (&+,u),

yov)r = 0.

The proof is readily carried out using Green’s formula as in the linear case [16]. 4. Relationship bet,ween the problem PM and the original physical problem If we can show that the solution of the problem ??)Mtakes is values in the interval [u,- T,, 0] we will have (u + T,) as the solution to the original problem. The physical problem has thus a unique solution in the interval [u,, T,]. 4.1. Positiveness in H’(0) We reca_ll some results on the positivity of functions in H’(a) (see [14, 161). Let E be a subset of a. u E H’(a) is said to be non-negative on E in the sense of H’(0) if there exists a sequence {u,} functions from %‘(fi) such that (i) u,rOonE; (ii) u,,, + u, strongly in H’(0). We set v+ = sup(v, 0) and v- = -inf(v, 0). The properties of v+ and v- are examined in detail in [16]. In particular we have the three properties: (1) v+ and v- are in H’(0); (2) sn Vv’ ,Vv-dJ2=0; (3) YOU’* you- = (you)’ . (you)- = 0. PROPOSITION

4.1. Let u be the solution of 9,. For almost every x in Q u(x) E [u, - T,, 01.

PROOF. We show that u is negative by proving that u+ = 0. The variational (Proposition 3.3) of 9, with v = u+ and u = u+ - u- yields

I IVu+l”dfl+ R

(d(You)+),

(y&)+)r

= 0 .

Since a(O) 2 0 (see Section 2.3) and since 4p is monotone,

0= 2

Therefore

I Ivu’l’ dfi + n

lVu’(*

M(Yo~+)

formulation

the Friedrichs inequality gives

- 40(o),(rou)+h- + (40(o),(you)+&

dR 2 const. - ((u+(J~~~~~.

u+ = 0. In the same way it can be shown that u 2 u, - Tr.

54

P. Witomski,

Heat equation

with nonlinear boundary

conditions

Thus, u + Tf satisfies the equations set up in section 1 and since physically the solution is unique on L?,, u + Tf must be that solution. 5. A constructive

method of solution for PM: lower and upper solutions

This method leads to an algorithm for calculating the solution which is enclosed by a set of lower and upper solutions. The existence of a solution to PM is shown by the convergence of the method.- _ LEMMA

Let u E H’(O)

5.1 (Positiveness lemma). Au=O,

and A E W’ such that

a;(Ayou + y,u) 2 0 .

OIYOU2 0,

Then u 2 0. PROOF.

We show that u- = 0. Writing Green’s formula with u and u- we find 0=-

IR

Au-u-da=

JIVu-I*dfi n

ZZ-

IR

Vu . Vu- dR - (ylu, you-)r

(Ayou

+ ylu, you-)r + A(you, you-h.

But (Toyou 2 0 implies that you- = a2y0u-. It follows that

and using Friedrichs’ inequality this becomes 0L

Therefore

In

da 2 const. * llu-~~&~.

IVu-1*

U- = 0.

In order to resolve YM numerically we will construct a sequence of iterations which are solutions of linear problems. A choice must be made for the linearisation procedure. 5.1. Linearisation of the problem 5% We use the following linearisation

9L”

(71 yowl

aT(Ayou,

=

0

in order to construct a monotone

,

+ ~1%)

= (T;(-&&-I))+

Ayou,-,

sequence of iterations.

P. Witomski, Heat equation with nonlinear boundary conditions

where

A, real and positive,

is chosen

such that

- cp(you) + cp(you) + A(you - yov)20.

u 2 2, 3

It suffices to take A greater than or equal to the Lipschitz shows that the correspondence u,_~ + u, is monotone.

REMARK instead

5.2. The introduction of A amounts of -&+,/a,, = uzP1 (see [16]).

constant

to linearising

since Lemma

in the form

-&+/a,,

5.1 clearly

= u&-~

5.2. Construction of the iterations The problem

.9%, has the form: given t E H-“*(r),

find u E H’(0,

A) such that

‘-Au=O, OJJL.

Ul

yo

=

0

,

a;r(Ay,,u + ylu) = u;(t),

A E

R’ .

THEOREM 5.3. The following two statements are equivalent: (i) u is solution of PL; (ii) u E H:,(n) and satisfies, for all v E H;,(a), the equality

I

R

Vu * Vv da + A(you, yov)j- = (t, yov)r .

The problem ?& has a unique solution which depends continuously on t. The proof follows quite simply from the Lax-Milgram lemma. The operator which associates with u E H:,(o), the solution of CP=,an element denoted by @. The following diagram indicates the pathway from one iteration

H’ Wn)+ 70 I H l’*(r) i

t E

H-l’*(r) is

to the next.

H’(0) @ I H-“‘(r) j* I

Recall from Section

2.1 that i is the injection

from H”2(r) into L2(r). It is obvious

that A is

56

P. Witomski, Heat equation with nonlinear boundary conditions

a continuous

operator from Hl(f2) into H’(n)

u1 2

uz a.e. j

and that A is increasing in the sense that

Au, 2 Au2 a.e. on a.

The sequence of the iterations is given by u, = Au,,_,. DEFINITION

5.4. We call u a lower (or upper) solution if Au 2 u (or Au 5 u).

PROPOSITION

5.5. (i) u. = u, - Tf is a lower solution of L?$+. (ii) v0 = 0 is ix2 upper solution of pL.

PROOF.

Put u1 = Au,. We have -A(u,-uo)=O, alyo(u, -

a;(hyo(ul-

uo)= -(UC- T,)>O, uo)+ Yl(Ul - uo))=

a-

cp(YoUo)) 2 0

3

since cp(u, - Tf) zz 0. It follows from Lemma 5.1 that u1 = Au0 2 uo. (ii) The proof is identical with ~(0) 2: 0. .5.3. The convergence of the iterations THEOREM 5.6. The sequence of iterations u, = Au,_,, v” = AC’ to u_(respectively a). The solution of 3% is given by g = V = u. PROOF.

converges strongly in H’(0)

We first show that u, converges. We have a sequence of measurable

functions such

that

Set u = sup, u,. The theorem of dominated convergence applied to (u - u,)* shows that un 3 e, strongly in L’(R). Suppose now that the sequence u, is bounded in H’(0). {Au,, n E N) is relatively compact in H’(Q). (A transforms the bounded set into a relatively compact set since j* is compact and cp is Lipschitzian.) Thus, there exists a sub-sequence u,, which converges strongly to u in H’(0) and we have u_= u. It remains to be shown that u, is bounded in H’(a). Use the variational form. For u, E Iif,

I

R

and

V v E H&(f2),

Vu, *Vu df2 + A(you,, ~ov)r = (a-l, you)r

where t,_, = - (P(~~~u,,_~)+ AY”u~_~,Hence on taking ZI= u,:

57

f? Witomski, Heat equation with nonlinear boundary conditions

Since you0 5 ~oyoul d - - - 5 you, I - * - I 0 the sequence ‘you, is bounded in 1;‘(r). It follows that ~~VU,,~(~~~~~ is bounded and, using Friedrich’s inequality, that u, is bounded in H’(0). Since A is continuous on H’(a) at the limit we find that _u= A_u and in the same way ~7= AC. _U- U is the solution of -A(t+j)=o, i

wo(g - 6) = 0,

1a4u u - ii therefore

- fi) = a:(-&w)

+ qY(yo6)).

satisfies

0 = IIW The monotony of Pp,.

of cp and Friedrichs’ inequality

imply that _u= 6 = u is the unique

solution

This method of upper and lower solutions yields an algorithm for the numerical calculation of an approximation of 9 M. This is developed in greater detail in the next section. ~~~~~ 5.7. The above theoretical results can be extended to a non-linear Lipschitzian operator cp from L*(T) into L*(r) by using the fixed point theory with compactness methods (Schauder, Schaeffer, as quoted in [16]).

6. The numerical method 6.1. Rotational invariance of the problem PM

Let (Ox, Oy, Oz) be ortho~ormai axes in R3. The domain a, the boundaries r1 and r2, have symmetry of revolution about the Oz-axis. We will show that the solution u of CPMalso has this symmetry.

NOTATION

6.3. Let 9&, be the rotational

sf3 : MC6 YI +-+M(&,

Y&

operator (02, 8):

Ze),

where x,=xcosO-ysinf?,

y,=xsinBfycos8.

9?!@ leaves each of $2, r1 and r2 invariant.

P. Witomski, Heat equation with nonlinear boundary conditions

58

Be is linear, continuous

and one to one; its inverse is

B2,l= SR, . We consider the mapping Ye from 9(a)’

into S@)

defined by

Ye is an isometry from %‘(jz) to %?‘(a) in the norm of H’(a); a linear isometry from H’(0) to H’(a) (see f16]).

it extends by density to become

6.2. u E H”(a) is said to be invariant under rotation if \d 0 E I%Y*u = u.

DEFINITION

A density argument and the rotational invariance properties can be used to show that the solution of CPMand 9% are invariant under rotation about the &-axis (see [16]). DEFINITION 6.3. The sub-space of Hk,(fl) tions will be denoted V, namely,

PROPOSITION

I

R

represented

by the rotationally

invariant func-

6.4. (i) u is a solution of 62% if dnd only if u E V and satisfies

Vu v Vu dS1 f (&ou),

you),- = 0,

\d z, E V.

(ii) u is a so~~t~o~of FL if and okay if EIE V and sutis~es

for all v in V.

The proof is obvious. It then becomes quite natural to introduce the following. 6.2. approximation

to the problems 9% by finite ele~~~entsof revolution

Given the symmetry of revolution of the probiem we seek the numerical solution of 2% in a finite dimensional sub-space of rotationally invariant functions. The sag-space V, of V Denote the following bouuued open set of R2 by

6.21.

w=nn{(x,y,2)ER3,X>0,y=0}.

We will suppose that w is an open polygonal set of IX’. We take a triangulation 19(d)

is the space of infinitely

differentiable

functions

on 6.

F,,, on W

P. Witomski, Heat equation with nonlinear boundary conditions

59

defined by: (i) V TEF,,, TCG and UTET,,T=h; reduces to one (ii) if T, and T2 belong to Jr,,r then either T, II T2 = 0 or their intersection side or to a vertex; (iii) if we set flh = r, fl6 and r~, = F2 n W, then r I,, and r2,, are junctions of sides of triangles. it is the largest diameter of the triangles The quantity h parameterises the triangulation: constructed. We now associate the following sub-space with Yh: W, = {u E %“((;j), V T E Y,,, uIT affine} , “Irh = {u E Wh, Uj~lh = 0).

We denote

by {bi} a basis of Y,,. We now construct

Vk using the basis of Y,, by the mapping

9 : Lip(O) + Lip(B) defined by .5&(x, y, z) = @(VT, z). Lip(O) (similarly fi) denotes the space of .Lipschitzian functions on (I, (similarly on 0). It is quite clear that the functions Ybi produce a space V, whose elements are invariant with respect to rotations about the Oz-axis. We also have dim V,, = dim Y,,. The problem PM is expressed in IX2 by introducing (see [ll])

L?(w)= { u H:(W)=

measurable,

{U E L:(W);

Iw

War,

u2r dr dz < +m

,

au/a.2 E L?(W)}.

With the scalar product

(u, u) H:(w) H:(w)

is a Hilbert

=

space.

PROPOSITION 6.5. The mapping Y extends into a linear continuous mapping from H:(o) H’(a) and gives

PROOF

to

(see [16]).

From a practical point of view the problem reduces to a classical solution in a domain of IX’. The integrals are taken with respect to the measure r dr dz. The scalar product Jn Vu . Vv d0 defines a norm on V which is equivalent to the norm induced by H’(a). We will use the notation

60

l? Witomski, Heat equation with nonlinear boundary conditions

rh: the projection of V on V,, in the sense of this scalar product; and ph: the injection of vh into V. In general we will omit ph in order to simplify the notation. 6.2.2. The problem PM(h) The numerical solution of the problem PM by finite elements method consists in:

v, ,

‘find u,, E v,, such that, for ah t.$,E vuh PM(h)’

9,(h)

’ vvh

dfi

+

(&youh),

yovh)T

=

0.

I R

is solved by linearisation

as explained in Section 5.2.

6.3. The algorithm for solving PM(h) The linearisation

of PM(h) leads to its reformulation

given w,, E v,,

find uh E v,,

as

such that, for all v,, E v,, ,

PI-(h) vuh

’ VVh

da

+

‘+(youh,

yovh)T

=

(--(p(yowh)

+

Ayowh,

yovh)T.

BL(h) has a unique solution Let A,, be the

the

projection

I

R

of

Vu - Vv da + A(you, y&r

PROPOSITION PROOF.

uh. operator which associates the Solution uh of 2%(h) with wh. We write & for V on vh in the sense of the modified scalar product:

(A > 0).

6.6. A/, = q,,Ap,,.

wh E vh. Aphvh satisfies

Let

YoVb= I V(APhWh)*VV da + A(Yo(AphWh),

(-&owh)

+

AYoWh,

3/o&-,

v

?J E

v.

0

Writing the orthogonality

I

v(q&‘hwh)

* vuh

property for the projections

d0

+

A(yo(q&Wh),

YOvh)f

R =

(-(P(yOwh)+

so, qh&hWh is a solution conclude that A,, = q&h.

AyoWh,

YOvhb,

of 9,(h).

v

vh

E

vh

The uniqueness

we obtain

=

. implies u,, =

&,Wh

=

q&hWh

and we

P. Witomski,

Heat equation with nonlinear boundary conditions

63.1. Numerical lower and upper solutions We Will take uh t0 be a nUmeriCal lower

(Ahuhs

A,&,, 2 u,,

The positiveness

uh)

(Upper)

SOhtiOn

61

if

.

is taken in the classical sense since vh C %“(fi).

6.3.2. Numerical algorithm The algorithm described uses the techniques of lower and upper solutions as in Section 5. We obtain an enclosure of the solution. The proof of the convergence shows the existence and the uniqueness of a solution of PM(h). The algorithm is as follows: (1) Determination of a lower solution UOhand of an upper solution II:. (2) Construction of the sequence &h = AhU+I)h and v: = &,vi-‘. We will show the convergence under the following hypothesis: Let @h be the operator from L’(f) to V, which associates with g the solution @h(g) of the equation (VV~E vh): (NN).

I

v@h(g)

- vvh

da+

A(?‘o@h(g),

yOvh)l-

=

(g,

yOvh)I-.

.

R

@h is supposed to be positive in the following sense

This hypothesis is related to the discrete maximum principle and the notion of non-negative triangulation (see [6]). We will develop this point later. 6.4. Convergence of the algorithm The proof of the convergence of the algorithm is based on the increasing nature of Ah which we show in the following lemma. LEMMA

6.7. Ah is an increasing operator from Vh to Vh, in the sense that ?.I;2 v; =$ A,vj, 2

&,?I;.

PROOF.

The choice of A allows us to write:

-cP(YoVi) + AYou;+ (P(YoV2h)(NN) 3

Ay04

2 0

a.e. on

@h(-cp(yod)+ Ayovt,+ (P(YOZ$ - hod) 2 0,

which together with the linearity of @h is equivalent to fi,V2,

-

r.

A,& 20.

62

l? Witomski, Heat equation with nonlinear boundary conditions

THEOREM 6.8. Under the hypothesis (NN) we have: (i) ux = 0 is a numerical upper solution, voh = u, - Tf is a ~~rner~~~~lower solution. (ii) The iterations u: = A,,uz-’ and v&, = &v(“-Qh converge towards the unique solution of g&h ).

PROOF.

(i) We recall the hypothesis made in Section 2.3 concerning Olcp(O)

and

4p.We have

cp(u,- TJi:O.

Further -q(O) i: 0 and (NN) above imply --Ah(O) 2 0. uz is therefore solution. Moreover it is easily verified that @&-$+&-

PD(k-

T,) 5 0

T,))=Ah(u,-

T,)-

a numerical

upper

@a-T,).

and (NN) imply further that

A,,(& - T,) 2 Z&- Tf . voh = ua- Tt is therefore a numerical lower solution. (ii) The sequences ui and v,h thereby generated nature of Ah we obtain

Since the sequences are bounded and monotone limu;t= ”

uh

and

lim&h=t& n

are monotone

and from the increasing

we have

when n--++m,

The uniqueness of the solution of PM(h) is obtained coercivity, and we have UJ,= ah.

as in the continuous

case using the

6.5. Sufficient conditions on the triangulation yh

We give the conditions on J5h which ensure that the hypothesis (NN) is satisfied. They should be compared to the conditions given in [6] in the case of the uniform convergence of the finite element method for a model Dirichlet problem. We will say that a triangulation J“h is non-negative if the following two hypotheses are satisfied: (Y),: For all h > 0 and all triangles T E yh, the angles of T are strictly less than $rr. (Y),: For all h > 0 and all triangles T having a side common with rzh, the angle opposite this side is less than 4~ - E (E > 0 independently of h). The following proposition ensures that the hypothesis (NN) is satisfied. PROPOSITION

6.9. Let G be the matrix of the linear system associated with 63$_(h).Let

63

P. Witomski, Heat equation with nonlinear boundary conditions

x = [Xl, . . . , XN(h)lE RN(hJbe such

that GX z 0. Under the hypothesis (Q

and (9x)2 and for h

sujjiciently small, then, X 2 0. The proof is long and rather technical; the details can be found in [16]. 6.6. Estimation of the error We now estimate the error between the theoretical solution of 9&h). (i) The theoretical solution u E V of PM satisfies Vu . Vu da + (~(you), y&r

= 0,

solution

of PM and the numerical

V v E V.

(ii) The numerical solution uh E V, of 9%(h) satisfies

I vu/,-V&I dfi + (&%%I),

yLlvh)l-

= 0,

v

vh

E vh .

R

(iii) We denote by u(h) the solution in V of

I Vu(h) * Vv do+

A(y&h),

yov>r = (-cp(yorhu) + Ayorhu,yovh-, V v E V.

R

PROPOSITION

PROOF.

6.10. There exists a constant K independent of h such that

We first show that a constant K1 exists such that

lb - u(WllN’(R) 5 K~llu-

fhuIIH1@)

.

We have

I Mu - u(h))l*do= R

(dwhu) - dyou), ydu - u(h))b + A(yo(u(h) - u), ydu - u(Wr + A(Yo(U - rhu), Yo(U - W)))r .

Since cp is Lipschitzian, with constant lower than A, we obtain

I

R

IV(u - u(h))12dfl~

~A//Yo(u -

Using the Friedrichs’ inequality we obtain

u(h))(l - (lyO(u- r&)1/ .

64

P. Witomski, Heat equation with nonlinear boundary conditions

(where c is a constant), whence

lb - Wll

H’(R)

s

Kll(u

-

d(H’(O)

.

Similarly we can show that a constant K2 also independent

where r,+(h) is the projection

of h exists such that

of u(h) on V,. The proof is therefore complete.

The evaluation of the error involved in replacing u by &, is therefore reduced to estimating the distance between u and vh. 6.7. A theorem of convergence DEFINITION 6.11. The triangulation have a lower bound $ > 0 (see [4]). THEOREM F?

yh is said to be regular if: V h > 0, all the angles of Fh

6.12. If the triangulation yh is regular then we have

IIu - rhuI(H’(o) = 0, tl u E v.

The usual results cannot be applied directly since vh has no polynomial basis. The problem is reduced in two dimensions. LEMMA

6.13. (see [16]). 9(a)

The proof of Theorem LEMMA that

rl V is dense in V for the norm H’(a).

6.12 also requires

6.14. If the triangulation yh is regular, then V u E g(a)

fy

Ill4-

UhllHl(fZ)

=

0

n V there exists vh + Vh such

.

PROOF. Let u E 9&?) f~ V; denotmg its trace on 0 by Tu, we have TuIrIk = 0. Since the triangulation is regular the projection ~~yh(Tu) on ?‘“,,converges strongly to Tu in H’(a) (see [I!!]).According to Proposition 6.5 we have

lb - %‘~TU)IIH~W, = /Iflu - +d-u)llH’(W = v%llTu -p~rh(Tu)lI~~~).

65

l? Witomski, Heat equation with nonlinear boundary conditions

Since 0 is bounded: towards u in H’(0).

I(- llH;(u, 5 KIJ . IIHltnj, .Yp,,(Tu) E V, and therefore

converges

strongly

The proof of Theorem 6.12 now follows directly. Given u E V, L?, is chosen in 9(fi) n V such that (1~- &(lnqnj < E. According to Lemma 6.14, there exists vh such that limh,ol)u - v,,I(~~~)= 0. Therefore limh-ollu - T,J&P(~) s E, V E > 0. The theorem is proved. 6.8. Numerical

experiments

We now give the results of numerical calculations derived from the experimental values obtained for the growth of a germanium crystal. The results after three iterations are shown in Fig. 3. We have used the values: initial lower solution

initial upper solution

T = Tf = 937°C T = Tf - u. = 573°C T = T, = 947°C T = Tf = 937°C

I

in the melt, in the crystal,

in the melt, in the crystal.

The numerical experiments were carried out with the finite element code DELTA installed at the Institute of Applied Mathematics at Grenoble (1.M.A.G).

Fig. 3. Triangulation of the domain. Isothermals in the melt with 1°C steps. Isothermals upper solutions; - ----- lower solutions. steps. -

[12]

in the crystal with u)“C

66

P. Witomski,

Heat equation with nonlinear boundary

conditions

References [I] [2] [3] [4]

J.P. Aubin, Approximation of Elliptic Boundary Value Problem (Wiley, New York, 1972). J.C. Brice, The Growth of Crystals from the Melt (North-Holland, Amsterdam, 1968). F.E. Browder, Probltmes Non Lineaires (Presses de I’Universite de Montreal, Montreal). P.G. Ciariet, Numerical analysis of the finite eIement method, Seminaire de Mathematiques superieure de Montreal, Universite de Montreal, 1975. [5] P.G. Ciarlet and P.A. Raviart, General Lagrange and Hermite interpolation in R” with applications to the finite element method, Arch. Rational and Mech. Anal. 46 (1972) 177-199. [6] P.G. Ciarlet and P.A. Raviart, Maximum principle and uniform convergence for the finite element method, Comput. Meths. Appl. Mech. Engrg. 2 (1973) 17-31. [7] G. Duvaut and J.L. Lions, Les Inequations en Mecanique et en Physique (Dunod, Paris, 1972). [8] R. Glowinski, J.L. Lions and R. Tremolikes, Analyse Numeriques des Inequations Variationnelles, Tome 2 (Dunod, Paris, 1976). [9] R.A. Laudise, The Growth of Single Crystals (Prentice-Hall, Englewood Cliffs, NJ, 1970). [In] J.L. Lions and E. Magenes, Problitmes aux limites Non Homogenes et Applications (Dunod, Paris, 1968). [ 111 J. Necas, Les Methodes Directes dans la Theorie des Equations aux Derivees Partiehes (Masson, Paris, 1967). 1121 A. Poncet, Autour de I’ecriture d’un code d’elements finis, Thesis, Universite Scientifique et Medicale de Grenoble, 1979. [ 131 I. Stakgold, Green’s Functions and Boundary Value Problem (Wiley, New York, 1979). [ 141 G. St~pacchia, Equations ~lliptiques du Second Ordre a Coefficients Discontinus {Presses Unive~itaires de Montreal, Montreal). [ 151 A. Weil, Elements des Cchanges thermiques (Gauthier-Villars, Paris, 1965). [I61 P. Witomski, Modelisation et etude numerique d’une probleme de croissance cristalhne, Thesis, Universitt Scientifique et Medicale de Grenoble, 1977. [ 171 G. Strang and G.J. Fix, An analysis of the Finite EIement Method (Prentice-ball, Englewood Cliffs, NJ, 1973).