You can read the original article at the World Congress Engineering web site.

* Abstract* — The quasilinear heat equation with thermal conductivity and heat capacity depending on the temperature field in three spatial dimensions is studied in application to the phase transition problem in permafrost ground. The conditions under which the alternating directions Douglas – Rachford finite difference scheme retains numerical stability are explicitly formulated. The comparison with the known analytical similarity solution to the Stefan problem in one spatial dimension is performed.

* Index Terms*—quasilinear heat equation, Stefan problem, finite differences, alternating directions scheme, numerical stability.

I. INTRODUCTION

Since the early formulations of the alternating directions implicit methods (ADI methods) [1], [2], they have been tremendously developed and found a vast number of applications [3], [4]. Nevertheless, serious difficulties are encountered with the use of these methods in application to problems with complex geometries [5] and/or nonlinear equations of mathematical physics [6].

Although the schemes of the ADI methods are proved to be efficient and economic with respect to time consumption and, in most cases, unconditionally stable, they exhibit some disadvantages:

1) Their finite-difference formulations permit to consider only rectangular spatial domains (due to commutativity requirements imposed on the factorized and split operators) [7];

2) The application of the ADI schemes to the problems with Neumann and Robin boundary conditions that are varying in time encounters serious problems due to the necessity of evaluation or approximation of these boundary conditions at the intermediate steps of the scheme [8];

3) When applied to the solution of nonlinear heat equations, the operators constituting an ADI scheme do not commute, thus leading to the loss of unconditional stability of the scheme [6].

The first of the above disadvantages can be overcome either by the use of finite elements methods in conjunction with operator splitting techniques, or by domain decomposition techniques. The second and the third disadvantages present an important problem for the successful application of the ADI scheme. To the best of our knowledge, no complete stability analysis for an ADI scheme applied to the nonlinear heat equation in a three-dimensional spatial domain is available in the literature, thus motivating this work.

Another motivation for the present work is the application of ADI scheme to the modeling of heat transfer in large scale environmental systems (e.g., large areas of permafrost ground) which, in the case of purely explicit finite-differences schemes, imposes stiff constraints on the time-step in order to guarantee the numerical stability. At the same time, implementation of implicit schemes can often lead to much greater computing expenses than that of explicit schemes, especially for the problems with rapidly changing coefficients in complex geometries and substantially nonhomogeneous meshes. Thus, in modeling of heat transfer in large scale systems the necessity of making an optimal choice between explicit and implicit schemes arises. In case of finite-elements method, applied to modeling of processes in permafrost ground, the analysis of numerical stability appears to be so complex, that the stability criterion is often established empirically [9].

In the present paper we discuss the application of the ADI Douglas – Rachford scheme to the solution of Stefan problem in porous permafrost ground. The paper is organized as follows: next section contains the problem formulation and some assumptions that will be used in the proof of numerical stability of the ADI scheme while section 3 exposes the proof itself. Section 4 presents some numerical results and is followed by Conclusions.

II. PROBLEM SETUP AND FORMALISM

*A. Quasilinear heat equation*

The equation describing heat transfer in a system with phase transition has the following form:

where _{} is the unknown function (temperature:_{} at every fixed _{} being a real linear space of positive valued functions), _{} is the heat capacity, _{} is thermal conductivity. By _{} we denote a bounded, open subset of an Euclidian space _{} with boundary _{}, the closure of _{} being denoted by _{}. Following the classification given in [10], we call Eq.(1) a quasilinear heat equation. The initial and boundary conditions are taken to be:

We consider the cases when _{} and _{}, thus Eq. (1) is uniformly parabolic and has a unique solution to the initial-boundary-value problem Eq. (3)–(2).[11].

*B. Stefan problem in permafrost ground*

Basing on results of [12], we formulate the model for the Stefan problem without explicitly invoking the front-tracking condition. This approach is justified by the observations that in case of explicit front-tracking models applied to Stefan problem in salted permafrost ground, there appears an overcooled region (frozen fringe zone) [13] (the analysis of this zone and related frost heave and cryogenic suction processes [9] are beyond the scope of the present work). Thus, we take into account the phase transition by introducing effective heat capacity _{} which incorporates the latent heat per unit mass _{}

where _{} and _{} are the values of volumetric heat capacities of the ground in thawed and frozen phases respectively, _{} is the fraction of frozen water, _{} is the smoothing parameter, _{} is the phase transition temperature, _{} is the density of water. The thermal conductivity is taken to be

where _{} and _{} are the values of thermal conductivities in thawed and frozen phases respectively.

*C. The linear space _{}*

