7 An Implicit Enthalpy Solution

Implicit

enthalpy

solution

for phase change problems:

V. R. Voller

Paper AMM 1065

An implicit enthalpy solution for phase change problems: with application to a binary alloy solidification V. R. Voller Mineral Resources Research Center, University of Minnesota, 56 East River Road, Minneapolis, MN55455, USA (Received January 1986; revised September 1986)

An enthalpy method is developed for analysis of one-dimensional phase change problems under heat conduction. This technique is based on a fixed space grid but a variable time step is used to ensure thatthephasefront isalwaysonanodepoint.TheresuIting numerical solution is implicit in nature and although iteration on the time step is needed at each time level, the tridiagonal matrix algorithm can be utilized within the iterations leading to efficient and accurate solutions. The node-jumping scheme is extended and a numerical solution to the well-known binary solidification model developed. This solution is implicit and performs well in respect to predicting discontinuities in the enthalpy and concentration fields and in predicting smooth non-oscillating temperature and concentration and histories. Keywords: enthalpy method, discontinuities, node-jumping

The purpose

of this paper is twofold:

(1) to introduced an implicit enthalpy scheme for fast and accurate solution to one-dimensional phase change problems (2) to modify this scheme and apply it to the solution of a binary alloy solidification problem Enthalpy methods have been a popular means of numerically solving phase change problems for some time.1’2 A major reason for this is that an enthalpy formulation removes the need to satisfy conditions on the moving boundary, which means that fixed domain solutions can be used. Direct applications of numerical enthalpy methods, however, may lead to results which have non-physical oscillations. Such behaviour has been noted and analysed by a number of authors.3-6 In order to overcome, or bypass, these oscillations remedial schemes have to be applied.2.7~s One of these is the so-

phase change,

binary

solidification,

called node-jumping scheme in which the space grid remains fixed but the time step changes ensuring that the phase change boundary is always on a node point. This scheme requires the application of an implicit time representation of the enthalpy formulation. Recently, Voller,’ developed an implicit enthalpy method which dovetails very neatly into the node-jumping technique. It is this combination which is presented and developed in this paper. In the first instance, the new scheme will be demonstrated on a simple test problem. After this, the scheme will be modified to deal with a problem of a binary alloy solidification, a problem which involves a phase change along with coupled heat and mass transfer. The major advantage of the node-jumping scheme when applied to this problem is that in the numerical scheme there is no need to account for coupling between enthalpy and concentration as is the case with previous numerical schemes. l&l2 0307-904X/87/02110-07/$03.00

110

Appl. Math. Modelling,

1987, Vol. 11, April

0

1987 Butterworth

& Co. (Publishers)

Ltd

Implicit

enthalpy

solution

for phase change problems:

where h = CT is the ‘sensible enthalpy’, is the diffusivity, and:

V. R. Voller K

_,I 3 = d(AH)

is a latent heat source term. The term AH in equation (4) is the latent heat evolved in the phase change. The functional relationship of AH with temperature, T, depends on the nature of the phase change. In the case of a single temperature phase change, refer to equation (2), then:

Control volumes

Liquid

Sol id

Liquid

AH=

0 N+l

a

AHN’

L/2

Solid

b Figure2

AHN=

0

T< T,,,

L

T> T,,,

The alternative enthalpy formulation given by equation (3) has proved to be useful in a number of respects. In particular, in seeking implicit numerical solutions the nonlinearity associated with the latent heat can be effectively dealt with via the source term, S.g

AHN+~ = L

Liquid

0

K/PC

(4)

at Figure 7

=

Node jumping scheme The scheme In the majority of physical situations, it is necessary to solve equation (3) via a numerical method. For a one-dimensional problem an appropriate method is that of finite difference control volumes. The region of interest is subdivided into a set of discrete nodes. With each node there is associated a volume such that values of variables at the node are representative of the values in the volume. Three adjacent control volumes are shown in Figure 1. From this figure, the following finite difference representation of equation (3) may be derived:

AHN+~ = L/2

One node jump: (a) time t; (b) time (t + ST)

