# One-Phase Stefan-Like Problems with Latent Heat Depending on the Position and Velocity of the Free Boundary and with Neumann or Robin Boundary Conditions at the Fixed Face

Authors: Julieta Bollati; Domingo A. Tarzia
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A one-phase Stefan-type problem for a semi-infinite material which has as its main feature a variable latent heat that depends on the power of the position and the velocity of the moving boundary is studied. Exact solutions of similarity type are obtained for the cases when Neumann or Robin boundary conditions are imposed at the fixed face. Required relationships between data are presented in order that these problems become equivalent to the problem where a Dirichlet condition at the fixed face is considered. Moreover, in the case where a Robin condition is prescribed, the limit behaviour is studied when the heat transfer coefficient at the fixed face goes to infinity.


---

## Body

## 1. Introduction

Stefan-like problems have attracted growing attention in the last decades due to the fact that they arise in many significant areas of engineering, geoscience, and industry [1–9]. The classical Stefan problem describes the process of a material undergoing a phase change. Finding a solution to this problem consists in solving the heat-conduction equation in an unknown region which has also to be determined, imposing an initial condition, boundary conditions, and the Stefan condition at the moving interface. For an account of the theory we refer the reader to [10].In the classical Stefan problem the latent heat is assumed to be constant. In this paper, we are going to consider a variable one. This assumption is motivated by the fact that it becomes meaningful in the study of the shoreline movement in a sedimentary basis [11], in the one-dimensional consolidation with threshold gradient [12], in the artificial ground-freezing technique [13], and in nanoparticle melting [14], among others [15–21].Many papers deal with a latent heat that depends on the position of the free boundary (size-dependent latent heat). In [18], a Stefan problem with a latent heat given as a function of the position of the interface L=φ(s(t)) has been considered. This hypothesis corresponds to the practical case when the influence of phenomena such as surface tension, pressure gradients, and nonhomogeneity of materials are taken into account. In [11], the shoreline movement in a sedimentary basin was studied, from where a one-phase Stefan problem with a latent heat arises that increases linearly with distance from the origin; i.e., L=γs(t) (with γ a given constant). The generalization to the two-phase problem was done in [19]. Also, in [20] a latent heat defined as a power function of the position, i.e., L=γsn(t) (with γ a given constant and n an arbitrary nonnegative integer), was considered. The extension to a noninteger exponent was done in [21] for flux and temperature boundary conditions, while the two-phase case was presented in [13]. In [22], a convective (Robin) condition was imposed for the one-phase case while for the two-phase case the analysis was done in [23].In [12], a one-dimensional consolidation problem with a threshold gradient was studied. This problem can be reduced to a one-phase Stefan problem where the latent heat can be expressed as L=γ/s˙(t), that is to say, rate-dependent latent heat. It must be noticed that the case considered in [12] is not properly a Stefan problem because the velocity of the moving boundary disappears, and it has to be treated as a free boundary problem with implicit conditions [24, 25].Recently, in [26] it was defined a generalized one-phase Stefan-like problem for a semi-infinite material x>0 with a latent heat given by L=γsβ(t)s˙δ(t) (with γ a given constant and β and δ arbitrary real constants), i.e., latent heat depending on the position and velocity of the moving boundary, taking a Dirichlet boundary condition. This paper intends to complete this model, by considering two new boundary conditions (Neumann and Robin conditions) at the fixed face x=0.In Section2 we present a problem (P) with a variable latent heat and a generalized boundary condition at the fixed face. We will obtain its exact solution, following the methodology given in [12, 21, 26], obtaining as immediate consequence the similarity solutions to two different problems: one with Neumann condition at x=0 and the other with a Robin one. Special cases will be treated in order to recover solutions recently reported in literature. Moreover, in Section 3, the equivalence between these problems and the problem with a Dirichlet condition considered in [26] will be proved under certain relationships between data. For the problem with a Robin boundary condition at the fixed face, the limit behaviour when the heat transfer coefficient goes to infinity will be also analysed in Section 4. This analysis will allow us to show that the Robin condition constitutes a generalization of the Dirichlet one, as happens in classical heat transfer problems [27]. Also, in Section 5, we will provide some plots and table of values in order to track the position of the free front and to show how the latent heat changes in time.


## 2. Formulation of the Problems and Exact Solution

### 2.1. Statement of the Problems

In this paper, the exact solution of two different free boundary problems is obtained. They will be defined as particular cases of the following problem(P) that consists in finding the function u=u(x,t) and the moving boundary x=s(t) such that(1)∂u∂tx,t=a2∂2u∂x2x,t,0<x<st,(2)ust,t=0,t>0,(3)-k∂u∂xst,t=Lst,s˙ts˙t,t>0,(4)s0=0,(5)k∂u∂x0,t=h0tλu0,t-u∞tα/2,t>0,where α, λ, h0, a2 (diffusivity), and k (conductivity) are nonnegative constants.The problem defined by specifyingλ=0 in (P) will be referred to asproblem (PN). In this case, condition (5) corresponds to the Neumann boundary condition:(6)k∂u∂x0,t=-q0tα-1/2,t>0,q0>0where a time dependent heat flux characterized by q0=h0u∞>0 is applied at the fixed face x=0. This flux is proportional to the power α-1/2 of time.The problem defined by specifyingλ=1 in (P) will be referred to asproblem (PR). In this case, condition (5) corresponds to the Robin boundary condition:(7)k∂u∂x0,t=h0tu0,t-u∞tα/2,t>0,h0>0where u∞ characterizes the bulk temperature at a large distance from the fixed face x=0 and h0>0 characterizes the heat transfer at the fixed face.Comparing these problems with respect to the classical Stefan problem, the new feature to be observed is that condition (3) at the free interface can be thought of as a generalized Stefan condition where the latent heat term L(s(t),s˙(t)) is not constant but rather a function of position and velocity of the moving boundary. Furthermore, in order to obtain a similarity type solution for problems (P), (PN), and (PR), L will be specified as(8)Lst,s˙t=γsβts˙δt,with γ, β, and δ being nonnegative given constants.


