On the Solutions of a Stefan Problem with Variable Latent Heat

Abstract

This paper is concerned with a reliable treatment of the Stephan problem with variable latent heat. The power series solutions have been presented. The method is examined for computational efficiency and accuracy. Also, an analytical solution based on a similarity variable is presented in the case when the Dirichlet boundary condition at the water-ice interface depends on time.

1. Introduction

The solution of the Stefan problem with variable latent heat consists of finding 𝑢 and the moving melt interface 𝑠 such that 𝜕𝑢=𝜕𝜕𝑡²𝑢𝜕𝑥²,0≤𝑥≤𝑠(𝑡), (1) which is the governing equation, subject to the following Neumann boundary condition at the left end of the domain describing the inlet heat flux: −𝜕𝑢()(𝑡)𝜕𝑥𝑡,0=ℎ,𝑡>0, (2) the Dirichlet boundary condition at the water-ice interface: 𝑢(𝑡,𝑠(𝑡))=𝑢^∗,𝑡>0, (3) the Stefan condition: 𝑠𝑑𝑠𝑑𝑡=−𝜕𝑢((𝑡(0)𝜕𝑥𝑡,𝑠)),𝑠=𝑏,𝑡>0, (4) and the initial temperature distribution: 𝑢(0,𝑥)=𝜑(𝑥),0<𝑥<𝑠(0)=𝑏. (5) On comparison of this condition with the one-phase Stefan condition, we observe that the latent heat term 𝑠 is not a constant but, rather, a linear function of position [1]. Problems such as this have been treated under the name of Stefan problems with variable latent heat [1–3].

An analytical solution for this Stefan problem is presented in [1], by introducing the similarity variables 𝜉=𝑥/√𝑡 and seeking the solution in the form 𝑢(𝑡,𝑥)=𝐹(𝜉) with 𝐹 being an unknown function. Accordingly, it is natural that the interface location 𝑠(𝑡) should be proportional to √𝑡; that is, 𝑠(𝑡)=𝐴√𝑡, where 𝐴 is a constant.

The purpose of this paper is to apply the Adomian decomposition method [4–25] to find the solution of (1), (2), and (5), that is, the temperature distribution of the water 𝑢(𝑡,𝑥), and then use a direct method to determine the position of the ice-water interface as a function of time. Also, an analytical solution based on a similarity variable is presented in the case when the Dirichlet boundary condition at the water-ice interface depends on time of the form 𝑢^∗(𝑡)=𝛼√𝑡, where 𝛼 is a constant.

2. The Adomian Decomposition Method

Based on the Adomian decomposition method, we write (1) in Adomian’s operator-theoretic notation as 𝐿_𝑥𝑥𝑢=𝜕𝑢,𝜕𝑡 (6) where 𝐿_𝑥𝑥=𝜕²𝜕𝑥². (7) We conveniently define the inverse linear operator as 𝐿⁻¹_𝑥𝑥(⋅)=^𝑥₀(⋅)𝑑𝑥𝑑𝑥. (8) Applying the inverse linear operator 𝐿⁻¹_𝑥𝑥 to (6) and taking into account that (𝜕𝑢/𝜕𝑥)(𝑡,0)=−ℎ(𝑥), we obtain 𝑢(𝑡,𝑥)=𝐴(𝑡)−ℎ(𝑡)𝑥+^𝑥₀𝜕𝑢𝜕𝑡𝑑𝑥𝑑𝑥, (9) where the unknown boundary condition 𝐴(𝑡)=𝑢(𝑡,0) will be determined.

