Application of Heat Balance Integral Methods to One-Dimensional Phase Change Problems

We employ the Heat Balance Integral Method (cid:3) HBIM (cid:4) to solve a number of thermal and phase change problems which occur in a multitude of industrial contexts. As part of our analysis, we propose a new error measure for the HBIM that combines the least-squares error with a boundary immobilisation method. We describe how to determine this error for three basic thermal problems and show how it can be used to determine an optimal heat balance formulation. We then show how the HBIM, with the new error measure, may be used to approximate the solution of an aircraft deicing problem. Finally we apply the new method to two industrially important phase change problems.


Introduction
Stefan problems arise in numerous industrial applications, such as the manufacture of steel, ablation of heat shields, contact melting in thermal storage systems, ice accretion on aircraft, and evaporation of water 1-5 . However, despite the plethora of applications, there remain very few analytical solutions to Stefan problems and in fact only one of these, the Neumann solution, is of practical use 1 . Consequently, solution methods tend to be numerical or approximate. In this paper we will focus on a particular form of approximate solution, namely, the heat balance integral method HBIM , and apply it to a problem related to aircraft deicing and then to two generic Stefan problems.
The heat balance integral method HBIM is a simple approximate technique originally developed for analysing thermal problems.  was the first to introduce the method, which was adapted from the Karman-Pohlhausen integral method 9 for analysing boundary layers, see 10 for a translated account of this work. However, since exact solutions exist for most standard fixed domain thermal problems, the HBIM has found its greatest use in Stefan problems, see 11, 12 for example. It has also been applied to problems such as thermal explosions, the Korteweg-de-Vries equation, microwave heating of grain, rewetting of surfaces, and boundary layer flows 13-17 . The standard heat balance integral method approximates solutions to the heat equation by first introducing a heat penetration depth, δ t , where for x ≥ δ the temperature change above the initial temperature is assumed to be negligible. An approximating function is then defined for the temperature, typically a polynomial, and by applying sufficient boundary conditions, all the unknown coefficients can be determined in terms of the unknown function δ. Finally, the governing heat equation is integrated for x ∈ 0, δ to produce a heat balance integral, leading to an ordinary differential equation for δ. The problem is therefore reduced to solving a first-order differential equation in time. An alternative approach to the HBIM is the refined integral method RIM , where the heat equation is integrated twice this was termed the RIM in 18 , although it had at least two previous titles, see 1 e.g. . A comparison of the two approaches, along with variations involving alternative approximating functions, is described in detail in Mitchell and Myers 12 . The HBIM has a number of well-known drawbacks. For example, it can only be applied to one-dimensional problems and the choice of approximating function is arbitrary, which is key to the method's accuracy. It is standard to use a polynomial function but even then there is debate over the power of the highest order term 12 . Recently, Mitchell and Myers 19,20 have developed a method where the exponent is determined during the solution process, producing significantly better results than standard models. It involves a combination of the conventional heat balance and refined integral methods, and is called the combined integral method CIM . The method has been applied to both thermal and Stefan problems 19, 20 . Another drawback is that without knowledge of an exact or numerical solution, there is no agreed method for measuring the accuracy. This question was addressed some time ago by Langford 21 who proposed using the least squares error when the approximate solution is substituted into the heat equation. This criterion was subsequently used to improve the accuracy of heat balance methods for thermal and Stefan problems in 22, 23 . Unfortunately, Langford's criterion is not a good measure of the error. For the most basic thermal problem, with a fixed temperature boundary condition, the error depends on time, E L ∝ t −3/2 . Despite the fact that the HBIM approximation may appear reasonable at small times, this measure shows that E L → ∞ as t → 0.
Indeed, Myers 22, 23 minimised Langford's error function to determine the exponent in the approximating function. In 22 the HBIM and RIM formulations were applied to the three standard fixed domain problems and in 23 the same formulations were applied to the classic Stefan problem. This method has the advantage that it is as easy to apply as Goodman's original method, since it simply involves using one of the exponents provided in 22, 23 , it can significantly improve accuracy and in general it is even more accurate than the CIM. However, there are situations where the exponent depends on time, for example with a Newton cooling boundary condition or a time dependent boundary temperature. In these cases, taking the initial value of the exponent will improve on Goodman's result but in general the CIM will provide a better approximation over large times. The relative accuracy of the various methods can be judged through the results provided in the following sections.
In this paper we will analyse various practical problems and along the way propose a new error measure. In the following section we set out the basic thermal problem and use it to illustrate Langford's error before going on to describe our new technique to determine International Journal of Differential Equations 3 a criteria which does not blow up as t → 0. In Section 3 different formulations that arise in the HBIM and the merits and weaknesses of each are discussed. We then examine constant flux and Newton cooling boundary conditions in Section 4. We show that for thermal problems, with either a fixed temperature or constant flux boundary condition, the new error definition is constant, whilst for a Newton cooling condition, the error is zero at t 0. In Section 5 we demonstrate how to apply the technique to a model for aircraft deicing systems. In Section 6 we consider two classic one-phase Stefan problems: the first involves a solid at the solidus melting due to a prescribed temperature at the boundary x 0, the second involves a semiinfinite subcooled solid acting to freeze the adjoining liquid layer. These two generic problems have been applied to diverse situations such as melting or production of ice, water filtration, freezing and thawing of foods, ice formation on pipe surfaces, solidification of metals and magma solidification 2, 5, 24-27 . In both cases the new error is shown to be constant, whereas Langford's definition is again singular at t 0.