For the formulation of finite difference scheme we introduce the following discretization procedure. Let the vector _{} have positive coordinates and _{} be the set of all points _{}, the indices _{} being integers. Two points _{} and _{} belonging to _{} are called neighbors if _{}. The points _{} all of whose neighbor belong to _{} are denoted by _{}. Following [14], we denote the points _{} with the property that at least one neighbor belongs to _{} by _{}. Thus, the full spatial mesh is _{}, _{} being the set of boundary points (outside of _{}). For the following we assume that _{} is a homogeneous (i.e., _{} are constants) cubic domain.

Let _{} and _{}. Then the time mesh is defined as follows:

The approximate solutions of Eq.(1), _{}, are defined on _{} and take their values in a real finite-dimensional linear space _{}, the dimension of which is equal to the number of points in _{}. A function _{} is called admissible if _{} and _{} on _{}. As was pointed out in [14], _{} is, in general, not known in _{} and thus the following assumption should be made: there exists a null sequence (that converges to zero) _{} of mesh spacings such that _{} and _{} always belongs to the sequence _{}. This assumption guarantees that an admissible function is uniquely specified in _{} in terms of the initial and boundary values of _{}.

For every _{} we define the first-order forward and backward difference operators:

where _{} being the number of points in _{} along _{}-direction in the Cartesian coordinates. For any suitably defined function _{} we set

Thus, _{} are the values of _{} at the fictitious intermediate nodes of the mesh _{} at _{}.

From the Eqs. (8)–(9) one obtains:

where

and _{} corresponds the value of _{} at _{} and _{} being defined as _{}. We also provide _{} with the _{}-inner product _{} and induced norm _{} on _{}:

The maximum norm is defined as following:

We define on _{} a modified _{}-inner product _{} and induced norm _{} by

where _{} and _{} correspond the values of _{} and _{} at _{} respectively. Note, that the so defined operator _{} implies that it is generally not true that _{} (see a discussion in the next subsection).

*D. Alternating direction finite difference scheme*

Let us now put in correspondence to Eq. (1) the following finite differences scheme:

where _{} are given by Eq. (15), _{}, _{} and _{} at _{}. Taking into account Eq. (8)–(11) and (15), Eq. (16) can be rewritten in the following form:

where _{} (same for _{}), _{} and _{} are the auxiliary functions that are used at the intermediate steps.Note that the indices _{} in the first equation of system (17) change in the following ranges _{}, _{}, while for the last two equations in (17) _{}, _{}, _{} and _{}, _{}, _{} respectively.

The scheme (17) has been described in [1] and is called the Douglas-Rachford scheme. It has been pointed out by Douglas in [1] that the scheme (17) with _{} defined by Eq. (11) is locally second-order correct in space, unless the boundary conditions lower the accuracy. From Eq. (16) one can see that Eq. (1) is approximated with the accuracy _{} in _{}. The boundary and initial conditions (3)–(2) are formulated on the mesh _{} in the following way:

where .

It is well known that the scheme (16) is unconditionally stable when _{} and _{} in Eq. (15) are constants (i.e.,_{} and _{}) and _{} are the finite-difference forms of the Laplace operators. This is readily seen from Eq.(16) and Eq. (14) if noticed that _{} are selfadjoint positive definite operators, _{}. In this case one can use the following Theorem proved in [15]:

**Theorem 1:** If operators _{} and _{} in Eq.(16) are self-adjoint, do not change in time, and _{} is positive definite, then the condition

is necessary and sufficient for the following estimate to hold:

where _{} is the norm induced by the _{}-inner product (14) in the case of constant _{} and _{}, _{} is the initial condition, and _{} is the numerical solution obtained at the _{}-th iteration.

Indeed, if _{ } and _{ } are constants, the operators _{ } defined by Eq. (15) are pairwise commutative and thus operator relation (19) holds for any small _{ }. Relation (20) means that the scheme (16) is numerically stable when _{} and _{} are constants.

The introduction of _{} and _{} that depend on _{} makes the operators _{} defined by (15) nonlinear, i.e. it is generally not true that _{}. This fact makes it impossible to prove numerical stability in a single norm, as it is done in Eq.(20). In order to proceed, we slightly modify the definition of operator _{} and make the following assumptions:

# Assumption 1:

In the ADI scheme (16) the operators _{} are acting on the vectors _{} in the following way:

This assumption leads to consideration of a linear finite difference scheme of type (16) for the heat equation with variable coefficients _{} and _{}, the forms of which depend on the temperature field _{} calculated at the previous time level. This kind of assumption is commonly used when considering nonlinear time dependent problems (e.g., [16]).

Let us rewrite (21) in the following equivalent form:

where

# Assumption 2:

The functions _{}, _{} and _{} (defined by (4)–(6), (9) and (24)) are mapping the elements _{} into a sequence space _{}, the elements of which _{} satisfy the following relations:

where the operators _{ } are defined similarly to (8); _{ } are real positive constants. _{ } is the Banach space with the maximum norm.

III. NUMERICAL STABILITY OF ADI SCHEME APPLIED TO THE NONLINEAR HEAT EQUATION

