Lobachevskii Journal of Mathematics Vol. 13, 2003, 25 – 38 \begin{matrix}c \\ ◯ \end{matrix}R.F. Kadyrov, E. Laitinen, and A.V. Lapin

R.F. Kadyrov, E. Laitinen, and A.V. Lapin

USING EXPLICIT SCHEMES FOR CONTROL PROBLEMS IN CONTINUOUS CASTING PROCESS

Abstract. In this article, we solve an optimal control problem of the cooling process in the steel continuous casting, which mathematical formulation is a coefficient identification problem for Stefan problem with prescribed convection. To minimize a cost function we use a gradient method, state and adjoint state problem being approximated by explicit mesh schemes with variable time steps. Presented numerical results show an advantage in calculations time of this approach in comparison with using implicit mesh schemes.

1. Introduction

A continuous casting process can be described by the Stefan problem with prescribed convection [1, 2]. The existence of a unique solution for this problem is proved in [3]. For numerical solving of this problem the fully implicit schemes or the implicit schemes with characteristics approximation of heat transfer operator are traditionally used (cf., e.g., [4, 5, 6] and bibliographies therein). These schemes are unconditionally stable, but SOR-types methods which are commonly used for its numerical solving have a slow convergence rate. More sophisticated methods are based on the domain decomposition [7, 8] and/or multigrid procedures [9]. Nevertheless, when using any of these approaches we get a system of nonlinear equations which have to be solved again by using SOR method.

On the other hand, a computational complexity of one step of an explicit scheme is the same as for one SOR-iteration. It is well known, that the explicit schemes with constant steps in time are only conditionally stable. This essentially restricts the field of their application in solving the applied problems. At the same time, in [10, 11, 12] the effective algorithms for solving both linear, and nonlinear non-stationary problems have been suggested, these algorithms being based on the explicit schemes with variable steps.

2000 Mathematical Subject Classification. Primary. 65M60,49J20.

Key words and phrases. Stefan problem, continuous casting process, optimization problem, finite element method, explicit scheme.

This work is supported by RFBR, project N 01-01-00068..

In this paper, the results of using the explicit schemes with variable steps for solving a 3D dynamic control of the secondary cooling process in continuous casting problem are presented.

The rest of the article is organized as follows. In the second section the boundary-value problem is posed in two formulations: a temperature formulation (unknown is the temperature field) and enthalpy formulation (unknown is enthalpy function). To solve the problem in temperature formulation we use an implicit scheme with FEM approximation in space variables, while for the problem in the enthalpy formulation an explicit Euler method with cycle of variable steps with the same FEM approximation in space is used. In the third section of the paper a cooling optimization problem for the continuous casting process is described. We construct a method for solving this problem which is based on the explicit approximations with variable steps of both direct and adjoint state problems. In the last section some numerical results are presented. The main conclusion of these numerical results is: the calculated control functions and temperature fields are very close for both methods, while the computational complexity of the explicit scheme is less than for the implicit one. This fact is very important for our applied problem which must be solved in real-time regime.

2. Boundary-value problem formulation

Let \overline{Ω} ∈ R^{2} be the rectangular domain [0, L_{x}] × [0, L_{y}] with round left-down corner of radius R ⩾ 0, \overline{Ω} = Ω ∪ ∂Ω. Let V = Ω × (0, L_{z}) be the 3D domain with boundary Γ = ∂V which consists of the following parts:

Γ_{0} = \overline{Ω} × \{0\} ,
Γ_{N} = {(x, y) ∈ ∂Ω : x = 0 ∨ y = 0} × (L_{M}, L_{z}] ∪ Ω × {L_{z}} ,
Γ_{S} = {(x, y) ∈ ∂Ω : x ≠ 0 ∧ y ≠ 0} × (0, L_{z}] ,
Γ_{M} = {(x, y) ∈ ∂Ω : x = 0 ∨ y = 0} × (0, L_{M}] ,

where 0 ⩽ L_{M} < L_{z}.

In our applied problem domain V represents a metallic bar (or melt), Ω is the quarter of z-cross section of a slab, L_{z} is the length of a slab and L_{M} is the length of primary cooling zone.

We find the temperature field T = T (x, y, z, t) for (x, y, z) ∈ \overline{Ω} and t ∈ (0, T_{f}] (T_{f} is the time of casting process ending).

