APPROXIMATION OF A SOLIDIFICATION PROBLEM

(1)

APPROXIMATION OF A SOLIDIFICATION PROBLEM

Rajae ABOULA¨ ICH

^∗

, Ilham HAGGOUCH

^∗

Ali SOUISSI

^∗∗

A two-dimensional Stefan problem is usually introduced as a model of solidification, melting or sublimation phenomena. The two-phase Stefan problem has been studied as a direct problem, where the free boundary separating the two regions is eliminated using a variational inequality (Baiocchi, 1977; Baiocchi et al., 1973; Rodrigues, 1980; Saguez, 1980; Srunk and Friedman, 1994), the en- thalpy function (Ciavaldini, 1972; Lions, 1969; Nochetto et al., 1991; Saguez, 1980), or a control problem (El Bagdouri, 1987; Peneau, 1995; Saguez, 1980). In the present work, we provide a new formulation leading to a shape optimization problem. For a semidiscretization in time, we consider an Euler scheme. Under some restrictions related to stability conditions, we prove an L

²

-rate of convergence of order 1 for the temperature. In the last part, we study the existence of an optimal shape, compute the shape gradient, and suggest a numerical algo- rithm to approximate the free boundary. The numerical results obtained show that this method is more efficient compared with the others.

Keywords: Stefan problem, free boundary, shape optimization, Euler method, finite-element method

1. Problem Statement

We consider the open bounded domain Ω ⊂

²

defined by Ω = x = (x

1

, x

2

) ∈

²

| 0 < x

1

< 1, 0 < x

2

< 1