In this section we make use of Assumptions 1 and 2 and obtain conditions under which the scheme (16) with Eq.(21) can be successfully applied to the solution of nonlinear heat equation (1). The main problem one has to consider is the following: obtain some relation between energy norms of operators _{ } that could guarantee numerical stability of (16). In order to proceed, let us accept the following definition of numerical stability given in [15]:

# Definition 1:

Let _{ } be a positive definite self-adjoint operator acting in a (real) Hilbert space *H*. Let _{} be a constant, such that _{ } for any _{ } from (7), _{} being a constant that does not depend on _{}. The scheme (16) is said to be stable in a Hilbert space _{} (with scalar product and norm defined by (14)) if for any _{} and for any sufficiently small _{} and _{}, for the solution of equation _{} the following relation holds:

where _{} is the identity operator.

In order to prove that the scheme (16) is stable in the above sense, we will need some auxiliary results proved in [15]:

# Lemma 1:

Let _{} be a Hilbert space (real or complex). Let _{} be a linear operator on _{}, such that _{}. Then the relation

(with _{}) is sufficient for

to hold.

# Lemma 2:

Let _{} be a time-dependent (for _{} from (7)), positive definite and self-adjoint operator. Then, relation

is equivalent to relation

where _{} is a constant that does not depend on _{}.

Now, based on the results of [15], we can formulate the main theorem:

# Theorem 2:

Let the operators _{} (defined by (21))satisfy the following condition for any _{}

where _{} is a constant. Then the relation

with _{} from (7), _{}) is sufficient for the following estimate to hold:

where _{} are the solutions of Eq.(16).

# Proof:

Operators _{} (defined by (21)) are positive definite and self-adjoint, thus permitting to define _{} and _{}. Let us rewrite Eq.(16) in the following form:

where

Now, using Lemma 1 and (34) we obtain:

Thus, for _{} and _{} one has:

Taking into account Lemma 2 and the first equation of (38), we obtain:

Hence,

And, finally, returning to the scalar products, one obtains (35).

Let us now examine the conditions under which relation (34) holds. One can see that the l.h.s of (34) contains the terms _{} and _{} (in addition to _{}) which are required for the inequality to hold. Taking into account Assumption 2, a straightforward calculation of a quantity _{} gives us:

where

the max() and min() being taken over all _{}.

From (44)–(48) it follows that for the inequality (34) to hold, the sufficient condition is:

or, more rigorously:

From (50), the sufficient stability criterion (that relates _{} and _{} and guarantees that (34) holds on every time level _{}) can be directly expressed as follows:

where _{} is a constant that can be deduced from (50). Note that (51) resembles the textbook expression for the stability criterion of a finite-differences scheme for a onedimensional heat equation with variable coefficients:_{}.

IV. NUMERICAL EXPERIMENTS

In this section we demonstrate the results of application of the ADI scheme (16) to the solution of a Stefan problem in a rectangular parallelepiped. We compare the obtained numerical results with the known analytical expressions of the similarity solutions for infinite domains [17]. Consider the following problem: a one-dimensional semiinfinite domain _{} contains the material in its liquid state in _{} and in its solid state in _{}, _{} being the point that separates the liquid and the solid phases. The temperature at the point _{} at time moment _{} is denoted by _{}. The phase transition temperature is assumed to be _{} degrees Celsius. We also assume zero heat fluxes at the ends of _{} and a piecewise constant initial condition:

The formulated problem is known to have the following solution [17]:

where _{} and the constant _{} is obtained by solving the equation:

where erfc() is the complementary error function, _{} is the latent heat of phase transition, _{} and _{} are the thermal diffusivities in the solid and the liquid phases respectively which are assumed to be constant in time and position.

In our numerical experiment we have considered the spatial domain with dimensions: _{}, _{} meters. The spatial mesh _{} has the following nodes numbers: 4х4х172 nodes. On all the faces of the parallelepiped the Neumann boundary conditions with zero heat fluxes have been imposed. The initial condition was taken to be (52) with _{} meters. Thermal conductivities and heat capacities of thawed and frozen phases that appear in (4)–(6) are taken to be: _{} _{}, _{}, _{}. Numerical stability for the problem with these parameters has been achieved when the time step _{} and spacial step _{} were subject to the following condition: _{}. This gave us the number of iterations: 6486.

Fig. 1. The dependence of temperature on the Z-coordinate: comparison of numerical results (blue stars and green triangles), obtained with the scheme (16), and analytical expression of the similarity solution (53) (red circles). Time moment 28 days. Inset shows the same curves at temperatures below Tph = 0 Celsius. Red vertical line inside the parallelepiped indicates the set of points, for which the results are presented.

