Time discretization of a nonlinear phase field system in general domains

This paper deals with the nonlinear phase field system \begin{equation*} \begin{cases} \partial_t (\theta +\ell \varphi) - \Delta\theta = f&\mbox{in}\ \Omega\times(0, T), \\[1mm] \partial_t \varphi - \Delta\varphi + \xi + \pi(\varphi) = \ell \theta,\ \xi\in\beta(\varphi)&\mbox{in}\ \Omega\times(0, T) \end{cases} \end{equation*} in a general domain $\Omega\subseteq\mathbb{R}^N$. Here $N \in \mathbb{N}$, $T>0$, $\ell>0$, $f$ is a source term, $\beta$ is a maximal monotone graph and $\pi$ is a Lipschitz continuous function. We note that in the above system the nonlinearity $\beta+\pi$ replaces the derivative of a potential of double well type. Thus it turns out that the system is a generalization of the Caginalp phase field model and it has been studied by many authors in the case that $\Omega$ is a bounded domain. However, for unbounded domains the analysis of the system seems to be at an early stage. In this paper we study the existence of solutions by employing a time discretization scheme and passing to the limit as the time step $h$ goes to $0$. In the limit procedure we face with the difficulty that the embedding $H^1(\Omega) \hookrightarrow L^2(\Omega)$ is not compact in the case of unbounded domains. Moreover, we can prove an interesting error estimate of order $h^{1/2}$ for the difference between continuous and discrete solutions.


