Estimates and Computations for Melting and Solidification Problems

J. M. Greenberg^{†}

Carnegie Mellon University Department of Mathematical Sciences Pittsburgh, PA 15213

1 Introduction

In this note we focus on rigorous estimates for melting and solidification problems. These estimates guarantee the convergence of solutions to nontrivial equilibrium patterns. We also present some interesting numerical simulations which demonstrate the equilibrium structures and the approach of the system to equilibrium. The novel feature of these calculations is the linking of the small parameter in the system, δ, to the grid spacing.

Our governing equations are an energy balance for the medium and an evolution equation for a nonconservative “order parameter,” p. In dimensionless variables the energy equation takes the form

T_{t} + (1 − p^{2})^{m}p_{t} − a^{2}ΔT = 0, (x, y) ∈ Ω  (1.1)

while the evolution equation for p may be written as

αδp_{t} − λ^{2}δ^{2}Δp = μT (1 − p^{2})^{m} + p (1 − p^{2}) , (x, y) ∈ Ω.  (1.2)

Here Δ is the two-dimensional Laplacian and T is the dimensionless temperature. The parameters α, a, λ, and μ are positive and order 1, 0 < δ << 1, and m is a nonnegative integer. The equations with m = 0 are the most frequently studied, the case where m = 2 has been investigated by Sekerka and coworkers [1-3] and Almgren and Almgren [4]. For definiteness we take m = 1 and restrict our attention to the situation where Ω is the unit square

Ω = {(x, y) | − 1/2 < x < 1/2 and − 1/2 < y < 1/2} .  (1.3)

On ∂Ω we assume that the normal derivatives of T and p vanish. Initial conditions for T and p are prescribed. Noting that at a local maximum (minimum) of p that Δp ≤ 0(Δp ≥ 0) it is easily checked that if the initial data for p satisfies −1 ≤ p(x, y, 0^{+}) ≤ 1, (x, y) ∈ \overline{Ω}, then the same estimates obtain for all future t > 0. We shall exploit this fact in what follows.

The system (1.1) and (1.2) is customarily referred to as the diffuse-interface description of melting and solidification. There is a corresponding sharp-interface description which is obtained formally from (1.2) by letting δ tend to 0^{+}. The derivation of this description may be found in Caginalp [5-6].

In the sharp-interface description p only takes on the values −1 and +1. Solid regions are those where p = −1 and liquid regions those where p = +1. In what follows

†This research was partially supported by the Applied Mathematical Sciences Program, U.S. Department of Energy and the Mathematics and Computer Science Division, Army Research Office.

Γ(t) = {(x, y) | x = x̂(s, t), y = ŷ(s, t), 0 ≤ s ≤ L}  (1.4)

will be a curve separating liquid and solid regions and n(s, t) will be the unit normal to Γ which points into the liquid region. The unit tangent, t, to Γ is chosen so that n(s, t) × t(s, t) = e_{3}. K is the curvature of Γ at (x̂, ŷ)(s, t), and

c(s, t) =^{def} (x̂_{t}, ŷ_{t}) · n(s, t)  (1.5)

is the speed at which Γ moves in the direction n. A singular perturbation analysis of (1.1) and (1.2), with m = 1, yields the following results:

1. T is continuous across Γ;

2. On Γ the Gibbs-Thomson relation

αc + λ^{2}δK = −2^{1/2}λμT_{M}  (1.6)

holds; and finally across Γ

\frac{4}{3}c = a^{2} (\frac{∂T^{−}}{∂n} − \frac{∂T^{+}}{∂n} ) .  (1.7)

In (1.6), T_{M}(s) = T (x̂(s, t), ŷ(s, t)) is the local melt temperature while in (1.7), \frac{∂T^{+}}{∂n} is the normal

derivative of T in the liquid region and \frac{∂T^{−}}{∂n} is the normal derivative of T in the solid region. In the sharp-interface description one solves the heat equation

T_{t} − a^{2}ΔT = 0  (1.8)

in the portion of Ω away from interfaces. At equilibrium, c = 0 and (1.6) implies that T_{M} = 0(δ). In what follows, we shall restrict our attention to initial data T (x, y, 0^{+}) = 0(δ^{\frac{1}{2}} ).

The sharp-interface description becomes unduly complicated in the presence of multiple interfaces since all must be tracked. Difficulties also arise when interfaces merge. For these, and a plethora of other reasons, we focus on (1.1) and (1.2). Our principal results consist of a sequence of a-priori inequalities satisfied by solutions of the system (1.1) and (1.2). These estimates guarantee regularity of the solutions and that the solutions converge to equilibria as t tends to infinity.

One point worthy of note is it is possible to obtain solutions of (1.1) and (1.2) which approximate those generated by the sharp interface descriptions. These are produced by solving (1.1) and (1.2) with initial data p^{δ}(x, y, 0^{+}) which approximate the jump discontinuous data of the sharp interface theory but transition smoothly from −1 to +1 over an interval of width of 0(δ) in the direction n. For such initial p’s (and initial temperature fields satisfying T (x, y, 0^{+}) = 0(δ^{\frac{1}{2}} )) we are guaranteed that terms like δ_{∫ ∫Ω}| ▽ p^{δ}|^{2}dxdy and \frac{1}{δ}∫ ∫ _{Ω} (1 − p^{2})^{2} dxdy are bounded independently of δ and t and it is these estimates which link the solutions of the two theories.

We also note there is an “intermediate” model where one replaces (1.1) (here with m = 1) with

T (x, y, t) ≡ T_{0}(t) = h_{0} − ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (p − p^{3}/3) (x, y, t)dxdy  (1.9)

and solves

αδp_{t} − λ^{2}δ^{2}Δp = (μT_{0} + p) (1 − p^{2}) , (x, y) ∈ Ω  (1.10)

and

\frac{∂p}{∂n} = 0 , (x, y) ∈ ∂Ω.  (1.11)

The constant h_{0} in (1.9) is given by

h_{0} = ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (T + p − p^{3}/3) (x, y, 0^{+})dxdy.  (1.12)