Thermal Problem with Fixed Boundary Temperature
We begin by illustrating the HBIM on a simple nondimensional thermal problem described by ∂u ∂t which has the exact solution To calculate the HBIM solution, we replace the boundary condition at infinity with where δ is sufficiently far from the boundary that the boundary temperature has a negligible effect. Two boundary conditions are required to replace the single condition at infinity due to the introduction of a new unknown, the position δ t , where δ 0 0. The appropriate approximating polynomial for u x, t is

4 International Journal of Differential Equations
The HBIM then proceeds by integrating the heat equation over x ∈ 0, δ , Substituting for u x, t , via 2.4 , and integrating lead to δ 2n n 1 t.

2.6
With traditional heat balance methods, the value of n is specified at the outset, typically n 2 or 3, and so the solution is given by 2.4 and 2.6 and the analysis is complete. Now, consider the function f x, t defined by The HBIM formulation and the least-squares error of Langford may then be written as If u x, t is an exact solution of the heat equation then both f and E L will be identically zero. However, with the HBIM the heat equation is only satisfied in an integral sense and f x, t / 0 for all x. That is, the heat equation is not satisfied by the polynomial approximation u x, t and the error E L > 0. This may be seen more clearly by considering the form of these functions. Noting the derivatives from 2.4 ∂u ∂t and E L 2n 4 − 7n 3 6n 2 2n − 1 2n n 1 4 n 1 2 2n − 3 2n − 1 2n 1 t −3/2 e L t −3/2 .

2.11
The above expressions highlight some of the problems inherent in applying heat balance methods. For example, at x 0, International Journal of Differential Equations For u to satisfy the second boundary condition in 2.3 we require n > 1 and so 2.12 indicates solutions of the form 2.4 will never satisfy the heat equation at x 0. At x δ, we find

2.13
With the restriction n > 1, this shows that the heat equation is only satisfied at x δ when n > 2. For n < 2 and n / 0, 1 , and 2.13 shows f δ, t → ∞, with n 2 we find f δ, t → −2/δ 2 . The expressions for E L and e L were first derived by Langford 21 . In Figure 1 we plot e L n as specified by 2.11 and also the equivalent curve for the RIM formulation. The expression in 2.11 is only valid for n > 3/2, below this the expression can predict unrealistic negative and zero values indicating that the integral should be evaluated more carefully in this region. However, provided we specify n > 3/2, the expression 2.11 is positive and there is no value of n that can make E L 0. Of course it should be expected that an approximate solution has some error, so perhaps the most worrying aspect of the expression for E L is that it is time-dependent and more specifically, as t → 0, the error E L → ∞ even though the approximate solution may appear accurate. Hence, while E L may be useful to represent relative errors, that is, for a given time, we can generate E L for different values of n to determine which leads to the most accurate approximation; it does not provide a meaningful measure of the error and is therefore not an appropriate parameter to quantify the error.

