An Accurate Approximation of the Two-Phase Stefan Problem with Coefﬁcient Smoothing

: In this work, we consider the heat transfer problems with phase change. The mathematical model is described through a two-phase Stefan problem and deﬁned in the whole domain that contains frozen and thawed subdomains. For the numerical solution of the problem, we present three schemes based on different smoothing of the sharp phase change interface. We propose the method using smooth coefﬁcient approximation based on the analytical smoothing of discontinuous coefﬁcients through an error function with a given smoothing interval. The second method is based on smoothing in one spatial interval (cell) and provides a minimal length of smoothing calculated automatically for the given values of temperatures on the mesh. The third scheme is a convenient scheme using a linear approximation of the coefﬁcient on the smoothing interval. The results of the numerical computations on a model problem with an exact solution are presented for the one-dimensional formulation. The extension of the method is presented for the solution of the two-dimensional problem with numerical results. the solution with increasing mesh parameters m and n to the reference solution on the ﬁne grid m = n = 400 using the corresponding method. We observe that both errors e f ( u T ) and e e ( u T ) are sensitive to the size of the mesh, where we obtain a large error for a very coarse grid n 1 = n 2 = 25. We obtain accurate results for n = 50, 100, and 100 using the h-smoothing method.


Introduction
One of the urgent problems of mathematical physics is the Stefan problem, which describes heat transfer with a phase change of matter. An example of physical processes with phase transitions is the problem of melting ice with a movable boundary between water and ice, the problem of melting a solid with an unknown boundary between the solid and liquid phases, and the heat transfer processes in permafrost areas. The Stefan problem is an initial-boundary value problem of a parabolic differential equation with discontinuous coefficients on the phase change interfaces. The phase change occurs for a given value of temperature (freezing point), where the energy balance on the interface is written in the form of the Stefan condition and represented as a sharp interface (a moving-boundary problem).
In the numerical solution of the Stefan problem, a major difficulty is associated with the tracking of the phase change interface evolution. The enthalpy method was presented in [1][2][3][4], where the concept of effective heat capacity was introduced through various methods of smoothing of discontinuous coefficients. The general approach of thermodynamic functions was presented in [5], where an online application was implemented to generate the first-order derivatives of the thermodynamic parameters for a closed system. The substantiation of the correctness and numerical methods for solving various topical applied problems of the Stefan type were summarized in the review paper [6] and the monographs [7][8][9]. The relaxed linearization methods for the solution of the heat transfer problem using finite element approximation were presented in [4]. An efficient enthalpy-based numerical approach with a trust region minimization algorithm was introduced in [10] with various smoothing functions. In [11], the mathematical model with the separate evolution of the phase change interface from the heat equation was presented with XFEMapproximation.
In this work, we consider an enthalpy method with different types of smoothing functions. We propose a smooth coefficient approximation-based on the analytical smoothing of discontinuous coefficients using the error function with a given smoothing interval. We study a grid and time step convergence to examine the accuracy of the method for different lengths of smoothing intervals and compared to the regular linear smoothing. The novel method with smoothing in one spatial interval is proposed. The method provides a minimal length of smoothing that is calculated automatically for the given temperature values on the mesh. We start with the presentation of the methods for the one-dimensional formulation of the heat equation and perform convergence studies against the grid and time-step sizes, where we calculate the errors between analytical and numerical solutions using the proposed smoothing methods. Then, we generalize the methods for finite element approximations and present the results for two-dimensional model problems. Triangular elements are used, but the formulation can be extended for any type of element.
The paper is organized as follows. In Section 2, we describe the problem formulation. In Section 3, we consider the problem in the one-dimensional domain and present a finite difference approximation. In order to construct an approximation, we present three schemes with different smoothing of discontinuous coefficients, which reduces the Stefan problem to a quasilinear problem for a parabolic equation with smooth coefficients. The generalization of the presented schemes is given in Section 4, where we present the finite element approximation for one-and two-dimensional problems. Numerical results for the model problem are presented in Section 5 for one-and two-dimensional formulations, where we investigate the accuracy of the proposed schemes with different types of smoothing. For the one-dimensional problem, we compare the numerical results with the exact solution. A result of the two-dimensional problem is presented for two geometries. The paper ends with the Conclusion.

