A discontinuous Galerkin method for mathematical modeling of ice melting at the interaction with the environment

A Stefan problem with a moving boundary in a heterogeneous two-phase medium is considered. For solving the Stefan problem, a discontinuous Galerkin method in a nonsymmetrical Inner Penalty formulation is used. We propose to split the problem solution into a discontinuous component in the phase transition zone and a continuous component outside the phase transition zone. A multiscale method allows to reduce a number of freedom degrees by using the discontinuous Galerkin method in the phase transition zone only. Mathematical modeling results of the ice melting at the interaction with the environment are given.


Introduction
A study of the ice cover and permafrost layer state is an important scientific problem, since it is aimed at preventing man-caused and environmental disasters. According to [1,2], two main factors of negative influence on the ice cover and permafrost layer are possible. The first factor is associated with natural phenomena, which can be seasonal in nature [1]. For example, anomalously warm and arid summer. The second factor is due to technogenic impacts [2] in the human activity process. From the physical point of view, both types of influence initialize the heat transfer process with different types of interface conditions and moving boundaries.
A solution of problems with a moving boundary is required in the study of physical processes with a phase transition or mixing of substances. Stefan's problem of the "solid-liquid" phase transition is a classic example of a problem with the moving boundary. This problem is of great practical importance in predicting a permafrost layer state and creating phase-change composite materials [1 -3].
The solid surface interacts with a medium whose temperature is higher than the phase transition temperature. As a result, the boundary between the liquid and solid is moved toward the liquid phase.
A modeling procedure of the heat transfer process in media with the phase transition is complicated by the search for the current position of the moving boundary [3]. Therefore, a computational scheme must be adaptable to the specifics of the problem and maximally independent of computational domain meshes.
To solve such a class of problems, computational schemes of finite differences methods (FDM), finite volume methods (FVM) and finite element methods (FEM) are used usually. FDM and FVM are successfully applied to solve one-dimensional multiphase Stefan problems, since they allow the use of such popular methods for tracking the phase transition front as catching the front node [4], enthalpy and phase field methods [5]. However, their application for solving the three-dimensional problems is complicated, since the FDM is sensitive to the orientation of the mesh [4], and the FVM requires the construction of special multilevel interpolation schemes [5]. 2 Modern computational schemes of FEM are implemented using explicit methods for tracking the phase transition front with fixed or moving meshes. An extended FEM (XFEM) uses the fixed meshes. XFEM is based on the introduction of additional freedom degrees in the front-moving zone [6]. The effectiveness of XFEM decreases when solving the multiphase Stefan problems, since the account of the phase transition boundaries becomes nontrivial at the functional level.
We propose to use a multiscale Discontinuous Galerkin Method (MDG) with a moving nonconformal mesh. In the MDG-method, the solution is represented as the sum of the discontinuous and continuous components. The discontinuous component exists in the phase transition zone only. The continuous component is defined in the entire computational domain. The strategy allows to construct a scheme of the discontinuous Galerkin method with the computational structure of the classical FEM [7,8].
The multiscale discontinuous Galerkin method in the nonsymmetrical variational "Inner Penalty" formulation for the two-phase three-dimensional Stefan problem using nonconformal tetrahedral moving meshes are used. Results of solving the problem of melting ice with constant heat exchange with the surrounding medium are presented. Figure 1 shows a computational domain. The computational domain consists of two phases: ice and water regions. Let Ω l is the liquid phase (water), and Ω s is the solid phase (ice). At the initial time, the medium is completely filled with ice. On the upper and lower boundaries S 1 the heat exchange occurs with the medium with the temperature of 1 °C and -5 °C, respectively. The side faces of S 2 are assumed to be thermally insulated. A front of the phase transition (melting) is denoted as ξ(t). At the initial instant of time, front is located in the plane Z = 0. The phase transition temperature is 0° C. The length of the computational domain is 25 m, width is 16 m, and height is 3 m. It is necessary to determine the behavior of the "ice-water" phase system over time.

Mathematical Model of the Stefan Problem
The temperature distribution in the each phase is described by the heat conductivity equations [3] ρ with initial conditions ρ -density, kg/m 3 ; c -specific heat, J/(kg·˚C); λ -thermal conductivity, W/(m·˚C); Т -temperature, ˚С; l -liquid phase, s -solid phase.
At the side faces S 2 at the upper and lower boundaries, the process of heat exchange with the medium is given as where β = 10 W/(m 2 ·˚С) -heat transfer coefficient.
On the boundary of the phase transition, ideal contact conditions for the temperature distribution are determined as and conjugation conditions for the heat flux are satisfied where L -phase transition heat, J / kg.

Variational Formulation of the Discontinuous Galerkin Method for the Stefan Problem
Let be a partition into disjoint sets Ω k . Let Г : k   be the set of the boundaries, be the set of the inner boundaries, and be the trace space. We introduce finite-dimensional subspaces on the set as , , where is a polynomial space of degree m.
where the index indicates belonging to Let : C D      and are defined as [8] : , where Г Г [ ] 0, { } : The nonsymmetrical formulation of the discontinuous Galerkin method for the problem (1) -(7) [8]: to find such where t  -discretization time step.

Experimental Results
To solve the problem, we choose adaptive nonconformal tetrahedral meshes. Figure 2 shows the 3D triangulation of the computational domain at the time moment t = 100 h. The combined solver BiCG + GMRES (5) is used to solve SLAE. BiCG is the oscillating method and allows you to quickly approach an exact solution. The GMRES method is a smoothing agent and ensures uniform convergence to the solution. Physical parameters are presented in the Table 1.  Figure 3 shows the dynamics of the position of the phase transition front over time depending on the time discretization step (from 10 min to 40 min).

Discussion of Results
Analyzing the graphs in Figure 3 and Figure 4 allows us to conclude that the results of modeling the process of ice melting when using the proposed computational scheme is correct. At the time t = 3000 min, the position of the phase transition front corresponds to the coordinate z = -0.052 m (see Figure  3), which is confirmed by the graphs in Figure 4. At the point z = -0.052 m, a "bend" of the curve is observed (see Figure 4), since at this point the physical properties of the medium is changed (the boundary of the phase transition of the ice into water).

Conclusion
Ice melting at the Interaction with the Environment is described by the two-phase Stefan problem. For solving the two-phase Stefan problem, we used the computational scheme of the multiscale discontinuous Galerkin method using the nonconformal tetrahedral meshes. The Stefan condition is taken into account in the weak form in the IP (Inner Penalty) variational formulation of the DGmethod. The mathematical modeling results of ice melting are obtained at the constant heat exchange with the media with the positive temperature depending on the time discretization step.