The enthalpy formulation Enthalpy, H, is defined as the sum of sensible and latent heats. For a melting/freezing problem under conduction, the enthalpy formulation is usually written as: where: p%

= V[K(VT)]

(1)

FP =

[Kdkv - hp)- K_(hp-

Ml

(7)

Sp=AG-AH;

where T is temperature, K conductivity and p density. The relationship between temperatures and enthalpy will depend on the nature of the phase change. For example for a pure material undergoing a melting/freezing phase change, all the associated latent heat is evolved at a fixed temperature, T = T,,, and the relationship is: H=

CT

T< T,,,

CT-I-L

T> T,,,

(2)

where C is specified heat and L is the latent heat. Vollerg has proposed an alternative formulation to that given in equation (1). Essentially, the latent heat and sensible heat terms are separated and equation (1) written as:

$ =V(KVh)

+

s

(3)

(8) Note that ( ) * indicates evaluation at the previous time level, upper case subscripts, e.g. P(oint) and W(est), indicate evaluation at a node point and lower case, e.g. w(est) and e(ast), indicate evaluation at a control volume face. The parameter 8 can take any value between 0 (fully explicit) and 1 (fully implicit). For a phase change problem, Voller9 has recently developed an efficient and stable implicit solution of equation (6) based on a fixed grid with fixed time step. Central to this solution technique is the interpretation of the discretized source term of equation (8). The value Sp gives the change in the nodal latent heat AH, over one time step and this change will be proportional to the amount of control volume which has changed state. In the case of a single temperature phase change in a one-dimensional region, this interpretation can be coupled with equation (6) to provide a highly efficient solution scheme. Assume that at time t the phase front is on node point N (see Figure 2(u)). Then half the control volume will be solid and have a latent heat value of 0. The other

Appl.Math.

Modelling,

1987,Vol.

11,April

111

Implicit

enthalpy

solution

for phase change problems:

V. R. Voller

Table 1 Node-jumping results Positionof

boundary

Temperature at0.5m

Time(h)

Numerical

True

Numerical

True

12.97 52.04 117.18 208.48 325.90 469.42 639.05 834.78

0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.00

0.124632 0.249622 0.374573 0.499617 0.624673 0.749707 0.874734 0.99976

1.79821 0.87151 0.32077 0.0 -1.89248 -3.19272 -4.1371 -4.85315

1.79216 0.886201 0.325863 0.00077 -1.89781 -3.20379 -4.15167 -4.86946

half will be liquid with a latent heat value of L. Hence, the representative nodal latent heat will be: AH; = L/2 Now, assume that the numerical time step t is such that at time t + 6t the phase front is located on the N + lth node (see Figure2(b)). Hence: AH,,,+, = L/2 and the change in latent heats in the time step in control volumes N and N + 1 will be: S, = AH; - AH, = L/2 and: S ,,,+, = AH;+I - AHN+, = L/2 If the time step was always chosen to ensure that the phase front always ‘jumped’ from one node to the next in a time step the form of the source term in equation (6) would always be known in advance. For example, if the phase front was on node N at the start of the time step and node N + 1 at the end: Sp =

L/2 [0

ifP=

NorP=

N+ 1

otherwise

This fact ensures that the system of equations resulting from equation 6 will be linear and tridiagonal and hence can be easily and efficiently solved by the tridiagonal matrix algorithm. The only unknown in the system is the time step required to jump from node to node. The value of the time step can be approximated using an iterative procedure, with the time step modified from one iteration to the next. The point of convergence for such a procedure being reached when the predicted value of h at node N + 1 is: h N+l

=

CT,

The node-jumping follows:

procedure

may be summarized

as

(1) phase front on node N (2) initial guess for time step (at), = previous value (3) set SN and SN+i = L/2, all other Si to zero (4) solve equation (6) for nodal distribution of sensible heat (5) update value of time step (via a discretized Newton method on the function hN+I - CT,,, = 0, for example) (6) repeat steps (4) and (5) until convergence, e.g. until [CT,,, - hN+,l < TOL.

112

Appl.

Math. Modelling, 1987,Vol.ll,April

Example problem A material in a liquid state and uniform temperature Ti = 2 > T,,, = 0 (the phase change temperature) fills the half space x > 0. At time t = 0, the temperature of the surface at x = 0 is lowered to T, = -10 so that, as time advances, a solid layer builds up adjacent to the wall. The material has the following thermal properties, conductivity K = 2, specific heat C = 2.5 x lo6 and latent heat L = 1O8. The solution of this problem via the node-jumping algorithm is compared with the known analytical solution in Table 1. Note that the step size used in these results was 6x = 0.0625 and for the first jump from x = 0 to x = 0.0625 a fully implicit mode (f3 = 1) was used, thereafter a Crank-Nicolson mode (0 = f) was adopted. Both the results for the movement of the boundary and temperature, given in Table 1, show that the method is accurate. In addition, the technique is very efficient, requiring somewhere in the order of two to five iterations at each time step. This means that accuracy can be taken to a fairly high level simply by using more nodes. It can be argued that the problem solved in Table I does not present a very severe test on the method. In the next section, however, the same node-jumping principal will be used and the technique modified to solve a binary alloy problem, a more complex problem of coupled heat and mass transfer.

Binary alloy problem The problem The binary alloy problem is that of solidification of a metal containing a small concentration of impurity (solute). The model of the system used here is identical to that proposed by Crowley and Ockendon.” Modification and application of the node-jumping scheme, however, will result in a more accurate and efficient solution than those previously presented. The developed solution will also be, to the author’s knowledge, the first implicit technique proposed for general treatment of the binary alloy problem (note the implicit scheme proposed by Whitei requires that the phase change temperature T, remains fixed; in general this will not be the case). An equilibrium phase diagram for a dilute binary alloy is shown in Figure 3 where the liquidus and solidus lines are assumed to be straight. Consider a small volume element of binary alloy in thermodynamic equilibrium at uniform temperature T and overall average concen-

Implicit

enthalpy

tration of solute c. If conditions are such that c + kzT > 0 (k, = (liquidus slope))‘) then the volume is in the stable liquid state. If c + k,T< 0 (k, = (solidus slope)-‘) then the element is in a stable solid state. When the values of c and T are such that they correspond to a point between the solidus and liquidus lines, the volume is undergoing the phase change and both solid and liquid phases, at different concentrations of solute, will be present separated by a plane front. If equilibrium is maintained, a discontinuity in concentration at the phase front will exist such that: T + q/k,

solution

for phase change problems:

In a similar onnodeN+

manner, 1:

after the jump

CN+I = (k, + M/2

V. R. Voller

when the front

Vv+,

These observations can be utilized in Crank-Nicolson finite difference representations of equations (10) and (11) to give for heat: hp-h;=;F;+;Fp+Sp

(14)

where A = &/6x2 and FP and Sp are defined (7) and (8), and for concentration:

= T + c2/k2 = 0

is

in equations

i.e.:

Cl

=

E,V,-E;V,“=;Gb+;G,

k, k2

This discontinuity makes the concentration similar to the enthalpy function H( 7) previously introduced. The enthalpy is a function of temperature T which is continuous at the boundary. Concentration can be made a function of the chemical activity, V,” defined by:

where GP = [D,.,( VW - VP) - D,(V, Under a node-jumping regime: Sp=

where i = 1, solid and i = 2, liquid. In this case at the interface, since ci = klcz/k2, the chemical activity is continuous, i.e. VI = V,. Assuming heat and mass transfer in one direction, the governing conservation equations are:

EP =

E; = (10)

(11) where D = (mass diffusion the relationships:

coefficient)

ki, along with

T+ V-CO

[h, h + L]H

T+V=O

h+L

T+ V>O

(12)

T+ V
k,Vc [k, V,

x

&VI

kzV

T+V=O T+ V>O

(13)

Application of node-jumping The basic principal involved in node-jumping is as already outlined. That is, the time step is chosen so that the phase front is always on a numerical node point. The treatment of the heat transfer is as before. As the phase front moves from node N to node N + 1 the latent heat lost from control volumes N and N + 1 is equal to L/2, with other latent heat losses equal to zero. For the mass transfer, the following observations are made. When the front is on node N the associated control volume is half solid and half liquid, hence the representative concentration is: c; = (k, + k,)/2

V;

L/2 [0

ifP=

k2

ci = k;V;

h

(15)

(9)

-c-2

(k,

NorP=

- V,)]

N-t 1

otherwise

+ k,)P

1

(16)

ifP>N+l ifP=N+l

kl

ifP
kz

ifP>N

(k2 + k,)/2

ifP=N

k,

ifP
(17a)

(17b)

When boundary conditions have been incorporated into equations (14) and (15), the result is two tridiagonal systems which can be solved when the time step is given. In order to close the solution, the correct time step (i.e. the time step which will accomplish the jump) needs to be found. As before, the correct time step can be found via an iterative technique, with the point of convergence reached when TN+, + V,,, = 0. The main steps then in applying the node-jumping technique to the binary alloy problem are: (1) (2) (3) (4) (5) (6) (7)

phase front on node N set (St), = previous value set S,v and SN+, = L/2 set EN+~ and Ei = (k, + k,)/2 solve equations (14) and (15) for h and V update time step repeat steps (5) and (6) until convergence, I%+, + V,v+,I
e.g.

An immediate advantage of the node-jumping solution of the binary alloy problem is the fact that there is no explicit numerical coupling required. In previous numerical techniques, l’ the heat and mass solutions are coupled. To achieve this, an ancillary relationship between temperature and concentration is required. In the node-jumping approach, the nodal latent heat values, equation (16) and the relationship between the nodal concentration, cp, and chemical activity VP, equations (17), are known a priori and no explicit coupling relationship is needed. The consistency between the heat and mass fields is maintained by ensuring that the

~ppl.

Math. Modelling,

1987,Vol.

1 l,ApriI 113

Implicit

enthalpy

solution

for phase change problems:

V. R. Voller 0.15

0.125 -

0.025 -

0 figure3

Equilibrium

diagram

I 0.2 X

I 0.1

0

Figure 6 4)

Concentration

I 0.3

0.4

profile at r = 0.02 (key as for Figure

2.0L

1.0 s

3

0

0.1

I 0.3

history at x = 0.1; (-

o 0.075 -

1analytic; (0)

Figure 4)

7

Enthalpy

profile at time r = 0.02 (key as for Figure

l

0.05 -

0.025 0

0 Figure 5

114

1

I

1

0.01

0.02

0.03

t

Cancellation history at X = 0.1 (key as for Figure 4)

Appl.

Math. Modelling,

0.4

X

t Figure 4 Temperature node jumping

I 0.2

1987,Vol.

11,April

Figure 8

Phase front movement

(key as for Figure41

implicit

enthalpy

final iterative time step is such that the sum of the predicted values in the control volume changing phase are close to zero (i.e. TN+, + V,,, = 0) which is the requirement of the heat mass coupling. Another feature of the node-jumping solution is that it is implicit. Previous techniques have, in the main, been explicit so that the coupling between temperature and concentration can be easily incorporated. The only other implicit technique known to the author, is that due to White.‘” This technique, however, appears to assume that the phase change occurs at a fixed value of temperature T,,, and chemical activity V,,,. This removes the need for a coupling relationship because the nature of a control volume, i.e. solid or liquid can be determined purely on the temperature value. Assuming a fixed temperature may not always be sound and some problems may exhibit variations in the phase change temperature and concentration as the phase change proceeds. The node-jumping technique determines the nature of a control volume from the sum T + V, an approach that will allow for variations.

An application To illustrate the application of the node-jumping binary alloy solution the test problem introduced by Crowley and Ockendon” is used. A binary mixture with an equilibrium diagram as in Figure 3, initially at uniform temperature Ti = 1 and concentration c, = 0.1 fills the half space x > 0. At time t = 0 the surface temperature is lowered to T, = -1, so that phase change can commence. It is assumed that there is no mass flow across the surface x = 0 and the properties of the alloy are: Kr = K,=

1

kZ = 2

k, = 1

Specific heat, C = 1

L=l

D,=D2=1 The latter expression implies that the mass diffusivities take the values 1 for the solid and f for the liquid. Rubinsteini4 has proposed an analytical similarity solution for the binary alloy problem. In certain cases, however,i5 the implicit assumption made by Rubinstein of a plane front breaks down and a mushy region occurs. For the test problem used here, no such breakdown occurs and the Rubinstein analytical solution reported by Crowley and Ockendon” can be used for comparison purposes. Solutions of the problem using a step of 6x = 0.0125 are compared in Figures 4-8. Figure 4 shows the temperature history of the point x = 0.1. The agreement between the analytical and numerical solutions is to within three significant figures. Furthermore, the oscillations noticeable in the results of Crowley and Ockendon (Figure 2(a) of reference 11) are not present in the nodejumping solution. Note, however, that in some later employing the remedial scheme of work, Crowley,ih Tacke,” has also removed the oscillations. Figure 5 gives the concentration history at x = 0.1. Again the numerical results are very close to the analytical solution and free of oscillations. The numerical values of the concentration may be manipulated to predict the jump discontinuity in the concentration history. The jump discontinuity will occur when the front is on the node located at x = 0.1 (node N say). At this time the nodal concentration at x = 0.1 is representative of

solution

for phase change problems:

the surrounding control volume, solid and half liquid, hence:

a volume

V. R. Voller which is half

c/v = (k, + k&v/2

(18)

The value cN is known, therefore V,,, can be calculated and the liquid and solid concentrations on either side of the front calculated as: c, = k,V,

cZ = k2V,\z

(19)

respectively. Both of these values, as derived from the numerical concentration when the front is on x = 0.1, are shown in Figure 5. The sharp manner in which the jump discontinuity is modelled by the node-jumping results is a notable feature. In the Crowley and Ockendon results (Figure 2(b) of reference 11) the numerical concentration history indicates a smearing of the discontinuity. Figures 6 and 7 show the enthalpy and concentration profiles at time t = 0.02. Once again, these results have high accuracy and discontinuities are well modelled by similar relationships to those given in equations (18)-(19). Additional information which is a natural consequence of the node-jumping method is the movement of the front. Numerical predictions are compared with analytical ones in Figure 8. As with all other results, the accuacy is good. Furthermore, note that other enthalpy-based methods will not provide a front position prediction quite so readily.

Conclusions In this paper, the enthalpy node-jumping technique for the solution of one-dimensional phase change problems first proposed by Voller and Cross’,’ has been updated. The major modification has been to incorporate the implicit numerical enthalpy technique reported by Voller.” This technique dovetails into the node-jumping concept and results in an efficient solution procedure for heat conduction problems which can utilize the tridiagonal matrix algorithm and requires few iterations per time step. The basic node-jumping technique has been extended and applied to a problem of binary alloy solidification. A problem which involves both heat and mass transfer. Enthalpy techniques have been applied to this problem before. In the main, however. these techniques are explicit in nature and require an additional relationship to account for the coupling between the temperature and concentration,” or if the technique is implicit not always physical assumptions are required.‘” The nodejumping technique presented here requires no additional relationships in that the numerical equations are solved in an uncoupled form. The heat/mass coupling is then accounted for by choosing a time step which moves the tront exactly one node space. In addition, the technique is implicit, yet still applicable to a general binary alloy problem in which the phase change temperature and concentration varies as the phase change proceeds. Results using the node-jumping scheme appear to be more accurate than previous methods retaining the expected physical behaviour especially with reference to discontinuities and temperature history. The efficiency of the technique also means that accuracy may be achieved on adding more node points to the solution.

Appl. Math. Modelling,

1987 Vol. 11,April

115

Implicit

enthalpy

solution

for phase change problems:

Extension of the application of the technique to problems which involve properties changing across the phase front should not present many problems. The major modification would be in the development of the finite difference equations, equations (6), (14)) and (15), and not in the solution technique itself. Problems in more dimensions or non-planar (i.e. mushy region) problems cannot currently be solved by the node-jumping approach. If a multi-dimensional problem is to be solved, then the shape of the phase change surface needs to be known a priori, so that the appropriate solid and liquid fractions in a control volume can be calculated. In the case of a mushy region, phase change difficulties arise because a control volume cannot be separated into two distinct phases (i.e. solid and liquid) which is a current reqmrement of node jumping.In fact, it should be noted that the node-jumping technique is only suitable for the binary alloy problem presented here in the circumstances that it retains a plane front and does not exhibit the mushy region behaviour noted by Wilson et al.15 Clearly more work is required to approach multi-dimensional and mush region problems via an enthalpy node-jumping or similar technique. In conclusion, the node-jumping enthalpy method presented in this paper should provide a tool for efficient analysis of one-dimensional heat conduction phase change problems. Furthermore, the success of the modification and application to the problem of binary alloy solidification indicates that the principle of node-jumping can cope with fairly sophisticated phase change problems. References 1 Crank, J. Free and moving Oxford, 1984, pp.

116

boundary

Appl. Math. Modelling,

problems,

1987,Vol.

Clarendon Press,

11, April

V. R. Voller 2 Voller, V.R. and Cross, M. Application of control volume enthalpy methods in the solution of Stefan problems, in Computational techniques in heat transfer, (eds. Lewis, R. W., et al.), Pineridge, Swansea, 1985 3 Price, P.H. and Slack, M.R. The effect of latent heat on numerical solutions of the heat flow equation, Br. J. Appl. Phys., 1954, 5, 245-2 75 4 Bonacina, C., Comini, G., Fasano, A. and Primicerio, M. Numerical solution of phase change problems, In?. J. Heat Mass Transfer, 1973, I&1825-I832 5 Voller, V.R., Cross, M. and Walton, P.G. Assessment of weak solution numerical techniques for solving Stefan problems, Numerical methods in thermal problems (eds. Lewis, R. W. and Morgan, K.), Pineridae Press. Swansea. 1979. pp. 172-181 6 Belly G.E. On the performance of the enthalpy method, Int. J. Heat Mass Transfer, 1982, 25,587-589 7 Voller, V.R. and Cross, M. Accurate solutions of moving boundary problems using the enthalpy method, ht. J. Heat Mass Transfer, 1981, 24,545-556 8 Tacke, K. Discretization of the enthalpy method for planar phase change, in Free boundary problems: application and theory, Vol. IV (eds. Bossavit, Damlamian and Fremond), Pitman, London, 1985 9 Voller, V.R. Implicit finite-difference solutions of the enthalpy formulation of Stefan problems, ZMA .I. Numer. Anal., 1985, 5, 201-214 10 Fix, G.J. Numerical methods for alloy solidification problems, in Moving boundary problems (eds. Wilson, D.G., et al.), Academic Press, New York, 1978 11 Crowley, A.B. and Ockendon, J.R. On the numerical solution of an alloy solidification problem, In?. J. Heat Mass Transfer, 1979, 22,941-947 12 Wilson, D.G., Solomon, A.D. andAlexiades, V. Amodel of binary alloy solidification, Int. J. Numer. Meth. Engin., 1984, 20,1067-1084 13 White, R.E. The binary alloy problem: existence, uniqueness and numerical approximation, SIAM J. Numer. Anal., 1985, 22, 205-244 14 Rubinstein, L.I. The Stefan problem, Trans. of Math. Monogr., 21,Am. Math. Sot., 1971 15 Wilson, D.G., Solomon, A.D. and Alexiades, V. A shortcoming of the explicit solution for the binary alloy problem, Left. in Heat Mass Transfer, 1982, 9,421-428 16 Crowley, A.B. Numerical solution of alloy solidification problem revisited, in Free Boundary Problems: application and theory, Vol. IV (eds. Bossavit, Damlamian and Fremond), Pitman, London, 1985