The Solution of Two-Phase Inverse Stefan Problem Based on a Hybrid Method with Optimization

The two-phase Stefan problem is widely used in industrial field.This paper focuses on solving the two-phase inverse Stefan problem when the interface moving is unknown, which is more realistic from the practical point of view. With the help of optimization method, the paper presents a hybrid method which combines the homotopy perturbation method with the improved Adomian decomposition method to solve this problem. Simulation experiment demonstrates the validity of this method. Optimization method plays a very important role in this paper, so we propose a modified spectral DY conjugate gradient method. And the convergence of this method is given. Simulation experiment illustrates the effectiveness of this modified spectral DY conjugate gradient method.


Introduction
The process of solidification/melting occurs in many industrial fields, such as steel-making continuous casting.In the continuous casting process [1,2], molten steel comes into slab by solidification.And this process is shown in Figure 1.Because of the heat symmetry, we only consider 1/2 of slab thickness, and we divide this domain into two regions by the freezing front, which is shown in Figure 2.  1 stands for the solid phase and  2 is taken by the liquid phase.The solidification problem of continuous casting is associated with Stefan problem, which is a model of the process with a phase change.
Many researches have focused on solving Stefan problems [3][4][5][6][7].Grzymkowski and Slota [8,9] applied the Adomian decomposition method (ADM) combined with optimization for solving the direct and the inverse one-phase Stefan problem, and Hetmaniok et al. [10] used the Adomian decomposition method (ADM) and variational iteration method to solve the one-phase Stefan moving boundary problems.In comparison to the studies on the one-phase Stefan problems, researches on two-phase Stefan problems [11] are much more scarce.Słota [12] applied homotopy perturbation method (HPM) to solve the two-phase inverse Stefan problem.However, in their formulation, these two-phase inverse Stefan problems were solved when the position of the moving interface is known.While the position of the moving interface is unknown and no temperature or heat flux boundary conditions are specified on one part of the boundary, this twophase inverse Stefan problem may be difficult to be solved.This problem is close to the problem which occurs in the continuous casting.Therefore, this paper focuses on solving this problem by a hybrid method with optimization which combines the homotopy perturbation method (HPM) with the improved Adomian decomposition method (IADM).
This paper is organized as follows.In Section 2, we give the description of the two-phase Stefan problem.And the hybrid method for solving this problem is introduced in Section 3. We give a modified spectral DY conjugate gradient method and the convergence of this method in Section 4. In Section 5, some numerical results are given to illustrate the validity of the hybrid method.The solidification process of continuous casting slab ( stands for the direction of heat transfer and  stands for casting direction).

The Description of Two-Phase Inverse Stefan Problem
In this section, the two-phase inverse Stefan problem is considered.The two domains are described (see Figure 2): and their boundaries with unknown function  = () are Γ 1 = {(, 0) ;  ∈ (0, ) ,  =  (0)} , Γ 2 = {(, 0) ;  ∈ (, ) ,  =  (0)} , In practice, the measurement of temperature on the boundary Γ 4 may not be easily obtained and the position of the moving interface described by means of function  = () is unknown.Therefore, the designed two-phase inverse Stefan problem consists in the calculation of the temperature distribution in the domains  1 and  2 , respectively, as well as in the reconstruction of the functions describing moving interface  = ().In order to solve this inverse problem, the additional information is given.In such a situation, the final time internal temperature at  =  * on the boundary Γ 6 is known, which can be measured in practice.The temperature distribution is determined by the functions  1 (, ) and  2 (, ), which are defined in the domains  1 and  2 , respectively, and they satisfy the following heat conduction equations: where  and  refer to time and spatial location and  1 and  2 are the thermal diffusivity in the solid phase and the liquid phase, respectively.On the boundaries Γ 1 and Γ 2 , the initial conditions are described in the following: On the boundary Γ 3 , it satisfies the Dirichlet boundary or Neumann boundary condition: On the boundary Γ 4 , it satisfies the Dirichlet boundary conditions: On the phase change moving interface (Γ 5 ), it satisfies the condition of temperature continuity and the Stefan condition: where  * is the phase change temperature.Hence, where  is the latent heat of fusion per unit volume and  1 and  2 are the thermal conductivity in the solid and liquid phase, respectively.The final time internal temperature on the boundary Γ 6 satisfies