### 2.2. Similarity Type Solutions

Before finding the similarity type solution to problem(P), the subsequent analysis will be necessary.Let us observe that if we use the similarity transformation presented in [20, 21](9)ux,t=tα/2φηwithη=x2at,then it is obtained that (10)∂2u∂x2x,t=tα/2-1φ″η14a2and∂u∂tx,t=tα/2-1α2φη-ηφ″η.Therefore, (1) is satisfied; i.e., ∂u/∂t(x,t)=a2∂2u/∂x2(x,t) if and only if(11)φ″η+2ηφ′η-2αφη=0.This second-order ordinary differential equation, known in literature as Kummer’s differential equation (see [28]), has a general solution that is given by(12)φη=tα/2C1M-α2,12,-η2+C2ηM-α2+12,32,-η2,with C1 and C2 being arbitrary constants. The function M(a,b,z) is called Kummer’s function or confluent hypergeometric function of the first kind and it is defined by the following series: (13)Ma,b,z=∑n=0∞anbnznn!,where b cannot be a nonpositive integer, and (a)n is the Pochhammer symbol defined by (14)a0=1,an=aa+1a+2⋯a+n-1.The detailed proof of the fact that the general solution of Kummer’s equation (11) can be written as (12) may be found in [26].The main properties of Kummer’s functionM(a,b,z) to be used throughout this paper can be found in [28] and they are stated in the following way:(15)Ma,b,0=1,(16)Ma,b,z=ezMb-a,b,-z,(17)e-z2=-2αz2M-α2+12,32,-z2M-α2+1,32,-z2+M-α2,12,-z2M-α2+12,12,-z2,(18)ddzMa,b,z=abMa+1,b+1,z,(19)ddzzb-1Ma,b,z=b-1zb-2Ma,b-1,z,(20)Ma,b,z≃ezza-bΓawhenz→∞,(21)M-n2,12,-z2=2n-1Γn2+1inerfc⁡z+inerfc⁡-z,n∈N,(22)zM-n2+12,32,-z2=2n-2Γn2+12inerfc⁡-z-inerfc⁡z,n∈N,where inerfc⁡· is the repeated integral of the complementary error function defined by(23)i0erfc⁡z=erfc⁡z=1-erf⁡z,erf⁡z=2π∫0ze-u2du,inerfc⁡z=∫z+∞in-1erfc⁡tdt.Now, we will look for the similarity solution to problem (P). In order to make the notation clearer, we will refer to the solution of problem (P) as the pair (u(x,t),s(t)) that satisfies (1)-(5).According to the previous analysis,u will satisfy (1) if it is written as(24)ux,t=tα/2C1M-α2,12,-η2+C2ηM-α2+12,32,-η2,with the similarity variable given by η=x/2at, where C1, C2 are constants to be determined so that u satisfies the rest of the conditions.Observe that from (2) it should be noted that φ defined by the transformation (9) has to satisfy φs(t)/2at=0 for all t>0. Therefore the moving boundary must adopt the following form:(25)st=2ξat,where ξ is a positive dimensionless coefficient to be determined.Hence, bearing in mind thatu is written as (24) and s as (25), finding the solution to problem (P) consists in determining the coefficients C1, C2, and ξ.The generalized boundary condition at the fixed face (5) and properties (15), (18)-(19) imply that(26)C2=2ah0kλC1-u∞.From condition (2), it can be deduced after some computations that(27)C1=u∞ξM-α/2+1/2,3/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.Therefore, replacing C1 in (26), we get(28)C2=-u∞M-α/2,1/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.Then, we have obtained C1 and C2 as functions of ξ.Finally, the Stefan-type condition given by (3) will give us an equation for ξ.Applying the derivation formulas (18)-(19) we claim that (29)∂u∂xst,t=tα-1/2u∞2a2αξ2M-α/2+1/2,3/2,-ξ2M-α/2+1,3/2,-ξ2-M-α/2,1/2,-ξ2M-α/2+1/2,1/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.Using relationships (16)-(17), the partial derivative of u is reduced to(30)∂u∂xst,t=-tα-1/2u∞2a1k/2ah0Mα/2+1/2,1/2,ξ2+λξMα/2+1,3/2,ξ2.Replacing (30) in (3) yields the following equality: (31)ku∞tα-1/22ak/2ah0Mα/2+1/2,1/2,ξ2+λξMα/2+1,3/2,ξ2=γ2βaβ+δ+1tβ-δ-1/2ξβ+δ+1,which makes sense if and only if α-1/2=β-δ-1/2, due to the fact that neither γ, ξ nor a depends on time. Thus, the similarity solution for problem (P) will exist if and only if(32)α=β-δ≥0and if ξ is a positive solution of the following equation:(33)ku∞γ2β+1aβ+δ+2fλz=zβ+δ+1,z>0,with(34)fλz=1k/2ah0Mα/2+1/2,1/2,z2+λzMα/2+1,3/2,z2.The notationfλ is adopted in order to emphasize the dependence of the solution to problem (P) on λ and therefore to obtain easily the solutions to problems (PN) and (PR).The task is now to prove the existence and uniqueness of solution to (33). From the relationships (15), (18), (19), and (20), we obtain that fλ satisfies(35)fλ′z=-k/2ah0α+1zMα/2+3/2,3/2,z2+λMα/2+1,1/2,z2k/2ah0Mα/2+1/2,1/2,z2+λzMα/2+1,3/2,z2<0,(36)fλ0+=2ah0k>0,(37)fλ+∞=0.We can deduce that the l.h.s. of (33) is a strictly decreasing function that goes from u∞h0/γ2βaβ+δ+1>0 to 0 when z increases from 0 to +∞, while the r.h.s. of (33), if β+δ+1>0, is a strictly increasing function that goes from 0 to +∞.In conclusion we obtain that ifβ+δ+1>0, we can ensure that (33) has a unique positive solution.It should be mentioned that due to restriction (32), i.e., α=β-δ≥0 and β+δ+1>0, we get that β≥max⁡δ,-1-δ.All the above analysis can be summarized in the following theorem.Theorem 1.
Letβ and δ be arbitrary real constants satisfying β≥max⁡(δ,-1-δ). Taking α=β-δ, there exists a unique solution (u,s) of a similarity type for problem (P), i.e., (1)-(5), which is given by (24) and (25), where C1 and C2 are given by formulas (27) and (28), respectively, and the dimensionless coefficient ξ is defined as the unique positive solution of (33).The solutions to problems(PN) and (PR) can be obtained as a consequence of Theorem 1, by fixing λ=0 or λ=1, respectively. As an immediate consequence we have the following results.Corollary 2 (caseλ=0).
Letβ and δ be arbitrary real constants satisfying β≥max⁡(δ,-1-δ). Taking α=β-δ, there exists a unique solution (uN,sN) of a similarity type forproblem (PN), i.e., (1)-(4), and (6) which is given by(38)uNx,t=tα/2C1NM-α2,12,-η2+C2NηM-α2+12,32,-η2,(39)sNt=2ξNat,where η=x/2at is the similarity variable. The coefficients C1N, C2N are defined by (40)C1N=2aq0ξNkM-α/2+1/2,3/2,-ξN2M-α/2,1/2,-ξN2,C2N=-2aq0k,and ξN is the unique positive solution of the equation ku∞/γ2β+1aβ+δ+2f0(z)=zβ+δ+1 that can be rewritten as(41)q0γ2βaβ+δ+1gz=zβ+δ+1,z>0,with(42)gz=1Mα/2+1/2,1/2,z2.Corollary 3 (caseλ=1).
Letβ and δ be arbitrary real constants satisfying β≥max⁡(δ,-1-δ). Taking α=β-δ, there exists a unique solution (uR,sR) of a similarity type forproblem (PR), i.e., (1)-(4) and (7), which is given by(43)uRx,t=tα/2C1RM-α2,12,-η2+C2RηM-α2+12,32,-η2,(44)sRt=2ξRat,where η=x/2at is the similarity variable. The coefficients C1R, C2R are defined by (45)C1R=u∞ξRM-α/2+1/2,3/2,-ξR2k/2ah0M-α/2,1/2,-ξR2+ξRM-α/2+1/2,3/2,-ξR2,C2R=-u∞M-α/2,1/2,-ξR2k/2ah0M-α/2,1/2,-ξR2+ξRM-α/2+1/2,3/2,-ξR2,and ξR is the unique positive solution of the following equation:(46)ku∞γ2β+1aβ+δ+2f1z=zβ+δ+1,z>0,with f1 defined by replacing λ=1 in fλ given in (34):(47)f1z=1k/2ah0Mα/2+1/2,1/2,z2+zMα/2+1,3/2,z2.Specifying different values forβ and δ in the above results, several solutions reported in literature can be recovered as a corollary. For instance, consider the following.Corollary 4.
The solution to the classical Stefan problem with a Neumann boundary condition at the fixed face can be recovered from Theorem1 by taking λ=0, β=δ=0.Takingβ=δ=0 and thus α=0, the latent heat L=γ is assumed to be constant like in the classical Stefan problem. In such case, fixing λ=0, the flux boundary condition is given by k∂u/∂x(0,t)=-q0/t. Moreover, the property M1/2,1/2,z=ez allows us to ensure that ξN is the unique solution to the following equation: (48)q0γae-z2=z,as was obtained in [29].Corollary 5.
The solutions provided in [20, 21] can be recovered from Theorem 1 by taking λ=0, β∈R+, and δ=0.Takingδ=0, we get that L=γsβ(t), i.e., a power function of the position. For such a case, by taking into account the fact that α=β we automatically obtain the solutions already presented in literature. It must be pointed out that if β is an integer, properties (21)-(22) should be applied.Remark 6.
It must be noticed that in case we want to recover a latent heat defined asL=γ/s˙(t), we have to set β=0, δ=-1, and thus α=1. However we cannot recover the solution given in [12], due to the fact that the boundary condition imposed at the fixed face (Dirichlet) does not agree with the boundary condition considered in problem (P).Corollary 7.
The solution to the classical Stefan problem with a Robin boundary condition at the fixed face can be recovered from Theorem1 by taking λ=1, β=δ=0 (See [27]).Corollary 8.
The solution to the Stefan problem studied in [22] can be recovered from Theorem 1 by taking λ=1, β∈R+, and δ=0.


