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

Authors: S. L. Mitchell; T. G. Myers
Journal: International Journal of Differential Equations
 (2012)

---

## Abstract

We employ the Heat Balance Integral Method (HBIM) 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.


---

## Body

## 1. 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. Goodman [6–8] 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, EL∝t-3/2. Despite the fact that the HBIM approximation may appear reasonable at small times, this measure shows that EL→∞ 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 a criteria which does not blow up ast→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 semi-infinite 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.


## 2. Thermal Problem with Fixed Boundary Temperature

We begin by illustrating the HBIM on a simple nondimensional thermal problem described by(2.1a)∂u∂t=∂2u∂x2,(2.1b)u(0,t)=1,(2.1c)u(x,0)=0,(2.1d)u⟶0asx⟶∞,which has the exact solution(2.2)u=erfc⁡(x2t).
To calculate the HBIM solution, we replace the boundary condition at infinity with(2.3)u(δ(t),t)=∂u∂x(δ(t),t)=0,
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(2.4)u(x,t)=(1-xδ)n.
The HBIM then proceeds by integrating the heat equation over x∈[0,δ],(2.5)∫0δ∂u∂tdx=∫0δ∂2u∂x2dx.
Substituting for u(x,t), via (2.4), and integrating lead to(2.6)δ=2n(n+1)t.
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 functionf(x,t) defined by(2.7)f(x,t)=∂u∂t-∂2u∂x2.
The HBIM formulation and the least-squares error of Langford may then be written as(2.8)∫0δf(x,t)dx=0,EL=∫0δf(x,t)2dx.
If u(x,t) is an exact solution of the heat equation then both f and EL 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 EL>0. This may be seen more clearly by considering the form of these functions. Noting the derivatives from (2.4)(2.9)∂u∂t=nxδtδ2(1-xδ)n-1,∂2u∂x2=n(n-1)δ2(1-xδ)n-2,
we find(2.10)f(x,t)=nδ2(1-xδ)n-2[xδt(1-xδ)-(n-1)],
and(2.11)EL=(2n4-7n3+6n2+2n-1)2n(n+1)4(n+1)2(2n-3)(2n-1)(2n+1)t-3/2=eLt-3/2.
The above expressions highlight some of the problems inherent in applying heat balance methods. For example, at x=0,(2.12)f(0,t)=-n(n-1)δ2.
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 atx=0. At  x=δ, we find(2.13)f(δ,t)=lim⁡x⟶δf(x,t)=-n(n-1)δ2lim⁡x⟶δ(1-xδ)n-2.
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 forEL and eL were first derived by Langford [21]. In Figure 1 we plot eL(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 EL=0. Of course it should be expected that an approximate solution has some error, so perhaps the most worrying aspect of the expression for EL is that it is time-dependent and more specifically, as t→0, the error EL→∞ even though the approximate solution may appear accurate. Hence, while EL may be useful to represent relative errors, that is, for a given time, we can generate EL 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.Figure 1
Plot ofeL(n) for n∈(1.5,5].The solid line denotes the HBIM formulation and the dashed line denotes the RIM formulation.


### 2.1. A New Error Measure

To define a new error measure that avoids the unrealistic behaviour ofEL, we use a boundary immobilisation technique. For thermal problems this simply means that we introduce a new coordinate(2.14)ξ=xδ(t),
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(2.15)∂U∂t-δtξδ∂U∂ξ=1δ2∂2U∂ξ2,
and the boundary conditions are(2.16)U(0,t)=1,U(1,t)=∂U∂ξ(1,t)=0.
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(2.17)F(ξ,t)=δ2∂U∂t-δδtξ∂U∂ξ-∂2U∂ξ2(2.18)=-δδtξ∂U∂ξ-∂2U∂ξ2,
after noting Ut=0. The HBIM formulation may then be written ∫01F(ξ,t)dξ, which leads to(2.19)δ=2n(n+1)t,
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:(2.20)EM=∫01F(ξ,t)2dξ=∫01[∂2U∂ξ2+δδtξ∂U∂ξ]2dξ.
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(2.21)EM=n2(2n4-7n3+6n2+2n-1)(2n-3)(2n-1)(2n+1).
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 EM=δ3EL.The RIM simply involves integrating the heat equation twice(2.22)∫01∫0ηF(ξ,t)dξdη=0,
where η is a dummy variable. After changing the order of integration this becomes ∫01(1-ξ)F(ξ,t)dξ=0 and then the HBIM formulation ∫01F(ξ,t)dξ=0 is substituted to give ∫01ξF(ξ,t)dξ=0. This leads to(2.23)δ=(n+1)(n+2)t.
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(2.24)EM=n(n-1)(10n4-29n3+5n2+32n+12)4(2n-3)(2n-1)(2n+1)
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 byEM.


## 2.1. A New Error Measure

To define a new error measure that avoids the unrealistic behaviour ofEL, we use a boundary immobilisation technique. For thermal problems this simply means that we introduce a new coordinate(2.14)ξ=xδ(t),
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(2.15)∂U∂t-δtξδ∂U∂ξ=1δ2∂2U∂ξ2,
and the boundary conditions are(2.16)U(0,t)=1,U(1,t)=∂U∂ξ(1,t)=0.
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(2.17)F(ξ,t)=δ2∂U∂t-δδtξ∂U∂ξ-∂2U∂ξ2(2.18)=-δδtξ∂U∂ξ-∂2U∂ξ2,
after noting Ut=0. The HBIM formulation may then be written ∫01F(ξ,t)dξ, which leads to(2.19)δ=2n(n+1)t,
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:(2.20)EM=∫01F(ξ,t)2dξ=∫01[∂2U∂ξ2+δδtξ∂U∂ξ]2dξ.
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(2.21)EM=n2(2n4-7n3+6n2+2n-1)(2n-3)(2n-1)(2n+1).
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 EM=δ3EL.The RIM simply involves integrating the heat equation twice(2.22)∫01∫0ηF(ξ,t)dξdη=0,
where η is a dummy variable. After changing the order of integration this becomes ∫01(1-ξ)F(ξ,t)dξ=0 and then the HBIM formulation ∫01F(ξ,t)dξ=0 is substituted to give ∫01ξF(ξ,t)dξ=0. This leads to(2.23)δ=(n+1)(n+2)t.
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(2.24)EM=n(n-1)(10n4-29n3+5n2+32n+12)4(2n-3)(2n-1)(2n+1)
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 byEM.


## 3. Defining the Exponentn

The classical HBIM involves specifyingn 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.


### 3.1. Local Matching

For this methodn 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 ux(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(3.1)ux(0,t)=-1πt=-nδ.
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 ux2/u2 at x=0. For a fixed boundary temperature, this is equivalent to matching ux 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 the temperature or gradient at a single point isno guarantee that the approximate solution will provide a good approximation throughout the domain. It is not even certain that it will provide a good solution to the heat equation even at the point of matching. For example, the discussion of the value of f(x,t), given after (2.12), (2.13), shows that n≈1.752 does not satisfy the heat equation at x=0 and leads to an infinite value of f at x=δ.


### 3.2. Global Matching

The method described in [22] involves leaving n unknown throughout the calculation. It is then determined by minimising EL 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 eL(n) in (2.11). Clearly this method may be easily modified to the present situation by choosing n to minimise EM(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 EL). 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.


### 3.3. Comparison of Results

In Figure2 we present the variation of EM 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 EM=0.784. The most accurate HBIM approximation requires n=2.008 and has EM=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 EM=0.8. The other common value chosen for the exponent, n=3, leads to EM=2.74,2.57 for HBIM and RIM, respectively. The solution obtained by matching at x=0 has errors EM=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 EM 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.Figure 2
EM(n) for the fixed temperature boundary condition, for HBIM (solid line) and RIM (dashed line): errors for n=1.752,2,3 are shown as diamonds and the minimum errors are shown as asterisks.


## 3.1. Local Matching

For this methodn 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 ux(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(3.1)ux(0,t)=-1πt=-nδ.
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 ux2/u2 at x=0. For a fixed boundary temperature, this is equivalent to matching ux 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 the temperature or gradient at a single point isno guarantee that the approximate solution will provide a good approximation throughout the domain. It is not even certain that it will provide a good solution to the heat equation even at the point of matching. For example, the discussion of the value of f(x,t), given after (2.12), (2.13), shows that n≈1.752 does not satisfy the heat equation at x=0 and leads to an infinite value of f at x=δ.


## 3.2. Global Matching

The method described in [22] involves leaving n unknown throughout the calculation. It is then determined by minimising EL 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 eL(n) in (2.11). Clearly this method may be easily modified to the present situation by choosing n to minimise EM(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 EL). 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.


## 3.3. Comparison of Results

In Figure2 we present the variation of EM 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 EM=0.784. The most accurate HBIM approximation requires n=2.008 and has EM=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 EM=0.8. The other common value chosen for the exponent, n=3, leads to EM=2.74,2.57 for HBIM and RIM, respectively. The solution obtained by matching at x=0 has errors EM=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 EM 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.Figure 2
EM(n) for the fixed temperature boundary condition, for HBIM (solid line) and RIM (dashed line): errors for n=1.752,2,3 are shown as diamonds and the minimum errors are shown as asterisks.


## 4. Constant Flux and Newton Cooling Boundary Conditions

### 4.1. 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 ux(0,t)=-1. This problem has the exact solution(4.1)u(x,t)=2tπe-x2/(4t)-xerfc⁡x2t.
The boundary immobilisation coordinate is again ξ=x/δ, but we now set u(x,t)=δU(ξ,t), chosen so that the polynomial approximation(4.2)U(ξ,t)=1n(1-ξ)n
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=𝒪(1). In the previous case u=1 at the boundary and no rescaling was necessary. The heat equation is given by(4.3)∂∂t(δU)-δtξ∂U∂ξ=1δ∂2U∂ξ2,
and since Ut=0, this reduces to(4.4)δδtU-δδtξ∂U∂ξ=∂2U∂ξ2.
The function F(ξ,t) is therefore given by rearranging  (4.4)(4.5)F(ξ,t)=δδtU-δδtξ∂U∂ξ-∂2U∂ξ2.
Again the HBIM is the integral of F over ξ∈[0,1], the RIM involves the double integral and the error EM is the integral of F2.The HBIM formulation givesδ=n(n+1)t and(4.6)EM=(n-1)(2n3-13n2+19n+10)4(2n-1)(2n+1)(2n-3).
Similarly the RIM gives δ=2(n+1)(n+2)t/3 and(4.7)EM=10n6-79n5+160n4+38n3-167n2-40n+249n2(2n-1)(2n+1)(2n-3).
In both cases the error is independent of time and EM=δEL. The corresponding expression using Langford’s method has an error EL∝t-1/2 and so, as in the constant temperature problem, it blows up as t→0.Results are shown in Figure3 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 EM=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, EM=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.Figure 3
EM(n) for the constant flux problem, for HBIM (solid line) and RIM (dashed line): errors for n=3,3.65,4 are shown as diamonds and the minimum errors are shown as asterisks.


### 4.2. Newton Cooling

Langford wisely avoided this problem, which does not lead to an analytical expression for the error. The boundary conditionu(0,t)=1 is now replaced by the cooling condition ux=u-1. Again this problem has an exact solution(4.8)u(x,t)=erfc⁡x2t-ex+terfc⁡x+2t2t.
Setting ξ=x/δ and u(x,t)=δU(ξ,t) leads to the polynomial profile(4.9)U=1n+δ(1-ξ)n.
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)δδt∫01Udξ+δddt∫01δUdξ=-∂U∂ξ|ξ=02δδt∫01ξUdξ+δddt∫01δξUdξ=U|ξ=0.
We define F(ξ,t) by (4.5) and then the error EM (=∫01F(ξ,t)2dξ) is(4.11)EM=∫01[nδ2δt(n+δ)2(1-ξ)n+δ3ntn+δ(1-ξ)nln⁡(1-ξ)-δ3nt(n+δ)2(1-ξ)n+nδ2δtn+δξ(1-ξ)n-1-n(n-1)δn+δ(1-ξ)n-2]2dξ.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 EM is time dependent, but there is no initial singularity (EM(0)=0).To avoid the difficulties of analysing a time-dependentn, in [22] the solution was calculated with n=n(0). This value was chosen since, due to the singular nature of EL, 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.Figure 4
EM against t for the cooling condition at small times.


## 4.1. 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 ux(0,t)=-1. This problem has the exact solution(4.1)u(x,t)=2tπe-x2/(4t)-xerfc⁡x2t.
The boundary immobilisation coordinate is again ξ=x/δ, but we now set u(x,t)=δU(ξ,t), chosen so that the polynomial approximation(4.2)U(ξ,t)=1n(1-ξ)n
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=𝒪(1). In the previous case u=1 at the boundary and no rescaling was necessary. The heat equation is given by(4.3)∂∂t(δU)-δtξ∂U∂ξ=1δ∂2U∂ξ2,
and since Ut=0, this reduces to(4.4)δδtU-δδtξ∂U∂ξ=∂2U∂ξ2.
The function F(ξ,t) is therefore given by rearranging  (4.4)(4.5)F(ξ,t)=δδtU-δδtξ∂U∂ξ-∂2U∂ξ2.
Again the HBIM is the integral of F over ξ∈[0,1], the RIM involves the double integral and the error EM is the integral of F2.The HBIM formulation givesδ=n(n+1)t and(4.6)EM=(n-1)(2n3-13n2+19n+10)4(2n-1)(2n+1)(2n-3).
Similarly the RIM gives δ=2(n+1)(n+2)t/3 and(4.7)EM=10n6-79n5+160n4+38n3-167n2-40n+249n2(2n-1)(2n+1)(2n-3).
In both cases the error is independent of time and EM=δEL. The corresponding expression using Langford’s method has an error EL∝t-1/2 and so, as in the constant temperature problem, it blows up as t→0.Results are shown in Figure3 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 EM=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, EM=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.Figure 3
EM(n) for the constant flux problem, for HBIM (solid line) and RIM (dashed line): errors for n=3,3.65,4 are shown as diamonds and the minimum errors are shown as asterisks.


## 4.2. Newton Cooling

Langford wisely avoided this problem, which does not lead to an analytical expression for the error. The boundary conditionu(0,t)=1 is now replaced by the cooling condition ux=u-1. Again this problem has an exact solution(4.8)u(x,t)=erfc⁡x2t-ex+terfc⁡x+2t2t.
Setting ξ=x/δ and u(x,t)=δU(ξ,t) leads to the polynomial profile(4.9)U=1n+δ(1-ξ)n.
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)δδt∫01Udξ+δddt∫01δUdξ=-∂U∂ξ|ξ=02δδt∫01ξUdξ+δddt∫01δξUdξ=U|ξ=0.
We define F(ξ,t) by (4.5) and then the error EM (=∫01F(ξ,t)2dξ) is(4.11)EM=∫01[nδ2δt(n+δ)2(1-ξ)n+δ3ntn+δ(1-ξ)nln⁡(1-ξ)-δ3nt(n+δ)2(1-ξ)n+nδ2δtn+δξ(1-ξ)n-1-n(n-1)δn+δ(1-ξ)n-2]2dξ.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 EM is time dependent, but there is no initial singularity (EM(0)=0).To avoid the difficulties of analysing a time-dependentn, in [22] the solution was calculated with n=n(0). This value was chosen since, due to the singular nature of EL, 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.Figure 4
EM against t for the cooling condition at small times.


## 5. 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,uxx=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(5.1)u=-1+αx,
and this provides the initial condition for the subsequent state.When the de-icing system is turned on, we require a heat flux atx=0 and the problem is governed by(5.2)∂u∂t=∂2u∂x2,0<x<H,∂u∂x|x=0=-1,u(H,t)=-1+αH,u(x,0)=-1+αx,
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)u(x,t)=-1+(1+α)H-x-∑n=0∞8H(1+α)π2(2n+1)2cos⁡[(2n+1)πx2H]×exp⁡(-(2n+1)2π2t4H2),
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 forz∈[0,δ] and the boundary condition at x=H  is replaced by two conditions u(δ,t)=-1+αδ, ux(δ,t)=α, where the second matches the gradient to the steady-state solution. We assume a temperature profile of the form(5.4)u(x,t)=a+b(1-xδ)+c(1-xδ)n,
which, after satisfying the boundary conditions, becomes(5.5)u=-1+αδ-αδ(1-xδ)+1n(1+α)δ(1-xδ)n.
The heat balance integral is(5.6)ddt∫0δudx-u|x=δδt=ux|x=δ-ux|x=0,
and this reduces to(5.7)δδt=n(n+1)2,
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 δ(5.8)U(ξ,t)=α-α(1-ξ)+1n(1+α)(1-ξ)n.
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 EM 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,t1)=0. Thus,(5.9)t1=n(n+1)(1+α)2,δ(t1)=n1+α.
The exact solution requires solving the nonlinear equation(5.10)1=(1+α)H-∑n=0∞8H(1+α)π2(2n+1)2exp⁡(-(2n+1)2π2t14H2).
Using parameter values from aircraft icing models [3, 30, 31], we find H=4.28, α=0.054 which then gives t1=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 t1, found from (5.9), are t1=0.7020 and t1=0.7129, respectively, giving errors of 0.76% and 0.79%. The CIM is slightly worse than both of these methods, with t1=0.7204 and an error of 1.89%.Finally, to demonstrate the accuracy of the heat balance solutions in Figure5 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.Figure 5
Solution profiles of the aircraft icing example before melting begins, at nondimensional timest=0 (dotted line) and t=0.1,0.2,0.3,0.4,0.5,0.6,0.7. The solid line denotes the exact solution and the dashed line denotes the HBIM solution (with n=3.535).