The Solution of Two-Phase Inverse Stefan Problem
3.1.The Homotopy Perturbation Method.The homotopy perturbation method was presented in [13,14]; it arose as a combination of two other methods: the perturbation method and the homotopy technique from topology.Homotopy perturbation method was applied on differential equations, partial differential equations, and the direct and inverse problems of heat transfer by many authors [15][16][17].
Homotopy perturbation method was used to find the exact or the approximate solutions of the linear and nonlinear integral equations [16,18], the two-dimensional integral equations [19][20][21], and the integrodifferential equations [22][23][24][25].This method enables seeking a solution of the following operator equation: where  is an operator,  is a sought function, and  is a known function on the domain of Ω.The operator  is presented in the form of sum: where  defines the linear operator and  is the remaining part of operator .So (18) can be written as We define a new operator  which is called the homotopy operator in the following: where  ∈ [0, 1] is an embedding parameter, V(, ) : Ω × [0, 1] → R, and  0 defines the initial approximation of a solution of (18); ( 21) is transformed into the following form: Therefore, For  = 0 the solution of operator equation (V, 0) = 0 is equivalent to solution of problem (V) − ( 0 ) = 0, whereas, for  = 1, the solution of operator equation (V, 1) = 0 is equivalent to solution of problem (20).Thus, a monotonous change of parameter  from zero to one is just that continuous change of the trivial problem (V) − ( 0 ) = 0 to the original problem.In topology, this is called homotopy.According to the HPM, we assume that the solution of ( 21) and ( 22) can be written as a power series in : setting  = 1, the approximate solution of ( 20) is obtained: in most cases, the series in ( 25) is convergent, and in [15] an additional remark on the convergence of this series can be found.

The Improved Adomian Decomposition
Method.The Adomian decomposition method was presented by Adomian [26].This method was used to solve a wide variety of linear and nonlinear problems.Several other researchers had developed a modification to the ADM [27].The modification usually involves a slight change and it is aimed at improving the accuracy of the series solution.
The Adomian decomposition method is able to solve a wide class of nonlinear operator equations [28] in the form where  :  →  is a nonlinear operator,  is a known element from Hilbert space  and  is the sought element from Hilbert space .The operator () can be written as where  is an invertible linear operator,  is a linear operator, and  is a nonlinear operator.The solution of ( 26) admits the decomposition into an infinite series of components the nonlinear term  is equated to an infinite series is the Adomian polynomials, which can be determined by Substituting ( 28), ( 29), (30), and (31) into operator equation ( 26) and using the inverse operator  −1 , the following equation can be obtained: where  * is the function dependent on the initial and boundary conditions.The improved Adomian decomposition method [27]: it is important to note that the improved Adomian decomposition method is based on the assumption that the function  =  * +  −1 () can be divided into two parts, namely,  1 and  2 .Under this assumption, the following equation can be obtained: (33) Accordingly, a slight variation is proposed on the components  0 and  1 .The suggestion is that only the part  1 is assigned to the component  0 , whereas the remaining part  2 is combined with other terms given in (33) to define  1 .Consequently, the recursive relation is written as follows: (34)

The Solution of the Two-Phase Inverse Stefan Problem.
In this section, we introduce the solution of the two-phase inverse Stefan problem.We split the two-phase inverse Stefan problem into two problems.The first is the determination of  1 (, ) in domain  1 and the moving interface () which is satisfying ( 7), ( 8), (10), and (13).Once the boundary  = () and the heat flux ( 2 (, )/)| =() have been obtained, the second problem for determining the temperature  2 (, ) in domain  2 is solved.The solution of  1 (, ) and moving interface () based on HPM: because the heat flux on the moving position is unknown in this two-phase inverse Stefan problem, the HPM is used to compute  1 (, ).In this section, the homotopy perturbation method is applied to solve the inverse problem.We make the homotopy map for (22): the solution to this equation, is sought in the form of a power series in  and it satisfies where  ∈ [0, 1]; by substituting (37) into (36), the following equation can be obtained: Next, by comparing the same powers of parameter  in (38), the following equation can be obtained: Equations ( 39) and (40) must be supplemented by the boundary conditions.Therefore, for (39), we set the following boundary conditions: and for (40) we set the conditions in the following ( ≥ 2): The above conditions are selected.So we give the -order approximate solution: because  1, are determined by the unknown function ().This function is derived in the form Thus we define a cost function  1 according to (10) and ( 17): Therefore,   are chosen according to the minimum of (47).
The moving interface () and the solution of  1 in domain  1 can be obtained.The boundary condition (13) on the boundary Γ 3 and the condition of temperature continuity (15) should be fulfilled for  ≥ 1.In this way, while looking for the solution of this problem in domain  1 , the initial approximation  1,1 will be found by solving the system of (39) with conditions (41) and (42).Conversely, functions  1, ,  = 1, 2, . .., are determined by the recurrent solution of (40) with conditions (42) and (44).Consequently, the moving interface () and the solution of  1 in domain  1 can be obtained by the homotopy perturbation method with optimization.According to the Stefan condition (16) on the phase change moving interface Γ 5 , the temperature  2 (, ) in domain  2 can be recovered by IADM with optimization method.
The solution of  2 (, ) based on IADM with optimization: in this problem, we consider the operator equations (27), where Mathematical Problems in Engineering 5 The inverse operator  −1 is given by the following equation: The boundary condition and the Stefan condition from ( 14) are used to obtain the following equation: Hence then the following recursive relation can be obtained: because  1 (, ) and () are determined by the above HPM with optimization method, ( 2 (, )/)| =() can be obtained by (16).An approximation solution is expressed in the form and in this inverse Stefan problem, () is unknown on the boundary Γ 4 and it is derived in the following form: where   ∈  and   () are the basis function.The coefficients   are selected to show minimal deviation of function (55) from conditions (11) and (15).Thus   are chosen according to the minimum of the following functional: where  2 ( 1 , . . .,   ) is a nonlinear function.The optimization method plays a very important role in solving this twophase inverse Stefan problem.Therefore, a modified spectral DY conjugate gradient method is presented.In Section 4, this optimization method is introduced.
Remark.The IADM with optimization method can be used to solve the direct and inverse one-phase Stefan problem which is designed in the literature [8,9].And this method is better than ADM with optimization method.In Section 5, we give two examples of one-phase Stefan problem to illustrate the effectiveness of the IADM with optimization method.