A New Error Measure
To define a new error measure that avoids the unrealistic behaviour of E L , we use a boundary immobilisation technique. For thermal problems this simply means that we introduce a new coordinate ξ x δ t , 2.14 which transforms the moving domain 0 < x < δ t to a fixed one 0 < ξ < 1. If we denote u x, t U ξ, t the heat equation 2.1a becomes and the boundary conditions are The polynomial approximation, 2.4 , is now U 1 − ξ n . A significant feature of this profile is that it is independent of δ, that is, the solution has a form that is constant in time. This is a consequence of the fact that the boundary transformation ξ x/δ t corresponds to the similarity solution transformation ξ x/ √ t. Equivalent to deriving equation 2.7 , we rearrange the heat equation, 2.15 , to define the function after noting U t 0. The HBIM formulation may then be written 1 0 F ξ, t dξ, which leads to δ 2n n 1 t, 2.19 and 2.6 is retrieved. We now propose that a new, better behaved error estimate may be obtained from the least-squares error of the boundary immobilised equation: Since δδ t n n 1 is independent of time this function clearly only varies with the value of n. Evaluating the integral then leads to That is, for this problem, we have now defined an error that is independent of time and, in particular, is not singular at t 0. It is related to the expression of Langford by E M δ 3 E L .
International Journal of Differential Equations 7 The RIM simply involves integrating the heat equation twice where η is a dummy variable. After changing the order of integration this becomes 1 0 1 − ξ F ξ, t dξ 0 and then the HBIM formulation 1 0 F ξ, t dξ 0 is substituted to give 1 0 ξF ξ, t dξ 0. This leads to δ n 1 n 2 t.

2.23
It should be noted that if the formulation of 2.22 is used without applying the HBIM integral then the method is referred to as the Alternative Refined Integral Method ARIM .
If we follow this approach to determine the error, we find which is once again independent of time. Now that we have a new error definition we will briefly outline some of the different ways to tackle heat balance problems and then compare the errors predicted by E M .

Defining the Exponent n
The classical HBIM involves specifying n at the start of a calculation. However, there exist a number of refinements which lead to different values of n. These can be categorised loosely as global or local matching methods and obviously they lead to various levels of accuracy.

Local Matching
For this method n is determined through an additional boundary condition, usually taken from the exact solution. For example, in the current problem the boundary condition specifies u 0, t and so, from the exact solution, one can also specify u x 0, t to determine n. If we equate the derivatives of the exact and HBIM expressions for u x, t and set x 0, we find Substituting for δ through 2.6 then shows n 2/ π − 2 ≈ 1.752. This is the approach suggested by Braga et al. 28 . The entropy minimisation method of Hristov 11 involves matching exact and approximate expressions for u 2 x /u 2 at x 0. For a fixed boundary temperature, this is equivalent to matching u x and so Hristov obtains the same n value.
There are two obvious drawbacks to this type of approach. Firstly, they require an exact solution to determine the extra condition or the series expansion, in which case the HBIM is redundant. Secondly, however the matching is carried out, matching a value of

Global Matching
The method described in 22 involves leaving n unknown throughout the calculation. It is then determined by minimising E L over n. This forces the solution to be close to the true solution over the whole domain, rather than at a single point. If we consider the example given above, then in the original error definition n is simply chosen by minimising e L n in 2.11 . Clearly this method may be easily modified to the present situation by choosing n to minimise E M n . Mitchell and Myers 19 take a similar approach to that of 22 in that they leave the exponent unknown and determine it as part of the solution process. In this case their additional equation comes from applying the RIM in addition to the HBIM, hence their method is termed the Combined Integral Method CIM .
In general the method of 22 will provide the most accurate solutions as measured by E L . This is perhaps clearest in Section 4.2, where the minimisation technique has an error half that of the CIM. However, in cases where n n t , due to the difficulty of implementation, the minimisation technique is used with n n 0 . The combined method deals more easily with n t and gives more accurate solutions in this situation.

Comparison of Results
In Figure 2 we present the variation of E M for the HBIM and RIM formulations for n ∈ 1.7, 2.5 . For larger n values both functions increase monotonically, so the only physically realistic minimum is shown on the graph. With the minimisation technique of 22 , this would then indicate that the most accurate approximation comes from using the RIM formulation with n 2.074, with an error E M 0.784. The most accurate HBIM approximation requires n 2.008 and has E M 0.7998. The best RIM and HBIM solutions are marked on the graph with an asterisk. The CIM method of 19, 20 requires the HBIM and RIM solutions to match and so can be viewed as the crossing point of the two curves in Figure 2. In this case the crossing occurs when n 2 and so the CIM matches the classical method of Goodman, with an error E M 0.8. The other common value chosen for the exponent, n 3, leads to E M 2.74, 2.57 for HBIM and RIM, respectively. The solution obtained by matching at x 0 has errors E M 1.145, 1.256 for the HBIM and RIM formulations respectively, that is, 60% higher than lowest possible error. As stated earlier, matching may give good local agreement but global agreement is not ensured. The local matching solutions are shown as diamonds on the graph. As noted in 22 the value of n that minimises the error is close to 2 and so the improvements in this case by taking n 2.008 or 2.074 from the classical approach of Goodman are small. Since the gradient of E M is small near n 2, any value close to 2 will lead to a reasonable approximation. Based on this observation, it may appear pointless to use any technique other than the classical HBIM with n 2. However, in subsequent sections we will consider different boundary conditions and then find that in general n 2 is a poor choice.

