A novel revision of Goodman’s profile and its application to a one-phase Stefan problem

Mechanics Research Communications 38 (2011) 456–462

Contents lists available at ScienceDirect

Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom

A novel revision of Goodman’s profile and its application to a one-phase Stefan problem O.P. Layeni a,c,∗ , A.M. Adegoke b a b c

Center for Research in Computational and Applied Mechanics, University of Cape Town, Rondebosch 7701, South Africa Department of Physics, Obafemi Awolowo University, Ile-Ife 220005, Nigeria Department of Mathematics, Obafemi Awolowo University, Ile-Ife 220005, Nigeria

a r t i c l e

i n f o

Article history: Received 6 April 2011 Received in revised form 27 May 2011 Available online 13 June 2011 MSC: 40c19 80A22

a b s t r a c t Due to its extreme simplicity the heat balance integral method (HBIM), introduced by Goodman for addressing transport problems, has attracted lots of attention in recent times. Specifically, the subtle objective of most of these has been that of significant improvement of its accuracy while retaining simplicity. A crucial factor in achieving this is the choice of test functions/profiles, of which holds little structure in current literature. It is shown, by elementary methods, how seemingly disjoint sets of earlier proffered profiles relate. Further, novel test profiles, due to pertinent generators, are introduced and shown to yield remarkably accurate approximate solutions to a benchmark Stefan (melting) problem. © 2011 Elsevier Ltd. All rights reserved.

Keywords: Heat balance integral method Stefan problem Moving front Generators

1. Introduction The heat balance integral method (HBIM) was introduced by Goodman (1958, 1964) to address transport problems. Suppose the approximate solutions to an initial boundary value partial differential equation is sought. The method roughly consists of the assumption of a functional profile with certain undetermined coefficients (which satisfies the prescribed boundary and initial conditions), integrating the differential equation to create an heat integral, and simultaneously obtaining both the undetermined coefficients and its solution. Due to the tractability of the method, it has found various application in the literature; thermistor model (Wood and Kutluay, 1995; Kutluay et al., 2006), approximation of temperature perturbation front in heat transfer (Kudinov et al., 2007, 2009), approximate solutions to heat transfer problems under convective boundary conditions (Braga et al., 2005; Roday and Kazmierczak, 2009a,b), application in the modeling of thermal protection systems for space vehicles (Hsiao and Chung, 1984; Braga et al., 2003), to mention a few. Various authors have since proffered modifications of the method. The main lines of improvement of the HBIM (for better

∗ Corresponding author at: Center for Research in Computational and Applied Mechanics, University of Cape Town, Rondebosch 7701, South Africa. E-mail addresses: [email protected] (O.P. Layeni), [email protected] (A.M. Adegoke). 0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.05.017

accuracy and tractability) are: the structural (which involves the framework of the integral method, see for instance Sadoun, 2006; Wood, 2001) and the prescriptive (which concerns suggesting effective approximate profile, for instance as in Mosally et al., 2002). This article concerns significant improvement in the latter. It is well known from literature on the HBIM, for instance those by Mosally et al. (2002), Hristov (2009), Mitchell and Myers (2010a) and Sadoun (2009), that test/approximate functions/profiles employed in the implementation of the technique are analytic over the interval of application, one of the most accurate being the Gaussian while the most occurent features the polynomial test profile of the form (for an (N + 1)-dimensional partial differential equation (pde), say,)

⎧ M    ⎪ ⎪ ⎨ (xi , t; aj ) ≈ (x, t) = aj (1 − y)j , ⎪ ⎪ ⎩y =

j=1

(1)

xk , ı(x1 , . . . , xk−1 , xk+1 , . . . , xN , t)

where k ∈ {2, . . ., N}, M ∈ Z+ , x = (x1 , . . . , xN )T , ı(·) is a ‘penetration depth’ functional and (x , t) is its exact solution. The (1 + 1)dimensional version of Eq. (1) which takes the form

(x, t; aj ) =

M 

aj 1 −

j=1

x ı(t)

j

,

is one which has been largely exploited in these articles.

(2)

O.P. Layeni, A.M. Adegoke / Mechanics Research Communications 38 (2011) 456–462

Though most practical cases are of spatial dimensions greater than 1, it is understood that these restrictions (to (1 + 1)-dimension) no doubt give significant inkling into the form of solutions to higher dimensional problems. This contribution is aimed at systematically improving the HBIM prescriptively. The applicability of our technique is illustrated and benchmarked with the (1+1)-dimensional one-phase Stefan (melting) problem ∂ ∂2  = ∂t ∂x2

 ≡ (0, ı(t)) × [0, T ]

