A meshless method for an inverse two-phase one-dimensional linear Stefan problem

We extend the meshless method of fundamental solutions proposed in [B.T. Johansson, D. Lesnic, and T. Reeve, A method of fundamental solutions for the one-dimensional inverse Stefan problem, Appl. Math. Model. 35 (2011), pp. 4367–4378] for the one-dimensional one-phase inverse Stefan problem to the two-phase change case. The implementation and analysis are more complicated and meaningful since one needs to handle composite layered materials. Furthermore, the inverse problem is ill-posed since small errors in the input data lead to large deviations in the solution. Therefore, the inverse problem is intractable to classical methods of linear inversion, and regularization is employed in our study. Numerical results obtained when the input data is either exact or contaminated with random noise are compared with the analytical solution, where available, or with the numerical results obtained by other methods, [D.D. Ang, A. Pham Ngoc Dinh, and D.N. Tranh, Regularization of an inverse two-phase Stefan problem, Nonlinear Anal. 34 (1998), pp. 719–731], otherwise.


Introduction
Heat conduction problems with phase change of state are of importance in several practical applications in metallurgy, ingot solidification, continuous casting of steel, scrap melting, ablation of metals using laser beams, welding of two materials, etc. [1]. All these involve solving (in general, numerically) some Stefan-type problems for the parabolic heat equation. The classical Stefan problem has the essential specific peculiarity of imposing simultaneously both the local thermodynamic equilibrium and the dynamical compatibility condition on the unknown free boundary [2].
The direct Stefan problem requires determining both the temperature and the moving free surface when the initial and boundary conditions, as well as the thermal properties of the medium involved, are known. In contrast, in inverse problems some part of the cause of a physical phenomenon is unknown and it is compensated for by the knowledge of some part of the effect. Consequently, inverse Stefan problems require determining some initial temperature and/or boundary conditions, and/or thermal properties from additional information, which may involve the partial knowledge of the free surface, its normal velocity or the interior temperature measured at some points inside the medium [3]. For example, in the process of recrystallization one has to solve the inverse Stefan problem which requires determining the temperature and the heat flux at the fixed surface that control the flatness of the crystallization front [4].
Prior to this study, much of the work on the inverse Stefan problem has been done for one-phase flow, see e.g. [5][6][7][8][9][10]. However, the literature on solving numerically the twophase inverse Stefan problem is much more scarce, see [11][12][13]. A reason for this, as it has been pointed out earlier in the literature, is that standard direct problem domain discretization methods such as the finite-element method and the finite-difference method are not straightforward to adjust to the case of internal (interface) boundaries. Instead, integral equation-based methods are more attractable to use for these type of problems. In [11], the problem was reduced to a Volterra integral equation of the second kind and solved by regularization; however, their numerically obtained solution appears rather unstable. We also mention that, more recently, [12] employed the alternating phase truncation method combined with genetic algorithms and Tikhonov regularization for identifying a piecewise constant convective heat transfer coefficient. However, this latter heat transfer coefficient identification problem will not be investigated herein. Instead, we shall concentrate on extending the numerical method developed by the authors in [5] for the one-dimensional one-phase inverse Stefan problem to the two-phase case. We shall also compare our numerical results to those obtained in [11] by using the different method briefly abstracted above.
The plan of this article is as follows. In Section 2, we give the mathematical formulation of the one-dimensional two-phase inverse Stefan problem and point out its illposedness. In Section 3, we describe the numerical meshless method which is employed to approximate the solution of the inverse problem and its regularization. Section 4 presents and discusses numerical results obtained for some benchmark test examples. Noise is added to the input data in order to investigate the stability of the numerical solution. In one example, the analytical solution is available, whilst in the second one we compare our numerical results to those obtained in [11], where a different method was employed. Finally, Section 5 presents conclusions and possible future work.