Constant Flux Boundary Condition
We now consider the same problem as in 2.1a -2.1d but with u 0, t 1 replaced by the constant flux condition u x 0, t −1. This problem has the exact solution The boundary immobilisation coordinate is again ξ x/δ, but we now set u x, t δU ξ, t , chosen so that the polynomial approximation is independent of δ. Note, for the fixed temperature boundary condition we set u U whereas here we set u δU. The difference is that now as t → 0, u → 0 and this scaling permits U O 1 . In the previous case u 1 at the boundary and no rescaling was necessary. The heat equation is given by and since U t 0, this reduces to The function F ξ, t is therefore given by rearranging 4.4 Similarly the RIM gives δ 2 n 1 n 2 t/3 and E M 10n 6 − 79n 5 160n 4 38n 3 − 167n 2 − 40n 24 9n 2 2n − 1 2n 1 2n − 3 .

4.7
In both cases the error is independent of time and E M δE L . The corresponding expression using Langford's method has an error E L ∝ t −1/2 and so, as in the constant temperature problem, it blows up as t → 0. Results are shown in Figure 3 for n ∈ 3, 4.5 . The minimum value of the errors occurs at n 3.535, 3.798 for the HBIM and RIM, respectively, and then E M 0.0097, 0.0124. For this situation matching the approximate and exact solutions at u 0, t 0 works well, leading to n 3.65, see 28 , with errors only slightly above the minimum values, E M 0.01, 0.014. We do not show the error for n 2 since this is significantly higher and its inclusion makes the other results less clear for the RIM the error increases by a factor 50 from the minimum . The other standard choice, n 3, gives an HBIM error twice that of the minimum value while the RIM error is quadrupled. The two curves cross at n 4, which is the CIM prediction, with a 50% increase over the minimum error.

Newton Cooling
Langford wisely avoided this problem, which does not lead to an analytical expression for the error. The boundary condition u 0, t 1 is now replaced by the cooling condition u x u − 1. Again this problem has an exact solution International Journal of Differential Equations

11
Setting ξ x/δ and u x, t δU ξ, t leads to the polynomial profile Note, if we assume a constant n / 0, then for small times δ → 0, and so n δ, which means that U tends to the constant flux solution, 4.2 . For large t, n δ and U behaves like the fixed boundary temperature solution. The two limits show a different dependence δ n and we must conclude n is also a function of time.
In this case, after integrating 4.3 , it follows that the HBIM and RIM formulations are given by

4.10
We define F ξ, t by 4.5 and then the error E M

4.11
Although this appears rather daunting, all terms can be integrated analytically. For the CIM it is relatively simple to calculate δ t and n t from the two first-order ODEs that result from 4.10 the initial conditions are δ 0 0, n 0 4, where the n 0 value comes from the constant flux solution . In this case since n n t , we find E M is time dependent, but there is no initial singularity E M 0 0 . To avoid the difficulties of analysing a time-dependent n, in 22 the solution was calculated with n n 0 . This value was chosen since, due to the singular nature of E L , the highest error occurred at t 0. In Figure 4 we compare the errors over time for the CIM solution and the HBIM and RIM solutions that use a constant n n 0 , for t ∈ 0, 0.05 . For small times the HBIM and RIM solutions, with n 3.535, 3.798, respectively, show the smallest error, but as the optimal value of n changes, the accuracy is lost. The CIM, which has a time-dependent n, quickly becomes more accurate than the RIM approximation and around t ≈ 0.09 it becomes more accurate than the HBIM.

