One-phase Stefan-like problems with a latent heat depending on the position and velocity of the free boundary, and with Neumann or Robin boundary conditions at the fixed face

In this paper, 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.


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], 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 non-homogeneity of materials are taken into account. In [11], the shoreline movement in a sedimentary basin was studied, from where arises a one-phase Stefan problem with a latent heat 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 = γs n (t) (with γ a given constant and n an arbitrary non-negative integer) was considered. The extension to a non-integer exponent was done in [21] for a 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 a 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)ṡ δ (t) (with γ a given constant and β and δ arbitrary real constants), i.e., a 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 Section 2 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], and [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 are 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: where α, λ, h 0 , a 2 (diffusivity) and k (conductivity) are non-negative constants. The problem defined by specifying λ = 0 in (P ), will be referred to as problem (P N ). In this case, condition (5) corresponds to the Neumann boundary condition: where a time dependent heat flux characterized by q 0 = h 0 u ∞ > 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 as problem (P R ). In this case, condition (5) corresponds to the Robin boundary condition: where u ∞ characterizes the bulk temperature at a large distance from the fixed face x = 0 and h 0 > 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 the condition (3) at the free interface can be thought of as a generalized Stefan condition where the latent heat term L(s(t),ṡ(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 ), (P N ) and (P R ), L will be specified as: with γ, β and δ non-negative given constants.

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 following similarity transformation presented in [20] and [21]: then it is obtained that Therefore, equation (1) is satisfied, i.e. ∂u ∂t (x, t) = a 2 ∂ 2 u ∂x 2 (x, t) if and only if This second order ordinary differential equation, known in literature as Kummer's differential equation (see [28]), has a general solution that is given by with C 1 and C 2 arbitrarily 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 where b cannot be a non-positive integer, and (a) n is the Pochhammer symbol defined by (a) 0 = 1, (a) n = a(a + 1)(a + 2) · · · (a + n − 1).
The detailed proof of the fact that the general solution of the Kummer's equation (10) can be written as (11) may be found in [26].
The main properties of the Kummer's function M (a, b, z) to be used throughout this paper can be found in [28] and they are stated in the following way: where i n erf c(·) is the repeated integral of the complementary error function defined by 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 equation (1) if it is written as: with the similarity variable given by η = x 2a √ t and where C 1 , C 2 are constants to be determined so that u satisfies the rest of the conditions. Observe that from (2) it should be that ϕ defined by the transformation (9) has to satisfy ϕ s(t) 2a √ t = 0 for all t > 0. Therefore the moving boundary must adopt the following form: where ξ is a positive dimensionless coefficient to be determined. Hence, bearing in mind that u is written as (20) and s as (21), finding the solution to problem (P ) consists in determining the coefficients C 1 , C 2 and ξ.
The generalized boundary condition at the fixed face (5), and properties (12), (15)- (16) imply that: From condition (2), it can be deduced after some computations that: Therefore, replacing C 1 in (22), we get Then, we have obtained C 1 and C 2 as functions of ξ. Finally, the Stefan-type condition given by (3) will give us an equation for ξ.
Applying the derivation formulas (15)-(16) we claim that: Using the relationships (13)- (14), the partial derivative of u is reduced to Replacing (25) in (3) yields to the following equality: , due to the fact that nor γ, ξ, or a depends on time. Thus, the similarity solution for problem (P ) will exist if and only if and if ξ is a positive solution of the following equation: with The notation f λ is adopted in order to emphasize the dependence of the solution to problem (P ) on λ, and therefore to obtain easily the solutions to problems (P N ) and (P R ).
The task is now to prove the existence and uniqueness of solution to equation (27). From the relationships (12), (15), (16) and (17), we obtain that f λ satisfies: We can deduce that the l.h.s. of equation (27) is a strictly decreasing function that goes from u∞h 0 γ2 β a β+δ+1 > 0 to 0 when z increases from 0 to +∞, while the r.h.s. of equation (27), if β + δ + 1 > 0, is a strictly increasing function that goes from 0 to +∞.
As a conclusion we obtain that if β + δ + 1 > 0, we can ensure that (27) has a unique positive solution.
All the above analysis can be summarized in the following theorem: Theorem 2.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. equations (1)-(5), which is given by (20) and (21), where C 1 and C 2 are given by the formulas (23) and (24) respectively, and the dimensionless coefficient ξ is defined as the unique positive solution of the equation (27).
The solutions to problems (P N ) and (P R ) can be obtained as a consequence of Theorem 2.1, by fixing λ = 0 or λ = 1, respectively. As an immediate consequence we have the following results.
Corollary 2.1. (case λ = 0) Let β and δ be arbitrary real constants satisfying β ≥ max(δ, −1 − δ). Taking α = β −δ, there exists a unique solution (u N , s N ) of a similarity type for problem (P N ), i.e. equations (1)-(4), and (6) which is given by and ξ N is the unique positive solution of the equation ku∞ γ2 β+1 a β+δ+2 f 0 (z) = z β+δ+1 , that can be rewritten as with Corollary 2.2. (case λ = 1) Let β and δ be arbitrary real constants satisfying β ≥ max(δ, −1 − δ). Taking α = β − δ, there exists a unique solution (u R , s R ) of a similarity type for problem (P R ), i.e. equations (1)-(4) and (7), which is given by where η = x 2a √ t is the similarity variable. The coefficients C 1R , C 2R are defined by and ξ R is the unique positive solution of the following equation with f 1 defined by replacing λ = 1 in f λ given in (28): Specifying different values for β and δ in the above results, several solutions reported in literature can be recovered as a corollary. For instance: 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) = − q 0 √ t . Moreover, the property M 1 2 , 1 2 , z = e z , allows us to ensure that ξ N is the unique solution to the following equation: as was obtained in [29].
Taking δ = 0, we get that L = γs β (t), i.e. a power function of the position. For such a case, by taking into account that α = β we automatically obtain the solutions already presented in literature. It must be pointed out that if β is an integer, properties (18)- (19) should be applied.
Remark 2.1. It must be noticed that in case we want to recover a latent heat defined as L = γ s(t) , we have to set β = 0, δ = −1, and thus α = 1. However we can not 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 2.5. The solution to the classical Stefan problem with a Robin boundary condition at the fixed face can be recovered from Theorem 2.1 by taking λ = 1, β = δ = 0 (See [27]).

Equivalence to the problem with Dirichlet condition
In [26], the unique similarity solution of a problem defined by equations (1)-(4) with a Dirichlet boundary condition at the fixed face characterized by u 0 > 0 was obtained, i.e.: The problem defined by conditions (1)-(4) and (40) will be referred to as problem (P D ) and its solution will be referred to as the pair (u D , s D ). According to [26], if β and δ are arbitrary real constants satisfying β ≥ max(δ, −1 − δ), taking α = β − δ, the unique solution to problem (P D ) is given by: where and ξ D is the unique positive solution of the following equation with In this section we will study conditions on the data of the problem (P ) that guarantee its equivalence with the problem (P D ). For equivalence it will be understood that both problems have the same solution.
Consider the problem (P ) with a given data λ, u ∞ , h 0 whose solution (u, s) is given by formulas (20), and (21) under the hypothesis that β and δ are arbitrary real constants with β ≥ max(δ, −1 − δ), and α = β − δ. Computing u(0, t) it is obtained that with ξ defined as the unique positive solution to equation (27). Suppose now that we fix and we solve the problem (P D ) obtaining (u D , s D ). Notice that the moving boundary s D is characterized by a dimensionless coefficient ξ D that will be the unique solution to equation (43), i.e. k γa β+δ+2 2 β+1 Notice that if we put z = ξ, the prior equation reduces to equation (27), meaning that z = ξ constitutes a solution to (45). Therefore, as the unique solution to (45) is given by ξ D , it results that ξ D = ξ. Then, it follows easily that C 1D = C 1 , C 2D = C 2 obtaining as consequence that the solution (u D , s D ) with the u 0 data given in function of λ, u ∞ , h 0 , coincides with the solution (u, s) of the problem (P ).
Conversely, consider the problem (P D ) with a given data u 0 whose solution (u D , s D ) is given by formulas (41) and (42) under the assumption that β and δ are arbitrary real constants, β ≥ max(δ, −1 − δ), and α = β − δ. Computing ∂u D ∂x (0, t) it is obtained that: Let us consider (P ) with the data h 0 given by , fixing λ and u ∞ such that λu 0 < u ∞ . The solution (u, s) of this problem can be obtained by (20) and (21). The free boundary s is characterized by a dimensionless coefficient ξ that is the unique solution of equation (27), i.e. satisfies: The prior equation has z = ξ D as a solution due to the fact that if we replace z by ξ D , it is obtained that equation (46) is equivalent to equation (43). As (46) has a unique solution given by ξ, we claim that ξ = ξ D . In addition, by some computations, it becomes C 1 = C 1D , C 2 = C 2D and so the solution (u, s) to problem (P ) given by a data h 0 in function of u 0 is equal to the solution (u D , s D ) to problem (P D ). Therefore the following theorem holds The coefficient ξ makes reference to the unique solution of equation (27) for problem (P ) which will coincide with the unique solution of (43) for problem (P D ).
As a consequence of the above result, by fixing λ = 0 and λ = 1, respectively we can obtain the following corollaries: Corollary 3.1. (case λ = 0) If β and δ are arbitrary real constants satisfying β ≥ max(δ, −1− δ) and α = β − δ, then the problem (P N ) defined by condition (1)-(4) and (6) is equivalent to problem (P D ) defined by (1)-(4) and (40), when the parameter q 0 in the problem (P N ) is related with the parameter u 0 in problem (P D ) by the following expression: The coefficient ξ N makes reference to the unique solution of equation (34) for problem (P N ) which will coincide with the unique solution of (43) for problem (P D ).
Corollary 3.2. (case λ = 1) If β and δ are arbitrary real constants satisfying β ≥ max(δ, −1− δ) and α = β − δ, then the problem (P R ) defined by condition (1)-(4) and (7) is equivalent to problem (P D ) defined by (1)-(4) and (40), when the parameters h 0 , u ∞ in the problem (P R ) are related with the parameter u 0 in problem (P D ) by the following expression: The coefficient ξ R makes reference to the unique solution of equation (38) for problem (P R ) which will coincide with the unique solution of (43) for problem (P D ).
4 Asymptotic behaviour when the coefficient h 0 → ∞ In this subsection we are going to analyse the behaviour of the problem (P R ) when the coefficient h 0 > 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 h 0 , we will rename it. Thus, we will consider u R (x, t, h 0 ) := u R (x, t) and s R (t) := s R (t, h 0 ) defined by equations (36)-(37). Let us define the problem (P D∞ ) defined by conditions (1)-(4) and the following condition of Dirichlet type at the fixed face x = 0 given by where u ∞ corresponds to the data of the problem (P R ). Notice that the solution (u D∞ , s D∞ ) to problem (P D∞ ) can be obtained from (41) and (42) replacing u 0 by u ∞ . Then we are going to state the following result: Theorem 4.1. If β and δ are arbitrary real constants satisfying β ≥ max(δ, −1 − δ), and α = β − δ the problem (P R ) converges to problem (P D∞ ) when h 0 tends to infinity, i.e.: In this context the term "convergence" means that: Proof. On the one hand, the free boundary solution to problem (P R ) is characterized by a dimensionless parameter ξ R (h 0 ) that is the unique solution to equation (38), i.e.
where C ∞ = ku ∞ γa β+δ+2 2 β+1 and f 1 (z, h 0 ) := f 1 (z) given by (39). On the other hand, the moving boundary s D∞ is characterized by a dimensionless parameter ξ D∞ which will be defined as the unique solution of the equation (43) replacing u 0 by u ∞ , i.e.
where f is defined by (44).
We are going to prove that when h 0 → ∞, the coefficient ξ R (h 0 ) 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 f 1 (z, h 0 ) is a strictly decreasing function in z as well but decreases from 2h 0 a k to 0 when z goes from 0 to +∞. After some computations, it can be seen that: Therefore it can be concluded that 0 < ξ R (h 0 ) < ξ D∞ , for all h 0 > 0. In addition, when h 0 → ∞ it can be easily seen that f 1 (z, h 0 ) → f (z) and so ξ R (h 0 ) → ξ D∞ . Once this equality has been proved, by taking the limit in the definitions of C 1R and C 2R one can obtain the required convergence for u R and s R .

Computational examples
In this section, we present and discuss some computational examples. From Theorem 2.1, the solution to problem (P ) is characterized by a dimensionless parameter ξ defined as the unique positive solution to equation (27). This equation can be rewritten as To solve the nonlinear equation F λ (z) = 0 we apply the following Newton's iteration formula: where z i is the value of z at the ith iteration step and 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 |z i+1 − z i | < 10 −10 .
In addition, given that the latent heat behaves as a function of the free front, we will plot L in order to show how it changes in time. Observe that Therefore it is deduced that the latent heat behaves as a power of time i.e L ∼ t p 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 2.1. Let us first analyse the problem with a Neumann boundary condition at the fixed face. From Corollary 2.1, the solution of the problem (P N ) is characterized by a dimensionless parameter ξ N defined as the unique solution of equation (34). This equation can be rewritten as (53) specifying λ = 0: where g is given by (35) and the dimensionless parameter Q is defined by: In Table 1 we present the computational results of ξ N for different values of Q. In Figure 1, we plot the coefficient ξ N that characterizes the free front s N , for different values of the parameters Q, β, δ.
In Figure 2, 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 ∼ t p 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.  Now, we turn to the problem with a Robin boundary condition at the fixed face. From Corollary 2.2, the solution of the problem (P R ) is characterized by a dimensionless parameter ξ R defined as the unique solution of equation (38). This equation can be rewritten as (53) fixing λ = 1: where f 1 is given by (39). Introducing the dimensionless parameter Ste, which constitutes a generalization of the Stefan number, and the generalized Biot number we get that f 1 can be rewritten f 1 (z) = In Table 2 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 that analysing h 0 → ∞ is equivalent to analysing Bi → ∞.  In Figure 3, we plot ξ R against Bi for different values of δ and β, fixing Ste=0.5.

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 βṡδ ). Using Kummer functions, exact solutions of similarity type were obtained for the cases when a 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 to mention 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.