Solution to the 1D Stefan problem using the unified transform method

In this paper the one-dimensional two-phase Stefan problem is studied analytically leading to a system of non-linear Volterra-integral-equations describing the heat distribution in each phase. For this the unified transform method has been employed which provides a method via a global relation, by which these problems can be solved using integral representations. To do this, the underlying partial differential equation is rewritten into a certain divergence form, which enables to treat the boundary values as part of the integrals. Classical analytical methods fail in the case of the Stefan problem due to the moving interface. From the resulting non-linear integro-differential equations the one for the position of the phase change can be solved in a first step. This is done numerically using a fix-point iteration and spline interpolation. Once obtained, the temperature distribution in both phases is generated from their integral representation.


Introduction
Many important heat conduction problems occurring in science, nature and engineering involve a phase change due to melting or freezing. Mathematically, these are special cases of moving phase boundary problems, in which the evolution of the interface itself is unknown and hence * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
where the author applied the approach on the heat conduction on a ring with two phases in order to provide an explicit solution for the temperature in both phases, where the interface position was assumed to be known. Another application of the UTM on the heat equation can be seen in [20], where the authors obtained an explicit solution for each layer of the onedimensional multi-layer heat equation, whereas once again the interface position was assumed to be known. In [21] the authors solved the heat conduction on networks of multiply connected rods using the UTM, where both the connectivity of the rods and their size are known, while the temperature and the heat fluxes represent the searched quantities. The variety of application possibilities for the UTM can be seen amongst others in [22], where the authors were able to simulate brain tumor growth within heterogeneous environments. The UTM rewrites the given partial differential operation (PDE) into a divergence form and applies Gauss' theorem and the Fourier transformation to obtain the general solution.
The objective of this work is twofold. First, we present analytical expressions for the temperature distributions and the interface position for the one-dimensional two-phase case. Second, we lay out a numerical scheme to solve the resulting integro-differential equations in a similar way as done in [23], where the authors used the UTM, to solve the two-phase Couette flow with transpiration through both walls. The solution is thus obtained in the form of integrals depending on initial and boundary values. The unknown values at the moving interface were determined by a system of linear Volterra-integral-equations, where these equations were solved by using a standard marching method (see e.g. [24]).
In the following this work proceeds by firstly describing the physical problem in section 2 including the underlying PDE system as well as the boundary and coupling condition in the Stefan problem. In section 3, the UTM is applied on the Stefan problem leading to a nonlinear integro-differentiated system in time for the temperature distributions in both phases and the interface. In section 4 the obtained integro-differential equations are evaluated numerically, thus leading to explicit solutions for the temperature distributions and the interface motion, which are presented in section 5. Finally, in the last section 6 the obtained analytical expressions, the numerical results and open problems are discussed.