This “intermediate” model is justified by two facts. The first is that solutions of (1.1) and (1.2) satisfying the boundary conditions \frac{∂T}{∂n} = \frac{∂p}{∂n} = 0, (x, y) ∈ ∂Ω satisfy the overall energy balance:

h_{0} ≡ ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (T + p − p^{3}/3) (x, y, t)dxdy  (1.13)

for all t ≥ 0 and the second is that the spatially varying portion of T decays to zero as t tends to infinity (for details see the comments following (2.17)).

2 Entropies and A-Priori Estimates

Our goal in this section is a set of inequalities of the form

\frac{∂η}{∂t} − div q ≤ G  (2.1)

satisfied by solutions of

T_{t} + (1 − p^{2}) p_{t} − a^{2}ΔT = 0 , (x, y) ∈ Ω,  (2.2)

αδp_{t} − λ^{2}δ^{2}Δp = (μT + p) (1 − p^{2}) , (x, y) ∈ Ω,  (2.3)

and

\frac{∂T}{∂n} = \frac{∂p}{∂n} = 0 , (x, y) ∈ ∂Ω.  (2.4)

Once again

Ω = {(x, y) | − 1/2 < x < 1/2 and − 1/2 < y < 1/2)} .  (2.5)

For (2.1) to be of value we require

0 ≤ η, (x, y) ∈ Ω and q · n = 0, (x, y) ∈ ∂Ω  (2.6)

and that either G ≤ 0 in Ω or G is a-priori bounded in L_{1}(Ω × (0, ∞)). Throughout we exploit the pointwise estimates

−1 ≤ p(x, y, t) ≤ 1 , (x, y) ∈ Ω and t ≥ 0  (2.7)

satisfied by solutions of (2.3).

Once again we record the overall energy balance satisfied by solutions of (2.2) – (2.4), namely the identity:

∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (T + p − p^{3}/3) (x, y, t)dxdy ≡ h_{0}.  (1.13)

This is obtained by integrating (2.2) over Ω and exploiting (2.4)_{1}.

Our first inequality is obtained by multiplying (2.2) by (\frac{μT}{δ} ) and (2.3) by \frac{p_{t}}{δ} and adding the resulting identities. This yields (2.1) with

η_{1} = \frac{μT^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})^{2}}{4δ} ,  (2.8)

q_{1} = \frac{μa^{2}}{δ} T ▽ T + λ^{2}δp_{t} ▽ p,  (2.9)

and

G_{1} = −αp_{t}^{2} − \frac{μa^{2}| ▽ T |^{2}}{δ}  (2.10)

and (2.8)-(2.10) yield

\begin{matrix}\frac{d}{dt}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{μT^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})}{4δ}^{2}) (x, y, t)dxdy \\ ≤ −∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (αp_{t}^{2} + \frac{μa^{2}}{δ} | ▽ T |^{2}) (x, y, t)dxdy. \end{matrix}  (2.11)

Once again we restrict our attention to initial data satisfying

∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{T^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})}{4δ}^{2}) (x, y, 0^{+})dxdy = 0(1)

independent of δ.

If we next multiply (2.2) by αT_{t} we obtain

\frac{αa^{2}}{2} \frac{∂}{∂t}| ▽ T |^{2} − αa^{2} div (T_{t} ▽ T) = −α (T_{t}^{2} + (1 − p^{2}) p_{t}T_{t})

and if we add this result to (2.1) with η_{1}, q_{1}, G_{1} we obtain

\frac{∂η_{2}}{∂t} − div q_{2} ≤ G_{2}

where

η_{2} = \frac{μT^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})}{4δ}^{2} + \frac{αa^{2}}{2} | ▽ T |^{2}  (2.12)

q_{2} = \frac{μa^{2}}{δ} T ▽ T + λ^{2}δp_{t} ▽ p + αa^{2}T_{t} ▽ T  (2.13)

and

G_{2} = −α (p_{t}^{2} + T_{t}^{2} + (1 − p^{2}) T_{t}p_{t}) − \frac{μa^{2}}{δ} | ▽ T |^{2} ≤ \frac{−α}{2} (p_{t}^{2} + T_{t}^{2}) − \frac{μa^{2}}{δ} | ▽ T |^{2}  (2.14)

and (2.12)-(2.14) imply that

\begin{matrix}\frac{d}{dt}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{μT^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})}{4δ}^{2} + \frac{αa^{2}}{2} | ▽ T |^{2}) (x, y, t)dxdy \\ ≤ −∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{α}{2} (p_{t}^{2} + T_{t}^{2}) + \frac{μa^{2}}{δ} | ▽ T |^{2}) (x, y, t)dxdy. \end{matrix}  (2.15)

If we now write

T (x, y, t) = T_{0}(t) + T_{1}(x, y, t) , (x, y) ∈ Ω  (2.16)

where

T_{0}(t) = ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} T (x, y, t)dxdy and ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} T_{1}(x, y, t)dxdy ≡ 0,  (2.17)

then (2.15) implies that if T and p are in H^{1}(Ω) at t = 0, then they are in H^{1}(Ω) for all t > 0 and that

T satisfies lim_{t→∞}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}T_{1}^{2}(x, y, t)dxdy = 0. ^{1}

If we make a similar decomposition of p, that is write

p(x, y, t) = p_{0}(t) + p_{1}(x, y, t)  (2.18)

where

p_{0}(t) = ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} p(x, y, t)dxdy and ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} p_{1}(x, y, t)dxdy = 0,  (2.19)

then (2.3) implies that

αδ \frac{dp_{0}}{dt} = ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}(μT + p)(1 − p^{2})(x, y, t)dxdy  (2.20)

and

αδ \frac{d^{2}p_{0}}{dt^{2}} = ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (μ (T_{t} (1 − p^{2}) − 2pT p_{t}) + (1 − 3p^{2}) p_{t}) (x, y, t)dxdy  (2.21)

||f||_{k}^{2} =^{def} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} ∑_{j=0}^{k} ( ∑_{\begin{matrix}α_{1}+α_{2}=j \\ α_{i}≥0 \end{matrix}} (\frac{∂^{j}f}{∂x^{α_{1}}∂y^{α_{2}}} )^{2}) (x, y)dxdy < ∞

