- Email: [email protected]
An enthalpy formulation of the conduction problem with phase changes
0045-7949/92 $5.00 + 0.00 0 1992 Pergamon Press Ltd
Compurers & Srrucrures Vol. 44, No. 5. pp. 983-987. 1992 Printed in Great Britain.
AN ENTHALPY FORMULATION OF THE CONDUCTION PROBLEM WITH PHASE CHANGES F. CESARI Department of Mechanical Engineering, University of Bologna, Viale Risorgimento 2, Bologna, Italy (Received 2 July 1991)
Abstract-In this paper a simple and effective algorithm is presented for the resolution of the heat conduction problems in the presence of phase changes. The peculiarity lies in the enthalpy formulation, which is able to treat the phase transition surface, using relatively coarse mesh and large time steps, in a simple manner.
1. INTRODUCTION
This paper deals with heat conduction that involves phase changes; the problem is often called the Stefan problem. The classical formulation of the Stefan problem has the temperature as the unknown variable. In this case the analytical and numerical difficulty lies in that the temperature gradients are discontinuous at the phase transition surface. Analytical, finite differences and finite element methods have been used from several authors, but for practical analysis the use of a finite element discretization appears to be the most attractive [l-3]. The objective of this paper is to present an alternate method for the solution of the Stefan problem. In the procedure, the heat conduction equilibrium equation is formulated with the enthalpy as the primary unknown variable, and the temperature is to be computed from the enthalpy [4]. The advantage lies in that the enthalpy formulation is able to treat the phase transition surface in a very simple and effective method. Relatively coarse mesh and large time steps are required. The differential equation to be solved is summarized and then the Newton-Raphson algorithm is presented to resolve the nonlinear system equation. Some sample solutions demonstrate the effectiveness of the proposed technique.
2. DIFFERENTIAL
U/SD
pici!?.p&g
( > =I),
(1)
where p, , ci , K, are the mass density, specific heat and conductivity coefficient, respectively. The boundary conditions are T,(O, t) = T, < T,
2
(L, t) = 0.
At the interface x = s(t) we find q(s(t),
t) = 7”’
(3)
where 1 is the fusion latent heat per unit mass. Let H be the enthalpy, defined for a pure substance by H(T)
EQUATION
Let us consider the solution of the classical Stefan problem, involving a change phase in the heat conduction of solids. We suppose that in a thin bar of liquid initially the temperature T,, is uniform, To > T,, where Tf is the freezing temperature. At time t = Of the temperature of the left-hand side is reduced to T, < Tf and maintained constant, while the righthand side is insulated. Then the liquid near the left-hand side becomes solid and a phase change front separates the solid and liquid region. Figure 1 shows us
the temperature T versus the coordinate x of the bar, where s(t) is the position of the front. The Fourier differential equation of the liquid (i = 1) and solid (i = s) phase is
= ps c, T, 1 p,c,(T
- Tf) + p,l + p,c, T,,
if T < T, if T > T,.
(5)
The definition of enthalpy enables us to define the enthalpy formulation of the heat conduction. Introducing (5) in (1), we obtain aHi ___ at 983
d
ax
(6)
F.
CESARI
where R is the nodes number of the monodimensional element, while N; are the shape functions and at the same time the test functions. For all test functions N; the weak equation for a element of 1 length is given by
s ‘aH’
o at
0
s(t*)
1(11)
‘a28’N’dx
-
so
ax*
= 0.
(12)
The substitution of eqn (11) into (12) and the integration by part yields the following
L
x
N’dx
Fig. 1. Phase front in a bar.
s I
Since H, K, p, c are discontinuous introduce the /?(H(T)) parameter
N’rN’A’dx
for T = T,,
+
‘B’Wj?‘dx s0
0
= 0,
(13)
where
I
K”H,
if H < p,c,T,
&=K,T,,
ifp,T,
Pscs
B=
(7)
Equation (13) becomes for a single element
CCke+ Kc/P = 0,
[$~~$-$~‘)KT,-PJ], ifH>H,,
where H, = pscs T,+ p,l and fi is a continuous tion of the enthalpy (Fig. 2). Introduce eqn (7) into (6)
aH(T)
--
~*/WW-N =
at
ax*
Now the initial condition
H(x, 0) = and the boundary
HOG’-o
condition
C=
.
>
7”)
(9)
kBdx s0
the conductivity matrix. For a plane problem with a two nodes element we obtain
(2) and for axisymmetric rl,
o
which coordinates
r2
(10) c’
Then the effect of the enthalpy formulation is to transform the equations (l)-(4) into a more compact form of (8)-(lo), with fl continuous in x.
=
1
rl +
r2/3
4 (rl +
r2)/3
1
rl + r2 K’=T
[
Since a*j/ax* may not exist in a classic sense, we need to define a weak solution of problem (8)-(10). In order to use Galerkin formulation of the finite element method, we introduce the semidiscretization of H and fi parameters for a single element
Equation assemblage
1
(6 + r2)/3 r2 +r1/3
-1
1
OF THE COUPLED SYSTEM
(14) for a single element becomes after Cfi+KjI
=0
(15)
“a :i/-/.
N;(x)H;
P$E,T,
“I = p,clTf*~,A
i=l H Hl Be=
i$,
N;(x)/%,
(11)
of nodes
1 .
-1
4. NUMERICAL SOLUTION
3. WEAK FORMULATION
i
capacity matrix
K'=
o
dx’
He=
N’Ndx 0
is the distributed
BP, t) = BoV, > T/l VW, t) =
I
where
func-
becomes
(14)
Fig. 2. Parameter /? versus enthalpy.
An enthalpy formulation of the conduction problem 1
985
T('K)
+h-O.2 *Tiworic l h11.0
242
233
t(ooc) 235 1.60
4.80
3.20
x(cm)
23a
‘
.m
8.m
6.40
.m
.m
1.w
2.40
3.20
4.00
Fig. 3. Solidification of a bar: T = T(f) with ten elements.
Fig. 5. Solidification of a bar: T = T(x) with 40 elements.
a nonlinear system of coupled differential first-order equations. Following the Euler emplicit time integration scheme, eqn (15) is rewritten at time t + At in the form
where P is a diagonal matrix with elements
which represents
c H(t + h) - H(r) h
c K
2
, ifHcH,
1P&s
FL* +Kj(t+h)=O,
”
(16)
if H, < H < HI
dH-
(21)
ifH>H,. where h = At is the time step. Equation (16) can be resolved using the NewtonRaphson algorithm. Let
F(H)=;H(t
+h)+K/?(t+h)-h-.
c H(t)
(17)
Multiplying K and P, we obtain a general nonsymmetric matrix. The solution of eqn (18) proceeds transforming /3 with system (7) and calculating the temperatures at the end of iteration process
c -
At the kth iteration the equilibrium
equation holds
H
ifH
I
if H, < H Q HI (22) aF
FAHk
= -F(Hk-‘)
H-P,~-PJ,T/
(18)
with Hk=Hk-‘+AHk.
(19)
+TJ,
P/Cl
ifH>H,.
The solution of the linear system equation is obtained with a standard subroutine, utilizing partial pivoting. For these reasons the temperature in the code must be expressed in absolute value (K).
The tangent matrix can be written 5. SAMPLE SOLUTIONS
g=; =
1
The solution of two classical test cases analysed by
(20) several authors is presented. The first case is the
+KP,
2.40
T('K)
1.92
1.44
.m
.44 x(cm)
.m
.a0
l.M
2.40
3.7.0
4.m
Fig. 4. Solidification of a bar: T = T(x) with ten elements.
.m
.m
.m
t.w
2.40
3.20
4.00
Fig. 6. Solidification of a bar: phase front in time.
F.
986
0.10
Phme Front(m)
I
0.08
a06 0.04
0.02
.oo 0.04
.M
0.08
Fig. 9. ~lidi~~tion
0.1 2
ais
of corner: phase front
a20 in time
for
x=y.
17, Fig. 7. Mesh for solution of solidification of corner region. 1
solution of the classical Stefan problem given previously, while the second is a bidimensional example, the solidification of corner region.
234367
8
9
10
IPC)
5.1. Solidification of a semifinite slab of liquid The uniform infinite slab of liquid is considered initially at temperature T, = 255.22 K (Fig. 1). At time t = O+ the temperaure of the left-hand side is reduced to T, = 230.20K. Let T,= 255.17K, p,c,= pscs = 1.8 J/&n3 K, p,A = 70.26 Jbrn’, K,= K, = 1.944 W/cm K Different numbers of elements and time steps have been used in the analysis. To choose an indicative time step, we know the eigenvalue of the monodimensional heat conduction problem, 1 = 12kj12, in the case of distributed capacity matrix and with k equal to the thermal diffusivity. Figure 3 shows the ten elements (I = 0.4 cm) solution when using h = 0.2 and 1 set to predict the temperature at x = 1 cm. The smaller time step gives oscillations at the first instants, after the solution is very accurately predicted. These oscillations are due to the latent heat being concentrated at a node which is absorbed or released before the nodal temperature can change. Figure 4, showing the temperature along the bar at three different times, shows these oscillations which disappear when a finer mesh is used (Fig. 5 and 40 elements). Figure 6 shows the freezing front position, using ten elements and h = 1 set, and
I t Fig. 10. Solidification of corner: thermal map at time t =O.l. ._ . . . _ -. ._ _ 40 elements and h = 0.1 sec. The 40-elements results are practical equal to those obtained theoretically. 5.2. Solidification of a corner region The corner of a uniform infinite medium carrying a liquid with initial temperature T, and freezing temperature T, is considered. Figure 7 shows the geometry of the plate with Dimensions 1 x 1 cm and relative spatial discretization with three-node finite elements. In this case the same equation (15) holds, where the distributed capacity matrix and conductivity matrix are 2 1 1 cc=+ 1 I 2 [1 1 2
I
l(%)
-17.6
TPC) -17.6 7,
-17.6
______~C
_____..____ ..T,_______-
-t7Y7 __._______._.__------____
__._- T ,""_"._~_----
-17.6
-17.Q
-16.1
\
J
-17Q -18.2 -t&O
-
-18.1 'I .m
\ .
, 004
,
, 0.08
.
. 0.12
-18.4 ,
. 0.10
.
‘ t(mC) 0.20
Fig. 8. Solidification of corner: T = T(r) at point B.
4.. .M
.m
J .
/
.
J”
.
.40
. .M
.
. .80
XiCi-8)
.
‘
t.m
Fig. 11. Solidification of corner: temperature at time I = 0.1 for x =y.
An enthalpy formulation of the conduction problem
Ke=-&
b: + c:
b,bz+c,cz
b,b2+cIc2
b:+c:
1b,b,+c,c,
1
b,b,+c,c, b,b, + czcg
b,b,+c,c,
b:+c:
. 1
A is the element area b,yz-ys,
b,=yJ-y,,
c,=xj-x~,
b,=y,-yy,
c*=x,-x,,
c, = x2 - x,
and xi, yi the node coordinates.
6. CONCLUSIONS
The previous procedure for the analysis of nonlinear transient problems with phase change effects has been applied to the study of a pure substance, showing simplicity and good performance, both with relatively coarse mesh and large time steps. The study of an alloy with a phase transition temperature interval greater than zero can be carried out in very simple manner, substituting the condition (5) with
if T < Tf
~3 c, T, H(T)
=
981
pscs T/+
~,a + p,c,(T
p,c,T,+
p,l + p,c,AT,+
if T,< T < T,+ AT,
- T/)9
p,c,(T - T,),
The
initial temperature of the plate is 17.61”C and at time t = O+ the temperature of the surface is reduced at T, = -18.33”C. Let T,= - 17.78”C, p,c, = p,c, = 1 J/cm3 “C, Pll= 0.139 J/cm’, K, = K, = 1 W/cm”C. Figures 8-11 give the results obtained for h = 0.02 set: T,,= -
(a) Figure 8 shows the variation temperature of node B (x = y = 0.5 cm) with time. (b) Figure 9 shows the thermal map at time 1 = 0.1 sec. (c) Figure 10 shows the phase front for x = y. (d) Figure 11 shows the temperature along x=y=OScmattimet=O.lsec. These results are very close to the theoretical ones. The same observations regarding the performance of the method can also be seen in the solution of the solidification of the bar of liquid.
if T > T,+ ATf,
where ATf is a given phase change temperature interval. The enthalpy formulation presented in this paper is applicable to nonsymmetric matrices, therefore is not readily incorporated into standard structural finite element programs, but in an ad hoc code. REFERENCES
G. Comini, S. Del Giudice, R. W. Lewis and 0. C. Zienkiewicz, Finite element solution of nonlinear heat conduction problems with special references to phase change. Int. J. Numer. Meth. Engng 8,613-624 (1974). Y. Ichikawa and N. Kikuchi, A one-phase multi-dimensional Stefan problem by variational inequalities. Inr. J. Numer. Merh. Engng 14, 1221-1239 (1979). Z. D. Rolph and K. J. Bathe, An efficient algorithm for analysis of nonlinear heat transfer with phase changes. Inf. J. Num. Merh. Engng 18, 119-134 (1980). R. E. White, An Introduction to the Finite Element Method with Applictions to Nonlinear Problems. John Wiley (1985).
Compurers & Srrucrures Vol. 44, No. 5. pp. 983-987. 1992 Printed in Great Britain.
AN ENTHALPY FORMULATION OF THE CONDUCTION PROBLEM WITH PHASE CHANGES F. CESARI Department of Mechanical Engineering, University of Bologna, Viale Risorgimento 2, Bologna, Italy (Received 2 July 1991)
Abstract-In this paper a simple and effective algorithm is presented for the resolution of the heat conduction problems in the presence of phase changes. The peculiarity lies in the enthalpy formulation, which is able to treat the phase transition surface, using relatively coarse mesh and large time steps, in a simple manner.
1. INTRODUCTION
This paper deals with heat conduction that involves phase changes; the problem is often called the Stefan problem. The classical formulation of the Stefan problem has the temperature as the unknown variable. In this case the analytical and numerical difficulty lies in that the temperature gradients are discontinuous at the phase transition surface. Analytical, finite differences and finite element methods have been used from several authors, but for practical analysis the use of a finite element discretization appears to be the most attractive [l-3]. The objective of this paper is to present an alternate method for the solution of the Stefan problem. In the procedure, the heat conduction equilibrium equation is formulated with the enthalpy as the primary unknown variable, and the temperature is to be computed from the enthalpy [4]. The advantage lies in that the enthalpy formulation is able to treat the phase transition surface in a very simple and effective method. Relatively coarse mesh and large time steps are required. The differential equation to be solved is summarized and then the Newton-Raphson algorithm is presented to resolve the nonlinear system equation. Some sample solutions demonstrate the effectiveness of the proposed technique.
2. DIFFERENTIAL
U/SD
pici!?.p&g
( > =I),
(1)
where p, , ci , K, are the mass density, specific heat and conductivity coefficient, respectively. The boundary conditions are T,(O, t) = T, < T,
2
(L, t) = 0.
At the interface x = s(t) we find q(s(t),
t) = 7”’
(3)
where 1 is the fusion latent heat per unit mass. Let H be the enthalpy, defined for a pure substance by H(T)
EQUATION
Let us consider the solution of the classical Stefan problem, involving a change phase in the heat conduction of solids. We suppose that in a thin bar of liquid initially the temperature T,, is uniform, To > T,, where Tf is the freezing temperature. At time t = Of the temperature of the left-hand side is reduced to T, < Tf and maintained constant, while the righthand side is insulated. Then the liquid near the left-hand side becomes solid and a phase change front separates the solid and liquid region. Figure 1 shows us
the temperature T versus the coordinate x of the bar, where s(t) is the position of the front. The Fourier differential equation of the liquid (i = 1) and solid (i = s) phase is
= ps c, T, 1 p,c,(T
- Tf) + p,l + p,c, T,,
if T < T, if T > T,.
(5)
The definition of enthalpy enables us to define the enthalpy formulation of the heat conduction. Introducing (5) in (1), we obtain aHi ___ at 983
d
ax
(6)
F.
CESARI
where R is the nodes number of the monodimensional element, while N; are the shape functions and at the same time the test functions. For all test functions N; the weak equation for a element of 1 length is given by
s ‘aH’
o at
0
s(t*)
1(11)
‘a28’N’dx
-
so
ax*
= 0.
(12)
The substitution of eqn (11) into (12) and the integration by part yields the following
L
x
N’dx
Fig. 1. Phase front in a bar.
s I
Since H, K, p, c are discontinuous introduce the /?(H(T)) parameter
N’rN’A’dx
for T = T,,
+
‘B’Wj?‘dx s0
0
= 0,
(13)
where
I
K”H,
if H < p,c,T,
&=K,T,,
ifp,T,
Pscs
B=
(7)
Equation (13) becomes for a single element
CCke+ Kc/P = 0,
[$~~$-$~‘)KT,-PJ], ifH>H,,
where H, = pscs T,+ p,l and fi is a continuous tion of the enthalpy (Fig. 2). Introduce eqn (7) into (6)
aH(T)
--
~*/WW-N =
at
ax*
Now the initial condition
H(x, 0) = and the boundary
HOG’-o
condition
C=
.
>
7”)
(9)
kBdx s0
the conductivity matrix. For a plane problem with a two nodes element we obtain
(2) and for axisymmetric rl,
o
which coordinates
r2
(10) c’
Then the effect of the enthalpy formulation is to transform the equations (l)-(4) into a more compact form of (8)-(lo), with fl continuous in x.
=
1
rl +
r2/3
4 (rl +
r2)/3
1
rl + r2 K’=T
[
Since a*j/ax* may not exist in a classic sense, we need to define a weak solution of problem (8)-(10). In order to use Galerkin formulation of the finite element method, we introduce the semidiscretization of H and fi parameters for a single element
Equation assemblage
1
(6 + r2)/3 r2 +r1/3
-1
1
OF THE COUPLED SYSTEM
(14) for a single element becomes after Cfi+KjI
=0
(15)
“a :i/-/.
N;(x)H;
P$E,T,
“I = p,clTf*~,A
i=l H Hl Be=
i$,
N;(x)/%,
(11)
of nodes
1 .
-1
4. NUMERICAL SOLUTION
3. WEAK FORMULATION
i
capacity matrix
K'=
o
dx’
He=
N’Ndx 0
is the distributed
BP, t) = BoV, > T/l VW, t) =
I
where
func-
becomes
(14)
Fig. 2. Parameter /? versus enthalpy.
An enthalpy formulation of the conduction problem 1
985
T('K)
+h-O.2 *Tiworic l h11.0
242
233
t(ooc) 235 1.60
4.80
3.20
x(cm)
23a
‘
.m
8.m
6.40
.m
.m
1.w
2.40
3.20
4.00
Fig. 3. Solidification of a bar: T = T(f) with ten elements.
Fig. 5. Solidification of a bar: T = T(x) with 40 elements.
a nonlinear system of coupled differential first-order equations. Following the Euler emplicit time integration scheme, eqn (15) is rewritten at time t + At in the form
where P is a diagonal matrix with elements
which represents
c H(t + h) - H(r) h
c K
2
, ifHcH,
1P&s
FL* +Kj(t+h)=O,
”
(16)
if H, < H < HI
dH-
(21)
ifH>H,. where h = At is the time step. Equation (16) can be resolved using the NewtonRaphson algorithm. Let
F(H)=;H(t
+h)+K/?(t+h)-h-.
c H(t)
(17)
Multiplying K and P, we obtain a general nonsymmetric matrix. The solution of eqn (18) proceeds transforming /3 with system (7) and calculating the temperatures at the end of iteration process
c -
At the kth iteration the equilibrium
equation holds
H
ifH
I
if H, < H Q HI (22) aF
FAHk
= -F(Hk-‘)
H-P,~-PJ,T/
(18)
with Hk=Hk-‘+AHk.
(19)
+TJ,
P/Cl
ifH>H,.
The solution of the linear system equation is obtained with a standard subroutine, utilizing partial pivoting. For these reasons the temperature in the code must be expressed in absolute value (K).
The tangent matrix can be written 5. SAMPLE SOLUTIONS
g=; =
1
The solution of two classical test cases analysed by
(20) several authors is presented. The first case is the
+KP,
2.40
T('K)
1.92
1.44
.m
.44 x(cm)
.m
.a0
l.M
2.40
3.7.0
4.m
Fig. 4. Solidification of a bar: T = T(x) with ten elements.
.m
.m
.m
t.w
2.40
3.20
4.00
Fig. 6. Solidification of a bar: phase front in time.
F.
986
0.10
Phme Front(m)
I
0.08
a06 0.04
0.02
.oo 0.04
.M
0.08
Fig. 9. ~lidi~~tion
0.1 2
ais
of corner: phase front
a20 in time
for
x=y.
17, Fig. 7. Mesh for solution of solidification of corner region. 1
solution of the classical Stefan problem given previously, while the second is a bidimensional example, the solidification of corner region.
234367
8
9
10
IPC)
5.1. Solidification of a semifinite slab of liquid The uniform infinite slab of liquid is considered initially at temperature T, = 255.22 K (Fig. 1). At time t = O+ the temperaure of the left-hand side is reduced to T, = 230.20K. Let T,= 255.17K, p,c,= pscs = 1.8 J/&n3 K, p,A = 70.26 Jbrn’, K,= K, = 1.944 W/cm K Different numbers of elements and time steps have been used in the analysis. To choose an indicative time step, we know the eigenvalue of the monodimensional heat conduction problem, 1 = 12kj12, in the case of distributed capacity matrix and with k equal to the thermal diffusivity. Figure 3 shows the ten elements (I = 0.4 cm) solution when using h = 0.2 and 1 set to predict the temperature at x = 1 cm. The smaller time step gives oscillations at the first instants, after the solution is very accurately predicted. These oscillations are due to the latent heat being concentrated at a node which is absorbed or released before the nodal temperature can change. Figure 4, showing the temperature along the bar at three different times, shows these oscillations which disappear when a finer mesh is used (Fig. 5 and 40 elements). Figure 6 shows the freezing front position, using ten elements and h = 1 set, and
I t Fig. 10. Solidification of corner: thermal map at time t =O.l. ._ . . . _ -. ._ _ 40 elements and h = 0.1 sec. The 40-elements results are practical equal to those obtained theoretically. 5.2. Solidification of a corner region The corner of a uniform infinite medium carrying a liquid with initial temperature T, and freezing temperature T, is considered. Figure 7 shows the geometry of the plate with Dimensions 1 x 1 cm and relative spatial discretization with three-node finite elements. In this case the same equation (15) holds, where the distributed capacity matrix and conductivity matrix are 2 1 1 cc=+ 1 I 2 [1 1 2
I
l(%)
-17.6
TPC) -17.6 7,
-17.6
______~C
_____..____ ..T,_______-
-t7Y7 __._______._.__------____
__._- T ,""_"._~_----
-17.6
-17.Q
-16.1
\
J
-17Q -18.2 -t&O
-
-18.1 'I .m
\ .
, 004
,
, 0.08
.
. 0.12
-18.4 ,
. 0.10
.
‘ t(mC) 0.20
Fig. 8. Solidification of corner: T = T(r) at point B.
4.. .M
.m
J .
/
.
J”
.
.40
. .M
.
. .80
XiCi-8)
.
‘
t.m
Fig. 11. Solidification of corner: temperature at time I = 0.1 for x =y.
An enthalpy formulation of the conduction problem
Ke=-&
b: + c:
b,bz+c,cz
b,b2+cIc2
b:+c:
1b,b,+c,c,
1
b,b,+c,c, b,b, + czcg
b,b,+c,c,
b:+c:
. 1
A is the element area b,yz-ys,
b,=yJ-y,,
c,=xj-x~,
b,=y,-yy,
c*=x,-x,,
c, = x2 - x,
and xi, yi the node coordinates.
6. CONCLUSIONS
The previous procedure for the analysis of nonlinear transient problems with phase change effects has been applied to the study of a pure substance, showing simplicity and good performance, both with relatively coarse mesh and large time steps. The study of an alloy with a phase transition temperature interval greater than zero can be carried out in very simple manner, substituting the condition (5) with
if T < Tf
~3 c, T, H(T)
=
981
pscs T/+
~,a + p,c,(T
p,c,T,+
p,l + p,c,AT,+
if T,< T < T,+ AT,
- T/)9
p,c,(T - T,),
The
initial temperature of the plate is 17.61”C and at time t = O+ the temperature of the surface is reduced at T, = -18.33”C. Let T,= - 17.78”C, p,c, = p,c, = 1 J/cm3 “C, Pll= 0.139 J/cm’, K, = K, = 1 W/cm”C. Figures 8-11 give the results obtained for h = 0.02 set: T,,= -
(a) Figure 8 shows the variation temperature of node B (x = y = 0.5 cm) with time. (b) Figure 9 shows the thermal map at time 1 = 0.1 sec. (c) Figure 10 shows the phase front for x = y. (d) Figure 11 shows the temperature along x=y=OScmattimet=O.lsec. These results are very close to the theoretical ones. The same observations regarding the performance of the method can also be seen in the solution of the solidification of the bar of liquid.
if T > T,+ ATf,
where ATf is a given phase change temperature interval. The enthalpy formulation presented in this paper is applicable to nonsymmetric matrices, therefore is not readily incorporated into standard structural finite element programs, but in an ad hoc code. REFERENCES
G. Comini, S. Del Giudice, R. W. Lewis and 0. C. Zienkiewicz, Finite element solution of nonlinear heat conduction problems with special references to phase change. Int. J. Numer. Meth. Engng 8,613-624 (1974). Y. Ichikawa and N. Kikuchi, A one-phase multi-dimensional Stefan problem by variational inequalities. Inr. J. Numer. Merh. Engng 14, 1221-1239 (1979). Z. D. Rolph and K. J. Bathe, An efficient algorithm for analysis of nonlinear heat transfer with phase changes. Inf. J. Num. Merh. Engng 18, 119-134 (1980). R. E. White, An Introduction to the Finite Element Method with Applictions to Nonlinear Problems. John Wiley (1985).