Stefan Problem
We consider the two-phase Stefan problem in the domain Ω. Let u(x, t) be an unknown function (temperature) and u * be a phase change temperature that occurs on the phase change interface ξ: where Ω 1 is a frozen zone and Ω 2 is a thawed zone: We have conductive heat transfer in Ω 1 and Ω 2 : On the phase change interface ξ, we have a continuous temperature and the jump of the heat flux, which is associated with the latent heat of fusion for water [1,2]: We write the system of Equations (1) and (2) as the two-phase Stefan problem in domain Ω: where δ is the Dirac delta function: and: The coefficients of the equation depend on the function u(x, t) and have a discontinuity of the first kind at u(x, t) = u * .
The main difficulty in the presented problem is associated with dynamic tracking of the phase change interface. Since the δ coefficient of Equation (3) at u(x, t) = u * becomes infinity, a numerical solution of this equation is possible after smoothing the function η [1,2,10]. Let the phase change not occur instantaneously, but take place in a small temperature range (slushy zone). The discontinuous coefficients η of Equation (3) are replaced by bounded continuous (smooth) functions of temperatureη(u): where:c The resulting Equation (4) is a quasilinear parabolic equation with positive, bounded, continuous (smooth) coefficients, which depends on the coefficient smoothing procedure.
We have the following initial condition: and boundary conditions: where Γ 1 ∪ Γ 2 = ∂Ω. Next, we present two novel approaches to the accurate approximation of the Stefan problem and consider a convenient approach with a linear approximation of the η function. First, we define a smooth approximation of coefficients using the error function and linear approximation. Then, we present a novel approach with smoothing in one spatial interval that provides a very accurate approximation of the problem and does not depend on the interval of a slushy zone, which is hard to define because it depends on the speed of phase change interface propagation and can vary in space and time.

Approximation of the One-Dimensional Problem
Let Ω be the one-dimensional computational domain (Ω = [0, L]) and T h be the uniform mesh: where h is the mesh size. For the time approximation, we use an implicit scheme with linearization from the previous time layer, where t m = mτ for m = 0, 1, . . . , M and τ = T/M is the time step size.
We have the following finite difference approximation of the problem (4)-(6) with smooth coefficients on the grid T h : where: for i = 1, . . . , N.
We consider the approximation (7) with the following initial conditions: Let Γ 1 and Γ 2 be the left and right boundaries, x = 0 and x = L. For the Dirichlet boundary condition on x = 0, we have: For the Neumann boundary condition on x = L, we have: for m = 1, 2, . . . , M.

Smoothing on the ∆ Interval
Let [u * − ∆, u * + ∆] be the smoothing interval (slushy zone). For the coefficients of Equations (7) and (8), we have:k Note that, in this approximation, coefficientsc m i ,k m i are defined on the mesh nodes, andk m i+0.5 is defined on the cell as the harmonic average of the values on nodes x i+1 and x i .
In schemes with smoothing on the ∆ interval, the coefficientsη andδ are defined in each node of the mesh and calculated using the given value of function u i on current node x i : We consider the following smoothing of discontinuous coefficients and define smooth functionsη as follows: • Smoothing on the ∆ interval using the error function: • Linear smoothing in the ∆ interval: We note that the accuracy of the presented schemes depends on the choice of the ∆ parameter. If we choose a very large ∆ interval, we overestimate the slushy zone interval and obtain large errors. In general, the parameter ∆ should depend on the spatial and time meshes in approximation as the parameter highly depends on the speed of the freezing front propagation. Other methods of such procedures were proposed in [10], where exponential-or sinusoidal-based functions can be used to smooth functionη on the ∆ interval. In [1,2], the authors presented smoothing using polynomials of a given degree.