|||f|||_{k}^{2} =^{def} ∫_{0}^{s_{∗}} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} ∑_{j=0}^{k} ( ∑_{\begin{matrix}α_{1}+α_{2}+α_{3}=j \\ α_{i}≥0 \end{matrix}} (\frac{∂^{j}f}{∂x^{α_{1}}∂y^{α_{2}}∂s^{α_{3}}} )^{2}) (x, y, s)dxdy) ds < ∞.

1Recall, H^{k}(Ω) consists of functions f satisfying

and H^{k}(Ω × (0, s_{∗})) consists of all functions f on Ω × (0, s_{∗}) satisfying

and (2.20) and (2.21) along with the fact that (2.15) implies that p_{t} and T_{t} are L_{2}(Ω × (0, ∞)) combine

to yield the result that \frac{dp_{0}}{dt} and \frac{d^{2}p_{0}}{dt^{2}} are in L_{2}(0, ∞) and that

lim_{t→∞}αδ \frac{dp_{0}}{dt} = lim_{t→∞} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}(μT + p)(1 − p^{2})(x, y, t)dxdy = 0.  (2.22)

With the estimate (2.15) we can glean weak information about the long-term behavior of solutions of (2.2)-(2.4). In what follows we let

(T^{t}, p^{t})(x, y, s) = (T, p)(x, y, t + s), (x, y) ∈ Ω and 0 ≤ s ≤ s_{∗}.  (2.23)

The inequality (2.15), together with (2.7), implies that for any s_{∗} the functions T^{t} and p^{t} are in H^{1}(Ω × (0, s_{∗})) with bounds which depend only on the H^{1}(Ω) norms of T (x, y, 0^{+}) and p(x, y, 0^{+}) and s_{∗} and that

lim_{t→∞} ∫_{0}^{s_{∗}} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{α}{2} ((T_{s}^{t})^{2} + (p_{s}^{t})^{2}) + \frac{μa^{2}}{δ} | ▽ T^{t}|^{2}) (x, y, s)dxdy) ds = 0.  (2.24)

Moreover, (2.3) implies that λ^{2}δ^{2}Δp^{t} ∈ L_{2}(Ω × (0, s_{∗})) with bounds which are independent of t. The boundedness of T^{t} and p^{t} in H^{1}(Ω × (0, s_{∗})) implies we can find an increasing sequence {t_{i}}_{i=1}^{∞} with lim_{i→∞}t_{i} = ∞ and funcitons T^{∞} and p^{∞} in H^{1}(Ω × (0, s_{∗})) and Δp^{∞} in L_{2}(Ω × (0, s_{∗})) with the following properties: (T^{∞}, p^{∞}) is the strong L_{2}(Ω × (0, s_{∗})) limit of (T^{t_{i}}, p^{t_{i}}) while T_{s}^{∞}, p_{s}^{∞}, ▽T^{∞}, ▽p^{∞}, and Δp^{∞} are the weak L_{2}(Ω × (0, s_{∗})) limits of the appropriate derivatives of T^{t_{i}} and p^{t_{i}}. Additionally, (2.24) implies

T_{s}^{∞} = p_{s}^{∞} = 0 and ▽ T^{∞} = 0 in Ω × (0, s_{∗})  (2.25)

and thus that T^{∞} is a constant and p^{∞} is independent of s. The energy balance (1.13) yields

T^{∞} = h_{0} − ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (p^{∞} − (p^{∞})^{3} /3) (x, y)dxdy  (2.26)

while (2.7) implies that −1 ≤ p^{∞} ≤ 1. (2.3) and (2.24) also imply that p^{∞} satisifes the equilibrium equation

λ^{2}δ^{2}Δp^{∞} + (μT^{∞} + p^{∞}) (1 − (p^{∞})^{2}) = 0, (x, y) ∈ Ω  (2.27)

and

\frac{∂p^{∞}}{∂n} = 0, (x, y) ∈ ∂Ω  (2.28)

where T^{∞} is given by (2.26) and finally (2.27) and (2.28) imply that p^{∞} is a critical point of the functional

N (p^{∞}) =^{def} \frac{μ (T^{∞})^{2}}{2δ} + ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{λ^{2}δ}{2} | ▽ p^{∞}|^{2} + \frac{(1 − (p^{∞})^{2})}{4δ}^{2}) (x, y)dxdy.  (2.29)

Additionally, the functional

η(t, s_{∗}) =^{def} \frac{1}{s_{∗}} ∫_{0}^{s_{∗}}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{μ(T^{t})^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p^{t}|^{2} + \frac{(1 − (p^{t})^{2})^{2}}{4δ} + \frac{αa^{2}}{2} | ▽ T^{t}|^{2}) (x, y, s)dxdyds (2.30)

is decreasing in t and has a limit η_{∞} satisfying

N (p^{∞}) ≤ η_{∞}.  (2.31)

To obtain stronger information we examine the time differentiated versions of (2.2)-(2.4), namely the equations:

T_{tt} + (1 − p^{2})p_{tt} − 2pp_{t}^{2} − a^{2}ΔT_{t} − 0, (x, y) ∈ Ω,  (2.32)

αδp_{tt} − λ^{2}δ^{2}Δp_{t} = (μT_{t} + p_{t})(1 − p^{2}) − 2p(μT + p)p_{t}, (x, y) ∈ Ω,  (2.33)

\frac{∂T_{t}}{∂n} = \frac{∂p_{t}}{∂n} = 0 , (x, y) ∈ ∂Ω.  (2.34)

We shall also exploit Lemma 1. Suppose f ∈ H^{1}(Ω). Then, for every ϵ > 0

∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} f^{4}(x, y)dxdy ≤ 4 [\frac{(1 + ϵ)}{ϵ} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} f^{2}(x, y)dxdy + ϵ ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} | ▽ f|^{2}(x, y)dxdy]^{2} .

We multiply (2.33) by p_{t} and make use of (2.7) to obtain

