Chemical Engineering Science. Vol. 41, No. Printed in Great Britain.
11, pp. 2767-2777,
$3.00 + 0.00 Journals Ltd.
ON THE DETAILED DYNAMICS OF COUPLED CONTINUOUS STIRRED TANK REACTORS SANG H. KIM and VLADIMIR HLAVACEKt Departmentof Chemical Engineering,State Universityof New York at Buffalo, Buffalo,N Y 14260, U.S.A. (Received 30 July 1984) Abstract-The dynamicbehaviourofcoupled paralleland serial arrays of cells is analysed. Phenomena such as synchrnization of oscillations at a common frequency, rhythm splitting, amplitude amplification, jump between oscillatory regimes and quasi-periodic oscillations of a torus type have been observed for an autocatalytic reaction with product inhibition. The possibility of steady-state operation driven by the coupling of two cells, each of which is oscillatory by itself, has been detected. The configuration of the cell arrays does not change the bifurcation pattern drastically.
A set of coupled cells with mutual mass and/or heat exchange can often be used as a prototype of model systems to describe heterogeneous catalytic reactors, surface diffusion on a matrix of solid catalysts, a packed bed reactor, and pattern formation in living organisms. The idea of coupled cells originated from Turing’s (1952) pioneering work early in 1952. Studies on the dynamic behaviour of coupled cells were initiated by many biologists; Goodwin (1969) proposed a phase-shift model for spatio-temporal organization in developing systems. Pavlidis (197 1) discovered the interesting phenomena of synchronization (phase locking) and rhythm splitting in a model system for circadian clocks. Burton and Canham (1973) analysed contact inhibition in cell division. Murray (1982) also investigated space instability for several models proposed to elucidate Turing’s pattern formation. Bar-Eli and Geiseler (1981) observed experimentally interesting features of the relative stability of coupled cells. Marek and Stuchl(l975) reported several important phenomena such as synchronization of oscillations at either a single or multiples of a common frequency, rhythm splitting, and amplitude amplification in two coupled oscillators for the Belousov-Zhabotinskii reaction. Neu (1979 a, b), using the singular perturbation technique, analysed synchronization and the pertinent bifurcation phenomena from this state. Marek and Schreiber (1982) showed the existence of a strange attractor in two coupled cells for a Brusselator model. Kumar et al. (1983) investigated the dynamic pattern of coupled cells in parallel for a chain of autocatalytic reactions with an exponential acceleration step. Jensen and Ray (1980) observed a single peak and complex multipeak oscillations on catalyst patches described by a so-called “fuzzy wire” model. It is known that a system of two coupled cells in nonoscillatory states may be driven into oscillations (Smale, 1974). The possibility of oscillations by couptTo whom correspondence should be addressed. 2767
ling two non-oscillatory cells was discussed by Tyson (1979). He examined two coupled cells, each described by an “Oregonator” model. Coupling the cells through a membrane, permeable only to bromous acid, can produce oscillations. As far as we know, no systematic studies have been reported on the dynamic behaviour of two or more coupled cells in different configurations, namely in parallel and in series (Fig. 1). In the present paper we attempt to comprehend the dynamic pattern of cell arrays in two different configurations. An autocatalytic reaction with product inhibition is considered to occur in the system. To achieve the goals of a detailed investigation, bifurcation analysis is adopted. The possible appearance of chaotic motion and the effect of the system size (number of cells) on the dynamic behaviour of coupled cells will also be discussed.
Consider the following chain of reactions with an autocatalytic step. This reaction exhibits a product inhibition step as seen in eqs (1) and (2): k3 r-------7
1 c3 -
dC, = klC1 -k3C1C3
= klC, - k2C2
-dC3 = k,C, - k3C3. dt
For a parallel configuration the mass balances in dimensionless form represent the differentialdifference equations given below [see eq. (3)]. Isothermal operation, an identical volume of every cell, the same feed composition to all cells, equivalence of the diffusive coupling rates for every component, and
Fig. 1. Schematic sketch of coupled cells: (a) in parallel and (b) in series.
constant physical properties are assumed throughout this paper. The governing equations are dut ~~=u~-~u,+Da,(u,-~cr,u,w,)+~(uz--u,) d7
u” - 74;+ Da,(ui+P(ui-,
3, . .
Da = Kk,, 4 V’
C; vi=G, al =-
, 3. ANALYSIS AND NUMERICAL
dw. I = w”-w~+Daj(a2ui-a3wi) dr i = 2,
vi = UP), wi = wt”.
The mass balance for coupled cells in series can be easily obtained in a similar way. Here we have assigned u, u, w the dimensionless concentrations of species Cr , Cz and Cj, respectively; aI, a2, a3, the dimensionless kinetic constants; Da, the Damkiihler number, p, the coupling coefficient for mass diffusion between the neighbouring cells, and T, the dimensionless time:
ui = nl”,
dvi = u0-zJ+Da,(ui-a2ui) dr +P(ui-r
dw, ---=wO-ww,+Dal(cr,o,-a0(3w1)+~(w2-wW1) d7
subject to the initial conditions
dr, ~=u~-uv,+Da,(u,-~a,v,)+~(u,-~,) d7
Numerical integration of a set of ordinary differential equations, eq. (3), was perf-ormed using Gear’s method for the integration of stiff systems of ordinary differential equations. The error of integration was controlled to six significant digits. The bifurcation analysis for tracing out the periodic as well as the stationary solution branches was performed using the software package AUTO (Doedel, 1980).
Dynamics of coupled CSTRs 3.1 A single ceil For a single CSTR, i.e. p = 0, eq. (3) is represented by a three-variable system: du -= ds
dv _=_ dr
dw -= dr
The type of steady state (us, us, w,) is determined by the eigenvalues of the following Jacobian matrix: afi _~_ au
evaluated at u = us, u = us and w = ws. We consider in this text only one component, C1, in the feed (i.e. u” = 1.0, v” = Wo = 0.0) and the kinetic parameters, cr,, CQ and aj, are fixed at a1 = a2 = cc3 = 0.2. Detailed analysis of the multiplicity and stability of steady states in a CSTR has been described in the literature (Kim and Hlavacek, 1986). The bifurcation diagram in Fig. 2 shows that a stable node, a stable focus and an unstable saddle-focus may exist in a CSTR. For higher values of the Damktihler number a stable limit cycle may occur around an unstable focus. Therefore, low flow rates favour oscillations while high flow rates support the steady-state operation in a CSTR. The steady-state concentrations and the periods of limit cycles in a CSTR are summarized in Table 1. The period of oscillations tends to increase with decreasing value of the Damkohler number (see Fig. 3). 15
__ - --__ I IO
Do Fig. 2. Characteristics of a steady state in a CSTR. Q, = a2 = CL3= 0.2.----, Stable steady states; ---, unstable steady states; A-..,stable periodic solution; 0, Hopf bifurcation point.
In the bifurcation diagrams the ordinate shows the measure of the concentration u, (u). The measure represents the steady-state value for stationary solution or a maximum value of a concentration u for periodic solution. This enables the stationary and periodic solutions to be plotted in the same graph. 3.2 Coupled cells The transient behaviour in coupled cells can be obtained by integration of eq. (3) for certain values of the initial conditions, @, t$) and wp)_ The initial conditions for the coupled ceil model were always chosen as an asymptotic solution to the uncoupled array of cells. Thus, to determine the initial conditions we have integrated eq. (6), representing a single cell in uncoupled state with initial conditions u = 3.0, v = w = 4.0. When the asymptotic regime was stabilized (oscillations or steady states), the algorithm switched to integration of the coupled cells. Based on the previous dynamic analysis for a single cell (cf. Fig. 2), the dynamic behaviour for coupled cells will be examined for different types of combinations, viz. “node-focus”, “oscillation-focus” and “oscillation-oscillation”. (a) Node-focus. The 3.2.1. Parallel conjiguration. bifurcation diagram for node-focus combination (i.e. Da, = 1.0, Daz = 0.1) reveals that the steady-state concentrations in each cell are different. They approach the same steady-state value with increasing value of interaction (see Fig. 4). In this text the subscript i in Da, and the number (1 or 2) in the bifurcation diagrams represent the cell number. (b) Oscillation focus. The typical bifurcation pattern for oscillation-focus combination, as shown in Figs 5 and 6, exhibits both steady and oscillatory states. The Hopf bifurcation point is detected for a very weak interaction in case 1 (Da, = 8.0, Do, = 1.0) while the steady state is reached for the whole range of coupling coefficient for case 2 (Da, = 8.0, Da, = 6.0) (see Figs 5 and 6). This fact reveals that the oscillation process in one cell can excite a stable steady state for a very weak interaction if the difference in the steady-state values of uncoupled cells is large. For both cases, the steadystate concentrations reached in both cells are different not only from each other but also from those of uncoupled cells. Amplitude amplification for very low values of the coupling coefficient for case 1 was detected. For case 1 we also observed a kind of jump phenomenon between oscillatory states during the transient operation. The trajectory suddenly jumped from an outer to an inner cycle cf. Fig. 7(a). The effect of the values of the coupling coefficient on the dynamic behaviour of coupled cells for case 1 reveals the existence of oscillations and stationary behaviour (see Fig. 8). (c) Oscillation-oscillation. Two different types of coupling are considered in order to investigate the effect of period difference on the dynamic behaviour of coupled cells (i.e. case 1: Da, = 8.0, Da2 = 16.0, case 2:
Table 1. Steady-state concentrations and periods of oscillation in a CSTR Period Da
Mode Stable node
0.89 1.00 2.00 6.00
4.8665 6.cKnIO 7.oOOO 2.9884
3.6939 5.00 10.0000 8.1501
0.5608 0.8333 2.857 1 4.4455
4.5597 4.6306 4.7347 4.8736
7.08t 8.00 10.00 16.00 tHopf
2.6552 2.4455 2.1306 1.6791
7.7803 7.5247 7.1020 6.3966
2.2709 2.0902 1.7843 1.2491
1.09 1.00 0.85 0.60
2.4552 4.6195 5.5234 6.0754
bifurcation point (HB).
Fig. 3. Period vs. Damkdhler number in a CSTR. = CL3= 0.2.
cxl = cyz
I 5 P
Fig. 5. Bifurcation diagram for two coupled cells in parallel. Da, = 8.0, Da, = 1.0. ~ > Stable steady states; ---, unstable steady states;. - - ., stable periodic solution; 0, Hopf bifurcation point.
Fig. 4. Bifurcation diagram for two coupled cells in parallel. Da, = 1.00, Da2 = 0.10. -, Stable steady states.
Da, = 8.0, Dal = 10.0). The bifurcation diagrams shown in Figs 9 and 10 reveal that the coupled cells may operate in a regime of quasi-periodic oscillations (torus), single peak oscillations of a limit cycle type or at a steady state. The type of global behaviour depends on the degree of interaction. The number of Hopf bifurcation points tends to increase with decreasing difference of periods between
Fig. 6. Bifurcation diagram for two coupled cells in parallel. Da, = 8.0, Da, = 6.0. -, Stable steady states.
Dynamics of coupled CSTRs
Fig. 7. Jump phenomena for two coupled cells in parallel. Da, = 8.0, Daz = 1.0, p = 0.02. Cell No. 1.
a; two oscillators. For case 1, the very weak interaction, p < 0.038, results in quasi-periodic oscillations while for the intermediate value of p, 0.308 -=I p -c 5.3 16, steadystate operation is established. For p > 5.316 and p
< 0.308, oscillations occur in both coupled cells. A large difference of amplitude of oscillations in both cells is observed in the low range of coupling coefficients. The amplitude of an oscillator with a long period is reduced to a low level. This small amplitude of oscillations may result from the interaction of the neighbouring oscillator with a short period. The period of oscillations increases insignificantly with increasing “distance” of coupling coefficients from the Hopf bifurcation point. Ruelle (1973) has shown that due to the interaction of two oscillating systems (periods of oscillations T, and T,), oscillation with the periods mT, and nT, can arise, where m and n are small numbers and n/m - Tl /Tz. Some periods of coupled oscillators for case 1 are presented in Table 2. If the difference between the frequency of oscillations at a given intensity of interaction (relatively high) is small, synchronization at the same frequency of oscillations in both cells can occur. This means that oscillations in the system with lower frequency are
s ;::. O-1 .oo
Periods u 0.10
0.20 0.30 5.50 6.00 10.00
1.2166 1.1878 1.1653 1.4602 1.4699 1.5076
1.1236 0.9500 1.1750 1.4350 1.4500 1.5076
t’denotes the condition after coupling.
$. 2 >
Table 2. Periods of oscillations vs. coupling coefficientsfor case 1 (Dal = 8.0, Dal = 16.0)
oi 8 ; 8 O-i
Fig. 8. Transientbehaviour of two coupled cells in parallel. Da, = 8.0, Daz = 1.0. (a) Cell No. 1. ~1= 0.20, (b)cell No. 2, p = 0.20; (c) cell No. 1. p = 10.0; (d) cell No. 2, p = 10.0.
and VLADIMIR HLAVACEK
Fig. 9. Bifurcation diagram for two coupled cells in parallel. Da, = 8.0, Da, = 10.0. -, Stable steady states; ---, unstable steady states; . . . ., stable periodic solution; 0 0 0, unstable periodic solution; + + +, quasi-periodic oscillations; 0, Hopf bifurcation point; I, bifurcation point to a torus.
Fig. 10. Bifurcation diagram for two coupled celis in parallel. Da, = 8.0,Da, = 16.0. ----, Stable steady states; ---. unstable steady states;. . ., stable periodic solution; 0 0 0, unstable periodic solution; + + + , quasi-periodic oscillations; 0, Hopf bifurcation point; m , bifurcation point to a torus.
driven by oscillations initially with higher frequency. This tendency is illustrated in Fig. 11. The original
period of oscillation in the first and second cells before coupling was T, = 2.09 and T, = 1.25, respectively. The periods of coupled cells for p = 10.0 were the same and equal to T; = FL = 1.5 1. Thus both oscillators are synchronized at the highest oscillation frequency. Several entrainment phenomena (or synchronization) such as strict entrainment, loose entrainment, pure free-running and jittery free-running entrainment were studied in detail by Nicolis et al. (1973). The possibilities of synchronization increase with increasing degree of interaction and decreasing difference between the periods of interacting oscillators. The synchronization regions in the coupling coefficient
domain are p 2 2.30 and ,U 2 10.0 for cases 2 and 1, respectively. It is an interesting feature of oscillatory systems that a coupled system of two cells operating at oscillatory states before coupling may be driven into steady-state operation. The region of coupling coefficient p to give steady-state operation is enlarged with increasing difference of the periods between two oscillators. The increasing number of cells in the system favours oscillatory regimes and eliminates steady-state operation. This tendency is evident after comparing Fig. 12 (four coupled cells in parallel) with Fig. 10. The stability of a periodic solution is determined by the Floquet multipliers, i.e. by eigenvalues of the monodromy matrix. In order to calculate the Floquet
Dynamics of coupled CSTRs
Fig. 11. Transient behaviour for two coupled cells in parallel. Da, = 8.0, DaZ = 16.0. p is: (a) 0.00, (b) 0.02, (c) Cell No. 1; ---, cell No. 2. 0.20, (d) 2.00, (e) 4.00, (f) 6.00, (g) 10.0. -,
multipliers, eq. (3) can be written in the following standard vector form. Here the vector X(t) has three components, u(t), u(t) and w(t): X(L) = f(X).
The monodromy matrix M (t) is given by the following equation:
Here X* denotes values for a limit cycle and I designates an identity matrix. The eigenvalues of the monodromy matrix M (T), calculated at a period T, are called the Floquet multipliers. One of these is always + 1 and the periodic solution is stable if the others lie inside the unit circle. When one of the Floquet multipliers leaves the unit circle, a limit cycle loses its stability and bifurcates to an another state. The nature of this bifurcation can often be analysed based upon the position on the unit circle where the multipliers
SANG H. KM
and VLADIMIR HLAVACEK
’ . __---.----&=_______ v 2 ,-
Fig. 12. Bifurcation diagram for four coupled cells in parallel. Da, = 8.0, Da, = 16.0, Da, = 16.0, Da, Stable steady states; ---, unstable steady states; . . . ., stable periodic solution; 0 0 0, = 16.0. p, unstable periodic solution; + + + , quasi-periodic oscillations; 0, Hopf bifurcation point; I ,bifurcation point to a torus.
cross. The case of crossing at - 1 implies perioddoubling bifurcation while crossing the unit circle at a pair ofcomplex points, i.e. a f /?i (a2 + f12 = l), leads to a quasi-periodic oscillation (torus). Leaving the unit circle at + 1 displays a bifurcation of periodic solution to another periodic branch. The Floquet multipliers at a point where a limit cycle becomes unstable and is replaced by a torus are illustrated in Table 3. Evidently, a* +/I2 = (0.8134)2 + (0.5817)’ x 1.0000 for case 1 verifies the occurrence of quasi-periodic motion. Figure 13 reveals the quasi-periodic oscillations of a torus type for case 2. When the conditions for synchronization are not fulfilled, the rhythm splitting phenomenon may occur. Rhythm (or period) for the slower oscillator is frequently split into two parts by the influence of the faster oscillator through the coupling. Such an effect is shown in Fig. 14. The original period of oscillations with the lower frequency, T, = 2.10 [cf. Fig. 14(a)], is split into two parts of oscillation, T; = 1.10 and TZ = T, - T; = 1.00.
of the interacting reactors in parallel. Complex oscillations, like a torus for the coupled parallel cells, have not been detected in the series configuration. For the sake of simplicity, we consider the situation where a reactor with Da = 8.0 is placed in front of the other cells. For an “oscillation-focus” combination, the region of oscillations is enlarged compared with that for a
3.2.2. Series conjiguration. The bifurcation pattern for coupled cells in series exhibits similar behaviour to that
multipliers in two coupled parallel
Case 1 (Dal
= 8.0. Da,
1.0000 -0.8134 +O.S817i 0.8611 x 10-s -0.8134AJ.5817i 0.6881 x 10-s 0.3977
Case 2 (Dn,
= 8.0, Da,
1.0000 0.6191 + 0.7853i 0.4218 x 1O-4 0.6191AI.7853i 0.769 1 x 10 5 0.7612
Fig. 13. (a)-(b)
Dynamics of coupled CSTRs
> =: d
parallel configuration (see Figs 15 and 16). This behaviour may result from the fact that the output oscillation concentration from the first cell becomes the forcing function fed to the second cell. For the interaction, similar be“oscillation-oscillation” haviour, except for the existence of quasi-periodic oscillations, to that in parallel is detected (see Figs 16 and 17). CONCLUSIONS
Fig. 13. Quasi-periodic oscillation (torus) for two coupled cells in parallel. Da* = 8.0, Da, = 10.0, p = 0.02. (a,~) Cell No. 1; (b,d) cell No. 2.
The dynamic behaviour of coupled parallel and serial cells in which an autocatalytic reaction with product inhibition occurs has been studied in the present paper. For a “node-focus” combination, both cells approach a steady state, which is different for each cell. For an “oscillation-focus” interaction, steady-state operation is obtained for the entire domain of values of coupling coefficients. For very weakly interacting cells, with a different steady state for a disconnected configuration, oscillations of a limit cycle type can be established. Amplitude amplification and jump
Fig. 14. Rhythm splitting for two coupled cells in parallel. Da, = 8.0, Da, = 16.0. Before coupling, (a) cell No. 1, (b) cell No. 2. After coupling o( = 0.29, (c) cell No. 1, (d) cell No. 2.
phenomena between two limit cycles has also been detected. For an “oscillation-oscillation” combination, single peak oscillations, quasi-periodic oscillations (torus) and steady-state operation result according to the values of the interaction coefficient cc. For very weak couplings (p -C 0.038 for Da, = 8.0, Da2 = 16.0), quasi-periodic oscillations have been observed while for intermediate values of interaction, i.e. 0.308 < p -z 5.316, steady-state operations have been established. Oscillations of a limit cycle type can be obtained for p < 0.308 and p > 5.316. Synchronization phenomena at the common frequency have been detected. Oscillations in the system with lower
SANG H. KIM
and VLADIMIR HLAVACEK
Fig. 15. Bifurcation diagram for two coupled cells in series. Da, = 8.0, Da, = 1.0. -, StabIe steady states; ---, unstable steady states; . . . ., stable periodic solution; 0, Hopf bifurcation point.
frequency are driven by the oscillations initially with a higher frequency. The increasing number of cells
I ::* l -
Fig. 16. Bifurcation diagram for two coupled cells in series. Da, = 8.0, Da, = 6.0. -, Stable steady states; ---, unstable steady states;. ., stable periodic solution; 0, Hopf bifurcation point.
favours oscillatory phenomena at the expense of steady-state operation. The occurrence of synchronization increases with increasing degree of interaction and decreasing difference of periods between the interacting oscillators. The bifurcation pattern calculated for coupled cells in series reveals behaviour similar to that of parallel arrays. This study has been motivated by our attempt to understand the dynamic behaviour of a distributed flow system, e.g. tubular flow reactors, interaction of isolated patches on the catalyst surface, and cooperation between separated biological entities. Acknowledgemenl-The bifurcation diagrams reproduced in this paper were calculated by the bifurcation package AUTO, which was provided by Dr. E. Doedel, Computer Science
Fig. 17. Bifurcation diagram for two coupled cells in series. Dal = 8.0, Da2 = 16.0. ~ 3Stable steady states; ---, unstable steady states; ., stable periodic solution; Q, Hopf bifurcation point.
Dynamics of coupled CSTRs Department,Concordia University,Montreal. His assistance is greatly appreciated. NOTATION A c1,
&I 4 4c t i-1,
V x Greek c!
Jacobian matrix defined by eq. (7) concentrations (or chemical species) involved in a reaction concentration of species C1 in a feed Damkijhler number in the i-th cell right-hand side defined by eq. (6) identity matrix reaction rate constants an integer satisfying n/m - T,/T2 monodromy matrix flow rate in the axial direction back-flow rate time periods of oscillation before coupling periods of oscillation after coupling measure of concentraaion u dimensionless concentrations of species C1, C2 and C3, respectively cell volume concentration vector
real part of a complex multiplier dimensionless kinetic constants defined
by eq. (5) imaginary part of a complex multiplier
REFERENCES Bar-Eli, K. and Geisler, W., 1981, Mixing and relative stabilities of pumped stationary states. J. phys. Chem. 85, 3461-3468. Burton, A. C. and Canham, P. B., 1973, The behavior of couuled biochemical oscillators as a model of contact inhibition of cellular division. .7. theor. Biol. 39, 555580. Doedel, E., 1980, AUTO: a program for the automatic bifurcation analysis of autonomous systems. Proc. 10th Manitoba Conf. Numerical Math. Computing, Winnipeg,
Canada. Goodwin, B. C., 1969, A phase-shift model for the spatial and
temporal organization of developing systems. J. theor. Biol. 2549-107. Jensen, K. F. and Ray, W. H., 1980, A new view of ignition, extinction, and oscillations on supported catalyst surfaces. Chem. Enana Sci. 35, 241-248. Kim, S. H. and Hlavacek, V., 1986, Qualitative properties of autocatalytic reactions occurring in a flow system. Z. fiir Nnturforsch. in press. B. D. and Kumar, V. R., Jayaraman, V. K., Kulkami, Doraiswamy, L. K., 1983, Dynamic behavior of coupled CSTRs operating under different conditions. Chem. Engng Sci. 38, 673686. Marek, M. and Schreiber, I., 1982, Strange attractors in 5D, 258-272. coupled reactionAiffusion cells. Physica Marek, M. and Stuchl, I., 1975, Synchronization in two interacting oscillatory systems. Biophys. Chem. 3, 241-248.
Murray, J. D., 1982, Parameterspace for Turing instabilityin reaction+tiffusionmechanisms:a comparison of models. .?. theor. Biol. 98. 143-163.
Ndi, J. C., 1979a, dhemical waves and the diffusive coupling of limit cycle oscillators. SIAM J. appl. Marh. 36, 509-515. Neu, J. C., 1979b, Coupled chemical oscillators. SIAM J. appi. Math. 37, 307-315. Nicolis, J. S., Galanes, G. and Protonotarios, E. N., 1973, A frequency entrainment model with relevance to systems displaying adaptive behavior. Int. J. Control l&1009-1027. Pavli&s; T.; 197 I, Populations of biochemical oscillators as cercadian clocks. J. theor. Biol. 33, 318-338. Ruelle, D., 1973, Some comments on chemical oscillations. Trans. N. Y. Acad. Sci. 35. 6&71 Smale, S., 1974, Lect. Math.‘Life Sci. 6, 17-26.
Turing. A. M., 1952, The chemical basis of morphoeenesis.
initial conditions inlet condition condition on a limit cvcle
Phil’ Trans. R. Sot. B237, 37-72. Tyson, J. J., 1979, Oscillations, bistability, and echo waves in models of the Belousov-Zhabotinskii reaction. Ann. N. Y. Acad. Sci. 316, 279-295.