A Modified Spectral DY Conjugate Gradient Method
Conjugate gradient methods are very efficient for solving the unconstrained optimization problem min ∈  ().Their iterative formula is given as follows: where step-size   is positive,   = ∇(  ), and   is a scalar.In addition,   is a step length which is calculated by the line search.  is a very important factor in conjugate gradient methods; the FR, PRP, DY, and CD methods [29] are used for choosing this scalar.Among them, the DY method is regarded as one of the efficient conjugate gradient methods.The scalar   is given by where  −1 =   −  −1 and ‖ ⋅ ‖ stands for Euclidean norm.Recently, Birgin and Martłnez [30] presented a spectral conjugate gradient method which is combining conjugate gradient method and spectral gradient method.And the direction   is obtained by where and   is regarded as spectral gradient.However, according to [31], the spectral conjugate gradient method can not guarantee producing the descent directions.
Based on the above analysis, in the following, we describe a modified spectral DY conjugate gradient method and give the convergence of this method with the Wolfe type line search rules.

The Modified Spectral DY Conjugate Gradient Method.
The new method is introduced: where   = 1 + ‖ −1 ‖ 2 /‖  ‖ 2 and   =  DY  .This method is called the MDY method.Algorithm 1.Consider the following.
Step 1.The following are given:  0 ∈   ,  > 0, and the max iteration step  max .
Step 3. Find   according to Wolfe type line search rule.

Convergence Analysis.
From the above, the MDY conjugate gradient method provides a descent direction for the objective function.In the following, we give the convergence of this algorithm with the Wolfe type line search rules.Firstly, we introduce some assumptions: (67) If we choose  = 1, then according to (65), we have So  DY 2 ≥ 0. For  > 1, according to (67), we have by induction.So we obtain the proof is completed.
Theorem 3. Let the sequences {  } and {  } be generated by the proposed algorithm (Algorithm 1) with Wolfe type line search rules; then where  > 0 is a constant.