\frac{αδ}{2} \frac{∂}{∂t}p_{t}^{2} − λ^{2}δ^{2} div (p_{t} ▽ p_{t}) ≤ −λ^{2}δ^{2}| ▽ p_{t}|^{2} + (\frac{μT_{t}^{2}}{2} + (1 + \frac{μ}{2}) p_{t}^{2}) + 2μ|T |p_{t}^{2}.  (2.35)

We note that (2.15) implies the term \frac{μT_{t}}{2}^{2} + (1 + \frac{μ}{2}) p_{t}^{2} is in L_{1}(Ω × (0, ∞)) and thus we confine our attention to 2μ|T |p_{t}^{2}. Schwarz’s inequality and Lemma 1 imply

\begin{matrix}2μ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (|T |p_{t}^{2}) (x, y, t)dxdy ≤ 2μ (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}T^{2}(x, y, t)dxdy)^{1/2} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{4}(x, y, t)dxdy)^{1/2} \\ ≤2μ (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}T^{2}(x, y, t)dxdy)^{1/2}((\frac{1+ϵ}{ϵ}) ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy+ϵ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}|▽p_{t}|^{2}(x, y, t)dxdy) . \end{matrix} (2.36)

Moreover, (2.11) guarantees that

(∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} T^{2}(x, y, t)dxdy)^{1/2} ≤ (\frac{2δ}{μ} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} η_{1}(x, y, 0^{+})dxdy)^{1/2}  (2.37)

and this guarantees that we may choose ϵ so that

2μϵ (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} T^{2}(x, y, t)dxdy)^{1/2} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} | ▽ p_{t}|^{2}(x, y, t) ≤ \frac{λ^{2}δ^{2}}{2} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} | ▽ p_{t}|^{2}(x, y, t)dxdy. (2.38)

With ϵ so chosen, (2.35) implies that

\begin{matrix}\frac{αδ}{2} \frac{d}{dt} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t) ≤ \frac{−λ^{2}δ^{2}}{2} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} | ▽ p_{t}|^{2}(x, y, t)dxdy + \frac{μ}{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}T_{t}^{2}(x, y, t)dxdy \\ + (1 + \frac{μ}{2} + 2^{3/2} (μδ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}η_{1}(x, y, 0^{+})dxdy)^{1/2} (1 + \frac{1}{ϵ})) ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy \end{matrix} (2.39)

and (2.39) then implies that if

p_{t}(x, y, 0^{+}) = \frac{1}{αδ} (λ^{2}δ^{2}Δp + (μT + p)(1 − p^{2})) (x, y, 0^{+})

is in L_{2}(Ω), then p_{t} and thus Δp are in L_{2}(Ω) independently of t and ▽p_{t} ∈ L_{2}(Ω × (0, ∞)).

The preceding estimate points up the advantage of the “intermediate” model discussed briefly in (1.9)-(1.12) over the full system. For the intermediate model the contankerous term 2μ|T |p_{t}^{2} in (2.35) would have been replaced by 2μ|T_{0}|p_{t}^{2} and this latter term is easily controlled since we know by our previous estimates that

|T_{0}| ≤ (\frac{2δ}{μ} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} η_{1}(x, y, 0^{+})dxdy)^{1/2}

and thus we would not have to invoke Lemma 1 and choose ϵ = 0 (δ^{3/2}). This latter choice, pumped

up the coefficient of ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy in (2.39) to be 0 (1/δ^{3/2}) which is far from optimal. If we

chose to add some multiple of (2.39) to (2.15) and wanted the right hand side of the resulting inequality to be negative definite our multiplier would have to be 0(δ^{3/2}), an unacceptably small number. For the “intermediate” model we could use an 0(1) multiplier (for details see (2.61)-(2.63)).

On the other hand, if we multiply (2.32) by (T_{t} + (1 − p^{2})p_{t})/a^{2} we find that

\begin{matrix}\frac{1}{2a^{2}} \frac{∂}{∂t} (T_{t} + (1 − p^{2})p_{t})^{2} = \frac{a^{2}}{2} \frac{∂}{∂t}(ΔT)^{2} \\ = div ((T_{t} + (1 − p^{2})p_{t}) ▽ T_{t}) − | ▽ T_{t}|^{2} − (1 − p^{2}) ▽ T_{t} · ▽p_{t} + 2pp_{t} ▽ T_{t} · ▽p \end{matrix}  (2.40)

and (2.40) implies that

\begin{matrix}\frac{a^{2}}{2} \frac{d}{dt}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}(ΔT)^{2}(x, y, t)dxdy ≤ \\ −\frac{1}{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}| ▽ T_{t}|^{2}(x, y, t)dxdy + ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}| ▽ p_{t}|^{2}(x, y, t)dxdy \\ +4∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}| ▽ p|^{2}(x, y, t)dxdy. \end{matrix}  (2.41)

The first term on the right hand side of (2.41) is in controllable by virtue of the minus sign and the second term is L_{1}(0, ∞) by virtue of (2.39). Thus we can confine our attention to the third term. To show this term is L_{1}(0, ∞) we note that

∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}|▽p|^{2}(x, y, t)dxdy≤(∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{4}(x, y, t)dxdy)^{1/2} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}|▽p|^{4}(x, y, t)dxdy)^{1/2}. (2.42)

Lemma 1 then guarantees that

(∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{4}(x, y, t)dxdy)^{1/2} ≤ \frac{(1 + \sqrt{5})}{2} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (p_{t}^{2} + | ▽ p_{t}|^{2}) (x, y, t)dxdy)  (2.43)

and that