in

(3)

with boundary and initial conditions ∂ dı = −ˇ , dt ∂x (ı, t) = 0,

x = ı(t), (4)

(0, t) = f (t) = / 0, (x, 0) = 0,

x > 0,

where  is the non-dimensional temperature, ı(t) is the location of the moving front, ˇ is the reciprocal of the Stefan number and f ∈ C1 [0, T]. An observation of the original Goodman’s approximate profile (for the above Stefan problem); M x j 

(x, t; aj ) = f (t)

aj

j=0

ı(t)

,

(5)

shall concretize our discussion by considering Stefan problem (3) and (4). Firstly, in our re-presentation of the Goodman’s series, we reflect on the form of the ‘spectral cubic’ (spherical) trial function (Layeni and Johnson, in preparation),



2 (x, t; 3, a ) = f (t)

1−

x ı(t)



+ a

x ı(t)





x2 2

ı(t)

−1

;

 ⎧ 1 ⎨ −∞, 2 × (0, 1) or (7)

⎪ x    a , ⊂



 ı(t) ⎪ ⎩ 12 , ∞ × 0, − 21 + 21 4 +aa ⊆ 12 , ∞ × (0, 1) . which though lacks a quadratic term, yields significantly more accurate HBIM results than the Taylor cubic due to Eq. (5). Consequently, we select A(j) = 2j + 1 for Eq. (6) and, up to the addition of a constant, propose a class of approximate profiles for Stefan problem (3) and (4) which are consequent of the following form. Heat Profile 1

⎛

1 (x, t; cj (p1 , p2 , · · · , pk )) = f (t) ⎝1 −

⎞ ∞ x 2j−1  ⎠ cj

j=1

ı(t)

;

(8)

∞

c = 1, and the pk are elements of a such that {cj }j≥1 ∈ B0 with j=1 j set of real parameters which cardinality are according as the desired application. 2.1. Revision

shows that a0 = 1 and there exists a convex combination of coefM a = 1. It is well known that this convexification is ficients j=1 j not sufficient to express the arbitrary constants aj , j ∈ {1, . . ., M} in terms of just one parameter when M > 2, that is when polynomials (in x/ı(t)) of degree greater than quadratics are employed. Hence, extra/artificial conditions need to be contrived from the dynamics of the problem. These, more often than not, worsen the accuracy of approximations obtained. A pertinent question then is whether there are means of retaining the simplicity of the Goodman profile (11), while simultaneously significantly increasing the accuracy of the HBIM and avoiding the injection of artificial conditions? This note proffers an answer in the affirmative. The approach taken is that of the deployment of an infinite series re-presentation of Goodman’s approximate profile in the form

⎧ ∞  ⎪ ⎪ ⎨ 1 (y(xi ), t; cj ) = f (t) cj yA(j) , ⎪ ⎪ ⎩y =

457

It is observed that profile (8) satisfies the initial and boundary condition of Stefan problem (3) and (4). The terms of the sequence {cj }n≥1 are expected to depend on parameters pk . With respect to the Stefan problem (3) and (4), we require one such parameter p1 = p, say. The next pertinent issue is the choice of the {cj (p)}j≥1 . It is well known that the linear heat equation (3) has a solution of the form (Rosenbloom and Widder, 1959), (x, t) =

∞  j=0

[j/2]  xj−2k k=0

(6)

j=0



where cj ∈ R, A : Z+ ∪ {0} → Z+ ∪ {0}, cj

j≥1

∈ B0 , the space of con-

vergent real sequences such that lim j→∞ cj = 0. It is observed that with this formulation, the problem of seeking optimal index of fine Taylor or Puiseux trial series, as considered for instance in Mitchell and Myers (2010b) is avoided. The utility of this prescription is evinced by highly accurate result, for Stefan problem (8) and (9) under both a constant boundary condition and the time-dependent condition f(t) = exp (t) − 1, which are atimes two orders of magnitude more accurate than the Gaussian. This note is arranged as follows: Section 2 provides a brief motivation for, and introduces the sequential approach taken, Section 3 details the numerical examples while Section 4 gives concluding statements. 2. Revision and formulation for a Stefan problem As earlier mentioned, the prescription in this note is one of the form (6). Without loss of generalization to higher dimensions, we