## 2.1. Statement of the Problems

In this paper, the exact solution of two different free boundary problems is obtained. They will be defined as particular cases of the following problem(P) that consists in finding the function u=u(x,t) and the moving boundary x=s(t) such that(1)∂u∂tx,t=a2∂2u∂x2x,t,0<x<st,(2)ust,t=0,t>0,(3)-k∂u∂xst,t=Lst,s˙ts˙t,t>0,(4)s0=0,(5)k∂u∂x0,t=h0tλu0,t-u∞tα/2,t>0,where α, λ, h0, a2 (diffusivity), and k (conductivity) are nonnegative constants.The problem defined by specifyingλ=0 in (P) will be referred to asproblem (PN). In this case, condition (5) corresponds to the Neumann boundary condition:(6)k∂u∂x0,t=-q0tα-1/2,t>0,q0>0where a time dependent heat flux characterized by q0=h0u∞>0 is applied at the fixed face x=0. This flux is proportional to the power α-1/2 of time.The problem defined by specifyingλ=1 in (P) will be referred to asproblem (PR). In this case, condition (5) corresponds to the Robin boundary condition:(7)k∂u∂x0,t=h0tu0,t-u∞tα/2,t>0,h0>0where u∞ characterizes the bulk temperature at a large distance from the fixed face x=0 and h0>0 characterizes the heat transfer at the fixed face.Comparing these problems with respect to the classical Stefan problem, the new feature to be observed is that condition (3) at the free interface can be thought of as a generalized Stefan condition where the latent heat term L(s(t),s˙(t)) is not constant but rather a function of position and velocity of the moving boundary. Furthermore, in order to obtain a similarity type solution for problems (P), (PN), and (PR), L will be specified as(8)Lst,s˙t=γsβts˙δt,with γ, β, and δ being nonnegative given constants.