\begin{matrix}(∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}| ▽ p|^{4}(x, y, t)dxdy)^{1/2} ≤ (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{x}^{4}(x, y, t)dxdy)^{1/2} + (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{y}^{4}(x, y, t)dxdy)^{1/2} \\ ≤ \frac{(1 + \sqrt{5})}{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (p_{x}^{2} + p_{y}^{2} + p_{xx}^{2} + p_{yy}^{2} + 2p_{xy}^{2}) (x, y, t)dxdy \\ ≤ \frac{(1 + \sqrt{5})}{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (p_{x}^{2} + p_{y}^{2} + 4(Δp)^{2}) (x, y, t)dxdy \end{matrix} (2.44)

and (2.43) and (2.44) together with (2.41) guarantee that if T (x, y, 0^{+}) and p(x, y, 0^{+}) are in H^{2}(Ω) then T (x, y, t) is in H^{2}(Ω) for all future times and additionally that ▽T_{t} ∈ L_{2}(Ω × (0, ∞)).

The estimates (2.39) and (2.44) also guarantee that if T (x, y, 0^{+}) and p(x, y, 0^{+}) are in H^{2}(Ω), then the resulting orbit is compact in H^{1}(Ω) and that the w-limit set is nonempty in H^{1}(Ω). The fact that ▽T_{t} ∈ L_{2}(Ω × (0, ∞)) further guarantees that the function T_{1} defined in (2.16) and (2.17) converges to zero strongly in H^{1}(Ω).

Our final estimates pertain to solutions generated by initial data T (x, y, 0^{+}) ∈ H^{2}(Ω) and p(x, y, 0^{+}) ∈ H^{3}(Ω). These are obtained by multiplying (2.33) by p_{tt}. The resulting inequality is

\begin{matrix}\frac{λ^{2}δ^{2}}{2} \frac{∂}{∂t}| ▽ p_{t}|^{2} + αδp_{tt}^{2} − λ^{2}δ^{2} div (p_{tt} ▽ p_{t}) = \\ μ (T_{t}(1 − p^{2}) − 2pT p_{t}) p_{tt} + (1 − 3p^{2})p_{t}p_{tt} \\ ≤ \frac{αδ}{2} p_{tt}^{2} + \frac{2μ^{2}}{αδ} (T_{t}^{2} + 4T^{2}p_{t}^{2}) + \frac{8}{αδ} p_{t}^{2}. \end{matrix}  (2.45)

Our preceding calculations imply that if T (x, y, 0^{+}) ∈ H^{2}(Ω) and p(x, y, 0^{+}) ∈ H^{3}(Ω), then

\frac{2μ^{2}}{αδ} (T_{t}^{2} + 4T^{2}p_{t}^{2}) + \frac{8}{αδ} p_{t}^{2}

is in L_{1}(Ω × (0, ∞)) and thus (2.45) implies that ▽p_{t} ∈ L_{2}(Ω) independently of t and that p_{tt} ∈ L_{2}(Ω×(0, ∞)). The former conclusion, together with (2.3), implies that p ∈ H^{3}(Ω) with bounds which are independent of t while the latter conclusion implies that

lim_{t→∞}αδp_{t} = lim_{t→∞} (λ^{2}δ^{2}Δp + (μT + p)(1 − p^{2})) = 0

and that this limit is obtained strongly in L_{2}(Ω).

We conclude this section with some observations about the “intermediate” model introduced in (1.9)-(1.12). Once again in that model one replaces (2.3) with

αδp_{t} − λ^{2}δ^{2}Δp = (μT_{0} + p)(1 − p^{2}), (x, y) ∈ Ω  (2.46)

and

\frac{∂p}{∂n} = 0 , (x, y) ∈ ∂Ω  (2.47)

where

T_{0}(t) = h_{0}−∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (p − p^{3}/3) (x, y, t)dxdy and h_{0} = ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (T + p − p^{3}/3) (x, y, 0^{+})dxdy. (2.48)

Solutions of this latter system also satisfy the pointwise estimates −1 ≤ p(x, y, t) ≤ 1 and a set of “entropy” inequalities similar to those satisfied by solutions of the full system. The most basic of these is

\frac{∂η̃_{1}}{∂t} − div q̃_{1} = G̃_{1}  (2.49)

where

η̃_{1} = \frac{μT_{0}^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})^{2}}{4δ}  (2.50)

q̃_{1} = λ^{2}δp_{t} ▽ p  (2.51)

and

G̃_{1} = −αp_{t}^{2}.  (2.52)

Moreover, (2.49) - (2.52) imply that

\begin{matrix}\frac{d}{dt}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{μT_{0}^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})}{4δ}^{2}) (x, y, t)dxdy \\ ≤ −α∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy. \end{matrix}  (2.53)

The time differentiated identities

αδp_{tt} − λ^{2}δ^{2}Δp_{t} = μ (Ṫ_{0}(1 − p^{2}) − 2pT_{0}p_{t}) + (1 − 3p^{2})p_{t},  (2.54)

Ṫ_{0} = − ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}(1 − p^{2})p_{t}(x, y, t)dxdy,  (2.55)

and

\frac{∂p_{t}}{∂n} = 0 , (x, y) ∈ ∂Ω.  (2.56)

also yields:

\begin{matrix}\frac{αδ}{2} \frac{∂p_{t}^{2}}{∂t} − λ^{2}δ^{2} div (p_{t} ▽ p_{t}) = \\ −λ^{2}δ^{2}| ▽ p_{t}|^{2} + μ (Ṫ_{0}(1 − p^{2}) − 2pT_{0}p_{t}) p_{t} + (1 − 3p^{2})p_{t}^{2} \\ ≤ −λ^{2}δ^{2}| ▽ p_{t}|^{2} + \frac{μ}{2} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy) \\ +2^{3/2}μ^{1/2}δ^{1/2} (∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}η̃_{1}(x, y, 0^{+})dxdy)^{1/2} p_{t}^{2} + (1 + \frac{μ}{2})p_{t}^{2} \end{matrix}  (2.57)

and

\begin{matrix}\frac{λ^{2}δ^{2}}{2} \frac{∂| ▽ p_{t}|^{2}}{∂t} − λ^{2}δ^{2} div (p_{tt} ▽ p_{t}) = \\ −αδp_{tt}^{2} + μ (Ṫ_{0}(1 − p^{2}) − 2pT_{0}p_{t})p_{tt} + (1 − 3p^{2})p_{t}p_{tt} \\ ≤ −\frac{αδ}{2} p_{tt}^{2} + \frac{1}{αδ} (2μ^{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy + (4 + 16μδ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}η̃_{1}(x, y, 0^{+})dxdy) p_{t}^{2}) . \end{matrix} (2.58)

