Symmetry analysis of equations for convection in binary mixture Текст научной статьи по специальности «Математика»

УДК 532.517

Symmetry Analysis of Equations for Convection in Binary Mixture

Ilya I.Ryzhkov*

Institute of Computational Modelling SB RAS Krasnoyarsk, 660036

Russia

Received 10.08.2008, received in revised form 10.10.2008, accepted 06.11.2008 The differential equations describing convection in binary mixture with Soret and Dufour effects are considered. The symmetry classification of these equations with respect to the constant parameters is made. It is shown that a generator producing equivalence transformations of constants is defined accurately up to a factor arbitrarily depending on these constants. The equivalence group admitted by the governing equations is calculated. Using this group, a transformation connecting the systems with and without Soret and Dufour terms is derived. In pure Soret case, it reduces to a linear change of temperature and concentration. The presence of Dufour effect requires an additional change of thermal diffusivity and diffusion coefficient. A scheme for reducing an initial and boundary value problem for Soret-Dufour equations to a problem for the system without these effects is proposed.

Keywords: Lie symmetry group, equivalence transformation, binary mixture, convection, Soret and Dufour effects

Introduction

It is well known that symmetry analysis provides a powerful tool for studying partial and ordinary differential equations [1]. This method is especially fruitful in application to the equations of physics and mechanics since many of them are derived on the basis of invariance principles. Symmetries give important information about qualitative properties of differential equations and provide a basis for their classification and simplification. Moreover, they can be effectively used for constructing exact solutions of non-linear problems.

This paper deals with symmetry analysis of equations describing convective motion in a binary mixture. Convection in binary mixtures has been an active field of research over the last decades [2-4]. The flow dynamics in mixtures is more complex than that of one-component fluids due to an interplay between advection and mixing, solute diffusion, and Soret and Dufour effects (also known as cross-effects). The Soret effect (or thermal diffusion) is a molecular transport of substance associated with a thermal gradient. It results in component separation in non-uniformly heated fluid. The situation when a concentration gradient induces heat transfer is called the Dufour effect. This effect is important only in gas mixtures.

We assume that the Oberbeck-Boussinesq approximation is valid. Then the deviations of temperature T and concentration C (expressed as a mass fraction) from their mean constant values To and Co are relatively small and the mixture density can be written in the form

p = po(1 - PtT - paC).

* e-mail: [email protected] © Siberian Federal University. All rights reserved

Here po is the mixture density at the mean values of temperature and concentration, Pt and pa are the thermal and concentration expansion coefficients respectively. It is supposed that C is the concentration of the lighter component, so pa > 0. The equations of motion have the form [2,5]

ut + (u ■ V)u = - — Vp + vAu - g(PtT + paC), po

Tt + u ■VT =(x + a2DN)AT + aDNAC, (1)

Ct + u ■ VC = DAC + aDAT, divu = 0,

where x = (xl,x2,x3) is the coordinate vector; u = (u1,u2,u3) is the velocity vector; p is the difference between actual and hydrostatic pressure; g = (0, 0, -g), where g is the gravitational acceleration; v is the kinematic viscosity; x is the thermal diffusivity; D is the diffusion coefficient. We suppose that all transport coefficients are constant and correspond to the mean values of temperature and concentration. The Soret effect is characterized by the parameter

a = C0(l - C0)DT/D,

where Dt is the thermal diffusion coefficient. The case a < 0 corresponds to positive Soret effect, when the lighter (heavier) component is driven towards the higher (lower) temperature region. In case of negative Soret effect we have a > 0, and the opposite situation is observed. Thermodynamic parameter

N

T (

cp\dC j T,p

>0

To, Co

is responsible for Dufour effect. Here cp is the specific heat at constant pressure, ^ is the chemical potential of the mixture. The expression in square brackets is taken at the mean values of temperature and concentration.

Symmetry properties of two-dimensional equations describing free convection in one-component fluid were studied in [6] (see also [7]). The case of stationary plane flows was considered in the earlier paper [8]. Symmetry analysis of equations for convection with Soret effect was first performed in [9]. Invariant solutions of rank 1 and 2 were classified in [10,11] (the first and second order optimal systems of subalgebras for the admissible Lie symmetry algebra were constructed). Invariant solutions describing convective motion in binary mixture with Soret effect in plane and cylindrical vertical layers have been recently studied in [12,13].

This paper is organized as follows. In section 1, symmetry classification of system (1) with respect to the nine constant parameters is performed. The one-parameter transformation subgroups produced by the generators of the admissible Lie symmetry algebras are described. Section 2 is devoted to studying the equivalence properties of the governing equations. It is shown that a generator producing equivalence transformations of constants is defined accurate to a factor arbitrarily depending on these constants. In section 3, we use the equivalence group to derive a transformation connecting the equations with and without Soret and Dufour terms. In pure Soret case, the number of transformed variables is minimized to temperature and concentration. The presence of Dufour effect also requires a change of thermal diffusivity and diffusion coefficient.

A scheme for reducing an initial and boundary value problem for equations (1) to a problem for the system without cross-effects is proposed.

1. Symmetry Properties of the Governing Equations

The equations of motion (1) contain nine constant parameters. Eight of them specify the properties of fluid, and one defines the gravitational acceleration. In this section, we consider the symmetry classification of system (1) with respect to these parameters. It is supposed that a, N, pT, Pc can take zero values (in this case, the corresponding terms in the equations are omitted). This approach allows us to study the symmetries of different models that can be obtained from system (1) by taking into account or neglecting Soret and Dufour effects or dependence of density on temperature and concentration. Note that the case of zero gravity (g = 0) is equivalent to the case Pt = Pc = 0 from symmetry analysis point of view. In what follows we assume that the parameters po, g, v, x, D are positive according to their physical meaning.

Let us introduce the following notation. If f (t, x) is an arbitrary function, then its derivatives are denoted as follows:

df _ , f d2f d2f_ . . _ . . dt = fu dx = fu dtdx = ^ dXdX = fij, i,j = 1,23 i -3

Using this notation, we rewrite equations (1) in coordinate form:

uj + v}u\ + u2u2 + u3u3 + p-1 pi — v (uj;1 + + U33) = 0, (2)

