Lobachevskii Journal of Mathematics Vol. 13, 2003, 15 – 24 \begin{matrix}c \\ ◯ \end{matrix}M.A.Ignatieva and A.V.Lapin

M.A.Ignatieva and A.V.Lapin

MIXED HYBRID FINITE ELEMENT SCHEME FOR STEFAN PROBLEM WITH PRESCRIBED CONVECTION

Abstract. We construct a mixed hybrid finite element scheme of lowest order for the Stefan problem with prescribed convection and suggest and investigate an iterative method for its solution. In the iterative method we use a preconditioner constructed by using ”standard” finite element approximation of Laplace operator on a finer grid.

The proposed approach develops the results of [1], where a spectrally equivalent preconditioner for the condensed matrix in mixed hybrid finite element approximation for linear elliptic equation was constructed.

1. Introduction

Stefan problem with prescribed convection serves as a mathematical model for the heat transfer and solidification process in the metal casting (see [2, 3]). Commonly used numerical methods of its solving are based on the implicit or semi-implicit mesh approximations in time variable with lowest order finite element approximation in space variables [4, 5].

Because in the applied problems both the temperature fields and fluxes are of practical interest, mixed and mixed hybrid finite element schemes appear as important method for its numerical solution.

Mixed and mixed hybrid finite element schemes are thoroughly investigated for the linear boundary-value problems (cf. [6, 7] and bibliography therein), while a few publications have concern with these methods for nonlinear problems, especially for free and moving boundary problems. We construct a mixed hybrid finite element scheme for a Stefan problem which is a case of moving boundary problem.

The main purpose of the article is to suggest an effective iterative algorithm to solve this finite element scheme. We construct and investigate an iterative

ditioners.

2000 Mathematical Subject Classification. Primary. 65M55, 65M60.

Key words and phrases. Mixed hybrid discretization, condensed matrices, variational inequalities, Stefan problem, iterative methods, spectrally equivalent precon-

This work was supported by RFBR, project No 01-01-00068.

process with a preconditioner, the rate of convergence of this algorithm does not depend on the mesh size. On the other hand, the implementation of the iterative method reduces to the solution of a ”standard” finite-dimensional variational inequality, which can be made by using any of coordinate or gradient relaxation methods.

2. Mathematical model

Let Ω ⊂ ℝ^{2} be a domain with piecewise smooth boundary ∂Ω = Γ_{D} ∪ Γ_{N}. We consider the following nonlinear problem: find (u(x, t), θ(x, t)) such that