The former follows from multiplying (2.54) by p_{t} and the latter from multiplying (2.54) by p_{tt}. Moreover, if we let

β_{2} = (1 + μ + 2^{3/2}μ^{1/2}δ^{1/2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}η̃_{1}(x, y, 0^{+})dxdy) ,  (2.59)

β_{3} = 4 + 2μ^{2} + 16μδ ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} η̃_{1}(x, y, 0^{+})dxdy,  (2.60)

η̃_{2} = η̃_{1} + \frac{α^{2}δ}{8β_{2}} p_{t}^{2} = \frac{μT_{0}^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})^{2}}{4δ} + \frac{α^{2}δ}{8β_{2}} p_{t}^{2},  (2.61)

q̃_{2} = q̃_{1} + \frac{αλ^{2}δ^{2}}{4β_{2}} p_{t} ▽ p_{t} = λ^{2}δp_{t} ▽ p + \frac{αλ^{2}δ^{2}}{4β_{2}} p_{t} ▽ p_{t},  (2.62)

\begin{matrix}G̃_{2} = G̃_{1} + \frac{α}{4β_{2}} (1 + \frac{μ}{2} + 2^{3/2}μ^{1/2}δ^{1/2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}η̃_{1}(x, y, 0^{+})dxdy) p_{t}^{2} \\ +\frac{αμ}{8β_{2}} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy − \frac{αλ^{2}δ^{2}}{4β_{2}} | ▽ p_{t}|^{2}, \end{matrix}  (2.63)

η̃_{3} = η̃_{2} + \frac{α^{2}λ^{2}δ^{3}}{8β_{3}} | ▽ p_{t}|^{2} = \frac{μT_{0}^{2}}{2δ} + \frac{λ^{2}δ}{2} | ▽ p|^{2} + \frac{(1 − p^{2})}{4δ} + \frac{α^{2}δ}{8β_{2}} p_{t}^{2} + \frac{α^{2}λ^{2}δ^{3}}{8β_{3}} | ▽ p_{t}^{2}|  (2.64)

q̃_{3} = q̃_{2} + \frac{α^{2}λ^{2}δ^{3}}{4β_{3}} p_{tt} ▽ p_{t} = λ^{2}δp_{t} ▽ p + \frac{αλ^{2}δ^{2}}{4β_{2}} p_{t} ▽ p_{t} + \frac{α^{2}λ^{2}δ^{3}}{4β_{3}} p_{tt} ▽ p_{t},  (2.65)

and

\begin{matrix}G̃_{3} = G̃_{2} + \frac{α}{4β_{3}} (4 + 16μδ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}η̃_{1}(x, y, 0^{+})dxdy) p_{t}^{2} \\ +\frac{α}{2β_{3}} μ^{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy − \frac{α^{2}δ^{2}}{8β_{3}} p_{tt}^{2} \end{matrix}  (2.66)

we obtain the inequalities

\frac{∂η̃_{2}}{∂t} − div q̃_{2} ≤ G̃_{2} and \frac{∂η̃_{3}}{∂t} = divq̃_{3} ≤ G̃_{3}  (2.67)

and these imply

\begin{matrix}\frac{d}{dt}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{μT_{0}^{2}}{2δ} + \frac{λ^{2}δ| ▽ p|^{2}}{2} + \frac{(1 − p^{2})^{2}}{4δ} + \frac{α^{2}δ}{8β_{2}} p_{t}^{2}) (x, y, t)dxdy \\ ≤ −\frac{3α}{4}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t)dxdy − \frac{αλ^{2}δ^{2}}{4β_{2}} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}| ▽ p_{t}|^{2}(x, y, t)dxdy, \end{matrix}  (2.68)

and

\begin{matrix}\frac{d}{dt}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} (\frac{μT_{0}^{2}}{2δ} + \frac{λ^{2}δ| ▽ p|^{2}}{2} + \frac{(1 − p^{2})^{2}}{4δ} + \frac{αδ}{8β_{2}} p_{t}^{2} + \frac{α^{2}λ^{2}δ^{3}}{8β_{3}} | ▽ p_{t}|^{2}) (x, y, t)dxdy \\ ≤ −\frac{α}{2}∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{t}^{2}(x, y, t) − \frac{αλ^{2}δ^{2}}{4β_{2}} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}| ▽ p_{t}|^{2}(x, y, t)dxdy \\ −\frac{α^{2}δ^{2}}{4β_{3}} ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}p_{tt}^{2}(x, y, t)dxdy. \end{matrix}  (2.69)

We leave it to the reader to supply the precise conclusions derivable from the inequalities (2.53), (2.63), and (2.69). The key point of these inequalities is that the parameter δ appears optimally.

3 Computational Experiments

In this section we present some numerical simulations of the system (2.2) - (2.4). We let N be an integer and set

dx = dy = 1/N.  (3.1)

Our approximate solutions will be piecewise constant on the cells

Ω_{i,j} = \{(x, y) |−\frac{1}{2} + \frac{(i − 1)}{N} < x < −\frac{1}{2} + \frac{i}{N} and − \frac{1}{2} + \frac{(j − 1)}{N} < y < −\frac{1}{2} + \frac{j}{N} \}  (3.2)

where 1 ≤ i, j ≤ N and thus the cell averages (T_{i,j}, p_{i,j}) may also be thought of as the point values at the cell centers

(x_{i}, y_{j}) = (−\frac{1}{2} + \frac{2i − 1}{2N}, −\frac{1}{2} + \frac{2j − 1}{2N}) .  (3.3)

We impose the boundary conditions

\frac{∂T}{∂n} = \frac{∂p}{∂n} = 0 , (x, y) ∈ ∂Ω  (3.4)

by insisting that the extended grid functions satisfy

\begin{matrix}& (T_{0,k}, p_{0,k}) = (T_{1,k}, p_{1,k}) , \\ & (T_{N+1,k}, p_{N+1,k}) = (T_{N,k}, p_{N,k}) , \\ & (T_{k,0}, p_{k,0}) = (T_{k,1}, p_{k,1}) , \\ and \\ & (T_{k,N+1}, p_{k,N+1}) = (T_{k,N}, p_{k,N}) \end{matrix} \}  (3.5)