zj j!

(9)

with

vj = j!

xk , ı(x1 , . . . , xk−1 , xk+1 , . . . , xN , t)

vj (x, t)

tk , (j − 2k)! k!

and

z ∈ C.

(10)

A reflection on this lends us to propose cj = a(p;j)(pj−1 /(j − 1)!), p ∈ R− , corresponding to which is the class of normalized generators ∞  j=1

a(p; j)

pj−1 . (j − 1)!

(11)

In what follows, we shall consider the simple cases wherein a(p;j) takes forms like exp (− p) or (1 − p)(j − 1)! which, respectively, lead to Poisson or geometric generators. It is clear that a Poisson generator induced by exp (− p) leads to a profile of the form

⎛

1 (x, t; cj ) = f (t) ⎝1 −

∞ 

exp(−p)

j=1

pj−1 (j − 1)!

⎞ x 2j−1 ⎠ ı(t)

,

(12)

while a geometric generator evoked by (1 − p)(j − 1)! yields

⎛

1 (x, t; cj ) = f (t) ⎝1 −

∞ 

⎞ x 2j−1 j−1 ⎠

(1 − p)p

j=1

In summary, we note the following:

ı(t)

.

(13)

458

O.P. Layeni, A.M. Adegoke / Mechanics Research Communications 38 (2011) 456–462

Theorem 1. The Poisson and geometric generators transform the sequential form (8), respectively, to the Gaussian approximating profile of Mosally et al. (2002)



f (t)

x exp ı(t)

1−

  p

−1 +



x2



p ∈ R− ∪ 0

,

2

ı(t)

(14)

Assumption 1. Suppose there exists an interval I ⊂ R− such that the power sequence {cj (p)}j≥1 ∈ B0 is such that for functions I × Z+ (p, j) → H, J(cj (p)) ∈ R/{0} defined by H(cj (p)) =

∞  2

and inverse multi-quadric (IMQ) of the form f (t)



x 1− ı(t)





x + ı(t)

x2 − ı(t)2



J(cj (p)) = 2 − ε = p−1 ∈ R− .

,

x2 − εı(t)2

1 + (ε − 1)

x ı(t)

−ε +

x2

.

x ı(t)

 −(1 + b) +

x2

−1 (17)

ı(t)2

x a+b ı(t)

 c+

x2

−1 

ı(t)2

,

a, b, c ∈ R, (18)



for Stefan problem (3) and (4). Corollary

1.

The

{pj /((−1 + exp(p))j!)}

j≥1

index

sequence

{cj }j≥1 =

yields the alternative Gaussian approximate

f (t)

2

2

(−1 + e(px )/(ı(t) ) )ı(t) 1− (−1 + ep )x



j−1 cj (p),

Let {cj (p)}j≥1 ∈ B0 be a power sequence satisfying

f (t)2

> 0,

(21)

then the HBIM gives the position of the moving front for Stefan problem (3) and (4) as

⎛ ⎞ ∞    (1 − 2j)cj ⎜ ⎟  t  ⎟ j=2 2 ⎜ ⎜ ⎟× ı(t) = f (s) ds,  ∞ ⎟ f (t) ⎜  0 ⎝ ⎠  2− −1 j cj

(22)

j=1

with ı(0) = limt→0 ı(t). Noting that

⎧ l = 0, ⎪ ⎨ c1∞ ∂1 (x, t; cj ) f (t)  × (x = l, t) = − (−1 + 2j)cj l = ı(t), ∂x ı(t) ⎪ ⎩

(23)

j=1

.

(19)

It is worth observing that the profiles in Theorem 1 and its corollary proffer one-parameter profiles which are analytic over R+ , and profile (19) satisfies Eq. (4)2 in the limiting sense. In applications, these can be trivially extended to suit arbitrary initial and boundary values considered. 2.2. HBIM formulation for Stefan problem We take the most effective approach in the application of the HBIM (Wood, 2001). This entails integrating both sides of Eq. (3) to the spatial variable x over the moving interval with respect  0, ı(t) to obtain ∂ ∂ d (x = ı(t)) − (x = 0) = dt ∂x ∂x

g(t) − g(0)

Proof. shifted

profile



,

Assumption 1. Suppose further that f ∈ C1 [0, T] is positive and monotone increasing with primitive g such that

t→0