Application to Deicing Systems
We now consider an industrial application where an ice layer in a cold environment is heated from below. This example is motivated by deicing systems, see for example 26,29 . Initially we will assume that the ice is in a steady state, determined by the ambient conditions. Energy is then applied to the lower surface; this represents switching on a de-icing device. The important issue here is the time taken for the ice to reach the melting temperature. After this, the thermal model is less important, since ice may then be free to slide off. Consequently we now study the process until the lower surface starts to melt.
The initial temperature of the ice is governed by a steady-state heat equation, u xx 0, subject to a fixed temperature at the base and a cooling condition where the ice is exposed to the ambient air flow. This leads to a temperature profile of the form u −1 αx, 5.1 and this provides the initial condition for the subsequent state. When the de-icing system is turned on, we require a heat flux at x 0 and the problem is governed by ∂u ∂t where we have taken our length-scale from the heat flux condition, rather than the ice thickness, hence H / 1. In fact this system is approximate in the sense that there is an assumption in line with the heat balance method that the temperature at x H is unaffected by the temperature at x 0. However, with this assumption we may find an exact solution using separation of variables

5.3
International Journal of Differential Equations 13 and this can be used to test our HBIM solution. Without the assumption we may still use a heat balance approach but lose the exact solution.
The HBIM system is the same as above, but the heat equation is defined for z ∈ 0, δ and the boundary condition at x H is replaced by two conditions u δ, t −1 αδ, u x δ, t α, where the second matches the gradient to the steady-state solution. We assume a temperature profile of the form which, after satisfying the boundary conditions, becomes

5.5
The heat balance integral is and this reduces to which is independent of α. Applying δ 0 0 gives δ n n 1 t. To carry through the new error analysis, we would again set ξ x/δ, but to allow a temperature close to −1 in the vicinity of x 0 at t 0, we choose u x, t δU ξ, t − 1. This gives a polynomial approximation independent of δ U ξ, t α − α 1 − ξ 1 n 1 α 1 − ξ n .

5.8
The new relation between u and U only differs from that of Section 4.1 by a constant. It then follows that the transformed heat equation is the same as that of Section 4.1, given by 4.4 and the function F ξ, t is defined by 4.5 . The obvious conclusion is that since this problem has a constant flux boundary condition, the values of n that minimises E M are those calculated previously, namely, n 3.535 for the HBIM and n 3.799 for the RIM. For the current problem the real quantity of interest is the time until melting begins. For the heat balance approach this may be found from 5.5 by setting u 0, t 1 0. Thus, The exact solution requires solving the nonlinear equation Using parameter values from aircraft icing models 3, 30, 31 , we find H 4.28, α 0.054 which then gives t 1 0.7073 eight terms in the sum is enough to obtain a value which does not change up to 10 −16 . The HBIM and RIM predictions of t 1 , found from 5.9 , are t 1 0.7020 and t 1 0.7129, respectively, giving errors of 0.76% and 0.79%. The CIM is slightly worse than both of these methods, with t 1 0.7204 and an error of 1.89%.
Finally, to demonstrate the accuracy of the heat balance solutions in Figure 5 we show temperature profiles before melting begins at times t 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7. The position of δ at each of these times is denoted by a circle. We only show the HBIM solution with n 3.535 because the RIM solution is very similar. It is quite clear from this figure that the approximate and exact solutions are extremely close giving further proof of the accuracy of this method.

One Phase Melting
We now extend the method to deal with one-phase Stefan problems. Well-known applications for this type of model include melting or production of ice, certain foods or metals. In Section 6.2 we study a model relevant to ablation of heat shields and laser drilling.

The Classic Stefan Problem
The basic nondimensional one-phase Stefan problem with a fixed temperature boundary condition is specified by ∂u ∂t where α satisfies the transcendental equation √ πβα erf α e α 2 1. The standard HBIM approximating function is where we assume a and m are constant 12, 20 . Then the Stefan condition in 6.2 may be integrated immediately to give s 2at/β. The function f x, t , defined in 2.7 , is given here as 6.5 From 6.5 we see that f 0, t can only be zero if a 1 since we require m > 1 : this specifies a linear temperature and in general provides a poor approximation. At x s we find that f s, t → ∞ if m < 2 and m / 0, 1 . For m ≥ 2, f s, t is finite but never zero. The error defined by Langford is as usual singular at t 0, E L ∝ t −3/2 . For the one-phase Stefan problem we propose a new error measure analogous to that discussed in Section 2.1. The only difference is that the boundary fixing coordinate is ξ x/s t , rather than x/δ t . This transforms the governing equations 6.
Note, to derive this equation, we have used the fact that from 6.8b ss t a/β and so it is clear that E M is independent of time whereas E L is once again singular at t 0 . However, it is clear that the error is dependent on the value of β and in fact we find the optimal m varies with β.
The HBIM and RIM formulations are given by The CIM is then determined from combining these equations with the Stefan condition 6.8b , which reduces to βss t a. Substituting for U it is easy to show that a and m will satisfy the same expressions as for the nontransformed system, that is, m solves the nonlinear equation For Stefan problems, with a known exact solution, there is an obvious error measure other than the accuracy of the temperature profile, namely, the value of the constant α in the expression s 2α √ t. It is interesting to note that when we consider the error in α, for β 1, the HBIM is less accurate than the RIM and CIM, whilst in Figure 6 the HBIM appears most accurate. Clearly satisfying the heat equation globally is no guarantee of accuracy in predicting the position of the melt front. For smaller values of β e.g., β < 0.1 , the HBIM is the most accurate according to both criteria.