Problem statement
The one-dimensional Stefan problem for two phases describes the temperature distribution in a fluid undergoing a phase change from its solid to its liquid phase due to an imposed temperature gradient. For this rod of the liquid phase on the interval [0, g(t)] is considered, which is separated by the moving phase boundary g(t) from the solid phase occupying the interval [g(t), L], as shown in figure 1. The boundary at x = 0 of the liquid phase and the boundary at x = L of the solid phase are non-moving. Due to the temperature gradients f 1 (t) on the left boundary of the liquid phase and f 2 (t) on the right boundary of the solid phase, the solid phase begins to melt, resulting in a moving phase boundary between the two phases. The heat equation in each phase for the temperature distributions T 1 (x, t) and T 2 (x, t) read as where T j t (x, t) stands for the time derivative of the temperature distribution for each phase, while T j x (x, t) and T j x,x (x, t) stand for the first order or second order spatial derivative, with j = 1, 2. The thermal diffusivity is given as α j = λ j ρ j c p, j with the thermal conductivity λ j , the density ρ j and the specific heat capacity c p, j of each phase. The initial temperature distributions and the phase boundary are given as Moreover the introduced temperature gradients f 1 (t), f 2 (t) are given as The melting condition on the moving phase boundary leads to two further conditions describing the heat distributions in each phase, reading as The Stefan-condition, which is based on the formulation in [15], finally leads to the last needed condition and reads as In the above system, the initial conditions (2a)-(2c) and the boundary conditions (3a) and (3b) are given, while the temperature distribution T( and the melting point position g(t) are the searched quantities.

Unified transform method
Using the UTM, as described in [25,26] equations for both the heat distributions for each phase and an equation for the interface condition can be derived, similar as done in [15]. Due to the extension to the full two-phase case as compared to the one-phase study done in [15] the subsequent numerical evaluation however differs significantly. Additionally due to the twophase nature, the spatial derivatives of the temperature at the interface T 1 y (g(s), s) and T 2 y (g(s), s) cannot be simply replaced by the temporal derivative of the interface as one in [15] for the onephase case. Thus, in this work the Stefan condition (5) is applied after general expressions for the temperature distributions in each phase T 1 (x, t) and T 2 (x, t) are obtained. For this purpose each phase is considered as a single phase individually.
The first step of the UTM is to rewrite the given PDE into a specific divergence form. The temperature distribution for each phase is multiplied by e −ikx+w j (k)t as in [26], with k ∈ C. Thus both heat equations (1a) and (1b) will be rewritten to the following conservative form: where w i (k), j = 1, 2, denotes the dispersion relation and is given by Applying Gauss's theorem on the divergence form of the heat equation on the considered for each phase as shown in figure  1, leads to a global relation for each domain separately General integral representations of the heat distributions T 1 (x, t) and T 2 (x, t) can be obtained by first multiplying (8) with e −w 1 (k)t and (9) with e −w 2 (k)t and then applying the inverse Fourier transformation, leading to

Eliminating the unknown Dirichlet conditions
The general integral representations (10) and (11) for each heat distribution still contain unknown Dirichlet boundary conditions (DBCs) T 1 (0, s) and T 2 (L, s). However, using the symmetry for the dispersion relation, i.e. w j (k) = w j (−k), the DBC can be eliminated. Firstly, the Dirichlet condition on the left boundary for the first phase is eliminated, by replacing k with −k in equation (8), leading to (12) Equation (12) is multiplied with e ikx−w 1 (k)t , resulting in Equation (13) is inserted into (10), leading to The unknown Dirichlet condition for the right boundary of the second phase is eliminated on a similar way, by once again replacing k with −k in equation (9), resulting in Equation (16) is inserted into equation (11), which leads to After changing the order of the integration in (14) and (17), i.e. integrating with respect to k first, the expressions can now be evaluated. The results for the temperature distributions are given as The terms for each integral representation in (14) and (17) containing T 1 (y, t) and T 2 (y, t) vanish due to Jordan's lemma (for details see e.g. [25,27]). Thus the temperature distributions for each phase could be obtained, using only known boundary and initial values and the temperature distributions.
In order to implement the Stefan condition (5) the heat fluxes T 1 y (g(t)) and T 2 y (g(t)), which are the spatial derivatives of each single phase (18) and (19), are evaluated at the moving interface. To evaluate these spatial derivatives for each phase at the interface, it is necessary to apply a lemma introduced by Friedman [28] while differentiating (18) and (19), where the lemma reads as with the fundamental solution of the heat equation being employed as the kernel K(x, t, g(s), s) as This lemma is being applied to the problem presented in this work, where presently the kernel is defined as This results in with i = 1, 2 representing the kernel for each phase. The sign of the additional term 1 2 √ α i ρ(t) depends on the direction in which the interface is being approached. The complete derivation of (23) can be found in the appendix. Thus by applying (23) the derivative of the temperature distribution at the interface of the first phase is given as which can be simplified to In an analogous way the derivative of the temperature distribution at the interface of the second phase can be obtained by once again applying (23) as which then can be simplified to Finally, the heat fluxes (25) and (27) can be used in order to use the Stefan condition (5) as (5) = α 2 (27) − α 1 (25), leading to a non-linear integro-differential equation for the evolution of the interface position g(t), which reads as 4α 2 (t−s) (g(t) − g(s))ds ds .
The above non-linear integro-differential equation for g(t) is solved numerically in the subsequent section.

Numerical methods
In order to solve the set of integro-differential equations for the temperature distributions ( (18) and (19)) and the interface evolution (28), a similar approach to the one used by Guenther [29] on a one-phase Stefan problem is employed. The author approximated the derivative of the interface g t (t) as a spline using polynomial expressions.
Since the unknown interface occurs as g(t) and its derivative g t (t) in the expressions of interest (18), (19) and (28) the first step is to express g(t) through its derivative and its initial position as g n (t) = g 0 + t n 0 g t (s)ds. Using this expression as well as the estimated heat fluxes T 1 x,n (g n (t)) and T 2 x,n (g n (t)) at the phase boundaries, the phase motion of the phase boundary g n (t) and the heat fluxes T 1 x,n (g n (t)) and T 2 x,n (g n (t)) for any t ∈ [0, t end ] can be obtained, whereas t end represents the end time. For each time step n the time is given within the interval t = [t 0 , t 1 , . . . , t n ] with n = 0, 1, . . . , N. Thus the size of the time step δt is defined as δt = t n+1 − t n . Therefore for any time, a fix-point iteration of iteration number l for the motion of the phase boundary can be used as A(g t,l , T 1 x,l (g l (t)), T 2 x,l (g l (t))) → (g t,l+1 , T 1 x,l+1 (g l+1 (t)), T 2 x,l+1 (g l+1 (t))).
The operator A represents the non-linear system of equations in (18), (19) and (28). Repeated execution of iteration leads to the next time step. The initial conditions are used for the first time step. The derivative of both the phase boundary g t,n and the heat fluxes T 1 x,n (g(t)) and T 2 x,n (g(t)) are constructed as spline approximations of the previous time steps and current estimates. Through multiple iteration l, the motion of the interface and the heat fluxes is calculated based on the previous time step. Numeric integrations of the spline structures lead to the phase boundary position and the temperature fields. By default the order for numeric integration and spline approximation is set to 3, resulting in a cubic approximation. Now all the expressions for a new estimate of the current time step can be calculated, since all the terms on the left-hand side in (29) are known. An abort criterion value is introduced as g l+1 n (t) − g l n (t) < δ . If the interface motion for one iteration step deviates from the previous estimate more than δ the iteration shall be started over using the new approximation as a base. For the next time step, the base assumptions are now the latest Table 1. Parameters for mimicking the one-phase Stefan problem.
Parameter is obtained. The latter numeric scheme was coded in MATLAB.

Results
The test case of a one-phase Stefan problem developed in [9] was used to prove the accuracy of the present solution scheme in section 5.1 by comparing with the results provided by the author. Furthermore, results for the variation of initial conditions from the symmetric case, in which both phases are identical, are presented in section 5.3, showing their influence on the motion of the phase boundary.

Verification studies
The one-phase Stefan problem has been studied already in various ways [9]. For comparison with the present scheme we consider the work by Asaithambi [9], who studied the evolution of the phase boundary for different classical configurations. Since the one-phase case, he studied, represents a reduced case of the one presented in this work, our scheme can be verified by comparing our results to the results in Asaithambi. The parameters used for the numeric scheme are chosen in accordance to the case presented by Asaithambi and are displayed in table 1.
Our results are presented alongside with those of Asaithambi in table 2, where the author is using the weak or Galerkin formulation of the initial boundary value problem in order to reduce it to a system of initial-value problems in ordinary differential equations. Various finitedifference marching methods were applied to solve the resulting initial-value problems. The results of four of these methods are presented in table 2. The first and second explicit method in [9] differ in the formulation of the finite-element approximation of the solution for the temperature, while the first and second implicit method differ in the size of the time steps with the second method having smaller time steps. The results of our numerical scheme match the values of Asaithambi within 0.5 percent.
A second test case of validating our numerical scheme is by computing a symmetric case with no interface evolution. For mirrored initial conditions in both phases and the same material parameters α 1 = α 2 , no movement of the phase boundary occurs. In figures 3 and 4 this case is marked by a solid line and from here on named the symmetric case. The parameters used for this numeric scheme are displayed in table 3, whereas t end represents the end time. Table 3. Parameters for the symmetric two-phase Stefan problem. Parameter 0.0001 0.05 0.000 01 Table 4. Physical and numerical parameters for the symmetric two-phase Stefan problem. Parameter

Temporal order of numerical convergence rate and numerical parameters
The numerical order of convergence and with it the choice of the time step δt has an impact on the accuracy of the obtained solution. Therefore, the influence is studied in the following. Table 4 shows the set of parameters used, within the present subsection. The numerical order of convergence and hence the time step influence is studied by choosing the time step as follows,

Parameter studies
In the following the influence of the initial temperature distribution T 1 0 (x) and the thermal diffusivity α 1 on the motion of the phase boundary and the temperature distribution is being studied. The same settings as for the symmetric two-phase Stefan problem in table 3 is used as a reference set-up, whereas each parameter is modified separately. The corresponding figures display the evolution of the phase boundary g(t) over the time.
First the influence of the initial temperature distribution is investigated. The initial temperature distribution for the first phase T 1 0 (x) is being increased and decreased in comparison to the reference case in order to visualize the effect on the moving interface g(t). In the balanced case T 1 0 (x) = 2x − 1 is used, while the increased and decreased cases consider Figure 3 shows the interface motion for the five initial temperature distributions. For an increased value of T 1 0 (x), compared to the symmetric case, the interface g(t) moves more towards the boundary of the second phase, whereas for lower temperature distributions T 1 0 (x) the interface g(t) moves more towards the boundary of the first phase, which is due to the fact that the increased initial temperature distribution leads to the melting of the solid phase resulting in the expansion of the liquid phase. For a decreased initial value the opposite is the case, since now the magnitude of the temperature of the solid phase is higher than for the liquid phase. It is also visible, that the interface moves faster for the increased initial distributions than for the cases using decreased initial distributions.
Next, the influence of α 1 is analysed by investigating the moving interface for an increased and a decreased value of α 1 in comparison to the reference case and results are presented in figure 4. An increased value of the heat conductivity α 1 results in the expansion of the first phase, which can be seen as the interface also moves more into the second phase towards its boundary. A faster motion of the interface g(t) can be observed. Decreasing the heat conductivity leads to the opposite observation. For every decreased heat conductivity, the interface reaches more into the first phase, leading to a retraction of this phase. This observation can be explained due to the fact that in the case in which α 1 > α 2 the first phase can transfer its temperature faster than the second phase, which then results in the melting of the ice and thus the expansion of the first phase, while in the example of α 1 < α 2 the opposite is the case.

Reconstruction of the temperature distributions
The temperature distributions in (18) and (19) can be computed for every time, once a solution for the motion of the phase boundary g(t) has been obtained. The set of temperature distributions is displayed over the domain for multiple time steps in figure 5 for the symmetric case as described in table 3. The time steps used are t = 0.002, 0.018, 0.034, 0.05. One can see that the temperature distribution in both phases stay fully symmetric over the time. Furthermore, with increasing time, one can observe that the temperatures of both phases equalize and simultaneously move towards T 1 (x) = T 2 (x) = 0.  The influence of the initial temperature distribution T 1 0 (x) is investigated. Figure 6 shows the temperature distribution for T 1 0 (x) = 0.75(2x − 1). It can be seen that a decreased initial temperature distribution as compared to the symmetric case T 1 0 (x) = (2x − 1) results in a temperature distribution for the first phase which decays faster towards T 1 (x) = 0, which is also to be expected, since the initial temperature distribution is already closer to T 1 0 (x) = 0. Furthermore, a kink in the temperature curve can be seen at the location of the phase boundary, since now the temperature distribution in both phases differ from each other. It can also be seen  that the temperature at the solid boundary of the second phase did not change compared to the symmetric case.
Lastly, the influence of an increased heat conductivity with α 1 = 4 may be taken from figure 7. It is visible that the first phase reaches T 1 (x) = 0 faster for an increased value of α 1 compared to the symmetric case α 1 = 1, which is to be expected, since the heat can now be transported faster. This results in a faster decaying temperature distribution. Overall the temperature shows values closer to zero for later times compared to the balanced state for an increased value of α 1 , while the motion of the interface towards the second phase is also clearly visible. Once again the temperature at the solid boundary of the second phase did not change compared to the symmetric case.

Conclusion and outlook
In this paper the unified transform method was applied to the one-dimensional Stefan problem for two phases in order to obtain a solution for the temperature distributions for each phase and, moreover, for the interface motion. In the second step the resulting non-linear integrodifferential equations were solved numerically. A combination of the fix-point iteration of the phase boundary motion at a specific time step and spline interpolation over those lead to the solution for the evolution of the phase boundary.
The obtained results for the motion of the interface for the one-phase case were compared to literature, which showed very good agreement, in order to test the implemented solver. As another test case for the solver symmetric initial conditions were employed, which retained successfully a symmetric solution for all times. Furthermore, the influence of modified initial parameters such as initial temperatures and thermal diffusivities on the evolution of the interface and the temperature distribution were investigated. The temperature distributions at any given point in time and space can be recovered from the solution of the phase boundary. The numerical implementation gives an order of convergence of about 1.5. Thus the obtained results prove that the two-phase Stefan problem can be solved using the presented approach without explicitly treating the phase boundary.
An obvious extension of this work would be the consideration of the two-dimensional and even the three-dimensional Stefan problem, which would have a closer applicability to physical science and engineering. Due to the importance of heat conduction problems in nature and engineering, this investigation could be useful for a variety of industrial processes. Thus an extension to even more phases and improving the numeric methods and thereby covering for a bigger scope of problems in engineering and science may be the focus of future studies. The present work further shows that the UTM can be applied on free boundary problems, such as the Stefan problem. Thus the application of the UTM on other free boundary problems e.g. the welding problem may be of interest in future works.
Taking the limit lim sup x→g(t)−0 of the equation above and using (42) and (43) leads to (45).
One can easily see that (46) approaches 0 for lim x→g(t) , and thus satisfies lim x→g(t) Adding (45) and (47) we get: lim sup Since the left side of (48) is independent of σ we conclude, on taking σ → 0, that it is zero, and the proof of the lemma is completed.