Ink J. Hear Mass Tramfer. Vol. 29, No. 2, PP. 285-291, 1986 Printed in Gnat Britain 0017-9310/86S3.00+0.00 Pergamon Press Ltd. The use of lumped capacitance in the finite-element solution of heat conduction problems with phase change Q. T. PHAM Meat Industry Research Institute of New Zealand (Inc.) P.O. Box 617, Hamilton, New Zealand zyxwvutsrqponmlk (Received 22 January 1985 and injnalform 22 July 1985) Abstract-In the finite-element solution of heat conduction problems with phase change, the singular behaviour of the specific heat near the phase-change point causes special problems. These problems can be overcome by lumping the thermal capacitance of the material at the nodes, to obtain a diagonal capacitance matrix. In this way the specific heat can be unambiguously calculated. The lumped-capacitance principle was combined with the explicit enthalpy method and Pham’s three-level enthalpy method. Tests against freezing problems show that for Stefan problems both resulting methods are more accurate than previous three-level distributed capacitance methods. zyxwvutsrqponmlkjihgfedcbaZYXWVUTS INTRODUCTION NUMEROUS methods have been proposed for the numerical solution of heat conduction processes with phase change, as reviewed for example by Fox [l], Furzeland [2] and Crank [3]. Of these methods, many can deal only with situations where a sharp phase- change boundary exists, and thus do not apply to the freezing of many materials of interest such as solutions, foodstuffs or alloys. Others are restricted to one- dimensional problems. This paper will concentrate on the so-called ‘weak solution’ [4] methods, which do not explicitly make use of the phase-change boundary and so are free of the restrictions mentioned above. A single partial differential equation is used to describe the heat transfer process throughout the material (rather than one equation for each phase). The differential equation can be written in one of two ways : c g = div [k grad (T’)] (1) or CYH at = div [k grad (T)]. Equation (1) is the basis of temperature methods, while (2) is the basis of enthalpy methods. Temperature methods Temperature methods have been used in conjunc- tion with both finite-difference [5--73 and finite- element schemes [S, 91. Their major drawback is that the specific heat c appears in the partial differential equation. Near the phase-change temperature, c changes extremely rapidly and may tend towards infinity, and the equation becomes highly non-linear. Pham [lo] reviews various ways proposed to overcome this problem. Of these, the most convenient (because of These modifications prevent jumping of the latent heat peak, but it cannot be said that they have a clear physical basis. Morgan et al. [13] suggest that c should, instead, be calculated from the temperature and enthalpy changes at the previous time levels, m and m+ 1. This constitutes a deviation from the three-level scheme (in which properties should be evaluated at the mid-level), and thus reasonably small time steps must still be used. Cleland et al. [14] use direct numerical integration over each element to calculate the specific heat. This method works only if latent heat release is gradual rather than sharp. its avoidance ofiteration andunconditional stability)is Lees’ [1 l] three-level scheme in which equation (1) is approximated by : c(T”‘+’ - Tm)/2At = {div [k grad (T”+‘)] +div [k grad (T”+‘)] +div [k grad (TM )]}/3 (3) and thermal properties are evaluated at the middle time level. Oscillations are often observed with non-linear problems and so the following ‘damping equation’ [ 121 is usually applied when updating at each time step : T” = (Tm+‘+Tm+‘+T’“)/3. (4) The singular behaviour of c still causes problems, and various other approximations have been resorted to. Comini et al. [8] assume that the enthalpy H follows the same distribution function as the temperature T and calculate the specific heat from : This equation was subsequently modified by Comini and Del Giudice [9] for the two-dimensional situation 285