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.