for 1 ≤ k ≤ N.

For any function u defined on the grid and satisfying (3.5) we let

Δ_{N}u_{i,j} = N^{2} (u_{i+1,j} + u_{i−1,j} + u_{i,j+1} + u_{i,j−1} − 4u_{i,j})  (3.6)

be the discrete laplacian applied to u. For functions u and v satisfying (3.5) we define

[u, v] = ∑_{i,j=1}^{N}u_{i,j}v_{i,j}/N^{2}  (3.7)

and note that

−[u, Δ_{N}u] = ∑_{i=1}^{N}∑_{j=1}^{N−1} (u_{i,j+1} − u_{i,j})^{2} + ∑_{j=1}^{N}∑_{i=1}^{N−1} (u_{i+1,j} − u_{i,j})^{2} .  (3.8)

Such grid functions also satisfy the discrete version of Lemma 1, namely the inequality

[u^{2}, u^{2}] = [u^{4}, 1] ≤ 4 (\frac{(1 + ϵ)}{ϵ}[u, u] − ϵ[u, Δ_{N}u])^{2}  (3.9)

for every ϵ > 0.

Finally, we replace (2.2) - (2.4) by the system

Ṫ_{i,j} + (1 − p_{i,j}^{2})ṗ_{i,j} − a^{2}Δ_{N}T_{i,j} = 0  (3.10)

and

αδṗ_{i,j} − λ^{2}δ^{2}Δ_{N}p_{i,j} = (μT_{i,j} + p_{i,j})(1 − p_{i,j}^{2})  (3.11)

where T_{i,j} and p_{i,j} satisfy the discrete boundary conditions (3.5) and the symbol · denotes differentiation with respect to t. Solutions of (3.10) and (3.11) obey the pointwise estimate that if −1 ≤ p_{i,j}(0^{+}) ≤ 1, 1 ≤ i, j ≤ N, then −1 ≤ p_{i,j}(t) ≤ 1, 1 ≤ i, j ≤ N for all t > 0. This system also satisfies estimates which are the exact analogue of those obtained in section (2). Corresponding to (2.8) we have the identity

[T + p − p^{3}/3, 1](t) ≡ h_{0}  (3.12)

while the analogues of (2.11) and (2.15) are the inequalities

\frac{d}{dt} ([(\frac{μT^{2}}{2δ} + \frac{(1 − p^{2})}{4δ}^{2}) , 1] − \frac{λ^{2}δ}{2} [p, Δ_{N}p]) (t) = −α[p_{t}, p_{t}](t) + \frac{μa^{2}}{δ} [T, Δ_{N}T](t),  (3.13)

and

\begin{matrix}\frac{d}{dt} ([\frac{μT^{2}}{2δ} + \frac{(1 − p^{2})^{2}}{4δ} , 1] − \frac{λ^{2}δ}{2} [p, Δ_{N}p] − \frac{αa^{2}}{2} [T, Δ_{N}T]) (t) \\ ≤ −\frac{α}{2} ([p_{t}, p_{t}] + [T_{t}, T_{t}]) (t) + \frac{μa^{2}}{δ} [T, Δ_{N}T] (t). \end{matrix}  (3.14)

We spare the reader the detailed analogues of (2.39), (2.41), and (2.45). Suffice it to say if we

replace terms like ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2}f(x, y, t)dxdy by [f, 1] (t) and terms like ∫_{−1/2}^{1/2} ∫_{−1/2}^{1/2} | ∇f |^{2} (x, y, t)dxdy

by − [f, Δ_{N}f] (t), then the resulting inequalities go through intact.

When doing simulations with the discrete system (3.10) and (3.11) we link the small parameter δ to the grid dx = dy = 1/N; specifically we set

δ = 1/N  (3.15)

and restrict our attention to initial data satisfying

N [\frac{μT^{2}}{2} + \frac{(1 − p^{2})^{2}}{4} , 1] (0^{+}) − \frac{λ^{2}}{2N} [p, Δ_{N}p] (0^{+}) − \frac{αa^{2}}{2} [T, Δ_{N}T] (0^{+}) = 0(1).  (3.16)

Interesting data satisfying (3.16) may be achieved in a variety of ways. Below we outline one such construction. We let Γ be a simple closed curve of finite arclength L(Γ) which is wholly contained in Ω, S_{c} be all cells Ω_{i,j} which are cut by Γ, S_{I} be all cells Ω_{i,j} which are interior to the region surrounded by Γ, and S_{E} be those cells Ω_{i,j} which are in Ω and exterior to Γ. For definiteness we choose