Mathematical formulation
Consider a one-dimensional two-phase Stefan model for a slab of thickness l 4 0 consisting partly of water and partly of ice separated by the moving interface x ¼ s(t). Mathematically, the problem consists of determining two functions u 1 and u 2 satisfying the one-dimensional linear heat equations subject to the initial conditions at t ¼ 0 the Stefan interface conditions along the curve x ¼ s(t), and a Dirichlet or Neumann boundary condition at x ¼ l, i.e. u 2 ðl, tÞ ¼ f l ðtÞ or k 2 @u 2 @x ðl, tÞ ¼ g l ðtÞ, 05 t T: ð7Þ Figure 1 depicts the two-phase inverse Stefan problem under investigation. We mention that the corresponding one-phase inverse Stefan problem investigated in [5] is given by Equations (1) and (3), and s 0 ðtÞ ¼ K 1 @u 1 @x ðsðtÞ, tÞ, 05 t T: In the above model, T 4 0 is a given final time of interest, i ¼ k i /( i c i ) are the thermal diffusivities of water (i ¼ 1) and ice (i ¼ 2), k i , i and c i are the thermal conductivities, densities and specific heats, respectively, K i ¼ k i /( 2 L) where L is the latent heat, and all the preceding physical quantities are constant, positive and given. The functions u 0 1 and u 0 2 represent the given initial data and satisfy the compatibility conditions at x ¼ s(0), namely where s(0) is assumed to be in the interval (0, l). We also require the compatibility condition at x ¼ l, namely in the case of the Dirichlet or Neumann data in (7), respectively. The Stefan boundary conditions (5) and (6) reflect, respectively, the facts that the temperature at the interface must be equal to the initial temperature u Ã (t) (usually taken to be zero in melting processes, for simplicity) and that the rate of melting is proportional to the rate of absorption of heat energy at the interface [14].
Denote by the two-phase rectangular solution domain (0, l) Â (0, T ] which is subdivided into two subdomains and D 2 T ¼ fðx, tÞ 2 Q T j sðtÞ 5 x 5 l, 0 5 t T g: The direct Stefan problem then requires finding the temperatures u i 2 C 2,1 ðD i T Þ, i ¼ 1, 2, and the interface s 2 C([0, T ]) \ C 1 ((0, T ]) satisfying s(t) 2 (0, l) for 0 t T, Equations (1)-(11) and a boundary condition at x ¼ 0, namely where f 0 , or g 0 , or E are given functions satisfying the compatibility conditions The well-posedness of the above direct two-phase Stefan problems were accomplished by Cannon and his co-workers in [15][16][17]. The infinite differentiability of the free surface function x ¼ s(t) has also been established in [18]. In contrast to the above direct nonlinear and well-posed Stefan problem, which requires finding the solution (u 1 , u 2 , s) satisfying Equations (1)-(13), the inverse linear and ill-posed Stefan problem requires finding the solution (u 1 , u 2 ) satisfying Equations (1)- (11) only, but the additional information is that the interface x ¼ s(t) is now known. The problems in which the additional information is the position of the moving interface are called design problems [12]. In particular, all the quantities defined in Equation (12) and satisfying (13) are unknown and need to be determined. The more difficult inverse nonlinear Stefan problem which requires finding the solution (u 1 , u 2 , s) satisfying Equations (1)-(5), (10) and the Cauchy data at x ¼ l (compare Equations (7)-(11)), namely with u 0 2 ðl Þ ¼ f l ð0Þ and k 2 @u 0 see [19], is deferred to a future work.
An initial attempt would be to split the two-phase inverse Stefan problem (1)-(11) into a direct problem for u 2 in D 2 T followed by a one-phase inverse Stefan problem for u 1 in D 1 T : More specifically, one can first solve, using for example the method of fundamental solutions (MFS) presented in [20], the direct well-posed problem for u 2 in D 2 T satisfying Equations (2), (4), (7), (11) and where and given s(t) 2 (0, l) for 0 t T, to determine the partial derivative This output can then be introduced into (6) to yield Then, Equations (1), (3), (19) and satisfying and given s(t) 2 (0, l) for 0 t T, form a one-phase inverse linear and ill-posed Stefan problem which has been previously solved using the MFS by the authors in [5]. The above splitting is useful in concluding the uniqueness of the solution of the twophase one-dimensional inverse Stefan problem (1)- (11). Namely, the direct problem given by Equations (2), (4), (7), (11), (16) and (17) has a unique solution u 2 , see [21], whilst the one-phase one-dimensional inverse Stefan problem given by Equations (1), (3), (19)-(21) also has a unique solution u 1 , see [22]. Moreover, uniqueness of a solution of this latter inverse problem also holds even if the initial condition (3) is not prescribed but is unknown, see [23].
However, when the data u 0 2 , u Ã , f l or g l are contaminated with noise the retrieval of the heat flux derivative (18) is numerically unstable, see [24,25]. Moreover, the inverse onephase Stefan problem for u 1 given by Equations (1), (3), (19)-(21) is ill-posed since small errors in the input data u 0 1 , u Ã , s cause large errors in the output data (12), see [6,9,26]. Thus, in order to obtain stable solutions, both the above problems need to be regularized; however, using Tikhonov regularization this involves choosing two different regularization parameters. It would seem better if we could have instead a non-splitting approach of solving Equations (1)-(11) simultaneously using a regularization procedure which involves only choosing a single regularization parameter. This combined approach has been used before by the authors for simultaneously determining a heat source and the initial temperature in the standard parabolic heat equation, see [27]. Therefore, in this article, we address finding the solution (u 1 , u 2 ) by solving Equations (1)-(11) simultaneously. We will also investigate the case when the initial condition (3) is not prescribed, see the corresponding one-phase case, which has been recently investigated by the authors in [28].

