Bifurcation analysis of a density oscillator using two-dimensional hydrodynamic simulation
aa r X i v : . [ n li n . PS ] M a y Bifurcation analysis of a density oscillator using two-dimensional hydrodynamicsimulation
Nana Takeda, Naoko Kurata, Hiroaki Ito, and Hiroyuki Kitahata ∗ Department of Physics, Chiba University, Chiba 263-8522, Japan (Dated: May 11, 2020)A density oscillator exhibits limit-cycle oscillations driven by the density difference of the two flu-ids. We performed two-dimensional hydrodynamic simulations with a simple model and reproducedthe oscillatory flow observed in experiments. As the density difference is increased as a bifurcationparameter, a damped oscillation changes to a limit-cycle oscillation through a supercritical Hopfbifurcation. We estimated the critical density difference at the bifurcation point and confirmed thatthe period of the oscillation remains finite even around the bifurcation point.
I. INTRODUCTION
Limit-cycle oscillations often appear in various systemswith energy gain and dissipation. Heartbeat, circadianrhythm, and firefly flashing are the famous examplesin biological systems [1–3]. Cyclic phenomena are alsofound in fluid systems such as geysers, thermohaline cir-culations, and the solar cycle [4, 5]. These limit-cycleoscillations in nature are often complex because manyfactors are cooperated, and thus it is difficult to under-stand the underlying mechanism. Therefore, it should bea good approach to first understand the essential mech-anism of ideal systems and then address more complexsystems.A density oscillator is an example of limit-cycle oscilla-tors in fluid systems, first reported by Martin in 1970 [6].The system consists of two different-density fluids putin two fixed containers. The inner and outer containersare for the higher- and lower-density fluids, respectively.The small hole at the bottom wall of the inner containerconnects the two fluids. In this system, the upstream ofthe lower-density fluid and the downstream of the higher-density fluid alternately occur through the hole in appro-priate conditions. Owing to the simplicity of the setup, adensity oscillator has been investigated mainly in exper-iments [7–23]. In previous studies, it has been reportedthat a density oscillator shows typical characteristics ofa limit-cycle oscillation such as orbital stability and syn-chronization among several oscillators [15–19]. In addi-tion, some studies proposed theoretical models by adopt-ing the Navier-Stokes equation for each of the upstreamand the downstream in the hole. Steinbock and coau-thors derived the theoretical description of the criticalwater level in the inner containers for the reversal of thedownstream in a two-dimensional system [20]. Kano andKinoshita focused on the intrusion of the different-densityfluid in the hole and proposed a model which explainedthe processes of the upstream, the downstream, and theswitching between them in a unified manner [21, 22].In spite of intensive studies on the dynamics of theflow, bifurcation structure by the change in parameters ∗ [email protected] such as the hole size, the density difference between thetwo fluids, and the viscosity of the fluids has not been in-vestigated in detail. To understand the dynamics of non-linear systems, it is important to interpret the dynam-ics qualitatively from the viewpoint of bifurcation, wherethe qualitative behavior of the system changes with a bi-furcation parameter. Our recent experiment suggestedthat the density oscillator shows a supercritical Hopf bi-furcation between the resting and oscillatory states de-pending on the density difference between the two fluids[23]. However, we could not definitely identify the bifur-cation class only from the experimental results becausethe measurement of a small-amplitude oscillation arounda bifurcation point suffered from relatively large error.In order to address the identification of the bifurca-tion structure, numerical simulation is useful since onecan accurately quantify the dynamics in various condi-tions. In general, hydrodynamic simulation for a densityoscillator has not been performed because it is difficultto treat the free surface which changes with time exceptfor few studies: Okamura and Yoshikawa carried out athree-dimensional hydrodynamic simulation for a densityoscillator with the free surface by adopting the volume offluid method [24]. In their study, they considered essen-tial factors in the system, e.g., the gradients of pressure,viscosity, and gravity, and proposed that the oscillationfollows the Rayleigh equation. The simulation param-eters were fixed to perform costly calculations, and thebifurcation structure remains as an open question. Us-ing such ordinary differential equations as the Rayleighequation, which were obtained by the reduction of the hy-drodynamic equation, mathematical analyses were per-formed and the type of bifurcation was discussed [25, 26].However, the bifurcation analysis directly based on thehydrodynamic equation has not been performed so far,which should be important to understand the actual dy-namics of the density oscillator and to discuss the validityof such reductions.In the present study, we carry out a two-dimensionalhydrodynamic simulation for a density oscillator. In thesimulation, we set the calculation area inside the fluidand associate the change in the water levels with thepressure at the calculationarea boundaries. Using thissimple model, we obtain the time series of the densityprofile and the water level. Then we investigate the de-tailed bifurcation structure depending on the density dif-ference between the two fluids. II. MODEL
We carried out a two-dimensional hydrodynamic sim-ulation for a density oscillator. For the incompressibleviscous fluid, we use the Navier-Stokes equation and theequation of continuity as governing equations, ρ (cid:20) ∂ v ∂t + ( v · ∇ ) v (cid:21) = −∇ p + µ ∇ v + ρ g , (1) ∇ · v = 0 , (2)where ρ ( r , t ) is the fluid density, v ( r , t ) =( v x ( r , t ) , v y ( r , t )) is the fluid velocity, p ( r , t ) is thepressure, µ is the fluid viscosity, and g is the accelerationof gravity. Here, r = ( x, y ) is a positional vector.For a miscible two-phase fluid, we define a normalizedconcentration c ( r , t ) (0 ≤ c ≤ c = 0 and c = 1 correspond to the lower- and higher-density fluids,respectively. The density is described using c as ρ ( r , t ) = ρ high + ( ρ low − ρ high )(1 − c ( r , t )) , (3)where ρ low and ρ high are the densities of the lower- andhigher-density fluids. The normalized concentration c satisfies an advection-diffusion equation as ∂c∂t + ∇ · ( c v ) = D ∇ c, (4)where D is a diffusion coefficient. We set D to be small,so that the fluids hardly mix with each other [23].Figure 1 shows the schematic drawing of a density os-cillator, where the calculation area is set inside the fluid.We set the origin of the Cartesian coordinates at thecenter of the hole and define the area for the hole as − a ≤ x ≤ a, − b ≤ y ≤ b . The calculation area is de-fined as − W ≤ x ≤ W, − b − H lower ≤ y ≤ b + H upper .Inside it, the areas − W ≤ x < − a, − b < y < b and a < x ≤ W, − b < y < b correspond to the bottom wallof the inner container. We assume that the pressures atthe upper and lower boundaries of the calculation area, p upper ( t ) and p lower ( t ), are determined by the water levelsin the inner and outer containers, y in and y out , as p upper ( t ) = ρ high g ( y in ( t ) − b − H upper ) , (5a) p lower ( t ) = ρ low g ( y out ( t ) + b + H lower ) , (5b)where g = | g | . Here, we set the atmospheric pressure tobe zero. The change in the water level is determined bythe amount of the fluid flowing through the hole per unittime Q ( t ), which is calculated as Q ( t ) = Z a − a v y ( x, b, t ) dx. (6) p upper ( t ) d in d out / 2 Q ( t ) xy d out / 2 p lower ( t ) y in ( t ) y out ( t ) (0,0) H upper H lower b FIG. 1. Schematic drawing of a density oscillator when anupstream occurs. The calculation area is denoted with thebroken rectangle, and the hatched area corresponds to thebottom wall of the inner container. p upper ( t ) and p lower ( t )are the pressures at the upper and lower boundaries of thecalculation area, respectively. Q ( t ) is the amount of the fluidflowing through the hole per unit time. y in ( t ) and y out ( t ) arethe water levels in the inner and outer containers, respectively. Due to the incompressibility, the change in the waterlevels can be directly obtained from Q as dy in dt = Qd in , (7a) dy out dt = − Qd out , (7b)where d in and d out are the widths of the inner and outercontainers, respectively, as shown in Fig. 1.In this way, we associate the change in the water levelswith the pressure at the calculation-area boundaries, andobtain a simple model, where the free surface need notbe directly dealt with. A non-slip boundary conditionfor the velocity, v = , and the Neumann boundary con-dition for the density, ∇ ⊥ ρ = 0, are set for the surfacesof the bottom wall, where ∇ ⊥ represents the derivativein the normal direction. The pressure at the surfacesof the bottom wall is determined so that it satisfies theNavier-Stokes equation in Eq. (1). The pressure at theupper and lower boundaries of the calculation area fol-lows Eqs. (5a) and (5b). The velocity, the pressure, andthe density at the other calculation-area boundaries alsofollow the Neumann condition. At the initial state, thetwo fluids are stationary ( v = ) and not mixed, wherethe concentration is set as c = 0 ( y < b ) and c = 1( y ≥ b ). The initial pressure of the lower- and higher-density fluids, p low , ( x, y ) and p high , ( x, y ), are given as p low , ( x, y ) = ρ low g ( y out (0) − y ) , (8a) p high , ( x, y ) = ρ high g ( y in (0) − y ) , (8b) t = 43.5 bottom wall (a)(b) (c) c = 1 c = 044.5 45.5 46.5 t = 60 61 62 63 t = 42 46 50 54 58 62 66 70 74 78 (d) (e) FIG. 2. Snapshots of density profile (a–c) and velocity field (d,e) for ∆ ρ = 0 .
2. (a) Typical time series of the oscillatory flow.(b,d) Detailed time series of downstream. (c,e) Detailed time series of upstream. The hatched area corresponds to the bottomwall of the inner container. where y out (0) and y in (0) are the initial water levels inthe outer and inner containers, respectively. y in (0) isdetermined by the gravitational equilibrium as y in (0) = b + ( ρ low /ρ high )( y out (0) − b ). In the simulation, we repre-sent the surface height of the fluid in the outer containerdefined as h ( t ) = y out ( t ) − y out (0) , (9)where the initial water level in the outer container y out (0)is not relevant to the behavior of the fluid because onlythe time derivatives of the water levels are included inthe governing equations in Eqs. (7a) and (7b).To solve the Navier-Stokes equation (Eq. (1)) and theequation of continuity (Eq. (2)), we used the Marker andCell method and calculated the velocity and the pressureon a staggered grid [27, 28]. We used an explicit methodfor the advection-diffusion equation (Eq. (4)) to calcu-late the time evolution of c . Here, we set the time step dt = 0 . dx = dy = 0 . W = 0 . H upper = H lower = 0 .
55. The simulation parameters were set as fol-lows: µ = 1 / , D = 0 . , g = 10 , a = 0 . , b = 0 . d in = d out = 4 .
8, and y out (0) = 10 .
05. To investigate thebifurcation structure, we changed the density difference∆ ρ (= ρ high − ρ low ) as a bifurcation parameter, where wefixed ρ low = 1 and varied ρ high . III. SIMULATION RESULTS
The snapshots of the density profile and the velocityfield for the density difference ∆ ρ = 0 . t = 43 .
5) and passes through the lower boundary of the calculationarea ( t = 46 . t = 60). Then theupstream grows and passes through the upper boundaryof the calculation area ( t = 63). The upstream becomesweaker and changes to the downstream again. In thisway, the limit-cycle oscillation consisting of the upstreamand the downstream occurs with a period of ∼ h . Figure 3(a) shows a relaxationoscillation in which h asymptotically approaches the twoheights alternately. As the density difference decreases,the waveform around the peaks changes as shown inFig. 3(b). The system does not exhibit a relaxation oscil-lation but exhibits a harmonic-like oscillation as shownin Fig. 3(c). In addition to the waveform change, theamplitude of the oscillation decreases, and the period in-creases. For the smaller density difference, the oscilla-tion gets damped as shown in Fig. 3(d). Figure 4 showsthe trajectories in a phase space of h and the flow ratethrough the hole Q , which corresponds to − d out dh/dt .For ∆ ρ = 0 . , . , and 0 .
05, the trajectory convergesto a closed orbit, suggesting that the system exhibitsa limit-cycle oscillation as shown in Figs. 4(a–c). For∆ ρ = 0 . ρ .To obtain the amplitude and the period of the os-cillation of h for each ∆ ρ , we detected the instanceswith Q = 0 and set the i -th instance as t i ( i =1 , , · · · ). t was set to be 0 since Q = 0 in the t h (a)(b)(c)(d) ttt h h h FIG. 3. Time series of h for (a) ∆ ρ = 0 .
2, (b) ∆ ρ = 0 .
1, (c)∆ ρ = 0 .
05, and (d) ∆ ρ = 0 . -20-12.5 (a) (b) h Q (c) (d) Q Q Q h h h FIG. 4. Trajectories in a phase space of h and Q for (a)∆ ρ = 0 . ≤ t ≤ ρ = 0 . ≤ t ≤ ρ = 0 .
05 (0 ≤ t ≤ ρ = 0 .
025 (0 ≤ t ≤ initial condition. The n -th period was defined as T n = t n − t n − and the n -th amplitude was definedas A n = [ | h ( t n − ) − h ( t n − ) | + | h ( t n − ) − h ( t n ) | ] / n = 1 , , · · · ). Then we defined the amplitude A andthe period T as A = ( A + A ) / T = ( T + T ) / A n and T n for n = 1 , , and 3 were not used toavoid the initial value dependence. Figure 5 shows theamplitude A and the period T of the oscillation of h for0 . ≤ ∆ ρ ≤ .
2. It should be noted that the systemfor ∆ ρ = 0 .
025 exhibits a damped oscillation as shownin Figs. 3(d) and 4(d), and thus the region of ∆ ρ plotted (a) (b)
080 0.2 A Δρ Δρ T FIG. 5. Bifurcation diagram of the density oscillator. (a) Am-plitude A and (b) period T of the oscillation of h dependingon ∆ ρ . Each inset shows the expanded area from ∆ ρ = 0 . ρ = 0 . in Fig. 5 reflects the behaviors for ∆ ρ in the both sidesof the bifurcation point. With a decrease in the den-sity difference, the amplitude steeply decreases to zerofor ∆ ρ . .
03, which suggests the existence of a bifurca-tion point ∆ ρ c as shown in Fig. 5(a). The period is largeespecially around the bifurcation point, while it clearlyremains finite and does not diverge as shown in Fig. 5(b).Since the resting state changes to the oscillatory statewith a finite period according to the increase in the den-sity difference as a bifurcation parameter, the bifurcationis classified into the supercritical Hopf bifurcation. IV. DISCUSSION
The simulation results show that a density oscilla-tor exhibits the supercritical Hopf bifurcation with thechange in the density difference as a bifurcation param-eter. In Fig. 5(a), we evaluated the amplitude from thefinite-time behaviors. In general, the convergence to afixed point or a limit cycle takes long time around thebifurcation point. Thus, the amplitude A calculatedonly from A and A could have a larger value thanthe amplitude of a limit cycle for ∆ ρ ≥ ∆ ρ c and 0 for∆ ρ ≤ ∆ ρ c . Here, we investigate the behavior of theamplitude and the period around the bifurcation pointin detail, and we estimate the critical density differenceat the bifurcation point ∆ ρ c . Figures 6(a–d) show thechanges in the surface height in the outer container h for∆ ρ = 0 . , . , , and 0 .
03 around the bifurcationpoint. Due to the slow convergence, we could not deter-mine whether the system shows a limit-cycle oscillationor a damped oscillation only from the apparent wave-forms within a finite duration. To judge the types of theoscillation, we check the time evolution of the n -th am-plitude A n for each ∆ ρ as shown in Fig. 6(e). A n shouldconverge to a nonzero value for a limit-cycle oscillationand to 0 for a damped oscillation as a long-time behav-ior. We also evaluate a damping rate | A n − A n +1 | /A n asshown in Fig. 6(f). It should converge to 0 for a limit-cycle oscillation and to a nonzero value for a dampedoscillation under the following assumption; if we assume tttt h (a)(b)(c)(d)
31 52 4 n | A n - A n + | / A n (f)(e) n A n
200 400200 400200 400200 400 0 2 4 6 8 1002460 h h h FIG. 6. Time series of h for (a) ∆ ρ = 0 .
03, (b) ∆ ρ = 0 . ρ = 0 . ρ = 0 . A n plotted against n for ∆ ρ = 0 .
03 ( • ; black), 0.029 ( ◦ ; red),0.028 ( (cid:4) ; green), 0.027 ( (cid:3) ; cyan), 0.026 ( N ; blue), and 0.025( △ ; magenta). (f) Damping rate | A n − A n +1 | /A n plottedagainst n for each ∆ ρ . The symbols correspond to those in(e). the amplitude of a linear damped oscillation representedby A ( t ) = A (0) e − γt , a damping rate is calculated as | A n − A n +1 | /A n = 1 − e − γT , where γ is a damping co-efficient. A damping rate converges to a nonzero valuebetween 0 and 1 because the period T has a constantvalue. The results shown in Figs. 6(e,f) suggest that thesystem exhibits a limit-cycle oscillation for ∆ ρ ≥ . ρ ≤ . ρ − ∆ ρ c ) / [29]. Figure 7(a) shows the squaredamplitude A plotted against ∆ ρ and the linear fit-ting by a least squares method with the three values∆ ρ = 0 . , . , and 0 . ρ for alimit-cycle oscillation close to the bifurcation point. For∆ ρ ≥ . A gradually deviate from thelinear fitting. The region of the linear fitting close to thebifurcation point means that the amplitude A increasesfrom 0 with the scaling of (∆ ρ − ∆ ρ c ) / . We estimate∆ ρ c = 0 . A = 0. Figure 7(b) shows the period T plottedagainst ∆ ρ − ∆ ρ c , where ∆ ρ c is estimated in Fig. 7(a).The slope of the period changes at around the bifurca-tion point ∆ ρ = ∆ ρ c , and the period has finite values for∆ ρ ≤ ∆ ρ c as well as ∆ ρ ≥ ∆ ρ c , which corresponds tothe behavior of the supercritical Hopf bifurcation.The simulation result agrees well with our previous ex-perimental result, where pure water and a sodium chlo-ride aqueous solution were adopted as the lower andhigher density solutions, respectively [23]. Actually, weobtained the bifurcation diagram in Fig. 5 qualitatively (a) (b) A Δρ Δρ - Δρ c T Δρ c =0.02880.5 -0.006 FIG. 7. Critical density difference at the bifurcation point∆ ρ c and the period around ∆ ρ c . (a) Squared amplitude A depending on ∆ ρ and linear fitting with the three values ∆ ρ =0 . , . , and 0 . T depending on ∆ ρ − ∆ ρ c . consistent with that obtained from the experiment, wherethe amplitude increased, and the period decreased, ac-cording to the increase in the density difference abovethe bifurcation point. The amplitude increased fromthe bifurcation point approximately with the scaling of(∆ ρ − ∆ ρ c ) / both in the simulation and in the experi-ment. Leaving the region close to the bifurcation point,these systems exhibited slightly different behaviors. InFig. 7(a) the values of the squared amplitude A grad-ually deviated downward from the linear fitting with anincrease in the density difference, whereas they deviatedupward in the experimental result. Despite the differ-ence in the behaviors far from the bifurcation point, itis notable that the simulation and experimental resultsexhibited the same bifurcation structure in the regionaround the bifurcation point.Here we discuss the validity of the two-dimensionalsimulation compared with the three-dimensional experi-mental system. It is expected that the fluid amount flow-ing through the hole and the resistance inside the holein the two-dimensional system are quantitatively differ-ent from those in the three-dimensional system. In spiteof these quantitative differences, bifurcation structures inthe both systems can be discussed in the same frameworkas follows. The existence of a critical density differenceis interpreted from the competition between the diffusionand advection of the solute. To evaluate the critical den-sity difference ∆ ρ c in the experiment, we considered theRayleigh number Ra = ( gL ∆ ρ/ρ ) / ( νD ), where L is acharacteristic length, and ν = µ/ρ is the kinematic vis-cosity of the fluid. We adopted Ra ∼ , which is typ-ical for Rayleigh-B´enard instability, and the estimatedvalue of ∆ ρ c agreed well to the experimental result [23].The critical density difference is estimated also by thesimulation parameters: g = 10, L ∼ . µ = 1 / D = 0 . L . The estimated value ∆ ρ c ∼ .
03 approx-imately agrees with the order of the estimation by thesimulation result in Fig. 7(a). The above discussion fromthe viewpoint of the transport phenomenon suggests thatthe bifurcations between the resting and oscillatory statesoriginate from the same mechanism, regardless of the dif-ference in dimensionality.In previous studies, the ordinary differential equa-tions reduced from a hydrodynamic equation were oftenadopted as simple mathematical models for the densityoscillators, and the detailed bifurcation structures for themathematical model have been investigated. For exam-ple, Aoki proposed a model with a nonlinear frictionalterm and claimed that the system has three fixed points.He investigated the stability of each point depending onthe bifurcation parameter, which accelerates the fluid ve-locity around the orifice, and discussed the bifurcationstructure [25]. Kenfack et al. also derived the piecewise-smooth ordinary differential equations and claimed thatthe system exhibits the nonconventional Hopf bifurcationby varying the resistance coefficient, determined by theviscosity and the geometry of the orifice, as a bifurca-tion parameter [26, 30]. It is difficult to directly discussthe correspondence between their models and our simu-lation result because their bifurcation parameters are es-sentially different from ours, i.e., the density difference.We guess that their results on the detailed bifurcationstructure may come from the assumption or approxima-tion in the reduction processes.In an actual density oscillator, the higher- and lower-density fluids mix, and the amplitude and the equilib-rium height changes with time. In our model, we as-sume that the density in the outside of the calculationarea is constant. Strictly, it should not be constant be-cause the different-density fluid goes through the calcula-tion boundaries. Considering that only the derivative of y out ( t ) is important and the absolute value of y out ( t ) doesnot matter owing to Eqs. (7a) and (7b), the size outsideof the calculation area can be set arbitrarily. Therefore,our model can be regarded as an ideal system in which the mixing effect between the higher- and lower-densityfluids is negligible independently of the choice of y out (0). V. CONCLUSION
We carried out a two-dimensional hydrodynamic sim-ulation for a density oscillator using a simple model asso-ciating the change in the water levels with the pressureat the calculation-area boundaries. Using this simula-tion, we clarified that a density oscillator shows super-critical Hopf bifurcation between a damped oscillationand a limit-cycle oscillation depending on the density dif-ference ∆ ρ . We confirmed that the amplitude increasesfrom 0 with the scaling of (∆ ρ − ∆ ρ c ) / , and the periodhas a finite value at the bifurcation point estimated fromthe scaling of the amplitude. The simulation of a densityoscillator can be useful for the investigation of other phe-nomenon in density oscillators such as synchronizationsand phase responses. The present study will contributeto further understanding of the nonlinear phenomenon inoscillators in fluidic systems, as well as a density oscilla-tor. ACKNOWLEDGMENTS
This work was supported by JSPS KAKENHI GrantsNo. JP19H00749, No. JP19K14675, and No.JP16H03949. It was also supported by Sumitomo Foun-dation (No. 181161), the Japan-Poland Research Co-operative Program Spatiotemporal patterns of elementsdriven by self-generated, geometrically constrained flows,and the Cooperative Research Program of NetworkJoint Research Center forMaterials and Devices (No.20194006). [1] S. H. Strogatz,
Sync: The Emerging Science of Sponta-neous Order (Hyperion, New York, 2003).[2] A. T. Winfree,
The Geometry of Biological Time (Springer, New York, 1980)[3] J. D. Murray,
Mathematical Biology I: An Introduction (Springer, Berlin, 1989).[4] P. Ball,
Flow, Nature’s patterns: A tapestry in three parts (Oxford University Press, New York, 2009).[5] C. Wunsch, Science , 1179 (2002).[6] S. Martin, Geophys. Fluid Dyn. , 143 (1970).[7] P.-H. Alfredsson, and T. Lagerstedt, Phys. Fluids , 10(1981).[8] K. Yoshikawa, S. Maeda, and H. Kawakami, Ferro-electrics , 281 (1988).[9] K. Yoshikawa, S. Nakata, M. Yamanaka, and T. Waki, J.Chem. Educ. , 205 (1989).[10] S. Upadhyay, A. K. Das, V. Agarwala, and R. C. Srivas-tava, Langmuir , 2567 (1992).[11] R. Cervellati and R. Sold`a, Am. J. Phys. , 543 (2001).[12] M. Ueno, F. Uehara, Y. Narahara, and Y. Watanabe,Jpn. J. Appl. Phys. Pattern formations & oscillatory phenomena ,edited by S. Kinoshita (Elsevier, Waltham, MA, 2013),pp. 119-163.[14] S. Maki, C. Udagawa, S. Morimoto, and Y. Tanimoto,Microgravity Sci. Technol. , 125 (2014).[15] K. Yoshikawa, N. Oyama, M. Shoji, and S. Nakata, Am.J. Phys. , 137 (1991).[16] S. Nakata, T. Miyata, N. Ojima, and K. Yoshikawa,Physica D , 313 (1998).[17] K. Miyakawa and K. Yamada, Physica D , 217(2001).[18] H. Gonz´alez, H. Arce, and M. R. Guevara, Phys. Rev. E , 036217 (2008).[19] M. Horie, T. Sakurai, and H. Kitahata, Phys. Rev. E ,012212 (2016).[20] O. Steinbock, A. Lange, and I. Rehberg, Phys. Rev. Lett. , 798 (1998).[21] T. Kano and S. Kinoshita, Phys. Rev. E , 046208(2007).[22] T. Kano and S. Kinoshita, Phys. Rev. E , 046217(2009). [23] H. Ito, T. Itasaka, N. Takeda, and H. Kitahata, EPL, , 18001 (2020).[24] M. Okamura and K. Yoshikawa, Phys. Rev. E , 2445(2000).[25] K. Aoki, Physica D , 187 (2000).[26] W. F. Kenfack, M. S. Siewe, and T. C. Kofane, ChaosSolitons Fract. , 321 (2018). [27] F. H. Harlow and J. E. Welch, Phys. Fluids , 2182(1965).[28] S. McKee, M. F. Tom´e, V. G. Ferreira, J. A. Cuminato,A. Castelo, F. S. Sousa, and N. Mangiavacchi, Comput.Fluids , 907 (2008).[29] S. H. Strogatz, Nonlinear dynamics and chaos: With ap-plications to physics, biology, chemistry, and engineering (Addison-Wesley, Reading, 1994).[30] W. F. Kenfack, M. S. Siewe, and T. C. Kofane, ChaosSolitions Fract.82