## 2.2. Similarity Type Solutions

Before finding the similarity type solution to problem(P), the subsequent analysis will be necessary.Let us observe that if we use the similarity transformation presented in [20, 21](9)ux,t=tα/2φηwithη=x2at,then it is obtained that (10)∂2u∂x2x,t=tα/2-1φ″η14a2and∂u∂tx,t=tα/2-1α2φη-ηφ″η.Therefore, (1) is satisfied; i.e., ∂u/∂t(x,t)=a2∂2u/∂x2(x,t) if and only if(11)φ″η+2ηφ′η-2αφη=0.This second-order ordinary differential equation, known in literature as Kummer’s differential equation (see [28]), has a general solution that is given by(12)φη=tα/2C1M-α2,12,-η2+C2ηM-α2+12,32,-η2,with C1 and C2 being arbitrary constants. The function M(a,b,z) is called Kummer’s function or confluent hypergeometric function of the first kind and it is defined by the following series: (13)Ma,b,z=∑n=0∞anbnznn!,where b cannot be a nonpositive integer, and (a)n is the Pochhammer symbol defined by (14)a0=1,an=aa+1a+2⋯a+n-1.The detailed proof of the fact that the general solution of Kummer’s equation (11) can be written as (12) may be found in [26].The main properties of Kummer’s functionM(a,b,z) to be used throughout this paper can be found in [28] and they are stated in the following way:(15)Ma,b,0=1,(16)Ma,b,z=ezMb-a,b,-z,(17)e-z2=-2αz2M-α2+12,32,-z2M-α2+1,32,-z2+M-α2,12,-z2M-α2+12,12,-z2,(18)ddzMa,b,z=abMa+1,b+1,z,(19)ddzzb-1Ma,b,z=b-1zb-2Ma,b-1,z,(20)Ma,b,z≃ezza-bΓawhenz→∞,(21)M-n2,12,-z2=2n-1Γn2+1inerfc⁡z+inerfc⁡-z,n∈N,(22)zM-n2+12,32,-z2=2n-2Γn2+12inerfc⁡-z-inerfc⁡z,n∈N,where inerfc⁡· is the repeated integral of the complementary error function defined by(23)i0erfc⁡z=erfc⁡z=1-erf⁡z,erf⁡z=2π∫0ze-u2du,inerfc⁡z=∫z+∞in-1erfc⁡tdt.Now, we will look for the similarity solution to problem (P). In order to make the notation clearer, we will refer to the solution of problem (P) as the pair (u(x,t),s(t)) that satisfies (1)-(5).According to the previous analysis,u will satisfy (1) if it is written as(24)ux,t=tα/2C1M-α2,12,-η2+C2ηM-α2+12,32,-η2,with the similarity variable given by η=x/2at, where C1, C2 are constants to be determined so that u satisfies the rest of the conditions.Observe that from (2) it should be noted that φ defined by the transformation (9) has to satisfy φs(t)/2at=0 for all t>0. Therefore the moving boundary must adopt the following form:(25)st=2ξat,where ξ is a positive dimensionless coefficient to be determined.Hence, bearing in mind thatu is written as (24) and s as (25), finding the solution to problem (P) consists in determining the coefficients C1, C2, and ξ.The generalized boundary condition at the fixed face (5) and properties (15), (18)-(19) imply that(26)C2=2ah0kλC1-u∞.From condition (2), it can be deduced after some computations that(27)C1=u∞ξM-α/2+1/2,3/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.Therefore, replacing C1 in (26), we get(28)C2=-u∞M-α/2,1/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.Then, we have obtained C1 and C2 as functions of ξ.Finally, the Stefan-type condition given by (3) will give us an equation for ξ.Applying the derivation formulas (18)-(19) we claim that (29)∂u∂xst,t=tα-1/2u∞2a2αξ2M-α/2+1/2,3/2,-ξ2M-α/2+1,3/2,-ξ2-M-α/2,1/2,-ξ2M-α/2+1/2,1/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.Using relationships (16)-(17), the partial derivative of u is reduced to(30)∂u∂xst,t=-tα-1/2u∞2a1k/2ah0Mα/2+1/2,1/2,ξ2+λξMα/2+1,3/2,ξ2.Replacing (30) in (3) yields the following equality: (31)ku∞tα-1/22ak/2ah0Mα/2+1/2,1/2,ξ2+λξMα/2+1,3/2,ξ2=γ2βaβ+δ+1tβ-δ-1/2ξβ+δ+1,which makes sense if and only if α-1/2=β-δ-1/2, due to the fact that neither γ, ξ nor a depends on time. Thus, the similarity solution for problem (P) will exist if and only if(32)α=β-δ≥0and if ξ is a positive solution of the following equation:(33)ku∞γ2β+1aβ+δ+2fλz=zβ+δ+1,z>0,with(34)fλz=1k/2ah0Mα/2+1/2,1/2,z2+λzMα/2+1,3/2,z2.The notationfλ is adopted in order to emphasize the dependence of the solution to problem (P) on λ and therefore to obtain easily the solutions to problems (PN) and (PR).The task is now to prove the existence and uniqueness of solution to (33). From the relationships (15), (18), (19), and (20), we obtain that fλ satisfies(35)fλ′z=-k/2ah0α+1zMα/2+3/2,3/2,z2+λMα/2+1,1/2,z2k/2ah0Mα/2+1/2,1/2,z2+λzMα/2+1,3/2,z2<0,(36)fλ0+=2ah0k>0,(37)fλ+∞=0.We can deduce that the l.h.s. of (33) is a strictly decreasing function that goes from u∞h0/γ2βaβ+δ+1>0 to 0 when z increases from 0 to +∞, while the r.h.s. of (33), if β+δ+1>0, is a strictly increasing function that goes from 0 to +∞.In conclusion we obtain that ifβ+δ+1>0, we can ensure that (33) has a unique positive solution.It should be mentioned that due to restriction (32), i.e., α=β-δ≥0 and β+δ+1>0, we get that β≥max⁡δ,-1-δ.All the above analysis can be summarized in the following theorem.Theorem 1.
Letβ and δ be arbitrary real constants satisfying β≥max⁡(δ,-1-δ). Taking α=β-δ, there exists a unique solution (u,s) of a similarity type for problem (P), i.e., (1)-(5), which is given by (24) and (25), where C1 and C2 are given by formulas (27) and (28), respectively, and the dimensionless coefficient ξ is defined as the unique positive solution of (33).The solutions to problems(PN) and (PR) can be obtained as a consequence of Theorem 1, by fixing λ=0 or λ=1, respectively. As an immediate consequence we have the following results.Corollary 2 (caseλ=0).
Letβ and δ be arbitrary real constants satisfying β≥max⁡(δ,-1-δ). Taking α=β-δ, there exists a unique solution (uN,sN) of a similarity type forproblem (PN), i.e., (1)-(4), and (6) which is given by(38)uNx,t=tα/2C1NM-α2,12,-η2+C2NηM-α2+12,32,-η2,(39)sNt=2ξNat,where η=x/2at is the similarity variable. The coefficients C1N, C2N are defined by (40)C1N=2aq0ξNkM-α/2+1/2,3/2,-ξN2M-α/2,1/2,-ξN2,C2N=-2aq0k,and ξN is the unique positive solution of the equation ku∞/γ2β+1aβ+δ+2f0(z)=zβ+δ+1 that can be rewritten as(41)q0γ2βaβ+δ+1gz=zβ+δ+1,z>0,with(42)gz=1Mα/2+1/2,1/2,z2.Corollary 3 (caseλ=1).
Letβ and δ be arbitrary real constants satisfying β≥max⁡(δ,-1-δ). Taking α=β-δ, there exists a unique solution (uR,sR) of a similarity type forproblem (PR), i.e., (1)-(4) and (7), which is given by(43)uRx,t=tα/2C1RM-α2,12,-η2+C2RηM-α2+12,32,-η2,(44)sRt=2ξRat,where η=x/2at is the similarity variable. The coefficients C1R, C2R are defined by (45)C1R=u∞ξRM-α/2+1/2,3/2,-ξR2k/2ah0M-α/2,1/2,-ξR2+ξRM-α/2+1/2,3/2,-ξR2,C2R=-u∞M-α/2,1/2,-ξR2k/2ah0M-α/2,1/2,-ξR2+ξRM-α/2+1/2,3/2,-ξR2,and ξR is the unique positive solution of the following equation:(46)ku∞γ2β+1aβ+δ+2f1z=zβ+δ+1,z>0,with f1 defined by replacing λ=1 in fλ given in (34):(47)f1z=1k/2ah0Mα/2+1/2,1/2,z2+zMα/2+1,3/2,z2.Specifying different values forβ and δ in the above results, several solutions reported in literature can be recovered as a corollary. For instance, consider the following.Corollary 4.
The solution to the classical Stefan problem with a Neumann boundary condition at the fixed face can be recovered from Theorem1 by taking λ=0, β=δ=0.Takingβ=δ=0 and thus α=0, the latent heat L=γ is assumed to be constant like in the classical Stefan problem. In such case, fixing λ=0, the flux boundary condition is given by k∂u/∂x(0,t)=-q0/t. Moreover, the property M1/2,1/2,z=ez allows us to ensure that ξN is the unique solution to the following equation: (48)q0γae-z2=z,as was obtained in [29].Corollary 5.
The solutions provided in [20, 21] can be recovered from Theorem 1 by taking λ=0, β∈R+, and δ=0.Takingδ=0, we get that L=γsβ(t), i.e., a power function of the position. For such a case, by taking into account the fact that α=β we automatically obtain the solutions already presented in literature. It must be pointed out that if β is an integer, properties (21)-(22) should be applied.Remark 6.
It must be noticed that in case we want to recover a latent heat defined asL=γ/s˙(t), we have to set β=0, δ=-1, and thus α=1. However we cannot recover the solution given in [12], due to the fact that the boundary condition imposed at the fixed face (Dirichlet) does not agree with the boundary condition considered in problem (P).Corollary 7.
The solution to the classical Stefan problem with a Robin boundary condition at the fixed face can be recovered from Theorem1 by taking λ=1, β=δ=0 (See [27]).Corollary 8.
The solution to the Stefan problem studied in [22] can be recovered from Theorem 1 by taking λ=1, β∈R+, and δ=0.


