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.

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 [15]. 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 [68] 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 [1317].

The standard heat balance integral method approximates solutions to the heat equation by first introducing a heat penetration depth, , where for 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 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, . Despite the fact that the HBIM approximation may appear reasonable at small times, this measure shows that as .

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 as . 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 . 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 , 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, 2427]. In both cases the new error is shown to be constant, whereas Langford’s definition is again singular at .

2. Thermal Problem with Fixed Boundary Temperature

We begin by illustrating the HBIM on a simple nondimensional thermal problem described by

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 , where . The appropriate approximating polynomial for is The HBIM then proceeds by integrating the heat equation over , Substituting for , via (2.4), and integrating lead to With traditional heat balance methods, the value of is specified at the outset, typically or 3, and so the solution is given by (2.4) and (2.6) and the analysis is complete.

Now, consider the function defined by The HBIM formulation and the least-squares error of Langford may then be written as If is an exact solution of the heat equation then both and will be identically zero. However, with the HBIM the heat equation is only satisfied in an integral sense and for all . That is, the heat equation is not satisfied by the polynomial approximation and the error . This may be seen more clearly by considering the form of these functions. Noting the derivatives from (2.4) we find and The above expressions highlight some of the problems inherent in applying heat balance methods. For example, at , For to satisfy the second boundary condition in (2.3) we require and so (2.12) indicates solutions of the form (2.4) will never satisfy the heat equation at . At , we find With the restriction , this shows that the heat equation is only satisfied at when . For (and ), and (2.13) shows , with we find .

The expressions for and were first derived by Langford [21]. In Figure 1 we plot as specified by (2.11) and also the equivalent curve for the RIM formulation. The expression in (2.11) is only valid for , 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 , the expression (2.11) is positive and there is no value of that can make . Of course it should be expected that an approximate solution has some error, so perhaps the most worrying aspect of the expression for is that it is time-dependent and more specifically, as , the error even though the approximate solution may appear accurate. Hence, while may be useful to represent relative errors, that is, for a given time, we can generate for different values of 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.

2.1. A New Error Measure

To define a new error measure that avoids the unrealistic behaviour of , we use a boundary immobilisation technique. For thermal problems this simply means that we introduce a new coordinate which transforms the moving domain to a fixed one . If we denote the heat equation (2.1a) becomes and the boundary conditions are The polynomial approximation, (2.4), is now . 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 corresponds to the similarity solution transformation .

Equivalent to deriving equation (2.7), we rearrange the heat equation, (2.15), to define the function after noting . The HBIM formulation may then be written , which leads to 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 is independent of time this function clearly only varies with the value of . 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 . It is related to the expression of Langford by .

The RIM simply involves integrating the heat equation twice where is a dummy variable. After changing the order of integration this becomes and then the HBIM formulation is substituted to give . This leads to 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 .

3. Defining the Exponent

The classical HBIM involves specifying at the start of a calculation. However, there exist a number of refinements which lead to different values of . 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 method is determined through an additional boundary condition, usually taken from the exact solution. For example, in the current problem the boundary condition specifies and so, from the exact solution, one can also specify to determine . If we equate the derivatives of the exact and HBIM expressions for and set , we find Substituting for through (2.6) then shows . This is the approach suggested by Braga et al. [28]. The entropy minimisation method of Hristov [11] involves matching exact and approximate expressions for at . For a fixed boundary temperature, this is equivalent to matching and so Hristov obtains the same 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 is no 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 , given after (2.12), (2.13), shows that does not satisfy the heat equation at and leads to an infinite value of at .

3.2. Global Matching

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

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 ). This is perhaps clearest in Section 4.2, where the minimisation technique has an error half that of the CIM. However, in cases where , due to the difficulty of implementation, the minimisation technique is used with . The combined method deals more easily with and gives more accurate solutions in this situation.

3.3. Comparison of Results