p_{i,j}(0^{+}) = \{ \begin{matrix}1, & if Ω_{i,j} ∈ S_{E} \\ p_{i,j}^{∗} ∈ [−1, 1] , & if Ω_{i,j} ∈ S_{c} \\ −1, & if Ω_{i,j} ∈ S_{I} \end{matrix}  (3.17)

and

T_{i,j} (0^{+}) = h_{0} − [(p − p^{3}/3) (0^{+}) , 1] = 0 (1/N^{1/2}) .  (3.18)

Such data models the situation where initially water surrounds ice. With such data it is easily checked that N [T^{2} (0^{+}) , 1] = 0(1), while N [(1 − p^{2})^{2}, 1] and −\frac{1}{N} [p, Δ_{N}p] are 0(L(Γ)) and thus the functional defined in (3.16) is 0(1) independently of N. For such data we note that δ^{3} [Δ_{N}p, Δ_{N}p] = \frac{1}{N^{3}} [Δ_{N}p, Δ_{N}p] is also 0(L(Γ)) independently of N. The last functional enters in the discrete analogue of (2.68).

The simulations we present were all run with the following parameters:

(α, λ, a, μ) = (1, 1, 1, 1) .  (3.19)

We chose N = 200

dx = dy = δ = 1/200  (3.20)

and

dt = dx/10.  (3.21)

With these parameters

m1 =^{def} dt/(dx)^{2} = 10  (3.22)

and

m2 =^{def} \frac{δdt}{(dx)^{2}} = \frac{1}{10}  (3.23)

and thus we were compelled to use an implicit integration scheme for (3.10). Our choice was ADI working with the fixed tri-diagonal matrix whose non-zero entries are

A(i, i) = \{ \begin{matrix}1 + m1 & if i = 1 and i = 200 \\ 1 + 2m1 & if 2 ≤ i ≤ 199 \end{matrix}  (3.24)

and

A(i, i + 1) = A(i + 1, i) = −m1, 1 ≤ i ≤ 199.  (3.25)

We used a first order explicit Euler scheme on (3.11).

Our first simulation was run with the following choice of initial data:

p_{i,j}(0^{+}) = \{ \begin{matrix}1, if \{101 + i ≤ j ≤ 200 and 1 ≤ i ≤ 99\} or \{302 − i ≤ j ≤ 200 and 102 ≤ i ≤ 200\} or \{1 ≤ j ≤ 100 − i and 1 ≤ i ≤ 99\} or \{1 ≤ j ≤ −101 + i and 102 ≤ i ≤ 200\}, \\ 0, if \{100 + i = j and 1 ≤ i ≤ 100\} or \{301 − i = j and 101 ≤ i ≤ 200\} or \{101 − i = j and 1 ≤ i ≤ 100\} or \{−100 + i = j and 101 ≤ i ≤ 200\}, \\ −1, if \{102 − i ≤ j ≤ 99 + i and 2 ≤ i ≤ 100\} or \{−99 + i ≤ j ≤ 300 − i and 101 ≤ i ≤ 199\}. \end{matrix}  (3.26)

This data satisfies [p, 1] = [p^{3}, 1] = 0. We chose h_{0} = .05 and T_{i,j}(0^{+}) ≡ .05, 1 ≤ i, j ≤ 200.

The initially superheated square block of ice converged to a circular block centered at (0, 0) with a radius consistent with the δ = 0^{+} asymptotics. Specifically, if we grant that the “equilibrium” curve separating the ice from the water is a circle of radius R, then the “Gibbs-Thomson” relation (1.6), with our choice of parameters, implies that

\frac{.005}{2^{1/2}R} = −T_{M}  (3.27)

while (3.12) implies that

T_{M} = .05 − [p − p^{3}/3, 1].  (3.28)

Noting that p_{i,j} ≃ −1 inside of the circle and p_{i,j} ≃ +1 outside yields

[p − p^{3}/3, 1] = \frac{2}{3} (1 − 2πR^{2})  (3.29)

and (3.27)-(3.29) imply that R satisfies

\frac{4π}{3} R^{2} + \frac{.005}{2^{1/2}R} = \frac{2}{3} − .05.  (3.30)

(3.30) has one positive solution close to zero and a second larger solution R = .3808 which is the approximate equilibrium radius we obtain in our simulation. This radius, when substituted into the Gibbs-Thomson relation (3.27) yields a melt temperature consistent with what we observe. The results of this simulation are shown in Figure 1.

Our final simulation starts with 16 regularly spaced circular blocks of ice surrounded by water. Specifically

p_{i,j}(0^{+}) = \{ \begin{matrix}−1, if & \sqrt{(i − 40m)^{2} + (j − 40n)^{2}} ≤ 15 \\ & and 1 ≤ m, n ≤ 4 \\ 1, & otherwise. \end{matrix}  (3.31)

For this data

[(p − p^{3}/3)(0^{+}), 1] = .2885.  (3.32)

We again chose h_{0} to be .05 and at t = 0 took the uniform temperature field

T_{i,j}(0^{+}) = h_{0} − [(p − p^{3}/3)(0^{+}), 1] = −.2385.  (3.33)

The solutions using this initial data converged to the same equilibrium as before but the approach was far more interesting. The solutions at times t = .05, .17, .24, .4, .6, .7, 2, 3.3, and 12.5 are shown in Figures 2-10.

We also ran both sets of data with the “intermediate” model and obtained virtually identical results.

Explanation of Graphics

In Figures 1-10 the upper left frame is a contour plot of p and p is 1 in a neighbourhood of the boundaries x = ±1/2 and y = ±1/2 and the upper right hand frame is a contour plot of T. In Figure 1 the lower left plot shows p, T and p − p^{3}/3 + T along the lines x = 0, y = 0, and the 45^{o} diagonal through the origin. These plots overlay one another because of the circular symmetry of the final state. In Figure 1 the bottom right plot again displays the same information. In Figures 2-10 the bottom left plot shows p, T, and p − p^{3}/3 + T along the 45^{o} diagonal through the origin and the bottom right plot shows the same quantities on the line x = 0.

Figure 1: 21 Figure 2:

22 Figure 3:

23 Figure 4:

24 Figure 5:

25 Figure 6:

26 Figure 7:

27 Figure 8:

28 Figure 9:

29 Figure 10:

30 References

1. S.L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, and G.B. Mc-Fadden, “Thermodynamically-Consistent Phase-Field Models, Physica D 69 (1993) 189-200.

2. G.B. McFadden, A.A. Wheeler, R.J. Brown, S.R. Coriell and R.F. Sekerka, “Phase-Field Models for Anisotropic Interfaces,” Phys. Rev. E 48 (1993) 2016-2024.

3. S.L. Wang and R.F. Sekerka, “Algorithms for Phase Field Computations of the Dendritic Operating State at Large Supercoolings,” J. Comp. Physics 127 (1996) 110-117.

4. A.S. Almgren and R.F. Almgren, “Phase Field Instabilities and Adaptive Mesh Refinement,” in Modern Methods for Modeling Microstructure in Materials, TMS/SIAM 1996.

5. G. Caginalp, “Stefan and Hele-Shaw Type Models as Asymptotic Limits of the Phase Field Equation,” Phys. Rev. A 39 (1989) 5887-5896.

6. G. Caginalp, “Phase Field Models and Sharp Interface Limits: Some Differences in Subtle Situations,” Rocky Mountain J. Math. 21 (1996), 2, 603-616.