which follows also derives directly from assuming a three parameter IMQ approximating heat profile of the form (x, t; a, b, c) = f (t)

Theorem 2.

lim



cj (p)

2

See Appendix A for illustration of this with cj being (1 − p)p−1+j , j(p − 1)2 p−1+j , p−1+j /(ep (− 1 + j)!) or pj /((−1 + ep ) j!).

(16)

Fix b = ε − 1, then the expression 1+b



sgn(H(cj (p))) = sgn(J(cj (p))) for all p ∈ I.

−1

ı(t)2

∞ 

 ∞ 

1

(15)

Proof. The proof follows straightaway. Alternatively, employing a symbolic mathematical portal like MATHEMATICA easily verifies that Eqs. (12) and (13) sum to, respectively, to Eqs. (14) and (15). The fact that Eq. (14) is an IMQ also follows simply from re-expressing the bracketed term as



d cj (p) − 2p dp



and



ı(t)

⎛ 1 (x, t; cj ) dx = ı(t)f (t) ⎝1 −

∞  cj

0

j=1

2j

⎞ ⎠,

(24)

we have, on observing Eq. (20) and taking some elementary steps, that the position of the moving front and the Neumann condition Eq. (4)1 are, respectively, given by

⎛ ⎞ ∞    (1 − 2j)cj ⎜ ⎟  t  ⎟ j=2 2 ⎜ ⎜ ⎟× ı(t) = f (s) ds,  ∞ ⎟ f (t) ⎜  0 ⎝ ⎠  2− −1 j cj

(25)

j=1

ı(t)

 dx.

(20)

and

0

Resolution of the time derivative of the heat integral (the right hand side of Eq. (20)), requires that the Stefan condition (that is the Neumann boundary condition Eq. (4)1 ) be applied. The approximate temperature profile, in the form of a series is substituted into Eq. (20), following which ı(t), and the relationship of its undetermined parameter(s) with the inverse Stefan number ˇ are determined. The only ingredients necessary for these are the boundary and initial boundary conditions together with the condition ı(0) = 0 (which is employed in a limiting sense in Theorem 2 below).

2

t

ˇ(f (t) − 2(

0

f (t)

(2 −

∞ 

∞ 

j−1 cj )

(−1 + 2j)cj



f (s) ds)f (t)) 4

=

j=1

j=1 ∞ 

2

≡ A(cj (p)).

(26)

(1 − 2j)cj

j=2

The assumptions on the power series {cj }n≥1 and the positivity of the integral guarantees that ı(t) makes sense and hence the proof. 

O.P. Layeni, A.M. Adegoke / Mechanics Research Communications 38 (2011) 456–462

459

Table 1 Generator–profile–sequence–restriction table. S. no.

Profile

Generator

cj (p)

A(cj (p))

Restrictions

P1

IMQ

Geometric

(1 − p)p−1+j

2(−3+p)p2 (−1+p)(−2p+(−1+p) log(1−p))

p ∈ (−3.92155, 0)

P2

–

”

j(− 1 + p)2 p−1+j

2p(6+(−3+p)p) −1+p2

p ∈ (−1.00000, 0)

Poisson

p−1+j ep (−1+j)!

”

pj (−1+ep )j!

P3

Gaussian

P4

Alt. Gaussian

−

2p(−1+ep (1+2p))

p ∈ (−1.25643, 0)

1+ep (−1+2p)

2(−1+ep (1−2p)+p)

p ∈ (−3.98124, 0)

−2+2ep ++[0,−p]+log(−p)

The relationship between the reciprocal Stefan number ˇ and the parameter p of a power series {cj (p)}, can be sieved from Eq. (26) for specific boundary condition f(t), following which the approximate moving front ı(t) positions and subsequently temperature profiles 1 (x, t;cj (p)) are determined. Table 1 shows profile–series–parameter relationships and restrictions for some selected power series. Suppose {cj (p)}j≥1 is as assumed in Theorem 2, it is easily

Absolute error in Θ

0.0015

0.0010

observed that a partial sum of Heat Profile 1 up to a cubic in x/ı(t), that is



f (t)

1−

 2 

cj (p)

x ı(t)

2j−1 



≡ f (t)

x x3 1 − c1 (p) − c2 (p) 3 ı(t) ı(t)

0.0005

(27)

j=1

0.1

2

