- Email: [email protected]
A note on the solution of a third-order eigenvalue problem arising in convective instability
A Note on the Solution bf a Third-order Eigenualue Problem arising in Convective Instability by
P.
G. DANIELS
Department of Mathematics, London EC1 V OHB, U.K. and
City University,
Northampton
Square,
M. WEINSTEIN
Ministry
of Defence, Armament
Development
Authority,
Haifa, Israel
ABSTRACT: Instability of a certain flow between heated vertical planes arises when a convective parameter y exceeds a critical value yC obtainedfrom the solution of a third-order eigenvalue problem. This problem involves a logarithmic singularity at the origin and an accurate determination of yC is achieved by a combination of series expansion and numerical integration techniques.
I. Zntroduction
An investigation of the instability of high Prandtl number flow between heated vertical planes has led to a simple criterion for the existence of stationary multiplecell convection : a convective parameter y based on the stratification and Rayleigh number of the flow must be greater than a critical value yCdetermined as the lowest eigenvalue of the system (1,2) d3f ---4x;y)f dx3 -f+o
=o
(I.11
(-;,
$0
(x=0).
(1.2)
Here CO= (dO/dx)/(d$/dx) where the basic non-dimensional horizontal temperature gradient de/dx and vertical fluid velocity - d$/dx are known in terms of elementary functions, with (1.3)
0 The Franklin Institute OOlCCIl32/92 $5.00+0.00
17
P. G. Daniels and M. Weinstein where do (X; y) = -4(D+ dY =(X;y)
cash Xcos X+ D_ sinh Xsin X),
(1.4)
= (D++D_)sinhXcosX+(D_-D+)coshXsinX,
X=yx
D,
and
(1.5)
= - (cash $y sin $y f sinh iy cos :y)/8d
(1.6)
where d = sinh2:y+sin2:y. These formulae constitute an exact solution of the Navier-Stokes-Boussinesq equations first discussed by Elder (3). Properties of dO/dXand dY/dXindicate that yCmust lie in the range ya < yC < y,, where ya z 4.73 is the first non-trivial zero of {tanh ty+ tan :y} and y,, z 7.85 is the first non-trivial zero of {tanh ty - tan $y}, these being the points at which the horizontal temperature gradient d@/dx and vertical shear d2$/dx2 reverse sign at x = 0, mid-way between the two planes x = k 4. Numerical solutions of (I. 1), (1.2) using a fourth-order Runge-Kutta scheme to within one step of the singularity in o at x = 0 have led to a prediction that yC z 6.30 (l), and while it is possible to conjecture that there may be an exact solution yC = 2n, this has not been established. Such a conjecture is suggested partly on physical grounds, being the point at which the centre-line stream function
(1.7) reaches a maximum as a function of y. However, the present work describes a refined numerical approach based on complete series expansions of the solution forfabout both x = - $ and x = 0 which confirms that yCis not 271 and is in fact given to six significant figures by yC = 6.29829.
II. Series Expansions of the Eigensolution It is convenient
to set f(x)
= F(X) so that
d3F --Q(X;y)F=O dX3 F=$
(x=
where R = (dO/dX)/(dY/dX) boundary conditions
(IT. 1)
$0
_ty);
(11.2)
(X=0)
and to note that 0 and Y satisfy the equations
d4Y ~ = d@ _.._ dX4 dX’ 78
(-+ydxdo)
d20 dX2=
dY
and
(11.3)
-4dx JCMd
of the Franklin Institute Pergamon Press plc
Third-order Eigenvalue dY Y =dx=o
o= *:,
Problem in Convective
(X=
i-:y>
Instability
(11.4)
with solutions given by (1.4) and (1.5) above. Although F = dY/dX is an exact solution of (II. 1) it does not satisfy dF/dX = 0 at either X = -47 or X = 0 and the relevant solutions of (II. 1) are found instead by constructing series expansions about these points. For the base functions dO/dX and dY’/dX, expansions about the origin proceed in even and odd powers of X respectively, with successive terms being determined from (11.3) to give
(11.5)
where a,, = -4D+,
b, = 2D_
and a, = -2b,_,/n b, = a,- ,/(2n+
1)2n(2n-
1)
(n = 1,2...).
(11.6)
Two independent series for F can then be constructed satisfying the requirement dF/dX = 0 at X = 0. The first of these F = F,, say, proceeds in even powers of X : F, = f
p,X*“.
(11.7)
The constant p, is arbitrary and substitution into (11.1) and use of (11.5) leads to the determination of the remaining coefficients pn from the formulae 92 = aopilbo,
P2 =
42124
n-2 qn
=
c
i= 0
QiP,_i& I - lg
h,y.,)/b,
\
(11.8)
pn = q,/2n(2n - 1) (2n - 2) The second solution F = F2, say, again proceeds additional set of logarithmic terms :
(n = 3,4...).
in even powers of X but with an
F2 = ,f {r,X2”+~,X2nln(-X)}.
(11.9)
?Z=O
Here r0 is arbitrary and the choice r, = 0 excludes the previously defined solution for F,. The remaining coefficients r,,, s, are determined from (11.1) and (11.5) as Vol. 329, No. I, PP. 77-83, Pnnted in Great Britain
1992
79
P. G. Daniels and M. Weinstein so = 0,
(alro-2b~sdlbo,
t2 =
(yg: a;s,-i-
u, =
t,= s,
=
roaoPb0,
SI =
I-
s2
u2 =
u2/24,
=
wllbo,
r2 = (t2-26s2)/24,
b;u.,>/bo,
y<
(n = 3,4..
~~‘air~_j~,-n~2bit~_,-2s,b~~,)i,,,, ,=O
;=
.).
(11.10)
I
u,/2n(2n - 1)(2n - 2),
r, = (t,-(12n2-
12n+2)sn)/2n(2n-
1) (2n-2)
It is of interest to note that the logarithmic part of F2 is proportional to F, (specifically (roao/2bop,)FI(X) In (-X)), a fact which can be established from the original differential equation as follows. Set F = F,(X) U(X) so that U satisfies d3U dF, d2U F1d~+3~~+3~~=0 and let dU/dX = U. be the solution Further, write dU/dX = U,(X)V(X) to
d2F, dU
(II. 11)
for which F - kX as X -+ 0 (i.e. U,, - koX-2). to obtain the third solution, corresponding
(11.12) Since U. and F, do not contain logarithmic terms in their expansions about X = 0 and Uo, and therefore U; ‘F; 3, contain no terms proportional to XP ‘, the expansion for V about X = 0 does not contain logarithmic terms. However, U. V can contain a term proportional to X- ’ in which case X
U. VdX
u=
(11.13)
s contains a single logarithmic term, K In (-X), in its expansion about X = 0. Thus, if the general solution for F contains a logarithmic term, it must be of the form KF,(X) In (-X). Confirmation of this follows from the observation that in (II. 10) and (11.8), s, = (aoro/2bo)(pn/p,), n = 1, 2. . Expansions for the base functions d@/dX and dY/dX about X = - fv are d@ -= dX
W c
(II. 14)
&nXn,
n-0
where 8 = X+ $y and Z. = -4(D+coshfycos$-D_ 6, = 2(D_coshiycos’,y-D, with the remaining 80
coefficients
sinh$ysin’,y), sinhjysinfy),
then determined
G, = 0, h; = -a,
(II. 15)
from (11.3) by the relations Journal
of the Franklin Pergamon
lnst~tute Press plc
Third-order Eigenvalue
Problem in Convective
ci,= 46_ ,/n, = Zn_J(n+
L, The series expansion where
for Fabout
l)n(n-
(n = 2,3. . .).
1) 1
J? = 0 satisfying
F = dF/dX
Instability
(II. 16)
= 0 there is F = F,(X)
F3 = f j@, n=O
(II. 17)
p”,, = p”, = 0 and fi2 is arbitrary. Substitution of (II. 17) and (11.14) into (II. 1) leads to the determination of the remaining coefficients j?n from the formulae 40 = 0,
4”, = p”&ib”,,
B+ I = kAn+
lM-
11,
I (11.18)
ZZZ.Numerical Method The solution of the singular eigenvalue problem (II. l)-(11.2) is now constructed in three regions, one where the series expansion (II. 17) is valid, giving F= AF3(X) another
where the expansions
(-4~
d Xb
X0),
(III. 1)
(11.7) and (11.9) are valid, giving
F= BF,(X)+CF,(X)
(X, < X< O),
(111.2)
and a third (X0 d X < X,) where the solution F3 is continued from X0 to X, using a fourth-order Runge-Kutta integration of the original differential equation. It is assumed that F,, F2 and F3 are uniquely fixed by the specification p , = r. = jj2 = 1. At X = X0 the values of F3, dFJdX and d*FJdX* determined by (11.17) provide the initial conditions for the Runge-Kutta computation and at X, a consistent overall solution is obtained provided that F, dFjdX and d2F/dX2 are continuous, requiring
(111.3) to vanish at X = X,. This must be achieved by suitable choice of the parameter y. It is of interest to note that the determinant D is a function of y but not of X,, the location at which it is evaluated. This is easily established using the fact that each F, (i = 1, 2, 3) satisfies Vol. 329, No. I, pp. 77-83, 1992 Printed in Great Britain
81
3E gg
$71 $2
2 : E ,g 4
40 40 40
400 400 400
M
N
0.49 0.375 0.3125
XOIY
0.49 0.49 0.375 0.375 0.375 0.375 0.3125 0.3125
200 200 200 200 400 400 400 400
5 20 20 40 20 40 20 40
-
- XII/Y
M
N
0.01 0.125 0.1875
-X,/Y
0.01 0.01 0.125 0.125 0.125 0.125 0.1875 0.1875
--X,/Y
TABLE I
-0.000105 -0.000126 -0.000125
D(6.298284)
Computations
- 0.932430 - 0.932430 -0.932425 -0.932425 -0.932425 -0.932425 - 0.932425 -0.932426
D(6.28) 0.088974 0.088974 0.088982 0.088982 0.088979 0.088979 0.088983 0.088982
D(6.30)
- 0.000094 - 0.000069 -0.000072
D(6.298285)
- 0.000008 - 0.000022 -0.000021
D(6.298286)
ofD near y = 6.298286
TABLE II
- 0.426388 -0.426388 -0.426385 -0.426385 - 0.426382 -0.426383 - 0.426384 -0.426385
D(6.29)
Computations of D war y = 6.30
0.000003 0.000030 0.000030
D(6.298287)
1.148328
0.613829 0.613829 0.613834 0.613834 0.613834 0.613834 0.613836 0.613835
0.000070 0.000086 0.000084
D(6.298288)
1.148325 1.148329 1.148327
1.148325
1.148326 1.148326
1.148328
D(6.32)
D(6.31)
P
!?
;f. G-
b g
a
Third-order Eigenvalue d’F, ~ -QFi dX3
Problem in Convective
Instability
= 0.
Multiplication of the equation for F, by {F3 dF2/dX- F2 dF,/dX} and the equations for F2 and F3 by similar expressions in cyclic rotation, integration of the resulting equations by parts and addition shows that D as defined by (111.3) is independent of x.
IV. Results
and Conclusions
The recurrence formulae defined in Section II were used to generate the functions F,, F2 and F3 on a computer up to terms in X2N and 2’” at the locations X = X0 and X = Xi. The Runge-Kutta integration was performed using A4 equal intervals of X in the region X0 d X < Xi. Tests were carried out to investigate the effects of each of the parameters N, X,,, Xi and M and some of the results are summarized in Table I. The constancy of the determinant D for different values of X, also provided a useful check. At low values of y the two series expansions provide excellent approximations to Facross the entire domain but as y increases above 5 the two radii of convergence in x decrease until when y = yC z 6.3 they appear to not quite overlap near X = -:y. Some computations were performed with X0 = X, = -:y and produced quite accurate results for yC at very high levels of truncation (N up to 90) but it became evident that the incorporation of the RungeeKutta region was necessary to produce definitive results. From these (Table II) it was concluded that the first zero of D occurs at y = yC = 6.29829 and no further zeros of D were identified in the range y d 7.7. Acknowledgement This work was partially supported by a Visiting Science and Engineering Research Council.
Fellowship
Research
Grant
from the
References (1) P. G. Daniels, “Convection in a vertical slot”, J. Fluid Mech., Vol. 176, pp. 419441, 1987. (2) P. G. Daniels, “Stationary instability of the convective flow between differentially heated vertical planes”, J. Fluid Mech., Vol. 203, pp. 525-540, 1989. (3) J. W. Elder, “Laminar free convection in a vertical slot”, J. Fluid Mech., Vol. 23, pp. 77-98, 1965.
Vol. 329, No. 1. pp. 77-83, Pnnted in Great Britain
1992
83
P.
G. DANIELS
Department of Mathematics, London EC1 V OHB, U.K. and
City University,
Northampton
Square,
M. WEINSTEIN
Ministry
of Defence, Armament
Development
Authority,
Haifa, Israel
ABSTRACT: Instability of a certain flow between heated vertical planes arises when a convective parameter y exceeds a critical value yC obtainedfrom the solution of a third-order eigenvalue problem. This problem involves a logarithmic singularity at the origin and an accurate determination of yC is achieved by a combination of series expansion and numerical integration techniques.
I. Zntroduction
An investigation of the instability of high Prandtl number flow between heated vertical planes has led to a simple criterion for the existence of stationary multiplecell convection : a convective parameter y based on the stratification and Rayleigh number of the flow must be greater than a critical value yCdetermined as the lowest eigenvalue of the system (1,2) d3f ---4x;y)f dx3 -f+o
=o
(I.11
(-;
$0
(x=0).
(1.2)
Here CO= (dO/dx)/(d$/dx) where the basic non-dimensional horizontal temperature gradient de/dx and vertical fluid velocity - d$/dx are known in terms of elementary functions, with (1.3)
0 The Franklin Institute OOlCCIl32/92 $5.00+0.00
17
P. G. Daniels and M. Weinstein where do (X; y) = -4(D+ dY =(X;y)
cash Xcos X+ D_ sinh Xsin X),
(1.4)
= (D++D_)sinhXcosX+(D_-D+)coshXsinX,
X=yx
D,
and
(1.5)
= - (cash $y sin $y f sinh iy cos :y)/8d
(1.6)
where d = sinh2:y+sin2:y. These formulae constitute an exact solution of the Navier-Stokes-Boussinesq equations first discussed by Elder (3). Properties of dO/dXand dY/dXindicate that yCmust lie in the range ya < yC < y,, where ya z 4.73 is the first non-trivial zero of {tanh ty+ tan :y} and y,, z 7.85 is the first non-trivial zero of {tanh ty - tan $y}, these being the points at which the horizontal temperature gradient d@/dx and vertical shear d2$/dx2 reverse sign at x = 0, mid-way between the two planes x = k 4. Numerical solutions of (I. 1), (1.2) using a fourth-order Runge-Kutta scheme to within one step of the singularity in o at x = 0 have led to a prediction that yC z 6.30 (l), and while it is possible to conjecture that there may be an exact solution yC = 2n, this has not been established. Such a conjecture is suggested partly on physical grounds, being the point at which the centre-line stream function
(1.7) reaches a maximum as a function of y. However, the present work describes a refined numerical approach based on complete series expansions of the solution forfabout both x = - $ and x = 0 which confirms that yCis not 271 and is in fact given to six significant figures by yC = 6.29829.
II. Series Expansions of the Eigensolution It is convenient
to set f(x)
= F(X) so that
d3F --Q(X;y)F=O dX3 F=$
(x=
where R = (dO/dX)/(dY/dX) boundary conditions
(IT. 1)
$0
_ty);
(11.2)
(X=0)
and to note that 0 and Y satisfy the equations
d4Y ~ = d@ _.._ dX4 dX’ 78
(-+ydxdo)
d20 dX2=
dY
and
(11.3)
-4dx JCMd
of the Franklin Institute Pergamon Press plc
Third-order Eigenvalue dY Y =dx=o
o= *:,
Problem in Convective
(X=
i-:y>
Instability
(11.4)
with solutions given by (1.4) and (1.5) above. Although F = dY/dX is an exact solution of (II. 1) it does not satisfy dF/dX = 0 at either X = -47 or X = 0 and the relevant solutions of (II. 1) are found instead by constructing series expansions about these points. For the base functions dO/dX and dY’/dX, expansions about the origin proceed in even and odd powers of X respectively, with successive terms being determined from (11.3) to give
(11.5)
where a,, = -4D+,
b, = 2D_
and a, = -2b,_,/n b, = a,- ,/(2n+
1)2n(2n-
1)
(n = 1,2...).
(11.6)
Two independent series for F can then be constructed satisfying the requirement dF/dX = 0 at X = 0. The first of these F = F,, say, proceeds in even powers of X : F, = f
p,X*“.
(11.7)
The constant p, is arbitrary and substitution into (11.1) and use of (11.5) leads to the determination of the remaining coefficients pn from the formulae 92 = aopilbo,
P2 =
42124
n-2 qn
=
c
i= 0
QiP,_i& I - lg
h,y.,)/b,
\
(11.8)
pn = q,/2n(2n - 1) (2n - 2) The second solution F = F2, say, again proceeds additional set of logarithmic terms :
(n = 3,4...).
in even powers of X but with an
F2 = ,f {r,X2”+~,X2nln(-X)}.
(11.9)
?Z=O
Here r0 is arbitrary and the choice r, = 0 excludes the previously defined solution for F,. The remaining coefficients r,,, s, are determined from (11.1) and (11.5) as Vol. 329, No. I, PP. 77-83, Pnnted in Great Britain
1992
79
P. G. Daniels and M. Weinstein so = 0,
(alro-2b~sdlbo,
t2 =
(yg: a;s,-i-
u, =
t,= s,
=
roaoPb0,
SI =
I-
s2
u2 =
u2/24,
=
wllbo,
r2 = (t2-26s2)/24,
b;u.,>/bo,
y<
(n = 3,4..
~~‘air~_j~,-n~2bit~_,-2s,b~~,)i,,,, ,=O
;=
.).
(11.10)
I
u,/2n(2n - 1)(2n - 2),
r, = (t,-(12n2-
12n+2)sn)/2n(2n-
1) (2n-2)
It is of interest to note that the logarithmic part of F2 is proportional to F, (specifically (roao/2bop,)FI(X) In (-X)), a fact which can be established from the original differential equation as follows. Set F = F,(X) U(X) so that U satisfies d3U dF, d2U F1d~+3~~+3~~=0 and let dU/dX = U. be the solution Further, write dU/dX = U,(X)V(X) to
d2F, dU
(II. 11)
for which F - kX as X -+ 0 (i.e. U,, - koX-2). to obtain the third solution, corresponding
(11.12) Since U. and F, do not contain logarithmic terms in their expansions about X = 0 and Uo, and therefore U; ‘F; 3, contain no terms proportional to XP ‘, the expansion for V about X = 0 does not contain logarithmic terms. However, U. V can contain a term proportional to X- ’ in which case X
U. VdX
u=
(11.13)
s contains a single logarithmic term, K In (-X), in its expansion about X = 0. Thus, if the general solution for F contains a logarithmic term, it must be of the form KF,(X) In (-X). Confirmation of this follows from the observation that in (II. 10) and (11.8), s, = (aoro/2bo)(pn/p,), n = 1, 2. . Expansions for the base functions d@/dX and dY/dX about X = - fv are d@ -= dX
W c
(II. 14)
&nXn,
n-0
where 8 = X+ $y and Z. = -4(D+coshfycos$-D_ 6, = 2(D_coshiycos’,y-D, with the remaining 80
coefficients
sinh$ysin’,y), sinhjysinfy),
then determined
G, = 0, h; = -a,
(II. 15)
from (11.3) by the relations Journal
of the Franklin Pergamon
lnst~tute Press plc
Third-order Eigenvalue
Problem in Convective
ci,= 46_ ,/n, = Zn_J(n+
L, The series expansion where
for Fabout
l)n(n-
(n = 2,3. . .).
1) 1
J? = 0 satisfying
F = dF/dX
Instability
(II. 16)
= 0 there is F = F,(X)
F3 = f j@, n=O
(II. 17)
p”,, = p”, = 0 and fi2 is arbitrary. Substitution of (II. 17) and (11.14) into (II. 1) leads to the determination of the remaining coefficients j?n from the formulae 40 = 0,
4”, = p”&ib”,,
B+ I = kAn+
lM-
11,
I (11.18)
ZZZ.Numerical Method The solution of the singular eigenvalue problem (II. l)-(11.2) is now constructed in three regions, one where the series expansion (II. 17) is valid, giving F= AF3(X) another
where the expansions
(-4~
d Xb
X0),
(III. 1)
(11.7) and (11.9) are valid, giving
F= BF,(X)+CF,(X)
(X, < X< O),
(111.2)
and a third (X0 d X < X,) where the solution F3 is continued from X0 to X, using a fourth-order Runge-Kutta integration of the original differential equation. It is assumed that F,, F2 and F3 are uniquely fixed by the specification p , = r. = jj2 = 1. At X = X0 the values of F3, dFJdX and d*FJdX* determined by (11.17) provide the initial conditions for the Runge-Kutta computation and at X, a consistent overall solution is obtained provided that F, dFjdX and d2F/dX2 are continuous, requiring
(111.3) to vanish at X = X,. This must be achieved by suitable choice of the parameter y. It is of interest to note that the determinant D is a function of y but not of X,, the location at which it is evaluated. This is easily established using the fact that each F, (i = 1, 2, 3) satisfies Vol. 329, No. I, pp. 77-83, 1992 Printed in Great Britain
81
3E gg
$71 $2
2 : E ,g 4
40 40 40
400 400 400
M
N
0.49 0.375 0.3125
XOIY
0.49 0.49 0.375 0.375 0.375 0.375 0.3125 0.3125
200 200 200 200 400 400 400 400
5 20 20 40 20 40 20 40
-
- XII/Y
M
N
0.01 0.125 0.1875
-X,/Y
0.01 0.01 0.125 0.125 0.125 0.125 0.1875 0.1875
--X,/Y
TABLE I
-0.000105 -0.000126 -0.000125
D(6.298284)
Computations
- 0.932430 - 0.932430 -0.932425 -0.932425 -0.932425 -0.932425 - 0.932425 -0.932426
D(6.28) 0.088974 0.088974 0.088982 0.088982 0.088979 0.088979 0.088983 0.088982
D(6.30)
- 0.000094 - 0.000069 -0.000072
D(6.298285)
- 0.000008 - 0.000022 -0.000021
D(6.298286)
ofD near y = 6.298286
TABLE II
- 0.426388 -0.426388 -0.426385 -0.426385 - 0.426382 -0.426383 - 0.426384 -0.426385
D(6.29)
Computations of D war y = 6.30
0.000003 0.000030 0.000030
D(6.298287)
1.148328
0.613829 0.613829 0.613834 0.613834 0.613834 0.613834 0.613836 0.613835
0.000070 0.000086 0.000084
D(6.298288)
1.148325 1.148329 1.148327
1.148325
1.148326 1.148326
1.148328
D(6.32)
D(6.31)
P
!?
;f. G-
b g
a
Third-order Eigenvalue d’F, ~ -QFi dX3
Problem in Convective
Instability
= 0.
Multiplication of the equation for F, by {F3 dF2/dX- F2 dF,/dX} and the equations for F2 and F3 by similar expressions in cyclic rotation, integration of the resulting equations by parts and addition shows that D as defined by (111.3) is independent of x.
IV. Results
and Conclusions
The recurrence formulae defined in Section II were used to generate the functions F,, F2 and F3 on a computer up to terms in X2N and 2’” at the locations X = X0 and X = Xi. The Runge-Kutta integration was performed using A4 equal intervals of X in the region X0 d X < Xi. Tests were carried out to investigate the effects of each of the parameters N, X,,, Xi and M and some of the results are summarized in Table I. The constancy of the determinant D for different values of X, also provided a useful check. At low values of y the two series expansions provide excellent approximations to Facross the entire domain but as y increases above 5 the two radii of convergence in x decrease until when y = yC z 6.3 they appear to not quite overlap near X = -:y. Some computations were performed with X0 = X, = -:y and produced quite accurate results for yC at very high levels of truncation (N up to 90) but it became evident that the incorporation of the RungeeKutta region was necessary to produce definitive results. From these (Table II) it was concluded that the first zero of D occurs at y = yC = 6.29829 and no further zeros of D were identified in the range y d 7.7. Acknowledgement This work was partially supported by a Visiting Science and Engineering Research Council.
Fellowship
Research
Grant
from the
References (1) P. G. Daniels, “Convection in a vertical slot”, J. Fluid Mech., Vol. 176, pp. 419441, 1987. (2) P. G. Daniels, “Stationary instability of the convective flow between differentially heated vertical planes”, J. Fluid Mech., Vol. 203, pp. 525-540, 1989. (3) J. W. Elder, “Laminar free convection in a vertical slot”, J. Fluid Mech., Vol. 23, pp. 77-98, 1965.
Vol. 329, No. 1. pp. 77-83, Pnnted in Great Britain
1992
83