2,12,22,32,-1 ( 2,2,2 \ n ro\

Ut + U Uj + U U2 + U U3 + p- p2 — V(Ujj + U22 + U33) =0, (3)

u3 + UjU3 + u2u3 + u3u3 + p-1 P3 — V (u31 + u22 + u33) — g(pT T + pc C) = 0, (4) Tt + u1T1 + u2T2 + u3T3 — (x + a2DN )(Tn + T22 + T33) —

— aDN (C11 + C22 + C33) = 0, (5)

Ct + u1C1 + u2C2 + u3C3 — D(Cn + C22 + C33) — aD(Tn + T22 + T33) = 0, (6)

U 1 + U2 + U3 = 0. (7)

To find the admissible Lie symmetry group, we calculate the corresponding Lie symmetry algebra of infinitesimal generators. The admissible generator for equations (1) is sought in the form

X=itdt++++nTdr+nCdC <8»

supposing that its coordinates depend on all dependent and independent variables (summation over i = 1, 2, 3 is assumed). To derive the determining equations, we need to apply the prolongation generator X to equations (2) - (7) and make the transition to the manifold given by this

system. However, the equations are not in involution, which makes it difficult to choose the external and internal variables. Let us supplement the system with its differential consequence [14]

(u 1)2 + (u2)2 + (u3)2 + 2(u2u2 + u1u? + u3u2)+ - 412 -

+P0 >11 + P22 + P33) - g(PTT3 + PcC3) = 0 (9)

obtained by differentiating equations (2), (3), (4) with respect to x^x^x3, respectively, and by using (7). When making the transition to the manifold, we also take into account the differential consequences from (7)

J'ti

+ U2 + = 0, uli + u2i + u3i =0, i =1,2,3. (10)

Equations (2) - (7), (9), (10) are in involution, and now it is easy to choose the external variables: «li, u2i, u3i, pii, T11, C11, «33, «13, «33, «33, «3. Note that T11 and C11 are expressed from (5) and (6). These equations are linear with respect to the above variables, and the solution always exists since the corresponding determinant xD is non-zero.

The determining equations are found by applying the prolongation generator X to the system and substituting the expressions for external variables in the obtained equations. After a considerable amount of calculations the solution is written in the form

Ct = 2c4t + c0, C1 = C4X1 + c1x2 + c2 x3 + f 1(t),

£2 = -c1x1 + c4x2 + c3x3 + f2(t), £3 = -c2x1 - c3x2 + c4x3 + f3(t),

n1 = —c4«1 + ci«2 + c2«3 + ft1(t), n2 = —ci«1 — c4«2 + c3«3 + ft2(t),

n3 = —c2«1 — c3«2 — c4«3 + ft3(t), (11)

np = Po(c5gPr x3 + cegPc x3 — fit (t)x1 — f2(t)x2 — f3 (t)x3) — 2c4p + f0(t),

nT = crT + cgC + c5, nC = cgC + cioT + c6.

Here co — cio are the group constants and fi(t), i = 0,1, 2, 3 are smooth arbitrary functions. The group constants are connected with the parameters of system (1) by the classifying equations

Pt (cr + 3c4) + pc cio = 0, pc (cg + 3c4) + Pt cg = 0,

aD(cg — cr) + (x — D + a2DN )cio = 0, (12)

(X — D)(cg — Ncio)=0, a(cg — Ncio)=0,

Pt c2 = 0, Pt c3 = 0, pc c2 = 0, pc c3 = 0.

Using formulae (11) and equations (12), we can find the Lie symmetry algebras admitted by the governing equations depending on the values of parameters. The results of symmetry classification are given in Table 1. The values of parameters a, Pt, Pc are specified in the first three columns. The basic generators are given in the fourth column, and the additional generators admitted in case of x = D are presented in the fifth column. Note that if a = 0, then N is considered to be zero. In case of a = 0, Pt = Pc = 0 the system admits the generator Ki, which turns into the generator K if N = 0. When a = 0, pT = 0, pc = 0, the generator K is admitted only if N = 0. Finally, when the parameters a, Pt,Pc are non-zero, the system admits the generator R3 in the special case N = No, a = pT/pc, where

N = Pc (Pt (X — D)+ PcaD) o Pt aD(PT — Pc a) . ( )

If No = 0, then the generator R3 turns into the generator Ri, which is admitted when Pt(x —

D)+ pc aD = 0.

1

Table 1. Symmetry classification of the governing equations.

a ßr ßc Generators x=D

0 0 0 1 3 1 3 Xo,Xij ,Hi,Ho,z,i ,1 , e ,e T 2,C2

0 0 =0 Xo, X12, Hi, Ho, Z2 ,U2,T1, T3 T2

0 =0 0 Xo, X12, Hi, H0, Z1, U1, C 1, C3 C2

0 =0 =0 Xo, X12, Hi, Ho, Z3, U1, U2 R1 , R2

= 0 0 0 Xo, Xjj, Hi, Ho, Z, R, T3, C3, K1

=0 0 =0 Xo, X12, Hi, Ho, Z3, U2, T3

=0 =0 0 Xo, X12, Hi, Ho, Z3, U1,C3, K (N = 0)

=0 =0 =0 Xo, X12, Hi, Ho, Z3, U1, U2

N = N0, a = ßr/ßc : R3

B i d i d i d i d

Xo = —, Xj = x ——t — xj ——t + u ——t — uj ——: Bt Bxj ox} Buj on1

B - - -

Hi(fi(t)) = fi(t) —. + fi(t) —. — poxifitt(t)i,j = 1 2, 3 (i < j),

B

' du1

Ho(f 0(t)) = f 0(t) dp, Z = 2^ + £

dt

i=i

i d i B x -77— — u —— dxi dui

B

— 2P

Bp

3 9 1 9 3 d 1 B

Ul = Pogx dp + ßr dT, U2 = Pogx dp + ßc de,

t i = T— T 2 = C-^- T3 = — дT, BT' BT'

(14)

ei = e

A C 2 = T A C 3 = A

BC, BC, BC'

R = T1 + C1, R1 = T1 — ßr C2, R2 = C1 — ßc T2,

, 1 ßc , 2 ßr ,

Z1 = Z - 3T1, Z2 = Z - 3 1, Z3 = Z - 3R,

R3 =

T ßr No C T — C

B ßT B

bt — ßC bc

K =

aT +(1 — x/D)C

B_ ~Ber

B

K1 = a.NC— +

aT +(1 — x/D — a2N )C

B

BC'

Now, let us describe one-parameter transformation subgroups that correspond to the generators presented in Table 1. These subgroups are obtained by solving the corresponding Lie equation for each generator.