(1) The boundary of Ω is written as ∂Ω. Time is denoted by t ∈ ]0, T [ , 0 < T < ∞. The field of temperature is θ : Ω × ]0, T [ −→

²

, so θ (x, t) is the temperature at point x ∈ Ω at time t ∈ ]0, T [. Let us denote by θ

^c

the fusion/solidification temperature.

At time t ∈ ]0, T [, the open bounded domain Ω ⊂ is partitioned as follows:

Ω = Ω

L

(t) ∪ Ω

S

(t) ∪ S (t) , (2)

∗ Universit´e Mohammed V-Agdal, Ecole Mohammadia d’Ing´enieurs, L.E.R.M.A, B.P. 765, Rabat, Morocco, e-mail: aboul,haggouch @emi.ac.ma

∗∗ Université Mohammed V-Agdal, Faculté des Sciences, Département de Mathématiques et d’Informatique, B.P. 1014, Rabat, Morocco, e-mail: souissi@fsr.ac.ma@emi.ac.ma

(2)

where

Ω

L

(t) = {x ∈ Ω | θ (x, t) > θ

^c

} , (3) Ω

S

(t) = {x ∈ Ω | θ (x, t) < θ

c

} , (4)

S (t) = {x ∈ Ω | θ (x, t) = θ

c

} . (5)

The interface S (t) separating the solid and liquid phases is a free boundary and is an unknown of the problem. The domain Ω is shown in Fig. 1.

x

2

x

1

1 1 0

S(t)

Ω

L

(t)

Ω

S

(t)

Fig. 1. The geometry of the domain.

We introduce the following notation:

Q = Ω × ]0, T [ , (6)

Q

L

= [

t∈]0,T [

Ω

L

t × {t} , (7)

Q

S

= [

t∈]0,T [

Ω

S

(t) × {t} , (8)

Σ = [

t∈]0,T [

S (t) × {t} . (9)

Let functions θ

0

∈ H

¹

(Ω) and θ

∂Ω

∈ L

²

(∂Ω) be given such that the following compatibility condition is satisfied:

θ

0

(x) = θ

∂Ω

(x) , a.e. x ∈ ∂Ω.

(3)

Thus the unknowns (θ, Q

L

) are the solutions of the following evolution free-boundary problem:

(P

1

)



 

 

 

 

Find θ and Q

L

such that

∂θ

∂t − div (C (θ) ∇θ) = 0 in Q

L

∪ Q

S

, C (θ) =

( c

s

if θ < θ

c

, c

l

if θ > θ

c

, θ (x, 0) = θ

0

(x) , x ∈ Ω,

θ (x, t) = θ

∂Ω

(x) , x ∈ ∂Ω, t ∈ ]0, T [ , θ (x, t) < θ

c

, (x, t) ∈ Q

^S

, θ (x, t) > θ

c

, (x, t) ∈ Q

^L

, θ (x, t) = θ

c

, (x, t) ∈ Σ, h C (θ) ∇θ·

^→

n i

S(t)

= λ V ·

^→

n, x ∈ S (t) , t ∈ ]0, T [ . Here λ is the latent heat of the material (a strictly positive coefficient), V signifies the velocity of the free boundary, c

s

means the diffusivity of the solid part, c

l

stands for the diffusivity of the liquid part (c

s

and c

l

are strictly positive coefficients), and n is the unitary normal to S (t) pointing towards Ω

L

(t).

The major difficulty in a direct problem lies in the fact that the moving boundary is utilized explicitly in the equation for the thermal state of the system. This difficulty is circumvented in Section 2, using the characteristic function of the liquid region. Such a formulation transforms the initial problem into a partial differential equation valid on the whole cavity occupied by the material. In Section 3 we re- call the regularization method proposed in (Humeau and Souza del Cursi, 1993). In Section 4 we discretize the regularization problem and study the convergence of the proposed scheme. Section 5 estabilishes an approximation of the free boundary for the obtained stationary problem using a shape optimization method. The existence of the free boundary, details of the shape gradient computation, and numerical results are provided.

2. Problem Reformulation

Let us introduce the following spaces:

H = L

²

(Ω) , V = H

¹

(Ω) , V

0

= H

₀¹

(Ω) , (10) with their usual scalar products

H = L

²

(0, T ; H) , B = L

²

(0, T ; V ) , B

⁰

= L

²

(0, T ; V

0

) , (11)

(4)

respectively. Let (·, ·) denote the scalar product on H corresponding to the norm k·k.

Consider the real-valued function

χ

QL

(x, t) =

( 1 if (x, t) ∈ Q

^L

,

0 otherwise. (12)

The problem equivalent to (P

1

) can be written as follows:

(P

2

)



 



 



Find θ ∈ B and Q

^L

⊂ Q such that

∂θ

∂t − ∇ · (C (θ) ∇θ) = −λ ∂

∂t χ

QL

in Q,

θ (x, t) = θ

∂Ω

(x) , x ∈ ∂Ω, t ∈ ]0, T [ , θ (x, 0) = θ

0

(x) , x ∈ Ω,

θ (x, t) > θ

c

, (x, t) ∈ Q

L

, θ (x, t) < θ

c

, (x, t) ∈ Q

S

, θ (x, t) = θ

c

, (x, t) ∈ Σ.

3. Problem Regularization

In order to overcome some numerical difficulties due to the discontinuity of the function C (θ), we consider the regularization method proposed in (Humeau and Souza del Cursi, 1993). Introduce φ : −→ given by

φ (β) = 3β − 2β

²

. (13)

Let ε > 0 be a fixed parameter. We set

C

ε

(β) =



 



 



C (β) if β / ∈ (θ

^c

− ε, θ

^c

+ ε) ,

c

s

φ ε − β + θ

^C

2ε

+ c

l

φ β + ε − θ

^C

2ε

otherwise.

Hence C

ε

can be considered as a Lipschitz continuous approximation of C.

(5)

β C

L

C

S

θ

c

− ε θ

c

θ

c

+ ε 0

C

ε

(β)

Fig. 2. Regularization of C.

The regularized problem associated with (P

2

) can be formulated as follows:

(P

3

)



 



 



Find θ

ε

∈ B and Q

L

⊂ Q such that

∂θ

ε

∂t − ∇ · (C

^ε

(θ

ε

) ∇θ) = −λ ∂

∂t χ

QL

in Q,

θ

ε

(x, t) = θ

∂Ω

(x) , x ∈ ∂Ω, t ∈ ]0, T [ , θ

ε

(x, 0) = θ

0

(x) , x ∈ Ω,

θ

ε

(x, t) > θ

c

, (x, t) ∈ Q

^L

, θ

ε

(x, t) < θ

c

, (x, t) ∈ Q

S

, θ

ε

(x, t) = θ

c

, (x, t) ∈ Σ.

For notational convenience, in what follows the index ε will be omitted, so θ

ε

and C

ε

will be denoted respectively by θ and C. Consider a function g ∈ H

¹

(Ω) such that

g (x) =

( θ

0

(x) if x ∈ Ω,

θ

∂Ω

(x) if x ∈ ∂Ω. (14)

Introduce the change of variables

u (x, t) = θ (x, t) − g (x) in Q. (15)

(6)

Then Problem (P

3

) can be written as follows:

(P

4

)



 



 



Find u ∈ B

⁰

and Q

L

⊂ Q such that

∂u

∂t − ∇ · (C (u + g) ∇ (u + g)) = −λ ∂

∂t χ

QL

in Q, u (x, 0) = 0, x ∈ Ω,

u + g > θ

c

, (x, t) ∈ Q

^L

, u + g < θ

c

, (x, t) ∈ Q

S

, u + g = θ

c

, (x, t) ∈ Σ.

Consider now the problem

(P

5

)



 



 



Find u ∈ B

0

such that

∂u

∂t − ∇ · C (u + g) ∇ (u + g) = −λ ∂

∂t χ

QL

in Q, u (x, 0) = 0, x ∈ Ω.

Let V be a closed subspace of B

0

defined by V =

v ∈ B

⁰

∂v

∂t ∈ H and v (x, 0) = v (x, T ) = 0, x ∈ Ω

(16) and equipped with the scalar product

(u, v)

_V

= (u, v)

_B₀

+ ∂u

∂t , ∂v

∂t

H

. (17)

The variational formulation associated with (P

5

) is as follows:

(P V

5

)



 



 



Find u ∈ B

⁰

such that

∂

∂t u, v

+ C (u + g) ∇ (u + g) , ∇v = −λ ∂

∂t χ

QL

, v

, ∀v ∈ V, u (x, 0) = 0, x ∈ Ω.

The existence and uniqueness of the solution to (P

5

) are established in (Haggouch, 1997; Humeau and Souza del Cursi, 1993), using the elliptic regularization method, cf. (Lions, 1969).

In the next section, we shall discretize Problem (P

5

) in time and then study the convergence of the proposed scheme.

4. Time Discretization

Consider a strictly positive integer N > 0 which implies the discretization step τ = T /N , and denote by t

n

the grid points of [0, T ] : t

n

= nτ, 0 ≤ n ≤ N. We define

Ω

ⁿ_L

= x ∈ Ω | θ (x, t

ⁿ

) > θ

c

(7)

and

χ

Ωⁿ_L

(x) = χ

ΩL

(x, t

n

) =

( 1 if θ (x, t

n

) > θ

c

, 0 otherwise.

Let

u

n

(x) ' u (x, t

n

) , u

0n

(x) = u

0

(x, t

n

) = θ

0

(x, t

n

) − g (x) . (18) Then the discretization of Problem (P

5

) can be written down as follows:

(P

n

)



 

 

Find (u

n+1

)

_{0≤n≤N −1}

⊂ V

0^N

such that u

n+1

− u

n

τ + ∇ · C (u

ⁿ

+ g) ∇ (u

ⁿ⁺¹

+ g) = −λ χ

_Ωⁿ⁺¹

L

− χ

^Ωⁿ_L

τ in Ω.

Note that we have a linear problem that changes with each value of n. That means that solution to (P

n

) requires N steps.

The variational problem associated with (P

n

) is as follows:

(P V

n

)



 

 

 

 

Find (u

n+1

)

_{0≤n≤N −1}

⊂ V

0^N

such that

u

n+1

− u

n

τ , v

+ C(u

n

+ g)∇(u

ⁿ⁺¹

+ g), ∇v

= −λ

χ

_Ωⁿ⁺¹

L

− χ

Ωⁿ_L

τ , v

, ∀v ∈ V

⁰

. Proposition 1. The function u

n

being a solution to (P V

n

) satisfies the following discrete a-priori estimates:

0≤n≤N

max ku

ⁿ

k ≤ C

¹

,

N −1

X

n=0

ku

n+1

− u

n

k

²

≤ C

2

,

τ

N

X

n=1

k∇u

ⁿ

k

²

≤ C

³

,

where C

1

, C

2

and C

3

are constants independent of τ . Proof. Setting k

1

= min(c

s

, c

l

) and k

2

= max(c

s

, c

l

), we have

k

1

≤ C (σ) ≤ k

²

. (19)

Choosing v = u

n+1

in (P V

n

) and applying the following elementary equality:

2p (p − q) = p

²

− q

²

+ (p − q)

²

, ∀ (p, q) ∈

²

, (20) we see that

ku

n+1

k

²

− ku

n

k

²

+ ku

n+1

− u

n

k

²

+ 2τ k

1

k∇u

n+1

k

²

+ 2τ k

1

∇ (g) , ∇u

ⁿ⁺¹

+ 2λ

χ

_Ωⁿ⁺¹

L

− χ

^Ωⁿ_L

, u

n+1

≤ 0. (21)

(8)

Moreover, applying the Young inequality to the last terms of (21) yields 2τ k

1

∇(g), ∇u

ⁿ⁺¹

≤ 2τk

¹

k∇gk k∇u

ⁿ⁺¹

k

≤ k

₁²

k∇gk

²

ε + ετ

²

k∇u

n+1

k

²

, ∀ ε > 0 (22) and

2λ χ

_Ωⁿ⁺¹

L

− χ

Ωⁿ_L

, u

n+1

≤ 2λ

χ

Ωⁿ⁺¹_L

− χ

Ωⁿ_L

ku

n+1

k

≤ 4λ √

mes Ω ku

ⁿ⁺¹

k

≤ 4λ

²

mes Ω

β + β ku

n+1

k

²

, ∀ β > 0. (23) Therefore, choosing arbitrary β and ε such that (1 − β) > 0 and (2k

1

− ετ) > 0, we get

(1 − β) ku

ⁿ⁺¹

k

²

− ku

ⁿ

k

²

+ ku

ⁿ⁺¹

− u

ⁿ

k

²

+ τ (2k

1

− ετ) k∇u

ⁿ⁺¹

k

²

≤ k

²₁

k∇gk

²

ε + 4λ

²

mes Ω β ≤ c.

Summing this inequality over n, 0 ≤ n ≤ p − 1, 1 ≤ p ≤ N, and using the fact that u

0

= 0, we get

(1 − β) ku

p

k

²

+

p−1

X

n=0

ku

n+1

− u

n

k

²

+ (2k

1

− ετ)

p

X

n=1

τ k∇u

n

k

²

≤ β

p−1

X

n=1

ku

n

k

²

+ pc.

By the discrete Gronwall inequality (Raviart and Girault, 1981), we obtain

(1 − β) ku

^p

k

²

+

p−1

X

n=0

ku

ⁿ⁺¹

− u

ⁿ

k

²

+ (2k

1

− ετ)

p

X

n=1

τ k∇u

ⁿ

k

²

≤ pce

^{β p}

.

Hence there exist constants C

1

, C

2

and C

3

independent of τ such that

0≤n≤N

max ku

n

k ≤ C

1

,

N −1

X

n=0

ku

n+1

− u

n

k

²

≤ C

2

,

N

X

n=1

τ k∇u

n

k

²

≤ C

3

,

which completes the proof.

Using these estimations and non-linear analysis, we can show the following the- orems (Humeau and Souza del Cursi, 1993):

Theorem 1. Problem (P

n

) has a unique solution u

n

in V

0

.

(9)

Consider t

_i+¹

2

= t

i

+ τ

2 , t

_i−¹

2

= t

i

− τ

2 , I

i

= h t

_i−¹

2

, t

_i+¹

2

i ∩ ]0, T [ ,

I

i

=

( 1 if t ∈ I

ⁱ

, 0 otherwise and

u

^N

(x, t) =

N

X

n=0

u

n

(x) I

n

(t) .

Theorem 2. The sequence u

^N

N >0

converges weakly in B

⁰

to a solution u to Problem (P

5

) as N → ∞.

5. Error Analysis

In this paragraph, we use a regularization of χ

ΩL

to obtain a continuous function χ

ΩL

such that ∂χ

ΩL

/∂t and ∂

²

χ

ΩL

/∂t

²

are regular. We consider, for example, the following regularization: For ε > 0, define

χ

^ε_Ω_L

u (x, t) =



 

 

χ

ΩL

u (x, t)

if u (x, t) / ∈ (0, ε) , ψ u (x, t)

ε

otherwise, (24)

where ψ is a function of C

¹

(Q). For notational simplicity, in place of χ

^ε_Ω_L

we will write χ

ΩL

.

In the following, we suppose that

∂u

∂t ∈ L

²

(0, T ; V

0

) , ∂

²

u

∂t

²

∈ L

²

0, T ; V

⁰

, ∂

²

χ

ΩL

∂t

²

∈ L

²

0, T ; V

⁰

. (25) Then we can deduce that

u ∈ C

⁰

(0, T ; V

0

) , ∂u

∂t ∈ C

⁰

(0, T ; H) , χ

ΩL

∈ C

⁰

(0, T ; V

0

) , ∂χ

ΩL

∂t ∈ C

⁰

(0, T ; H) .

(26)

First of all, we show that the consistency error is of order one and that the proposed

scheme is stable.

(10)

5.1. Estimate of the Consistency Error

Let χ

ΩL

(t) and u (t) be two functions defined on Ω by χ

ΩL

(t) : x −→ χ

^ΩL

u (x, t)

and

u (t) : x −→ u (x, t) .

We define the consistency error ε

n

∈ V

⁰

by hε

ⁿ

, υi = 1

τ u (t

n+1

) − u (t

ⁿ

) , v + C u (t

n

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v

+ λ

τ χ

ΩL

(t

n+1

) − χ

^ΩL

(t

n

) , v, (27) where V

⁰

is the dual space of V , and h·, ·i denotes the duality product between V and V

⁰

.

Lemma 1. (Lions, 1969) When n = 2, all the elements ϕ of V

0

satisfy kϕk

L⁴(Ω)

≤ 2

¹⁴

kϕk

¹²

k∇ϕk

¹²

.

Suppose that a constant M > 0 exists such that for all t ∈ [0, T ] we have

(H)

∇ u (t) + g

_L4(Ω)

< M.

Consider

∂

∂t u (t

n+1

) , v + C u (t

n+1

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v

+ λ ∂

∂t χ

ΩL

(t

n+1

) , v = 0. (28) Calculating the difference of (27) and (28), we get

hε

ⁿ

, υi = u (t

n+1

) − u (t

ⁿ

)

τ − ∂u

∂t (t

n+1

) , v

+ λ χ

ΩL

(t

n+1

) − χ

^ΩL

(t

n

)

τ − ∂

∂t χ

ΩL

(t

n+1

) , v

− C u (t

n+1

) + g − C u (t

ⁿ

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v . By Taylor’s formula with integral remainder, defined by

1 τ f (t

n+1

) − f (t

n

) , v

= ∂

∂t f (t

n+1

) , v

+ 1 τ

Z

tn+1

tn

(t − t

n+1

) ∂

²

∂t

²

f (t) , v

dt (29)

(11)

for f = u and f = χ

Ω_L

, we obtain

hε

n

, vi = 1 τ

Z

tn+1

tn

(t − t

n+1

) ∂

²

u

∂t

²

(t) , v

dt

+ λ τ

Z

tn+1

t_n

(t − t

n+1

) ∂

²

∂t

²

χ

ΩL

(t) , v

dt

− C u (t

ⁿ⁺¹

) + g − C u (t

ⁿ

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v . (30)

But

C u (t

n+1

) + g − C u (t

ⁿ

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v

≤

C u (t

n+1

) + g − C u (t

ⁿ

) + g ∇ u (t

ⁿ⁺¹

) + g k∇vk

≤ max

σ∈

C

⁰

(σ)

u (t

n+1

) − u (t

ⁿ

) ∇ u (t

ⁿ⁺¹

) + g k∇vk

≤ max

σ∈

C

⁰

(σ)

u (t

n+1

) − u (t

ⁿ

)

_L4(Ω)

∇ u (t

ⁿ⁺¹

) + g

_L4(Ω)

k∇vk .

From Lemma 1 and the hypothesis (H), we deduce that

C u (t

n+1

) + g − C u (t

ⁿ

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v

≤ c

1

u (t

n+1

) − u (t

n

)

1 2

∇ u (t

n+1

) − u (t

n

)

1 2

k∇vk ,

where

c

1

= 2

¹⁴

M max

σ∈

C

⁰

(σ)

(31)

Since the embedding H

¹

(Ω) → L

²

(Ω) is compact and u (t

n

) ∈ H

0¹

(Ω), we obtain

C u (t

n+1

) + g − C u (t

ⁿ

) + g ∇ u (t

ⁿ⁺¹

) + g , ∇v

≤ c

2

∇ u (t

n+1

) − u (t

n

)

k∇vk . (32) Let

u (t

n+1

) − u (t

n

) = Z

tn+1

tn

∂u (t)

∂t dt. (33)

(12)

Then we have

Z

tn+1

tn

∂u (t)

∂t dt

2

V0

= Z

Ω

∇ Z

tn+1

tn

∂u (t)

∂t dt

²

dx =

Z

Ω

Z

tn+1

tn

∇ ∂u (t)

∂t dt

²

dx

≤ Z

Ω

"

τ Z

tn+1

tn

∇ ∂u (t)

∂t

2

dt

# dx

≤ τ Z

tn+1

tn

"

Z

Ω

∇ ∂u (t)

∂t

2

dx

# dt

≤ τ Z

tn+1

tn

∂u (t)

∂t

2

V0

dt. (34)

By the definition of the norm k · k

V⁰

:

kε

ⁿ

k

V⁰

= sup

v∈V

hε

ⁿ

, vi kvk

V

, (35)

from (30) we get

kε

n

k

V⁰

≤ 1 τ

Z

tn+1

tn

(t − t

n+1

)

∂

²

u

∂t

²

V⁰

dt

+ λ τ

Z

tn+1

tn

(t − t

ⁿ⁺¹

)

∂

²

∂t

²

χ

ΩL

(t)

_V0

dt

+ c

2

τ Z

tn+1

tn

∂u (t)

∂t

2

V0

dt

!

¹2

. (36)

But

1 τ

Z

tn+1

tn

(t − t

n+1

)

∂

²

u (t)

∂t

²

_V0

dt

≤ 1 τ

Z

tn+1

tn

(t − t

ⁿ⁺¹

)

²

1

2

Z

tn+1

tn

∂

²

u (t)

∂t

²

2

V⁰

dt

!

¹2

≤ τ

Z

tn+1

tn

∂

²

u (t)

∂t

²

2

V⁰

dt

!

¹2

(37)

(13)

and

1 τ

Z

tn+1

tn

(t − t

ⁿ⁺¹

)

∂

²

∂t

²

χ

ΩL

(t)

_V0

dt

≤ 1 τ

Z

tn+1

tn

(t − t

n+1

)

²

1

2

Z

tn+1

tn

∂

²

∂t

²

χ

ΩL

(t)

2

V⁰

dt

!

¹2

≤ τ

Z

tn+1

t_n

∂

²

∂t

²

χ

ΩL

(t)

2

V⁰

dt

!

¹2

. (38)

Substituting (37) and (38) into (36), we get

kε

ⁿ

k

V⁰

≤ τ Z

tn+1

tn

∂

²

u (t)

∂t

²

2

V⁰

dt

!

¹2

+ λ τ Z

tn+1

tn

∂

²

∂t

²

χ

ΩL

(t)

2

V⁰

dt

!

¹2

+ c

2

τ Z

tn+1

tn

∂u (t)

∂t

2

V0

dt

!

¹2

. (39)

Then

kε

ⁿ

k

²V⁰

≤ cτ

"

Z

tn+1

tn

(

∂

²

u (t)

∂t

²

2

V⁰

+

∂

²

∂t

²

χ

ΩL

(t)

2

V⁰

+

∂u (t)

∂t

2

V0

) dt

#

. (40)

Summing this inequality over n, 0 ≤ n ≤ p − 1, 1 ≤ p ≤ N, we obtain

τ

p−1

X

n=0

kε

ⁿ

k

²V⁰

≤ cτ

²

"

∂

²

u (t)

∂t

²

2

L²

(

^{0,T ;V}⁰

) +

∂

²

∂t

²

χ

ΩL

(t)

2

L²

(

^{0,T ;V}⁰

) +

∂u (t)

∂t

2

L²(0,T ;V0)

dt

#

. (41)

Using (25), we get the following result:

Proposition 2. Suppose that u satisfies the regularity conditions (25) and that there exists a constant M > 0 such that for all t ∈ [0, T ] we have k∇ (u (t) + g)k

L⁴(Ω)

<

M . Then the consistency error is of order 1, i.e.

τ

p−1

X

n=0

kε

n

k

²V⁰

!

1 2

≤ cτ, 1 ≤ p ≤ N, (42)

where c > 0 is a constant independent of τ .

(14)

5.2. Stability of the Scheme For all n ∈ {0, . . . , N − 1}, we have

1 τ (u

n+1

− u

ⁿ

, v) + C (u

n

+ g) ∇ (u

ⁿ⁺¹

+ g) , ∇v + λ

τ

χ

_Ωⁿ⁺¹

L

− χ

^Ωⁿ_L

, v

= 0. (43)

Set

e

ⁿ

= u (t

n

) − u

ⁿ

. (44)

Proposition 3. If u

0

= u (t

0

), the following stability criterion holds for 1 ≤ p ≤ N:

ke

^p

k

²

+

p−1

X

n=0

e

ⁿ⁺¹

− e

ⁿ

2

+ τ k

1 p

X

n=1

k∇e

ⁿ

k

²

≤ 2τ k

1

p−1

X

n=0

kε

^∗n

k

²V⁰

!

exp(cT ),

where c > 0 is a constant independent of τ and k

1

= min(c

s

, c

l

).

Proof. Subtracting (27) from (43) yields e

ⁿ⁺¹

− e

ⁿ

, v

+ τ k

1

∇e

ⁿ⁺¹

, ∇v

+ τ C u (t

n

) + g − C (u

ⁿ

+ g) ∇ u (t

ⁿ⁺¹

) + g , ∇v

≤ τ hε

ⁿ

, vi , ∀v ∈ V

⁰

, 0 ≤ n ≤ N − 1. (45) Taking v = e

n+1

in (45) and applying the inequality

C (α) − C (β) ≤ max

σ∈

C

⁰

(σ)

|α − β| , (46)

we see that

e

ⁿ⁺¹

2

− ke

ⁿ

k

²

+

e

ⁿ⁺¹

− e

ⁿ

2

+ 2τ k

1

∇e

ⁿ⁺¹

2

≤ 2τ ε

n

, e

ⁿ⁺¹

+ 2τ max

σ∈

C

⁰

(σ)

e

ⁿ

∇ u (t

n+1

) + g , ∇e

ⁿ⁺¹

. (47) The right-hand side is estimated as follows:

2 ε

n

, e

ⁿ⁺¹

≤ 2 kε

ⁿ

k

V⁰

e

ⁿ⁺¹

_V

≤ 2

k

1

kε

ⁿ

k

²V⁰

+ k

1

2 ∇e

ⁿ⁺¹

2

. (48)

(15)

Using Lemma 1 and the assumption (H), we get 2 e

ⁿ

∇ u (t

n+1

) + g , ∇e

ⁿ⁺¹

≤ 2 ke

ⁿ

k

L⁴(Ω)

∇ u (t

n+1

) + g

L⁴(Ω)

∇e

ⁿ⁺¹

≤ 2 M 2

¹⁴

ke

ⁿ

k

¹²

k∇e

ⁿ

k

¹²

∇e

ⁿ⁺¹

≤ c

¹

1 ε ke

ⁿ

k k∇e

ⁿ

k + ε ∇e

ⁿ⁺¹

2

≤ c

¹

1 ε

1 2δ ke

ⁿ

k

²

+ δ

2 k∇e

ⁿ

k

²

+ ε

∇e

ⁿ⁺¹

2

, (49) where ε > 0 and δ > 0 are arbitrary. Hence we have the estimate

2τ max

σ∈

C

⁰

(σ)

e

ⁿ

∇ u (t

n+1

) + g , ∇e

ⁿ⁺¹

≤ k

1

4 τ

k∇e

ⁿ

k

²

+

∇e

ⁿ⁺¹

2

+ τ c

2

ke

ⁿ

k

²

. (50) Substituting (48) and (50) into (47), we get

e

ⁿ⁺¹

2

− (1 + τc

2

) ke

ⁿ

k

²

+

e

ⁿ⁺¹

− e

ⁿ

2

+ 5k

1

4 τ

∇e

ⁿ⁺¹

2

≤ k

1

4 τ k∇e

ⁿ

k

²

+ 2τ

k

1

kε

n

k

²V⁰

. (51)

If we sum this last relation over n, 0 ≤ n ≤ p − 1, 1 ≤ p ≤ N, and use the fact that e

0

= 0, then we obtain

ke

^p

k

²

+

p−1

X

n=0

e

ⁿ⁺¹

− e

ⁿ

2

+ τ k

1 p

X

n=1

k∇e

ⁿ

k

²

≤ τc

²

p−1

X

n=1

ke

ⁿ

k

²

+ 2τ k

1

p−1

X

n=0

kε

ⁿ

k

²V⁰

.

The discrete Gronwall inequality gives

ke

^p

k

²

+

p−1

X

n=0

e

ⁿ⁺¹

− e

ⁿ

2

+ τ k

1 p

X

n=1

k∇e

ⁿ

k

²

≤ 2τ k

1

p−1

X

n=0

kε

n

k

²V⁰

!

exp(c

2

pτ )

≤ 2τ k

1

p−1

X

n=0

kε

n

k

²V⁰

!

exp(c

2

T ). (52)

(16)

Combining Propositions 2 and 3, we derive immediately the following result:

Theorem 3. Under the assumptions of Proposition 2, there exists a constant c > 0 independent of τ such that

0≤n≤N

max ku (t

n

) − u

n

k + τ

p−1

X

n=0

ku (t

n

) − u

n

k

²V

!

1 2

≤ cτ. (53)

6. Formulation in the Framework of Shape Optimization

Under some regularity of the free boundary, we shall expose a new formulation of the semi-discrete problem associated with (P

3

), using the shape optimization techniques.

The existence results for an optimal domain and the shape gradient are presented. For the computation of the gradient, we suggest the material derivative (Zol´esio, 1981) and the duality methods (Lions, 1968).

The partial differential equation in Problem (P

3

) is approximated as θ

n+1

− θ

ⁿ

τ −∇· C (θ

ⁿ

) ∇θ

ⁿ⁺¹

= −λ χ

_Ωⁿ⁺¹

L

− χ

Ωⁿ_L

τ , ∀n ∈ {0, . . . , N − 1} . (54) We introduce functions a

n

and b

n

defined as follows:

a

n

(x) = τ C θ

n

(x) , x ∈ Ω, b

n

(x) = λχ

Ωⁿ_L

+ θ

n

(x) , x ∈ Ω.

Thus we get the following semi-discretized problem associated with (P

3

):

(P

6

)



 

 

 

 

Find (θ

n+1

)

_{0≤n≤N −1}

⊂ V

^N

and Ω

ⁿ⁺¹_L

0≤n≤N −1

⊂ (O

ad

)

^N

such that θ

n+1

> θ

c

in Ω

ⁿ⁺¹_L

,

θ

n+1

< θ

c

in Ω

ⁿ⁺¹_S

, θ

n+1

= θ

c

on S

ⁿ⁺¹

,

and θ

n+1

is the solution of the problem

P Ω

ⁿ⁺¹_L



 



 



Find θ

n+1

∈ V such that

θ

n+1

− ∇ · (a

ⁿ

∇θ

ⁿ⁺¹

) = −λχ

Ωⁿ⁺¹_L

+ b

n

in Ω,

θ

n+1

= θ

∂Ω

on ∂Ω,

where O

^ad

is the set of admissible domains that will be defined later.

For notational convenience, we eliminate the indices n, and the sequences

θ

n

, Ω

ⁿ_S

, Ω

ⁿ_L

, S

ⁿ

, a

n

and b

n

are denoted by θ, Ω

S

, Ω

L

, S, a and b, respectively.

(17)

Introduce the cost functional J (Ω

L

) ≡ J Ω

L

, θ (Ω

L

)

= 1 2 Z

Ω

χ

ΩS

h θ (Ω

L

) − θ

^c

+

i

2

dx + 1 2 Z

Ω

χ

ΩL

h θ

c

− θ (Ω

^L

)

+

i

2

dx. (55) The optimal shape design problem is formulated as follows:

(P

^op

)







ΩL

min

∈Oad

J (Ω

L

) such that θ (Ω

L

) is the solution of P (Ω

L

) . By setting

F =

Ω

L

, θ (Ω

L

) | Ω

L

∈ O

ad

and θ (Ω

L

) is the solution to P (Ω

L

) , (56) the optimization problem can be written as

(P

op

) minimize J (Ω

^L

) | (Ω

L

, θ (Ω

L

)) ∈ F . (57) The new formulation we propose makes use of the regularity of the free boundary.

The existence of the latter was proved in (Baiocchi et al., 1973; Lions, 1969; Saguez, 1980). As in (Lions, 1969), we assume that the free boundary is of measure zero on Ω and, moreover, we suppose that it is defined by a curve described by the equation x

2

= α (x

1

), where α is a regular function. Then there exists a solution in F such that J (Ω

^L

) = 0. We deduce easily the equivalence of Problems (P

6

) and (P

^op

). In the next section, we shall study Problem (P

^op

).

6.1. Existence Result

The existence of an optimal solution to (P

op

) requires the choice of an adequate topology on the admissible domain, permitting to obtain the compactness of O

ad

and the lower semicontinuity of J .

The set of admissible functions which parameterize the free boundary S is defined as follows:

U

^ad

= n

α ∈ C [0, 1]

α (x

1

) − α (x

¹

)

≤ k |x

¹

− x

¹

| ∀ x

^1,

x

1

∈ [0, 1] , α (0) = c

1

, α (1) = c

2

and α (x

1

) < c

3

, x

1

∈ [0, 1] o

. The constants k c

1

, c

2

and c

3

are chosen in such a way that U

^ad

is not empty. U

^ad

is equipped with the following norm:

kαk

∞

= max

0≤x1≤1

α (x

1

)

, α ∈ U

^ad

. (58)

We define α

n

n→∞

α in [0, 1] ⇐⇒ kα

ⁿ

− αk

∞_n→∞

→ 0, (59)

(18)

and the convergence in U

^ad

by

‘α

n _n→∞

→ α’ in U

ad

⇐⇒ α

n

n→∞

α in [0, 1] . (60)

The different regions can be characterized by

Ω

L

(α) = x ∈ Ω | x

²

> α (x

1

) , (61)

Ω

S

(α) = x ∈ Ω | x

²

< α (x

1

) , (62)

S (α) = x ∈ Ω | x

²

= α (x

1

) , (63)

and

χ

ΩL(α)

=

( 1 if x

2

> α (x

1

) , 0 otherwise.

We consider O

ad

as the set of the admissible domains

O

^ad

= Ω (α) ⊂ Ω | α ∈ U

^ad

, (64)

where

Ω (α) = Ω

L

(α) ∪ Ω

^S

(α) ∪ S (α) .

We require O

ad

to be equipped with the appropriate topology and convergence defined by

‘Ω

n

→

n→∞

Ω’ ⇐⇒ ‘α

n

→

n→∞

α’ in U

ad

, (65)

where Ω

n

= Ω (α

n

) and Ω = Ω (α).

Let α ∈ U

^ad

. For any Ω

L

(α) we consider the following boundary-value problem:

P (α)



 



 



Find θ (α) ∈ V such that

θ (α) − ∇ · a∇θ (α) = −λ χ

^ΩL(α)

+ b in Ω,

θ (α) = θ

∂Ω

on ∂Ω.

The cost functional is given by J (α) ≡ J α, θ (α)

= 1 2 Z

Ω

χ

ΩS(α)

h θ (α) − θ

c

+

i

2

dx + 1 2 Z

Ω

χ

ΩL(α)

h θ

c

− θ (α)

+

i

2

dx. (66)

We set

F =

α, θ (α)

α ∈ U

^ad

and θ (α) is the solution to P (α) , (67)

(19)

endowed with the topology defined by the following convergence:

‘ Ω

L

(α

n

) , θ (α

n

)

n→∞

→ Ω

L

(α) , θ (α) ’

⇐⇒







Ω

L

(α

n

) →

n→∞

Ω

L

(α) in O

ad

, θ (α

n

) *

n→∞

θ (α) (weakly) in B. (68) Thus the optimization problem is written as

P

op

(α) minimize J (α)

α, θ (α) ∈ F .

Using the approach presented in (Haslinger and Neittaanm¨ aki, 1988), we have the following results, and the details of the proofs are given in (Haggouch, 1997).

Proposition 4. Let θ

n

= θ (α

n

) be the solutions of P (α

n

) , α

n

∈ U

ad

and Ω

ⁿ_L

= Ω

L

(α

n

). Then there exist a subsequence of {(α

n

, θ

n

)} (again denoted by {(α

n

, θ

n

)}) and elements α ∈ U

ad

, θ ∈ V such that

‘Ω

ⁿ_L_n→∞

→ Ω

L

(α)’ in O

ad

and θ

n

*

n→∞

θ (weakly) in B Moreover, θ solves P (α).

Proposition 5. The function α → J (α) is continuous on U

ad

. Using these propositions, we establish our next theorem.

Theorem 4. There exists at least one solution to Problem P

op

(α) , α ∈ U

ad

. Proof. We define q by

q = inf

α∈Uad

J (α) . (69)

Let Ω

ⁿ_L

= Ω

L

(α

n

) , α

n

∈ U

^ad

be a minimizing sequence, i.e.

n→∞

lim J (α

ⁿ

) = q, (70)

and θ (α

n

) be the solution to Problem P (α

n

).

Proposition 2 implies that there exist a subsequence

α

nj

, θ α

nj

⊂ {(α

ⁿ

, θ (α

n

))} and an element {(α

^∗

, θ (α

^∗

))} ∈ F such that

α

nj

j→∞

α

^∗

in [0, 1] (71)

and

θ α

nj

*

j→∞

θ (α

^∗

) (weakly) in V. (72)

(20)

By Proposition 3, we have

j→∞

lim J α

nj

= J (α

^∗

) . (73)

The uniqueness of the limit implies

α∈U

inf

ad

J (α) = J (α

^∗

) . (74)

6.2. Numerical Approximation of the Free Boundary 6.2.1. Existence of the Gradient

Consider a vector field W, defined on [0, β] × U with values in

²

, U being an open neighborhood of Ω and β > 0. Let W ∈ C [0, β] , D

^k

U,

²

, k ≥ 1. We transform Ω into Ω

τ

through the function T

τ

defined by

T

^τ

(X) = x (τ, X) , (75)

where x (τ, X) is the unique solution to the differential equation

P



 

  d

dτ x (τ, X) = W τ, x (τ, X) , x (τ, 0) = X.

We suppose that this transformation makes the domain Ω invariant and preserves the functional spaces, i.e.

APPROXIMATION OF A SOLIDIFICATION PROBLEM