One Phase, Subcooled Stefan Problem
We now discuss a Stefan problem which involves a one-phase semi-infinite, subcooled material. An application occurs when considering whether ice melts or water freezes when hot water is thrown over cold ice 32, 33 . This can be formulated as a Stefan problem involving a constant heat source term in the condition at the moving boundary. A similar industrial process is that of ablation, where a mass is removed from an object by vapourisation or other similar erosive processes 2, 34 .
The dimensionless model is defined as ∂u ∂t If the constraintṡ t ≥ 0 holds then for early times instead of the conditions in 6.17a and 6.17b we have ∂u ∂x −1, at x 0. 6.19 International Journal of Differential Equations The process then involves two stages, and Mitchell 34 has successfully applied the CIM to this problem. In this paper we consider the single stage process described by 6.15 -6.18 . The introduction of a subcooled region requires an analysis of the temperature there, and consequently we introduce the heat penetration depth δ. For the system above, δ is defined by replacing 6.18 with u δ t , t −1, 6.20a δ t , t 0, 6.20b To apply a heat balance integral method, we would use an approximating polynomial of the form which then automatically satisfies the boundary conditions 6.20a and 6.17a . However, to calculate E M we first need to define the boundary immobilisation coordinate which transforms the moving domain s t < x < δ t into a fixed one 0 < ξ < 1. If we denote u x, t U ξ, t then the system 6.15 -6.18 becomes The polynomial approximation, 6.21 , is now Although the profile is independent of s and δ we cannot conclude that U t 0 in 6.23 . As discussed in 33 , the conditions 6.16 and 6.17a imply that the gradient u x is infinite as t → 0 . Thus at small times the gradient term in 6.17b must be balanced with the left hand side, and this implies that s, δ ∼ √ t. However, as time increases, the constant term in 6.17b becomes important indicating a shift in solution form and so the possibility that m is time dependent.
International Journal of Differential Equations

19
The error measure here is defined as The HBIM formulation can be found from integrating 6.23 once with respect to ξ in the usual way. This gives For the RIM, since U ξ 1, t 0, we choose to first integrate over ξ, 1 . Then, on integrating again between 0, 1 and changing the order of integration, we find

6.29
These are coupled with the Stefan condition in 64 and this allows us to determine s, δ, and m. We know that s 0 δ 0 0 and the initial condition for m is found by considering the limit as t → 0 . Setting s −2μ √ t, δ 2λ √ t, with μ, λ > 0, and substituting into 64 , 6.29 leads to three algebraic equations to solve for μ, λ, m 0 m 0 . This leads to chosen to minimise E M at t 0 and the CIM. The CIM is far more accurate here because m is allowed to vary with time. In the right plot we compare the constant value of m used in the HBIM and RIM which are superimposed on this plot with the time-dependent m from the CIM. We see that m varies significantly with time and thus using the time dependent value is much more accurate. As β increases, the variation in m with time is less pronounced in the CIM and then the HBIM and RIM solutions improve.

Conclusions
In this paper we have shown how the heat balance method, in various guises, may be applied to thermal and Stefan problems. We illustrated the method on a range of industrially important problems. The application of the method to further problems, we believe, is obvious.
The HBIM has been criticised for a lack of accuracy, hence we developed an error measure, which does not require knowledge of an exact solution and verified the accuracy of our solutions. We noted that if the exponent of the approximating polynomial is not specified at the outset, then methods for determining the exponent fall into two categories, namely, local and global matching methods. By defining the function f u t − u xx , it was also shown that the standard polynomial approximation will never satisfy the heat equation at x 0 and it is only satisfied at x δ when n > 2. This indicates that a more general approximation should be used. However, our goal in this paper was to analyse a number of physically important problems rather than develop new forms of HBIM, and hence we have not investigated alternative forms for the temperature.