Introduction and results
In the present contribution we address a nonlinear phase field system in a general domain Ω ⊆ R N and discuss it using a time discretization procedure, which turns out to be useful and efficient for the approximation and the existence proof. Moreover, we show an error estimate of order 1/2 in suitable norms for the difference between continuous and discrete solutions.
Our analysis moves from the consideration of the following simple version of the phasefield system of Caginalp type (cf. [8,22]; one may also see the monographs [7,23,43]): in Ω × (0, T ), (1.2) where Ω is the three-dimensional domain in which the evolution takes place, T is some final time, θ refers to the relative temperature around some critical value, that is taken to be 0 without loss of generality, and φ is the order parameter. Moreover, c 0 , ℓ, κ and η are positive constants, f is a source term and F ′ is the derivative of a smooth double-well potential F , whose prototype reads F reg (r) = 1 4 (r 2 − 1) 2 , r ∈ R. (1. 3) In our approach, we aim to consider other kinds of potentials, which are by now widely known and extensively used in the mathematical literature, that are the logarithmic potential and the double obstacle potential: F log (r) = (1 + r) ln(1 + r) + (1 − r) ln(1 − r) − c 1 r 2 if r ∈ (−1, 1), (1.4) F obs (r) = I(r) − c 2 r 2 , r ∈ R. (1.5) Here, the coefficient c 1 in (1.4) is larger than 1, so that F log admits a double well, and c 2 in (1.5) is an arbitrary positive constant, whereas the function I in (1.5) is the indicator function of [−1, 1], i.e., it takes the values 0 or +∞ according to whether or not r belongs to [−1, 1]. Note that F log can be extended by continuity to the closed interval [−1, 1], but its derivative F ′ log turns out to be singular as the variable approaches −1 from the right or +1 from the left. On the other hand, F obs is even a non-smooth potential and for it it is no longer possible to consider the derivative, but we have to deal with the subdifferential of I. Hence, in order to be as general as possible, in our investigation we allow the potential F to be just the sum where β is a convex function that is allowed to take the value +∞ somewhere, and π is a smooth perturbation which may be concave. In such a case, β is supposed to be proper and lower semicontinuous, so that its subdifferential is well defined and can replace the derivative which might not exist. Let us point out that, in the case of a multivalued subdifferential like ∂I, equation (1.2) becomes a differential inclusion.
The system (1.1)-(1.2) is of course complemented by initial conditions like θ(0) = θ 0 and φ(0) = φ 0 , and suitable boundary conditions. Concerning the latter, we set the usual homogeneous Neumann condition for both θ and φ, that is, where ∂ ν denotes differentiation with respect to the outward normal of ∂Ω. Indeed, by these conditions we are assuming that there is no flow exchange with the exterior of Ω.
Equations (1.1)-(1.2) yield a system of phase field type. Such systems have been introduced (cf. [7,8,43]) in order to include phase dissipation effects in the dynamics of moving interfaces arising in thermally induced phase transitions (the reader may also see [6,13,34,38,39]). In our framework, we are actually considering the following form for the total free energy: where c 0 and ℓ stand for the specific heat and latent heat coefficients, respectively, with a terminology motivated by earlier studies (see [21]) on the Stefan problem; let us also refer to the monography [23], which deals with phase change models as well. In this connection, it is worth to introduce the enthalpy e by where δG δθ denotes the variational derivative of G with respect to θ, so that (1.7) yields e = c 0 θ + ℓϕ. Then the governing balance and phase equations are given by where q denotes the thermal flux vector, f represents some heat source and δG δϕ is the variational derivative of G with respect to ϕ. Hence (1.9) reduces exactly to (1.2) along with the homogeneous Neumann boundary condition for ϕ. Moreover, if we assume the classical Fourier law q = −κ ∇θ, then (1.8) is nothing but the usual energy balance equation of the Caginalp model. Moreover, (1.1) follows from (1.8) and the Neumann boundary condition for θ is a consequence of the no-flux condition q · n = 0 on the boundary. We notice that the above phase field system has received a good deal of attention in the last decades [1,9,12,15,19,35] and can be deduced as a special gradientflow problem (cf., e.g., [41] and references therein).
Let us point out that questions related to the well-posedness, long-time behaviour of solutions and optimal control problems have been investigated for the Caginalp system (1.1)-(1.2) and for some variation or extension of this phase field system. Without any sake of completeness, we mention the contributions [2,10,11,14,18,20,22,25,29,36,37,40] for various qualitative analyses and [5,16,17,27,28] for some related control problems.
For the sake of simplicity, in the sequel of the paper we simply let c 0 = κ = η = 1 by observing that our treatment can be easily extended to the case of coefficients different from 1. On the other hand, we keep the parameter ℓ in both equations (1.1) and (1.2) since ℓ plays a role in the estimates.
The case of unbounded domains has the difficult mathematical point that compactness methods cannot be applied directly (related discussions can be found, e.g., in [24,[30][31][32][33]). It would be interesting to construct an applicable theory for the case of unbounded domains and to set assumptions for the case of unbounded domains by trying to keep some typical features in previous works, that is, in the case of bounded domains. By considering the case of unbounded domains, it may be possible to make a new finding which is not given or known in the case of bounded domains. Also, the new finding would be useful for other studies of partial differential equations.
In this paper we consider the initial-boundary value problem on a general domain for the nonlinear phase field system (P) where Ω is a bounded domain or an unbounded domain in R N with smooth bounded boundary ∂Ω (e.g., Ω = R N \ B(0, R), where B(0, R) is the open ball with center 0 and radius R > 0) or Ω = R N (in this case, no boundary condition should be prescribed) or Ω = R N + . Moreover, we let ℓ > 0 and deal with the following conditions (A1)-(A4): (A1) β ⊂ R × R is a maximal monotone graph with effective domain D(β) and β(r) = ∂ β(r), where ∂ β denotes the subdifferential of a proper lower semicontinuous convex function β : R → [0, +∞] satisfying β(0) = 0.
Please note that the three functions in (1.3)-(1.5) actually satisfy (A1) and (A2), indeed we have that Let us define the Hilbert spaces respectively, and with the related Hilbertian norms. Moreover, we use the notation This paper is organized as follows. In Section 2 we introduce a time discretization of (P), set precisely the approximate problem, and state the main theorems. Section 3 contains the proof of the existence for the discrete problem. In Section 4 we establish uniform estimates for the approximate problem and pass to the limit. Section 5 show error estimates between solutions of (P) and solutions of the approximate problem.

Time discretization and main results
We will prove existence of solutions to (P) by employing a time discretization scheme. More precisely, we will establish existence for (P) by passing to the limit in the problem for a.a. t ∈ (nh, (n + 1)h), n = 0, ..., N − 1, we can rewrite (P) n as the problem in Ω.
as the reader can check directly by using the definitions (2.2) and (2.3).
We define solutions of (P) as follows.
is called a solution of (P) if (θ, ϕ, ξ) satisfies Now the main results read as follows.
Then there exists a unique solution of (P).

Existence of discrete solution
In this section we will prove Theorem 2.1.

Error estimates
In this section we will prove Theorem 2.3. for all h ∈ (0, h 1 ).
Remark 5.1. The above argument can be used also to the check that f h → f strongly in L 2 (0, T ; H) (5.7) as h ց 0, in the case when f is just in L 2 (0, T ; H). Indeed, by density, for all ε > 0 there exists a function g ∈ W 1,1 (0, T ; H) such that f − g L 2 (0,T ;H) < ε. for all h ∈ (0, h) (cf. (5.5) and (5.6)). Hence this gives a proof of (5.7).