## 3. Equivalence to the Problem with Dirichlet Condition

In [26], the unique similarity solution of a problem defined by (1)-(4) with a Dirichlet boundary condition at the fixed face characterized by u0>0 was obtained; i.e.,(49)u0,t=tα/2u0>0,t>0.The problem defined by conditions (1)-(4) and (49) will be referred to asproblem (PD) and its solution will be referred to as the pair (uD,sD).According to [26], if β and δ are arbitrary real constants satisfying β≥max⁡(δ,-1-δ), taking α=β-δ, the unique solution to problem (PD) is given by(50)uDx,t=tα/2C1DM-α2,12,-η2+C2DηM-α2+12,32,-η2,(51)sDt=2ξDat,where (52)C1D=u0,C2D=-u0M-α/2,1/2,-ξD2ξDM-α/2+1/2,3/2,-ξD2,and ξD is the unique positive solution of the following equation:(53)ku0γaβ+δ+22β+1fz=zβ+δ+1,z>0,with(54)fz=1zMα/2+1,3/2,z2.In this section we will study conditions on the data of the problem (P) that guarantee its equivalence with the problem (PD). For equivalence it will be understood that both problems have the same solution.Consider the problem (P) with given data λ, u∞, h0 whose solution (u,s) is given by formulas (24) and (25) under the hypothesis that β and δ are arbitrary real constants with β≥max⁡(δ,-1-δ), and α=β-δ. Computing u(0,t) it is obtained that (55)u0,t=u∞ξM-α/2+1/2,3/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2tα/2,with ξ defined as the unique positive solution to (33). Suppose now that we fix (56)u0=u∞ξM-α/2+1/2,3/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2,and we solve the problem (PD) obtaining (uD,sD). Notice that the moving boundary sD is characterized by a dimensionless coefficient ξD that will be the unique solution to (53); i.e., (57)kγaβ+δ+22β+11zMα/2+1,3/2,z2u∞ξM-α/2+1/2,3/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2=zβ+δ+1.Notice that if we putz=ξ, the prior equation reduces to (33), meaning that z=ξ constitutes a solution to (57). Therefore, as the unique solution to (57) is given by ξD, ξD=ξ results. Then, it follows easily that C1D=C1, C2D=C2 obtaining as a consequence that the solution (uD,sD) with the u0 data given in function of λ,u∞,h0 coincides with the solution (u,s) of the problem (P).Conversely, consider the problem (PD) with a given data u0 whose solution (uD,sD) is given by formulas (50) and (51) under the assumption that β and δ are arbitrary real constants, β≥max⁡(δ,-1-δ), and α=β-δ. Computing ∂uD/∂x(0,t), the following is obtained: (58)∂uD∂x0,t=-tα-1/2u0M-α/2,1/2-ξD22aξDM-α/2+1/2,3/2,-ξD2.Let us consider (P) with the data h0 given by (59)h0=-ku0M-α/2,1/2,-ξD22aξDM-α/2+1/2,3/2,-ξD2λu0-u∞,fixing λ and u∞ such that λu0<u∞. The solution (u,s) of this problem can be obtained by (24) and (25). The free boundary s is characterized by a dimensionless coefficient ξ that is the unique solution of (33), i.e., satisfies (60)ku∞γ2β+1aβ+δ+21k/2ah0Mα/2+1/2,1/2,z2+λzMα/2+1,3/2,z2=zβ+δ+1,z>0.The prior equation has z=ξD as a solution due to the fact that if we replace z by ξD, it is obtained that (60) is equivalent to (53). As (60) has a unique solution given by ξ, we claim that ξ=ξD. In addition, by some computations, it becomes C1=C1D, C2=C2D and so the solution (u,s) to problem (P) given by a data h0 in function of u0 is equal to the solution (uD,sD) to problem (PD). Therefore the following theorem holds.Theorem 9.
Ifβ and δ are arbitrary real constants satisfying β≥max⁡(δ,-1-δ) and α=β-δ, then the problem (P) defined by condition (1)-(5) is equivalent to problem (PD) defined by (1)-(4) and (49), when the parameters λ,u∞, and h0 in the problem (P) are related to the parameter u0 in problem (PD) by the following expression:(61)u0=u∞ξM-α/2+1/2,3/2,-ξ2k/2ah0M-α/2,1/2,-ξ2+λξM-α/2+1/2,3/2,-ξ2.The coefficient ξ makes reference to the unique solution of (33) for problem (P) which will coincide with the unique solution of (53) for problem (PD).As a consequence of the above result, by fixingλ=0 and λ=1, respectively, we can obtain the following corollaries.Corollary 10 (caseλ=0).
Ifβ and δ are arbitrary real constants satisfying β≥max⁡(δ,-1-δ) and α=β-δ, then the problem (PN) defined by conditions (1)-(4) and (6) is equivalent to problem (PD) defined by (1)-(4) and (49), when the parameter q0 in the problem (PN) is related to the parameter u0 in problem (PD) by the following expression:(62)u0=2aq0kξNM-α/2+1/2,3/2,-ξN2M-α/2,1/2,-ξN2.The coefficient ξN makes reference to the unique solution of (41) for problem (PN) which will coincide with the unique solution of (53) for problem (PD).Corollary 11 (caseλ=1).
Ifβ and δ are arbitrary real constants satisfying β≥max⁡(δ,-1-δ) and α=β-δ, then the problem (PR) defined by conditions (1)-(4) and (7) is equivalent to problem (PD) defined by (1)-(4) and (49), when the parameters h0, u∞ in the problem (PR) are related to the parameter u0 in problem (PD) by the following expression:(63)u0=u∞ξRM-α/2+1/2,3/2,-ξR2k/2ah0M-α/2,1/2,-ξR2+ξM-α/2+1/2,-ξR2.The coefficient ξR makes reference to the unique solution of (46) for problem (PR) which will coincide with the unique solution of (53) for problem (PD).


