Stability and convergence of two discrete schemes for a degenerate solutal non-isothermal phase-field model

ESAIM: M2AN 43 (2009) 563–589 ESAIM: Mathematical Modelling and Numerical Analysis DOI: 10.1051/m2an/2009011 www.esaim-m2an.org STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES FOR A DEGENERATE SOLUTAL NON-ISOTHERMAL PHASE-FIELD MODEL ∗ Francisco Guill´ en-Gonz´ alez1and Juan Vicente Guti´ errez-Santacreu1 Abstract. We analyze two numerical schemes of Euler type in time and C0finite-element type with P1-approximation in space for solving a phase-field model of a binary alloy with thermal properties. This model is written as a highly non-linear parabolic system with three unknowns: phase-field, solute concentration and temperature, where the diffusion for the temperature and solute concentration may degenerate. The first scheme is nonlinear, unconditionally stable and convergent. The other scheme is linear but conditionally stable and convergent. A maximum principle is avoided in both schemes, using a truncation operator on the L2projection onto the P0finite element for the discrete concentration. In addition, for the model when the heat conductivity and solute diffusion coefficients are constants, optimal error estimates for both schemes are shown based on stability estimates. Mathematics Subject Classification. 35Q72, 35K65, 65M12, 65M60. Received November 15, 2007. Revised May 21, 2008 and October 1st, 2008. Published online April 30, 2009. 1. Introduction 1.1. The model The phase-field method provides a mathematical description for free-boundary problems associated to physical processes with phase transitions. It postulates the existence of a function, called the phase-field, whose value identifies the phase at a particular point in space and time. The method is particularly suitable for cases with complex growth structures occurring during phase transitions. The mathematical model studied in this work describes the solidification process occurring in a binary alloy with temperature-dependent properties. It is based on a highly non-linear parabolic system of partial differential equations with three dependent variables: phase-field, solute concentration and temperature. Moreover, the temperature and concentration equation have nonlinear degenerate diffusivity. Keywords and phrases. Phase-field models, diffuse interface model, solidification process, degenerate parabolic systems, backward Euler schemes, finite elements, stability, convergence, error estimates. ∗This work has been partially supported by DGI-MEC (Spain), Grant MTM2006–07932 and CGCI MECD-DGU Brazil/Spain, Grant 117/06. 1Dpto. E.D.A.N., University of Sevilla, Aptdo. 1160, 41080 Sevilla, Spain. [email protected]; [email protected] Article published by EDP Sciences c EDP Sciences, SMAI 2009 564 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU Let Ω ⊆Rd(d= 2 or 3) be a bounded domain with boundary Γ. Denote by [0,T] the time interval (T>0). We use the notation Q=Ω×(0,T), Σ = Γ ×(0,T)andn(x) is the outwards unit normal vector to Ω at the point x∈Γ. After some physical simplifications [5], we consider the following differential problem, related to a phase-field model of a binary alloy with thermal properties [1]: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ αε2φt−ε2Δφ=1 2(φ−φ3)+β(θ−θAc−θB(1 −c)) in Q, CVθt+l 2φt=∇·[K1(φ)∇θ]inQ, ct=∇·[K2(φ)(∇c+Mc(1 −c)∇φ)] in Q. (1.1) This model is completed with the Neumann boundary conditions ∂φ ∂nΣ=0,(K1(φ)∇θ)·nΣ=0,(K2(φ)∇c)·nΣ=0 (1.2) and the initial conditions φ(x,0) = φ0(x),θ(x,0) = θ0(x),c(x,0) = c0(x)x∈Ω.(1.3) The unknowns for this problem are: φ:Q→R(phase-field) is the state variable characterizing the different phases so that φ= 1 represents the liquid phase and φ=−1 represents the solid phase, θ:Q→Ris the temperature of the material, c:Q→[0,1] (concentration) represents the fraction of one of the two materials in the mixture. The parameter α>0 is the relaxation scaling; the parameter βis given by β=ε[s]/3σ,where ε>0 is the measure of the interface width, σthe surface tension and [s] the entropy density difference between phases; θA,θBare the melting temperatures of each of the two materials in the alloy; CV>0 is the specific heat; l>0 the latent heat; K1≥0 the thermal conductivity; K2≥0 the solute diffusivity; M∈Ris a constant related to the slopes of solid and liquid lines. We will assume that K1=K1(φ)andK2=K2(φ) are two globally Lipschitz continuous functions satisfying 0≤K1(r)≤b1,0≤K2(r)≤b2∀r∈R, with b1,b 2>0. In this sense, the problem is singular with respect to the temperature and concentration when K1(φ)=0orK2(φ) = 0, respectively. As physically the diffusion of material in the solid phase can be considered close to zero [5]; this leads to a degenerate solute diffusion. Such a phenomenon is included in this model, assuming that K2(φ)=0ifφ=−1. On the other hand, although the heat conductivity is nonzero in both solid and liquid phases, we also consider a degenerate diffusion for the temperature with the aim of considering a more general model. Moreover, this may help the development of numerical methods for systems with similar characteristics. The phase-field model for solidification (1.1) is used to treat phenomena such as crystal growth and the fusion of materials. Now we introduce the definition of weak solutions similar to that given in [1,12] which take into account the heat and solute degenerate diffusivity, respectively. Moreover, the maximum principle for the concentration equation says us that 0 ≤c≤1inQif 0 ≤c0≤1inΩ. Definition 1.1. Atriplet(φ, θ, c) is called a weak solution of (1.1)–(1.3)in(0,T)if: (1) φ∈L2(0,T;H2(Ω)) ∩L∞(0,T,H1(Ω)),φ t∈L2(Q),φ(0) = φ0,∂φ ∂n=0a.e. on Σ, (2) θ∈L∞(0,T;L2(Ω)),θ t∈L2(0,T,H1(Ω)),θ(0) = θ0, (3) c∈L∞(0,T;L2(Ω)),c t∈L2(0,T;H1(Ω)),c(0) = c0,0≤c≤1a.e. in Q, (4) J1:= ∇(K1(φ)θ)−θ∇K1(φ)∈L2(Q), (5) J2:= ∇(K2(φ)c)−c∇K2(φ)∈L2(Q), STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 565 verifying αε2φt−ε2Δφ=1 2(φ−φ3)+β(θ−θAc−θB(1 −c)) a.e. in Q, CVT 0 θt,ηdt+l 2T 0φt,ηdt+T 0J1,∇η=0, T 0 ct,ηdt+T 0J2,∇ηdt+MT 0K2(φ)c(1 −c)∇φ, ∇η=0, for each η∈L2(0,T;H1(Ω)). If, in addition, K1,K 2≥b0>0, then θ,c ∈L2(0,T;H1(Ω)) and J1=K1(φ)∇θ, J2=K1(φ)∇c. Here and in what follows, ·,·denotes the inner product in L2(Ω) and ·,· denotes the duality between H1(Ω)and H1(Ω). 1.2. Known results In [1], the existence of weak solutions of problem (1.1)–(1.3) but with a constant solute diffusivity (K2>0) is obtained via the introduction of a regularized problem approximating the degenerate thermal conductivity K1by a strictly positive, regular function followed by the derivation of suitable aprioriestimates and the application of compactness arguments. More concretely, the following existence result was established in [1]. Theorem 1.2. Let Ωbe an open bounded domain of Rd,d=2or 3, with smooth boundary Γ. Assume φ0∈H1+γ(Ω) with 1/2<γ≤1such that ∂φ0 ∂n=0on Γ,θ0∈L2(Ω) and c0∈H1(Ω) such that 0≤c0≤1a.e. in Ω.Then,thereexists(φ, θ, c)aweaksolutionof(1.1)–(1.3)(with K1>0a constant) in (0,T). In addition, in [1] the authors say that the hypothesis φ0∈H1+γ(Ω) with 1/2<γ≤1 is not essential, and the result holds for φ0∈H1(Ω). Scheid [12] proved the existence of weak solutions, by using a similar methodology to [1], for the following isothermal phase-field model of a binary alloy αε2φt−ε2Δφ=F1(φ)+cF2(φ)inQ, ct=∇·[D1(φ)(∇c+D2(c, φ)∇φ)] in Q, (1.4) which has a degenerate solute diffusivity D1(φ)≥0. The main difficulty of model (1.4) is the treatment of the nonlinear term involving D1(φ)D2(c, φ)∇φin the concentration equation. Moreover, the maximum principle for the phase-field variable gives −1≤φ≤1 under the assumptions that the above nonlinearities F1(φ)and F2(φ)vanishwhenφ=−1andφ=1. Error estimates of nonlinear numerical schemes for isothermal phase-field models related to binary alloys are given in [10] for a model as (1.4), and, in [4], considering anisotropic diffusion for the phase-field equation and a more general right-hand side in the phase-field equation, changing the terms F1(φ)+cF2(φ)consideredin [10]byS(c, φ) being a bounded, Lipschitz function. In [7], optimal error estimates are given for a fully discrete nonlinear numerical scheme of a more simplified phase-field model than (1.1) without the concentration (that is one material is only considered) and with constant thermal conductivity K1>0, paying special attention on the dependency of the parameter ε. Stability estimates independent of εare proved for ksmall enough with respect to ε,andα, β are constants depending on ε.Itis also shown some error bounds depending only on a lower polynomial order for 1/ε. Moreover, error estimates are used to establish the convergence of the fully discrete scheme to solutions of the sharp interface limits under different scaling hypotheses in its coefficients. In [2], a time-discrete nonlinear scheme is proposed for a phase-field problem again without the concentration variable and replacing in the equation for the temperature the term l 2φtby the more general term l 2f(θ,φ)t, where fis a generic function satisfying some adequate properties. Convergence of this semi-discrete in time scheme is proved, obtaining the existence and regularity of solutions for the limit problem. 566 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU 1.3. Main results of the paper In this work we will consider two numerical schemes in order to approximate problem (1.1) using continuous P1-finite elements for the tree variables (φ, θ, c). Since a maximum principle cannot be verified in general by the discrete concentration, we introduce a truncation operator on the L2projection onto P0, in order to guarantee aL∞bound for some terms in the discrete concentration equation. A similar idea of truncation, but without the L2projection onto P0, has been used in [8]fora2DNavier-Stokes model with mass diffusion. First of all, we will present in Section 2the nonlinear numerical scheme (2.1)–(2.3) which will be unconditionally stable and convergent. Theorem 1.3 (unconditionally stable, convergent nonlinear scheme).Assume φ0∈H1(Ω),θ 0∈L2(Ω) and c0∈L2(Ω) such that 0≤c0≤1a.e. in Ω. Let Ωbe such that the H2-regularity for the Neumann problem (3.19)holds. Let Thbe a regular, quasi-uniform family of a polyhedral domain Ω. Then, there exists a convergent subsequence of functions φh,k,θ h,k and ch,k associated to scheme (2.1)–(2.3)(see Def. 3.5) towards a weak solution (φ, θ, c)of problem (1.1)–(1.3)in (0,T), as (h, k)→0in the following sense: θh,k →θ, ch,k →c, in L∞(0,T;L2(Ω))-weak∗, φh,k →φ, in L∞(0,T;H1(Ω))-weak∗,and in L2(0,T;H1(Ω))-strong. Second, we construct the linear numerical scheme (6.1)–(6.3) which will be conditionally stable and convergent. Theorem 1.4 (conditionally stable, convergent linear scheme).Assume the assumptions of Theorem 1.3 and the constraint (S) lim (h,k)→0 k h=0. Then, there exists a convergent subsequence of functions φh,k,θh,k and ch,k associated to scheme (6.1)–(6.3) (see Def. 3.5) towards a weak solution (φ, θ, c)of problem (1.1)–(1.3)in (0,T),as(h, k)→0in the same sense of Theorem 1.3. At this point, it is well to point out that, in particular, the two previous theorems provide the existence of weak solutions of problem (1.1)–(1.3) under hypotheses on the data weaker than those imposed in [1](see Thm. 1.2). To be more precise, the hypothesis on c0is relaxed from c0∈H1(Ω) imposed in [1]toc0∈L2(Ω) as was considered in [12]. Recall that in [12] the isothermal case is considered and in [1] there is not degeneration in the solute diffusivity. Finally, assuming both the heat conductivity and solute diffusion coefficients are constants, error estimates of order O(k+h) are shown towards a regular enough continuous solution to (1.1)–(1.3). Theorem 1.5. Under hypotheses of Theorem 1.3 (respectively Thm. 1.4), admitting that K1,K 2are two positive constants and that there exist a continuous solution (φ, θ, c)to (1.1)–(1.3)which verifies the regularity assumptions (7.2), then the discrete solution of scheme (2.1)–(2.3)(respectively (6.1)–(6.3)) satisfies for all n<N en+1 φ2 H1(Ω) +1 4ε2en+1 φ4 L4(Ω) +CV|en+1 θ|2+|en+1 c|2+αk n  l=1  el+1 φ−el φ k 2 +K1k n  l=1 |∇el+1 θ|2+K2k n  l=1 |∇el+1 c|2≤Ch2+k2, STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 567 where Cis a constant independent of hand k, and the errors are denoted by en+1 φ=φn+1 h−φ(tn+1), en+1 θ=θn+1 h−θ(tn+1)and en+1 c=cn+1 h−c(tn+1). The rest of the paper is described as follows. In Section 2the nonlinear scheme (2.1)–(2.3)ispresented, obtaining its unconditionally stability in Section 3. In Section 4some necessary compactness results are proved, passing to the limit in Section 5and concluding the proof of Theorem 1.3. In addition, a conditionally stable and convergent linear scheme is studied in Section 6giving an outline of the proof of Theorem 1.4. Finally, Section 7is devoted to studying optimal error estimates for both schemes. 2. A nonlinear scheme In what follows, let us consider a uniform partition tn=nk of the time interval [0,T]withk=T/N the time step, let Ω ⊂Rd(d= 2 or 3) be a domain with polyhedral boundary and Thbe a family of triangulations of Ω with Ω=  K∈Th K.Hereh:= max K∈Th hKwith hKthe diameter of K.LetXhbe the finite element subspace of H1(Ω) furnished by globally continuous, piecewise linear functions, that is, Xh={xh∈C0(Ω): xh|K∈P1(K),∀K∈T h,}· A first idea to approximate equation (1.1)1is ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ αφn+1 h−φn h k,x h+∇φn+1 h,∇xh+1 2ε2(φn+1 h)3,x h =1 2ε2φn h,x h+β ε2θn h−θAcn h−θB(1 −cn h),x h,∀xh∈Xh. But, we add φn+1 h,x hto the left-hand side and φn h,x hto the right-hand side which will cancel each other in the limit as (h, k) go to zero. The reason why we introduce these terms is to get stability constants only of polynomial order with respect to εavoiding exponential dependence. Concretely, since β=O(ε), we will get stability constants depending on 1/ε (see Rem. 3.2 below). Then we propose the following scheme to approximate problem (1.1)–(1.3): Initialization:Let(φ0 h,θ0 h,c 0 h)∈Xh×Xh×Xhbe suitable approximations of (φ0,θ 0,c 0). Step n+1: Given(φn h,θn h,c n h)∈Xh×Xh×Xh. Find φn+1 h∈Xhas a solution of the problem: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ αφn+1 h−φn h k,x h+∇φn+1 h,∇xh+φn+1 h,x h+1 2ε2(φn+1 h)3,x h =1 2ε2+1 φn h,x h+β ε2θn h−θAcn h−θB(1 −cn h),x h,∀xh∈Xh. (2.1) Find θn+1 h∈Xhand cn+1 h∈Xhas solutions of the decoupled variational problems: CVθn+1 h−θn h k,x h+Kh 1(φn+1 h)∇θn+1 h,∇xh=−l 2φn+1 h−φn h k,x h,∀xh∈Xh,(2.2) ⎧ ⎪ ⎨ ⎪ ⎩cn+1 h−cn h k,x h+Kh 2(φn+1 h)∇cn+1 h,∇xh =−MKh 2(φn+1 h)[P0cn h]T(1 −[P0cn h]T)∇φn h,∇xh,∀xh∈Xh. (2.3) 568 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU Here Kh 1=K1+g1(h), Kh 2=K2+g2(h), where giare positive functions to be chosen later, and P0is the L2orthogonal projector onto X0 h,whereX0 his the finite element space of piecewise constant functions, and [·]T is a truncation operator defined from X0 hinto X0 has follows: Given xh∈X0 h,then[xh]T∈X0 hsuch that ∀K∈T h,[xh]T|K=⎧ ⎨ ⎩ xh|Kif xh|K∈[0,1], 0ifxh|K<0, 1ifxh|K>1. Since (2.2)and(2.3) are quadratic linear systems, it is easy to check the existence and uniqueness of solutions. On the other hand, (2.1) is a discrete nonlinear variational problem and its existence and uniqueness can be proved as follows: We define J(φh)= α 2kΩ |φh|2+1 2Ω|∇φh|2+|φh|2+1 8ε2Ω |φh|4−Ω gφ h,(2.4) where g=α kφn h+1 2ε2+1 φn h+β ε2θn h−θAcn h−θB(1 −cn h). Clearly, Jis a strictly convex functional on Xh, then the minimum problem min φh∈Xh J(φh) has a unique solution characterized by its Euler equation (2.1). We will denote by Cgeneric positive constants always independent of the discretization parameters hand k. 3. APRIORI estimates and weak convergences Let us add and subtract the term 1 2ε2φn+1 h,x hto the left-hand side of (2.1) in order to rewrite (2.1) with respect to the so-called Ginzburg-Landau function f(φ)= 1 2ε2φ2−1φwhich has the potential function F(φ)= 1 8ε2φ2−12,thatis,f(φ)=∇φF(φ). Then, (2.1) is rewritten as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αφn+1 h−φn h k,x h+∇φn+1 h,∇xh+φn+1 h,x h +f(φn+1 h),x h+1 2ε2φn+1 h−φn h,x h=φn h,x h +β ε2θn h−θAcn h−θB(1 −cn h),x h∀xh∈Xh. (3.1) It is easy to check that if we select φ0 h=Ihφ0,θ0 h=Ihθ0and c0 h=Ihc0as initial approximations, where Ih is an interpolation operator into Xhsatisfying stability properties in the L2,L4and H1norms, it follows that there exists a constant C2>0 (independent of ε) such that 2CVβ lε2|θ0 h|2+|c0 h|2+1+φ0 h2 H1(Ω) +1 4ε2Ω (|φ0 h|2−1)2≤C2 ε2·(3.2) For instance, this is true when Ihis the L2-projector onto Xh,orIhis the Cl´ement or Scott-Zhang regularization operator. Let us denote by |·|the L2(Ω)-norm and by · H1(Ω) the H1(Ω)-norm. With such a notation we establish the following stability result. STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 569 Lemma 3.1. Assume φ0∈H1(Ω),θ0∈L2(Ω) and c0∈L2(Ω) such that 0≤c0≤1a.e. in Ω. Then, for each k such that βk ε2is sufficiently small, the discrete solution of scheme (2.1)–(2.3)satisfies the following estimates: (i)max 0≤n≤Nφn h2 H1(Ω) ≤C, (ii) N−1  n=0 φn+1 h−φn h2 H1(Ω) ≤C, (iii)k N−1  n=0  φn+1 h−φn h k 2 ≤C, (iv)max 0≤n≤N|θn h|2≤C, (v) N−1  n=0 |θn+1 h−θn h|2≤C, (vi)k N−1  n=0 |Kh 1(φn+1 h)∇θn+1 h|2≤C, (vii)max 0≤n≤N|cn h|2≤C, (viii) N−1  n=0 |cn+1 h−cn h|2≤C, (ix)k N−1  n=0 |Kh 2(φn+1 h)∇cn+1 h|2≤C, where C>0depends on εand the data (φ0,θ 0,c 0)but is independent of (h, k). Proof. Let xh=4β lε 2kθn+1 hand xh=2kc n+1 hbe test functions in (2.2)and(2.3), respectively. Now by using the identity (a−b, 2a)=|a|2−|b|2+|a−b|2and bounding adequately the right-hand side, we have 2CVβ lε2|θn+1 h|2−|θn h|2+|θn+1 h−θn h|2+4β lε2k|Kh 1(φn+1 h)∇θn+1 h|2 =−2β ε2kφn+1 h−φn h k,θn h+(θn+1 h−θn h) ≤− 2β ε2kφn+1 h−φn h k,θn h+α 2k φn+1 h−φn h k 2 +2β2 αε4k|θn+1 h−θn h|2, (3.3) |cn+1 h|2−|cn h|2+|cn+1 h−cn h|2+k|Kh 2(φn+1 h)∇cn+1 h|2≤Ck|∇φn h|2.(3.4) By choosing βk ε2sufficiently small to control the last term on the right-hand side of (3.3), this inequality reduces to 2CVβ lε2|θn+1 h|2−|θn h|2+1 2|θn+1 h−θn h|2+4β lε2k|Kh 1(φn+1 h)∇θn+1 h|2≤ −2β ε2kφn+1 h−φn h k,θn h+α 2k φn+1 h−φn h k 2 ·(3.5) Next, take xh=2kφn+1 h−φn h kas a test function in (2.1), it follows that αk φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω)+2 f(φn+1 h),φ n+1 h−φn h +1 ε2|φn+1 h−φn h|2≤2β ε2kθn h,φn+1 h−φn h k+Ck|φn h|2+Cβ2 αε4k(|cn h|2+1),(3.6) where Cis a constant independent of ε. 570 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU Now, using again the identity (a−b, a)=1 2(|a|2−|b|2+|a−b|2) twice, we can rewrite 2f(φn+1 h),φ n+1 h−φn h=1 2ε2Ω(φn+1 h)2−1(φn+1 h)2−(φn h)2+(φn+1 h−φn h)2 =2ΩF(φn+1 h)−F(φn h)+ 1 8ε2((φn+1 h)2−(φn h)2)2 +1 2ε2Ω(φn+1 h)2−1|φn+1 h−φn h|2.(3.7) Note that the negative term on the right-hand side of (3.7) can be absorbed by the last term on the left-hand side of (3.6). This property can be summarized as 1 2ε2(φn+1 h)3−φn hφn+1 h−φn h≥F(φn+1 h)−F(φn h) (3.8) which is an appropriate discrete version of the equality f(φ)φt=F(φ)t. Next, if we add up (3.4), (3.5)and(3.6), the first term on the right hand side of (3.5)and(3.6) disappears, and we get 2CVβ lε2|θn+1 h|2−|θn h|2+1 2|θn+1 h−θn h|2+4β lε2k|Kh 1(φn+1 h)∇θn+1 h|2 +(|cn+1 h|2+1)−(|cn h|2+1)+|cn+1 h−cn h|2+k|Kh 2(φn+1 h)∇cn+1 h|2 +α 2k φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω) +2ΩF(φn+1 h)−F(φn h)≤Ckφn h2 H1(Ω) +Cβ2 αε4k(|cn h|2+1). Finally, by summing over nthe discrete Gronwall lemma and the initial bound (3.2) provide the desired estimates, and this completes the proof.  Remark 3.2. Since the initial estimates have order O(ε−2), see (3.2), then the stability estimates obtained in Lemma 3.1 are of order O(ε−2eCβ2ε−4) for the variables (φ, √β εθ,c). As βis of order O(ε), the order reduces to O(ε−2eCε−2). In particular, these estimates would be independent of εif βwere of order O(ε2) and considering an initial bound (3.2) independent of ε. Furthermore, if we truncate the discrete concentration cn hin (2.1)as made in (2.3), that is replacing β ε2θn h−θAcn h−θB(1 −cn h),x hby β ε2θn h−θA[P0cn h]T−θB(1 −[P0cn h]T),x h, this modified scheme has stability estimates of order O(ε−2+β2ε−4). Remark 3.3. Observe that in this nonlinear scheme, we have used a first-order semi-implicit approximation of the Ginzburg-Landau function f(φ), which provides a stationary problem to solve in each time step, identified with the critical point of a convex functional (see (2.4)). Moreover, this approximation verifies the property (3.8). For instance, if we use the first-order implicit approximation f(φn+1 h), then the associated stationary problem (of Allen-Cahn type) is related to the critical points of a non-convex functional and the property (3.8)innot verified, because a negative term appears on the right-hand side. To be more concrete, it follows that f(φn+1 h)φn+1 h−φn h≥F(φn+1 h)−F(φn h)−1 4ε2(φn+1 h−φn h)2. STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 571 Consider the linear operator Lh:Xh→Xhdefined as: Lhφh,x h=∇φh,∇xh+φh,x h∀xh∈Xh.(3.9) Then, the discrete phase-field equation (2.1) can be rewritten as: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩φn+1 h−φn h k,x h+1 αLhφn+1 h,x h+1 2αε2(φn+1 h)3,x h =1 2αε2+1 αφn h,x h+β αε2θn h−θAcn h−θB(1 −cn h),x h,∀xh∈Xh. (3.10) Taking xh=Lhφn+1 has a test function in (3.10) and using the estimates of Lemma 3.1, the following result can be established. Corollary 3.4. Under the hypotheses of Lemma 3.1, it holds k N−1  n=0 |Lhφn+1 h|2≤C. On the other hand, since K1(·)≤b1and K2(·)≤b2,from(vi)and(ix) of Lemma 3.1 we also have k N−1  n=0 |Kh 1(φn+1 h)∇θn+1 h|2≤C, k N−1  n=0 |Kh 2(φn+1 h)∇cn+1 h|2≤C. Definition 3.5. We define φh,k (respectively  φh,k) as the piecewise constant functions in time taking values φn+1 hon (tn,t n+1] (respectively φn h). Analogously, we define θh,k, θh,k,andch,k,ch,k. Moreover, we define  φh,k,  θh,k,ch,k ∈C0([0,T]; Xh) as the piecewise linear functions in time such that  φh,k(tn)=φn h, θh,k(tn)=θn h, ch,k(tn)=cn h, respectively. An easy consequence of the previous definition, Lemma 3.1 and Corollary 3.4 is the following result. Lemma 3.6. Under the hypotheses of Lemma 3.1, the following estimates hold: {θh,k}h,k,{ θh,k}h,k,{ θh,k}h,k is bounded in L∞(0,T;L2(Ω)),(3.11) {ch,k}h,k,{ch,k}h,k,{ch,k}h,k is bounded in L∞(0,T;L2(Ω)),(3.12) {φh,k}h,k,{ φh,k}h,k,{ φh,k}h,k is bounded in L∞(0,T;H1(Ω)),(3.13) {Kh 1(φh,k)∇θh,k}h,k is bounded in L2(0,T;L2(Ω)),(3.14) {Kh 2(φh,k)∇ch,k}h,k is bounded in L2(0,T;L2(Ω)),(3.15) d dt φh,kh,k is bounded in L2(0,T;L2(Ω)).(3.16) {Lhφh,k}h,k is bounded in L2(0,T;L2(Ω)).(3.17) 578 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU 5. Passing to the limit In order to pass to the limit in the discrete concentration equation, we will use the following result, which is easy to prove because equation (5.1) satisfies the maximum principle: Lemma 5.1. The following two systems are equivalent: ct=∇·(K2(φ)[∇c+MT1 0c(1 −T1 0c)∇φ]) in Q, (5.1) and 0≤c≤1,c t=∇·(K2(φ)[∇c+Mc(1 −c)∇φ]) in Q. To pass to the limit in scheme (2.1)–(2.3), we rewrite the scheme as follows: Taking xh=ηn+1 h∈Xha suitable approximation at time tn+1 of any function η∈C0([0,T]; C∞ c(Ω)) such that η(T) = 0 (clearly ηN h=0) as a test function in (2.1), (2.2)and(2.3), multiplying by k, summing over nand denoting the function ηh,k similarly to Definition 3.5, one arrives at ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ αε2T 0d dt φh,k,η h,k+ε2T 0∇φh,k,∇ηh,k+ε2T 0φh,k − φh,k,η h,k +1 2T 0(φh,k)3− φh,k,η h,k−βT 0 θh,k −θAch,k −θB(1 −ch,k),η h,k=0, CVT 0d dt θh,k,η h,k+l 2T 0d dt φh,k,η h,k+T 0Kh 1(φh,k)∇θh,k,∇ηh,k=0, ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ T 0d dtch,k,η h,k+T 0Kh 2(φh,k)∇ch,k,∇ηh,k +MT 0Kh 2(φh,k)[P0ch,k]T(1 −[P0ch,k]T)∇ φh,k,∇ηh,k=0. (5.2) By applying all the convergences already obtained, there are no additional difficulties in passing to the limit obtaining that (φ, θ, c) is a weak solution of (1.1). In particular, taking (h, k)→0 in the discrete equation for the concentration cand using (4.7), we arrive at the limit equation (5.1); hence 0 ≤c≤1andT1 0c=c. Finally, the discrete phase-field equation is verified pointwise in Qthanks to the strong regularity of φ. The proof of Theorem 1.3 is finished. Remark 5.2. AsmentionedinRemark3.2,taking[P0cn h] instead of cn hin the discrete equation for the phase field variable (2.1) provides better stability estimates with respect to ε. Nevertheless we find that there is a limit function ϑ∈L∞(0,T;L2(Ω)) such that [P0cn h] tends to ϑweakly* in L∞(0,T;L2(Ω)), but it is not clear how to identify ϑwith the limit function c.Byusing(4.14)andthat0≤c≤1a.e. Q (owing to the limit in (5.2)can be taken as before), we can only deduce that ϑ=ca.e. in  Q, Qbeing defined in the proof of Proposition 4.3. 6. A conditionally stable, convergent linear scheme In this section we study a more explicit scheme, where the nonlinear discrete approximation (2.1)of(1.1)1 is considered completely in the previous step time, resulting a linear (and decoupled) scheme. Contrary to the previous nonlinear scheme, now to obtain stability we will impose a constraint on the discrete parameters. Recall the definition of the Ginzburg-Landau function f(φ)= 1 2ε2(φ2−1)φassociated to the potential function F(φ)= 1 8ε2(φ2−1)2. We propose the following linear scheme: Initialization:Let(φ0 h,θ0 h,c 0 h)∈Xh×Xh×Xhbe suitable approximations of (φ0,θ 0,c 0). Step n+1: Given(φn h,θn h,c n h)∈Xh×Xh×Xh. STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 579 Find φn+1 h∈Xhas a solution of the problem: αφn+1 h−φn h k,x h+∇φn+1 h,∇xh+φn+1 h,x h=−f(φn h),x h+φn h,x h +β ε2θn h−θAcn h−θB(1 −cn h),x h,∀xh∈Xh.(6.1) Find θn+1 h∈Xhand cn+1 h∈Xhas solutions of the decoupled variational problems: CVθn+1 h−θn h k,x h+Kh 1(φn+1 h)∇θn+1 h,∇xh=−l 2φn+1 h−φn h k,x h,∀xh∈Xh,(6.2) ⎧ ⎪ ⎨ ⎪ ⎩cn+1 h−cn h k,x h+Kh 2(φn+1 h)∇cn+1 h,∇xh =−MKh 2(φn+1 h)[P0cn h]T(1 −[P0cn h]T)∇φn h,∇xh,∀xh∈Xh. (6.3) The conditional stability of scheme (6.1)–(6.3) will be obtained by induction on the time step n. First of all, we establish the following result which provides a basic recursive inequality. Lemma 6.1. Assume the constraint: (S) lim (h,k)→0k/h =0. If there exists a constant Cd>0(independent of h,kand n,butdependentonε) such that φn h2 H1(Ω) +2CVβ lε2|θn h|2+|cn h|2+1≤Cd,(6.4) then the following inequalities hold for (h, k)sufficiently small (independent of n), ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω) +2CVβ lε2|θn+1 h|2−|θn h|2+1 2|θn+1 h−θn h|2 +2ΩF(φn+1 h)−F(φn h)+ 1 8ε2((φn+1 h)2−(φn h)2)2 +α 2k φn+1 h−φn h k 2 +2β lε2k|Kh 1(φn+1 h)∇θn+1 h|2 ≤R1k|φn h|2+R2 β2 ε4k(|cn h|2+1), (6.5) |cn+1 h|2−|cn h|2+|cn+1 h−cn h|2+k|Kh 2(φn+1 h)∇cn+1 h|2≤R1k|∇φn h|2,(6.6) where R1and R2are positive constants independent of h,k,n,andε. Proof. Firstly, we consider xh=2kφn+1 h−φn h kin (6.1) and bound the right-hand side 3 2αk  φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω) +2 f(φn h),φ n+1 h−φn h≤2β ε2kθn h,φn+1 h−φn h k+C αk|φn h|2+Cβ2 αε4k(|cn h|2+1).(6.7) 580 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU Now, we handle the last term on the left-hand side of (6.7) as follows: 2φn+1 h−φn h,f(φn h)=1 ε2φn+1 h−φn h,((φn+1 h)2−1)φn h+1 ε2φn+1 h−φn h,((φn h)2−(φn+1 h)2)φn h:= I1−I2. Next,wecontinuerewritingI1as follows: I1=1 2ε2Ω ((φn+1 h)2−1)((φn+1 h)2−(φn h)2−(φn+1 h−φn h)2) =1 4ε2Ω((φn+1 h)2−1)2−((φn h)2−1)2+((φn+1 h)2−(φn h)2)2 +k2 2ε2Ω (1 −(φn+1 h)2) φn+1 h−φn h k 2 ·(6.8) The term I2is bounded as I2≤C1 ε2k2φn h2 L∞(Ω)  φn+1 h−φn h k 2 +1 8ε2Ω ((φn+1 h)2−(φn h)2)2. Therefore, we get from (6.7) and the previous computations 3 2αk φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω) +2ΩF(φn+1 h)−F(φn h)+ 1 8ε2((φn+1 h)2−(φn h)2)2+k2 2ε2 φn+1 h−φn h k 2 ≤2β ε2kθn h,φn+1 h−φn h k+C αk|φn h|2+Cβ2 ε4αk(|cn h|2+1) +C1 ε2αkφn h2 L∞(Ω) +φn+1 h2 L∞(Ω)k φn+1 h−φn h k 2 ≤2β ε2kθn h,φn+1 h−φn h k+C αk|φn h|2+Cβ2 ε4k(|cn h|2+1) +C1 ε2 k hφn h2 H1(Ω) +φn+1 h2 H1(Ω)k φn+1 h−φn h k 2 , (6.9) where in the last line the inverse estimate xhL∞(Ω) ≤Ch −1/2xhH1(Ω) has been used. Now we are looking for the bound φn+1 hH1(Ω) ≤C1where C1>0 depends on the constant Cdof hypothesis (6.4) but it will be independent of n. It will be carried out by bounding φn+1 hH1(Ω) in terms of φn hH1(Ω), |θn h|and |cn h|and using hypothesis (6.4). Indeed, taking again xh=2kφn+1 h−φn h kas a test function in (6.1), STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 581 but now bounding directly the term depending on f(φn h) on the right-hand side, we get αk φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω)≤Cβ2 ε4αk|θn h|2 +C αk|φn h|2+Cβ2 ε4αk(|cn h|2+1)+Ck|f(φn h)|2 ≤Cβ2 ε4αk|θn h|2+C αk|φn h|2+Cβ2 ε4αk(|cn h|2+1) +C1 ε4kφn h6 H1(Ω) +C1 ε4kφn h2 H1(Ω). In particular, by using hypothesis (6.4), the previous inequality says us φn+1 h2 H1(Ω) ≤φn h2 H1(Ω) +Ck Cdβ ε2+Cd+Cdβ2 ε4+C3 d ε4+Cd ε4≤φn h2 H1(Ω) +C1(ε)k with C1(ε) independent of h,kand n. Thus, by using the previous estimate in (6.9) and again hypothesis (6.4), we get 3 2αk φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω) +2ΩF(φn+1 h)−F(φn h)+ 1 8ε2((φn+1 h)2−(φn h)2)2 ≤2β ε2kθn h,φn+1 h−φn h k+C αk|φn h|2+Cβ2 αε4k(|cn h|2+1)+Ck h 1 ε2Cd+C1(ε)kk φn+1 h−φn h k 2 · By taking into account the constraint (S), in particular lim (h,k)→0 k h 1 ε2Cd+C1(ε)k= 0, so for any (h, k)small enough such that Ck h 1 ε2Cd+C1(ε)k≤1 2α, the last term on the right-hand side can be absorbed, and remains αk φn+1 h−φn h k 2 +φn+1 h2 H1(Ω) −φn h2 H1(Ω) +φn+1 h−φn h2 H1(Ω)+2ΩF(φn+1 h)−F(φn h) +1 8ε2((φn+1 h)2−(φn h)2)2≤2β ε2kθn h,φn+1 h−φn h k+C αk|φn h|2+Cβ2 αε4k(|cn h|2+1).(6.10) On the other hand, take xh=4β lε2kθn+1 hin (6.2) to arrive at inequality (3.5), that is 2CVβ lε2|θn+1 h|2−|θn h|2+1 2|θn+1 h−θn h|2+4β lε2k|Kh 1(φn+1 h)∇θn+1 h|2≤ −2β ε2kφn+1 h−φn h k,θn h+α 2k φn+1 h−φn h k 2 ·(6.11) Consequently, it suffices to add up (6.11)and(6.10)toget(6.5). Finally, inequality (6.6) is easily obtained by testing (6.3)bycn+1 hand bounding adequately as in the proof of Lemma 3.1. 582 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU On the other hand, we turn our attention to the initial bound (3.2) which in particular verifies hypothesis (6.4) imposed in Lemma 6.1. It is very important in order to guarantee a correct induction argument. Now, we are in position to give the following stability result. Lemma 6.2. Under the hypotheses of Lemma 6.1, the discrete solution of scheme (6.1)–(6.3)satisfies the following estimates: (i)max 0≤n≤Nφn h2 H1(Ω) ≤C, (ii) N−1  n=0 φn+1 h−φn h2 H1(Ω) ≤C, (iii)k N−1  n=0  φn+1 h−φn h k 2 ≤C, (iv)max 0≤n≤N|θn h|2≤C, (v) N−1  n=0 |θn+1 h−θn h|2≤C, (vi)k N−1  n=0 |Kh 1(φn+1 h)∇θn+1 h|2≤C, (vii)max 0≤n≤N|cn h|2≤C, (viii) N−1  n=0 |cn+1 h−cn h|2≤C, (ix)k N−1  n=0 |Kh 2(φn+1 h)∇cn+1 h|2≤C, where C>0is independent of (h, k)and depends on the data (φ0,θ 0,c 0),α, β and ε. Proof. Obviously, if we let (6.5)and(6.6) hold for n=0, ..., N −1, we get all the statements of this lemma by adding (6.5)and(6.6) and applying the discrete Gronwall lemma. Therefore, it suffices to prove that (6.5) and (6.6) hold for n=0, ..., N −1. Let us consider Cd=e (R1+β2 ε4R2)TC2/ε2with C2>0 given in (3.2)andR1,R 2given in Lemma 6.1.Asthe initial approximations hold hypothesis (6.4)forn= 0, inequalities (6.5)and(6.6) are satisfied for n=0. The final induction step can be easily seen by assuming that inequalities (6.5)and(6.6) hold for l=0, ..., n−1. Then, adding up (6.5)and(6.6)from0ton−1, one has 2CVβ lε2|θn h|2+|cn h|2+1+φn h2 H1(Ω) +2Ω F(φn h)≤2CVβ lε2|θ0 h|2+|c0 h|2+1+φ0 h2 H1(Ω) +2Ω F(φ0 h)+ k n−1  l=0 R1φl h2 H1+R2 β2 ε4(|cl h|2+1) . Now, the discrete Gronwall lemma and (3.2) yield 2CVβ lε2|θn h|2+|cn h|2+1+φn h2 H1(Ω) +2Ω F(φn h)≤e(R1+β2 ε4R2)(n−1) k2CVβ lε2|θ0 h|2+|c0 h|2 +1+φ0 h2 H1(Ω) +2Ω F(φ0 h)≤e(R1+β2 ε4R2)TC2/ε2:= Cd. Then, we find that hypothesis (6.4) is satisfied. Therefore, in view of Lemma 6.1, inequalities (6.5)and(6.6) hold.  Note that the stability estimates obtained in Lemma 6.2 are of order O((1/ε2)e(β2/ε4)) for the variable (φ, √β εθ,c) as in the nonlinear scheme, but now an adequate constraint for (h, k) small enough (depending exponentially on 1/ε) is necessary (recall that in the nonlinear scheme, only βk ε2small enough was imposed). To finish the proof of Theorem 1.4 it is necessary to prove the convergence of the linear scheme (6.1)–(6.3). But, as the argument for this is similar to that developed for the nonlinear scheme (2.1)–(2.3), it is left to the reader.  STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 583 7. Error estimates for the non-degenerate case In this section we deal with the error analysis of both linear and nonlinear scheme. The presence of the truncation operator applied to the piecewise constant operator P0makes nonstandard this error analysis and this particular truncation is responsible of order O(h) in error estimates, although higher-order finite elements were considered. In order to be able to guarantee a sufficient regular solution of problem (1.1)–(1.3) we assume the nondegenerate case. For simplicity, we assume that K1and K2are positive constants, providing in particular standard Neumann boundary conditions in (1.2). Let {Th},0<h≤1, be a regular, quasi-uniform family of subdivisions of a polyhedral domain Ω ⊂Rm, m= 2 or 3, whose boundary Γ is such that the problem −Δu+u=fin Ω,∂u ∂n= 0 on Γ (7.1) holds the stability property uH2(Ω) ≤|f|,foreachf∈L2(Ω). Recall that, both previous hypotheses are also assumed in Theorem 1.3 to prove the convergence. Define the global error en φ=φn h−φ(tn), en θ=θn h−θ(tn), and en c=cn h−c(tn). These errors are decomposed into a discrete error e·,d, and an interpolation error e·,i as follows en φ,d =φn h−P1 hφ(tn),e n φ,i =P1 hφ(tn)−φ(tn), en θ,d =θn h−P0 hθ(tn),e n θ,i =P0 hθ(tn)−θ(tn), en c,d =cn h−P0 hc(tn),e n c,i =P0 hc(tn)−c(tn), where P1 h:H1→Xhis the H1-projection operator defined as ψ−P1 hψ,xh+∇ψ−∇P1 hψ,∇xh=0 ∀xh∈Xh, and P0 h:L2→Xhis the L2-projection operator defined as ψ−P0 hψ,xh=0 ∀xh∈Xh. Finally, let us recall some approximation properties of P1 hand P0 hto be used later on (see [6], Prop. 1.134, p. 73): ψ−P0 hψH1(Ω) ≤ChψH2(Ω),ψ−P1 hψH1(Ω) ≤ChψH2(Ω), |ψ−P0ψ|≤ChψH1(Ω),|ψ−P1 hψ|≤ChψH1(Ω). In particular, from the last inequality, one has k|δtψ(tn+1)−δtP1 hψ(tn+1)|2≤Ch2tn+1 tn ψt2 H1(Ω), whereweuseδtto denote the discrete backward Euler time derivative, that is δtψ(tn+1)=ψ(tn+1)−ψ(tn) k· Note that, for the H1-interpolation error of P0 h, a quasi-uniform family of finite elements must be assumed and, for the L2-interpolation error of P1 h, a duality argument is required where the elliptic H2-regularity for the Elliptic-Neumann problem (7.1)isimposed. 584 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU Throughout the section we assume a regular solution of (1.1)–(1.3). Concretely, one assumes φ∈L∞(0,T;W1,∞(Ω)) ∩L2(0,T;H2(Ω)),φ t∈L2(0,T;H1(Ω),φ tt ∈L2(0,T;L2(Ω)), θ∈L∞(0,T;H1(Ω)) ∩L2(0,T,H2(Ω)),θ t∈L2(0,T;L2(Ω)),θ tt ∈L2(0,T;H1(Ω)), c∈L∞(0,T;H1(Ω)) ∩L2(0,T,H2(Ω)),c t∈L2(0,T;L2(Ω)),c tt ∈L2(0,T;H1(Ω)). (7.2) 7.1. Errorestimatesforthenonlinearscheme We now state our error estimates for the fully discrete nonlinear scheme (2.1)–(2.3). If we compare the exact problem with the scheme and use the equality a3−b3=(a−b)3+3ab(a−b), then the error equations are given by ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αδten+1 φ,d ,x h+∇en+1 φ,d ,∇xh+en+1 φ,d ,x h+1 2ε2(en+1 φ,d )3,x h =−αδten+1 φ,i ,x h−1 2ε2(en+1 φ,i )3,x h−3 2ε2P1 h(φ(tn+1))φ(tn+1)en+1 φ,i ,x h −3 2ε2φn+1 hP1 h(φ(tn+1))en+1 φ,d ,x h+1 2ε2en φ,d +en φ,i,x h +en φ,d +en φ,i,x h+β ε2en θ,d,x h−β ε2(θA+θB)en c,d,x h+Rn+1 φ,x h, (7.3) ⎧ ⎪ ⎨ ⎪ ⎩ CVδten+1 θ,d ,x h+K1∇en+1 θ,d ,∇xh =−K1∇en+1 θ,i ,∇xh−l 2δten+1 φ,x h+Rn+1 θ,x h, (7.4) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δten+1 c,d ,x h+K2∇en+1 c,d ,∇xh=−K2∇en+1 c,i ,∇xh −MK 2[P0cn h]T(1 −[P0cn h]T)∇en φ,∇xh −MK 2(1 −[P0cn h]T)en cT∇φ(tn),∇xh +MK 2c(tn)en cT∇φ(tn),∇xh+Rn+1 c,x h, (7.5) where δten+1 =(en+1 −en)/k,en cT=[P0cn h]T−c(tn)and Rn+1 φ=α ktn+1 tn (tn−t)φtt(t)dt+1 2ε2+1 tn+1 tn φt(t)dt +β ε2tn+1 tn θt(t)dt+β ε2(θA+θB)tn+1 tn ct(t)dt, (7.6) Rn+1 θ=CV ktn tn (tn−t)θtt(t)dt−l 2ktn+1 tn (tn−t)φtt(t)dt, Rn+1 c=1 ktn+1 tn (tn−t)ctt(t)dt+∇·1−c(tn+1)c(tn+1)tn+1 tn ∇φt(t) +1−(c(tn)+c(tn+1)tn+1 tn ct(t)∇φ(tn). Theorem 7.1. Under the assumptions of Theorem 1.3, if the solution (φ, θ, c)of (1.1)–(1.3)satisfies that (φ(t),θ(t),c(t)) ∈H2(Ω)3for all t∈[0,T]and the regularity given in (7.2), then the following error estimates hold for ksmall enough: max 0≤n≤Nen+1 φ2 H1(Ω) +|en+1 θ|2+|en+1 c|2+k N−1  n=0 |δten+1 φ|2+|∇en+1 θ|2+|∇en+1 c|2≤Ch2+k2,(7.7) STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 585 where the constant C>0depends on the exact solution, but is independent of (h, k). Proof. Set xh=2kδten+1 φ,d ∈Xhas a test function in (7.3), and using the equalities a3(a−b)=1 2a2(a2−b2+(a−b)2)=1 4(a4−b4+(a2−b2)2)+1 2a2(a−b)2 for a=en+1 φand b=en φ,weget ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ α2kδten+1 φ,d  2+en+1 φ,d 2 H1(Ω) −en φ,d2 H1(Ω) +en+1 φ,d −en φ,d2 H1(Ω) +1 4ε2en+1 φ,d 4 L4(Ω) −en φ,d4 L4(Ω) +|(en+1 φ,d )2−(en φ,d)2|2+2|(en+1 φ,d )(en+1 φ,d −en φ,d)|2 =−α2kδten+1 φ,i ,δ ten+1 φ,d −1 ε2k(en+1 φ,i )3,δ ten+1 φ,d  −3 ε2kP1 h(φ(tn+1))φ(tn+1)en+1 φ,i ,δ ten+1 φ,d −3 ε2kφn+1 hP1 h(φ(tn+1))en+1 φ,d ,δ ten+1 φ,d  +1 ε2ken φ,d +en φ,i,δ ten+1 φ,d +2ken φ,d +en φ,i,δ ten+1 φ,d  +2β ε2ken θ,d,δ ten+1 φ,d −2β ε2(θA+θB)ken c,d,δ ten+1 φ,d +2kRn+1 φ,δ ten+1 φ,d := 9  i=1 Ii. (7.8) Now, we must bound each term on the right-hand side of (7.8). We just focus on the terms I2,I3and I4: I2≤Cλ1 ε4αken+1 φ,i 6 L6(Ω) +λkα|δten+1 φ,d |2≤Cλ1 ε4αkh 6φ(tn+1)6 H2(Ω) +λkα|δten+1 φ,d |2, I3≤3 ε2kP1 h(φ(tn+1))L6(Ω)φ(tn+1)L6(Ω)en+1 φ,i L6(Ω)|δten+1 φ,d | ≤Cλ 1 ε4αkφ(tn+1)4 H1(Ω)en+1 φ,i 2 H1(Ω) +λkα|δten+1 φ,d |2 ≤Cλ 1 ε4αkh 2φ(tn+1)4 H1(Ω)φ(tn+1)2 H2(Ω) +λkα|δten+1 φ,d |2, I4≤3 ε2kφn+1 hL6(Ω)P1 h(φ(tn+1))L6(Ω)en+1 φ,d L6(Ω)|δten+1 φ,d | ≤Cλ 1 ε4αkφ(tn+1)2 H1(Ω)en+1 φ,d 2 H1(Ω) +λkα|δten+1 φ,d |2. In the last line, the stability estimate φn+1 hL6(Ω) ≤CgiveninLemma3.1 has been applied. Note that the previous inequalities hold for any λ>0andCλ>0 are different constants of order O(1/λ). Finally, the constants Cλbounding I2and I3are independent of ε, but the constant Cλrelated to I4depends on εvia the bound of Lemma 3.1. The remainder of the terms can be bounded more easily. Thus, by choosing λsmall enough, we get ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αkδten+1 φ,d  2+en+1 φ,d 2 H1(Ω) −en φ,d2 H1(Ω) +en+1 φ,d −en φ,d2 H1(Ω) +1 4ε2en+1 φ,d 4 L4(Ω) −en φ,d4 L4(Ω) +|(en+1 φ,d )2−(en φ,d)2|2+2|en+1 φ,d (en+1 φ,d −en φ,d)|2 ≤Ckh 2h4φ(tn+1)6 H2(Ω) +φ(tn+1)4 H1(Ω)φ(tn+1)2 H2(Ω) +φ(tn)2 H1(Ω) +Ck 1 ε4αβ2|en c,d|2+β2CV|en θ,d|2+en φ,d2 H1(Ω) +φ(tn+1)2 H1(Ω)en+1 φ,d 2 H1(Ω) +Ch 2tn+1 tn φt(s)2 H1(Ω) +Ck 2tn+1 tn|φt(s)|2+|θt(s)|2+|ct(s)|2+|φtt(s)|2ds, (7.9) where C>0 are different constants independent of (h, k) and independent of the exact solution (φ, θ, c). 586 F. GUILL´ EN-GONZ´ ALEZ AND J.V. GUTI´ ERREZ-SANTACREU We now test (7.4)with2ken+1 θ,d and bound the right-hand side ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ CV|en+1 θ,d |2−|en θ,d|2+|en+1 θ,d −en θ,d|2+2K1k|∇en+1 θ,d |2 =−2K1k∇en+1 θ,i ,∇en+1 θ,d −lkδten+1 φ,d +δten+1 φ,i ,e n+1 θ,d +2kRn+1 θ,e n+1 θ,d  ≤Ckh 2θ(tn+1)2 H2(Ω) +K1k|∇en+1 θ,d |2+Ch 2tn+1 tn φt(s)2 H1(Ω)ds +α 2k|δten+1 φ,d |2+CC Vk|en+1 θ,d |2+Ck 2tn+1 tn (θtt(s)2 H1(Ω)+φtt(s)2 H1(Ω))ds. (7.10) Let us now take xh=2ken+1 c,d as a test function into (7.5), ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ |en+1 c,d |2−|en c,d|2+|en+1 c,d −en c,d|2+2K2k|∇en+1 c,d |2 =−2K2k∇en+1 c,i ,∇en+1 c,d −2MK 2k[P0cn h]T(1 −[P0cn h]T)∇(en φ,d +en φ,i),∇en+1 c,d  −2MK 2k(1 −[P0cn h]T)en cT∇φ(tn),∇en+1 c,d  +2MK 2kc(tn)en cT∇φ(tn),∇en+1 c,d +2kRn+1 c,e n+1 c,d . We firstly bound the truncated error |en cT|2=|[P0cn h]T−P0c(tn)+P0c(tn)−c(tn)|2 ≤C|[P0cn h]T−P0c(tn)|2+|P0c(tn)−c(tn)|2 ≤C|P0cn h−P0c(tn)|2+h2|∇c(tn)|2 ≤C|en c,d|2+|en c,i|2+h2|∇c(tn)|2 ≤C|en c,d|2+h2c(tn)2 H1(Ω), (7.11) where in the last line we have used the stability property |P0ψ|≤|ψ|. Note that the interpolation error |P0c(tn)−c(tn)|appearing in (7.11)isonlyoforderO(h), independent of the finite element approximation. By virtue of (7.11), we bound ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ |en+1 c,d |2−|en c,d|2+|en+1 c,d −en c,d|2+K2k|∇en+1 c,d |2 ≤Ckh 2c(tn+1)2 H2(Ω) +Ck|∇en φ,d|2+Ckh 2∇φ(tn)2 L∞(Ω)c(tn)2 H1(Ω) +Ckh 2φ(tn)2 H2(Ω) +Ck∇φ(tn)2 L∞(Ω)|en c,d|2 +Ck 2tn+1 tn |∇φt(s)|2ds +Ck 2∇φ(tn+1)2 L∞tn+1 tn |ct(s)|2ds+Ck 2tn+1 tn ctt(s)2 H1(Ω)ds. (7.12) Again, C>0 are different constants independent of (h, k) and independent of the exact solution (φ, θ, c). STABILITY AND CONVERGENCE OF TWO DISCRETE SCHEMES 587 By adding (7.9), (7.10)and(7.12) and applying the generalized discrete Gronwall lemma, we establish the following estimate for all n<N and for ksmall enough: en+1 φ,d 2 H1(Ω) +1 4ε2en+1 φ,d 4 L4(Ω) +CV|en+1 θ,d |2+|en+1 c,d |2+αk n  l=0 δtel+1 φ,d  2 +K1k n  l=0 |∇el+1 θ,d |2+K2k n  l=0 |∇el+1 c,d |2≤exp(C1T)C2h2+C3k2,(7.13) where C1= C1+φ2 L2(0,T;W1,∞(Ω)) 1−Ckmax{1,φ2 L∞(0,T;H1(Ω))} C2=Cφ4 L∞(0,T;H1(Ω))φ2 L2(0,T;H2(Ω)) +φ2 L2(0,T;H2(Ω)) +φt2 L2(0,T;H1(Ω)) +θ2 L2(0,T;H2(Ω)) +φ2 L2(0,T;W1,∞(Ω))c2 L∞(0,T;H1(Ω)) +c2 L2(0,T;H2(Ω)), and C3=Cφt2 L2(0,T;L2(Ω)) +θt2 L2(0,T;L2(Ω)) +ct2 L2(0,T;L2(Ω)) +φtt2 L2(0,T;L2(Ω)) +θtt2 L2(0;T;H1(Ω))+φtt2 L2(0,T;H1(Ω))+φt2 L2(0,T;H1(Ω)) +φ2 L∞(0,T;W1,∞(Ω))ct2 L2(0,T;L2(Ω))ds+ctt2 L2(0,T;H1(Ω)), with C>0 independent of (h, k) and independent of the exact solution (φ, θ, c). The same estimates are obtained for the total errors by using the interpolation errors; hence (7.7)canbe deduced.  Remark 7.2. Observe that, to obtain error estimates in the previous theorem, the monotony property (analogous to (3.8)) 1 2ε2((en+1 φ,d )3−en φ,d)(en+1 φ,d −en φ,d)≥F(en+1 φ,d )−F(en φ,d) is not used, because the corresponding initial term is Ω F(e0 φ,d)= 1 8ε2Ω ((e0 φ,d)2−1)2=O1 ε2· Therefore, by using this property, the error estimates of order O(k+h) does not hold, because in this case the final bound remains of order O(k+h+1/ε2). 7.2. Errorestimatesforthelinearscheme As the derivation of the error equations for θn+1 hand cn+1 his exactly the same as in the previous nonlinear scheme, we only treat in this section the error equation for φn+1 h,whichisgivenby ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αδten+1 φ,d ,x h+∇en+1 φ,d ,∇xh+en+1 φ,d ,x h+1 2ε2(en φ,d)3,x h =−αδten+1 φ,i ,x h−1 2ε2(en φ,i)3,x h−3 2ε2P1 h(φ(tn))φ(tn)en φ,i,x h −3 2ε2φn hP1 h(φ(tn))en φ,d,x h+1 2ε2en φ,d +en φ,i,x h +en φ,d +en φ,i,x h+β ε2en θ,d −(θA+θB)en c,d,x h+Rn+1 φ+Sn+1 φ,x h (7.14)