Analysis and optimal control of some solidification processes

DISCRETE AND CONTINUOUS doi:10.3934/dcds.2014.34.3985 DYNAMICAL SYSTEMS Volume 34, Number 10, October 2014 pp. 3985–4017 ANALYSIS AND OPTIMAL CONTROL OF SOME SOLIDIFICATION PROCESSES Roberto C. Cabrales Departamento de Ciencias B´asicas, Universidad del B´ıo-B´ıo Campus Fernando May, Casilla 447 Chill´an, Chile Gema Camacho and Enrique Fern´ andez-Cara Departamento de Ecuaciones Diferenciales y An´alisis Num´erico Universidad de Sevilla, Aptdo. 1160 41080 Sevilla, Spain Abstract. In this paper we consider a mathematical model that describes the solidification of a binary alloy. We prove some existence and uniqueness results for a regularized problem, depending on a small parameter ε. We also analyze the behavior of the regularized solutions as ε→0. Then, we consider some associated optimal control problems. We prove existence and optimality results and we present and discuss some iterative methods. 1. Introduction. The mathematical model. The solidification of metals is one of the most difficult problems to model and analyze in engineering. The complex nature of the processes is an obstacle to obtain solutions in a simple way. Indeed, to accomplish this task, we have to use sophisticated mathematical models based on systems of nonlinear partial differential equations. However, by using correctly the laws governing solidification processes and appropriate mathematical and numerical techniques, it is possible to gain insight in the phenomena, obtain products of good quality and, also, diminish the related costs. This paper deals with the theoretical analysis and optimal control of a model of solidification for a binary alloy. The plan is the following. In the remainder of this Section we present the considered model and also a regularized version. Then, in Section 2we study the existence and uniqueness of the weak solution to the regularized problem, as well as its behaviour as the regularizing parameter εgoes to zero. Section 3deals with a first optimal control problem for the regularized solidification problem; here, the control variable is a heat source. We deduce the existence of optimal controls and we present the associated optimality system; some iterative algorithms are also introduced and discussed. Finally, in Section 4, we consider and solve a second optimal control problem where the solidification time plays a crucial role. 2010 Mathematics Subject Classification. Primary: 93C20, 35Q93, 49K20, 35A01, 35A02; Secondary: 76D05. Key words and phrases. Nonlinear PDEs, solidification models, regularization, weak solutions, optimal control. The first author is partially supported by Universidad del B´ıo-B´ıo grant GI121909/C (Chile), FONDECYT grant 1140074 (Chile) and by MICINN grant MTM2012-32325 (Spain). The second and third authors were supported by MICINN grant MTM2010-15992 (Spain). 3985 3986 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA 1.1. The example of a cavity. In order to present the basic ideas of solidification processes, we consider a situation where the mould is modelled as a rectangular cavity Ω filled with an incompressible and diluted binary alloy, initially at liquid state, with uniform temperature and composition, under the influence of gravity. The components of the mould are the solute and the solvent; together, they form the melting. The variable that indicates the proportion of solute in the melting is called the concentration. The alloy is deposited in the mould in such a way that, at time t= 0, the temperature of the left side of Ω is instantaneously dropped and kept under the cooling point, while the other sides of Ω are kept thermically isolated (see Figure 1.1). This is the origin of a temperature gradient in the alloy. The coupled action of gravity forces, convection, diffusion and reaction may introduce changes in the alloy density. On the other hand, the region next to the left side of Ω experiences a phase change, passing from the liquid to the solid state. The corresponding interface is called the solidification front. In general, it may have a highly irregular geometry and may exhibit dendritic forms for rapid solidification. The formation of dendrites is a consequence of the instability of growth of the solidification front. Remark 1. Solidification is a multi-scale problem, varying from meters to micrometers (and even nanometers). Because of this, it is very difficult to present a model involving all the relevant variables and phenomena occurring at each scale. In this paper, we will consider a model that can capture the dynamics of the phenomena at the macroscopic (mould) scale. To describe appropriately the very rich and complex dynamics of dendritic growth, it is necessary to consider other models (see for instance [2], [15] and the references therein). On the other hand, in order to avoid some technical (and maybe nontrivial) difficulties, latent heat effects have been neglected. After a while, the alloy is solidifies. If we analyze the composition of the resultant material, we will probably find variations of the solute. These variations are produced at two different scales: from the macroscopic viewpoint, the solidification front acts like a filter on the solute, that is rejected from the solid to the liquid state (this phenomenon is called macro-segregation); at another (microscopic) scale, the solute is partially trapped by the dendrites of the solidification front, generating highly concentrated solute grains (this is the micro-segregation phenomenon). Figure 1. A schematic view of the cavity problem, with details of the solidification front. 1.2. The solidification problem. In the sequel, Ω is a connected, open and bounded open set in Rd(d= 2,3) with Lipschitz-continuous boundary ∂Ω, such ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 3987 that ∂Ω = ΓD∪ΓN, ΓD∩ΓN=∅and ΓN6=∅. Let us assume that T > 0 and let us set Q:= Ω ×(0, T). Since it is very difficult to work with the highly irregular geometry of solidification front, it is convenient to perform an averaging process on the main variables of the problem in a such way that the front is replaced by a smooth phase change zone called the mushy or moisture zone, where we have coexistence of the liquid and solid states. Then, an acceptable mathematical model for the solidification problem is the following (see for example [3,4,5]): ct+∇ · (vcl(c, θ)−D∇c)=0,(1) θt+∇ · (vθ−χ∇θ) = hθ,(2) vt+∇ · (v⊗v−ν∇v) + Fi(c, θ)v+∇p=Fe(c, θ),(3) ∇ · v= 0.(4) In (1)–(4), hθis a known function and, for any aand b,a⊗bstands for the tensor whose (i, j) component is aibj;cand clare respectively the concentration and the liquid concentration of solute of the binary alloy; θis the temperature of the moisture; vis the velocity of the liquid and pis the pressure; D,χand νare positive constants, respectively denoting the diffusion coefficient of the solute, the thermal conductivity and the kinematic viscosity; Fiand Feare functions associated to the internal and external forces acting on the system (1)–(4). In general, Fiand Feare bounded functions depending on cand θ. By introducing the solid fraction fs(the function denoting the proportion of solid in the mould) and using the Carman-Kozeny approximation to model the effects in the mould as a porous media and the Boussinesq approximation to model the thermic and solutal stresses, we get the following expressions: Fi=M0 fs(c, θ)2 (1 −fs(c, θ))3,(5) Fe=g(1 + βccl(c, θ) + βθθ).(6) Here, M0denotes a positive constant depending on the material, gis the gravity force, βθand βcare constants representing the thermal and solute expansion coefficients and one has c=csfs+ (1 −fs)cl, where csis the solid concentration of solute, given by cs=rclwith 0 < r < 1. 1.3. The phase diagram. Assume that t≥0. By using the solid fraction fs, we can identify the following sets: Ωs(t) = {x∈Ω : fs(x, t) = 1},Ωl(t) = {x∈Ω : fs(x, t) = 0}, Ωm(t) = {x∈Ω : fs(x, t)∈(0,1)}. This allows the solid and non-solid regions to be defined as Qs={(x, t) : x∈Ωs(t), t ∈(0, T)} and Qml ={(x, t) : x∈Ωm(t)∪Ωl(t), t ∈(0, T)}, respectively; here, “ml” stands for “moisture-liquid”. The equations (1) and (2) must hold in the whole cylinder Q, with v= 0 in Qs. The equations (3) and (4) must hold in Qml. 3988 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA The sets Qml and Qscan be computed from the phase diagram by using the values of cand θ(see Figure 2). In this diagram, θfand θedenote the fusion and eutectic temperature, respectively. By simplicity, we assume that the solid and liquid curves θ=ηs(c) and θ=ηl(c) are straight lines; they respectively indicate the couples (c, θ) where the solid alloy begins to melt and the couples where the liquid begins to solidify. Specifically, we have: ηl(c) = θf−(θf−θe)c ce and ηs(c) = θf−(θf−θe)c ca , whee caand ceare characteristic concentration values. Also, we can identify in the phase diagram the following sets: L={(c, θ) : c > ceor θ > ηl(c)},(7) M={(c, θ):0≤c≤ce,max(ηs(c), θe)≤θ≤ηl(c)},(8) S={(c, θ):0≤c≤ce, θ < max(ηs(c), θe)}.(9) Figure 2. A typical phase diagram for a binary alloy The temperature and the concentration can be used to determine the solid fraction fsthrough the following equalities: fs(c, θ) =                0,if c≥ce, 0,if 0 ≤c≤ce, θ > ηl(c), ηl(c)−θ ηl(c)−max (ηs(c), θe),if 0 ≤c≤ce,and max (ηs(c), θe)≤θ≤ηl(c), 1,otherwise. (10) Taking into account this definition of fs, we can deduce that the liquid concentration is given by: cl(c, θ) =          ce,if c≥ce, c, if 0 ≤c≤ce, θ > ηl(c), (1 −fs)c, if 0 ≤c≤ce,max (ηs(c), θe)≤θ≤ηl(c), 0,otherwise. (11) In the sequel, we will always assume that hθ=h1ωfor some (small) non-empty open set ω⊂Ω (1ωis the characteristic function of ω). In practice, this means ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 3989 that we only act on the system through the right hand side of the transport-heat equation (2), by imposing a heat source localized in ω×(0, T). In fact, it would have been more realistic to assume that the action is performed through (a part of) the boundary, for instance by imposing a non-zero Dirichlet condition on θ. However, this leads to a more complex analysis involving technical difficulties and we have preferred to consider the present situation and deal with distributed or internal controls. 1.4. The regularized problem. Notice that the interfaces that separate the solid and the non-solid region are not known a priori; this is a serious difficulty. Moreover, the Carman-Kozeny term is singular in Qs, because the solid fraction fsreaches 1 there. Accordingly, we will introduce a parameter ε∈(0,1] and we will replace the function Fi(c, θ) defined in (5) by the following: Fε i=M0 fs(c, θ)2 (1 −fs(c, θ) + ε)3.(12) This leads to the so called regularized problem: ct+∇ · (vcl(c, θ)−D∇c)=0,(13) θt+∇ · (vθ−χ∇θ) = h1ω,(14) vt+∇ · (v⊗v−ν∇v) + Fε i(c, θ)v+∇p=Fe(c, θ),(15) ∇ · v= 0.(16) Now, the four equations (13)–(16) must hold in the whole set Q. The system is completed with boundary conditions for c,θand von ∂Ω×(0, T): (D∇c)·n= 0 on ∂Ω×(0, T),(17) (χ∇θ)·n= 0 on ΓN×(0, T),(18) θ= 0 on ΓD×(0, T),(19) v=0on ∂Ω×(0, T).(20) Finally, we add initial conditions: c(x,0) = c0(x); θ(x,0) = θ0(x); v(x,0) = v0(x) in Ω.(21) Let us now justify the introduction of (13)–(16). Assume that, as ε→0, the associated solutions to (13)–(16) remain bounded. For small ε, the term Fε ibecomes very large in the solid part of the domain, where fs≈1; thus, in this set the velocity field is at least O(ε3) and converges to zero. Contrarily, in any open set O ⊂⊂ Qml, Fε iremains bounded and converges in some sense to Fi. This gives the motion equation (3) in the limit. This explanation clarifies the role of Fε iin the regularized problem: Fε iis a coefficient that blows up in the solidified part of the domain and remains bounded in the other part (there, it is related to dissipation). These properties of Fε iwill be established rigorously when d= 2 in Section 2.3. 2. Existence and uniqueness results. In this Section, we analyze the existence of weak solutions to the regularized solidification problem. The problem is given by equations (13)–(16) and the boundary and initial conditions are given by (17)– (21). We will also prove a uniqueness result when d= 2. Finally, we will show that, at least when d= 2, the solutions to the regularized problems converge (in an appropriate sense) towards a weak solution to the original problem. 3990 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA The proofs of the results that follow are more or less standard. Maybe, the most interesting points are the following: •The estimates 0 ≤c≤ce. •The energy estimates of cε,θεand vε, that are independent of ε. •The fact that cεand θεare actually strong solutions if c0and θ0are regular enough and •The estimates of cεand θεin a space of H¨older-continuous functions in Q when d= 2. Some previous works concerning this topic are [6,10,9]. 2.1. Preliminaries and weak formulation. In the sequel, we consider the usual Sobolev spaces, given by Wm,p(Ω) = {f∈Lp(Ω) : ∂αf∈Lp(Ω) ∀|α| ≤ m} for all m≥1 and 1 ≤p≤+∞. When p= 2, we write Hm(Ω) instead of Wm,2(Ω); we will denote by Hm 0(Ω) the closure of C∞ 0(Ω) in Hm(Ω). Sometimes, we will also need to work with the Sobolev spaces Wr,p(Ω), where r∈R+; for the definitions and main properies, see for instance [1]. If there is no confusion, we write Lpinstead of Lp(Ω), Hminstead of Hm(Ω), etc. For any τ∈[0,1), we will denote by C0,τ,τ/2(Q) the space of functions ϕ∈C0(Q) that are H¨older-continuous in the following sense: sup (x,t),(x0,t0)∈Q |ϕ(x, t)−ϕ(x0, t0)| |x−x0|τ+|t−t0|τ/2<+∞. In general, we will denote by k·kXthe norm of the normed space X. If Xis a Banach space and 1 ≤q≤+∞, we denote by Lq(0, T;X) the Banach space of the X-valued (classes of) functions defined on the interval [0, T ] that are Lq-integrable in the sense of Bochner. On the other hand, C0([0, T]; X) and C0 w([0, T]X) will respectively denote the spaces of continuous and weakly continuous functions f: [0, T]7→ X. Sometimes we write Lq(X) instead of Lq(0, T;X), etc. Spaces of Rd-valued functions, as well as their elements, are usually denoted by bold faced letters. We will consider the space V(Ω) = {v∈C∞ 0(Ω) : ∇ · v= 0 in Ω }. Let Hand Vrespectively stand for the closures of V(Ω) in L2(Ω) and H1(Ω). It is then possible to show that H={v∈L2(Ω) : ∇ · v= 0 in Ω,v·n= 0 on ∂Ω}, V={v∈H1 0(Ω) : ∇ · v= 0 in Ω }, where nis the outward unit normal vector to ∂Ω, see for instance [19]. In the sequel, we denote by Ca generic positive constant depending only on Ω and the other data of the problem. As usual, it may have different values in different expressions. We will sometimes emphasize that the constants may have different values by putting C1, C2, . . . We will emphasize that Cdepends on a specific data D(a set, a parameter, . . . ), by writing C(D). ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 3991 For the weak formulation of the regularized solidification problem described by (13)–(21), we will also need the following space: Hθ(Ω) = {ψ∈H1(Ω) : ψ= 0 on ΓD}. Notice that the function clin (11) satisfies cl∈W1,∞ loc (R2).(22) This is a consequence of the fact that fsis Lipschitz-continuous. Moreover, cl≡ce for c≥ceand cl≡0 for c≤0. From now on, we will assume that h∈L2(ω×(0, T)),(23) c0∈L∞(Ω),0≤c0(x)≤cea.e. in Ω,(24) θ0∈L2(Ω),(25) v0∈H,(26) g∈L∞(Ω).(27) Definition 2.1. We say that (c, θ, v) is a weak solution of (13)–(21) if c∈L2(0, T;H1(Ω)) ∩L∞(0, T;L2(Ω)),(28) θ∈L2(0, T;Hθ(Ω)) ∩L∞(0, T;L2(Ω)),(29) v∈L2(0, T;V)∩L∞(0, T;H),(30) the following relations hold in (0, T ) hct, ϕi+ZΩ v· ∇cl(c, θ)ϕ+ZΩ D∇c· ∇ϕ= 0,(31) hθt, ψi+ZΩ v· ∇θψ +ZΩ χ∇θ· ∇ψ=Zω hψ, (32) hvt,wi+ZΩ (v· ∇)v·w+ZΩ ν∇v· ∇w+ZΩ Fε i(c, θ)v·w=ZΩ Fe·w,(33) for all ϕ∈H1(Ω), ψ∈ Hθ(Ω) and w∈Vand c|t=0 =c0, θ|t=0 =θ0,v|t=0 =v0.(34) Let us assume that (23)–(27) holds. It can be shown that any (c, θ, v) satisfying (28)–(30) and (31)–(33) also satisfies c∈C0([0, T]; L2(Ω)) and, furthermore, θ∈ C0([0, T]; L2(Ω)) and v∈C0([0, T]; H) if d= 2 and θ∈C0 w([0, T]; L2(Ω)) and v∈C0 w([0, T]; H) if d= 3; the argument is given below, in the proof of Theorem 2.2. In particular, the initial conditions (34) make sense. Another property satisfied by cis the following. Let us introduce the quantity cg:= ZΩ c0(x), i.e. the initial total amount of solute in the moisture. By using (24), we see that 0≤cg≤ce|Ω|<+∞. Furthermore, if we take ϕ= 1 in (31), a short computation shows that d dtZΩ c(x, t) = −Z∂Ω (v·n)cldΓ=0 ∀t∈(0, T). 3992 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA Thus, ZΩ c(x, t) is constant in time and ZΩ c(x, t) = cg=ZΩ c0(x)∀t∈(0, T). In other words, the amount of solute in the moisture does not depend on time. 2.2. On the existence and uniqueness of solution of the regularized problem. In the proof of existence of a weak solution to the regularized problem, we will use the following result: Proposition 1. Assume that (c, θ, v)is a weak solution to (13)–(21). Then 0≤c(x, t)≤cea.e. in Q. Proof. Let us denote by ϕ+(resp. ϕ−) the positive (resp. negative) part of ϕ. For each t∈[0, T ], let us choose ϕ=c−(·, t) in (31). Since c=c+−c−, taking into account the definition of cl, we easily get: −1 2 d dtkc−k2 L2−Dk∇c−k2 L2= 0. Integrating in time, we have −1 2kc−(·, t)k2 L2−DZt 0 k∇c−k2 L2ds =−1 2k(c0)−k2 L2= 0. Therefore, kc−(·, t)kL2= 0 for all t, that is, c≥0. Now, for each t∈[0, T], we choose ϕ= (c−ce)+(·, t) in (31). By using an argument similar to the previous one, we obtain that 1 2k(c(·, t)−ce)+k2 L2+DZt 0 k∇(c−ce)+k2 L2ds =1 2k(c0−ce)+k2 L2= 0 for all t∈[0, T], whence c≤ce. Our first result concerns the existence of a solution to (13)–(21): Theorem 2.2. There exists at least one weak solution (c, θ, v)to the regularized problem (13)–(21), with 0≤c≤ce. Proof. We will introduce Galerkin approximations. To this end, we will consider three “special” bases Bc={ϕk(x) : k∈N},Bθ={ψk(x) : k∈N}and Bv={wk(x) : k∈N}, respectively in H1(Ω), Hθ(Ω) and V. It will be assumed that they are orthogonal for the scalar product in H1and orthonormal for the scalar product in L2. Obviously, this is the case if the ϕkare the eigenfunctions of the Neumann Laplacian in Ω, that is,    −∆ϕk=λkϕk, x ∈Ω, ∇φk·n= 0, x ∈∂Ω, kφkkL2= 1, and similar definitions hold for ψkand wk. Let us fix m∈N. We consider the m-dimensional spaces Sm c,Sm θand Sm v, respectively spanned by the first mfunctions of Bc,Bθand Bv. ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 3993 For any t, we define the approximations cm,θmand vmas follows: cm(x, t) = m X k=1 λk,m(t)ϕk(x), θm(x, t) = m X k=1 ξk,m(t)ψk(x), vm(x, t) = m X k=1 σk,m(t)wk(x). The coefficients λk,m(t), ξk,m(t) and σk,m(t) are computed in such way that hcm t, ϕki+ZΩ vm· ∇cl(cm, θm)ϕk+ZΩ D∇cm· ∇ϕk= 0,(35) hθm t, ψki+ZΩ vm· ∇θmψk+ZΩ χ∇θm· ∇ψk=Zω hψk,(36) hvm t,wki+ZΩ (vm· ∇)vm·wk+ZΩ ν∇vm· ∇wk +ZΩ Fε i(cm, θm)vm·wk=ZΩ Fe(cm, θm)·wk, (37) for all k= 1, . . . , m and cm(·,0) = Pc,m(c0), θm(·,0) = Pθ,m(θ0),vm(·,0) = Pv,m(v0).(38) Here, Pc,m :H1(Ω) 7→ Sm c,Pθ,m :Hθ(Ω) 7→ Sm θand Pv,m :V7→ Sm vare the orthogonal projectors. For each m, (35)–(37) is an ordinary differential system for the unknowns λk,m(t), ξk,m(t) and σk,m(t). It is complemented with the initial conditions (38). This initial value problem has a local in time solution (cm, θm,vm) defined in some interval [0, tm). In order to prove that tm=T, we need to show some a priori estimates. If we multiply equations (35), (36) and (37) by λk,m(t), ξk,m(t) and σk,m(t), respectively, and we sum over k, we obtain (using that vm∈V): 1 2 d dtkcmk2 L2+Dk∇cmk2 L2=−ZΩ vm· ∇cl(cm, θm)cm 1 2 d dtkθmk2 L2+χk∇θmk2 L2=Zω hθm 1 2 d dtkvmk2 L2+νk∇vmk2 L2=ZΩ Fe(cm, θm)·vm−ZΩ Fε i(cm, θm)|vm|2 Now, we have to estimate the terms in the right-hand side of these inequalities. By using (22) and H¨older and Young inequalities, we can bound the first term as follows: ZΩ vm· ∇cl(cm, θm)cm=−ZΩ cl(cm, θm)vm· ∇cm ≤D 2k∇cmk2 L2+C1kvmk2 L2. (39) 4000 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA To this end, we first notice that Qml can be written as a countable union of open sets of the form O × I, with O ⊂ Ω, I ⊂(0, T) (an open interval), sup O×I fs(c, θ)<1.(72) Indeed, we first have Qml =[ n≥1 Q(n) ml ,with Q(n) ml ={(x, t)∈Q:fs(c(x, t), θ(x, t)) <1−1 n}; then, any Q(n) ml can be written as a countable non-decreasing union of compact sets K(n) m; finally, each K(n) mcan be covered by a finite union of sets O × Isatisfying (72). Consequently, it will suffice to check that, for any such set O × I, we have vε→va.e. in O × I(73) at least for a subsequence. Thus, let Oand Ibe given, satisfying (72). Let us introduce the space V(O), with V(O) = {w∈H1 0(O) : ∇ · w= 0 in O } (74) and let us denote by V(O)0the dual of V(O). Let uεdenote the restriction to O × Iof the field vε. We obviously have kuεkL2(I;H1(O)) ≤C. (75) On the other hand, for any w∈V(O) and any t∈I, one has: huε t,wi=ZΩ vε t·w=ZΩ(vε· ∇)w·vε+ν∇vε· ∇w−(Fε e−Fε ivε)·w ≤C(O × I)1 + kvεk2 L4+k∇vεkL2kwkV(O), where we have used (69), (70) and the fact that, for any sufficiently small ε,Fε iis uniformly bounded in O × I. These estimates show that kuε tkV(O)0≤C(O × I)1 + kvεk2 L4+k∇vεkL2 and, consequently, kuε tkL2(I;V(O)0)≤C(O × I).(76) From (75) and (76), we deduce from well known compactness results that uε belongs to a compact set in L2(O × I). Therefore, we get (73) at least for a subsequence, as desired. This proves (71). From (70) and (71), we obtain: vε→vstrongly in Lp(Qml)∀p∈[1,4). On the other hand, (69) and (71) give (vε· ∇)vε→(v· ∇)vweakly in L4/3(Qml).(77) We can now take limits in the equations satisfied by cε,θεand vε. Thus, using (66)–(77), we find that cand θsatisfy (31), (32) and the first two equalities of (34). ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 4001 Also, using (69), (70) and (68), the following is found for every w∈C∞ 0(Q) with Supp w⊂Qml and ∇ · w≡0: ZZQ (−v·wt+ (v· ∇)v·w+ν∇v· ∇w+ (Fi(c, θ)v−Fe(c, θ)) ·w) = lim ε→0ZZQ (−vε·wt+ (vε· ∇)vε·w+ν∇vε· ∇w+ (Fε ivε−Fε e)·w) = 0. In order to check that vt=0 =v0in Ωml(0), we argue as follows. Let O ⊂⊂ Ωml(0) be a non-empty open set and let us denote again by V(O) the space (74). Let ρ > 0 be such that O × [0, ρ]⊂Qml ∪(Ω × {0}); the definition of Qml and the fact that cand θare uniformly continuous prove that such a ρexists. Let us denote by uε(resp. u) the restriction of vε(resp. v) to O × (0, ρ) and let u0 be the restriction of v0to O. Then kuεkL2(0,ρ;H1(O)) +kuεkL∞(0,ρ;L2(O)) ≤C and an argument similar to the proof of (76) shows that kuε tkL2(0,ρ;V(O)0)≤C(O × (0, ρ)). Hence, uεcan be regarded as a continuous V(O)0-valued function and, also, as a weakly continuous L2(O)-valued function. Furthermore, kuεkC0([0,ρ];V(O)0)≤C(O × (0, ρ)) and k(uε, ϕ)L2(O)kC0([0,ρ]) ≤C(O × (0, ρ), ϕ)∀ϕ∈L2(O). This shows that ucan also be regarded as a L2(O)-valued weakly continuous function and uε|t=0 converges weakly in L2(O) to u|t=0. But we also have uε|t=0 =u0. Therefore, u|t=0 =u0in Oand, since Ois arbitrary, the desired initial condition is satisfied by v. Let us finally prove that v= 0 in int Qs. To this purpose, we will view again (vε, pε) as a weak solution to (15) satisfying (62). From De Rham’s lemma, we have ∇pε∈W−1,∞(0, T ;H−1(Ω)) and k∇pεkW−1,∞(0,T ;H−1(Ω)) ≤C. Indeed, we can write (15) in the form hSε,ϕi= 0 ∀ϕ∈ V(Ω), where Sε:= vε t+∇ · (vε⊗vε−ν∇vε) + Fε i(cε, θε)vε−Fe(cε, θε) is uniformly bounded in W−1,∞(0, T;L2(Ω)) ∩L2(0, T;H−1(Ω)). Accordingly, S can be viewed as a distribution in L2(Ω; W−1,∞(0, T)) that vanishes when it is applied to any ϕin V(Ω), see for instance [18]; therefore, Sεmust be a gradient, that is, Sε=−∇pεfor some ∇pεuniformly bounded in H1(Ω; W−1,∞(0, T)) ∼ = W−1,∞(0, T;H−1(Ω)). Thus, it can be assumed that ∇pε→ ∇pweakly-∗in W−1,∞(0, T ;H−1(Ω)).(78) Let ϕ∈C∞ 0(Q) be given, with K:= Supp ϕ⊂int Qs. Notice that fs(cε, θε)→1 and ε3Fε i→1 uniformly in Kas ε→0.(79) 4002 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA Also, 0 = ε3ZZQ (vε t+ (vε· ∇)vε−ν∆v+Fε ivε+∇pε−Fε e)·ϕ =ε3ZZQ (vε t+ (vε· ∇)vε+∇pε−Fε e)·ϕ +νε3ZZQ ∇vε· ∇ϕ+ε3ZZQ Fε ivε·ϕ Taking limits as ε→0 and using (65), (69), (78) and (79), we deduce that ZZQ v·ϕ= 0. Since this must hold for all ϕ∈C∞ 0(Q) with Supp ϕ⊂int Qs, we find that v vanishes in int Qs. This shows that (c, θ, v) solves the solidification problem in the sense of definition 2.4 and ends the proof. Remark 2. The regularity assumptions on c0and θ0have been needed to ensure the uniform H¨older-continuity of cεand θεand thus the compactness of these functions in C0(Q). On the other hand, the last passage to the limit in the previous proof relies, among other things, on the fact that ZZQ vε t·ϕ=−ZZQ vε·ϕt and consequently converges as ε→0. Remark 3. Observe that, if d= 3 and (cε, θε,vε) is a weak solution to the regularized problem for every ε, we can only prove that vεis uniformly bounded in L10/3(Q), see [8,19]. But this is not enough to obtain (63) and (64) and we cannot deduce that cεand θεbelong to a compact set of C0(Q); we would need a uniform estimate in L5/2(Q), see [13]. It is at this point that the argument begins to fail in the three-dimensional case. Consequently, the analog of theorem 2.5 in the three-dimensional case is open. 3. Some first optimal control problems for the regularized solidification model. Let us consider some first (standard) control problems for the regularized solidification system (13)–(21). The underlying goal is to govern the growth of a solidification front (and its geometrical shape) by imposing a heating mechanism that acts on the non-empty subset ω⊂Ω. A mathematical formulation can be obtained under the form of an optimal control problem. More precisely, we will try to minimize an appropriate cost function subject to (13)–(21) and some additional constraints on the control h. In this Section, we will propose several choices for the cost functional (see the definitions (81) and (82) below). In each case, our purpose will be to achieve three main tasks: •To prove the existence of an optimal control. •To characterize the optimal controls or, at least, to obtain necessary conditions for optimality. •To provide iterative algorithms for the computation of optimal controls. ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 4003 A very interesting question is what happens when ε→0, i.e. which is the behavior of the solutions to the optimal control problems when we take limits in (13)–(21) and we approach the original solidification model. This is a difficult question. It will not be analyzed in this paper, but will be the objective of a forthcoming work. Let J=J(h, c, θ, v) be a cost function. The considered optimal control problem will have the following general form: (Find (h∗, c∗, θ∗,v∗)∈ E such that J(h∗, c∗, θ∗,v∗) = min (h,c,θ,v)∈E J(h, c, θ, v),(80) where Eis a non-empty set that will be specified later. In practice, any (h, c, θ, v)∈ E will be assumed to satisfy (13)–(21). Furthermore, it will be realistic to consider appropriate constraints on h, i.e. to impose that h belongs to a set of admissible controls Uad ⊂L2(ω×(0, T)). Let α,β,γbe nonnegative constants and let us assume that N > 0. We will consider the following two possible choices of J, that seem reasonable: First Choice: Let cd, θd∈L2(Q) and vd∈L2(Q) be given. We set J(h, c, θ, v) = α 2ZZQ |c−cd|2+β 2ZZQ |θ−θd|2+γ 2ZZQ |v−vd|2 +N 2ZZω×(0,T ) |h|2 (81) Second Choice: Let ce, θe∈L2(Ω) and ve∈L2(Ω) be given. We now set J(h, c, θ, v) = α 2ZΩ |c(x, T)−ce(x)|2+β 2ZΩ |θ(x, T)−θe(x)|2 +γ 2ZΩ |v(x, T)−ve(x)|2+N 2ZZω×(0,T ) |h|2 (82) As usual, these objective functions correspond to the goal of having (c, θ, v) close to a desired state in the cheapest feasible way. 3.1. The existence of optimal controls. Let us fix a non-empty set Uad ⊂ L2(ω×(0, T)). Let us set E={(h, c, θ, v) : h∈ Uad,(c, θ, v) solves (13)–(21)}. Notice that Eis a non-empty subset of L2(ω×(0, T ))×E, where Eis the energy space for the weak solutions to (13)–(21), i.e. the space of triplets (c, θ, v) satisfying (28)–(30). Of course, Eis a Banach space for the norm kckL2(H1)+kckL∞(L2)+kθkL2(Hθ)+kθkL∞(L2)+kvkL2(V)+kvkL∞(H). Let us consider the control problem (80). The following result holds: Theorem 3.1. Assume that Uad is weakly closed in L2(ω×(0, T)) and the following hypotheses are satisfied: 1. Either Jis a coercive functional, that is, J(hn, cn, θn,vn)→+∞if (hn, cn, θn,vn)∈ E,khnkL2(ω×(0,T )) →+∞, or Uad is a bounded set. 4004 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA 2. Jis sequentially weakly-∗lower semicontinuous, that is, if (hn, cn, θn,vn)→ (h, c, θ, v)weakly-∗in L2(ω×(0, T)) ×E, then lim inf n→+∞J(hn, cn, θn,vn)≥J(h, c, θ, v). Then, there exists at least one solution to (80), that is, an optimal control and its associated state. Proof. The argument is standard and well known. However, for completeness, we will give a sketch. Let {(hn, cn, θn,vn)}be a minimizing sequence for (80). Then all the controlstates (hn, cn, θn,vn) belong to Eand the hnare uniformly bounded in L2(ω× (0, T)). Indeed, if this were not the case, it could be assumed that khnkL2(ω×(0,T )) → +∞, whence we should have J(hn, cn, θn,vn)→+∞, which is an absurd. From the estimates in the proof of theorem 2.2, we easily get 0≤cn≤ce,k(cn, θn,vn)kE≤C, and kcn tkL2(0,T ;(H1(Ω))0)+kθn tkLσ(0,T ;(Hθ(Ω))0)+kvn tkLσ(0,T ;V0)≤C, where σ= 2 if d= 2 and σ= 4/3 if d= 3. Consequently, at least for a new minimizing (sub)sequence, we have the weak and/or weak-∗convergence of hn,cn,θn,vn,cn t,θn tand vn tin appropriate spaces and the strong convergence of cn,θnand the components of vnin L2(Q). Thus, as in the proof of theorem 2.2, we deduce that the limit (h∗, c∗, θ∗,v∗) belongs to E, i.e. h∗∈ Uad (here we use that Uad is weakly closed) and (c∗, θ∗,v∗) solves (13)–(21) for h=h∗. We also have lim inf n→+∞J(hn, cn, θn,vn)≥J(h∗, c∗, θ∗,v∗). Hence, (h∗, c∗, θ∗,v∗) solves (80) and the proof is achieved. Remark 4. The hypotheses in this result are satisfied when Jis given by (81) or (82) and Uad ⊂L2(ω×(0, T)) is non-empty, closed and convex. Typical examples are the following: Uad =L2(ω×(0, T)) Uad ={h∈L2(ω×(0, T)) : |h| ≤ Ra.e. } Uad ={h∈L2(ω×(0, T)) : h= m X j=1 hj(x)1Ij(t), hj∈L2(ω)} where R > 0 and the Ij⊂(0, T) are disjoint intervals. Consequently, in these cases, there exists at least one solution to (80). 3.2. The optimality system. Next, we will deduce the optimality system associated to the previous control problems. We will need a regularity assumption on the optimal control-state (h∗, c∗, θ∗,v∗): The set of points (x, t)where θ∗=ηl(c∗)or max(ηs(c∗), θe)is negligible. (83) In the case of the cost functional (81), the following holds: ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 4005 Theorem 3.2. Assume that Uad ⊂L2(ω×(0, T )) is non-empty, closed and convex and Jis given by (81). Let (h∗, c∗, θ∗,v∗)be an optimal solution to (80)satisfying (83)and assume that (c∗, θ∗,v∗)is the unique state associated to h∗. Then, there exists (φ, ψ, w)∈Esuch that one has:                    c∗ t+v∗· ∇c∗ l−D∆c∗= 0, θ∗ t+v∗· ∇θ∗−χ∆θ∗=h∗1ω, v∗ t+ (v∗· ∇)v∗−ν∆v∗+ (Fε i)∗vε+∇pε=F∗ e, ∇ · v∗= 0, (D∇c∗)·n= 0,v∗= 0 on ∂Ω×(0, T), (χ∇θ∗)·n= 0 on ΓN×(0, T), θ∗= 0 on ΓD×(0, T ), c∗(x,0) = c0(x); θ∗(x,0) = θ0(x); v∗(x,0) = v0(x)in Ω. (84)                                −φt−(Dccl)∗v∗· ∇φ−D∆φ= ((DcFe)∗−(DcFε i)∗v∗)·w +α(c∗−cd), −ψt−v∗· ∇ψ−χ∆ψ= ((DθFe)∗−(DθFε i)∗v∗)·w +(Dθcl)∗v∗· ∇φ+β(θ∗−θd), −wt−ν∆w−(v∗· ∇)w+ (∇v∗)tw+ (Fε i)∗w+∇π =−φ∇c∗ l−ψ∇θ∗+γ(v∗−vd), ∇ · w= 0, (D∇φ)·n= 0,w= 0 on ∂Ω×(0, T), (χ∇ψ)·n= 0 on ΓN×(0, T), ψ = 0 on ΓD×(0, T), φ(x, T) = 0; ψ(x, T) = 0; w(x, T )=0 in Ω. (85) ZZω×(0,T ) (ψ+Nh∗)(h−h∗)≥0∀h∈ Uad, h∗∈ Uad.(86) Notice that the same conclusion is obtained if (h∗, c∗, θ∗,v∗) is only a local minimizer of J. Proof. Let us take h=h∗+am with a∈R+(small), m∈L2(ω×(0, T)) and h∗+am ∈ Uad. Let (c, θ, v) be a state associated to h. We can then write (c, θ, v)=(c∗, θ∗,v∗) + a(z, y, u) + a(z0 a, y0 a,u0 a), with            zt+v∗· ∇((Dccl)∗z+ (Dθcl)∗y) + u· ∇c∗ l−D∆z= 0, yt+v∗· ∇y+u· ∇θ∗−χ∆y=m1ω, ut+ (v∗· ∇)u+ (u· ∇)v∗−ν∆u+ (Fε i)∗u +((DcFε i)∗z+ (DθFε i)∗y)v∗+∇q= (DcFe)∗z+ (DθFe)∗y, ∇ · u= 0, (87) and            z0 a,t +v∗· ∇((Dccl)∗z0 a+ (Dθcl)∗y0 a) + u0 a· ∇c∗ l−D∆z0 a=−Za, y0 a,t +v∗· ∇y0 a+u0 a· ∇θ∗−χ∆y0 a=−Ya+ ((DcFε i)∗z0 a, u0 a,t + (v∗· ∇)u0 a+ (u0 a· ∇)v∗−ν∆u0 a+ (Fε i)∗u0 a +(DθFε i)∗y0 a)v∗+∇q= (DcFe)∗z0 a+ (DθFe)∗y0 a−Ua, ∇ · u0 a= 0. (88) Here, we have used the following notation: c∗ l=cl(c∗, θ∗),(Dccl)∗=Dccl(c, θ)|c=c∗,θ=θ∗,(Fε i)∗=Fε i(c∗, θ∗), . . . The functions (z, y, u) and (z0 a, y0 a,u0 a) must satisfy the same boundary conditions than (c∗, θ∗,v∗) and homogeneous initial conditions at t= 0. 4006 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA We have Za=∇ · ((cl(c, θ)−c∗ l)(u+u0 a) + Z0 a), where Z0 a=1 a(cl(c, θ)−c∗ l)−((Dccl)∗(z+z0 a)+(Dθcl)∗(y+y0 a))v∗ := [D1,a(z+z0 a) + D2,a(y+y0 a)] v∗ Similar expressions hold for Yaand Ua: Ya=a∇ · ((y+y0 a)(u+u0 a)) and Ua=a∇ · ((u+u0 a)⊗(u+u0 a)) + (Fε i(c, θ)−(Fε i)∗)(u+u0 a) + U0 a, where U0 a=βc1 a(cl(c, θ)−c∗ l)−((Dccl)∗(z+z0 a)+(Dθcl)∗(y+y0 a))g +1 a(Fε i(c, θ)−(Fε i)∗)−((DcFε i)∗(z+z0 a)+(DθFε i)∗(y+y0 a))v∗ := [D1,a(z+z0 a) + D2,a(y+y0 a)] g+ [D3,a(z+z0 a) + D4,a(y+y0 a)] v∗. Here, D1,a and D2,a (resp. D3,a and D4,a) denote appropriate combinations of the partial derivatives of cl(resp. Fε i). For instance, D1,a =Z1 0 Dccl(c∗+sa(z+z0 a), θ∗+sa(y+y0 a)) ds −(Dccl)∗.(89) Let us see that (z, y, u),(z0 a, y0 a,u0 a)∈E, with k(z, y, u)kE≤CkmkL2(ω×(0,T )),k(z0 a, y0 a,u0 a)kE→0 as a→0+.(90) Taking into account that (h∗, c∗, θ∗,v∗)∈ E, 0 ≤c≤ceand cl∈W1,∞(R2), from the usual energy estimates for linear parabolic systems, we easily deduce the first part of (90): by multiplying the first, second and third equations in (87) respectively by z,yand u, integrating in space and setting e:= (z, y, u) and µ:= min(D, χ, ν), we find that d dtkek2 L2+µk∇ek2 L2≤1 + k∇θ∗k2 L2+k∇v∗k2 L2kek2 L2+km1ωk2 L2 for all t∈(0, T). This leads to the desired estimate. Also, k(z0 a, y0 a,u0 a)kEis bounded, independently of a. More precisely, using the expressions of Za,Z0 a,Ya,Uaand U0 aand introducing e0 a:= (z0 a, y0 a,u0 a), we now have d dtke0 ak2 L2+µk∇e0 ak2 L2 ≤1 + k∇θ∗k2 L2+k∇v∗k2 L2+k∇yk2 L2+k∇uk2 L2ke0 ak2 L2+Ma for all t∈(0, T), where Ma=Pa+Qa+Ra+Sa, Pa:= CZΩ |cl(c, θ)−c∗ l|2|u|2, Qa:= CZΩ |Fε i(c, θ)−(Fε i)∗|2|u|2,(91) Ra:= CZΩ|D1,a|2|z|2+|D2,a|2|y|2(92) ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 4007 and Sa:= CZΩ|D3,a|2(|∇z|2+|∇v∗|2) + |D4,a|2(|∇y|2+|∇v∗|2).(93) First, since Ma≤C(1 + k∇zk2 L2+k∇yk2 L2+k∇vk2 L2) for all aand t, we deduce that k(z0 a, y0 a,u0 a)kEis bounded by a constant depending on kmkL2(ω×(0,T )). Secondly, observe that this yields uniform bounds for z0 a,t and y0 a,t respectively in Lσ(0, T; (H1(Ω))0) and Lσ(0, T; (Hθ(Ω))0). Indeed, it suffices to consider the first two equations in (88) and notice that Zaand Yaare uniformly bounded in these spaces. As a consequence, at least for a subsequence, we have that z0 aand y0 a converge strongly in L2(Q). In view of (91) and Lebesgue’s theorem, we find that ZT 0 Pa(t)dt →0 and ZT 0 Qa(t)dt →0.(94) On the other hand, we also have ZT 0 Ra(t)dt →0 and ZT 0 Sa(t)dt →0.(95) Indeed, the uniqueness of (c∗, θ∗,v∗) and the regularity assumption (83) imply that, for instance, D1,a →0 a.e. in Q, since clis C1in a neighborhood of any (c, θ) with θ6=ηl(c) and θ6= max(ηs(c), θe). The same is true for D2,a,D3,a and D4,a. Consequently, Lebesgue’s theorem also leads to (95). From (94) and (95), we find that ZT 0 Ma(t)dt →0. Since ke0 ak2 L2(t) + µZt 0 k∇e0 ak2 L2(s)ds ≤CZt 0 Ma(s)ds ∀t∈[0, T], we deduce that, at least for a subsequence, the second part of (90) is fulfilled. Since this argument can be applied to any subsequence of {(z0 a, y0 a,u0 a)}, the convergence must hold for the whole sequence. This proves (90). By hypothesis, J(h, c, θ, v)−J(h∗, c∗, θ∗,v∗)≥0. Dividing by aand taking limits as a→0+, we see that ZZQ (α(c∗−cd)z+β(θ∗−θd)y+γ(v∗−vd)·u) + NZZω×(0,T ) h∗m≥0.(96) Let us introduce the linear (adjoint) system (85). From classical arguments, it is clear that (85) possesses at least one weak solution (φ, ψ, w)∈E. Furthermore, a straighforward integration by parts yields the following identity: ZZQ (α(c∗−cd)z+β(θ∗−θd)y+γ(v∗−vd)·u) = ZZω×(0,T ) ψ m. This, together with (96), gives the inequality ZZω×(0,T ) (ψ+Nh∗)m≥0. Since this must hold for any mof the form m=h−h∗with h∈ Uad, we find (86). This ends the proof. 4008 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA For the second choice of the cost functional, given by (82), a very similar reult can be obtained: Theorem 3.3. Assume that Uad ⊂L2(ω×(0, T )) is non-empty, closed and convex and Jis given by (82). Let (h∗, c∗, θ∗,v∗)be an optimal solution to (80)satisfying (83)and assume that (c∗, θ∗,v∗)is the unique state associated to h∗. Then, there exists (φ, ψ, w)∈Esuch that one has (84),                                −φt−(Dccl)∗v∗· ∇φ−D∆φ= ((DcFe)∗−(DcFε i)∗v∗)·w, −ψt−v∗· ∇ψ−χ∆ψ= ((DθFe)∗−(DθFε i)∗v∗)·w +(Dθcl)∗v∗· ∇φ, −wt−ν∆w−(v∗· ∇)w+ (∇v∗)tw+ (Fε i)∗w+∇π =−φ∇c∗ l−ψ∇θ∗, ∇ · w= 0, (D∇φ)·n= 0,w= 0 on ∂Ω×(0, T), (χ∇ψ)·n= 0 on ΓN×(0, T), ψ = 0 on ΓD×(0, T), φ(x, T) = α(c∗(x, T)−ce(x)); ψ(x, T ) = β(θ∗(x, T)−θe(x)) in Ω, w(x, T) = γ(v∗(x, T)−ve(x)) in Ω. (97) and (86). Again, notice that the same holds if (h∗, c∗, θ∗,v∗) is a local minimizer. Remark 5. Observe that theorems 3.2 and 3.3 do not assert that h7→ J(h, c, θ, v) is differentiable at h∗. In fact, nothing indicates that this function is well defined, since in general a control hclose to h∗can have several associated states. Nevertheless, we have been able to express the variation of Jat (h∗, c∗, θ∗,v∗) in the direction determined by min the form ZZω×(0,T ) (ψ+Nh∗)m, where ψsolves, together with φand w, the adjoint system (85) or (97). For this reason, we can interpret (ψ+Nh∗)|ω×(0,T )as the “gradient” of h7→ J(h, c, θ, v) at h∗. Remark 6. In the most simple case, Uad =L2(ω×(0, T)), and (86) means that h=−1 Nψin ω×(0, T).(98) More generally, since Uad is a closed convex set of L2(ω×(0, T)), (86) is equivalent to h=Pad −1 Nψ|ω×(0,T ),(99) where Pad :L2(ω×(0, T)) 7→ Uad is the orthogonal projector. 3.3. Some iterative algorithms. We will now propose some iterates to compute the solution to the previous optimal control problems. For simplicity, we will only refer to the case where Jis given by (81) and, consequently, the optimality system is (84)–(86). The adaptation to the case (82) is straightforward and will not be given. The following algorithms rely on the ideas in the proof of theorem 3.2. Specifically, we notice that, if (h, c, θ, v)∈ E and the couple (c, θ) is “regular” in the sense ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 4009 of (83), then for any m∈L2(ω×(0, T)), any small a > 0 and any state (c0, θ0,v0) associated to h0=h+am, one has J(h0, c0, θ0,v0) = J(h, c, θ, v) + aZZω×(0,T ) (ψ+Nh∗)m+a O(a), where ψis, together with φand w, the unique solution to                                −φt−(Dccl)v· ∇φ−D∆φ= ((DcFe)−(DcFε i)v)·w +α(c−cd) −ψt−v· ∇ψ−χ∆ψ= ((DθFe)−(DθFε i)v)·w +(Dθcl)v· ∇φ+β(θ−θd) −wt−ν∆w−(v· ∇)w+ (∇v)tw+ (Fε i)w+∇π =−φ∇cl−ψ∇θ+γ(v−vd) ∇ · w= 0 (D∇φ)·n= 0,w= 0 on ∂Ω×(0, T ) (χ∇ψ)·n= 0 on ΓN×(0, T ), ψ = 0 on ΓD×(0, T) φ(x, T) = 0; ψ(x, T ) = 0; w(x, T) = 0 in Ω (100) and O(a)→0 as a→0+. The first proposed algorithm is the following: Algorithm 1 a. Choose h0∈ Uad; b. Then, for given n≥0 and hn∈ Uad, do until convergence: 1. Solve (84) with h=hn, to obtain (cn, θn,vn); 2. Solve (100) with (c, θ, v)=(cn, θn,vn), to obtain (φn, ψn,wn); 3. Set dn= (ψn+Nhn)|ω×(0,T )and find ρnsuch that jn(ρn) = infρ>0jn(ρ). Here, jn(ρ) is the value of Jat any (hn−ρdn, cn(ρ), θn(ρ)vn(ρ)), where (cn(ρ), θn(ρ),vn(ρ)) is a state associated to hn−ρdn; 4. Set hn+1 =Pad(hn−ρndn). Table 1. The optimal step gradient method with projection. Let us assume that (84) possesses exactly one weak solution (c, θ, v) for each h∈ Uad (this is the case if d= 2) and that all the (h, c, θ, v)∈ E satisfy (83). Then algorithm 1must be viewed as a classical optimal step gradient method. Since (84) is nonlinear and we have to solve this system by using an iterative scheme, it is reasonable to introduce a variant where we perform mixed loops. This is described in Table 2. Remark 7. A natural choice of the convergence criteria can be kkn+1 −knkL2(ω×(0,T )) < κkkn+1kL2(ω×(0,T )), for κsmall enough. Notice however that this can be not completely significative for linear and not superlinear convergence. Consequently, this should be followed by an additional test where we check whether the necessary optimality conditions are satisfied. On the other hand, since the numerical computation of ρncan be 4016 R. C. CABRALES, G. CAMACHO AND E. FERN´ ANDEZ-CARA Proof of theorem 4.4.Let us introduce S∈R,m∈L2(ω×(0, T)) and a∈R+, such that T:= ˆ T+aS ∈[0, T0], h := ˆ h+am ∈ Uad,kθ(·, T)−θek2 L2=δ2(129) and let (c, θ, v) be a state associated to h. In view of (129), 0≤Φ(T, h)−Φ( ˆ T, ˆ h) = a ˆ TS +NZZω×(0,T0) ˆ h m! +a2 2"S2+NZZω×(0,T0) |m|2# for any small a > 0, whence ˆ TS +NZZω×(0,T0) ˆ hm ≥0.(130) As in the proof of theorem 3.2, we can write (c, θ, v) = (ˆc, ˆ θ, ˆ v) + a(z, y, u) + a(z0 a, y0 a,u0 a) with (z, y, u) and (z0 a, y0 a,u0 a) solving systems respectively similar to (87) and (88), together with homogenous boundary and initial conditions. Arguing as we did in that proof, we see that (z, y, u),(z0 a, y0 a,u0 a)∈E0and k(z0 a, y0 a,u0 a)kE0→0 as a→0+. Moreover, 0 = kθ(·, T)−θek2 L2−δ2=k(θ(·, T)−ˆ θ(·,ˆ T)) + (ˆ θ(·,ˆ T)−θe)k2 L2−δ2 =kθ(·, T)−ˆ θ(·,ˆ T)k2+ 2(θ(·, T)−ˆ θ(·,ˆ T),ˆ θ(·,ˆ T)−θe)L2 Taking into account that θ(·, T)−ˆ θ(·,ˆ T) = ay(·,ˆ T) + aˆ θt(·,ˆ T)S+a O(a) in L2(Ω) where O(a)→0, we easily deduce that −(ˆ θ(·,ˆ T)−θe,ˆ θt(·,ˆ T))L2S= (ˆ θ(·,ˆ T)−θe, y(·,ˆ T))L2.(131) Let us introduce λ∈Rwith −(ˆ θ(·,ˆ T)−θe,ˆ θt(·,ˆ T))L2λ=ˆ T(132) and let (φ, ψ, w)∈ˆ Ebe the solution to (116). Thanks to (114), λis well defined. Furthermore, ˆ TS =−(ˆ θ(·,ˆ T)−θe,ˆ θt(·,ˆ T))L2λS = (λ(ˆ θ(·,ˆ T)−θe), y(·,ˆ T))L2 = (ψ(·,ˆ T), y(·,ˆ T))L2 and, using the equations and boundary and initial conditions satisfied by (φ, ψ, w) and (z, y, u), we see after some integrations by parts that ˆ TS =ZZω×(0,T0) ψ m. In view of (130), this yieldsZZω×(0,T0) (ψ+Nˆ h)m≥0. ANALYSIS, OPTIMAL CONTROL AND SOLIDIFICATION 4017 On the other hand, since ˆ T∈(0, T0) and λis given by (132), the equality (118) is trivially satisfied. Consequently, the couple ( ˆ T, ˆ h), the associate state (ˆc, ˆ θ, ˆ v), the multiplier λ∈Rand the adjoint state (φ, ψ, w) satisfy (115)–(119). This ends the proof.  Remark 12. The optimality system can be used to deduce iterative algorithms for the computation of an optimal ( ˆ T, ˆ k, ˆc, ˆ θ, ˆ v). This will be the goal of a forthcoming paper. REFERENCES [1] R. A. Adams, Sobolev Spaces, Pure and Applied Mathematics, 65, Academic Press, New York-London, 1975. [2] W. J. Boettinger, J. A. Warren, C. Beckermann and A. Karma, Phase-field simulation of solidification,Annual Review of Materials Research,32 (2002), 163–194. [3] J. Ni and C. Beckermann, A volume-averaged two-phase model for solidication transport phenomena, Metallurgical Transactions B,22 (1991), 349–361. [4] W. D. Bennon and F. P. Incropera, A Continuum Model for Momentum, Heat and Species Transport in Binary Solid-Liquid Phase Change Systems: I. Model Formulation,Int. J. Heat Mass Transfer,30 (1987), 2161–2170. [5] W. D. Bennon and F. P. Incropera, A Continuum Model for Momentum, Heat and Species Transport in Binary Solid Liquid Phase Change Systems: II. Application to Solidification in a Rectangular Cavity,Int. J. Heat Mass Transfer,30 (1987), 2171–2187. [6] Ph. Blanc, L. Gasser and J. Rappaz, Existence for a stationary model of binary alloy solidification, RAIRO Mod´el. Math. Anal. Num´er.,29 (1995), 687–699. [7] J. L. Boldrini and G. Planas, A tridimensional phase-field model with convection for phase change of an alloy,Discrete Contin. Dyn. Syst.,13 (2005), 429–450. [8] P. Constantin and C. Foias, Navier-Stokes Equations, Chicago Lectures in Mathematics, University of Chicago Press, Chicago, IL, 1988. [9] M. Dur´an, E. Ortega-Torres and J. Rappaz, Weak solution of a stationary convection-diffusion model describing binary alloy solidification processes,Math. Comput. Modelling,42 (2005), 1269–1286. [10] M. Dur´an and E. Ortega-Torres, Existence of weak solutions for a non-homogeneous solidification model, Electron. J. Differential Equations,2006 17 pp. [11] I. Ekeland and R. Temam, Convex analysis and Variational Problems, Society for Industrial and Applied Mathematics, Philadelphia, 1999. [12] R. Glowinski, J.-L. Lions and J. He, Exact and Approximate Controllability for Distributed Parameter Systems: A Numerical Approach, Encyclopedia of Mathematics and its Applications, 117, Cambridge University Press, Cambridge, 2008. [13] O. A. Ladyzenskaya, V. A. Solonnikov and N. N. Uraltzeva, Linear and Quasilinear Equations of Parabolic Type, Translations of Mathematical Monographs, 23, American Mathematical Society, Providence, R.I., 1967. [14] G. Planas and J. L. Boldrini, A bidimensional phase-field model with convection for change phase of an alloy,J. Math. Anal. Appl.,303 (2005), 669–687. [15] J. C. Ramirez and C. Beckermann, Examination of binary alloy free dendritic growth theories with a phase-field model,Acta Materialia,53 (2005), 1721–1736. [16] J.-C. Saut and B. Scheurer, Unique continuation for some evolution equations,J. Differ. Equations,66 (1987), 118–139. [17] J. Simon, Compact sets in the space Lp(0, T ;B),Annali di Matematica Pura ed Applicata, 146 (1986), 64–96. [18] J. Simon, A constructive proof of a theorem of G. de Rham, C. R. Acad. Sci. Paris S´er. I Math.,316 (1993), 1167–1172. [19] R. Temam, Navier Stokes Equations. Theory and Numerical Analysis, Reprint of the 1984 edition, AMS Chelsea Publishing, Providence, RI, 2001. Received October 2012; revised January 2013. E-mail address:[email protected] E-mail address:[email protected] E-mail address:[email protected]