Fig. 1 shows the numerical results along the line _{} meters (i.e., along the central line ofthe parallelepiped) together with the corresponding analytical result obtained with (53). One can see both a qualitative and quantitative agreement between numerical and analytical results, the accuracy being within 0.5 degree Celsius. As expected, a better agreement is achieved at greater values of the smoothing parameter _{} in (5), which is responsible for the account of sharpness of the phase transition.

A good agreement between numerical and analytical results is also observed for the time evolution of temperature at a given point. Thus, Fig. 2 shows the time dependence of temperature at the point with coordinates _{} meters.

Fig. 2. The dependence of temperature on time t at the point _{} meters: comparison of numerical results, obtained with the scheme (16), (blue stars and green triangles) and analytical expression of the similarity solution (53) (red circles)._{} days.

V. CONCLUSIONS

We have shown that under the Assumptions 1 and 2, the uniform numerical stability (Definition 1) can be established for the Douglas – Rachford scheme (17) (which approximates the equation (1)), the general form of the sufficient stability criterion being (51). The used approach consists in reduction of a quasilinear equation to an equation with variable coefficients (subject to the constraints (25)–(27)), the forms of which change after each iteration. The comparison of obtained numerical results and an analytical solution (Fig.(1) - (2)) suggests that the ADI method can be successfully applied to the solution of Stefan problems in permafrost ground.

ACKNOWLEDGMENT

REFERENCES

[1] J. Douglas, “Alternating direction methods for three space variables,” Numerische Mathematik, vol. 4, no. 1, pp. 41–63, Dec. 1962.

[2] N. N. Yanenko, The method of fractional steps;: The solution of problems of mathematical physics in several variables, 1st ed. Springer-Verlag, 1971.

[3] G. I. Marchuk, Splitting and Alternating Direction Methods, ser. Finite Difference Methods. North Holland: Elsevier Science Publishers B.V., 1990, vol. 1, pp. 199–462.

[4] K. J. in ’t Hout and B. D. Welfert, “Unconditional stability of secondorder ADI schemes applied to multi-dimensional diffusion equations with mixed derivative terms,” Applied Numerical Mathematics, vol. 59, no. 3-4, pp. 677–692, Mar. 2009.

[5] P. Angot, J. Keating, and P. D. Minev, “A direction splitting algorithm for incompressible flow in complex geometries,” Computer Methods in Applied Mechanics and Engineering, vol. 217-220, pp. 111–120, 2012.

[6] E. Hansen and A. Ostermann, “Dimension splitting for quasilinear parabolic equations,” IMA Journal of Numerical Analysis, vol. 30, no. 3, pp. 857–869, 2010.

[7] G. Birkhoff and R. S. Varga, “Implicit alternating direction methods,” Trans. Amer. Math. Soc., vol. 92, pp. 13–24, 1959.

[8] B. Bialecki and R. Fernandes, “An alternating-direction implicit orthogonal spline collocation scheme for nonlinear parabolic problems on rectangular polygons,” SIAM Journal on Scientific Computing, vol. 28, no. 3, pp. 1054–1077, 2006.

[9] H. R. Thomas, P. Cleall, Y. C. Li, C. Harris, and M. Kern-Luetschg, “Modelling of cryogenic processes in permafrost and seasonally frozen ground,” G´eotechnique, vol. 59, no. 3, pp. 173–184, 2009.

[10] O. A. Ladyˇzenskaja, V. A. Solonnikov, and N. N. Ural’ceva, Linear and quasi-linear equations of parabolic type, ser. Translations of Mathematical Monographs. Providence, RI: American Mathematical Society, 1968, vol. 23.

[11] S. Kamenomostskaya, “On the Stefan problem.” Mat. Sb., N. Ser., vol. 53(95), no. 4, pp. 489 –514, 1961.

[12] A. A. Samarskii and P. N. Vabishchevich, Mathematical Modelling, Volume 1, Computational Heat Transfer, volume 1 ed. Wiley, 1996.

[13] L. Bronfenbrener and E. Korin, “Two-phase zone formation conditions under freezing of porous media,” Journal of Crystal Growth, vol. 198-199, pp. 89–95, 1999.

[14] M. Lees, “Alternating direction and semi-explicit difference methods for parabolic partial differential equations,” Numerische Mathematik, vol. 3, no. 1, pp. 398–412, 1961.

[15] A. A. Samarskii and A. V. Gulin, Stability of Difference Schemes. Moskow: Nauka, 1973.

[16] J. E. Dendy, “Alternating direction methods for nonlinear Time- Dependent problems,” SIAM J. Numer. Anal., vol. 14, no. 2, pp. 313–326, 1977.

[17] E. Javierre, C. Vuik, F. J. Vermolen, and S. van der Zwaag, “A comparison of numerical models for one-dimensional stefan problems,” Journal of Computational and Applied Mathematics, vol. 192, no. 2, pp. 445–459, 2006.