together with the condition j=1 cj (p) = c1 (p) + c2 (p) = 1, is no different from the spherical profile Eq. (7). The structure of Heat Profile 1 permits its accuracy to monotonically increase with increasing partial sums, which consequently lifts the burden of seeking optimal indices (which is a prevalent theme in the application of Goodman’s profiles and its type in the literature) while implementing the HBIM. 3. Numerics for the Stefan problem A comparison of the accuracy of the results due to test profiles Pj , j ∈ {1, 2, 3, 4} is detailed in the following. We shall consider cases both constant and time-dependent boundary conditions f(t) = 1 and f(t) = exp(t) − 1 respectively.

∞ 

(2 −

∞ 

j−1 cj )

j=1

ˇ=



ı(t)





2

(28)

(1 − 2j)cj

establishes the relationship between parameter p and the reciprocal Stefan number ˇ. A log of these relationship for profiles Pj is furnished in Table 2. A comparison of the L2 errors of moving front approximation ı(t) due to the class of profiles Pj introduced above is given in Table 3 while illustrative absolute error plots are provided in Figs. 1 and 2. The L2 errors of the temperature distribution estimates are over the domain [0, ı(t)) × [0, T] and in the sense



T



!

ı(t)



T



||2 dx dt − 0

∂ d ∂t



dx = (x = ı(t), t) − (x = 0, t) − ı(t)

0

T = 1 is taken in Table 3.



0

ı(t)

d (x, t) dx − dt



∂ (x = 0), ∂x

(29)

ı(t)

x(x, t) dx = −(x = 0, t) 0

∂ (x = ı(t), t). ∂x

(30)

Absolute Error in Θ 0.000035

j=2

!

x

0

−ı(t)

,

∞

x

In this time dependent boundary value case, we shall determine approximate position ı(t) of the moving front through the refined integral method (RIM). The RIM was introduced by Sadoun (2006) and it has the advantage of being able to estimate the position of the front and temperature distribution more accurately than the HBIM in a reasonable neighbourhood of the interface while being less accurate at estimating the temperature field elsewhere. We integrate Eq. (3) twice to have

d ı(t) dt

(−1 + 2j)cj

j=1

0.4

3.2. f(t) = exp(t) − 1

leading to

In this case, Eq. (26) which reads off as

0.3

Fig. 1. HBIM  absolute error plots at ˇ = 1, t = 0.1; Plain line - P1, Dashed line - P2, Thick line - P3, Dotted line- P4.

0

3.1. f(t) = 1

0.2

0

0.000025 0.00002 0.000015 0.00001 5. x 10- 6

ı(t)

|app |2 dx dt . 0

0.00003

0.02

0.04

0.06

0.08

0.10

0.12

0.14

x

Fig. 2. HBIM  absolute error plots at ˇ = 10, t = 0.1; Plain line - P1, Dashed line P2, Thick line - P3, Dotted line- P4.

460

O.P. Layeni, A.M. Adegoke / Mechanics Research Communications 38 (2011) 456–462

Table 2 ˇ − p relationship and selected (ˇ, p) coordinates for profiles Pj when f(t) = 1. S. no.

Profile

ˇ(p)

p, ˇ = 1

p, ˇ = 10

P1

IMQ

2p(1+p)−(−1+p2 ) log(1−p) 2(−3+p)p2

−0.130248

−0.016178

P2

–

(1+p)(1+3p) − 2p(6+(−3+p)p)

−0.061755

−0.008035

P3

Gaussian

−0.243123

−0.032068

−0.117369

−0.015962

P4

(1+ep (−1+2p))(−2+2ep ++[0,−p]+log(−p)) 2(−1+ep )(−1+ep (1−2p)+p)

−1 +

Alt. Gaussian

1 2p

−

2 −1+ep (1+2p)

Table 3 L2 error comparison of profiles Pj s HIBM moving front and temperature approximations, respectively, over [0, 1] and [0, ı(t)] × [0, 1]. ˇ

Spectral cubic

P1 (IMQ)

P2

P3 (Gaussian)

P4 (Alt. Gaussian)

1 10

HBIM: ı(t) errors 4.3439 × 10−3 2.7749 × 10−3

3.2328 × 10−4 2.5983 × 10−6

9.1826 × 10−4 4.1881 × 10−6

2.1083 × 10−3 1.0938 × 10−5

1.2761 × 10−3 6.4099 × 10−6

1 10

HBIM: (x, t) errors 1.6914 × 10−5 7.5987 × 10−7

1.0759 × 10−5 9.8403 × 10−8