## 4. Asymptotic Behaviour When the Coefficienth0→∞

In this subsection we are going to analyse the behaviour of the problem (PR) when the coefficient h0>0 which characterizes the heat transfer coefficient at the fixed face x=0 tends to infinity. Due to the fact that the solution of this problem depends on h0, we will rename it. Thus, we will consider uR(x,t,h0)≔uR(x,t) and sR(t)≔sR(t,h0) defined by (43)-(44).Let us define the problem (PD∞) defined by conditions (1)-(4) and the following condition of Dirichlet type at the fixed face x=0 given by(64)u0,t=tα/2u∞,where u∞ corresponds to the data of the problem (PR). Notice that the solution (uD∞,sD∞) to problem (PD∞) can be obtained from (50) and (51) replacing u0 by u∞.Then we are going to state the following result.Theorem 12.
Ifβ and δ are arbitrary real constants satisfying β≥max⁡(δ,-1-δ) and α=β-δ the problem (PR) converges to problem (PD∞) when h0 tends to infinity; i.e., (65)limh0→+∞⁡PR=PD∞.In this context the term “convergence” means that (66)limh0→+∞⁡ξRh0=ξD∞,limh0→+∞⁡sRt,h0=sD∞t,∀t>0,limh0→+∞⁡uRx,t,h0=uD∞x,t,∀t>0,x>0.Proof.
On the one hand, the free boundary solution to problem (PR) is characterized by a dimensionless parameter ξR(h0) that is the unique solution to (46); i.e., (67)f1z,h0=zβ+δ+1C∞,z>0,where C∞=ku∞/γaβ+δ+22β+1 and f1(z,h0)≔f1(z) given by (47). On the other hand, the moving boundary sD∞ is characterized by a dimensionless parameter ξD∞ which will be defined as the unique solution of (53) replacing u0 by u∞; i.e., (68)fz=zβ+δ+1C∞,z>0,where f is defined by (54).
We are going to prove that whenh0→∞, the coefficient ξR(h0) converges to the coefficient ξD∞. We know that zβ+δ+1/C∞ is a strictly increasing function that goes from 0 to +∞ when z increases from 0 to +∞; f is a strictly decreasing function that goes from +∞ to 0 and f1(z,h0) is a strictly decreasing function in z as well but decreases from 2h0a/k to 0 when z goes from 0 to +∞. After some computations, it can be seen that (69)fz-f1z,h0=k2ah0fzMα/2+1/2,1/2,z21/fz+k/2ah0Mα/2+1/2,1/2,z2>0,z>0.Therefore it can be concluded that 0<ξR(h0)<ξD∞, for all h0>0. In addition, when h0→∞ it can be easily seen that f1(z,h0)→f(z) and so ξR(h0)→ξD∞. Once this equality has been proved, by taking the limit in the definitions of C1R and C2R one can obtain the required convergence for uR and sR.