Xi : t = t + a;

Xj : x® = x® cos a — xj sin a, xj = x® sin a + xj cos a; H®(/® (t)) : a?4 = x® + a/®(t), U = U + a/t®(t),

p = P - - aPo /«(t) (x 4 + a/®(t)/2) ; Hq(/0(t)) : p = P + a/0(t);

Z : 5 = e2at, X = eax®, ul = e-aw®, pT = e- 2ap ;

Ui : 5 = 3 = p + apogx . , T = T + aP-i;

U2 : 5 = = p + apogx . , C = : C + aP-i ;

Ti : T = eaT ; T2 : T = T + aC ; T3: T = T + a ;

Ci : C = eaC ; C2 : C = C + aT ; C3: C = C + a

R : T = = eaT, C = eaC ;

Ri : T = eaT, C = = C + PtP-i(1 - e°)T ;

R2 : c = eaC, t = = T + PcP-i(1 - e°)C ;

R3 : C = C +

Pt(PcT - PtNqC) [1 - exp((1 + NqPT/PC)a)]

t

(15)

NqP2 + PC

Pt Nq(Pt T + Pc C) + Pc (Pc T - Pt NqC) exp((1 + NqPT/PC )a) .

NqP2 + PC

K : C = (l - e(i-*/D) ^ T + e(1-*/D) ° C ; y - D V y

K : T =

y - D

(s1 es2° - s2 esia) T + (esia - es2°) aNC

si - S2

C =

(esia - es2a) aT + (Si esi° - s2 es2°) C

si - S2

Here a is a real parameter (every subgroup has its own parameter!), and the value a = 0 corresponds to the identity transformation. The variables that are not mentioned in the above formulae remain unchanged.

The transformation generated by Ki contains two parameters si and S2 that are roots of the quadratic equation

s2 + (a2N + y/D - 1)s - a2N = 0.

This equation always has two different real roots except for the case a = 0, y = D, in which Ki = 0. Finally, the transformations generated by Zi, Z2, Z3 are obtained by extending the one-parameter subgroup corresponding to the generator Z by the following transformations:

Zi : T = e-3aT ; Z2 : C = e-3aC ; Z3 : T = e-3aT, C = e-3aC.

The Lie symmetry algebra admitted by equations (1) in every case presented in Table 1 generate the transformation group that acts in the space of dependent and independent variables and leaves the equations invariant. The basic generators of the admissible algebra generate one-parameter subgroups of this group. Every transformation of the admissible group can be represented as a composition of transformations belonging to one-parameter subgroups.

Note that the generators Xo,X®j, H®,Ho are admitted by many models of continuum mechanics, while the generators Ui, U2, T®, C®, R, R®, Z, Z®, K, Ki (i = 1, 2,3) are specific for the equations of convection with Soret and Dufour effects. When the Dufour effect is not allowed for (N = 0), the results presented here agree with those reported in [9]. It can be shown that symmetry classification of two-dimensional equations (1) is obtained by restricting the action of the admissible Lie algebras and the corresponding Lie groups to the space of variables t, x1, x3, u1, m3, p, T, C. The case of plane Oberbeck-Boussinesq equations corresponds to C = 0, a = 0, pT = 0, pC =0, and the admissible generators are X0, X12, H1, H3, H0, Z1, U1. It is consistent with the previous results [6,7].

It should be noted that the equivalence transformations of parameters are not taken into account when performing the symmetry classification (they will be calculated in the next section). These transformations are usually used to simplify the arbitrary elements (parameters or functions) entering into the equations [1]. When arbitrary elements are constants, the aim is to set as many constants as possible to zero or to unity. However, the so-obtained equations do not possess the necessary physical parameters, and using them for constructing physically meaningful solutions with a help of symmetries is not convenient. In contrast to this approach, the classification presented here shows the dependence of symmetry properties on physical effects incorporated in the model.

In this section we also calculate the admissible reflections of dependent and independent variables for system (1). Let us introduce new variables

5 = (-1)ao t, x® = (-1)ai x®, M® = (-1)Y «®, i = 1,2,3, p = (-1)74 p, T = (-1)Y5 T, C = (-1)Y6 C,

where ao,..., a3,71,..., 76 € {0,1}, and substitute them into equations (2) - (7). Multiplying equation (2) by (- 1)Yl gives

(-1)a0U1 + (-1)ai+Yl M1«! + (_1)«2+Y2u2M1 + (-1)^3+73 u3u3 = = ( 1)ai +Yl +Y4 p0-1pi + v(uii + u22 + m!3).

This equation remains unchanged if all the powers of (-1) are zero. Applying a similar procedure to the other equations, we obtain the conditions

ao = ai + 71 = a2 + 72 = a3 + 73 = 74 = 73 + 75 = 73 + 76 =0 (16)

for the invariance of system (1). As there are 10 parameters satisfying seven equations (16), the three parameters a\,a2, a can be considered arbitrary. Successively setting one of them to unity and the remaining ones to zero, we find the admissible reflections

di : x1 = — x1, u1 = — u1; : x2 = — x2, u2 = — u2;

d3 : x3 = —x3, u3 = —u3, T = —T, C = —C. (17)

These reflections are discrete transformations of the governing equations. Note that the transformation d3 is physically meaningful since T and C are the deviations from the mean values, and these deviations can be negative.

2. Equivalence Transformations

Let us now turn to studying the equivalence transformations of the governing equations. First of all, we prove a statement that will be widely used in the further analysis. Consider a system of differential equations

E(x, u,u,... ,u,c) =0.

(18)

Here x = (x1,... ,xn) and u = (u1,..., um) are vectors of independent and dependent variables respectively, c = (c1,..., c q) is a vector of constant parameters, u is a collection of derivatives of

r

u with respect to x of order r = 1, ... k:

u = iuJi1...ir I j = l,...,m; h,...,ir = l,...,n},

where the following notation is used

j

ui1...ir

d r uj

dxi1 ... dxir

Suppose that system (18) admits a one-parameter group of equivalence transformations

x = x(x, u, c, a), u = u(x, u, c, a), c = c(c, a)

with the group parameter a. This group preserves the differential structure of system (18) and acts on the parameters c only. Since the transformed parameters x should be constant magnitudes, the corresponding transformation can only depend on the original parameters c (but not on x and u). The infinitesimal generator of the equivalence group can be written in the form

d

d

d