1.0769 × 10−4 1.3306 × 10−7

1.4421 × 10−5 3.5314 × 10−7

9.9532 × 10−6 2.0891 × 10−7

Applying (4)1 and (4)3 to Eq. (30), we obtain d dt



ı(t)

x(x, t) dx = (x = 0, t) − ˇı(t) 0

dı(t) , dt

Proof.

(31)

and subsequently (upon reflection of the time dependent boundary condition) d dt



ı(t)

x(x, t) dx = f (t) − ˇı(t) 0

dı(t) . dt

(32)

Form (32) together with the Neumann boundary condition Eq. (4)1 establish the relationship between the reciprocal Stefan number ˇ and the sequence parameter p in time. Further, these determine the position of the moving front ı(t) and subsequently the temperature distribution across the domain. Theorem 3. Suppose there exists an interval I ⊂ R− such that {cj (p)}j≥1 ∈ B0 is a power sequence for which

Form (32) gives the ordinary differential equation in ı(t)2

  ⎧ ˇ d d 2 ⎪ + G(cj ; j) (ı(t)2 ) = f (t), (ı(t) ) f (t)G(cj ; j) + ⎪ ⎪ 2 dt dt ⎪ ⎨ ⎛ ⎞ ∞ ⎪ 1  cj ⎪ ⎝ ⎠ , ı(0) = 0, ⎪ ⎪ ⎩ G(cj ; j) = 2 − 2j + 1 j=1

which solution is Eq. (33). For ı(t) so defined to make sense, the ∞ restriction on series (c (p))/(2j + 1) is necessary.  j=1 j The Neumann boundary condition ∂ dı(t) = −ˇ dt ∂x

x = ı(t)

at

which clearly yields

 f (t)

∞ 

ˇ + f (t) − 2f (t)

 cj 1+2j

+

 t 0

f (t) dt

j=1

∞  cj (p) j=1

1 ≤ 2 2j + 1

 for all

(34)

p ∈ I.

∞ 

ˇ + f (t) − 2f (t)



∞ 

−1 + 2

 cj

d f (t) dt

1+2j j=1

2 cj 1+2j

j=1 ∞  (−1 + 2j)cj

Suppose further that f ∈ C1 [0, T] is positive and monotone increasing, then the RIM gives the position of the moving front for Stefan problem (3) and (4) as

  t  2  ⎞ ⎞ ı(t) =  ⎛⎛ f (t) dt.  ∞ 0   cj ⎠ f (t) + ˇ⎠  ⎝⎝1 − 2 2j+1

(33)

j=1

= f (t)

ˇ

,

(35)

j=1

relates ˇ and p in time via the boundary condition f(t). Appendix B lists the resultants of Eq. (35) for profiles Pj . In the following, we consider the specific case wherein f(t) = exp (t) − 1. The exact solution for this when ˇ = 1 is given as (x, t) = exp(t − x) − 1,

ı(t) = t.

(36)

Table 4 Comparison of exact, enthalpy methods and Pj RIM moving front approximations for the one-phase Stefan problem with non-constant boundary condition f(t) = exp (t) − 1 when ˇ = 1. t

Exact

Enthalpy

P1 (IMQ)

P2

P3 (Gaussian)

P4 (Alt. Gaussian)

0.25 0.5 0.75 1 N

t

ı(t) 0.2500 0.5000 0.7500 1.0000 – –

0.2505 0.5006 0.7507 1.0007 100 0.001

0.2499 0.4992 0.7474 0.9944 – –

0.2499 0.4991 0.7470 0.9930 – –

0.2499 0.4990 0.7466 0.9917 – –

0.2499 0.4991 0.7469 0.9927 – –

O.P. Layeni, A.M. Adegoke / Mechanics Research Communications 38 (2011) 456–462

Table 4 details a comparism of moving front position ı(t) approximations obtained via the sequential RIM approach with both the exact solution and that obtained through the enthalpy method (as given in Caldwell and Kwan, 2004).

461

where  is the Euler–Gamma constant and [a, x] the incomplete Gamma function defined by



∞

t a−1 exp(−t) dt.

x

4. Conclusion