Define the solution 𝑢(𝑡,𝑥) by an infinite series of components in the form 𝑢(𝑡,𝑥)=^∞_𝑛=0𝑢_𝑛(𝑡,𝑥). (10) Consequently, the components 𝑢_𝑛 can be elegantly determined by setting the recursion scheme 𝑢₀(𝑡)(𝑡)𝑢=𝐴−ℎ𝑥,_𝑛+1=^𝑥₀𝜕𝑢_𝑛𝜕𝑡𝑑𝑥𝑑𝑥,𝑛≥0, (11) for the complete determination of these components. In view of (11), the components 𝑢₀(𝑡,𝑥),𝑢₁(𝑡,𝑥),𝑢₂(𝑡,𝑥),… are immediately determined as 𝑢₀(𝑡)(𝑡)𝑢=𝐴−ℎ𝑥,₁=𝐴^′(𝑡)𝑥²2!−ℎ^′(𝑡)𝑥³,𝑢3!₂=𝐴^′′(𝑡)𝑥⁴4!−ℎ^′′(𝑡)𝑥⁵,5!⋯. (12) Consequently, the solution is readily found to be 𝑢(𝑡,𝑥)=^∞_𝑛=0𝐴^(𝑛)(𝑡)𝑥^2𝑛()!−2𝑛^∞_𝑛=0ℎ^(𝑛)(𝑡)𝑥^2𝑛+1()!,2𝑛+1 (13) obtained by substituting (12) into (10).

We remark here that the unknown boundary condition 𝐴(𝑡) can be easily determined by using the initial condition equation (5). Substituting 𝑡=0 into (13) and using the Taylor expansion of 𝜑(𝑥) lead to ^∞_𝑛=0𝐴^(𝑛)(0)𝑥^2𝑛()!−2𝑛^∞_𝑛=0ℎ^(𝑛)(0)𝑥^2𝑛+1()!=2𝑛+1^∞_𝑛=0𝜑^(𝑛)(0)𝑥^𝑛.𝑛! (14) Equating the coefficients of like power of 𝑥 in both sides of (14) and taking into account that the compatibility conditions ℎ^(𝑛)(0)=−𝜑^(2𝑛+1)(0),𝑛≥0, (15) yield 𝐴^(𝑛)(0)=𝜑^(2𝑛)(0),𝑛≥0. (16) Thus 𝐴(𝑡)=^∞_𝑛=0𝜑^(2𝑛)(0)𝑡^𝑛(𝑛)!. (17) Accordingly, the solution equation (13) is completely determined by defining the function 𝐴(𝑡).