The MFS
We approximate the solution (u 1 , u 2 ) of the heat equations (1) and (2) using the MFS for the transient heat conduction in layered materials, recently developed by the authors in [29,30]. Let where H the Heaviside function be the fundamental solutions of the heat equations (1) and (2). In the MFS, we approximate u 1 and u 2 by linear combinations of fundamental solutions (22) of the form In expressions (23) and (24), the source points ð y ðnÞ j , j Þ for n ¼ 1, 2 and j ¼ 1, M, are located outside their corresponding solution domains in the following way: where A representation of the domains, boundaries and possible placement of source points is shown in Figure 2. Figure 2. Representation of the two-phase problem with domains D 1 T and D 2 T , unknown boundary data (Â), and source points (ÁÁÁÁÁ) and (---) placed externally to the domains D 1 T and D 2 T , respectively.
The MFS approximations (23) and (24) are motivated by the denseness properties shown in [5,29] in the standard L 2 -spaces on the respective boundaries of the solution domains. Alternatively, one could choose points below the bottom surface at points with coordinates ð2j À 1Þsð0Þ 2M , Àh and ð2j À 1Þðl À sð0ÞÞ as suggested in [31] (under-determined situation), or at both these latter points and the top half of the lateral points (25) and (26), as suggested in [32] (over-determined situation). However, in this article we employ only the more natural determined situation given by the selection (25) and (26), and for which we know denseness properties of the approximation due to the proofs provided in [5,29].

Determination of the coefficients in the MFS
We now discuss how to determine the coefficients c ðnÞ j for n ¼ 1, 2, and j ¼ 1, M: We shall collocate the initial conditions (3) and (4), the interface conditions (5) and (6) and the boundary condition (7) at certain collocation points. For the ease of implementation we select a uniform distribution of collocation points given by and collocate Equations (3)-(11) as follows: Note that it is not important whether we impose (31) at i ¼ 0 or not. In total, via (23) and (24), Equations (29)-(32) form a system of (4M 2 þ 4 þ 2K) linear equations with 2M ¼ 8M 1 unknown coefficients c ¼ ðc ðnÞ j Þ j¼1,M, n¼1, 2 : A necessary condition for a unique solution is that Based on (23) and (24) we can write the system of Equations (29)-(32) in matrix form as with the obvious notation. As pointed out in [25,33] the matrix A will have a high condition number both because of the ill-conditioning of the numerical MFS used and the ill-posedness of the mathematical problem. In order to handle this lack of stability, we consider the Tikhonov regularization method, which solves the modified stabilized system of linear equations instead of the ill-conditioned system of Equations (34), where the superscript tr denotes the transpose of a matrix, I is the identity matrix and 4 0 is a regularization parameter to be chosen. The now well-conditioned linear system of equations (35) can be solved using any classical method such as the Gaussian elimination method. The regularization parameter can, for example, be chosen according to the L-curve criterion, [34]. We note that solving the system of Equations (35) yields an approximation to the vector of coefficients c which in turn, via the MFS approximants (23) and (24), produces the desired missing boundary data (12).