when p ∈ (− 3.98124, 0), sgn(H(cj (p))) = + 1 = sgn(J(cj (p))). In this article, we have presented an applicable formalism for profile selection for the HBIM, in the approximate solution of two-dimensional moving boundary partial differential equations, reflecting on the fundamental solutions to their state equations. This gives a framework and basis for application and extension to higher dimensional problems. The utility of our characterization is demonstrated in its application to a Stefan melting problem. The HBIM approximations obtained via the suggested profiles are remarkable, with the inverse multi-quadric HBIM being at least of two orders of magnitude more accurate than those in current literature. In fact the implementation of the Gaussian profile re-introduced herein, albeit with a different motivation, is easier implemented than as in Mosally et al. (2002). The reliability of our approach is further elicited in the timedependent boundary condition case, where accurate RIM moving front position estimates obtained are comparable with those due to spatio-temporal discretization methods in the literature. An immediate area of extension of application of these is in the approximate solutions to steady two dimensional laminar flows of Non-Newtonian fluids, as a potent replacement for the Pohlhaussen’s integral method, the penetration depth functional ı in this case determining the boundary layer thickness.

Appendix B. In the cases wherein profiles P1 , P2 , P3 and P4 are considered, Eq. (35) respectively relates ˇ and p in terms of f(t) and its derivatives according as Eqs. (B.1)–(B.4) viz: −

√ −1 √ f (t)(2r 5/2 ˇ − (1 + r)((−2 + r) r − 2(−1 + r) tanh ( r))f (t)) (−1 + r)ˇ



√ −1 √ r 3/2 ((−2 + r) r − 2(−1 + r) tanh ( r))(

f (t) dt)

(B.1)

d f (t) dt

0

+

√ −1 √ r 3/2 ˇ − ((−2 + r) r − 2(−1 + r) tanh ( r))f (t)

= 0,

√ 2 −1 √ 2 r 3/2 (r 3/2 ˇf (t) + ( r(−1 + 2r) + (−1 + r) tanh ( r))f (t) ) √ √ 2 2 −1 (r 3/2 ˇ + ( r(−1 + 2r) + (−1 + r) tanh ( r))f (t)) =

+

