Chaos, Solitons and Fractals 28 (2006) 1046–1054 www.elsevier.com/locate/chaos
Fractal structure of iterative time profiles of a tubular chemical reactor with recycle M. Berezowski *, K. Bizon Silesian University of Technology, Faculty of Mathematics and Physics, Institute of Mathematics, 44–100 Gliwice, ul. Kaszubska 23, Poland Accepted 24 August 2005
Abstract The scope of the paper is the theoretical analysis of the time rate in which a chemical reactor reaches a stable stationary state or stable temperature and concentration oscillations of the fluid flux. The method used for the analysis is based on the so-called iterative time profiles, demonstrating a chaotic and fractal nature of some of the profiles. The results were presented in the form of two and three-dimensional graphs. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction The dynamics of a homogeneous tubular chemical reactor without dispersion and with external feedback was discussed, for example, in [1–5], where it was demonstrated that the concentration and temperature of the fluid flux may be constant in the steady state or may oscillate in the form of rectangular time wave. The oscillations may be more or less complex, depending on the periodicity of the time series of the above-mentioned variables of a given state. Due to an asymptotic nature of the changes occurring in the reactor, the time rate required for the system to reach a steady state is, from the theoretical point of view, infinitely long. In practice, however, it may be assumed that the reactor is in the steady state when the values of its variables are relatively close to the values of the theoretical state. Thus, the time rate required for the system to reach the steady state is limited and depends on the values of the reactor parameters as well as on the values of temperature and concentration of the fluid flux at the initial moment. The reactor subjected to the theoretical analysis is a homogeneous tubular chemical reactor without dispersion and with external mass recycle [1,3, 5]. The time required for the reactor to reach the steady state was illustrated by means of various types of two and threedimensional iterative time profiles [6]. The profiles were demonstrated to have a chaotic nature and fractal structure.
2. The model of the reactor and its iterative time profiles The differential mass balance model of a chemical reactor with recycle [1,3] was analyzed in the following form: *
Corresponding author. Fax: +48 32 237 28 64. E-mail address:
[email protected] (M. Berezowski).
0960-0779/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2005.08.117
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
1047
Notations Da f n a b H d c k n
Damko¨hler number recycle coefficient order of reaction conversion degree coefficient related to enthalpy of reaction dimensionless temperature dimensionless heat transfer coefficient dimensionless activation energy eigenvalue dimensionless position coordinate along reactor
Subscripts 0 initial condition H heat exchanger s fixed point
dakþ1 ðnÞ ¼ ð1 f Þ/½akþ1 ðnÞ; Hkþ1 ðnÞ dn dHkþ1 ðnÞ ¼ ð1 f Þf/½akþ1 ðnÞ; Hkþ1 ðnÞ þ d½HH Hkþ1 ðnÞg dn akþ1 ð0Þ ¼ f ak ð1Þ; Hkþ1 ð0Þ ¼ f Hk ð1Þ bH . /ða; HÞ ¼ Dað1 aÞn exp c 1 þ bH
ð1Þ ð2Þ ð3Þ ð4Þ
The purpose of the analysis was to determine the impact of the initial values of conversion degree a0(1) and dimensionless temperature H0(1) on the number of discrete time steps N required for the reactor to reach a stable steady state, that is stable stationary point or stable periodic orbit [6]. It was assumed that the steady state for the reactor is the state when the distance between the trajectory generated by the model (1)–(4) and the fixed point as(1), Hs(1) is not longer than the one determined by the following condition: að1Þ as ð1Þ Hð1Þ Hs ð1Þ þ 100% < e. as ð1Þ Hs ð1Þ
ð5Þ
The analysis was based on the catastrophic set shown in Fig. 1 [1,3], generated for the following values of the parameters: Da = 0.15, n = 1.5, c = 15, b = 2, d = 3. Crossing of the lines SN, FB or HB of this set leads to the change in the multiplication factor (saddle-node bifurcation), generation of jump oscillations (flip bifurcation) or generation of quasiperiodic oscillations (Hopf bifurcation), respectively. Thus, right to the zone marked by the FB line, the temperature and concentration of the fluid flux are constant in the steady state, practically for any value of the recycle coefficient within the range of 0 < f < 1. Accordingly, what we are dealing with are fixed stationary points. But, within the boundaries of the zone marked by the FB line, the temperature and concentration oscillate in the steady state in a jump way [3]. Accordingly, we are dealing with fixed points of oscillation solutions. By analyzing the eigenvalues a of the model (1)–(5), it is possible to determine the ‘‘power’’ by means of which a given fixed point attracts trajectories. The most powerful attractor is the point for which all the eigenvalue modules are equal to zero, whereas the least powerful point is the one for which the module of at least one eigenvalue is equal to 1. The eigenvalues of k-periodic solution may be derived from the following relation: y k ¼ f k
R1 k k Y Y J dn e 0 j y 0 ¼ f k M j y 0 ¼ f k Ak y 0 j¼1
ð6Þ
j¼1
where J j is Jacobi matrix of the right sides of Eqs. (1) and (2), whereas y is the state vector of the linear approximation of the equations. The monodromy matrix M j may be derived from the following equation:
1048
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
1
f
0.9
SN
0.8
SN HB
0.7 0.6 0.5 0.4
FB
0.3
FB 0.2 0.1 0
ΘH – 0.09
–0.08
–0.07
–0.06
–0.05
–0.04
–0.03
–0.02
–0.01
0
Fig. 1. Catastrophic set of the reactor model.
dM j ¼ J jM j dn
ð7Þ
assuming, however, that M j ðn ¼ 0Þ is an identity matrix. The searched eigenvalues of the model (1)–(5) are the eigenvalues of matrix f k Ak . Furthermore, assuming that HH = 0.001 (right to the zone marked by FB), the bifurcation graph shown in Fig. 2 was designated, from which it may be indicated that the model of the reactor has a stationary solution within the full variation range of 0 < f < 1. Assuming that k = 1, the dependence between the modules of eigenvalues ki (i = 1, 2) and the recycle coefficient was shown in Fig. 3. As seen in the figure, the maximal eigenvalue occurs for f = 0.427. Accordingly, it is for this value of the recycle coefficient that the stationary state of the reactor is a weak attractor of the trajectories. This means a long way for the trajectory to lead the reactor to the state of dynamical equilibrium. In such case the number of recurrence steps N is big. To substantiate the above deliberation, the tree of the recurrence solutions is shown in Fig. 4, where, for a given value of f, different values of N were marked, corresponding to different initial values of the conversion degree 0 < a0(1) < 1. The maximal value corresponds to the maximal eigenvalue from Fig. 3. The calculations were based on the assumption that H0(1) = 0.2.
0.24
Θ s(1)
0.23 0.22 0.21 0.2 0.19 0.18 0.17 0.16 0.15
f
0.14 0
0.2
0.4
0.6
Fig. 2. Bifurcation diagram, HH = 0.001.
0.8
1
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
1
1049
| λi| | λ1 |
0.9 0.8
| λ1|=|λ2|
0.7 0.6 0.5 0.4 0.3 0.2
| λ2 |
0.1 0
f 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 3. The eigenvalues of the model, HH = 0.001.
1200
N
1000
800
600
400
200
f 0 0.25
0.3
0.35
0.4
0.45
0.5
0.55
0.6
Fig. 4. Tree of iterative time profiles, HH = 0.001.
Assuming that HH = 0.002 (the internal zone marked by the line FB), the bifurcation diagram shown in Fig. 5 was designated, where, apart from stationary solutions, twoperiodic oscillation solutions in the range of 0.3795 < f < 0.4695 may be observed. Assuming that k = 2, the dependence between the modules of eigenvalues ki (i = 1, 2) and the recycle coefficient was shown in Fig. 6. As seen in the figure, the maximal value from Fig. 3 corresponds (with big accuracy) to the minimum from Fig. 6. Because the minimum lies in the above-mentioned variation range of parameter f, it concerns an oscillation solution. For the value of the recycle coefficient f = 0.427 the twoperiodic orbit is a strong attractor of the trajectories. This means a short way for the trajectory to lead the reactor to the stable temperature and concentration oscillations. In such case the number of recurrence steps N is small. At the ends of the variation range of parameter f,
1050
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
Θ s(1)
0.24 0.23 0.22 0.21 0.2 0.19 0.18 0.17 0.16 0.15 0.14
f 0
0.2
0.4
0.6
0.8
1
Fig. 5. Bifurcation diagram, HH = 0.002.
| λ1 |
1
| λ1 | |λ 1 |=|λ2|
0.8
0.6
0.4
0.2
|λ2 |
f
0 0
0.2
0.4
0.6
0.8
1
Fig. 6. The eigenvalues of the model, HH = 0.002.
that is, for f = 0.3795 and f = 0.4695, the eigenvalue is equal to 1. Accordingly, at this particular position the fixed point is the least powerful attractor of the trajectories. To substantiate the above deliberation, another tree of recurrent solutions is shown in Fig. 7, where, the minimum corresponds to the minimal eigenvalue from Fig. 6 in the range of 0.3795 < f < 0.4695. The two maximums that are clearly seen in the figure designate the range of oscillations FB, equal to the range presented in Figs. 5 and 6. By decreasing the value of the dimensionless temperature of the cooling medium the model generates another FB bifurcation, which leads to the doubling of the period of the present oscillation solution. In consequence, for HH = 0.012, fourperiodic oscillations in the range of 0.381 < f < 0.430 are generated. To substantiate the above deliberation, another tree of recurrent solutions is shown in Fig. 8, where, the three minimums correspond to strong attraction by a fourperiodic attractor, and the four maximums to FB bifurcation points—which are weak attractors for the trajectories. Decreasing further the value of HH other FB bifurcations are obtained, generating the successive doubling of the period [3]. Each of the bifurcations offers (k 1) minimums, i.e. (k 1) strong attractors of the trajectories and k maximums, i.e. k weak attractors.
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
7000
1051
N
6000
5000
4000
3000
2000
1000
f
0 0.25
0.3
0.35
0.4
0.45
0.5
0.55
0.6
Fig. 7. Tree of iterative time profiles, HH = 0.002.
7000
N
6000
5000
4000
3000
2000
1000
0 0.15
f 0.2
0.25
0.3
0.35
0.4
0.45
Fig. 8. Tree of iterative time profiles, HH = 0.012.
0.5
0.55
0.6
1052
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054 N
1200
N
900
1100
800
1000
700
900
600
800
500
700 400 600 300
500
200
400
100
300 200
a
α (1) 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
b
α0(1) 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 9. (a) Iterative time profile, HH = 0.001, f = 0.427, H0(1) = 0.2. (b) Iterative time profile, HH = 0.002, f = 0.427, H0(1) = 0.2.
N
900 800 700 600 500 400 300 200 100 0
α0(1) 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 10. Iterative time profile, HH = 0.012, f = 0.427, H0(1) = 0.2.
Nj+1
900 800 700 600 500 400 300 200 100 0
Nj 0
100
200
300
400
500
600
700
800
900
Fig. 11. Poincare´ cross-section, HH = 0.012, f = 0.427, 0.4 < a0(1) < 0.42.
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
1053
For example, assuming that HH = 0.001, f = 0.427 and H0(1) = 0.2 an iterative time profile demonstrating the impact of the initial value of the conversion degree on the number of recurrence steps N required for the reactor to reach a stable stationary state is shown in Fig. 9a with the accuracy of e = 0.001%. For a0(1) = 0.559 the visible peak directed downwards certifies that the initial condition of the reactor is located closest to the stationary attractor. Likewise, in Fig. 9b the iterative time profile for HH = 0.002, f = 0.427 and H0(1) = 0.2 demonstrates the impact of the initial value of the conversion degree a0(1) on the number of recurrence steps N required for the reactor to reach a stable limit cycle. For a0(1) = 0.546 the visible peak directed upwards certifies that the initial condition of the reactor is located farthest from the twoperiodic attractor. In view of an insignificant discrepancy in the values of HH in both of the above-mentioned cases, the peaks are also insignificantly displaced from each other, yet directed oppositely. The change
Fig. 12. Three-dimensional iterative time profile, HH = 0.012, f = 0.427.
Fig. 13. Fragment of the profiles from Fig. 12.
1054
M. Berezowski, K. Bizon / Chaos, Solitons and Fractals 28 (2006) 1046–1054
in the direction of the peak is a result of crossing the bifurcation line FB on the graph from Fig. 1. This regularity also concerns the next bifurcations, leading to the successive doubling of the period of the oscillation solutions. In Fig. 9b the visible peak directed downwards [a0(1) = 0.048] certifies that the initial condition of the reactor is located closest to the twoperiodic attractor. Further decrease in the value of parameter HH leads to the occurrence of successive FB bifurcations, which, in consequence, leads to the generation of more and more peaks on the iterative time graph. The number of the peaks may increase very rapidly. In Fig. 10 the iterative time profile for HH = 0.012, f = 0.427 and H0(1) = 0.2 has a chaotic nature in the range of 0.4 < a0(1) < 0.42. This chaos is revealed by a great sensitivity of this part of the graph to the infinitesimal change in the position of a fourperiodic attractor, as well as by a disordered cloud of the points on Poincare´ cross-section (Fig. 11). The number of iteration steps N corresponding to the jth and (j + 1)th peak from Fig. 10 is marked on the coordinates of the cross-section. Moreover, the graph in Fig. 10 was found to have a fractal nature, which was demonstrated on three-dimensional iterative graphs (Figs. 12 and 13). Other exemplary fractal diagrams of iterative time profiles of the discussed reactor are to be found in [7].
3. Concluding remarks The scope of the paper was a theoretical analysis of the dynamics of a homogeneous chemical reactor without dispersion and with mass recycle. Particular attention was given to the impact of the process parameters and the initial values of the concentration and temperature of the fluid flux on the time rate required for the reactor to reach steady states. In the course of the research the so-called iterative time profiles were designated and their chaotic and fractal nature was shown. The results are useful from the view point of process technology, and interesting from the theoretical point of view. Due to substantial and significant mass effects and thermal effects caused by the chemical reaction, the discussed method makes it possible to optimize the design of the reactor system, whereas the results give an insight into the dynamics of chemical reactors.
References [1] Berezowski M. Spatio-temporal chaos in tubular chemical reactors with the recycle of mass. Chaos, Solitons & Fractals 2000;11:1197–204. [2] Berezowski M, Grabski A. Chaotic and non-chaotic mixed oscillations in a logistic system with delay and heat-integrated tubular chemical reactor. Chaos, Solitons & Fractals 2002;14:97–103. [3] Jacobsen EW, Berezowski M. Chaotic dynamics in homogeneous tubular reactors with recycle. Chem Eng Sci 1998;53:4023–9. [4] Jacobsen EW, Berezowski M. Dynamics of heat-integrated homogeneous tubular reactors. In: 5th IFAC symposium on dynamics and control of process systems. Corfu, Greece, June 8–10, 1998. [5] Luss D, Amundson NR. Stability of loop reactors. AIChE J 1966;13:279–90. [6] Peitgen H-O, Ju¨rgens H, Saupe D. Fractals for the classroom. Part 2: Complex systems and Mandelbrot set. New York: Springer; 1992. [7] http://zeus.polsl.gliwice.pl/~mberez/galstud/galstud.html.