Once the function 𝑢(𝑡,𝑥) is obtained, we can rewrite the Stefan condition equation (4) in terms of the known function 𝑢(𝑡,𝑥) including the Dirichlet and Neumann boundary conditions. For this, integrating (1) with respect to 𝑥 from 0 to 𝑥 and taking into account that −(𝜕𝑢/𝜕𝑥)(𝑡,0)=ℎ(𝑡), we obtain 𝜕𝑢(𝑡)+𝜕𝑥=−ℎ^𝑥₀𝜕𝑢𝜕𝑡𝑑𝑥. (18) Thus (4) can be replaced by 𝑠𝑑𝑠(𝑡)−𝑑𝑡=ℎ₀^𝑠(𝑡)𝜕𝑢(0)𝜕𝑡𝑑𝑥,𝑠=𝑏,𝑡>0. (19) Using the following Leibniz rule for differentiation under the integral sign: 𝑑𝑑𝑡₀^𝑠(𝑡)𝑢(𝑡,𝑥)𝑑𝑥=₀^𝑠(𝑡)𝜕𝑢((𝑡𝜕𝑡𝑑𝑥+𝑢𝑡,𝑠))𝑑𝑠,𝑑𝑡 (20) and taking into account that 𝑢(𝑡,𝑠(𝑡))=𝑢^∗, we obtain ₀^𝑠(𝑡)𝜕𝑢𝑑𝜕𝑡𝑑𝑥=𝑑𝑡₀^𝑠(𝑡)𝑢(𝑡,𝑥)𝑑𝑥−𝑢^∗𝑑𝑠.𝑑𝑡 (21) Substituting (21) into (19), we obtain 𝑠𝑑𝑠𝑑𝑡−𝑢^∗𝑑𝑠+𝑑𝑑𝑡𝑑𝑡₀^𝑠(𝑡)𝑢(𝑡,𝑥)𝑑𝑥=ℎ(𝑡). (22) Since (𝑑𝑠²/𝑑𝑡)=2𝑠(𝑑𝑠/𝑑𝑡), we have 𝑑𝑠𝑑𝑡²−2𝑢^∗𝑠+2₀^𝑠(𝑡)𝑢(𝑡,𝑥)𝑑𝑥=2ℎ(𝑡). (23) Applying the inverse linear operator 𝐿_𝑡⁻¹(⋅)=∫^𝑡₀(⋅)𝑑𝑡 to (23), we obtain 𝑠²−2𝑢^∗𝑠+2₀^𝑠(𝑡)𝑢(𝑡,𝑥)𝑑𝑥=𝑏²−2𝑏𝑢^∗+2^𝑡₀ℎ(𝑡)𝑑𝑡+2^𝑏₀𝑢(0,𝑥)𝑑𝑥. (24) It follows that 𝑠²−2𝑢^∗𝑠+2₀^𝑠(𝑡)𝑢(𝑡,𝑥)𝑑𝑥=𝑏²−2𝑏𝑢^∗+2^𝑡₀ℎ(𝑡)𝑑𝑡+2^𝑏₀𝜑(𝑥)𝑑𝑥. (25) Substituting the solution equation (13) of the heat equation into (25), we obtain 𝑠²−2𝑢^∗𝑠+2₀^{𝑠(𝑡)∞}_𝑛=0𝐴^(𝑛)(𝑡)𝑥^2𝑛()!−2𝑛^∞_𝑛=0ℎ^(𝑛)(𝑡)𝑥^2𝑛+1()!(𝑡),2𝑛+1𝑑𝑥=𝑔 (26) where 𝑔(𝑡)=𝑏²−2𝑏𝑢^∗+2^𝑡₀ℎ(𝑡)𝑑𝑡+2^𝑏₀𝜑(𝑥)𝑑𝑥. (27) Thus 𝑠²−2𝑢^∗𝑠+𝐹(𝑡,𝑠(𝑡))=𝑔(𝑡), (28) where 𝐹(𝑡,𝑠(𝑡))=2^∞_𝑛=0𝐴^(𝑛)(𝑡)𝑠^2𝑛+1()!−2𝑛+1^∞_𝑛=0ℎ^(𝑛)(𝑡)𝑠^2𝑛+2()!.2𝑛+2 (29) We now can determine the shoreline 𝑠(𝑡) with time by solving the nonlinear equation (28). In order to demonstrate the feasibility and efficiency of this method, we consider the following case: ℎ(𝑡)=𝑒^−𝑡. (30) If we choose 𝜑(𝑥)=𝑒^𝑥, then a simple calculation leads to 𝐴(𝑡)=𝑒^𝑡, and (28) becomes 𝑠²−2𝑢^∗𝑠+𝑒^𝑡∞_𝑛=0𝑠^𝑛(𝑡),𝑛!=𝑔 (31) where 𝐹(𝑡,𝑠(𝑡))=𝑒^𝑡+𝑠.

Let us write 𝑠(𝑡) and 𝑔(𝑡) in series forms 𝑠(𝑡)=∑^∞_𝑛=0𝑠_𝑛𝑡^𝑛 and 𝑔(𝑡)=∑^∞_𝑛=0𝑔_𝑛𝑡^𝑛. Thus 𝑠(𝑡)=^∞_𝑛=0𝑠_𝑛𝑡^𝑛=^∞_𝑛=01𝑠𝑛!^𝑛𝑡^𝑛𝑛. (32) To compute [𝑠_𝑛𝑡^𝑛]^𝑛 we need the following theorem [4].