We define the enthalpy function H(T) and the Kirchoff’s temperature K(T) by

H (T) = ρ (∫_{0}^{T} c (ξ) dξ + L (1 − f_{s} (T))), K (T) = ∫_{0}^{T} k (ξ) dξ ,

Figure 1. Domain and boundaries

where ρ, c(T), k(T) and L are density, heat capacity, thermal conductivity and latent heat, f_{s} (T) is the solid fraction at temperature T, such that

f_{s} (T) = \{ \begin{matrix}1, T > T_{s}, \\ 0, T < T_{l}. \end{matrix}

Here T_{s} is the solidification and T_{l} is the melting temperature. On Figures 2, 3 the graphs of enthalpy function and Kirchhoff’s temperature for steel are presented. Note that because of the monotonicity and continuity of enthalpy function H = H(T), we can define the inverse function T = T (H).

Figure 2. Enthalpy function Figure 3. Kirchhoff’s temperature function

The continuous casting process can be modelled by the following boundaryvalue problem:

\{ \begin{matrix}\frac{∂H(T)}{∂t} + v \frac{∂H(T)}{∂z} − ΔK (T) = 0 in V × (0, T_{f}] , \\ T = T_{0} on Γ_{0} × (0, T_{f}] , \\ \frac{∂K (T)}{∂n} + h (T − T_{w}) + σε (T^{4} − T_{ext}^{4}) = 0 \\ on Γ_{N} × (0, T_{f}] , \\ \frac{∂K (T)}{∂n} = Q on Γ_{M} × (0, T_{f}] , \\ T (x, y, z, 0) = T^{0} in V. \end{matrix}  (1)

Here n is the unit vector of outward normal to ∂V, σ > 0 is the Stefan-Bollzmann constant, ε > 0 is the emissivity, v > 0 is the casting speed, h = h (x, y, z, t) and Q = Q (x, y, z) are known cooling functions, T_{w} = T_{w} (x, y, z) and T_{ext} = T_{ext} (x, y, z) are known temperatures of cooling liquid and environment, T_{f} is the time of casting process ending.

Besides of the temperature-enthalpy formulation (1) we use the following enthalpy formulation:

\{ \begin{matrix}\frac{∂H}{∂t} + v \frac{∂H}{∂z} − ΔK (H) = 0 in V × (0, T_{f}] , \\ H = H(T_{0}) on Γ_{0} × (0, T_{f}] , \\ \frac{∂K (H)}{∂n} + h (T (H) − T_{w}) + σε (T^{4} (H) − T_{ext}^{4}) = 0 \\ in Γ_{N} × (0, T_{f}] , \\ \frac{∂K (H)}{∂n} = Q in Γ_{M} × (0, T_{f}] , \\ H (x, y, z, 0) = H^{0} on V. \end{matrix}  (2)

Unknown function in (2) is the enthalpy function H(·), by K(H) we mean here K(T (H)).

3. Mesh problem and its solving

Let for some ε > 0 we know τ = τ (ε), such that the following condition holds:

\frac{∂H}{∂t} + v \frac{∂H}{∂z} = \frac{H − H̃}{τ} + ω (t) , with ∥ω (t) ∥ ⩽ ε, ∀t : 0 ⩽ t < t + τ ⩽ T_{f},

where H̃ = H (x, y, z − vτ, t − τ). In fact, this means that for this τ the local ε-approximation of the term ∂H/∂t + v∂H/∂z in a norm ∥ · ∥ is achieved.

The weak formulation of the semi-discrete scheme to problem (1) for a fixed time level can be written as follows:

\begin{matrix}find H ∈ H^{1}(V), H |_{Γ_{0}}= H(T_{0}) : \\ \frac{1}{τ} ∫_{V} (H (T) − H̃)τηdx + ∫_{V} ∇K (T) · ∇ηdx+ \\ + ∫_{Γ_{N}} \{h (T − T_{w}) + σε (T^{4} − T_{ext}^{4})\} ηdx = 0, \\ ∀η ∈ V^{0} = \{v ∈ H^{1} (V) : η|_{Γ_{0}} = 0\} . \end{matrix}  (3)

Here H^{1} (V) is the usual Sobolev space.

Let Ω be partitioned into a set of quadrilateral finite elements, this triangulation being topologically equivalent to rectangular one. The triangulation of Ω and partitioning of [0, L_{z}] form the triangulation of 3D domain V into a set of prismatic 3D finite elements. Now we can write the finite element (FEM) discretization of problem (3). Let V_{h} be a finite element approximation of the space H^{1}(V), based on the polylinear elements. Further w_{h} stands for a V_{h}-approximation of function w (for a continuous function w w_{h} is its V_{h}-interpolant). Let

V_{h}^{T_{0}} = {v_{h} ∈ V_{h} (V) : v_{h} = T_{0h} on Γ_{0}} ,
V_{h}^{0} = {v_{h} ∈ V_{h} (V) : v_{h} = 0 on Γ_{0}} .

Approximation of problem (3) by finite element method is defined as

\begin{matrix}\frac{1}{τ} ∫_{V} (H_{h} (T_{h}) − H̃_{h})ηdx + ∫_{V} ∇K_{h} (T_{h}) · ∇ηdx+ \\ + ∫_{Γ_{N}} \{h (T_{h} − T_{w,h}) + σε ((T^{4})_{h} − (T_{ext}^{4})_{h})\} ηdΓ_{N} = 0, \\ ∀η ∈ V_{h}^{0}. \end{matrix}  (4)

Discrete problem (4) is equivalent to the system of nonlinear algebraic equations (hereafter we use the same notations for continuous and mesh functions):

\frac{H (T) − H̃}{τ} + AK (T) + B (h) (T − T_{w}) + B (T^{4} − T_{ext}^{4}) = 0.  (5)

Here: A is the discrete Laplas operator, in current FEM realization A is the 11diagonal matrix; B, B (h) are the diagonal matrices with sparse main diagonal.

In enthalpy formulation problem (5) becomes

\begin{matrix}\frac{H − H̃}{τ} + AK (H) + B (h) (T (H) − T_{w}) + \\ +B (T (H)^{4} − T_{ext}^{4}) = 0. \end{matrix}  (6)

Below we consider the problem in enthalpy statement (6). Let

f (H) = AK (H) + B (h) (T (H) − T_{w}) + B (T (H)^{4} − T_{ext}^{4}) .

Then an explicit Euler method can be written as

H^{k+1} = H^{k} − τ^{k+1}f (H^{k}) , H^{0} = H (T^{0}) , k = 0, 1, 2, …  (7)

Here τ^{k} are the time steps, H^{k} = H (t = ∑_{i=1}^{k} τ^{i}) is the numerical solution at

the time level t.

The Frechet derivative of the mapping f (H) is given by

J(H) = A ∘ K′(H) + B (h) ∘ T′(H) + 4B ∘ (T^{3}(H)T′(H)) .

Let (λ, ξ) is the eigen pair of operator J, then

(K′ + A^{−1}B (h) T′ + 4A^{−1}BT^{3}T′) ξ = λA^{−1}ξ.

Hereafter we omit the argument H of functions T, T′ and K′. Executing some technical estimates we get the following inequalities:

(K′ξ, ξ) + (A^{−1}B (h) T′ξ, ξ) + 4 (A^{−1}BT^{3}T′ξ) ⩽
⩽ (L_{K}^{H}λ_{max}^{A}) (A^{−1}ξ, ξ) + 4 (A^{−1}BT^{3}T′ξ, BT^{3}T′ξ) (A^{−1}ξ, ξ)^{1/2} +
+ (A^{−1}B (H) T′ξ, B (H) T′ξ)^{1/2} (A^{−1}ξ, ξ)^{1/2} ⩽
⩽ (L_{K}^{H}λ_{max}^{A} + cond (A)^{1/2} L_{T}^{H} ∥B (h)∥ +
+4 cond (A)^{1/2} L_{T}^{H}L_{T^{3}}^{H}) (A^{−1}ξ, ξ) .

Now, maximum eigenvalue of operator J can be estimated by

λ_{max} ⩽ L_{K}^{H}λ_{max}^{A} + cond (A)^{1/2} L_{T}^{H} ∥B (h)∥
+4 cond (A)^{1/2} L_{T}^{H}L_{T^{3}}^{H}.

Here:

• λ_{max}^{A} is the maximum eigenvalue of the operator A;

• cond (A) = \frac{λ_{max}^{A}}{λ_{min}^{A}} is the condition number of operator A;

• L_{K}^{H} = sup_{H∈[H(0),H(T_{s})]} (K′ (H));

• L_{T}^{H} = sup_{H∈[H(0),H(T_{s})]} (T′ (H));

• L_{T^{3}}^{H} = sup_{H∈[H(0),H(T_{s})]} (T^{3} (H)).

Let τ^{k} = τ_{0}, then we have the explicit Euler method with constant step. Approximation condition in this case is defined by following inequality: τ_{0} ⩽ τ. Stability condition can be written as τ_{0} ⩽ cou, where Courant number cou = 2/λ_{max}. As it is known, for large value of λ_{max} the last condition is too restrictive and too many steps in time can be necessary when using an explicit scheme.

Let {τ^{k}}_{k=1}^{N} be a cycle of time steps in (7), N be the length of the cycle. The total step of the cycle is l_{N}(τ_{max}) = ∑_{k=1}^{N} τ^{k}, where τ_{max} is the maximum step of the cycle. The approximation condition of this scheme can be written as

τ_{max} ⩽ τ(ε)  (8)

The relaxed stability condition of this so-called super-time-stepping scheme is

|∏_{k=1}^{N} (1 − τ_{k}λ)| < 1 ∀λ : (λ_{min} ⩽ λ ⩽ λ_{max}),  (9)

where λ_{min} and λ_{max} is the minimum and maximum eigenvalues of the operator J.

To obtain an optimal scheme, a maximum superstep size must be found as maximum l_{N}(τ(ε)) when (8) and (9) hold. An optimal cycle is found by using the optimal properties of Chebyshev polinomials in [10] for the cases N = 2^{p} and N = 2^{p}3^{q}. The construction of this scheme is based on the using the following parameters: τ_{max} ⩽ τ(ε) – the maximum step of the cycle, cou – Courant number of the operator.

This scheme is described for the operators with real spectrum. In our problem operator J is symmetrizable, thus it has only real spectrum.

In our case we have parameter τ_{f} ≫ τ(ε) and we must calculate solution for every time t = kτ_{f}, k = 0, 1, 2, … It means that we need the time steps cycle such that l_{N} ≈ τ_{f}. Obviously, function g(x) = |τ_{f} − l_{N}(x)|, 0 < x ⩽ τ(ε) is unimodal for every fixed N, and it means that for all fixed ξ > 0 we can easily find τ : g(τ) ⩽ ξ.

Values N = 2^{p} and τ were defined experimentally. Using explicit scheme with cyclic set of steps we get some advantage in computing complexity in comparison with the scheme with a constant step, but it has not a theoretical N-order advantage, because we have condition l_{N} ≈ τ_{f}. In the next section we also compare explicit schemes with variable steps with standard implicit method (see Table 1). The solution of nonlinear system arising in the implicit scheme was made by Gauss-Seidel method with preliminary Newton linearization. In column ”Number of iterations” of Table 1 we present: for the implicit scheme – average number of the iterations required for the solving of corresponding system of the nonlinear equations at one time level, for the explicit scheme – number of steps in a cycle. In column ”Time of iteration” the average time of one iteration in milliseconds is recorded.

4. Optimization problem

The secondary cooling region is divided into M cooling zones, in every one the heat transfer coefficient h_{k}, k = \overline{1, M} have to be controlled. In optimal control method our aim is to minimize a goal function which is constructed by means of metallurgical cooling criteria. By these criteria we have two observation lines on the slab surface (along the corner and center lines of slab surface): l_{1} = ((x, y, z) ∈ Γ : x = y = R(1 − 1/\sqrt{2}), z ∈ [L_{M}, L_{z}]), l_{2} = ((x, y, z) ∈ Γ : x = L_{x}, y = 0, z ∈ [L_{M}, L_{z}]) and given target temperatures \{t_{k}^{tar}, k = \overline{1, M}\} for starting point of each cooling zone along this lines. Let T^{tar}(x), x ∈ l = l_{1} ∪ l_{2} is the linear interpolation of {t_{k}^{tar}}. Thus, the goalfunction can be defined as

I(T) = \frac{1}{2} ∫_{l} (T^{tar} − T)^{2}dx.

Let J(T) is the mesh approximation of I(T), vector J′(T) is the derivative of the goal-function. We can formulate the optimal control problem in the following way:

Find h^{∗} ∈ U such that:

J(h^{∗}) = min_{h∈U} J(h),  (10)

where U = \{u = (u_{1}, u_{2}, …, u_{M}) : h_{min} ⩽ u_{i} ⩽ h_{max}∀i = \overline{1, M}\} is the set of constraints for control variables and T = T (h) is defined via solution of mesh scheme (6) with the set of variable time steps. Thus, we have a nonlinear function minimization problem with constraints. To solve problem (10) we use a gradient method.

To calculate the gradient of the goal-function we construct the adjoint state problem by using the Lagrange function

L = J(T^{N})+

+ ∑_{k=1}^{N−1} (\frac{H^{k+1} − H^{k}}{τ_{k}} + AK(T^{k}) + B(h)(T^{k} − T_{w}) + B ((T^{k})^{4} − T_{ext}^{4}) , λ^{k+1}).

Here we use the inner product (·, ·) in the finite dimensional space of vectors. Variation gives

L′_{y}v = (J′(T^{N})) + ∑_{k=1}^{N−1} (\frac{v^{k+1} − v^{k}}{τ_{k}} + A (K′(T^{k})(T^{k})′_{H}v^{k}) , λ^{k+1})+

+ ∑_{k=1}^{N−1} (B(h)(T^{k})′_{H}v^{k} + 4B(T^{k})^{3}(T^{k})′_{H}v^{k}, λ^{k+1}) = 0.

Obviously,

∑_{k=1}^{N−1} \frac{v^{k+1} − v^{k}}{τ_{k}} λ^{k+1} = ∑_{k=1}^{N−1} v^{k} (\frac{λ_{k}}{τ_{k}} − \frac{λ^{k+1}}{τ_{k+1}} ) + \frac{v^{N}λ^{N}}{τ_{N}} , λ^{1} = 0

From here we derive the following system of the linear algebraic equations, which is the mesh adjoint state problem:

\{ \begin{matrix}\frac{λ^{N}}{τ_{N}} + J′(T^{N})(T^{N})′_{H} = 0 \\ \frac{λ^{k}}{τ_{k}} − \frac{λ^{k+1}}{τ_{k+1}} + (T^{k})′_{H}K′(T^{k})A^{∗}λ^{k+1} + (T^{k})′_{H}B^{∗}(h)λ^{k+1}+ \\ +4(T^{k})′_{H}(T^{k})^{3}B^{∗}λ^{k+1} = 0, k = N − 1, …, 1. \end{matrix}  (11)

As it is shown in [10] scheme (11) is stable and stable to round-off errors. Now the descend direction is defined by

(J′(h))_{j} = (B(h_{j})(T − T_{w}), λ^{1}),

where: h_{j} = (0, …, 0, \overbrace{1, …, 1}^{i_{k}^{j}}, 0, …, 0)^{T} , j = \overline{1, M}; i_{k}^{j} - all mesh points in the

j-th cooling zone. Vector (−J′(h)) is the descent direction. The projection p = (p_{1}, p_{2}, …, p_{M}) of the vector −J′(h) to constrains set U can be calculated as

p_{i} = \{ \begin{matrix}0, (h_{i} = h_{min} ∧ J′_{j} > 0) ∨ (h_{i} = h_{max} ∧ J′_{j} < 0), \\ −J′_{i}, otherwize. \end{matrix}

Then we consider the one-dimensional minimization problem:

ρ_{opt} = arg min_{0⩽ρ⩽ρ_{max}} J(h + ρp),  (12)

where maximum descent step ρ_{max} can be found from: start point h, current descent direction p and constraints set U. In this work we solve (12) by the Brent algorithm of one-dimension minimization and by the following algorithm: Step 1.

Let h(ρ) = h + pρ

ρ_{0} = \frac{ρ_{max}}{2^{d}}, where d > 0 is the arbitrary parameter of algorithm

If J(h(ρ_{0})) > J(h(0)) Then ρ_{opt} = 0 End

i = 1

GoTo Step 2 Step 2.

ρ_{opt} = ρ_{i−1}

ρ_{i} = 2ρ_{i−1}

If (ρ_{i} > ρ_{max}) or (J(h(ρ_{i})) > J(h(ρ_{i−1}))) End

i = i + 1

GoTo Step 2. The solutions of (10) using both one-dimensional minimization algorithms are the same, while the second algorithm is much simpler than Brent algorithm.

The results of numerical experiments (optimal control and temperature fields calculated by using both explicit and implicit schemes) are presented on Figures 4 – 9.

34 R.F. KADYROV, E. LAITINEN, AND A.V. LAPIN

5. Numerical results

Direct solution.

Method Number of iterations

Time of one iteration (ms)

\begin{matrix}Mesh: 11x17x221 \\ T_{f} = 1000 \end{matrix} Explicit scheme  8  116
τ_{f} = 5  Implicit scheme  24  256
\begin{matrix}Mesh: 8x16x150 \\ T_{f} = 1000 \end{matrix} Explicit scheme  8  56
τ_{f} = 5  Implicit scheme  19  107
\begin{matrix}Mesh: 8x16x150 \\ T_{f} = 900 \end{matrix} Explicit scheme  4  28
τ_{f} = 3  Implicit scheme  14  73

Table 1. Number of iterations and time of calculation for solution of the state problem.

Optimization.

Figure 4. Optimal cooling control results by implicit scheme

EXPLICIT SCHEMES FOR CONTROL PROBLEMS 35 Figure 5. Optimal cooling control results by explicit scheme

Optimal and target temperatures Figure 6. Temperature by implicit scheme along corner line 36 R.F. KADYROV, E. LAITINEN, AND A.V. LAPIN

Figure 7. Temperature by explicit scheme along corner line Figure 8. Temperature by implicit scheme along center line EXPLICIT SCHEMES FOR CONTROL PROBLEMS 37

Figure 9. Temperature by explicit scheme along center line

References

[1] Laitinen E., Neittaanmäki P. On numerical simulation of the continuous casting problem //J. Eng. Math. – 1987. – V. 22. – P. 335–354.

[2] Louhenkilpi S., Laitinen E., Nieminen R. Real time simulation of heat transfer in continuous casting // Metall. Trans. B. – 1993. – V. 24B. – P. 685–693.

[3] Rodrigues J. F., Yi F. On a two-phase continuous casting Stefan problem with nonlinear flux // Euro J. App. Math. – 1990. – V. 1. – P. 259–278.

[4] Chen Z., Jiang L. Approximation of a two phase continuous casting problem // J. Part. Diff. Equations. – 1998. – V. 11. – P. 59–72.

[5] Laitinen E., Lapin A. Implicit approximation and iterative solution for the continuous casting problem // Preprint, July, 1999, Dep. of Math. Sci., University of Oulu. – 1999.

– 17 p.

[6] Laitinen E., Lapin A., Pieskä J. Mesh approximation and iterative solution of the continuous casting problem // in: ENUMATH 99. Ed. by P. Neittaanmäki, T. Tiihonean and P. Tarvainen. World Scientific, Singapore. – 2000. – P. 601–617.

[7] Laitinen E., Saranen J., Pieskä J., Lapin A. Comparison of domain decomposition methods for solving continuous casting problem // in: Domain Decomposition Methods DDM-13. Ed. by N. Debit, M. Garbey, R. Hoppe, J. Periaux, D. Keyes, Y. Kuznetsov. CIMNE, Barselona. – 2002. – P. 411–418.

[8] Laitinen E., Lapin A., Pieskä J. Asinchronous domain decomposition methods for solving continuous casting problem // J. of Comp. and Appl. Math. – 2003. – V. 154. – P. 393–413.

[9] Dautov R.Z., Ignatieva M.A. Multigrid solution of a two-phase Stefan problem with prescribed convection // Proceedings of Russian-Finnish Workshop: Numerical methods for continuous casting and related problems, V. 9. – Kazan, DAS Publisher, 2001, – pp. 3-11.

[10] Lebedev V.I. Functional analysis and calculus mathematics. – M.: Science, 2000. – 290 P.

38 R.F. KADYROV, E. LAITINEN, AND A.V. LAPIN

[11] Lebedev V.I. Explicit difference schemes for solving stiff systems of ODEs and PDEs with complex spectrum//Russ. J. Numer. Anal. Math. Modelling. – 1998. – V. 13. – No 2. – P. 107–116.

[12] Lebedev V.I. Explicid difference schemes for solving stiff schemes with complex or partitioned spectrum//JNM and MP. – 2000. – V. 40. – No 12. – P. 1801–1812.

Kazan State University, Russia University of Oulu,Finland Kazan State University, Russia E-mail address: alapin@ksu.ru

Received October 15, 2003