Numerical solution of two-dimensional (2D) nonlinear heat conductivity problem on moving grids

For the numerical solution of two-dimensional (2D) nonlinear heat conduction problem a transition to a moving coordinate system is used. A local speed of the latter is chosen in such a way that the considered process is close to stationary. Such coordinate system guarantees a concentration of grid nodes in a range of the solution features. This allows for a quality improvement of the obtained numerical solution with a small number of computational grid nodes. The developed numerical algorithm is tested on the heat wave propagation problem which has an exact solution.


Introduction
It is known that the quality of numerical solution often strongly depends on the choice of the computational grid used in nonlinear problems. Presence or occurrence of regions where the solution is either with high gradients, or changes rapidly, or loses smoothness or even continuity -demands special measures in order to keep the result on acceptable level. On the other hand, there are a lot of problems where mentioned peculiarities are of great interest, for example, fronts of shock-waves and boundary layers in gas dynamics [1] or phase-change fronts in relevant applications. The latter should be carefully considered while dealing with cryosurgery operations modeling because of a high risk for the patient's health [2].
The case of elliptic equations, where the optimal grid is built before the computation, seems to be the most studied [3,4]. In general case, the position and the moment of peculiarity occurrence is not known, thus the development of a universal algorithm for grid adaptation is practically important. It should be noted that such approach must not have any side parameters fitting it to the current problem.
One of the most promising and elegant approaches is associated with the transition to the moving coordinate system where the considered processes is close to stationary -quasistationarity principle [5]. Such approach automatically thickens the computational grid for the peculiarities without introducing any side parameters. In this paper we propose and study in detail a multidimensional generalization of this approach.

The transition to the moving coordinate system
In this paper we consider the following model of the initial-boundary problem for the heat conduction equation: The considered approach is based on the moving coordinate system introduction: We introduce the following designations for convenience: Here we will call Ψ α β a local deformation matrix, while Q α is associated with the moving coordinate system local speed. The differentiation operators of (1) can be rewritten in terms of new basis ∂ ∂τ , ∂ ∂y β as following: Thus, if the following compatibility conditions are satisfied, the governing equation in (1) equals to the following: In this way, the introduction of the moving coordinate system transforms the governing equation (1) to the pair of (6), which describes the considered physical process in new coordinates, and (5), responsible for the evolution of the coordinate lines geometry. It should me mentioned, that due to (5) the equation (6) can be reintroduced in conservative form:

Quasi-stationarity principle
According to the principle of quasi-stationarity the coordinate system where the considered process is close to the stationary should be chosen. We propose the following multidimensional modification of this approach. Let us consider (6) where ∂u ∂τ = 0: If we septate the convection and diffusion parts and neglect the diffusion terms proportional to ∂ 2 u ∂y γ ∂y β , the (8) will become D β ∂u ∂y β = 0, and it should be satisfied for every solution u. Then D β = 0. Therefore, the optimal coordinate system's speed is:

Numerical approximation of the governing equation and compatibility conditions
For the numerical solution of the considered problem the uniform orthogonal grid on I 2 in yspace was introduced. The functions u, Q α was calculated in the nodes of the grid, Ψ α β -on the faces. During the numerical experiments the governing equation in the conservative form (7) was numerically integrated via an explicit finite difference scheme.
Due to "Ψ-Q" conditions Thus, if "Ψ-Ψ" conditions are satisfied simultaneously, they will be satisfied at every moment of time. In this way it is preferable to use the non-deformed coordinate grid x α (y β , 0) = y α at the initial moment, when Ψ α β ≡ δ α β and "Ψ-Ψ" conditions are satisfied wittingly. It is important to use such approximation of compatibility conditions that the grid analogue of this fact would be saved.
According to the definition of the functions Q α and Ψ α The transitions between the neighbouring nodes on the current layer and the neighbouring states of the current node in the physical x-space must be commutative. It leads to the approximation of compatibility conditions in the following form: for "Ψ-Ψ" condition and  for "Ψ-Q" condition. One can directly check the following Proposition. Due to (14), the grid analogue of "Ψ-Ψ" condition in form of (13) does not depend on time.
This circumstance allowed for the consideration of only "Ψ-Q" compatibility condition during the numerical integration of the problem.

Numerical results
A program realization of the proposed numerical algorithm was tested on the "heat-wave" problem. It is known that the equation allows for an exact solution in travelling waves [6]: Thus, via rotation of the coordinate system one can obtain a good two-dimensional (2D) test problem with non-smooth solution. The following parameter values were used: p = 2, D = 0.25. The results of numerical experiments for two different moments in time are shown Figure 1 and Figure 2. On the left of both results the numerical grid in x-space is shown. On the right of both results the solution profile with error isolines are presented.

Conclusion
A multidimensional generalization of the quasi-stationarity principle for the dynamic adaptation of the computational grids for the numerical solution of some parabolic problems was proposed.
The numerical algorithm was tested on the heat wave propagation problem, which has an exact solution.