Theorem 1. If 𝑢=∑^∞_𝑛=0𝑎_𝑛𝑡^𝑛 and 𝐹(𝑢)=∑^∞_𝑛=0𝑏_𝑛𝑡^𝑛 are convergent, then 𝐹(𝑢)=∑^∞_𝑛=0𝛽_𝑛𝑡^𝑛 is convergent, where 𝛽₀=𝑏₀+^∞_𝑛=0𝑏_𝑛𝐴^𝑛₀𝑎₀,𝛽₁=^∞_𝑛=0𝑏_𝑛𝐴^𝑛₁𝑎₀,𝑎₁,𝛽₂=^∞_𝑛=0𝑏_𝑛𝐴^𝑛₂𝑎₀,𝑎₁,𝑎₂,…, (33) and 𝐴^𝑛_𝑚(𝑎₀,𝑎₁,…,𝑎_𝑚) are the Adomian polynomials.

Using this theorem with the given formula 𝑏_𝑛=1/𝑛!, we see that 𝑠(𝑡)=^∞_𝑛=0𝑠_𝑛𝑡^𝑛=^∞_𝑛=0𝛽_𝑛𝑡^𝑛,𝛽₀=^∞_𝑛=01𝑠𝑛!^𝑛0,𝛽₁=^∞_𝑛=11𝑛!𝑛𝑠₁𝑠₀^𝑛−1,…. (34)(35) Consequently, we obtain the recurrence relations for the coefficients 𝑐_𝑛−2𝑢^∗𝑠_𝑛+𝛾_𝑛=𝑔_𝑛,𝑛≥0, (36) where 𝑐_𝑛=∑^𝑛_𝑘=0𝑠_𝑘𝑠_𝑛−𝑘 and 𝛾_𝑛=∑^𝑛_𝑘=0(1/𝑘!)𝛽_𝑛−𝑘.

This solution, of course, is 𝑠(𝑡)=^∞_𝑛=0𝑠_𝑛𝑡^𝑛. (37)

3. Similarity Solution