Numerical results and discussion
In this section, we apply the MFS outlined in the previous section to the inverse two-phase Stefan problem (1)- (11). Preliminary investigations indicate that the MFS parameter h, giving the distance between the boundary and the source points, should be chosen neither too small nor too large and in this section we present the numerical results obtained with a typical value of h ¼ 2. Other values of h in the range (0.5, 4) produced similar results.
In the first example we have an analytical solution available. For this example we also investigate the case when the initial condition (3) for the function u 1 is not given and is part of the unknowns. In the second example, since an analytical solution is not available, we compare our numerical results with those obtained in [11] where a different method was used.
One can also check that the condition (6) is satisfied. Finally, the boundary condition at x ¼ l ¼ 1 is taken as The analytical solution of the problem (1), (2), (6), (36), (37), (39) and (40) is then given by, [29,35], which can be verified by direct substitution. Of particular interest is to recover the data (12) at the inaccessible (hostile) boundary x ¼ 0, given by As always with inverse and improperly posed problems, behaviour of the numerical solution in the presence of noise is an important part of the stability testing process. Noise could be added in the initial data (36), (37) or the moving boundary (38), see [12], but this case is not pursued herein. In our example, considering only the Dirichlet boundary condition in (40) given by we add to it random additive noise generated as where N(0, 2 ) represents the normal distribution with mean zero and standard deviation where represents the percentage of noise level. It was assumed that the initial data (36), (37) and the moving boundary (38) would be known fairly accurately, and the measurement error would concentrate mostly in the boundary temperature (45) or the heat flux (49) part of the data.  From these figures it can be seen that more accurate approximations can be obtained by increasing the number of collocation and source points. Figures 4(a) and (b) show the plots of the exact solution and exact normal derivative, respectively, in comparison with the MFS approximations for various noise levels added in the input data (45). From these figures one can see that the error in the approximation is of the same order as the noise level everywhere except for time t close to the final time T ¼ 1. This conclusion is consistent with previous studies on the inverse heat conduction problems, see [36][37][38], which showed that if one wishes to determine a stable and accurate approximation for the boundary temperature and the heat flux in (10) for all times closer to the upper final time t ¼ T, then one needs to consider and use the input Cauchy data (12) over a slightly future time interval (0, T þ r), where r 4 0 is related to the concept of Beck's future temperatures, see [39].
In Figure 5(a) the exact solution and the MFS approximation for the best, average and least accurate approximations after 10 trials are plotted, and in Figure 5(b) the absolute error is plotted for the noise level ¼ 5%.
Next, we look at the case when the initial condition (36) is not imposed. In this case, Equation (29) is collocated only for n ¼ 2 and this together with (30)-(32) then forms a system of (4M 2 þ 4 þ K) linear equations with 2M ¼ 8M 1 unknown coefficients c ¼ ðc ðnÞ j Þ j¼1,M, n¼1, 2 . A necessary condition for a unique solution is that compare with (33). To further highlight the flexibility of the proposed MFS approximation, in the remaining figures in this example we apply the Neumann boundary data in (40), instead of the Dirichlet data (44), given by @u 2 @x ð1, tÞ ¼ g 1 ðtÞ ¼ À We further add to it random additive noise generated as where ¼ Â max  initial data at t ¼ 0 for 0 x s(0), respectively. All plots in this example have shown the flexibility of the MFS to handle different and noisy boundary data and produce accurate and stable results.

Example 2
In this example, taken from [11], an analytical solution is not available. The input data are given by Noise is added in the zero Neumann data (52) as We also tested adding noise to the initial data (51), but this appeared to be less influential on the obtained approximation. Note that, in the similar example in [11], the initial data (3) was additionally imposed as However, it is unlikely that a solution to this inverse problem exists. Moreover, the plot presented in [11] is not clearly explained and since their numerically obtained solution appears not to match the given initial condition at x ¼ t ¼ 0, this seems to indicate the lack of a solution. Therefore, we relaxed the problem and took the initial data u 0 1 to be an additional unknown function to be determined. Since, as it is highlighted at the end of Section 2, the existence and uniqueness of the solution in D T 2 is guaranteed, from the given Cauchy data on the interface s(t) it should be possible, due to analytic continuation [22], to determine a solution in D T 1 provided that there are no additional redundant boundary or initial conditions imposed. This is further indicated by the numerical results below, in which a stable numerical solution with the expected physical properties is obtained. The inverse problem given by Equations (1), (2), (5), (6), (51) and (53) is then solved. In Figure 8(a) and (b), we plot the MFS approximations when T ¼ 1 and T ¼ 7, respectively, for different levels of noise applied in (53). The plots show that the temperature decays to zero as t increases, which is expected due to the diffusive nature of heat transfer. Figures 9(a)   plots of the MFS approximations for the initial data (3). In both plots we observe that stable and accurate results are obtained. Finally, for completeness, Figure 10 shows the MFS approximation for all (x, t) 2 [0, /2] Â [0, 1].

Conclusion
In this article, a regularized MFS has been developed for solving an inverse two-phase onedimensional linear Stefan problem. The numerical results show that the method is accurate and stable with respect to noise added in the input data. Future work will concern extending the MFS developed in this study to three-phase and two-dimensional inverse Stefan problems, see [26] and [40], respectively.