- Email: [email protected]
Travelling waves and chaos in the Kolmogorov-Spiegel-Sivashinsky model
0020-7225192 $5.00 + 0.00 Copyright @ 1992 Pergamon Press plc
Int. J. Engng Sci. Vol. 30, No. 5, pp. 593-610, 1992 Printed in Great Britain. All rights reserved
TRAVELLING WAVES AND CHAOS IN THE ISOLMOGOROV-SPIEGEL-SIVASHINSKY MODEL G. ONAL
and E. S. SUHUBl
Istanbul Technical University, Faculty of Sciences, Maslak 80626, Istanbul, Turkey Abstract-Kolmogorov-Spiegel-Sivashinsky (KSS) equation has been derived to model large-scale structures arising in two-dimensional turbulence [Null. Phys. (Proc. Suppi.) 2, 453 (19t37)j. In the present study, we have obtained its periodic, quasi-periodic and solitary wave solutions analytically to a certain degree of approximation by discarding the linear damping term. Normal form analysis succeeds to capture Hopf, Neimark and infinite-period bifurcations before the onset of chaos. Melnikov analysis has been carried out to identify the homoclinic bifurcation. Transient spatiotemporal chaos has been observed, by constructing the Poincart-section and it has been identified by computing Lyapunov Characteristic Exponents (LCE).
1.
INTRODUCTION
The spontaneous formation of large-scale structures in a viscous fluid driven at small scales is one of the challenging problems of hydrodynamic turbulence. The latter ensues from the corresponding small-scale flow as a manifestation of the long-wave instability [2]. Kolmogorov flow and large-scale turbulent solar convection possess these features [l, 21. The former is a two-dimensional viscous fluid flow subjected to a unidirectional external force field periodic in one of the space coordinates [2]. KSS equation is capable of modeling large-scale structures in flows mentioned above. A specific form of the KSS equation can be written as $* + BQ, + Y& + (a - 3W,2MJXX+ #XXXX =0
(1.1)
where (Y, fi, y and 6 are physical positive constants. Sivashinsky [2] derived the equation (1.1) (with LY= 2) for a unidirectional large-scale flow parallel to a direction where the effective viscosity of the secondary flow is negative in the context of Kolmogorov flow. In this case, +(x, t) is the resealed large-scale stream function obtained after subtracting the mean periodic field component, x is the resealed preferential direction of negative viscosity and E is the resealed time variable. Spiegel et al. [3] derived an equation similar to (1.1) in multiple-scales analysis of a compressible solar convection layer zone. Their equation also includes a hydrodynamic blow-up term and a term carrying non-Boussinesqian compressibility effects. Circumventing the latter terms leads to equation (1.1) with #(x, t) representing the large-scale horizontal fluctuations of the temperature field after subtraction of the mean field and x being a large-scale horizontal variable [l]. The role played by each term in equation (1.1) has been expounded in references [ 1,4] (see also [5]). In the present study, linear damping term will be annulled by setting the coefficient /3 = 0. This is merely due to the fact that travelling wave solution found in our previous work [5] does not remain bounded for & = 0. In the same work, we were able to obtain the following dynamical system by resorting to group theoretical means: zi1= u2 l&=l+ iiig= l&q . a4
=
-j?u~
+
clip
-
cm3
+
3&&
- yu$
(1.2)
and z=xct, and c is the celerity of waves. Solutions to this dynamical where u 1=#(z) system are known to be travelling waves which they retain their shape at constant celerity. 593
G. iiNAL
594
and E. S. .$UHUBi
When c = 0 and /3 = 0, eigenvalues of the Jacobian derivative gradient matrix of the right hand side of (1.2) become a double zero and a pair of complex conjugate pure imaginary numbers. This situation introduces an extra degeneracy to the problem. The latter can be circumvented by the following change of variables u2 = ul, u3 = u2 and uq = us. Thus, the dynamical system which we are concerned becomes,
ti,
=
CV1 -
&V2 + %J:V2
-
J’V:.
(1.3)
Our main objective here is to study the local properties of the dynamical system (1.3) in the neighborhood of an equilibrium point. In Section 2, we identify the stability of the equilibrium points by employing the linear stability analysis. Normal form analysis has been performed in Section 3. In Section 4, periodic, quasi-periodic and solitary wave solutions to (1.3) are constructed. Behavior of the solutions obtained in Section 3 under the perturbation of higher order terms are investigated by resorting to Melnikov analysis in Section 5. In Section 6, transition to chaos has been discussed and numerically demonstrated by computing Poincare sections and LCEs.
2. LINEAR
STABILITY
Dynamical system (1.3) has two equilibrium
ANALYSIS
points located at
v” = (0,0,O)T and
v1 = (c/y, 0, O)=.
(2.1)
The Jacobian derivatives of dynamical system (1.3) at equilibrium points are:
(2.2)
where subscripts 0 and 1 indicate that they are obtained characteristic polynomial corresponding to Jo is found to be
at v” and v’, respectively.
i13+&-c=o.
The
(2.3)
The eigenvalues are then given by, (2.4) where II = 3[(27c2 + 4~u~)l’~+ 3fi 2 = 3ti
~1~‘~- 31’36u3cr,
[(27c2 + 4&3)“2 + 3fi
A = [(27c2 + 4~~)“~ + 3ti
~1~‘~+ 35’662’3u,,
~]~‘~3”~6”~.
(25)
When the wave celerity vanishes, then the eigenvalues simplify to A0= 0,
A1= -i(urf2,
A-, = id”.
(2.6)
While c > 0, eigenvalue )Lobecomes a positive real number, and real parts of A, and A_, are both negative. This situation reverses when c < 0 indicating that Hopf bifurcation occurs in the subspace spanned by the eigenvectors corresponding to 3L1and il-,. Sign change in & points out
Travelling waves and chaos in the KSS model
595
that simple bifurcation occurs in the corresponding subspace. Similar behavior has been observed near the equilibrium point v*, although it differs in sign from the latter. Existence of Hopf bifurcation in a system ferrets out the location of periodic orbits [6]. Our aim will be to determine these orbits in the next two sections.
3. NORMAL
FORM
ANALYSIS
The singularity at c = 0 has codimension two (for general codimension two bifurcations see [6]). The normal form analysis has a paramount importance at codimension two bifurcations, and our task is now to perform this analysis. Brjuno normal form analysis is general in character [7] and has the advantage of being amenable to computer applications. Thus, it has been preferred to apply it to the present problem. The prominent features of this analysis has been exploited in [5,7]. The similarity transformation 11 S= 0 ia -a2 [0
1 -iu _a2
1
(3.1)
casts the original system (1.3) into f, = Lc& + L& + a+w,
+ bvl,,,gw+..
(3.2)
Here, a = u112 and Y, i, I, m, p =0, 1, -1. In (3.2) we have complied to the notation introduced by Brjuno [7]. The relations between the original dynamical vector v and the new vector x are defined to be x = S’v and v = Sx. The matrix Lo appearing the right hand side (r.h.s.) of (3.2) is a diagonal matrix of eigenvalues 3L,. The matrix L” is given by -2 1
-2 1
-2 1
[ 1
1
1
v=-$
The coefficients of quadratic below: a Olm
=
--
and cubic monomials
Y 2,
a
b CKlOl=;b 6.
b 1111
a 11~
=
a-th
b 0011
=
a
%m~,
26. 1,
=
(3.3)
appearing
on the r.h.s. of (3.2) are given
5 b 0111
ho-m= born
b-l = hem h-l-l = ho-l,
b- limp -bun, -
(I,
=
1 .
= a36. 1,
boll-1
= booI
ho-, = -ho, b,.,.,-,= -h,,,
m, p = 0, 1, -1).
(3.4)
The coefficients a,+, and b,,,,,,, are symmetric with respect to their subscripts. According to fundamental theorem of Brjuno [7] there exists a normalizing transformation XV=YV + ~%lmYIY??l+ WLrnpY,Yt?lYp + ~Y”IrnprYrYl?lYpYr
(35)
which casts the initial dynamical system (3.2) (with c = 0) into normal form equations ?V = &Yv + %‘vI~Y,Y,
+ ~kz,Y/Y~Yp
+ =%,rY/Y,Y,Yr
(3.6)
G. ijNAL
596
and E. S. $LJHUBi
which involves only the resonant monomials satisfying the resonance condition: +.
(A, Q) =&q, where q,, are the exponents 91,
. . + A,q,
= 0
(3.7)
of monomials appearing in (3.6) with the following property
. . . , qv-1,
. . . , 4v
qv+1,
2
0,
qv
2
-1,
iqi’-l.
[7] (3.8)
j=l
In (3.5) and (3.6) coefficients of monomials are also assumed to be symmetric with respect to their subscripts. Substituting (3.5) into (3.2), having used equation (3.6) in this identity so obtained and retaining terms up to fifth order a new identity can be attained. After symmetrizing the coefficients involved in the latter, following identity can be found [7,8]:
a vmr4lY~YmYpYr +i
(bvjpdm+ bvjrdkp+ bvjlmaj,r + bvjmpd+ bvjlp$mr + bvjmrdp)Y~YmYpYr (3.9)
Summation convention applies to repeated be useful in the forthcoming calculations: A vim=
1 0
if A, = A, + A, ifA,#&+A,’ A vlmpr
1 =
(0
indices in identity (3.9). Following symbols will
A vrmp
=
1
if A, = A, + a, + A,
0
if A, # il, + A, + hp
if A, = h, + il, + A, + il, if A, # AI+ A, + LP + A, ’
(3.10)
Comparing the terms with equal powers at both sides of the identity (3.9), using the symbols introduced above and noting that only the resonant monomials should appear in the normal form equations, complex valued coefficients of the monomials involved in (3.5) and (3.6) can be obtained [7]: (3.11-l)
(3.11-2)
Traveliing waves and chaos in the KSS model
tit
+
~,&I
+
%TtdLp
+
(y,lp&r
597
+
%mr
dIp )}
(3.11-6)
The coefficients (3.11-1 to 6) have been computed on IBM 4381 and they are given in the Appendix. Quadratic, cubic and quartic resonant monomials which will contribute to normal form equations are given in the following in accordance with 191:
(Z)>
(Y~l)?
(“E’)
(Y+)?
(Y~~-l)y
(,I,)?
(,I,)
(,,;J (3.12)
Here, we will only deal with terms of up to fourth order. In Section 5, fourth order terms will be taken into consideration. By utilizing the coefficients given in the Appendix the normal form equations of (3.2) can be written as Y 2 Yo= -;;zYo -2&i
+ O(4)
Yi = iaYl+ -$YoY, f .
Y-l = -my-, ES 30:5-c?
Y
-t pay-1
3( y” - &X4) 2a5
-
iY”,Yt-
3( y2 - csa”) 22
y2 +
g&a4 iY;Y-r + O(4) h5
iYiY_I +
y2 + 9&a4 iy1y2_1+ O(4). 6~’
(3.13)
G. tiNAL
598
and E. S. $UHUBi
Next, we consider the case in which c is a small constant different from zero. Since the normal form analysis of the nonlinear terms has been already performed, transformation sought for the linear term due to wave celerity should not produce changes in nonlinear terms of the normal form equations. Furthermore, this transformation is also desired to allow linear terms to be located on the main diagonal of the normal form equations to alleviate the further analysis. Thus, the normalizing transformation in question becomes (3.14)
x = Y + V(Y) + G(c)Y
where, q(y) stands for the nonlinear terms and G(c) is of order O(c). Normal form equations then take the form [5]
where perturbation
jt = LOy+ L=‘y+ g(y)
(3.15)
L” = LOG- GL” + L”.
(3.16)
matrix L” is
Here, g(y) stands for the nonlinear terms and Lo is a diagonal matrix with eigenvalues A,. Since the perturbation matrix L” is desired to be a main diagonal matrix, upon granting this restriction, the relation (3.16) allows us to determine the elements of the both matrices, namely L”’ and G: Lc’=$
:I
;J,
G=$
[ -ii
-;
ii].
(3.17)
Finally we are led to the following normal form equations: .
c
Y
Yo=~Yo-;;Tyo
2
-SYiY-i+
O(4)
It is easily seen that normal form equations (3.18) has the following property
where - denotes the complex conjugate. variables: YO= y0,
Y, =y-1
(3.18-1)
The last relation enables us to introduce
the new
Then we see that the normal form equations
(3.19)
y-, = re+.
Y, = reip,
(3.18) can be rewritten as (3.20-l)
+-&+Y
,2
2u2
@=a+
3(y2 - Sa2) 2a5
YC
(3.20-2)
w
y2 + 9Sa4 r2 3a5
’
Notice that (3.20-l to 3) are free from the angular variable implying that the KSS equation has an internal symmetry near the equilibrium point v” for /3 = 0. Fortunately this property allows
599
Travelling waves and chaos in the KSS model
us to study the dynamical problems on the phase plane. Second step in ~odimension two bifurcation theory is to eliminate the terms of order O(ly, rl*) in (3.20-3) and study the orbits of the phase plane. This will be treated in the next section.
4. PERIODIC,
QUASI-PERIODIC
AND SOLITARY
WAVES
Equilibrium points of the planar system (3.20-1, 3.20-2) are Y” =
(0,w=,
c >T C,o
yL
Y
and
Here, y = (yO, r)= is the dynamical variable vector. Linear stability analysis reveals that the equilibrium point y” is a center with pure imaginary eigenvalues. Equilibrium point y” and y1 are saddle points with real eigenvalues of opposite sign. Phase plane in which we are interested consists of limit cycles around the center and homoclinic loop joining the saddle points. First integral of the planar system (3.20-1, 3.20-2) can be found by standard means: H = i yor2 - y&’ - r4.
(4. I)
The level curves are then defined by the relation ~(y~, r) = S = const. The parameter the value of
S takes
c4
on the center and zero on the homoclinic loop. Hence for a given value of wave celerity and for different values of y, there exist infinite number of limit cycles. Since the frequency of rotation around the center is smaller than the frequency about the y. axis (c/fia’
and a are anticipated as solutions to (1.3). Expressing y. in terms of r and S via first integral (4.1) and then substituting the latter in the planar system (3.20-1, 3.20-2), following equation can be found: 2
(
l/2
_r4+Cr2_~
i=*tr a*
(4.2)
>
4Y2
which yields immediately r1
I r
dp [@Z
_
p2)(p”
_
r;)]l”
=
*$
(*
+
cJ
(4.3)
where (4.4)
600
SUHUBi
and E. S.
G. ONAL
It is a simple matter to convert this integral to an elliptic integral of the first kind [lo]. Having done this, solution to planar system is obtained, after some algebra as y,=~_42r,snecne dn 8 2Y
r = rl dn 8,
e = $ (2 + c,)~,
sgn 8,
(4.5)
where rf - rf 2 rI
4’
and sn 8, cn 0, dn 0 are Jacobian elliptic functions. An expression for the period can also be found immediately by replacing the lower limit of integral (4.3) with r, and then integrating so that T(S) =
z
K(q)
1
where K(q) is a complete elliptic integral. Following limits would be useful in regard to understanding the behavior of solution (4.5) in the limit of homoclinic loop and the center: (4.7) Homoclinic loop is made up of two heteroclinic orbits. One of them emanates from the saddle point y’, follows its unstable manifold and ends on the saddle point y” following its stable manifold. Along with the lines mentioned in the periodic case, solution associated with this heteroclinic orbit is found to be: r=$sech8,y,=&(l-sgnBtanh8). We are also interested
e=$(z+C2).
in the values taken by the latter at z ---, w lim y. = f,
lim y. = 0,
,?---m
lim r = 0.
*--
(4.9)
z- *m
The other heteroclinic orbit completes rest of the homoclinic loop by emanating from y” and ending at y*. Similarly solution leading to heteroclinic orbit can be written as
Y yo = 1 +
r = 0,
(4.10)
e(cla*)r-d,
where d, is an arbitrary constant. If z ---, i-m, limits of this solution become lim y. = 0, L---m
lim y. = c. Z--rY
(4.11)
Original dynamical variable v can be expressed in terms of the variables yo, r and Q, by employing the normalizing transformation (3.5) together with the coefficients in the Appendix as follows: v=y,+ +&
(
2c-
3r
>
'Yy, ;;Jsinv+
r2 sin
2q +
[
-
9( y2 - 6a4)
19y* - 186a4 9a6
2a6
yZ+
135 y2 + 276a4 18a6 2
y,r2cos2q-2$yi-
r*+
1
2 r cos Q,
y2 + 96a4 36a6 r3cos3g,
(4.12)
601
Travelling waves and chaos in the KSS model
where v = vl, rp = az + Co and C,, is an arbitrary following expression can be found:
+2y2q4r: sn2 8 cn2 e sgn e a6 dn2 8
Substituting
(4.5) into (4.12)
q2r1 1
sn e cn e 3 y2q2rf sn 0 cn 8 dn 8 sin rp + sgn e a3 dn 0
+ _9(Y2 - 6a4) c2 + 135~’ + 276a4
r?dn28+2
18a6
8y2a6
constant.
1 1
q2rIsn8cn8 dn 8
I
q2r:sn8cn8
r,dntIcosg,
cos~+Lr~dn2Bsin2~ 3a3 *
2a6
19~’ - 186a4 cdn0 + rf dn 6 cos 29~ - sgn 8q2rI sn e cn e 9a6 2Y - y2 + 96a4 r: dn3 tI cos 391. 36a6
(4.13)
Here sgn 8 stems from (4.2). At the equilibrium point 3 of the planar system, values taken by the parameters are
s=c4
q=o
64 y4’
and
Using these values in (4.14) and assuming that (cl << 1, (4.13) simplifies to v=;
cl (
1 2+/zcoscp.
(4.14)
)
Thus, the equilibrium point y” of the planar system (3.20-1, 3.20-2) corresponds to a limit cycle in three-dimensional phase space (other coordinates being it and ti). Hence the predictions of the linear stability analysis are verified. While s<-
C4 64f
solution (4.13) and its first and second derivatives form a torus in three-dimensional phase space. This merely indicates that Neimark bifurcation is occurring. These are illustrated by
“3
J0
"1
,..’ /
.^
_.-.\ ,_.-
x0’
------+
(’
,,,:/
-b
._: ‘\.__._.
Hopf bifurcation
_.
-
”
Neimark
bifurcation
Fig. 1. Bifurcations from the equilibrium point and periodic orbits.
602
G. tiNAL
and E. S. SUHUBi
plotting the solution (4.13) and (4.14) in Fig. 1. Thus we have identified the first two bifurcations before the onset of turbulence. Hopf bifurcation followed by a Neimark bifurcation is known to be “Basic Sequence” [ll]. The latter is a prelude to the scenarios given in Refs [12-141. Right after the Neimark bifurcation infinite-period bifurcations occur. This can be discerned by considering the solution (4.13) together with expression (4.6). Langford [1.5] has demonstrated the existence of infinite-period bifurcations numerically in the context of codimension two bifurcations, The same author also expounded that Neimark bifurcations result from the interactions between simple and Hopf bifurcations. As it can be easily seen from (4.13) elliptic functions are modulated by trigonometric ones. In a sense, elliptic and trionometric functions play the role of large and small scale eddies, Amplitude of the smaller eddies are controlled by the parameter y which appear in front of the convective term of equation (1.1). This solution implies the existence of energy transfer from large eddies to smaller ones. Naturally, the coefficient 6 pertaining to negative viscosity term appears in solution (4.13). This can explain the propagation of travelling waves without dissipation on the viscous fluid. All of these agree well with interpretations given in Refs [l, 61. Substituting solution (4.8) into expression (4.12), solitary wave solution to the system (1.3) can be easily found. After having done this and taking limits we obtain lim u = 0. 2-m
(4.15)
When
first quantity in (4.15) vanishes. At this value of the wave celerity, solitary wave with many humps should appear (see Fig. 2). Similarly, solitary wave with single hump can also be obtained from (4.10) as it can be seen from Fig. 3. We have not pursued the analysis further to see whether these are also solitons or not.
Fig. 2. Solitary wave with humps plotted for the parameter values of y = 0.5, CY= 0.75, 6 = 0.25 and c = -0.459.
Travelling waves and chaos in the KSS model
V
603
-0 y
] -0.2 -
-0.4 -
-0.6-
,
I
I I I I,
-30
I,,
,
,
,
,
,
,
,
,
-20
,
,
I,,
,
,
,
,
,,
-10
,
,
,
,
,
,
,
+
,
,
0
,
,
,
,
,
10
,
,
,
,
,
,
20
Fig. 3. Solitary wave with a single hump plotted for the parameter values of y = 1, (Y= 0.25, 6 = 0.5 and c = -0.088.
5. MELNIKOV
ANALYSIS
What happens to the trajectories of the planar system (3.20-1,3.20-2) under the perturbation of higher order terms? In order to be able to provide an answer to this question, fou~h-order normal form analysis together with Melnikov analysis have to be carried out. By making use of the results obtained Section 3 and the coefficients given in the Appendix, normal form equations are found as
@=
a + WIYO, r12).
(5.1-3)
Resealing the variables and the celerity in (5.1-1 to 3) with a small parameter the following manner r= 6,
yo = C&
v,
2’ = Z/E
c = v&J,
E (0 < E << 1) in
(5.2)
yields * = f(w, p) + &g(W)*
(5.3)
Here, dot denotes the derivative taken with respect to z’ and w = (u,
~m~nents
ZI)T,
f=
vm,
f2w)=,
g =
tfim%
g*(w))‘.
of the vector field have the form fi(w) =$
(-J4U + 2yuu),
gr(w) = UrUV3+ b&,
f*(w) = - (/z?J - yu2 - 2yu2) :2
(5.4-l)
g,(w) = u2u4 + b2u2v2+ c2u4
(W-2)
604
G. fiNAL
and E. S. $UHUBi
where b = y(82y2 + 816a4)
a = _Y(17Y2 - 246a4) 1 4aX ’ b = 9Y(Y2 - 26a4) 2
2aX
1 c2= -
’
l&&8
y(27ly*+
y(4y2 - 36a4) ’
‘* =
ax
3786a4)
18aX
(5.4-3)
’
Melnikov function is proportional to the distance between stable and unstable manifolds inherent in homoclinic and heterclinic connections [9]. Following Salam [16], Melnikov integral for dissipative systems can be written as MC%) = [= f(q”(z - zoo))A g(q”(z - r,))exp[ -I--CC
tr D&q’(s))
ds] dz
(5-5)
Here, q” stands for the unperturbed trajectories, A denotes the wedge product (f/\g= fig, -figI) and tr D,o is the trace of Jacobian derivative along the unperturbed trajectory. Unperturbed periodic orbits of the phase plane corresponding to (5.3) can be easily obtained from (4.5) by replacing u with r, y. with u and c with p. Trace of the Jacobian derivative taken along the periodic orbits of the unperturbed system (5.3 with E = 0) is tr D,of(s) = sgn 0
Yq2rI sn 8 cn 8 dn8 ’
8 = 4 r,s. a
(5.6)
After having evaluated the exponential expression in (5.5), we can convert (5.5) into a line integral. Then by making use of the Green’s identity, we are led to the following integral u3 (2b, + 6,) II v + (4az + 2u,) ; u3} dv] du
(5.7)
where 1, =;+$
(
4y2s II2
+44y=u*--
u2 1 ’
P*-4y2u*__
12=g-$
(
4y2s 1’2 u* 1
’
rI and rz are obtained from (4.4) by replacing c with ~1,When the integral in braces is evaluated we obtain Z, + (B - 2A)Z, - 2ASZs3
(5.8)
where the following quantities and integrals used A=-
2a2 + al
2 Z, =
I
’
r’ (p2u2 - 4 y2u4 - 4~~s)“~ du,
B = 2b1 + b2
Z, =
(5.9-l)
” u2(p2u2 - 4Y2u4 - 4Y2S)“* du, IQ
*I *u2 - 4y*u4 - 4Y*S)l” Z, = r2 (P du. I r2 u2
(5.9-2)
After these integrals have been evaluated by making use of the integrals Melnikov function for this case is found to be 1270
gP2Y2SK+
270~’ - 355 r2SE 24
1
I
given in [lo],
(5.10)
605
Travelling waves and chaos in the KSS model
where E(q) and K(q) elliptic integrals of the first kind and second kind, respectively. For a fixed value of ~1, Melnikov function increases while the parameter S approaches to zero. Equilibrium point
(u, v)‘=
(*y, gT
of the unperturbed system (5.3 with E = 0) is a center. However, higher order monomials this equilibrium point switches to UC-
CL
ti
y + ’ 2304a6y3
2fi
under the perturbation
of
(179y* - 9186a4)p3 + O(s*)
ti P ~ (71~~ - 2976a4)p3 ’ = G + E 288a6 y3
+ O(E*).
(5.11)
Linear stability analysis near the equilibrium point reveals that it is indeed an unstable focus. Both analyses agree on the destruction of torus under the perturbation of higher order terms. Similar calculations for the homoclinic loop showed that Melnikov function has the form (5.12) Since Melnikov function does not have a simple zero in both cases, transversal intersection of stable and unstable manifolds does not occur as it has been anticipated. Thus it can be inferred that the homoclinic bifurcation should occur as a result of the perturbation of higher order terms [9].
6. TRANSITION
TO CHAOS
The consequence of the homoclinic bifurcation in the planar system and the existence of horseshoe mapping and chaotic dynamics heralded from it can best be viewed by constructing Silnikov’s example for the present problem. Following Guckenheimer [9], let us construct a cylinder of radius p. and height of {I in the neighborhood of the origin which is a saddle point (see Fig. 6.5.2 of [9]). L inearized system corresponding to (1.3) is obtained as
(6.1) where
By utilizing the similarity transformation (6.1): v=Sr,
(3.1), the following relations are defined in obtaining
r = s-iv,
r = (5, 9, 0’.
Solutions of (6.1) rules the flow near the saddle point:
[I [ E
77 (z) = I;
euZ[cos /3z E(O) - sin /3z n(O)] euZ[cos /3z E(O) + sin /3z n(O)]
e”’f (0)
. I
(6.2)
G. ONAL and E. s. SUHUBi
606
When trajectories pass by the origin, 15‘1 increases decreases. Let us define the surfaces Za and Z, by
and radial coordinate
p = (6” + $)“’
z:o = ((5 11, 5‘) I 5” + r2 = P:,>0-c 5 < 5‘1) 2, = ((5, rl, I;) I E2 + v2 = Pit, c = c-1‘0).
(6.3)
It has been assumed that Co and Ci are in the neighborhood U where the flow is linear [9]. The mapping 0 :X0+ Zi can be found by substituting the time of flight z = A-’ In(~l/~(0)) into (6.2). Thus, we obtain Vz = ((51/f)““](cos Here o = /In-’ ln(
w)5’ - (sin
whl, (1;JC;)“‘*[(sin ~15 + (~0so)rl, fd
By introducing the new variables can be obtained [9]:
(6.4)
f = p0 cos Y and rl = p. sin v, the
@(v, 5) = (~~(5;~/0”‘” cos(v + w), ~~(5;~/0”” sin(v + 0))
(6.5)
0 maps a vertical segment Y = constant of &, to a logarithmic spiral surrounding the 5; axis and lies in Xi. Stretching of the vertical segment, wrapping it around the spiral can easily be seen by calculating the determinant of Jacobian derivative [9]:
(6.6) Stretching ratio s = (plpo)/( 1;/5;1) can be found from the last quantity: s=-
Cl 5
1+0/A .
(>
(6.7)
Stability of the homoclinic orbit in the neighborhood of the saddle point is being determined by the stretch ratio and it is stable for s < 1 [6, 91. While k > 0 and o < 0, linear stability condition is found to be
14’ 14. When the situation is reversed,
linear stability is determined
14 ’ 14.
(6.8) by (6.9)
For the system (1.3) la/Al = l/2, thus it has a homoclinic orbit which is stable according to (6.9). For the right travelling waves near the v” and for the left travelling waves near the VI, homoclinic orbit can be observed. What follows from the analysis is that the vertical segment {o = constant, < E (0, &,)} on Z,, are mapped into a spiral of maximum radius po(l;l/l;)“” (C&J and the lengths of these segments are stretched by an amount which is unbounded as &,-,O [9]. Let us now consider a region A = {(p, w, 5) 1p = po, joI< q, 0 < I; -=cQ} near a point on the surface X0. The region A will be stretched under the mapping 0, yet its maximum radius will shrink when the condition (6.8) is fulfilled. According to Silnikov’s theorem when the condition is satisfied the region A fl (q!~0 0) (A) contains a horseshoe [9]. Horseshoe map has the properties of folding, stretching and mixing. Cantor like invariant set of the horseshoe map contains (a) a countable set of periodic orbits of all periods; (b) an uncountable set of nonperiodic motions; (c) a dense orbit; (d) the periodic orbits are all of saddle type and they are dense in the invariant set [9]. Thus at the parameter values close to the one which homoclinic orbit appears, chaos should prevail. Analytical results obtained hitherto will be verified by computing Poincare surface of section and LCEs. Initial conditions needed to construct surface of section obtained from the planar
607
Traveliing waves and chaos in the KSS model
system for different values of the parameter S (0 c S < c4/64y4). With these initial conditions dynamical system (1.3) has been integrated by employing fourth-order Runge-Kutta for the parameter values of c = -0.5, y = 1, cy = 2 and S = 0.5. Poincare surface of section can be seen from Fig. 4. Poincare map corresponding to the latter has a fixed point located at the center of section which has been obtained by taking
s=c4
64y4
This fixed point corresponds to a limit cycle in three dimensional phase space. By varying the parameter S closed curves su~ounding the fixed point appear. These are orbits having two incommensurate frequencies, which is so called quasiperiodicity, and they correspond to tori of different sizes in the phase space. One of the salient features of these orbits is their winding number. It is the ratio of the number of turns around the torus in each of the two directions [17]. Since the frequencies of nonlinear waves depend on the amplitudes, the winding number of orbits varies ~ntinuously for different initial conditions. Here, parameter f of the planar system is controlling the winding number. For the rational winding numbers orbits close on themselves after a finite number of iterations. For the irrational values of the winding number, the orbit covers the entire closed curve densely. RAM theorem (after Kolmogorov, Arnold and Moser) predicts that the orbits with sufficiently irrational winding number persist under small but finite ~~urbation. These orbits are also called RAM curves 191. Thus we infer that the dynamical system (1.3) is locally conservative. Hitherto all the properties of analytical solution (4.13) agress very well with the numerically obtained Poincare section. Formation of separatrix with seven hyperbolic fixed points and islands inside of it appearing in Fig. 4 can be understood by considering Moser twist map along with the Poincare-Birkhoff theorem [9]. When the frequency ratio of large scales (elliptic functions) and small scales (trigonome~c ~n~ions) is rational, perturbation of higher order terms switched on by the
E!
-0.4
y
-0.8
1
-1.0
I,,,,> -0.6
f,,,,,,,l,,,,,,,, -0.5
/ / -0.4
/
/ ,,
,,,,,(
-0.3
/i,,,
/,
-0.2
Vl
Fig. 4. Poincart? surface of section.
,
,,,,,,,,/,,,,
-0.1
,,,
II -0
608
G. ijNAL and E. S. SUHUBi
06
;;
-
04-
> E
3
E -I
0.2
-
-0
-
400
800
1200
1600
2000
2
Fig. 5. Positive LCE.
magnitude of celerity has an important effect on the periodic orbits. In other words resonance occurs and energy is exchanged between large and small scales. Direction of the energy transfer is determined by the relative phases and the parameters y and 6. Hyperbolic fixed points of the perturbed map and orbits connecting them from the separatrix (resonant curve). Islands inside of the separatrix are formed around the elliptic fixed points of the perturbed map. As a result 7: 1 resonant curve and islands appear on the surface of section. As a principle for every rational winding number there corresponds an island structure. Very close to the parameter value S = 0, celebrated chaotic trajectory showed up as having 0.6
-
I
Gi > 5
0.4
-
0.2
-
-O-
cu w Y
-0.2
-
-0.4
-
-06
1,,,,,,,,,,,,,,,,~,,,,,,,,,~,,,,,,,,,,,, -400
0
,,,,,,,,, 400
800 L
Fig. 6. Zero LCE.
1200
,,,,,,,,,, 1600
2cQo
Travelling waves and chaos in the KSS model
609
L
Fig. 7. Negative LCE.
fuzzy crown like (outer orbit in Fig. 4) structure. Homoclinic points resulting from the transversal intersections of the stable and unstable manifolds near the hyperbolic fixed points are prominent. It is worth noting that chaos sets in before the limit of homoclinic orbit. This is accounted for the parameter values taken here (i.e. c + u”/fi). Occurrence of homoclinic bi~rcation on the onset of the chaos is evident from Fig. 4. The chaotic trajectory observed here is a transient one and it escapes from the attractor after certain period of computer integration time. The closer the values of parameter S taken to the value of S for which separatrix appears the more time elapses before the trajectory escapes. Computation of LCEs are performed to check whether the chaotic orbit survives on the strange attractor or not. The LCEs are the average exponential rates of the convergence or divergence of the nei~borhood trajectories in the phase space fI8). If one of the LCEs is positive and their sum is negative, thus the chaotic trajectory lies on the strange attractor [19]. As seen from Figs 5-7, the LCEs for the chaotic trajectory are A1= 0.0415 bit/s. Since the sum of these exponents is negative (--0.00002) and one of them is positive, chaotic trajectory lies on the strange attractor. Hausdorff dimension of the latter is calculated by Kaplan-Yorke conjecture 1191 and is found to be approximately 2.9976. REFERENCES B. NICOLAENKO, &cl. Phys. {Proc. Suppl..) 2,453 (1987). G. I. SIVASHINSKY, Physicu 17D,243 (1980). M. C. DEPASSIER and E. A. SPIEGEL, A.rrron. I; &t(3), 4% (1981). D. K. CAMPBELL, Los AIqos Sci. Spec. iss. 218 (1987). G. UNAL and E. S. SUHUBI, Znt. J. Engng Sci. 30,579 (1992). J. GUCKENHEIMER, SIAM J. Math. Anal. xS(l), 1 (1984). V. M. STARZHINSKI, Applied Methods in the Theory of Nonlinear Oscillations. Mir, Moscow (1980). G. UNAL, A group theoretical approach to turbulence (Ph.D. thesis in Turkish). Istanbul Technical University, Istanbul, Turkey (1991). J. GUC~N~I~R and P. HOLMES, Nonlinear Osc~fatio~, ~y~mical Systems and Bi~rcation of Vector Piehis. Springer, New York (1983). I. S. GRADSHTEYN and I. M. RYZHIK, Tables of integrals Series und Products. Academic Press, New York (1971). R. H. ABRAHAM and J. E. MARSDEN, Foundations of Mechanics. Benjamin-Cummings, Reading MA (1978).
610
G. ANAL and E. S. SUHUBI
[12] [13] [14] [15]
L. D. LANDAU and E. M. LIFSHITZ, Fluid Mechanics. Pergamon Press, Oxford (1959). D. RUELLE and F. TAKENS, Commun. M&z. Phys. 20, 167 (1971). S. NEWHOUSE, D. RUELLE and F. TAKENS, Commun. M&h. Phys. 64, 35 (1978). W8:j.LANGFORD, in Chaos, Fractals and Dynamics (Edited by P. FISCHER). Marcel Dekker, New York
[16] [17] [18] [19]
F. M. SALAM, SZAM J. Appl. Math. 47(2), 233 (1987). J. M. GREENE, Physica 18D, 427 (1986). G. BENETTIN, L. GALGANI, A. GIORGILLI and J. M. STRELCYN, Mecca&a 15,9 (1980). A. WOLF, J. B. SWIFT, H. L. SWINNEY and J. A. VASTANO, Physica 16D, 285 (1985). (Received 31 August 1991; accepted 11 September
1991)
APPENDIX The coefficients of resonant monomials contributing which are not vanishing will be reproduced here:
K
,oo,= Y’ F,’
y2 - 96a4
h4
K,,,., = -51, 52oooo=
K-,-,-I, = -K,,,.,,
%x1-1 =
82yZ+ 816a4
216us
Q”.,.,,, = -
to the normal form equations are given below. Only the ones
18a
4yz - 36a4 aR yp
Q-,0,x,-, = f&u,,,
yj
271 y2 - 3786a4 108~s ”
K.,w-I = -K,nn, Q,,,,
= -
Q,“.,.,,
a,,=--lli
a loo=&i,
’
ab,, = “lx&,
a10.1= (r,,/z
~11,= -I,,
al-l-1 = %0/3, =11--l
&n-1=
=
o-,00= a-,-,-*
-@I,.,,
13y* + 96a4 Q
a3
F
=
w-,-,0
8*-1-m= B ,--1-1-1 =-q$Y, B- ,110 --B ,-,-*o~
~,I., = ~ICKI
ooo,
-(Y,fJo,
B
L-111 = B,-,-,-,r
a- ,,, =
Y2 -___
r%&,,
5 y2 - 96a4
pr_,,=-$>
=
’
= B”,,&,j
IL-11
T
-ff1-I-,
&14
u6
= BIX,“?
27a6
are:
cu,.,.,= acd2
ff_ 10, = -o,o-,?
-o,m,
= Q,“,,.,
52
And the coefficients of monomials appearing in the normalizing transformation %o,=,s’,y.
17y2 - 24Sa4 l&, y,
= -8,,,o>
B-,00,=-/%c#J,, /L,,,-,
= -7yy&:ba4’
8-I-1-10
= 8,,1”P
Pm,-,-,-,
= /$I,,.
G42)
Int. J. Engng Sci. Vol. 30, No. 5, pp. 593-610, 1992 Printed in Great Britain. All rights reserved
TRAVELLING WAVES AND CHAOS IN THE ISOLMOGOROV-SPIEGEL-SIVASHINSKY MODEL G. ONAL
and E. S. SUHUBl
Istanbul Technical University, Faculty of Sciences, Maslak 80626, Istanbul, Turkey Abstract-Kolmogorov-Spiegel-Sivashinsky (KSS) equation has been derived to model large-scale structures arising in two-dimensional turbulence [Null. Phys. (Proc. Suppi.) 2, 453 (19t37)j. In the present study, we have obtained its periodic, quasi-periodic and solitary wave solutions analytically to a certain degree of approximation by discarding the linear damping term. Normal form analysis succeeds to capture Hopf, Neimark and infinite-period bifurcations before the onset of chaos. Melnikov analysis has been carried out to identify the homoclinic bifurcation. Transient spatiotemporal chaos has been observed, by constructing the Poincart-section and it has been identified by computing Lyapunov Characteristic Exponents (LCE).
1.
INTRODUCTION
The spontaneous formation of large-scale structures in a viscous fluid driven at small scales is one of the challenging problems of hydrodynamic turbulence. The latter ensues from the corresponding small-scale flow as a manifestation of the long-wave instability [2]. Kolmogorov flow and large-scale turbulent solar convection possess these features [l, 21. The former is a two-dimensional viscous fluid flow subjected to a unidirectional external force field periodic in one of the space coordinates [2]. KSS equation is capable of modeling large-scale structures in flows mentioned above. A specific form of the KSS equation can be written as $* + BQ, + Y& + (a - 3W,2MJXX+ #XXXX =0
(1.1)
where (Y, fi, y and 6 are physical positive constants. Sivashinsky [2] derived the equation (1.1) (with LY= 2) for a unidirectional large-scale flow parallel to a direction where the effective viscosity of the secondary flow is negative in the context of Kolmogorov flow. In this case, +(x, t) is the resealed large-scale stream function obtained after subtracting the mean periodic field component, x is the resealed preferential direction of negative viscosity and E is the resealed time variable. Spiegel et al. [3] derived an equation similar to (1.1) in multiple-scales analysis of a compressible solar convection layer zone. Their equation also includes a hydrodynamic blow-up term and a term carrying non-Boussinesqian compressibility effects. Circumventing the latter terms leads to equation (1.1) with #(x, t) representing the large-scale horizontal fluctuations of the temperature field after subtraction of the mean field and x being a large-scale horizontal variable [l]. The role played by each term in equation (1.1) has been expounded in references [ 1,4] (see also [5]). In the present study, linear damping term will be annulled by setting the coefficient /3 = 0. This is merely due to the fact that travelling wave solution found in our previous work [5] does not remain bounded for & = 0. In the same work, we were able to obtain the following dynamical system by resorting to group theoretical means: zi1= u2 l&=l+ iiig= l&q . a4
=
-j?u~
+
clip
-
cm3
+
3&&
- yu$
(1.2)
and z=xct, and c is the celerity of waves. Solutions to this dynamical where u 1=#(z) system are known to be travelling waves which they retain their shape at constant celerity. 593
G. iiNAL
594
and E. S. .$UHUBi
When c = 0 and /3 = 0, eigenvalues of the Jacobian derivative gradient matrix of the right hand side of (1.2) become a double zero and a pair of complex conjugate pure imaginary numbers. This situation introduces an extra degeneracy to the problem. The latter can be circumvented by the following change of variables u2 = ul, u3 = u2 and uq = us. Thus, the dynamical system which we are concerned becomes,
ti,
=
CV1 -
&V2 + %J:V2
-
J’V:.
(1.3)
Our main objective here is to study the local properties of the dynamical system (1.3) in the neighborhood of an equilibrium point. In Section 2, we identify the stability of the equilibrium points by employing the linear stability analysis. Normal form analysis has been performed in Section 3. In Section 4, periodic, quasi-periodic and solitary wave solutions to (1.3) are constructed. Behavior of the solutions obtained in Section 3 under the perturbation of higher order terms are investigated by resorting to Melnikov analysis in Section 5. In Section 6, transition to chaos has been discussed and numerically demonstrated by computing Poincare sections and LCEs.
2. LINEAR
STABILITY
Dynamical system (1.3) has two equilibrium
ANALYSIS
points located at
v” = (0,0,O)T and
v1 = (c/y, 0, O)=.
(2.1)
The Jacobian derivatives of dynamical system (1.3) at equilibrium points are:
(2.2)
where subscripts 0 and 1 indicate that they are obtained characteristic polynomial corresponding to Jo is found to be
at v” and v’, respectively.
i13+&-c=o.
The
(2.3)
The eigenvalues are then given by, (2.4) where II = 3[(27c2 + 4~u~)l’~+ 3fi 2 = 3ti
~1~‘~- 31’36u3cr,
[(27c2 + 4&3)“2 + 3fi
A = [(27c2 + 4~~)“~ + 3ti
~1~‘~+ 35’662’3u,,
~]~‘~3”~6”~.
(25)
When the wave celerity vanishes, then the eigenvalues simplify to A0= 0,
A1= -i(urf2,
A-, = id”.
(2.6)
While c > 0, eigenvalue )Lobecomes a positive real number, and real parts of A, and A_, are both negative. This situation reverses when c < 0 indicating that Hopf bifurcation occurs in the subspace spanned by the eigenvectors corresponding to 3L1and il-,. Sign change in & points out
Travelling waves and chaos in the KSS model
595
that simple bifurcation occurs in the corresponding subspace. Similar behavior has been observed near the equilibrium point v*, although it differs in sign from the latter. Existence of Hopf bifurcation in a system ferrets out the location of periodic orbits [6]. Our aim will be to determine these orbits in the next two sections.
3. NORMAL
FORM
ANALYSIS
The singularity at c = 0 has codimension two (for general codimension two bifurcations see [6]). The normal form analysis has a paramount importance at codimension two bifurcations, and our task is now to perform this analysis. Brjuno normal form analysis is general in character [7] and has the advantage of being amenable to computer applications. Thus, it has been preferred to apply it to the present problem. The prominent features of this analysis has been exploited in [5,7]. The similarity transformation 11 S= 0 ia -a2 [0
1 -iu _a2
1
(3.1)
casts the original system (1.3) into f, = Lc& + L& + a+w,
+ bvl,,,gw+..
(3.2)
Here, a = u112 and Y, i, I, m, p =0, 1, -1. In (3.2) we have complied to the notation introduced by Brjuno [7]. The relations between the original dynamical vector v and the new vector x are defined to be x = S’v and v = Sx. The matrix Lo appearing the right hand side (r.h.s.) of (3.2) is a diagonal matrix of eigenvalues 3L,. The matrix L” is given by -2 1
-2 1
-2 1
[ 1
1
1
v=-$
The coefficients of quadratic below: a Olm
=
--
and cubic monomials
Y 2,
a
b CKlOl=;b 6.
b 1111
a 11~
=
a-th
b 0011
=
a
%m~,
26. 1,
=
(3.3)
appearing
on the r.h.s. of (3.2) are given
5 b 0111
ho-m= born
b-l = hem h-l-l = ho-l,
b- limp -bun, -
(I,
=
1 .
= a36. 1,
boll-1
= booI
ho-, = -ho, b,.,.,-,= -h,,,
m, p = 0, 1, -1).
(3.4)
The coefficients a,+, and b,,,,,,, are symmetric with respect to their subscripts. According to fundamental theorem of Brjuno [7] there exists a normalizing transformation XV=YV + ~%lmYIY??l+ WLrnpY,Yt?lYp + ~Y”IrnprYrYl?lYpYr
(35)
which casts the initial dynamical system (3.2) (with c = 0) into normal form equations ?V = &Yv + %‘vI~Y,Y,
+ ~kz,Y/Y~Yp
+ =%,rY/Y,Y,Yr
(3.6)
G. ijNAL
596
and E. S. $LJHUBi
which involves only the resonant monomials satisfying the resonance condition: +.
(A, Q) =&q, where q,, are the exponents 91,
. . + A,q,
= 0
(3.7)
of monomials appearing in (3.6) with the following property
. . . , qv-1,
. . . , 4v
qv+1,
2
0,
qv
2
-1,
iqi’-l.
[7] (3.8)
j=l
In (3.5) and (3.6) coefficients of monomials are also assumed to be symmetric with respect to their subscripts. Substituting (3.5) into (3.2), having used equation (3.6) in this identity so obtained and retaining terms up to fifth order a new identity can be attained. After symmetrizing the coefficients involved in the latter, following identity can be found [7,8]:
a vmr4lY~YmYpYr +i
(bvjpdm+ bvjrdkp+ bvjlmaj,r + bvjmpd+ bvjlp$mr + bvjmrdp)Y~YmYpYr (3.9)
Summation convention applies to repeated be useful in the forthcoming calculations: A vim=
1 0
if A, = A, + A, ifA,#&+A,’ A vlmpr
1 =
(0
indices in identity (3.9). Following symbols will
A vrmp
=
1
if A, = A, + a, + A,
0
if A, # il, + A, + hp
if A, = h, + il, + A, + il, if A, # AI+ A, + LP + A, ’
(3.10)
Comparing the terms with equal powers at both sides of the identity (3.9), using the symbols introduced above and noting that only the resonant monomials should appear in the normal form equations, complex valued coefficients of the monomials involved in (3.5) and (3.6) can be obtained [7]: (3.11-l)
(3.11-2)
Traveliing waves and chaos in the KSS model
tit
+
~,&I
+
%TtdLp
+
(y,lp&r
597
+
%mr
dIp )}
(3.11-6)
The coefficients (3.11-1 to 6) have been computed on IBM 4381 and they are given in the Appendix. Quadratic, cubic and quartic resonant monomials which will contribute to normal form equations are given in the following in accordance with 191:
(Z)>
(Y~l)?
(“E’)
(Y+)?
(Y~~-l)y
(,I,)?
(,I,)
(,,;J (3.12)
Here, we will only deal with terms of up to fourth order. In Section 5, fourth order terms will be taken into consideration. By utilizing the coefficients given in the Appendix the normal form equations of (3.2) can be written as Y 2 Yo= -;;zYo -2&i
+ O(4)
Yi = iaYl+ -$YoY, f .
Y-l = -my-, ES 30:5-c?
Y
-t pay-1
3( y” - &X4) 2a5
-
iY”,Yt-
3( y2 - csa”) 22
y2 +
g&a4 iY;Y-r + O(4) h5
iYiY_I +
y2 + 9&a4 iy1y2_1+ O(4). 6~’
(3.13)
G. tiNAL
598
and E. S. $UHUBi
Next, we consider the case in which c is a small constant different from zero. Since the normal form analysis of the nonlinear terms has been already performed, transformation sought for the linear term due to wave celerity should not produce changes in nonlinear terms of the normal form equations. Furthermore, this transformation is also desired to allow linear terms to be located on the main diagonal of the normal form equations to alleviate the further analysis. Thus, the normalizing transformation in question becomes (3.14)
x = Y + V(Y) + G(c)Y
where, q(y) stands for the nonlinear terms and G(c) is of order O(c). Normal form equations then take the form [5]
where perturbation
jt = LOy+ L=‘y+ g(y)
(3.15)
L” = LOG- GL” + L”.
(3.16)
matrix L” is
Here, g(y) stands for the nonlinear terms and Lo is a diagonal matrix with eigenvalues A,. Since the perturbation matrix L” is desired to be a main diagonal matrix, upon granting this restriction, the relation (3.16) allows us to determine the elements of the both matrices, namely L”’ and G: Lc’=$
:I
;J,
G=$
[ -ii
-;
ii].
(3.17)
Finally we are led to the following normal form equations: .
c
Y
Yo=~Yo-;;Tyo
2
-SYiY-i+
O(4)
It is easily seen that normal form equations (3.18) has the following property
where - denotes the complex conjugate. variables: YO= y0,
Y, =y-1
(3.18-1)
The last relation enables us to introduce
the new
Then we see that the normal form equations
(3.19)
y-, = re+.
Y, = reip,
(3.18) can be rewritten as (3.20-l)
+-&+Y
,2
2u2
@=a+
3(y2 - Sa2) 2a5
YC
(3.20-2)
w
y2 + 9Sa4 r2 3a5
’
Notice that (3.20-l to 3) are free from the angular variable implying that the KSS equation has an internal symmetry near the equilibrium point v” for /3 = 0. Fortunately this property allows
599
Travelling waves and chaos in the KSS model
us to study the dynamical problems on the phase plane. Second step in ~odimension two bifurcation theory is to eliminate the terms of order O(ly, rl*) in (3.20-3) and study the orbits of the phase plane. This will be treated in the next section.
4. PERIODIC,
QUASI-PERIODIC
AND SOLITARY
WAVES
Equilibrium points of the planar system (3.20-1, 3.20-2) are Y” =
(0,w=,
c >T C,o
yL
Y
and
Here, y = (yO, r)= is the dynamical variable vector. Linear stability analysis reveals that the equilibrium point y” is a center with pure imaginary eigenvalues. Equilibrium point y” and y1 are saddle points with real eigenvalues of opposite sign. Phase plane in which we are interested consists of limit cycles around the center and homoclinic loop joining the saddle points. First integral of the planar system (3.20-1, 3.20-2) can be found by standard means: H = i yor2 - y&’ - r4.
(4. I)
The level curves are then defined by the relation ~(y~, r) = S = const. The parameter the value of
S takes
c4
on the center and zero on the homoclinic loop. Hence for a given value of wave celerity and for different values of y, there exist infinite number of limit cycles. Since the frequency of rotation around the center is smaller than the frequency about the y. axis (c/fia’
and a are anticipated as solutions to (1.3). Expressing y. in terms of r and S via first integral (4.1) and then substituting the latter in the planar system (3.20-1, 3.20-2), following equation can be found: 2
(
l/2
_r4+Cr2_~
i=*tr a*
(4.2)
>
4Y2
which yields immediately r1
I r
dp [@Z
_
p2)(p”
_
r;)]l”
=
*$
(*
+
cJ
(4.3)
where (4.4)
600
SUHUBi
and E. S.
G. ONAL
It is a simple matter to convert this integral to an elliptic integral of the first kind [lo]. Having done this, solution to planar system is obtained, after some algebra as y,=~_42r,snecne dn 8 2Y
r = rl dn 8,
e = $ (2 + c,)~,
sgn 8,
(4.5)
where rf - rf 2 rI
4’
and sn 8, cn 0, dn 0 are Jacobian elliptic functions. An expression for the period can also be found immediately by replacing the lower limit of integral (4.3) with r, and then integrating so that T(S) =
z
K(q)
1
where K(q) is a complete elliptic integral. Following limits would be useful in regard to understanding the behavior of solution (4.5) in the limit of homoclinic loop and the center: (4.7) Homoclinic loop is made up of two heteroclinic orbits. One of them emanates from the saddle point y’, follows its unstable manifold and ends on the saddle point y” following its stable manifold. Along with the lines mentioned in the periodic case, solution associated with this heteroclinic orbit is found to be: r=$sech8,y,=&(l-sgnBtanh8). We are also interested
e=$(z+C2).
in the values taken by the latter at z ---, w lim y. = f,
lim y. = 0,
,?---m
lim r = 0.
*--
(4.9)
z- *m
The other heteroclinic orbit completes rest of the homoclinic loop by emanating from y” and ending at y*. Similarly solution leading to heteroclinic orbit can be written as
Y yo = 1 +
r = 0,
(4.10)
e(cla*)r-d,
where d, is an arbitrary constant. If z ---, i-m, limits of this solution become lim y. = 0, L---m
lim y. = c. Z--rY
(4.11)
Original dynamical variable v can be expressed in terms of the variables yo, r and Q, by employing the normalizing transformation (3.5) together with the coefficients in the Appendix as follows: v=y,+ +&
(
2c-
3r
>
'Yy, ;;Jsinv+
r2 sin
2q +
[
-
9( y2 - 6a4)
19y* - 186a4 9a6
2a6
yZ+
135 y2 + 276a4 18a6 2
y,r2cos2q-2$yi-
r*+
1
2 r cos Q,
y2 + 96a4 36a6 r3cos3g,
(4.12)
601
Travelling waves and chaos in the KSS model
where v = vl, rp = az + Co and C,, is an arbitrary following expression can be found:
+2y2q4r: sn2 8 cn2 e sgn e a6 dn2 8
Substituting
(4.5) into (4.12)
q2r1 1
sn e cn e 3 y2q2rf sn 0 cn 8 dn 8 sin rp + sgn e a3 dn 0
+ _9(Y2 - 6a4) c2 + 135~’ + 276a4
r?dn28+2
18a6
8y2a6
constant.
1 1
q2rIsn8cn8 dn 8
I
q2r:sn8cn8
r,dntIcosg,
cos~+Lr~dn2Bsin2~ 3a3 *
2a6
19~’ - 186a4 cdn0 + rf dn 6 cos 29~ - sgn 8q2rI sn e cn e 9a6 2Y - y2 + 96a4 r: dn3 tI cos 391. 36a6
(4.13)
Here sgn 8 stems from (4.2). At the equilibrium point 3 of the planar system, values taken by the parameters are
s=c4
q=o
64 y4’
and
Using these values in (4.14) and assuming that (cl << 1, (4.13) simplifies to v=;
cl (
1 2+/zcoscp.
(4.14)
)
Thus, the equilibrium point y” of the planar system (3.20-1, 3.20-2) corresponds to a limit cycle in three-dimensional phase space (other coordinates being it and ti). Hence the predictions of the linear stability analysis are verified. While s<-
C4 64f
solution (4.13) and its first and second derivatives form a torus in three-dimensional phase space. This merely indicates that Neimark bifurcation is occurring. These are illustrated by
“3
J0
"1
,..’ /
.^
_.-.\ ,_.-
x0’
------+
(’
,,,:/
-b
._: ‘\.__._.
Hopf bifurcation
_.
-
”
Neimark
bifurcation
Fig. 1. Bifurcations from the equilibrium point and periodic orbits.
602
G. tiNAL
and E. S. SUHUBi
plotting the solution (4.13) and (4.14) in Fig. 1. Thus we have identified the first two bifurcations before the onset of turbulence. Hopf bifurcation followed by a Neimark bifurcation is known to be “Basic Sequence” [ll]. The latter is a prelude to the scenarios given in Refs [12-141. Right after the Neimark bifurcation infinite-period bifurcations occur. This can be discerned by considering the solution (4.13) together with expression (4.6). Langford [1.5] has demonstrated the existence of infinite-period bifurcations numerically in the context of codimension two bifurcations, The same author also expounded that Neimark bifurcations result from the interactions between simple and Hopf bifurcations. As it can be easily seen from (4.13) elliptic functions are modulated by trigonometric ones. In a sense, elliptic and trionometric functions play the role of large and small scale eddies, Amplitude of the smaller eddies are controlled by the parameter y which appear in front of the convective term of equation (1.1). This solution implies the existence of energy transfer from large eddies to smaller ones. Naturally, the coefficient 6 pertaining to negative viscosity term appears in solution (4.13). This can explain the propagation of travelling waves without dissipation on the viscous fluid. All of these agree well with interpretations given in Refs [l, 61. Substituting solution (4.8) into expression (4.12), solitary wave solution to the system (1.3) can be easily found. After having done this and taking limits we obtain lim u = 0. 2-m
(4.15)
When
first quantity in (4.15) vanishes. At this value of the wave celerity, solitary wave with many humps should appear (see Fig. 2). Similarly, solitary wave with single hump can also be obtained from (4.10) as it can be seen from Fig. 3. We have not pursued the analysis further to see whether these are also solitons or not.
Fig. 2. Solitary wave with humps plotted for the parameter values of y = 0.5, CY= 0.75, 6 = 0.25 and c = -0.459.
Travelling waves and chaos in the KSS model
V
603
-0 y
] -0.2 -
-0.4 -
-0.6-
,
I
I I I I,
-30
I,,
,
,
,
,
,
,
,
,
-20
,
,
I,,
,
,
,
,
,,
-10
,
,
,
,
,
,
,
+
,
,
0
,
,
,
,
,
10
,
,
,
,
,
,
20
Fig. 3. Solitary wave with a single hump plotted for the parameter values of y = 1, (Y= 0.25, 6 = 0.5 and c = -0.088.
5. MELNIKOV
ANALYSIS
What happens to the trajectories of the planar system (3.20-1,3.20-2) under the perturbation of higher order terms? In order to be able to provide an answer to this question, fou~h-order normal form analysis together with Melnikov analysis have to be carried out. By making use of the results obtained Section 3 and the coefficients given in the Appendix, normal form equations are found as
@=
a + WIYO, r12).
(5.1-3)
Resealing the variables and the celerity in (5.1-1 to 3) with a small parameter the following manner r= 6,
yo = C&
v,
2’ = Z/E
c = v&J,
E (0 < E << 1) in
(5.2)
yields * = f(w, p) + &g(W)*
(5.3)
Here, dot denotes the derivative taken with respect to z’ and w = (u,
~m~nents
ZI)T,
f=
vm,
f2w)=,
g =
tfim%
g*(w))‘.
of the vector field have the form fi(w) =$
(-J4U + 2yuu),
gr(w) = UrUV3+ b&,
f*(w) = - (/z?J - yu2 - 2yu2) :2
(5.4-l)
g,(w) = u2u4 + b2u2v2+ c2u4
(W-2)
604
G. fiNAL
and E. S. $UHUBi
where b = y(82y2 + 816a4)
a = _Y(17Y2 - 246a4) 1 4aX ’ b = 9Y(Y2 - 26a4) 2
2aX
1 c2= -
’
l&&8
y(27ly*+
y(4y2 - 36a4) ’
‘* =
ax
3786a4)
18aX
(5.4-3)
’
Melnikov function is proportional to the distance between stable and unstable manifolds inherent in homoclinic and heterclinic connections [9]. Following Salam [16], Melnikov integral for dissipative systems can be written as MC%) = [= f(q”(z - zoo))A g(q”(z - r,))exp[ -I--CC
tr D&q’(s))
ds] dz
(5-5)
Here, q” stands for the unperturbed trajectories, A denotes the wedge product (f/\g= fig, -figI) and tr D,o is the trace of Jacobian derivative along the unperturbed trajectory. Unperturbed periodic orbits of the phase plane corresponding to (5.3) can be easily obtained from (4.5) by replacing u with r, y. with u and c with p. Trace of the Jacobian derivative taken along the periodic orbits of the unperturbed system (5.3 with E = 0) is tr D,of(s) = sgn 0
Yq2rI sn 8 cn 8 dn8 ’
8 = 4 r,s. a
(5.6)
After having evaluated the exponential expression in (5.5), we can convert (5.5) into a line integral. Then by making use of the Green’s identity, we are led to the following integral u3 (2b, + 6,) II v + (4az + 2u,) ; u3} dv] du
(5.7)
where 1, =;+$
(
4y2s II2
+44y=u*--
u2 1 ’
P*-4y2u*__
12=g-$
(
4y2s 1’2 u* 1
’
rI and rz are obtained from (4.4) by replacing c with ~1,When the integral in braces is evaluated we obtain Z, + (B - 2A)Z, - 2ASZs3
(5.8)
where the following quantities and integrals used A=-
2a2 + al
2 Z, =
I
’
r’ (p2u2 - 4 y2u4 - 4~~s)“~ du,
B = 2b1 + b2
Z, =
(5.9-l)
” u2(p2u2 - 4Y2u4 - 4Y2S)“* du, IQ
*I *u2 - 4y*u4 - 4Y*S)l” Z, = r2 (P du. I r2 u2
(5.9-2)
After these integrals have been evaluated by making use of the integrals Melnikov function for this case is found to be 1270
gP2Y2SK+
270~’ - 355 r2SE 24
1
I
given in [lo],
(5.10)
605
Travelling waves and chaos in the KSS model
where E(q) and K(q) elliptic integrals of the first kind and second kind, respectively. For a fixed value of ~1, Melnikov function increases while the parameter S approaches to zero. Equilibrium point
(u, v)‘=
(*y, gT
of the unperturbed system (5.3 with E = 0) is a center. However, higher order monomials this equilibrium point switches to UC-
CL
ti
y + ’ 2304a6y3
2fi
under the perturbation
of
(179y* - 9186a4)p3 + O(s*)
ti P ~ (71~~ - 2976a4)p3 ’ = G + E 288a6 y3
+ O(E*).
(5.11)
Linear stability analysis near the equilibrium point reveals that it is indeed an unstable focus. Both analyses agree on the destruction of torus under the perturbation of higher order terms. Similar calculations for the homoclinic loop showed that Melnikov function has the form (5.12) Since Melnikov function does not have a simple zero in both cases, transversal intersection of stable and unstable manifolds does not occur as it has been anticipated. Thus it can be inferred that the homoclinic bifurcation should occur as a result of the perturbation of higher order terms [9].
6. TRANSITION
TO CHAOS
The consequence of the homoclinic bifurcation in the planar system and the existence of horseshoe mapping and chaotic dynamics heralded from it can best be viewed by constructing Silnikov’s example for the present problem. Following Guckenheimer [9], let us construct a cylinder of radius p. and height of {I in the neighborhood of the origin which is a saddle point (see Fig. 6.5.2 of [9]). L inearized system corresponding to (1.3) is obtained as
(6.1) where
By utilizing the similarity transformation (6.1): v=Sr,
(3.1), the following relations are defined in obtaining
r = s-iv,
r = (5, 9, 0’.
Solutions of (6.1) rules the flow near the saddle point:
[I [ E
77 (z) = I;
euZ[cos /3z E(O) - sin /3z n(O)] euZ[cos /3z E(O) + sin /3z n(O)]
e”’f (0)
. I
(6.2)
G. ONAL and E. s. SUHUBi
606
When trajectories pass by the origin, 15‘1 increases decreases. Let us define the surfaces Za and Z, by
and radial coordinate
p = (6” + $)“’
z:o = ((5 11, 5‘) I 5” + r2 = P:,>0-c 5 < 5‘1) 2, = ((5, rl, I;) I E2 + v2 = Pit, c = c-1‘0).
(6.3)
It has been assumed that Co and Ci are in the neighborhood U where the flow is linear [9]. The mapping 0 :X0+ Zi can be found by substituting the time of flight z = A-’ In(~l/~(0)) into (6.2). Thus, we obtain Vz = ((51/f)““](cos Here o = /In-’ ln(
w)5’ - (sin
whl, (1;JC;)“‘*[(sin ~15 + (~0so)rl, fd
By introducing the new variables can be obtained [9]:
(6.4)
f = p0 cos Y and rl = p. sin v, the
@(v, 5) = (~~(5;~/0”‘” cos(v + w), ~~(5;~/0”” sin(v + 0))
(6.5)
0 maps a vertical segment Y = constant of &, to a logarithmic spiral surrounding the 5; axis and lies in Xi. Stretching of the vertical segment, wrapping it around the spiral can easily be seen by calculating the determinant of Jacobian derivative [9]:
(6.6) Stretching ratio s = (plpo)/( 1;/5;1) can be found from the last quantity: s=-
Cl 5
1+0/A .
(>
(6.7)
Stability of the homoclinic orbit in the neighborhood of the saddle point is being determined by the stretch ratio and it is stable for s < 1 [6, 91. While k > 0 and o < 0, linear stability condition is found to be
14’ 14. When the situation is reversed,
linear stability is determined
14 ’ 14.
(6.8) by (6.9)
For the system (1.3) la/Al = l/2, thus it has a homoclinic orbit which is stable according to (6.9). For the right travelling waves near the v” and for the left travelling waves near the VI, homoclinic orbit can be observed. What follows from the analysis is that the vertical segment {o = constant, < E (0, &,)} on Z,, are mapped into a spiral of maximum radius po(l;l/l;)“” (C&J and the lengths of these segments are stretched by an amount which is unbounded as &,-,O [9]. Let us now consider a region A = {(p, w, 5) 1p = po, joI< q, 0 < I; -=cQ} near a point on the surface X0. The region A will be stretched under the mapping 0, yet its maximum radius will shrink when the condition (6.8) is fulfilled. According to Silnikov’s theorem when the condition is satisfied the region A fl (q!~0 0) (A) contains a horseshoe [9]. Horseshoe map has the properties of folding, stretching and mixing. Cantor like invariant set of the horseshoe map contains (a) a countable set of periodic orbits of all periods; (b) an uncountable set of nonperiodic motions; (c) a dense orbit; (d) the periodic orbits are all of saddle type and they are dense in the invariant set [9]. Thus at the parameter values close to the one which homoclinic orbit appears, chaos should prevail. Analytical results obtained hitherto will be verified by computing Poincare surface of section and LCEs. Initial conditions needed to construct surface of section obtained from the planar
607
Traveliing waves and chaos in the KSS model
system for different values of the parameter S (0 c S < c4/64y4). With these initial conditions dynamical system (1.3) has been integrated by employing fourth-order Runge-Kutta for the parameter values of c = -0.5, y = 1, cy = 2 and S = 0.5. Poincare surface of section can be seen from Fig. 4. Poincare map corresponding to the latter has a fixed point located at the center of section which has been obtained by taking
s=c4
64y4
This fixed point corresponds to a limit cycle in three dimensional phase space. By varying the parameter S closed curves su~ounding the fixed point appear. These are orbits having two incommensurate frequencies, which is so called quasiperiodicity, and they correspond to tori of different sizes in the phase space. One of the salient features of these orbits is their winding number. It is the ratio of the number of turns around the torus in each of the two directions [17]. Since the frequencies of nonlinear waves depend on the amplitudes, the winding number of orbits varies ~ntinuously for different initial conditions. Here, parameter f of the planar system is controlling the winding number. For the rational winding numbers orbits close on themselves after a finite number of iterations. For the irrational values of the winding number, the orbit covers the entire closed curve densely. RAM theorem (after Kolmogorov, Arnold and Moser) predicts that the orbits with sufficiently irrational winding number persist under small but finite ~~urbation. These orbits are also called RAM curves 191. Thus we infer that the dynamical system (1.3) is locally conservative. Hitherto all the properties of analytical solution (4.13) agress very well with the numerically obtained Poincare section. Formation of separatrix with seven hyperbolic fixed points and islands inside of it appearing in Fig. 4 can be understood by considering Moser twist map along with the Poincare-Birkhoff theorem [9]. When the frequency ratio of large scales (elliptic functions) and small scales (trigonome~c ~n~ions) is rational, perturbation of higher order terms switched on by the
E!
-0.4
y
-0.8
1
-1.0
I,,,,> -0.6
f,,,,,,,l,,,,,,,, -0.5
/ / -0.4
/
/ ,,
,,,,,(
-0.3
/i,,,
/,
-0.2
Vl
Fig. 4. Poincart? surface of section.
,
,,,,,,,,/,,,,
-0.1
,,,
II -0
608
G. ijNAL and E. S. SUHUBi
06
;;
-
04-
> E
3
E -I
0.2
-
-0
-
400
800
1200
1600
2000
2
Fig. 5. Positive LCE.
magnitude of celerity has an important effect on the periodic orbits. In other words resonance occurs and energy is exchanged between large and small scales. Direction of the energy transfer is determined by the relative phases and the parameters y and 6. Hyperbolic fixed points of the perturbed map and orbits connecting them from the separatrix (resonant curve). Islands inside of the separatrix are formed around the elliptic fixed points of the perturbed map. As a result 7: 1 resonant curve and islands appear on the surface of section. As a principle for every rational winding number there corresponds an island structure. Very close to the parameter value S = 0, celebrated chaotic trajectory showed up as having 0.6
-
I
Gi > 5
0.4
-
0.2
-
-O-
cu w Y
-0.2
-
-0.4
-
-06
1,,,,,,,,,,,,,,,,~,,,,,,,,,~,,,,,,,,,,,, -400
0
,,,,,,,,, 400
800 L
Fig. 6. Zero LCE.
1200
,,,,,,,,,, 1600
2cQo
Travelling waves and chaos in the KSS model
609
L
Fig. 7. Negative LCE.
fuzzy crown like (outer orbit in Fig. 4) structure. Homoclinic points resulting from the transversal intersections of the stable and unstable manifolds near the hyperbolic fixed points are prominent. It is worth noting that chaos sets in before the limit of homoclinic orbit. This is accounted for the parameter values taken here (i.e. c + u”/fi). Occurrence of homoclinic bi~rcation on the onset of the chaos is evident from Fig. 4. The chaotic trajectory observed here is a transient one and it escapes from the attractor after certain period of computer integration time. The closer the values of parameter S taken to the value of S for which separatrix appears the more time elapses before the trajectory escapes. Computation of LCEs are performed to check whether the chaotic orbit survives on the strange attractor or not. The LCEs are the average exponential rates of the convergence or divergence of the nei~borhood trajectories in the phase space fI8). If one of the LCEs is positive and their sum is negative, thus the chaotic trajectory lies on the strange attractor [19]. As seen from Figs 5-7, the LCEs for the chaotic trajectory are A1= 0.0415 bit/s. Since the sum of these exponents is negative (--0.00002) and one of them is positive, chaotic trajectory lies on the strange attractor. Hausdorff dimension of the latter is calculated by Kaplan-Yorke conjecture 1191 and is found to be approximately 2.9976. REFERENCES B. NICOLAENKO, &cl. Phys. {Proc. Suppl..) 2,453 (1987). G. I. SIVASHINSKY, Physicu 17D,243 (1980). M. C. DEPASSIER and E. A. SPIEGEL, A.rrron. I; &t(3), 4% (1981). D. K. CAMPBELL, Los AIqos Sci. Spec. iss. 218 (1987). G. UNAL and E. S. SUHUBI, Znt. J. Engng Sci. 30,579 (1992). J. GUCKENHEIMER, SIAM J. Math. Anal. xS(l), 1 (1984). V. M. STARZHINSKI, Applied Methods in the Theory of Nonlinear Oscillations. Mir, Moscow (1980). G. UNAL, A group theoretical approach to turbulence (Ph.D. thesis in Turkish). Istanbul Technical University, Istanbul, Turkey (1991). J. GUC~N~I~R and P. HOLMES, Nonlinear Osc~fatio~, ~y~mical Systems and Bi~rcation of Vector Piehis. Springer, New York (1983). I. S. GRADSHTEYN and I. M. RYZHIK, Tables of integrals Series und Products. Academic Press, New York (1971). R. H. ABRAHAM and J. E. MARSDEN, Foundations of Mechanics. Benjamin-Cummings, Reading MA (1978).
610
G. ANAL and E. S. SUHUBI
[12] [13] [14] [15]
L. D. LANDAU and E. M. LIFSHITZ, Fluid Mechanics. Pergamon Press, Oxford (1959). D. RUELLE and F. TAKENS, Commun. M&z. Phys. 20, 167 (1971). S. NEWHOUSE, D. RUELLE and F. TAKENS, Commun. M&h. Phys. 64, 35 (1978). W8:j.LANGFORD, in Chaos, Fractals and Dynamics (Edited by P. FISCHER). Marcel Dekker, New York
[16] [17] [18] [19]
F. M. SALAM, SZAM J. Appl. Math. 47(2), 233 (1987). J. M. GREENE, Physica 18D, 427 (1986). G. BENETTIN, L. GALGANI, A. GIORGILLI and J. M. STRELCYN, Mecca&a 15,9 (1980). A. WOLF, J. B. SWIFT, H. L. SWINNEY and J. A. VASTANO, Physica 16D, 285 (1985). (Received 31 August 1991; accepted 11 September
1991)
APPENDIX The coefficients of resonant monomials contributing which are not vanishing will be reproduced here:
K
,oo,= Y’ F,’
y2 - 96a4
h4
K,,,., = -51, 52oooo=
K-,-,-I, = -K,,,.,,
%x1-1 =
82yZ+ 816a4
216us
Q”.,.,,, = -
to the normal form equations are given below. Only the ones
18a
4yz - 36a4 aR yp
Q-,0,x,-, = f&u,,,
yj
271 y2 - 3786a4 108~s ”
K.,w-I = -K,nn, Q,,,,
= -
Q,“.,.,,
a,,=--lli
a loo=&i,
’
ab,, = “lx&,
a10.1= (r,,/z
~11,= -I,,
al-l-1 = %0/3, =11--l
&n-1=
=
o-,00= a-,-,-*
-@I,.,,
13y* + 96a4 Q
a3
F
=
w-,-,0
8*-1-m= B ,--1-1-1 =-q$Y, B- ,110 --B ,-,-*o~
~,I., = ~ICKI
ooo,
-(Y,fJo,
B
L-111 = B,-,-,-,r
a- ,,, =
Y2 -___
r%&,,
5 y2 - 96a4
pr_,,=-$>
=
’
= B”,,&,j
IL-11
T
-ff1-I-,
&14
u6
= BIX,“?
27a6
are:
cu,.,.,= acd2
ff_ 10, = -o,o-,?
-o,m,
= Q,“,,.,
52
And the coefficients of monomials appearing in the normalizing transformation %o,=,s’,y.
17y2 - 24Sa4 l&, y,
= -8,,,o>
B-,00,=-/%c#J,, /L,,,-,
= -7yy&:ba4’
8-I-1-10
= 8,,1”P
Pm,-,-,-,
= /$I,,.
G42)