Numerical Simulation Experiments
In order to illustrate the validity of the IADM with optimization method, we introduce the one-phase inverse Stefan problem in Example 1. From the literature [9], we know the inverse Stefan problem consists in the calculation of the temperature distribution in the domain and in the reconstruction of the temperature distribution on the boundary, when the position of the function which describes the moving interface is known.And then, in Example 2, we solve the twophase inverse Stefan problem based on the hybrid method which combines the Adomian decomposition method with the homotopy perturbation method.
Example 1.In this example, the one-phase inverse Stefan problem [9] is taken for the following values of parameters: Then the exact solution of such formulated one-phase Stefan problem can be found from the following functions: The following equation, is taken as basis function.The approximate solutions are compared with the exact solution and the values of the absolute errors are calculated from where ℎ() ∈ {(), ()} (() stands for the boundary condition on Γ 1 in the literature [9] and () is the moving interface),  is the domain of the temperature distribution, ℎ  () is the exact value of function ℎ(), ℎ  () is the approximate value of function ℎ(),   (, ) is the exact distribution of the temperature in the domain, and   (, ) is the reconstructed distribution of temperature.Percentage relative errors are calculated from The comparison of the ADM and the IADM for reconstructing distribution of the temperature on boundary Γ 1 in the literature [9] is shown in Figure 3 for  = 1,  = 2 and  = 2,  = 3.From Figure 3, the method based on the improved Adomian decomposition method (IADM) is better than the Adomian decomposition method (ADM).Example 2. This example is introduced to illustrate the effectiveness of the hybrid method which combines the Adomian decomposition method with the homotopy perturbation method, which is presented in the previous sections.Here, we consider an example of the two-phase inverse Stefan problem [12], in which  (93) First, we use HPM with optimization method to solve  1 (, ) and ().As the initial approximation  1,0 , this fulfills the initial condition: By using the function  1,0 in (39) and the conditions (41) and (42), the following system of equations can be obtained: With the conditions  1,1 (0, ) =  3/10 ( (/10)−1 ) , this system of equations can be solved.And the solution of this system is shown in the following: 1,1 (, ) = −  3/10−/5 + 1 −  (3+)/10  ()  +  (+3)/10 .(97) Next, we designate the functions  1, for  ≥ 2 in the following: and the conditions are shown: By using (98) and the boundary conditions (99), we derive  1, (, ) for  ≥ 2; according to (45), the approximate solution  1 (, ) is obtained.Here, we define basis functions according to (46): We use the modified spectral DY (MDY) conjugate gradient method to solve function (47).The solution of moving interface () is shown in Figure 4.So the temperature  1 (, ) is derived in domain  1 .Values of the error in reconstruction of the position of the moving interface () and the temperature distribution  1 (, ) are shown in Table 1.
In order to verify the validity of the modified spectral DY (MDY) conjugate gradient method, we choose  = 2,  = 2,  = 10 −4 , the initial point  0 = (2.2,2.5), and the max step number  max = 100.The comparison results of DY and MDY are shown in Figure 5 and Table 2. Comparing to DY conjugate gradient method, the MDY conjugate gradient method is slightly better.
The moving interface () and  1 (, ) in domain  1 are obtained by the homotopy perturbation method with optimization.Next, with the help of condition ( 16), we can determine the temperature  2 (, ) in domain  2 by the improved Adomian decomposition method (IADM) with optimization.Without loss of generality, we choose the basis functions The exact and the reconstructed distribution of the temperature on boundary Γ 4 are shown in Figure 6 for different number of basis functions   ().And the errors in reconstruction of the temperature () (on boundary Γ 4 in Figure 2) and the temperature distribution  2 (, ) are shown in Table 3.The inversion results of the function () match the exact value very well.
The measured temperature   1 () at the final time  * on the boundary Γ 6 in Figure 2 is perturbed by the random error.So in the next section we consider the influence of measurement errors on the results.We set the measured The ADM (m = 2, n = 1) The IADM (m = 2, n = 1) temperature   1 () with an error of 1% and 3%.We can obtain the moving interface () and (); the results are shown in Figures 7 and 8.
In Figures 7 and 8, the exact and the reconstructed distribution of the temperature on moving interface () and the boundary Γ 4 are shown, and the errors in reconstruction of the function describing the moving interface and boundary conditions are compiled in   method with optimization is stable with the input errors.
When the input data are contaminated by the errors, the error of the reconstructed boundary conditions does not exceed the error of the input data.

Conclusion
Based on the homotopy perturbation method (HPM) with the improved Adomian decomposition method (IADM), a hybrid method with optimization is proposed to solve the two-phase inverse Stefan problem.This problem is more realistic from the practical point of view.The simulation experiment is used to verify this method.Experimental results match the exact value very well.Furthermore, this investigated method can also be used for solving one-phase direct Stefan problems.Optimization method is very crucial in this paper, so a modified spectral DY conjugate gradient method is    presented.And we give the convergence of this MDY method.Simulation experiment illustrates the validity of this optimization method.

2 Mathematical
Figure 1: The solidification process of continuous casting slab ( stands for the direction of heat transfer and  stands for casting direction).

Figure 2 :
Figure 2: The domain formulation of the two-phase Stefan problem.

Figure 5 :Figure 6 :
Figure 5: The cost function  1 changes with the increasing of iteration numbers.

Figure 8 :
Figure 8: The reconstructing distribution of the temperature () on boundary Γ 4 in Figure 2.

Table 1 :
Values of the error in reconstruction of the position of the moving interface and the temperature distribution  1 (, ).

Table 4 .
The obtained reconstruction results illustrate that the functions () and () are reconstructed very well.From the results, the hybrid

Table 2 :
The comparison results of DY and MDY.

Table 4 :
Values of the error in reconstruction the moving interface () and the temperature () on boundary Γ 4 in Figure2.