Smoothing in One Spatial Interval
In order to construct a scheme without the estimation of the smoothing interval parameter ∆, we present a novel method to an accurate approximation of the two-phase Stefan problem by choosing the minimal possible interval of the slushy zone for the current mesh.
Let a phase change occur on one spatial interval [ In this scheme,η andδ are defined on cell [x i , x i+1 ] and approximated using values of the functions u i and u i+1 . For the coefficients of Equations (7) and (8), we have: where coefficientc m i is defined on each node x i and coefficientsη m i+0.5 andδ m i+0.5 are calculated taking into account the presence or absence of the phase change interface in the temperature interval [u i , u i+1 ].
For the accurate approximation of the two-phase Stefan problem, we present a smoothing in one spatial interval:η At the calculation of the phase change interface position, we use a linear interpolation of the desired function within the indicated temperature interval. Further, we determine the value of the corresponding coefficient for the difference equation.

Generalization of Methods for Finite Element Approximation
In this section, we generalize the presented algorithm for the solution of the two-dimensional problems on unstructured meshes [12]. We use a finite element method and set: For the heat problem (4)-(6), we have the following variational formulation of the problem: find u m ∈ V such that: where: Let T h be the mesh for the domain Ω: where N c is the number of grid cells.
We replace the space V by a finite dimensional subspace V h , V h ⊂ V. We suppose that space V h is composed of piecewise linear basis functions φ i on the grid cell, i = 1, . . . N and N = dim(V h ) is the number of nodes (vertices) for linear basis functions. We have the following decomposition of u h in the basis of V h : Problem (18) is replaced by the following discrete problem: We have the following discrete system in the matrix form for U = (u 1 , ..., u N ) T : where S m and A m are global mass and stiffness matrices: Here, S m,i and A m,i are local matrices in cell K i : In order to obtain a similar form of discrete system as in the previous section with finite difference approximation, we suppose that coefficientk i is defined in each cell (i = 1, . . . , N c ) and coefficientc j is defined on mesh node (j = 1, . . . , N). We use a trapezoidal integration to obtain a diagonal mass matrix.

Smoothing on the ∆ Interval
Similarly to the previous section,η andδ are defined in each node of the mesh and calculated using the given value of function u i on current node x i : with the following approximations: • Smoothing on the ∆ interval using the error function: • Linear smoothing on the ∆ interval: Let n K be the number of vertices in cell K, where n K = 2 for the one-dimensional problem and n K = 3 for the two-dimensional problem with triangular cells. For the stiffness matrix, we have: For the mass matrix with trapezoidal integration, we have: Here, coefficientsc m j ,k m j are defined on the mesh nodes, andk m K i is defined on the cell as the harmonic average of the values on the cell nodes.

Smoothing in One Mesh Cell
Similar to smoothing in one spatial interval for the one-dimensional problem, we construct smoothing in one cell of the grid for the general case, where a cell can be of any shape (interval, triangle, tetrahedron, quadrilateral, etc.).
Let a phase change occur on one mesh cell . . , n K (n K is the number of vertices of mesh cell). We defineη andδ on cell K i using a similar approximation by the values of the function u K i j in cell K i (j = 1, . . . , n K ). For the stiffness matrix, we have: For the mass matrix, we have: where u K i ,m j is the value of function u on node j in cell K i and j, l = 1, . . . , n v . Here,c K i ,m j is defined on mesh vertex j in cell K i using the corresponding value of function u, andk m K i is defined on cell K i . The coefficientη m K i is calculated by taking into account the presence or absence of the phase change interface in cell K i .
For the one-dimensional problem with n K = 2, we have cell Therefore, we have: Note that this approximation is equal to the smoothing presented in the previous section for the one-dimensional formulation using the finite-difference method.
Before the generalization of the smoothing in one mesh cell, we present a brief explanation of the concept that we used in this method. Let V K i + and V K i − be the volume of thawed and frozen zones in cell K i , V K i = |K i | the volume of the cell, and The coefficient η is the unfrozen water concentration that can be represented as a fraction of the thawed zone from the cell volume: . In order to calculate the volume of the thawed zone in mesh cell V K i + , we use a linear approximation of the temperature in cell.
In the one-dimensional case, we have: where: Therefore, the location of the phase change interface in cell K i is the following (see Figure 1): and therefore, we obtain:η for the one-dimensional problem when the phase change occurs in cell In the approximation construction for the two-dimensional case, we use a similar approach, where we find an interface position on each edge between two nodes and use it to calculate a fraction of the unfrozen zone in cell K i . For example, for the case with u 1 > u * , u 2 < u * , and u 3 < u * (see Figure 1), we have:η where the unfrozen zone volume and the volume of the cell are calculated as the area of corresponding triangles. Approximations of theη for other cases can be calculated in a similar way. In Figure 1, we present an illustration of thawed and frozen zones on mesh cell K i . Finally, for the two-dimensional formulation with triangular cells, we have: Here, the values ofδ K i ,2 andδ K i ,3 are calculated in a similar way. The presented method can be extended to other types of elements.
We obtained an approximation of unfrozen water concentration with one mesh cell as a function of cell temperatures on nodes. The proposed methods work with unstructured meshes that can handle complex computational domains.

Numerical Results
We present the results for the proposed difference schemes on one-and two-dimensional examples. We solve the two-phase Stefan problem, which serves as a mathematical model for the pore water freezing/thawing process with input data specified in the international SI system: The phase change temperature is u * = 0 • C. In the calculations, we use a linearized difference scheme, where nonlinear coefficients are calculated using a solution from the previous time layer.
We present numerical results for the three scheme: • erf-smoothing: the scheme with smoothing using the error function on the ∆ interval. • delta-smoothing: the scheme with linear smoothing on the ∆ interval. • h-smoothing: the scheme with smoothing in one spatial interval (mesh cell).
The first and third approximations are based on the definition of a slushy zone interval with parameter ∆.

One-Dimensional Problem
The computational domain size L = 8 m. We simulate for T = 10 7 s with initial condition: The boundary conditions are the following: where we consider two test cases with g = −5 • C and g = −15 • C. For a one-dimensional problem, an exact solution of the two-phase Stefan problem is defined using the phase change interface location ξ = ξ(t): where γ is the coefficient of proportionality (m/s 0.5 ).
The exact solution has the following form [13]: Using a phase change interface condition (2), we find γ as the root of the transcendental equation: where a 1 = √ k 1 /c 1 , a 2 = √ k 2 /c 2 , and erf is the error function [13][14][15][16]. The value of γ is calculated using the iterative secant method with the accuracy of < 10 −12 . We have γ = 0.00023897230346 for g = −5 and γ = 0.0004188066281859222 for g = −15.
In Figure 2, we present exact and numerical solutions at the final time T = 10 7 s ≈ 115 days and a phase change interface location ξ(t). We present the results of the numerical solution of the problem using three types of smoothing. Numerical calculations were performed with n = 200 (mesh by space) and m = 200 for the time approximation. We consider the case with g = −5 • C and g = −15 • C. The smoothing parameter is chosen constant ∆ = 0.3 throughout the entire computation time t ∈ (0, T].
In the left figure, a solution at the final time is presented, where we see that the first derivative of the function has a discontinuity of the first kind with respect to spatial variable x at u = u * , which agrees with the Stefan condition. In the right figure, we depict a phase change interface location by time.
In numerical calculations, the phase change interface is found through the linear interpolation method in the interval, where the sought function has different signs at the ends (u m i · u m i+1 ) < 0. We observe that the numerical solution both for the temperature at the final time and for the phase change interface location is accurate for smoothing using error functions (erf-smoothing) and for smoothing in one spatial interval (h-smoothing). We note that the novel method based on smoothing in one spatial interval is more accurate than the scheme with approximation using the ∆ parameter, that is, in general, it should be carefully chosen for each test case and depends on the mesh size, number of time steps, and speed of the phase change interface propagation.   Next, we investigate the influence of the mesh size and the number of time steps on the accuracy of the presented methods. We consider n = 25, 50, 100, 200 for the spatial mesh and m = 50, 100, 200 for the time approximation. We perform simulations and calculate errors for two smoothing parameters ∆ = 0.3 and ∆ = 0.6 in the erf-smoothing and delta-smoothing methods. We calculate the relative L 2 errors for a solution and phase change interface: where y, ξ are the numerical solutions, u(x, t), ξ e (t) are the exact solutions, andũ, ξ f are the solutions on the very fine grid with n = m = 800.
In Tables 1-3, we present the relative errors for the three methods: erf-smoothing, h-smoothing, and linear delta-smoothing. We present relative errors in % and investigate the influence of the boundary condition, smoothing interval ∆, and the space and time mesh parameters on the accuracy of the numerical solution. In Table 1, the relative errors for the scheme with smoothing using an error function are presented (erf-smoothing) for ∆ = 0.3 and ∆ = 0.6, where errors for g = −5 are presented on the left and for g = −15 presented on the right. In Table 2, the results for the scheme with smoothing in one spatial interval (h-smoothing) are presented. Relative errors for the scheme with linear smoothing on the ∆ interval are presented in Table 3 (h-smoothing) for ∆ = 0.3 and ∆ = 0.6. From the errors between the exact and numerical solutions (e e (u T ) and e e (γ)), we observe that the erf-smoothing and linear delta-smoothing methods are highly dependent on parameter ∆. The error decreases along with the spatial mesh refining as the number of time layers increases. Table 1. One-dimensional problem with erf-smoothing. Relative L 2 error for the temperature at the final time and phase change interface location with respect to the exact solution and a very fine grid solution.   The errors between a very fine grid solution (m = n = 800) and a numerical solution for different parameters m and n (e f (u T ) and e f (γ)) demonstrate the convergence of the solution depending on m and n. We observe that for the coarse parameters m and n, we can obtain results with a large error. For example, we have e f (u T ) = 25% and e f (γ) = 70% for ∆ = 0.3, g = −15, and m = n = 50 for the erf-smoothing method. For larger ∆, we can obtain more stable results. However, for a finer mesh and more time steps, the parameter ∆ = 0.6 is too large to obtain good approximation. For example, we have e f (u T ) = 3% for ∆ = 0.6 and e f (u T ) = 1% for ∆ = 0.3 in the test problem with g = −5 and m = n = 200. A new method with smoothing in one spatial interval (h-smoothing) does not have parameter ∆ and provides very good results with small errors and convergence depending on parameters m and n. This method automatically determines a sufficient smoothing parameter.
A graphical illustration of the errors from Tables 1-3 for the presented methods is shown in Figures 3 and 4. We depict errors for m = 100 and m = 200 with different spatial mesh parameters n and ∆. The h-smoothing scheme provides accurate results on a coarse spatial grid at large time steps. We observe that the erf-smoothing method is less sensitive to the parameter ∆ and provides sufficiently good results with a good choice of the parameter ∆. By comparing the results for g = −5 and g = −15 in Figures 3 and 4, we see that for a larger temperature gradient, we obtain worse results. In the erf-smoothing and delta-smoothing methods, we should use a finer mesh and more time steps for better approximation with a good choice of parameter ∆. The scheme with smoothing in one spatial interval provides an accurate numerical solution with the minimal smoothing interval, which is different in each spatial interval and time layer.

Two-Dimensional Problem
In this section, we present numerical results for a two-dimensional heat problem. We consider the two-phase Stefan problem with similar parameters as in the previous section. We consider computational domain Ω = [L 1 × L 2 ] with size L 1 = L 2 = 2 m and simulate for T = 10 6 s with initial condition u 0 = 5 • C. As boundary conditions, we set g = −5 • C and g = −15 • C for x ∈ Γ 1 = {x ∈ ∂Ω : x = 0, y = 0} and zero flux on other boundaries, ∂Ω/Γ 1 . The computational mesh is uniform with triangular cells with size n = n 1 × n 2 . Numerical calculations are performed with m = 200 for the time approximation. We present the results of the numerical solution of the problem using three types of smoothing schemes for n = n 1 × n 2 , n 1 = n 2 = 25, 50, 100, and 100.
We calculate relative L 2 errors for the solution at the final time: where y is the numerical solution andũ is the solution on the fine grid with n = m = 400 using the current method. Here, for u f ,hs (x, T), we use a reference solution with the h-smoothing method. We compare the results with h-smoothing to show that this method provides an accurate solution and does not depend on the parameter ∆ as in the previous examples for the one-dimensional problem.
In Figure 5   In Tables 4-6, we investigate the influence of the mesh size (n) on the accuracy of the presented methods. Time step number m is fixed, m = 200. We consider n = 25, 50, 100, 200 for the spatial mesh. In Table 4, we present the relative errors for the erf-smoothing scheme with ∆ = 0.15, 0.3, and 0.6. In the left table, we present the errors for boundary condition g = −5. The table at the right demonstrates the results for g = −15. We present relative errors in % and investigate the influence of mesh size on the accuracy of the numerical solution. In Table 5, the results for the h-smoothing scheme are presented for g = −5 and −15. Relative errors for delta-smoothing are shown in Table 6. The first error e e (u T ) is an error between the current solution and a reference solution using the h-smoothing scheme on mesh m = n = 400. We observe that in order to obtain good results using the erf-smoothing scheme, we should take smaller ∆ for g − 5 than for the test case with g = −15. For example, we have 6.5% of the error for ∆ = 0.6 and 1.2% of the error for ∆ = 0.15 in the test case with g = −5. Table 6 shows the optimal value of ∆ = 0.3 for g = −5 and ∆ = 0.6 for g = −15. For the h-smoothing scheme with smoothing in one spatial interval, we obtain good results for both test cases with g = −5 and g = −15.
The second error e f (u T ) illustrates the convergence of the solution with increasing mesh parameters m and n to the reference solution on the fine grid m = n = 400 using the corresponding method. We observe that both errors e f (u T ) and e e (u T ) are sensitive to the size of the mesh, where we obtain a large error for a very coarse grid n 1 = n 2 = 25. We obtain accurate results for n = 50, 100, and 100 using the h-smoothing method. Table 4. Two-dimensional problem with erf-smoothing. Relative L 2 error for the temperature at the final time with respect to the finest grid solution using the presented method (e f (u T )) and h-smoothing (e e (u T )).  Table 5. Two-dimensional problem with h-smoothing. Relative L 2 error for the temperature at the final time with respect to the finest grid solution using the presented method (e f (u T ) = e e (u T )).  Table 6. Two-dimensional problem with delta-smoothing. Relative L 2 error for the temperature at the final time with respect to the finest grid solution using the presented method (e f (u T )) and h-smoothing (e e (u T )).

Geometry with an Unstructured Grid
The main advantage of the finite element method is that it provides the possibility of problem solving in complex geometries using the construction of the unstructured mesh with triangular elements.
We consider the computational domain with size L 1 = 6 m with a rough boundary and two circle perforations. In Figure 6, we depict the computational domain and unstructured mesh. In the left figure, we present the computational domain. The unstructured grid with triangular elements is presented in the right figure. We depict a rough top boundary Γ 1 in green color, where we set the Dirichlet boundary conditions: u(x, t) = g, x ∈ Γ 1 . Similar to the previous results, we consider g = −5 and g = −15. On the circle perforation, we set the Robin-type boundary conditions: where we set α = 20. In Figure 6, the first circle is indicated in red color with g 1 = 10, and the second circle is indicated in orange color with g 2 = 15. Circle perforations illustrate a pipe influence to the temperature distribution. The pipes contain fluid with temperature g 1 (first) and g 2 (second) with heat transfer parameter α. The radius of the first circle is r 1 = 0.1, and the radius of the second one is r 2 = 0.2. On other boundaries (blue color), we set the Neumann boundary condition: We simulate for T = 10 6 s with initial condition u 0 = 2 • C. A computational mesh is fixed and contains n vertices, n = 18,242.
The numerical solution at the final time for g = −5 • C and g = −15 • C is presented in Figures 7 and 8 We calculate similar errors for the solution at the final time (e e (u T ) and e f (u T )). The error e f (u T ) is calculated using the solution with m = 400 using the current method. The error e f (u T ) is calculated using the reference solution with the h-smoothing method. In Tables 7-9, we investigate the influence of the number of time steps m and ∆ on the accuracy of the presented methods. We consider m = 50, 100, 200 and ∆ = 0.15, 0.3, 0.6.
In Table 7, we present relative errors for the erf-smoothing scheme with ∆ = 0.15, 0.3, and 0.6. In the left table, we present errors for boundary condition g = −5. In the right table, results for g = −15 are shown. In Table 8, the results for the h-smoothing scheme are presented for g = −5 and −15. The relative errors for delta-smoothing are shown in Table 9 for different values of m and ∆. We observe good results for the erf-smoothing scheme, when we take ∆ = 0.15 for g = −5 and ∆ = 0.3 for g = −15. In the delta-smoothing scheme, we have good results for ∆ = 0.3 for g = −5 and g = −15. We obtain an accurate solution for the scheme with h-smoothing for m = 50, 100, and 200 for g = −5 with less then one percent of error. The scheme with h-smoothing works very well for g = −15 with less than one percent for m = 200.    We present and investigate three schemes: smoothing on the ∆ interval using the error function (erf-smoothing), linear smoothing on the ∆ interval (delta-smoothing), and smoothing in one mesh cell (h-smoothing). The first algorithm (erf-smoothing) is based on the analytical spreading of the point source function and smoothing of discontinuous coefficients. The algorithm with smoothing in one mesh cell (h-smoothing) provides a minimal length of the smoothing interval that is calculated automatically for the given temperature values on the mesh at a given time layer. The third scheme (delta-smoothing) is a convenient scheme based on smoothing using linear approximation on the ∆ interval. We study the proposed schemes on one-dimensional and two-dimensional model problems. We consider a problem with two values of the boundary conditions on different spatial and time meshes. We present errors for the methods with smoothing on the ∆ interval for different values of the smoothing interval parameter ∆, which is fixed and given as a constant value. Numerical results show that the smoothing scheme with the minimal smoothing interval shows its computational efficiency and provides good results.

Conclusions
In this work, we consider the two-phase Stefan problem, which describes the heat transfer problems with phase change. In a two-phase formulation, the mathematical model is defined in the whole domain, which contains frozen and thawed zones.
The main problem is associated with the sharp phase change interface. For the numerical solution of the problem, we present and investigate schemes with smoothing in one spatial interval (cell) and smoothing on the ∆ interval using the linear function and the error function. We observe that smoothing with the error function provides better results with smaller errors than linear smoothing on the ∆ interval. However, the novel smoothing on one spatial interval (cell) gives better results, does not depend on the parameter ∆, and can automatically smooth the coefficient in one mesh cell. We generalize the presented methods for the two-dimensional formulation of the problem and present results for two geometries with structured and unstructured meshes.