## 6. 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 Section6.2 we study a model relevant to ablation of heat shields and laser drilling.


### 6.1. The Classic Stefan Problem

The basic nondimensional one-phase Stefan problem with a fixed temperature boundary condition is specified by(6.1)∂u∂t=∂2u∂x2,0<x<s(t),(6.2)u(0,t)=1,u(s,t)=0,βdsdt=∂u∂x|x=s,
and this has exact solution(6.3)u(x,t)=1-erf⁡[x/(2t)]erf⁡(α),s(t)=2αt,
where α satisfies the transcendental equation πβαerf(α)eα2=1. The standard HBIM approximating function is(6.4)u=a(1-xs)+(1-a)(1-xs)m,
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)f(x,t)=xsts2[a+m(1-a)(1-xs)m-1]-m(m-1)(1-a)s2(1-xs)m-2.
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, EL∝t-3/2.For the one-phase Stefan problem we propose a new error measure analogous to that discussed in Section2.1. The only difference is that the boundary fixing coordinate is ξ=x/s(t), rather than x/δ(t). This transforms the governing equations (6.1), (6.2) into(6.6)∂U∂t-stsξ∂U∂ξ=1s2∂2U∂ξ2,0<ξ<1,(6.7)U(0,t)=1,
(6.8a)U(1,t)=0,(6.8b)βsst=-∂U∂ξ|ξ=1, and  (6.6) is equivalent to (2.15) with δ replaced by s. The approximating profile is
(6.9)U=a(1-ξ)+(1-a)(1-ξ)m.
Since a is constant, it follows that Ut=0 and so (equivalent to (2.18))
(6.10)F=-sstξ∂U∂ξ-∂2U∂ξ2.
The error is defined as EM=∫01F2dξ or
(6.11)EM=∫01[m(m-1)(1-a)(1-ξ)m-2-a2βξ-m(1-a)aβξ(1-ξ)m-1]2dξ.
Note, to derive this equation, we have used the fact that (from (6.8b)) sst=a/β and so it is clear that EM is independent of time (whereas EL 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(6.12)sst∫01Udξ=∂U∂ξ|ξ=1-∂U∂ξ|ξ=0,2sst∫01ξUdξ=1+∂U∂ξ|ξ=1.
The CIM is then determined from combining these equations with the Stefan condition (6.8b), which reduces to βsst=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(6.13)β(2m-3)(2+m)(1+m)(2m2+5m-6)+12(m-2)=0,
and a satisfies(6.14)a=6(2-m)(m-1)(2m2+5m-6),
further details are given in [20].In Figure6 we present the variation of EM for the HBIM and RIM formulations for m∈[1.5,2.1] and two values of β. With β=1 we see that the most accurate approximation comes from using the RIM formulation with m=1.793, which has an error of 0.021. The most accurate HBIM approximation requires m=1.786 and has a slightly smaller error of EM=0.020. The CIM predicts m=1.546, with a much higher error of EM=0.113. Similarly, when β=10, we find that m=1.761,1.765,1.506 and EM=3.8×10-4,4.3×10-4,3.0×10-3 for the HBIM, RIM, and CIM formulations, respectively. Also shown are solutions for the HBIM and RIM with m=2 and the solution found by matching ux at x=0 to the exact value (found from (6.3)).EM(m) for HBIM (solid line) and RIM (dashed line), minimum errors shown as asterisks. (a) β=1, diamonds show errors for m=1.5531, m=1.5514 (obtained by matching ux at x=0), m=1.5463 (CIM) and m=2.(b) β=10, diamonds show errors for m=1.5061, m=1.5060 (matching ux at x=0), m=1.5056 (CIM) and m=2.
(a)(b)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.


### 6.2. 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(6.15)∂u∂t=∂2u∂x2,s(t)<x<∞,(6.16)u(x,0)=-1,s(0)=0,
(6.17a)u(s,t)=0,(6.17b)βdsdt=1+∂u∂x|x=s,   (6.18)u⟶-1asx⟶∞.
If the constraint ṡ(t)≥0 holds then for early times instead of the conditions in (6.17a) and (6.17b) we have(6.19)∂u∂x=-1,atx=0.
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(6.20a)u(δ(t),t)=-1,(6.20b)(δ(t),t)=0,(6.20c)δ(0)=0.To apply a heat balance integral method, we would use an approximating polynomial of the form(6.21)u=-1+(δ-xδ-s)m,s<x<δ,
which then automatically satisfies the boundary conditions (6.20a) and (6.17a). However, to calculate EM we first need to define the boundary immobilisation coordinate(6.22)ξ=x-sδ-s,
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(6.23)∂2U∂ξ2=(δ-s)2∂U∂t-(δ-s)[st+(δt-st)ξ]∂U∂ξ,
with(6.24)U(0,t)=0,U(1,t)=-1,∂U∂ξ(1,t)=0,βst=1+1δ-s∂U∂ξ|ξ=0.
The polynomial approximation, (6.21), is now(6.25)U=-1+(1-ξ)m.
Although the profile is independent of s and δ we cannot conclude that Ut=0 in (6.23). As discussed in [33], the conditions (6.16) and (6.17a) imply that the gradient ux 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.The error measure here is defined as(6.26)EM=∫01F(ξ,t)2dξ,F(ξ,t)=∂2U∂ξ2-(δ-s)2∂U∂t+(δ-s)[st+(δt-st)ξ]∂U∂ξ.
The HBIM formulation can be found from integrating (6.23) once with respect to ξ in the usual way. This gives(6.27)(δ-s)(δt-st)∫01Udξ+(δ-s)2ddt∫01Udξ+(δ-s)δt=-∂U∂ξ|ξ=0.
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.28)2(δ-s)(δt-st)∫01ξUdξ+(δ-s)2ddt∫01ξUdξ+(δ-s)st∫01Udξ+(δ-s)δt=1.
After substituting the profile (6.25) into the integral equations (6.27) and (6.28) we have the pair of ODEs(6.29)δt+mst-δ-sm+1mt=m(m+1)δ-s,2δt+mst-(2m+3)(δ-s)(m+1)(m+2)mt=(m+1)(m+2)δ-s.
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)=m0. This leads to(6.30)βμ=m02(λ+μ),m0(m0+1)2(λ+μ)=-m0μ+λ,(m0+1)(m0+2)2(λ+μ)=-m0μ+2λ.
We can eliminate μ and λ to give an explicit expression for m0 in terms of β:(6.31)m0=β+9β2+8β2(1+β).
Note, only the positive root of the quadratic equation ensures m0>0. In fact, to ensure m0>1 requires β>0.5: for β<0.5, the CIM breaks down (note, most applications of this model involve values of β>1 [32–34]).In Figure7 we have plotted s and m against t for β=1. Comparisons are shown between a numerical solution, using the Keller box scheme [35], the HBIM and RIM (with m chosen to minimise EM 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.Plots ofs against t (a) and m against  t (b) when β=1.
(a)(b)


## 6.1. The Classic Stefan Problem

The basic nondimensional one-phase Stefan problem with a fixed temperature boundary condition is specified by(6.1)∂u∂t=∂2u∂x2,0<x<s(t),(6.2)u(0,t)=1,u(s,t)=0,βdsdt=∂u∂x|x=s,
and this has exact solution(6.3)u(x,t)=1-erf⁡[x/(2t)]erf⁡(α),s(t)=2αt,
where α satisfies the transcendental equation πβαerf(α)eα2=1. The standard HBIM approximating function is(6.4)u=a(1-xs)+(1-a)(1-xs)m,
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)f(x,t)=xsts2[a+m(1-a)(1-xs)m-1]-m(m-1)(1-a)s2(1-xs)m-2.
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, EL∝t-3/2.For the one-phase Stefan problem we propose a new error measure analogous to that discussed in Section2.1. The only difference is that the boundary fixing coordinate is ξ=x/s(t), rather than x/δ(t). This transforms the governing equations (6.1), (6.2) into(6.6)∂U∂t-stsξ∂U∂ξ=1s2∂2U∂ξ2,0<ξ<1,(6.7)U(0,t)=1,
(6.8a)U(1,t)=0,(6.8b)βsst=-∂U∂ξ|ξ=1, and  (6.6) is equivalent to (2.15) with δ replaced by s. The approximating profile is
(6.9)U=a(1-ξ)+(1-a)(1-ξ)m.
Since a is constant, it follows that Ut=0 and so (equivalent to (2.18))
(6.10)F=-sstξ∂U∂ξ-∂2U∂ξ2.
The error is defined as EM=∫01F2dξ or
(6.11)EM=∫01[m(m-1)(1-a)(1-ξ)m-2-a2βξ-m(1-a)aβξ(1-ξ)m-1]2dξ.
Note, to derive this equation, we have used the fact that (from (6.8b)) sst=a/β and so it is clear that EM is independent of time (whereas EL 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(6.12)sst∫01Udξ=∂U∂ξ|ξ=1-∂U∂ξ|ξ=0,2sst∫01ξUdξ=1+∂U∂ξ|ξ=1.
The CIM is then determined from combining these equations with the Stefan condition (6.8b), which reduces to βsst=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(6.13)β(2m-3)(2+m)(1+m)(2m2+5m-6)+12(m-2)=0,
and a satisfies(6.14)a=6(2-m)(m-1)(2m2+5m-6),
further details are given in [20].In Figure6 we present the variation of EM for the HBIM and RIM formulations for m∈[1.5,2.1] and two values of β. With β=1 we see that the most accurate approximation comes from using the RIM formulation with m=1.793, which has an error of 0.021. The most accurate HBIM approximation requires m=1.786 and has a slightly smaller error of EM=0.020. The CIM predicts m=1.546, with a much higher error of EM=0.113. Similarly, when β=10, we find that m=1.761,1.765,1.506 and EM=3.8×10-4,4.3×10-4,3.0×10-3 for the HBIM, RIM, and CIM formulations, respectively. Also shown are solutions for the HBIM and RIM with m=2 and the solution found by matching ux at x=0 to the exact value (found from (6.3)).EM(m) for HBIM (solid line) and RIM (dashed line), minimum errors shown as asterisks. (a) β=1, diamonds show errors for m=1.5531, m=1.5514 (obtained by matching ux at x=0), m=1.5463 (CIM) and m=2.(b) β=10, diamonds show errors for m=1.5061, m=1.5060 (matching ux at x=0), m=1.5056 (CIM) and m=2.
(a)(b)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.


## 6.2. 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(6.15)∂u∂t=∂2u∂x2,s(t)<x<∞,(6.16)u(x,0)=-1,s(0)=0,
(6.17a)u(s,t)=0,(6.17b)βdsdt=1+∂u∂x|x=s,   (6.18)u⟶-1asx⟶∞.
If the constraint ṡ(t)≥0 holds then for early times instead of the conditions in (6.17a) and (6.17b) we have(6.19)∂u∂x=-1,atx=0.
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(6.20a)u(δ(t),t)=-1,(6.20b)(δ(t),t)=0,(6.20c)δ(0)=0.To apply a heat balance integral method, we would use an approximating polynomial of the form(6.21)u=-1+(δ-xδ-s)m,s<x<δ,
which then automatically satisfies the boundary conditions (6.20a) and (6.17a). However, to calculate EM we first need to define the boundary immobilisation coordinate(6.22)ξ=x-sδ-s,
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(6.23)∂2U∂ξ2=(δ-s)2∂U∂t-(δ-s)[st+(δt-st)ξ]∂U∂ξ,
with(6.24)U(0,t)=0,U(1,t)=-1,∂U∂ξ(1,t)=0,βst=1+1δ-s∂U∂ξ|ξ=0.
The polynomial approximation, (6.21), is now(6.25)U=-1+(1-ξ)m.
Although the profile is independent of s and δ we cannot conclude that Ut=0 in (6.23). As discussed in [33], the conditions (6.16) and (6.17a) imply that the gradient ux 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.The error measure here is defined as(6.26)EM=∫01F(ξ,t)2dξ,F(ξ,t)=∂2U∂ξ2-(δ-s)2∂U∂t+(δ-s)[st+(δt-st)ξ]∂U∂ξ.
The HBIM formulation can be found from integrating (6.23) once with respect to ξ in the usual way. This gives(6.27)(δ-s)(δt-st)∫01Udξ+(δ-s)2ddt∫01Udξ+(δ-s)δt=-∂U∂ξ|ξ=0.
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.28)2(δ-s)(δt-st)∫01ξUdξ+(δ-s)2ddt∫01ξUdξ+(δ-s)st∫01Udξ+(δ-s)δt=1.
After substituting the profile (6.25) into the integral equations (6.27) and (6.28) we have the pair of ODEs(6.29)δt+mst-δ-sm+1mt=m(m+1)δ-s,2δt+mst-(2m+3)(δ-s)(m+1)(m+2)mt=(m+1)(m+2)δ-s.
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)=m0. This leads to(6.30)βμ=m02(λ+μ),m0(m0+1)2(λ+μ)=-m0μ+λ,(m0+1)(m0+2)2(λ+μ)=-m0μ+2λ.
We can eliminate μ and λ to give an explicit expression for m0 in terms of β:(6.31)m0=β+9β2+8β2(1+β).
Note, only the positive root of the quadratic equation ensures m0>0. In fact, to ensure m0>1 requires β>0.5: for β<0.5, the CIM breaks down (note, most applications of this model involve values of β>1 [32–34]).In Figure7 we have plotted s and m against t for β=1. Comparisons are shown between a numerical solution, using the Keller box scheme [35], the HBIM and RIM (with m chosen to minimise EM 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.Plots ofs against t (a) and m against  t (b) when β=1.
(a)(b)


## 7. 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 functionf=ut-uxx, 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.



---

*Source: 187902-2012-04-23.xml*