## 5. Computational Examples

In this section, we present and discuss some computational examples.From Theorem1, the solution to problem (P) is characterized by a dimensionless parameter ξ defined as the unique positive solution to (33). This equation can be rewritten as Fλ(z)=0 with(70)Fλz=ku∞γ2β+1aβ+δ+2fλz-zβ+δ+1=0,z>0.To solve the nonlinear equation Fλ(z)=0 we apply the following Newton iteration formula:(71)zi+1=zi-FλziFλ′zi,where zi is the value of z at the ith iteration step and(72)Fλ′z=-ku∞γ2β+1aβ+δ+2k/2ah0α+1Mα/2+3/2,3/2,z2+λMα/2+1,1/2,z2k/2ah0Mα/2+1/2,1/2,z2+λξMα/2+1,3/2,z2-β+δ+1zβ+δ.We have implemented a MATLAB program to compute the dimensionless coefficientξ for different values of the parameters. The stopping criterion used is the boundedness of the absolute error zi+1-zi<10-10.In addition, given that the latent heat behaves as a function of the free front, we will plotL in order to show how it changes in time. Observe that(73)L=γsβts˙δt=γ2ξatβξatδ=γ2βaβ+δξβ+δtβ-δ/2.Therefore it is deduced that the latent heat behaves as a power of time, i.e., L~tp with p<1 if β-δ<2, p=1 for β-δ=2 and p>1 in case β-δ>2. It must be pointed out that in all cases β should be β≥max⁡{δ,-1-δ} in order to meet the hypothesis of Theorem 1.Let us first analyse the problem with a Neumann boundary condition at the fixed face. From Corollary2, the solution of the problem (PN) is characterized by a dimensionless parameter ξN defined as the unique solution of (41). This equation can be rewritten as (70) specifying λ=0:(74)F0z=Q2βgz-zβ+δ+1=0,z>0,where g is given by (42) and the dimensionless parameter Q is defined by(75)Q=q0γaβ+δ+1.In Table 1 we present the computational results of ξN for different values of Q.Table 1
Computational results ofξN.
Q = 0.1 Q = 0.2 Q = 0.3 Q = 0.4 Q = 0.5 δ = 0, β = 0 0.0990 0.1927 0.2777 0.3531 0.5237 β = 1 0.2138 0.2912 0.3453 0.3875 0.4225  δ = - 1 / 2, β = 0 0.0100 0.0398 0.0879 0.1496 0.2172 β = 1 0.1319 0.2016 0.2543 0.2970 0.2952  δ = 1, β = 1 0.3534 0.4357 0.4904 0.5321 0.5661 β = 3 0.3838 0.4323 0.4627 0.4851 0.5031In Figure1, we plot the coefficient ξN that characterizes the free front sN, for different values of the parameters Q, β, and δ.Figure 1
Plot ofξN against Q for different values of β and δ.In Figure2, we can observe graphically what can be analytically deduced, in the sense that when δ=0, β=1, we obtain that the latent heat behaves as power of time, i.e., L~tp with p=1/2<1. In the case that δ=1, β=3, we obtain that p=1 and for δ=1, β=4, the power becomes p=3/2>1.Figure 2
Plot ofL against time for different values of β and δ assuming Q=0.5, a=1, and γ=1.Now, we turn to the problem with a Robin boundary condition at the fixed face. From Corollary3, the solution of the problem (PR) is characterized by a dimensionless parameter ξR defined as the unique solution of (46). This equation can be rewritten as (70) fixing λ=1:(76)F1z=Ste2β+1f1z-zβ+δ+1=0,z>0,where f1 is given by (47). Introducing the dimensionless parameter Ste, which constitutes a generalization of the Stefan number, and the generalized Biot number(77)Ste=u∞kγaβ+δ+2,Bi=h0ka,we get that f1 can be rewritten as f1(z)=1/1/2BiMα/2+1/2,1/2,z2+zMα/2+1,3/2,z2.In Table2 we present the computational results of ξR for different values of Bi, β, and δ, fixing Ste=0.5. Observe that the last column of the table intends to show that when Bi increases, ξR becomes closer to ξD which is the dimensionless parameter of the free front to the problem with Dirichlet condition. This convergence is in agreement with the prior section, taking into account the fact that analysing h0→∞ is equivalent to analysing Bi→∞.Table 2
Computational results ofξR.
Ste=0.5 Bi=1 Bi=10 Bi=50 Bi=100 ξ D δ = 0, β = 0 0.2926 0.4422 0.4601 0.4625 0.4648 β = 1 0.3490 0.4485 0.4617 0.4635 0.4652  δ = - 1 / 2, β = 0 0.1430 0.3375 0.3617 0.3648 0.3680 β = 1 0.2701 0.3837 0.3994 0.4015 0.4036  δ = 1, β = 1 0.4736 0.5514 0.5609 0.5621 0.5634 β = 3 0.4615 0.5181 0.5260 0.5270 0.5281In Figure3, we plot ξR against Bi for different values of δ and β, fixing Ste=0.5.Figure 3
Plot ofξR against Bi for different values of β and δ, with Ste=0.5.


## 6. Conclusions

In this paper two different one-phase Stefan-like problems were studied for a semi-infinite material. The main feature of both problems resides in the fact that a variable latent heat depending on the power of the position and the rate of change of the moving boundary is considered (L=γsβs˙δ). Using Kummer functions, exact solutions of similarity type were obtained for the cases when Neumann or Robin boundary conditions are imposed at the fixed face.In addition, the necessary and sufficient relationships between the data of the two problems in order to obtain an equivalence with the problem with a Dirichlet condition are obtained.For the problem with Robin boundary condition, the limit behaviour of the solution when the heat transfer coefficient at the fixed face goes to infinity was analysed, obtaining as a result the convergence to the solution of a Stefan problem with a Dirichlet boundary condition.This paper constitutes a mathematical generalization of the classical one because it can be obtained by fixing the parametersβ=1, δ=0. Also, the results obtained when a latent heat is considered as a linear or power function of the position of the free boundary can be recovered.We have provided tables and plots in order to show how the free front evolves in each case for specific values of the parameters.It is worth mentioning that finding exact solutions is meaningful not only to understand better the physical processes involved but also to verify the accuracy of numerical methods that solve Stefan problems.



---

*Source: 4960391-2018-08-05.xml*