t √ 2 −1 √ f (t) dt)((d f (t))/(dt))) r 3/2 (( r(−1 + 2r) + (−1 + r) tanh ( r))( 0

√ 2 2 −1 √ (r 3/2 ˇ + ( r(−1 + 2r) + (−1 + r) tanh ( r))f (t)) (1 + 3r)f (t) , ˇ − rˇ

2p5/2 f (t) +

Acknowledgment

t

(B.2) (B.2)

√ √ (1 + 2p)((−1 + p) p + G( p))f (t)2 ˇ

t √ √ f (t) p3/2 ((−1 + p) p + G( p))( 0 f (t) dt) d dt =− √ √ p3/2 ˇ + ((−1 + p) p + G( p))f (t)

The first author holds an NRF SARChI postdoctoral fellowship at the University of Cape Town, funded through the South African Research Chair in Computational Mechanics.

(B.2)

(B.3)

and Appendix A.

√ √ √ √ 2 (−1 + ep ) p((−1 + ep ) pˇf (t) + ((1 + ep ) p − 2ep G( p))f (t) ) √ √ √ 2 p p p ((−1 + e ) pˇ + ((1 + e ) p − 2e G( p))f (t))

Illustration for Assumption 1.

t √ √ √ (−1 + ep ) p(((1 + ep ) p − 2ep G( p))( f (t) dt)(d f (t))/(dt)) 0 = √ √ √ 2 ((−1 + ep ) pˇ + ((1 + ep ) p − 2ep G( p))f (t))

(i) cj (p) = (1 − p)p−1+j . H(cj (p)) = −p,



J(cj (p)) = 2 + −1 +

1 p



(A.1)

log(1 − p).

when p ∈ (− 3.92155, 0), sgn(H(cj (p))) = + 1 = sgn(J(cj (p))). (ii) cj (p) = j(− 1 + p)2 p−1+j .



H(cj (p)) = 2 +

.

−1 + exp(−p) , p

(A.3)

J(cj (p)) = exp(−p) (−1 + exp(p) + 2p) . when p ∈ (− 1.25643, 0), sgn(H(cj (p))) = + 1 = sgn(J(cj (p))). (iv) cj (p) =

pj . (−1+ep )j!

H(cj (p)) =



1 − p + exp(p) −2 + exp(p) + p − 2p2

J(cj (p)) = 2 +

(−1 + exp(p))

2

 + [0, −p] + log(−p) , −1 + exp(p)

x

exp(y2 ) dy.

(B.5)

0

References

when p ∈ (− 1, 0), sgn(H(cj (p))) = + 1 = sgn(J(cj (p))). ep (−1+j)!

(B.4)

where G is the Dawson integral defined by

(A.2)

J(cj (p)) = 1 + p. (iii) cj (p) =

(1 + ep (−1 + 2p))f (t) , (−1 + ep )ˇ

G(x) = exp(−x2 )

H(cj (p)) = p(−2 + 3p),

p−1+j

+

 , (A.4)

Braga, W.F., Mantelli, M.B.H., Azevedo, J.L.F., June 2005. Analytical solution for onedimensional semi-infinite heat transfer problem with convection boundary condition. In: AIAA Thermophys. Conference. Braga, W.F., Mantelli, M.B.H., Azevedo, J.L.F., June 2003. Approximate analytical solution for one-dimensional ablation problem with time-variable heat flux. In: AIAA Thermophys. Conference. Caldwell, J., Kwan, Y.Y., 2004. Numerical methods for one-dimensional Stefan problems. Commun. Numer. Meth. Eng. 20, 535–545. Goodman, T.R., 1958. The heat-balance integral and its application to problems involving a change of phase. Trans. ASME 80, 335–342. Goodman, T.R., 1964. Application of integral methods to transient nonlinear heat transfer. Adv. Heat Transfer 1, 51–122. Hristov, J., 2009. The heat-balance integral method by a parabolic profile with unspecified exponent: analysis and benchmark exercises. Therm. Sci. 13 (2), 27–48. Hsiao, J.S., Chung, B.T.F., January 1984. A heat balance integral approach for twodimensional heat transfer in solids with ablation. In: AIAA Aerospace Sciences Meeting.

462

O.P. Layeni, A.M. Adegoke / Mechanics Research Communications 38 (2011) 456–462

Kudinov, V.A., Stefanyuk, E.V., Antimonov, M.S., 2007. Analytical solutions of the problems of heat transfer during liquid flow in plane-parallel channels by determining the temperature perturbation front. J. Eng. Phys. Thermophys. 80 (5), 1038–1049. Kudinov, V.A., Stefanyuk, E.V., Antimonov, M.S., 2009. Analytical solution method for heat conduction problems based on the introduction of the temperature perturbation front and additional boundary conditions. J. Eng. Phys. Thermophys. 82 (3), 537–555. Kutluay, S., Wood, A.S., Esen, A., 2006. A heat balance integral solution of the thermistor problem with a modified electrical conductivity. Appl. Math. Model. 30 (4), 386–394. Layeni, O.P., Johnson, J.V. Chebyshev spectral polynomial heat balance integral methods with application to a one-phase Stefan problem, in preparation. Mitchell, S.L., Myers, T.G., 2010a. Application of standard and refined heat balance integral methods to one-dimensional Stefan problems. SIAM Rev. 52 (1), 57–86. Mitchell, S.L., Myers, T.G., 2010b. Improving the accuracy of heat balance integral methods applied to thermal problems with time dependent boundary conditions. Int. J. Heat Mass Transfer 53, 3540–3551.

Mosally, F., Wood, A.S., Al-Fhaid, A., 2002. An exponential heat integral balance method. Appl. Math. Comput. 130, 87–100. Roday, A.P., Kazmierczak, M.J., 2009a. Analysis of phase-change in finite slabs subjected to convective boundary conditions: Part I. Melting. Int. Rev. Chem. Eng. 1 (1), 87–99. Roday, A.P., Kazmierczak, M.J., 2009b. Analysis of phase-change in finite slabs subjected to convective boundary conditions: Part II. Freezing. Int. Rev. Chem. Eng. 1 (1), 100–108. Rosenbloom, P.C., Widder, D.V., 1959. Expansions in terms of heat polynomials and associated functions. Trans. Amer. Math. Soc. 92, 220–266. Sadoun, N., 2006. On the refined integral method for the one-phase Stefan problem with time-dependent boundary conditions. Appl. Math. Model. 30, 531–544. Sadoun, N., 2009. On the Goodman heat-balance integral method for Stefan likeproblems: further considerations and refinements. Therm. Sci. 13 (2), 81–96. Wood, A.S., Kutluay, S., 1995. A heat balance integral model of the thermistor. Int. J. Heat Mass Transfer 38 (10), 1831–1840. Wood, A.S., 2001. A new look at the heat balance integral method. Appl. Math. Model. 25, 815–824.