We begin our approach by introducing the similarity variable 𝑥𝜉=√𝑡 (38) and look for solutions of (1) in the form 𝑢(𝑡,𝑥)=𝑡^𝑝𝜇(𝜉(𝑡,𝑥)), (39) where the number 𝑝 and the function 𝜇 are to be determined. Substituting (39) into (1) we find that 𝑡^𝑝−1𝑑𝑝𝜇−²𝜇𝑑𝜉²−𝜉2𝑑𝜇𝑑𝜉=0 (40) and so 𝑑²𝜇𝑑𝜉²+𝜉2𝑑𝜇𝑑𝜉=𝑝𝜇. (41) This is also difficult to solve for arbitrary values of 𝑝 but for special values we can do something. Define 𝐹(𝜉)=𝜇𝑒^𝜉4/4, and thus (41) can be transformed into 𝑑²𝐹𝑑𝜉²−𝜉2𝑑𝐹=𝑑𝜉2𝑝+12𝐹. (42) This has the general solution 𝐹(𝜉)=𝐶𝜉𝑒^𝜉2/4+𝐷−√𝜋𝜉𝑒^−𝜉2/4𝜉erf2−2, (43) where 𝐶 and 𝐷 are constants, provided 𝑝=1/2. Hence 𝜇(𝜉)=𝐶𝜉+𝐷−√𝜋𝜉𝑒^−𝜉2/2𝜉erf2−2𝑒^−𝜉2/4. (44) This gives a full solution for 𝑢(𝑡,𝑥) in the form 𝑢(𝑡,𝑥)=𝐶𝑥+𝐷−√𝜋𝑥𝑒^{−𝑥2/2𝑡}𝑥erf2√𝑡−2√𝑡𝑒^{−𝑥2/4𝑡}. (45) Taking the spatial derivative of the solution 𝑢 given by (45) 𝜕𝑢()√𝜕𝑥𝑡,𝑥=𝐶+𝐷𝜋𝑒^{−𝑥2/2𝑡}𝑥erf2√𝑡𝑥²𝑡+𝑥−1√𝑡𝑒^{−𝑥2/4𝑡}−𝑒^{−3𝑥2/4𝑡}. (46) From the Neumann boundary condition at 𝑥=0, we obtain 𝜕𝑢()𝜕𝑥𝑡,0=𝐶. (47) It follows that ℎ(𝑡)=𝑞 must be a constant. Hence 𝐶=−𝑞, (48) and from the Dirichlet condition we get −𝑞𝑠+𝐷−√𝜋𝑠𝑒^{−𝑠2/2𝑡}𝑠erf2√𝑡−2√𝑡𝑒^{−𝑠2/4𝑡}=𝑢^∗(𝑡). (49) Let 𝜆=𝑠/2√𝑡, and thus (49) becomes −2𝑞𝜆√𝑡+𝐷−2√𝜋𝜆𝑒^−2𝜆2erf(𝜆)−2𝑒^−𝜆2√𝑡=𝑢^∗(𝑡). (50) The analytical solution based on a similarity variable [1] requires that 𝑢^∗=0. This means that 𝑢^∗ does not depend on time. In this treatment we will assume that 𝑢^∗ is a function of 𝑡, that is, 𝑢^∗(𝑡)=𝛼√𝑡, where 𝛼 is a constant, which needs to be specified to realize an explicit solution. Thus −2𝑞𝜆+𝐷−2√𝜋𝜆𝑒^−2𝜆2erf(𝜆)−2𝑒^−𝜆2=𝛼. (51) Since 𝐷 in (51) is a constant, it follows that 𝜆 must also be constant. Thus 𝐷=𝛼+2𝑞𝜆−2√𝜋𝜆𝑒^−2𝜆2erf(𝜆)−2𝑒^−𝜆2,𝑠(𝑡)=2𝜆√𝑡. (52) The Stefan condition at the free boundary 𝑥=𝑠(𝑡) is 𝑠𝑑𝑠𝑑𝑡=−𝜕𝑢(),𝜕𝑥𝑡,𝑠 (53) where 𝑠𝑑𝑠𝑑𝑡=2𝜆²,𝜕𝑢()√𝜕𝑥𝑡,𝑠=−𝑞+𝐷𝜋𝑒^−2𝜆2erf(𝜆)4𝜆²𝑒−1+2𝜆^−𝜆2−𝑒^−3𝜆2. (54)(55) By putting (53) and (54) in the Stefan condition, we get the following transcendental equation: 2𝜆²+𝐷√𝜋𝑒^−2𝜆2erf(𝜆)4𝜆²−1+2𝜆𝑒^−𝜆2−𝑒^−3𝜆2=𝑞. (56) Thus we have proved the following theorem.

Theorem 2. For the Stefan problem with variable latent heat 𝜕𝑢=𝜕𝜕𝑡²𝑢𝜕𝑥²,0≤𝑥≤𝑠(𝑡),𝜕𝑢()𝑢((𝑡√𝜕𝑥𝑡,0=−𝑞,𝑡>0,𝑡,𝑠))=𝛼𝑠𝑡,𝑡>0,𝑑𝑠𝑑𝑡=−𝜕𝑢((𝑡𝜕𝑥𝑡,𝑠)),𝑡>0. (57) The solution is given by 𝑢()−√𝑡,𝑥=−𝑞𝑥+𝐷𝜋𝑥𝑒^{−𝑥2/2𝑡}𝑥erf2√𝑡−2√𝑡𝑒^{−𝑥2/4𝑡},𝑠(𝑡)=2𝜆√𝑡, (58)(59) where 𝐷=𝛼+2𝑞𝜆−2√𝜋𝜆𝑒^−2𝜆2erf(𝜆)−2𝑒^−𝜆2,𝑞=2𝜆²+𝐷√𝜋𝑒^−2𝜆2erf(𝜆)4𝜆²𝑒−1+2𝜆^−𝜆2−𝑒^−3𝜆2. (60)(61)

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.