In Figure 2 we present the variation of for the HBIM and RIM formulations for . For larger 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 , with an error . The most accurate HBIM approximation requires and has . 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 and so the CIM matches the classical method of Goodman, with an error . The other common value chosen for the exponent, , leads to for HBIM and RIM, respectively. The solution obtained by matching at has errors 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 that minimises the error is close to 2 and so the improvements (in this case by taking or 2.074) from the classical approach of Goodman are small. Since the gradient of is small near , 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 . However, in subsequent sections we will consider different boundary conditions and then find that in general is a poor choice.

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 replaced by the constant flux condition . This problem has the exact solution The boundary immobilisation coordinate is again , but we now set , chosen so that the polynomial approximation is independent of . Note, for the fixed temperature boundary condition we set whereas here we set . The difference is that now as , and this scaling permits . In the previous case at the boundary and no rescaling was necessary. The heat equation is given by and since , this reduces to The function is therefore given by rearranging  (4.4) Again the HBIM is the integral of over , the RIM involves the double integral and the error is the integral of .

The HBIM formulation gives and Similarly the RIM gives and In both cases the error is independent of time and . The corresponding expression using Langford’s method has an error and so, as in the constant temperature problem, it blows up as .

Results are shown in Figure 3 for . The minimum value of the errors occurs at for the HBIM and RIM, respectively, and then . For this situation matching the approximate and exact solutions at works well, leading to , see [28], with errors only slightly above the minimum values, . We do not show the error for 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, , gives an HBIM error twice that of the minimum value while the RIM error is quadrupled. The two curves cross at , which is the CIM prediction, with a 50% increase over the minimum error.

4.2. Newton Cooling

Langford wisely avoided this problem, which does not lead to an analytical expression for the error. The boundary condition is now replaced by the cooling condition . Again this problem has an exact solution Setting and leads to the polynomial profile Note, if we assume a constant , then for small times , and so , which means that tends to the constant flux solution, (4.2). For large , and behaves like the fixed boundary temperature solution. The two limits show a different dependence and we must conclude is also a function of time.

In this case, after integrating (4.3), it follows that the HBIM and RIM formulations are given by We define by (4.5) and then the error () is

Although this appears rather daunting, all terms can be integrated analytically. For the CIM it is relatively simple to calculate and from the two first-order ODEs that result from (4.10) (the initial conditions are , , where the value comes from the constant flux solution). In this case since , we find is time dependent, but there is no initial singularity ().

To avoid the difficulties of analysing a time-dependent , in [22] the solution was calculated with . This value was chosen since, due to the singular nature of , the highest error occurred at . In Figure 4 we compare the errors over time for the CIM solution and the HBIM and RIM solutions that use a constant , for . For small times the HBIM and RIM solutions, with , respectively, show the smallest error, but as the optimal value of changes, the accuracy is lost. The CIM, which has a time-dependent , quickly becomes more accurate than the RIM approximation and around it becomes more accurate than the HBIM.

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, , 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 and this provides the initial condition for the subsequent state.

When the de-icing system is turned on, we require a heat flux at and the problem is governed by where we have taken our length-scale from the heat flux condition, rather than the ice thickness, hence . 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 is unaffected by the temperature at . However, with this assumption we may find an exact solution using separation of variables 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 and the boundary condition at is replaced by two conditions , , 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 The heat balance integral is and this reduces to which is independent of . Applying gives .