\{ \begin{matrix}\frac{∂θ}{∂t} + v \frac{∂θ}{∂x_{1}} − Δu = 0, & in Ω, t > 0, \\ u = z(x, t) & on Γ_{D}, t > 0, \\ \frac{∂u}{∂n} = g(x, t) & on Γ_{N}, t > 0, \\ θ(x, t) ∈ H(u(x, t)) & in Ω, t > 0, \\ θ(x, 0) = θ_{0}(x) & in Ω̄, \end{matrix}  (1)

where n is the unit vector of outward normal, v = const > 0, z, g and θ_{0} are given functions. We consider the case when the graph of H : ℝ^{1} → ℝ^{1} monotonically increases and contains a vertical segment and suppose that the function H has a single values at points on Γ_{D}.

Problem (1) can serve as a simplified model of continuous casting process, where u is the temperature of casting metal, θ(u) is the enthalpy function and v is the casting speed in x_{1} direction. The enthalpy function has a mentioned above property for example in the case of copper casting.

The existence and uniqueness of a weak solution for problem (1) are studied in [8, 9].

3. Semi-discretization

First, we introduce the semi-discretization of problem (1) using the characteristics of the first order differential operator and constant steps τ in time. Namely, if (x_{1}, x_{2}, t) is the point at the time level t we use the following approximation:

(\frac{∂}{∂t} + v \frac{∂}{∂x_{1}} ) H ≈ \frac{1}{τ} (H(x_{1}, x_{2}, t) − H(x̃_{1}, x_{2}, t − τ)) , x̃_{1} = x_{1} − vτ.

If x_{1} − vτ < 0 then we put x̃_{1} = 0.

After semi-discretization problem (1) on each time level can be formally written in the pointwise form as

\{ \begin{matrix}−Δu + P u ∋ f & in Ω, \\ u = z & on Γ_{D}, \\ \frac{∂u}{∂n} = g & on Γ_{N}, \end{matrix}  (2)

where P u = H(u)/τ is the multivalued maximal monotone nonlinear operator and the right-hand side f = H̃(u)/τ also arises due to semi-discretization.

Because of maximal monotonicity of function H and the property Dom(H) = ℝ the operator P u is the subdifferential of convex continuous function φ(u). The weak formulation of problem (2) can be written in the form of variational inequality: find u(x) ∈ H^{1}(Ω), u(x) = z(x) on Γ_{D} such that

∫_{Ω} ∇u∇(q − u)dx + φ(q) − φ(u) ⩾ ∫_{Ω} f(q − u)dx + ∫_{Γ_{N}} g(q − u)dΓ (3)

∀q ∈ H^{1}(Ω), q(x) = z(x) on Γ_{D}.

It is well known that (3) has a unique solution [10].

4. Mixed hybrid formulation of the problem

Now, let v = ∇u be so called flux function, then we get the following mixed formulation of (2):

\begin{cases}v − ∇u = 0 & in Ω, \\ div v − P u ∋ −f & in Ω, \\ u = z & on Γ_{D}, \\ v · n = g & on Γ_{N}. \end{cases}  (4)

Let H(div, Ω) = {w ∈ L_{2}(Ω)^{n} : div w ∈ L_{2}(Ω)} with the norm ∥w∥^{2} = ∫_{Ω}(|w|^{2} + | div w|^{2})dx, H(div, Ω) = {w ∈ H(div, Ω) : w · n ∈ L_{2}(∂Ω)} with the norm ∥w∥_{H}^{2} = ∥w∥^{2} + ∫_{∂Ω}(w · n)^{2}dΓ and subspaces H_{N}(div, Ω) = {w ∈ H(div, Ω) : w · n̄ = g a. e. on Γ_{N}}, H_{N}^{0}(div, Ω) = {w ∈ H(div, Ω) : w · n = 0 a. e. on Γ_{N}}. Now by a weak solution of problem (4) we mean a triple (u, v, σ) ∈ L_{2}(Ω) × H_{N} × L_{2}(Ω), such that

\begin{cases}∫_{Ω} v · w dx + ∫_{Ω} u div wdx − ∫_{Γ_{D}} z(w · n)dΓ = 0 ∀w ∈ H_{N}^{0}, \\ ∫_{Ω} div v qdx − ∫_{Ω} σ(x)qdx = − ∫_{Ω} fqdx ∀q ∈ L_{2}(Ω), \\ σ(x) ∈ P u(x) for a. e. x ∈ Ω. \end{cases}  (5)

Note, that by construction if u is solution of (5) then (u, ∇u, σ) is a solution of problem (5).

Let Ω̄ = ⋃_{i=1}^{m} ē_{i} be a partitioning of the domain into m nonoverlapping subdo-

mains, where e_{i} has a piecewise smooth boundary. Hereafter we suppose that the parts Γ_{D} and Γ_{N} of the boundary ∂Ω are composed by the whole sides ∂e_{i}. The common sides of elements e_{i} and e_{j} we denote by Γ_{ij}, where i ≠ j, i, j = \overline{1, m}. We suppose that ∂Ω consists of s segments, which we denote as Γ_{1}, …, Γ_{s}. Let the intersections ∂e_{i} with Γ_{j}, j = \overline{1, s} be denoted by Γ_{i,m+j}, i = \overline{1, m}, j = \overline{1, s} and let Γ_{i,m+j} for j = \overline{1, s_{1}}, s_{1} < s compose Γ_{N}, while Γ_{i,m+k}, i = \overline{1, m}, k = \overline{s_{1} + 1, s} compose Γ_{D}.

Let further v = (v_{1}, …, v_{m}), u = (u_{1}, …, u_{m}). Then system (4) can be written in the following form

\begin{cases}v_{i} − ∇u_{i} = 0, div v_{i} − σ_{i} = −f_{i}, σ_{i} ∈ P u_{i} in e_{i}, \\ u_{i} − u_{j} = 0, v_{i} · n_{ij} − v_{j} · n_{ij} = 0 on Γ_{ij}, i, j = \overline{1, m}, \\ u_{i} = z on Γ_{i,m+k}, k = \overline{s_{1} + 1, s}, \\ v_{i} · n_{i} = g on Γ_{i,m+j}, j = \overline{1, s_{1}}, \end{cases}  (6)

where n_{i} is the outward normal vector to ∂e_{i} and n_{ij} is the unit normal vector to Γ_{ij} directed from e_{i} to e_{j}. Let us note that in semi-discretized problem the flux function v · n is continuous in Ω though in initial differential problem flux may have a jump.

On the basis of (6) we define a weak mixed hybrid formulation of problem

(2). Let U = ∏_{1⩽i⩽m} L_{2}(e_{i}), V = ∏_{1⩽i⩽m} H(div, e_{i}) and Λ = ∏_{i>j} L_{2}(Γ_{ij}). We

introduce the bilinear forms M : V × V → ℝ, B : U × V → ℝ, C : Λ × V → ℝ and D : U × U → ℝ by the following equalities

M(v, w) = ∑_{i=1}^{m} ∫_{e_{i}} v_{i} · w_{i}dx, B(u, w) = ∑_{i=1}^{m} ∫_{e_{i}} u_{i} div w_{i}dx,

C(λ, w) = ∑_{i=1}^{m} ( ∑_{j=i+1}^{m} ∫_{Γ_{ij}} λ_{ij}(w_{j} − w_{i}) · n_{ij}dΓ − ∑_{j=1}^{s_{1}} ∫_{Γ_{i,m+j}} λ_{i,m+j}(w_{i} · n_{i,m+j})dΓ),

D(q, σ) = ∑_{i=1}^{m} ∫_{e_{i}} σ_{i}q_{i}dx

and linear functionals

F (q) = ∑_{i=1}^{m} ∫_{e_{i}} fq_{i}dx, ℓ(μ) = ∑_{i=1}^{m} ∑_{j=1}^{s_{1}} ∫_{Γ_{i,m+j}} gμ_{i,m+j}dΓ,

r(w) = ∑_{i=1}^{m} ∑_{k=s_{1}+1}^{s} ∫_{Γ_{i,m+k}} z(w_{i} · n_{i})dΓ.

The weak mixed hybrid formulation is as follows: find (u, v, λ, σ) ∈ U ×V × Λ × U such that

\begin{cases}M(v, w) + B(u, w) + C(λ, w) = r(w) ∀w ∈ V, \\ B(q, v) + D(q, σ) = −F (q) ∀q ∈ U, \\ C(μ, v) = −ℓ(μ) ∀μ ∈ Λ, \\ σ(u) ∈ P̄u, \end{cases}  (7)

where P̄u = (H(u_{1})/τ, …, H(u_{m})/τ) for u ∈ U.

Proposition 1. The problems (3) and (7) are equivalent in the following sense:

if u is the solution of (3), then (u, v, λ, σ) with u_{i} = u|_{e_{i}} – restriction of u to e_{i}, v_{i} = ∇u_{i} a. e. in e_{i}, λ_{ij} = u_{i} a. e. on Γ_{ij} and σ ∈ P̄(u) is the solution of (7);

backwards, if (u, v, λ, σ) is the solution of problem (7), then the function u : u|_{e_{i}} = u_{i} is a solution to problem (3).

5. Approximation

Let further Ω be a polygonal domain and τ_{h} = {e_{1}, e_{2}, …, e_{m}} be its conforming triangulation [11]. We assume that all e_{i} are convex polygons.

Let U_{h}, V_{h} and Λ_{h} be finite element subspaces of the corresponding spaces

U_{h} = ∏_{1⩽i⩽m} U_{ih}, V_{h} = ∏_{1⩽i⩽m} V_{ih}

and refer to Λ_{h} as the space of vector-functions with constant components. Here V_{ih} is a finite dimensional subspace in H(div, e_{i}) of dimension n_{i} consisting of vector-functions v_{i} ∈ H(div, e_{i}) such that v_{i} · n_{ij} is a constant on the corresponding interface Γ_{ij} and U_{ih} is a subspace of U consisting of vectorfunctions such that each its coordinate is a constant on each subdomain e_{i}. The value of n_{i} is equal to the total number of interfaces Γ_{ij} belonging to the

boundary of e_{i}, i = \overline{1, m}. Thus, the dimension of V_{h} is equal to n̂ = ∑_{i=1}^{m} n_{i}, the

dimension of U_{h} is equal to m, and the dimension of Λ_{h} is equal to ň where ň is the total number of interfaces.

The finite element approximation of (7) reads as follows: find (u_{h}, v_{h}, λ_{h}, σ_{h}) ∈ U_{h} × V_{h} × Λ_{h} × U_{h} satisfying the following relations:

\begin{cases}M(v_{h}, w) + B(u_{h}, w) + C(λ_{h}, w) = r(w) ∀w ∈ V_{h}, \\ B(q, v_{h}) + D(q, σ(u_{h})) = −F (q) ∀q ∈ U_{h}, \\ C(μ, v_{h}) = −ℓ(μ) ∀μ ∈ Λ_{h}, \\ σ_{h} ∈ P̄u. \end{cases}  (8)

Let now v_{ij} be the degrees of freedom for vector-function v_{ih}, associated with Γ_{ij}, u_{i} be the degrees of freedom for the function u_{h}, associated with e_{i}, and λ_{ij} be the degrees of freedom for λ_{h}, associated with Γ_{ij}, j > i. The algebraic formulation of (8) is: to find (v̄, ū, λ̄, σ̄) such that

\begin{cases}M̂(v̄, w̄) + B̂(ū, w̄) + Ĉ(λ̄, w̄) = r̂(w̄), \\ B̂(q̄, v̄) + D̂(q̄, σ(ū)) = −F̂(q̄), \\ Ĉ(μ̄, v̄) = −ℓ̂(μ̄), \\ σ̄ ∈ P̄ū, \end{cases}  (9)

∀(w̄, q̄, μ̄) ∈ ℝn̂ × ℝ^{m} × ℝň. Here

M̂(v̄, w̄) = ∑_{i=1}^{m}(M_{i}v̄_{i}, w̄_{i}), v̄_{i}, w_{i} ∈ ℝ^{n_{i}},

B̂(ū, w̄) = ∑_{i=1}^{m} u_{i}(− ∑_{j=1}^{i−1} w_{ij}|Γ_{ij}| + ∑_{j=i+1}^{m+s} w_{ij}|Γ_{ij}|),

Ĉ(λ̄, w̄) = ∑_{i=1}^{m}( ∑_{j=i+1}^{m} λ_{ij}(w_{ji} − w_{ij})|Γ_{ij}| − ∑_{i=1}^{s_{1}} λ_{i,m+j}w_{i,m+j}|Γ_{i,m+j}|),

D̂(q̄, σ̄) = ∑_{i=1}^{m} σ_{i}q_{i}mes(e_{i})

are the bilinear forms on ℝn̂×ℝn̂, ℝ^{m}×ℝn̂, ℝň×ℝn̂ and ℝ^{m}×ℝn̂, respectively, and

r̂(w̄) = ∑_{i=1}^{m} ∑_{j=s_{1}+1}^{s} z_{i}(w_{i} · n_{i})dΓ, ∫_{Γ_{i,m+j}} zdΓ, i = \overline{1, m}, j = \overline{s_{1} + 1, s},

F̂(q̄) = ∑_{i=1}^{m} f_{i}q_{i}, f_{i} = ∫_{e_{i}} fdx, i = \overline{1, m},

ℓ̂(μ̄) = ∑_{i=1}^{m}(∑_{j=1}^{s_{1}} g_{ij}μ_{i,m+j}), g_{ij} = ∫_{Γ_{i,m+j}} gdΓ, i = \overline{1, m}, j = \overline{1, s_{1}}

are linear forms defined on ℝň, ℝ^{m} and ℝň.

The entries m_{k,l}^{(i)} of matrices M_{i} are defined by the standard way using the L_{2}(e_{i}) scalar products of the nodal basis functions of the subspaces V_{ih}.

In matrix-vector form problem (9) can be written as follows:

A (v̄ ; ū ; λ̄) + (\begin{matrix}0 \\ −P̄(ū) \\ 0 \end{matrix}) ∋ ( r̄ ; −f̄ ; −ḡ) , A = (\begin{matrix}M & B^{T} & C^{T} \\ B & 0 & 0 \\ C & 0 & 0 \end{matrix} ) .  (10)

We rewrite the system as

\begin{cases}Mv̄ + B^{T} ū + C^{T} λ̄ = r̄, \\ Bv̄ − P̄ū ∋ −f̄, \\ Cv̄ = −ḡ \end{cases}  (11)

and eliminate v̄ from the system. After that we obtain a reduced system which can be written in block form by

(\begin{matrix}BM^{−1}B^{T} & BM^{−1}C^{T} \\ CM^{−1}B^{T} & CM^{−1}C^{T} \end{matrix} ) (ū ; λ̄) + (\begin{matrix}P̄ & ū \\ 0 \end{matrix}) ∋ (\begin{matrix}f + BM^{−1}r̄ \\ g + CM^{−1}r̄ \end{matrix}) .

Using notations

S = (\begin{matrix}BM^{−1}B^{T} & BM^{−1}C^{T} \\ CM^{−1}B^{T} & CM^{−1}C^{T} \end{matrix} ) , μ = (ū ; λ̄) ,

F = (\begin{matrix}f + BM^{−1}r̄ \\ g + CM^{−1}r̄ \end{matrix}) , P̃ = (\begin{matrix}P̄ & 0 \\ 0 & 0 \end{matrix}) ,

we finally obtain the following finite-dimensional problem:

Sμ + P̃μ ∋ F.  (12)

Note that Schur complement matrix S is a symmetric and positive definite matrix [1], while P̃ is a maximal monotone operator. Owing to this fact problem (12) has a unique solution, whence problem (11) has also a unique solution.

6. Iterative method

For the sake of simplicity we analyze only the case of rectangular meshes. We suppose that Ω is the unit square and construct mesh with step h in both directions which defines the partitioning of Ω into elements e_{i}.

To obtain preconditioner for (12) we construct finer grid in Ω with the step h/2. We denote by τ_{h/2} the new partitioning of Ω. Let further W_{h/2} be the piecewise finite element subspace of H^{1}(Ω) and A be the stiffness matrix, corresponding to the approximation in this subspace of Laplace operator with Dirichlet boundary conditions on Γ_{D}. Let the nodes of τ_{h/2} consist of two groups: the first group contains the nodes of τ_{h}, while the second one contains all others (called by the fictitious ones). Then matrix A can be represented in the corresponding block form:

A = (\begin{matrix}A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}) .

It is shown in [1] that the operator S is spectrally equivalent to the Schur complement S_{A}:

S_{A} = A_{11} − A_{12}A_{22}^{−1}A_{21}

with constants of equivalence which do not depend on mesh size:

α(S_{A}μ, μ) ⩽ (Sμ, μ) ⩽ β(S_{A}μ, μ) ∀μ.

In the case under consideration we get α = 1, β = 6.

This observation allows us to use the matrix S_{A} as a preconditioner in iterative process for solving (12):

S_{A} \frac{μ^{n+1} − μ^{n}}{τ} + Sμ^{n} + P̃μ^{n+1} ∋ F.  (13)

The following statement is valid.

Proposition 2. The iterative method (13) converges for any τ ∈ (0, 2/β) and for τ = 2/(α + β) the following estimate holds:

(S_{A}(μ^{n+1} − μ), μ^{n+1} − μ)^{1/2} ⩽ \frac{β − α}{β + α}(S_{A}(μ^{n} − μ), μ^{n} − μ)^{1/2}.

To avoid the explicit calculation of S_{A} on each step of process (13) we use the following trick. We complete the system (12) by the equations in fictitious nodes, so that the algebraic size of resulting system

\begin{cases}Sμ + P̃μ ∋ F, \\ 0 = 0 \end{cases}

is equal to number of fine grid nodes. After we write iterative process with operator A as a preconditioner for this system using block representation of A:

\begin{cases}A_{11} \frac{μ^{n+1} − μ^{n}}{τ} + A_{12} \frac{μ̃^{n+1} − μ̃^{n}}{τ} + Sμ^{n} + P̃μ^{n+1} ∋ F, \\ A_{21} \frac{μ^{n+1} − μ^{n}}{τ} + A_{22} \frac{μ̃^{n+1} − μ̃^{n}}{τ} + 0 = 0, \end{cases}  (14)

where μ̃ are fictitious components. Eliminating from the second equation fictitious values we can verify that it is an equivalent form of process (13) for solving (12) with S_{A} as preconditioner.

For a fixed n problem (14) is the finite element variational inequality with positive definite and symmetric matrix A, so, it can be solved, for example, by using any coordinate relaxation or gradient relaxation method.

On each step of constructed method we need to calculate F − Sμ. It is easy to do if we take into account the identity

F − Sμ^{n} = (\begin{matrix}Bv̄^{n} \\ Cv̄^{n} \end{matrix})

with v̄^{n} = M^{−1}(r̄ − B^{T} ū^{n} − C^{T} λ̄^{n}).

Thus, to solve (14) we get the the following Algorithm

(1) Define μ^{0} = (ū^{0}, λ̄^{0}).

(2) For n ⩾ 0 on each element e_{i} calculate v̄^{n} by formula

v̄_{i}^{n} = M_{i}^{−1}(r̄_{i} − B_{i}^{T} ū_{i}^{n} − C^{T} λ̄_{i}^{n}).

(3) Calculate μ^{n+1} = (ū^{n+1}, λ̄^{n+1}) by solving (14).

(4) n := n + 1 goto step 2.

7. Numerical results

In numerical test we take Ω = (0, 1) × (0, 1), Γ_{D} = {(x_{1}, x_{2}) ∈ ∂Ω : x_{2} = 0}∪{(x_{1}, x_{2}) ∈ ∂Ω : x_{2} = 1}, Γ_{N} = ∂Ω\Γ_{D}. We solve problem in time interval [0,1] using constant time step τ = 0.05 and various grids in space to compare number of iterations. As H we take the function

H(u) = \{ \begin{matrix}0.5u & if u < 0, \\ [0, 1] & if u = 0, \\ u + 1 & if u > 0. \end{matrix}

On the top of square Ω we put g = −3 in the boundary condition, on bottom we use value g = 0, on the left side z = 2, and on the right side z = −1. As initial condition we take u_{0} = −1.

To solve inequality with matrix A on each step of iterative process (14) we use SOR-method, the stopping criterion was ∥μ^{k+1} − μ^{k}∥ ⩽ 10^{−13} and the stopping criterion of outer process was ∥μ^{n+1} − μ^{n}∥ ⩽ 10^{−12}. As iterative parameter in outer process we take τ = 2/7, which we get using estimates of equivalence of matrices S and S_{A}.

In the second column of table 1 we show the number of iteration which we need to solve problem on the first time level (on the next levels the number of iterations becomes smaller) and in the third column – the total number of iteration, which is equal to the sum of iterations on all time steps. From the table we can see that the number of iterations does not depend on mesh size.

Figure

1. Temperature distribution on t = 0.0125 Figure

2. Temperature distribution on t = 0.0750

Figure

3. Temperature distribution on t = 0.0225 Figure

4. Temperature distribution on t = 1

On the figures we show temperature distributions at several time levels.

Grid size Iter 1 Total Iter
11 × 11  81  1195
21 × 21  82  1183
41 × 41  83  1189
81 × 81  83  1189
161 × 161 83  1171

Table 1. Number of iterations

References

[1] Kuznetsov Yu. A. Spectrally equivalent preconditioners for mixed hybrid discretizations of diffusion equations on distorted meshes // J. Numer. Math. – 2003. – V. 11 – P.61–74.

24 M.A.IGNATIEVA AND A.V.LAPIN

[2] Louhenkilpi S., Laitinen E., Nieminen R. Real-time simulation of heat transfer in continuous casting//Metallurgical Trans. B. – 24B. – 1996. – P. 685–693.

[3] Mäkelä M., Männikkö T., Schramm H. Applications of nonsmooth optimization methods to continuous casting of steel//Math.Ins., Univ. Bayreuth. – Rep. No 421 – 1993. – 24 p.

[4] Chen Z., Jiang L. Approximation of a two phase continuous casting problem // J. Part. Diff. Equations. – 1998. – V. 11. – P. 59–72.

[5] Laitinen E., Lapin A., Pieskä J. Mesh approximation and iterative solution of the continuous casting problem // in: ENUMATH 99. Ed. by P. Neittaanmäki, T. Tiihonean and P. Tarvainen. World Scientific, Singapore. – 2000. – P. 601–617.

[6] Brezzi F., Fortin M. Mixed and hybrid finite element methods. – New-York: Springer Verlag, 1991.

[7] Roberts J.E., Thomas J.M. Mixed and hybrid methods//Numer. Anal. – 1991. – V. II. – P. 523–639.

[8] Rodrigues J.F., Yi F. On a two-phase continuous casting Stefan problem with nonlinear flux//Euro. J. Appl. Math. – 1990. – No 1. – P. 259-278.

[9] Rodrigues J.F. Variational methods in the Stefan problem//Lect. notes in math. – 1994. – Springer Verlag. – P. 149–212.

[10] Lions J.L. Quelques méthodes de résolution des problèmes aux limites non linéaires. – Paris: Dunod Gauthier-Villars, 1969.

[11] Ciarlet P.G., Lions J.L. Handbook of numerical Analysis. Volume II: Finite element methods (Part I). – Amsterdam: Elsevier, 1991. – 928 p.

Kazan State University, Russia E-mail address: alapin@ksu.ru

Received October 15, 2003