x = C(x,u,c)dx + (x,u,c)duj + T (c)del,

duj

where i = 1,... ,n, j = 1, .. .m, l = 1,.. . ,q. Further, we suppose that indexes i, j, l possess the above values and summation over a repeating index is assumed.

Statement 1. If system (18) admits the generator X, then this system also admits the generator FX with an arbitrary function F(c).

Proof. Suppose the condition of the statement is satisfied. Then the following relation holds

X E|e=o = 0, (19)

where X is the k-th prolongation of the generator X : k

x = e—+»j—+Ti— + vC - *

r=1 -1...-

The symbol |e-o means that the equality XE = 0 is satisfied on the manifold given by equations (18). Let us introduce the operator of full differentiation with respect to x®:

d ■ d ■ d

D® = dX® + ujdU- + ••• + + ■■■ ■

The coordinates of the prolongation generator are given by the formula [15]

= Dn ■ ■ ■ Dir (nj - ru) + e<..,ri. (20)

Consider a generator

FX = Ff -X- + F„j A + FT 1A dx® dwj dc1

The coordinates of the k-th prolongation (FX) are calculated by substituting F£® and Fnj

k

instead of and n- in formula (20). Then it is obvious that the relation (FX) = FX holds.

kk

Now, using (19), we have

(Ff )EIe=0 = F (X E )IE=0 = 0,

from which it follows that the generator FX is admitted by system (18). □

According to statement 1, a generator producing the group of equivalence transformations of constants is defined accurate to a factor arbitrarily depending on these constants. It should be noted that this fact was established earlier for the equivalence groups admitted by some equations, namely, Navier-Stokes equations [16] and modified Burgers equation [17]. The general proof is given here for the first time.

Now, we proceed to calculating the equivalence group of system (1). The corresponding infinitesimal generator is written as

X = X + na— + nN A + „^ A + + + A + „A + — + ng —,

' da 1 dN ' dpT ' dpc ' dx ' dD 1 dv ' dpo ' dg'

where X is given by (8). The coordinates of XX are assumed to depend on all independent and dependent variables as well as nine parameters entering the equations. These parameters are considered as additional variables that do not depend on t, x®, u®,p, T, C, i = 1, 2, 3. To ensure this, the conditions for zero derivatives of a, N, Pt, Pc, X, D, v, po, g with respect to the above variables are added to equations (1). The coordinates of the prolongation generator are calculated by the formulae obtained in [18]. Applying the prolongation generator to the system and making a transition to the corresponding manifold, we obtain the determining equations. It

follows from them that the coordinates of X that correspond to the parameters do not depend on t,xi,ui,p,T, C. The solution of the determining equations is written as

^ = (2c4 + v-1)t + co,

e1 = (c4 + VU v-1)x1 + cix2 + c2x3 + f 1(t),

e2 = -c1x1 + (c4 + v-1)x2 + csx3 + f 2(t), (21)

e3 = -c2x1 - csx2 + (c4 + v-1)x3 + f3(t), n1 = —c4u1 + c1u2 + c2u3 + ft(t), rj2 = —c1u1 — c4u2 + c3u3 + ft(t), n3 = —c2u1 — c3u2 — c4u3 + ft3(t),

rp = Po(c5gpTx3 + ceg^cx3 — ft\(t)x1 — f2 (t)x2 — f3 (t)x3)+ + (rP0 P-1 — 2c4)p + f 0(t),

rT = crT + cgC + c5, rC = cgC + c1oT + c6.

The coordinates that correspond to the parameters are given by

ra = (cg — c7)a + cga2 + c1o(xD-1 — l), (22)

arN = 2(c7 — cg)aN + (l — xD-1 — o2N)(c9 + c1oN), (23)

r^T = —(3c4 + c7 + rv v-1 + rgg-1)pT — c1o^c, (24)

= —(3c4 + cg + rv v-1 + rgg-1)Pc — cq^t , (25)

rx = xrv v-1 + (Cq — c1oN )ax, (26)

rD = D rvv-1 — (Cq — c1oN)aD. (27)

The following conditions must be satisfied:

[3t c2 = PT c3 = ^c c2 = Pc c3 = 0. (28)

In the presented formulae the quantities ci, i = l,..., 10, the functions fj(t), j = 0, l, 2, 3, and the coordinates rv, rp°, rg arbitrarily depend on the parameters a, N, Pt, Pc, X, D, v, po, g. It follows from statement 1 that without loss of generality we can consider ci as arbitrary constants and choose undefined coordinates as rv = cuv, rPo = c12po, ng = c13g, where en, c12, c13 are arbitrary constants. It is convenient to eliminate c4 from the expressions for r^T, by introducing new group constants X = 3c4 + c7, cg = 3c4 + cg.

It follows from (21) and (28) that the governing equations admit the generators X0, X12, Hi, Ho, Z3, U1, U2 (see (14)). If pT = pc = 0, then X13 and X23 are also admitted. Note that in case of pT = 0 (pc = 0) the generator T3 (C3) is admitted instead of U1 (U2). The above generators do not produce equivalence transformations of parameters. To obtain the generators producing equivalence transformations, we successively set one of the constants cy, cg, cq — c13 to unity and the remaining ones to zero. Note that if some parameters are equal to zero, then

the corresponding coordinates of the generator XX should be also set to zero. In what follows we consider the following three cases.

Table 2. Equivalence transformations.

Generator

Transformation

E1 TJ = s-1T, a = sa, NJ = s-2 N, Pt = spT

E2 (7 = s-1C, a = s-1a, NJ = s2N, pc = spc

E3 TJ = T + (1 - s)a-1 C, J = •s-1*^ DJ = sD,

a = s-1a, N = sN + (s - - 1)(XD-1 - s)a- 2,

pc = = Pc + (s - 1)a-1pT

Ea

C = C + aT, a = a — a + a*(D(1 + aNaj)

X = X(1 + aNa)--1

D = D(1 + aNa),

N = aDN(D(a - a)(1 + aNa) + ax) 1, pT = Pt - apC

E5 t = eat, x" = eax" , Pt = e-a Pt , pc = e-aPc,

J = - eaX, DJ = eaD, J = eav

Eg J = eap, Jo = eapo

Er Pt = eapT, pc = ea 'Pc, J = e-ag

1. Let a = 0 and N = 0, so the Soret and Dufour effects are allowed for. Then it follows from (22) - (27) that

na = (c8 - c7)a + cga2 + cw(xD-1 - 1), nN = 2(C7 - C8)N + (a-1 - x(aD)-1 - aN) (cg + ciqN),

(29)

n

Pt

-(cr + C11 + ci3)pT - cioPc,

nPC = -(C8 + C11 + C13)Pc - cqpt, nx = C11X + (cg - C10N)ax, 'qd

C11D - (cg - C10N)aD.

The admissible generators are written as

d d d d E1 = TdT - + 2NdN - pT ~3pr,

= CA + aA +

a dT da

D - x

a2D

N

d

dN

d d d d E = Cdc + - 2Ndm - pcdpC,

-PT— + x- - dA

a dpc dx dD

E = TdC +

X -1

d_

da

D - x aD

- aN

d d d d NdN - Pc WT - aNXdX + aNDdD,

1

д i д 2 d 3 д д д д д д

E = + x + x + x дхз - втдвТ -всдве + хдХ + +

„ д д т^ п д п д д

Еб = Рдр + Р0 дР0, E = вт двт + вс две - g д£ •

The corresponding equivalence transformations are given in Table 2. In the presented formulae a is a group parameter and s = e-a. Every transformation has its own independent group parameter. The omitted variables remain unchanged. In what follows the transformation corresponding to a generator will be denoted by the same symbol as this generator for simplicity.

Remark 1. The generator E3 is obtained by setting C9 = a-1 while the remaining group constants are zero. This choice allows to present the corresponding transformation in a convenient form.

Table 3. Additional equivalence transformations.

Generator Transformation

Eii Ei2 T = T = s-iT, а = sa, вт = ввт s-iT, вТ = «вт

E2i E22 (5 = C = s-iC, а = s-ia, вс = ввс s-iC, вС = «вс

Ei T = T + aC, вС = вс - авт

Ei E2 C = (5 = C + aT, а = а + a(xD-i — 1), вт = вт — авс C + aT, вт = вт — авс

Table 4. Admissible equivalence transformations.

N

Equivalence transformations x = D

= 0 = 0 Ei, , E3, E4, E5, E6, E7 = 0 0 El, E2, E41, E5, E6, E7

0 0 E2, E5, E6, E7 E3, E4

2. Let a = 0 and N = 0, so the Dufour effect is neglected. Then = 0, and from (23) it follows that (x — d)cg = 0. If x = D, we have cg = 0. Otherwise, the equality = , where and are given by (26) and (27), leads to cg = 0. The coordinates na, and become

na = (c8 — c7)a + cio(xD-1 — 1), nx = ciix, = ciiD,

while and are given by (29) with cg = 0. In this case the system admits the generators

д д д x д д д да - втдвт, E = сдс + ада - всдве'

Ei = т— - а— - вт

а

e4=Tdc+(xD-1 -1) da- pc dpT,

and E5, E6, Er from (30). The transformations generated by E}, E2, E4 are presented in Table 3.

3. Finally, when a = N = 0, the Soret and Dufour effects are not allowed for. Setting na = r/N = 0 in (22) and (23), we have cg(x - D) = c10(x - D) = 0. The admissible generators are

d d d d

e2 = TdT - Pt dpT, E2 = CdC - Pc Wc,

and E5, E6, Er from (30). If x = D, then the additional generators

d „ d _2 _ d „ d

E3 = CdT - Pt dpc, E = TdC - Pc dPT

are admitted. The transformations produced by E2, E2, E3, E4 are given in Table 3.

The three cases considered are summarized in Table 4, where the values of a and N are presented in the first and second columns and the admissible transformations are given in the third and fourth columns.

To complete the study of equivalence transformations, let us find the admissible reflections of parameters. Applying the scheme similar to the one described in Section 1, we find the following reflections

¿4 : t = -1, ul = -u", J = -X, D = -D, v = -v, i = 1, 2, 3;

d5 : pt = -p, Jo = -po; ¿6 : Pt = -Pt, Pc = -pc, tt = -g; (31)

dr : T = -T, Pt = - Pt , J = -a; ¿8 : C = - C, pc = - pc, J = -a

in addition to ¿1, ¿2, ¿3 from (17). Note that transformations ¿4 and ¿5 are physically meaningless.

Remark 2. The parameter s in transformations E1 and E2 is positive. Applying the reflections ¿r and ¿8, respectively, we can consider the above transformations with negative s also. The same is true for El, E2 and E2, E2.

3. Influence of Soret and Dufour Effects on Temperature and Concentration Fields

3.1. Transformation of the Governing Equations

As it was shown in the previous section, the governing equations admit an extensive equivalence group. In what follows, we show that with a help of this group system (1) can be transformed to the system with a = 0 and N = 0. In other words, a transformation connecting the equations with and without Soret and Dufour terms is derived.

Consider system (1) with a = 0 and N > 0. To vanish the Dufour terms in the heat equation, we apply the transformation E3 (see Table 2) with the parameter s that satisfies the quadratic equation

Ds2 - (a2DN + x + D)s + x = 0. (32)

Note that this parameter should be positive since s = e-a. When x > 0, D > 0, and N > 0, equation (32) always has two different positive roots except for the case N = 0, x = D, where these roots coincide (s = 1). As a result, we have

~ 1 — s T = T +-C,

Pc

s - 1

-Pt + Pc ,

X = s ^

D = sD,

a = s 1a,

NV = 0

(the other variables are unchanged). Now, we can use the transformation E4i (see Table 3) with the parameter a = aD(D — x)- to vanish the Soret term. It is admitted since a = 0 and N = 0. The transformed variables are written as

— 1 — s T = T +-C,

_ s _ 1

Pc = -Pt + Pc ,

- aDsT + (Ds - x) c — (x - Ds) pt + aDsPc C =-—-, pt = —

Ds2 - X

X - Ds2

X = s ^

D = sD,

= N = 0.

So, the successive action of transformations E3 and E4 makes a and N vanish. Finally, let us apply the transformations Ej2 and E| (see Table 3) with the parameters and

respectively (they are admitted since a = N = 0, see Table 4 and Remark 2). Then system (1) becomes

ut + (u ■ V)u = - — Vp + vAu - g(PTT' + PcC'), P0

Tt' + u ■ VT' = x 'AT',

Ct + u ■ VC' = D 'AC', divu = 0,

(33)

where

T' =

C' =

Pt(x - Ds) + PcaDs Pt(x - Ds2)

Pt (s - 1) + Pc a Pca(Ds2 - x)

t + !-s c

aDsT + (Ds - x)C

(34)

X ' = s 1 X,

D ' = sD.

(35)

So, systems (1) and (33) are connected by transformation (34), (35). Note that substituting the expressions for T', C', x ', D ' into (33) does not give the heat and mass transfer equations in the form (1). After changing the variables, it is also necessary to make two linear combinations of equations for T and C with appropriate coefficients (the systems connected by this transformation are, in fact, identical).

Transformation (34), (35) is invertible if the corresponding Jacobian

J

xD(Pt (X - Ds) + Pc aDs)(Pr (s - 1) + Pc a) aPTPc(X - Ds2)

is finite and non-zero. First of all, x — Ds2 = 0 when N > 0 (the case N = 0 will be considered below). It follows from the fact that equation (32) does not have the roots s = x/D for

positive N. Then J = 0 if and only if

, poa xpT

s = 1---— and s =

Pt ^ D(Pt - Po a)'

It can be shown that equation (32) does not have these roots if and only if N = No, where N0 is given by (13).

Since equation (32) has two different roots, there are two ways of transforming system (1) into system (33). The resulting equations will differ in the values of x' and D'.

If the Dufour effect is not allowed for (N = 0), then equation (32) has the roots s =1 and s = x/D. In case s = 1, formulae (34) become

' Pt(x - D) + PoaD T aD T + C

T =-TT~!-TT\-T, C = n-T + C, (36)

Pt(x - D) D - x

while x' = x and D ' = D. This transformation is invertible if and only if x = D and a =

Pt(D - x)(PoD)-1.

The inverse transformation for (34), (35) has the form

= Pt(x 's - D ') + Po (s - 1) '

Pt (x's - D') + Po aD' + Pt (s - 1) + Po a ' ( )

C =_ptaD'_T' + Po a C'

Pt (x 's - D ')+ po aD ' Pt (s - 1) + Po a '

x = sx', D = s-1D '. (38)

Here s is a root of the equation

x's2 - (x' + D ')s + D '(a2N + 1) = 0, (39)

which is obtained from (32) by substituting the expressions for x and D from (38). This equation has two positive roots if

a2N < (x' - D ')2 . (40)

- 4x'D' y '

Otherwise, the roots are complex numbers. The above transformation is invertible if and only if

' = D ' =, Po a = D '(Pt - Po a)

x = D , s = 1--, s = -.

x ^ Pt , ^ Ptx'

In case N = 0 equation (39) has the roots s = 1 and s = D '/x'. If we put s = 1, then formulae (37) become

T = pT(x- D) T' C = pTaD T' + C' (41)

Pt (x - D) + Po aD , Pt (x - D) + Po aD + , ()

and the inverse transformation for (36) is found.

The obtained formulae can be used for transforming solutions of system (1) into solutions of system (33) and vice versa. Suppose that the functions u, p, T', C' depending on x' and D' provide a solution of equations (33) (the dependence on the other parameters is not important here). If we introduce new variables T and C by formulae (37) and express x' and D' in terms

of x and D according to (35), then the functions u, p, T, C will provide a solution of equations (1). Note that s should be chosen as a root of equation (32).

Now, let the functions u, p, T, C depending on x, D, a, N form a solution of system (1). The corresponding solution of system (33) can be easily obtained by setting a = N = 0. On the other hand, we can introduce new variables T' and C' by formulae (34) and transform the parameters according to (38). In this case, s is a root of equation (39). The functions so obtained provide a solution to system (33), but it is defined for all x' and D' satisfying condition (40).

If the Dufour effect is neglected (N = 0), then the solutions are transformed by formulae (36) and (41).

3.2. Dimensionless Formulation

In practice, it is convenient to present a solution of a concrete problem in dimensionless form. Non-dimensional equations and boundary conditions are also used in numerical simulation. In this subsection, we derive the formulae connecting systems (1) and (33) in dimensionless variables. Using these formulae, we will show that an initial and boundary value problem for Soret-Dufour equations can be reduced to a problem for the system without cross-effects.

Let us introduce the scales of length L, time L2/v, velocity v/L, pressure p0v2/L2, temperature ©, and concentration PtO/Pc. Non-dimensional governing equations have the form

ut + (u • V)u = —Vp + Au + Gr(T + C)e,

Tt + u •VT = Pr-iAT — ^QSc-i(AC — ^AT), (42)

Ct + u • VC = Sc- (AC — ^AT), divu = 0,

where e = (0,0,1). The system contains five non-dimensional parameters: the Grashof number Gr = g^TOh3/v2, the Prandtl number Pr = v/x, the Schmidt number Sc = v/D, the separation ratio ^ = —aPc/Pr, which characterizes the Soret effect, and the Dufour number Q = NpT/PC > 0. The cases ^ > 0 and ^ < 0 correspond to positive and negative Soret effect, respectively.

Equations (33), where the Soret and Dufour effects are not allowed for, in dimensionless form are written as

ut + (u • V)u = —Vp + Au + Gr(T' + C')e,

T/ + u • VT' = Pr'— AT', (43)

Ct + u •VC' = Sc'-iAC', divu = 0,

where Pr' = v/x' and Sc' = v/D'.

The transformation connecting systems (42) and (43) is readily obtained from (34), (35) by a formal change of parameters

PT ^ Gr, PC ^ Gr, x ^ Pr_i, D ^ Sc_i,

a ^ —N ^ Q, x' ^ Pr'_i, D' ^ Sc'_i.

As a result, we have

T ' = C ' =

Sc - sPr(^ + 1)

Sc — s2Pr

^ + 1 — s

T + i-! C

(44)

^(Sc — s2 Pr)

s^ Pr T + (Sc — sPr) C

Pr' = sPr, Sc' = s 1Sc. where s is a root of the quadratic equation

Pr s2 — (^2Q Pr + Pr + Sc)s + Sc = 0.

(45)

(46)

Equation (46) always has two different positive roots except for the case Q = 0, Pr = Sc, in which these roots are equal (s = 1). We suppose that ^ = 0 (otherwise the above transformation is identical). So, there are two ways of transforming system (42) into system (43). The resulting equations will differ in the values of Pr' and Sc'.

The inverse transformation exists if and only if

Pr(^ + 1) - Sc

Q=

Pr ^ +1)

Note that Sc - s2Pr = 0 since equation (46) does not have the roots s = Sc/Pr for positive Q.

If the Dufour effect is not allowed for (Q = 0), then the roots are s = 1 and s = Sc/Pr. In the former case formulae (44) become

T' = Sc — Pr(^ + 1) T

C ' =J^ T + C,

Sc — Pr Sc — Pr

where the conditions Sc = Pr and Sc = Pr(^ + 1) must be satisfied. The inverse transformation for (44), (45) is written as

sSc' — Pr' . T' + s — 1 C'

(47)

T

C

sSc' — Pr'(^ + 1) s — 1 — ^ ' ^Pr' T' + ^ C'

Pr'(^ +1) — sSc' ^ +1 — s ,

(48)

Pr = s-1Pr', Sc = sSc', (49)

where s is a root of the equation

Sc's2 - (Pr' + Sc')s + Pr'(^2Q + 1) = 0. (50)

This equation has two positive roots if

*2Q < ^.

Otherwise, the roots are complex numbers. So, system (43) can be transformed into two systems (42) with different values of Prandtl and Schmidt numbers. The transformation is invertible if and only if

Pr'(^ + 1)

Pr' = Sc', s = ^ + 1, s =

Sc'

In case of Q = 0, equation (50) has the roots s = 1 and s = Pr'/Sc'. Taking s = 1, we obtain the inverse transformation for (47):

T = Sc — Pr(- + 1) T', C = Pr(- + 1) — Sc T' + C'. (51)

It should be noted that a linear change of temperature and concentration, which vanishes the Soret and Dufour terms, was found in [19] without using symmetry methods. The momentum equation was written in terms of vorticity, so the pressure term was excluded. The transformation also involved changing time, velocity, and several dimensionless parameters (Prandtl number, Lewis number, and temperature and concentration Rayleigh numbers). In this paper, the equivalence properties of governing equations allowed us to minimize the number of transformed variables to temperature, concentration, thermal diffusivity, and diffusion coefficient (or the Prandtl and Schmidt numbers in dimensionless formulation). In pure Soret case, only linear change of temperature and concentration is required. A detailed study of conditions, under which the derived transformation is invertible, has been also performed.

3.3. On Initial and Boundary Value Problems for Soret—Dufour Equations

The formulae obtained in previous section can be used for reducing an initial and boundary value problem for Soret-Dufour equations (42) to a problem for system (43), which has a simpler form due to the absence of Soret and Dufour terms. The following scheme of three steps is proposed:

1. Boundary and initial conditions for system (42) are transformed to those for system (43) by using formulae (45), (48) or (51) (in pure Soret case).

2. The obtained problem is solved numerically or, if possible, analytically. Note that if the Dufour effect is present (Q = 0), then the values of Pr' and Sc' are required. They can be found from (45).

3. To obtain the solution of the original problem, the temperature and concentration fields are transformed according to (45), (48) or (51) (if there is no Dufour effect).

Let us consider these steps in more detail. Suppose that we have a problem on finding the functions u, p, T, C that should satisfy equations (42) with — = 0, Q > 0 and some initial and boundary conditions. These conditions can be considered as relations on the unknown functions and their derivatives that should be satisfied on some manifolds in the space of t, x:

Fj (u,p, T, C) =0, j = 1,...,n. (52)

jV ' Q, (t,x)=0 V '

In general, Fj is a differential operator that acts on the unknown functions.

The described problem can be reduced to the problem on finding the functions u, p, T', C' that satisfy equations (43) and the transformed imposed conditions. The latter are obtained from (52) by expressing T and C in terms of T' and C' according to the formulae (see (45), (48))

TVT' n' ) — s*Sc (si) — Pr (si) T" i s* — 1 n' (K?\

T(T,C,sij = si Sc'(si) — Pr'(si) (- +1) T + si — 1 — - C, (53)

C(T', C', s4) = „ , / Pr'(si)-— T' + —t— C',

y , , tJ Pr'(si) (t + 1) - si Sc'(si) t + 1 - si ,

Pr'(si) = siPr, Sc' (si) = s-1Sc. (54)

Here si is a root of the quadratic equation (46), i = 1, 2. For t = 0, Q > 0, and positive numbers Pr and Sc, this equation always has two different positive roots except for the case Q = 0, Pr = Sc, where si = 1. In this case, transformation (53) is not invertible, so these values of parameters are not considered below.

We see that the original initial and boundary value problem for equations (42) can be reduced to two problems for system (43) that correspond to different roots of the quadratic equation. Let u(si), p(si), T'(si), C'(si) be a solution of the problem

ut + (u • V)u = -Vp + Au + Gr(T' + C')e,

Tt' + u • VT' = Pr'(si)-1AT',

Ct' + u •VC' = Sc'(si)-1AC', (55)

divu = 0,

Fj( u, p, T (T ',C',si), C (T ',C ',si)) where T, C, Pr', Sc' are given by (53), (54).

= 0,

fi, (t,x)=0

Statement 2. The solutions of problem (55) corresponding to the roots si and s2 are connected by the following relations:

u(si) = u(s2), p(si)= p(s2), T'(si) = C'(s2), C' (si) = T'(s2). (56)

Proof. Consider problem (55) in case of i = 1. Taking into account that sis2 = Sc/Pr (see (46)) and using (53), (54), we have

( , , i) Sc — siPri' + 1) + si — 1 — ' s2 — 1 T ' + Sc — s2Pr C ' =

s2 — 1 — ' Sc — s2Pr(' + 1) s2 Sc'(s2) — Pr'(s2) C' + s2 — 1 t' = T(C', T', s2). (57)

s2 Sc'(s2) - Pr'(s2)(t + 1) s2 - 1 - t

The relation C(T', C',s1) = C(C',T',s2) can be established in the same way. Now, using the equalities Pr'(s1) = s1Pr = s-1Sc = Sc'(s2), Sc'(s1) = Pr'(s2), and the relations obtained above, we can rewrite problem (55) for i = 1 in the form

ut + (u • V)u = -Vp + Au + Gr(C' + T')e,

C' + u • VC' = Pr'(s2)-1AC',

Tt' + u •VT' = Sc'(s2)-1AT', (58)

divu = 0,

Fj ( u, p, T (C ',T '.sa), C (C',T ',s2 ) )

= 0.

Q (t,x) = 0

Here the second and third equations from (55) are written in a reversed order. Comparing problems (58) and (55) for i = 2, we see that the latter can be obtained from the former by replacing T' with C' and vice versa. It follows that the solutions of these problems have the same relation, so equalities (56) hold true. □

As soon as the solution of problem (55) is found, the solution of the original problem for Soret-Dufour equations is given by formulae (53). It should not depend on the choice of the root si. The equality of velocity vectors and pressures corresponding to different roots directly follows from Statement 2. Using (56) and (57), we have for temperature fields

T(T'(si), C'(si), si) = T(C'(si), T'(si),s2) = T(T'(s2),C'M,s2).

The equality of concentrations corresponding to the roots si and s2 can be established in the same way.

Note that we can also use dimensional equations and boundary conditions when reducing a problem for Soret-Dufour equation to a problem for the system without these effects (see section 3.1). In this case, it is easy to show that the resulting solution does not depend on the choice of the root as well.

Here we do not discuss the variety of initial and boundary conditions on the velocity, pressure, temperature, and concentration fields (see [2-4,20]). However, it would be useful to consider an example of such condition and to see how it is changed under the transformations derived above. Suppose that the Soret effect is taken into account. Then the mass diffusion flux is driven by concentration and thermal gradients:

j = —po(DVC + aDVT).

It is often assumed that the mass diffusion flux through a rigid or free boundary r is zero: j • n|r= 0, where n is the unit normal vector to r. In dimensionless variables introduced in the previous section, this relation can be written as

dC ^dT dn dn

= 0.

r

The corresponding condition for system (43) is obtained by applying formulae (48):

= 0.

dC' Sc'(^ + 1 - s) dT'

I dn Pr'(^ + 1) - sSc' dn

The other possible initial or boundary conditions for system (42) can be easily transformed to those for system (43) by applying (48), (49) or (51) (in pure Soret case).

Conclusion

The equations describing convection in binary mixture with Soret and Dufour effects are considered. The symmetry classification of this system with respect to the nine constant parameters is made. The one-parameter transformation subgroups produced by the generators of the admissible

Lie symmetry algebras are described. It is shown that a generator producing equivalence transformations of constants is defined accurate to a factor arbitrarily depending on these constants. The equivalence group admitted by the governing equations is calculated. This group is used to derive a transformation connecting the systems with and without Soret and Dufour terms. In pure Soret case, it reduces to a linear change of temperature and concentration. The presence of Dufour effect requires an additional change of thermal diffusivity and diffusion coefficient. It is shown that there are two ways of transforming the systems into each other. They correspond to different values of parameter, which solves the quadratic equation. The transformation is presented in dimensional and non-dimensional forms. The conditions, under which the transformation is invertible, are investigated. A scheme for reducing an initial and boundary value problem for Soret-Dufour equations to a problem for the system without cross-effects is proposed. These results are useful for analytical and numerical modelling of convection in binary systems with Soret and Dufour effects.

This work is supported by the INTAS YS Project 06-1000014-6257 and the Russian Foundation for Basic Research Project 08-01-00762-a. The author is grateful to Prof. Victor K. Andreev for helpful comments and discussion.

References

[1] L.V.Ovsyannikov, Group Analysis of Differential Equations, New York, Academic Press, 1982.

[2] G.Z.Gershuni, E.M.Zhukhovitskii, Convective Stability of Incompressible Fluids, Jerusalem, Keter, 1976.

[3] G.Z.Gershuni, E.M.Zhukhovitskii, A.A.Nepomnyashchy, Stability of Convective Flows, Moscow, Nauka, 1989 (in Russian).

[4] J.K.Platten, J.C. Legros, Convection in Liquids, Berlin, Springer, 1984.

[5] L.D.Landau, E.M. Lifshitz, Fluid Mechanics, Oxford, Pergamon Press, 1987.

[6] O.N.Goncharova, Group Classification of the Free Convection Equations, Continuum dynamics: collection of papers, Novosibirsk, 79(1987), 22-35 (in Russian).

[7] V.K.Andreev, O.V.Kaptsov, V.V.Pukhnachov, A.A.Rodionov, Applications of Group-Theoretical Methods in Hydrodynamics, Kluwer Academic Publishers, 1998.

[8] V.L.Katkov, Exact Solutions of Certain Convection Problems, J.Appl. Math. and Mech., 32(1968), no. 3, 489-495.

[9] V.K.Andreev, I.I.Ryzhkov, Symmetry Classification and Exact Solutions of the Thermal Diffusion Equations, Differential equations, 41(2005), no. 4, 538-547.

[10] I.I.Ryzhkov. On the Normalizers of Subalgebras in an Infinite Lie Algebra, Communications in Nonlinear Science and Numerical Simulation, 11(2006), no. 2, 172-185.

[11] I.I.Ryzhkov, Symmetry Analysis of the Thermal Diffusion Equations in the Planar Case, Proceedings of 10th International Conference on MOdern GRoup ANalysis, Larnaca, Cyprus, (2005),182-189.

12] I.I.Ryzhkov, Invariant Solutions of the Thermal Diffusion Equations for a Binary Mixture in the Case of Plane Motion, J. Appl. Mech. Tech. Phys, 47(2006), no. 1, 79-90.

13] I.I.Ryzhkov, On Double Diffusive Convection with Soret Effect in a Vertical Layer between Co-Axial Cylinders, Physica D: Nonlinear phenomena, 215(2006), 191-200.

14] A.F.Sidorov, V.P.Shapeev, V.P.Yanenko, The Method of Differential Connections and its Applications to Gas Dynamics, Novosibirsk, Nauka, 1984 (in Russian).

15] N.H.Ibragimov, Transformation Groups Applied to Mathematical Physics, Dordrecht, Rei-del, 1985.

16] N.H.Ibragimov, G.Unal, Lie Groups in Turbulence, Lie Groups and their Applications, 1(1994), no. 2, 98-103.

17] V.F.Kovalev, V.V.Pustovalov, Lie Algebra of Renormalization Group Admitted by Initial Value Problem for Burgers Equation, Lie Groups and their Applications, 1(1994), no. 2, 104-120.

18] S.V.Meleshko, Generalization of the Equivalence Transformations, Nonlinear mathematical physics, 3(1996), no. 1-2, 170-174.

19] E.Knobloch, Convection in Binary Fluids, Physics of Fluids, 23(1980), no. 9, 1918-1920.

20] S.Hollinger, M.Liike, Influence of the Dufour Effect on Convection in Binary Gas Mixtures, Physical review E, 52(1980), no. 1, 642-656.