To carry through the new error analysis, we would again set , but to allow a temperature close to −1 in the vicinity of at , we choose . This gives a polynomial approximation independent of The new relation between and 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 is defined by (4.5). The obvious conclusion is that since this problem has a constant flux boundary condition, the values of that minimises are those calculated previously, namely, for the HBIM and 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 . Thus, The exact solution requires solving the nonlinear equation Using parameter values from aircraft icing models [3, 30, 31], we find , which then gives (eight terms in the sum is enough to obtain a value which does not change up to ). The HBIM and RIM predictions of , found from (5.9), are and , respectively, giving errors of 0.76% and 0.79%. The CIM is slightly worse than both of these methods, with 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 . The position of at each of these times is denoted by a circle. We only show the HBIM solution with 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.

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 Section 6.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 and this has exact solution where satisfies the transcendental equation . The standard HBIM approximating function is where we assume and are constant [12, 20]. Then the Stefan condition in (6.2) may be integrated immediately to give . The function , defined in (2.7), is given here as From (6.5) we see that can only be zero if (since we require ): this specifies a linear temperature and in general provides a poor approximation. At we find that if (and ). For , is finite but never zero. The error defined by Langford is as usual singular at , .

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 , rather than . This transforms the governing equations (6.1), (6.2) into and (6.6) is equivalent to (2.15) with replaced by . The approximating profile is Since is constant, it follows that and so (equivalent to (2.18)) The error is defined as or Note, to derive this equation, we have used the fact that (from (6.8b)) and so it is clear that is independent of time (whereas is once again singular at ). However, it is clear that the error is dependent on the value of and in fact we find the optimal 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 . Substituting for it is easy to show that and will satisfy the same expressions as for the nontransformed system, that is, solves the nonlinear equation and satisfies further details are given in [20].

In Figure 6 we present the variation of for the HBIM and RIM formulations for and two values of . With we see that the most accurate approximation comes from using the RIM formulation with , which has an error of 0.021. The most accurate HBIM approximation requires and has a slightly smaller error of . The CIM predicts , with a much higher error of . Similarly, when , we find that and for the HBIM, RIM, and CIM formulations, respectively. Also shown are solutions for the HBIM and RIM with and the solution found by matching at to the exact value (found from (6.3)).

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 . It is interesting to note that when we consider the error in , for , 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., ), 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    If the constraint holds then for early times instead of the conditions in (6.17a) and (6.17b) we have 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

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 we first need to define the boundary immobilisation coordinate which transforms the moving domain into a fixed one . If we denote then the system (6.15)–(6.18) becomes with The polynomial approximation, (6.21), is now Although the profile is independent of and we cannot conclude that in (6.23). As discussed in [33], the conditions (6.16) and (6.17a) imply that the gradient is infinite as . Thus at small times the gradient term in (6.17b) must be balanced with the left hand side, and this implies that , . However, as time increases, the constant term in (6.17b) becomes important indicating a shift in solution form and so the possibility that is time dependent.

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 , we choose to first integrate over . Then, on integrating again between and changing the order of integration, we find After substituting the profile (6.25) into the integral equations (6.27) and (6.28) we have the pair of ODEs These are coupled with the Stefan condition in (64) and this allows us to determine , , and . We know that and the initial condition for is found by considering the limit as . Setting , , with , and substituting into (64), (6.29) leads to three algebraic equations to solve for . This leads to We can eliminate and to give an explicit expression for in terms of : Note, only the positive root of the quadratic equation ensures . In fact, to ensure requires : for the CIM breaks down (note, most applications of this model involve values of [3234]).

In Figure 7 we have plotted and against for . Comparisons are shown between a numerical solution, using the Keller box scheme [35], the HBIM and RIM (with chosen to minimise at ) and the CIM. The CIM is far more accurate here because is allowed to vary with time. In the right plot we compare the constant value of used in the HBIM and RIM (which are superimposed on this plot) with the time-dependent from the CIM. We see that varies significantly with time and thus using the time dependent value is much more accurate. As increases, the variation in with time is less pronounced in the CIM and then the HBIM and RIM solutions improve.

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 function , it was also shown that the standard polynomial approximation will never satisfy the heat equation at and it is only satisfied at when . 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.

Acknowledgments

S. Mitchell acknowledges the support of the Mathematics Applications Consortium for Science and Industry (MACSI, http://www.macsi.ul.ie/) funded by the Science Foundation Ireland Mathematics Initiative Grant 06/MI/005. The research of T. Myers was supported by a Marie Curie International Reintegration Grant Industrial applications of moving boundary problems Grant no. FP7-256417 and Ministerio de Ciencia e Innovación Grant MTM2010-17162.