- Email: [email protected]
Chaos in the lid-driven square cavity
Accepted Manuscript Chaos in the lid-driven square cavity Salvador Garcia
PII: DOI: Reference:
S0378-4754(17)30186-6 http://dx.doi.org/10.1016/j.matcom.2017.04.010 MATCOM 4457
To appear in:
Mathematics and Computers in Simulation
Received date : 30 June 2015 Revised date : 10 April 2017 Accepted date : 26 April 2017 Please cite this article as: S. Garcia, Chaos in the lid-driven square cavity, Math. Comput. Simulation (2017), http://dx.doi.org/10.1016/j.matcom.2017.04.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Chaos in the lid-driven square cavity Salvador Garcia Instituto de Ciencias F´ısicas y Matem´ aticas Universidad Austral de Chile Casilla 567, Valdivia Chile
Abstract As the Reynolds number Re −→ ∞ the primary eddy switches from clockwiseto counterclockwise-rotating, a switching repeatedly reversing as the time t −→ ∞. The primary eddy is the eddy leading the dynamic. Marking the genesis of chaos, the switching first is observed at Re = 307, 000 at the upper side of the total kinetic energy range. But then, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range and unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. At extreme Reynolds numbers 1 Re = O( ) = O(1015 ), where eps is the machine epsilon in double precieps sion, as the time t −→ ∞, an interesting dynamic still endures where the switching repeatedly reverses or where it stays put, the primary eddy rotating solely counterclockwise. Keywords: Navier-Stokes equations, lid-driven square cavity flow, chaos 1. Introduction As the Reynolds number Re −→ ∞ the primary eddy switches from clockwise- to counterclockwise-rotating, a switching repeatedly reversing as the time t −→ ∞. The primary eddy is the eddy leading the dynamic. Marking the genesis of chaos, the switching first is observed at Re = 307, 000 at Email address: [email protected] (Salvador Garcia) Preprint submitted to Mathematics and Computers in Simulation
April 10, 2017
the upper side of the total kinetic energy range. But then, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range and unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. At extreme Reynolds numbers 1 Re = O( ) = O(1015 ), where eps is the machine epsilon in double precieps sion, as the time t −→ ∞, an interesting dynamic still endures where the switching repeatedly reverses or where it stays put, the primary eddy rotating solely counterclockwise. Chaos appears as a systematic deterioration of a stationary configuration, while always obeying a definite set of principles [26, 27] and while in the end reversing somehow such a stationary configuration. The domain is the unit square cavity. And the viscous incompressible flow is governed by the two-dimensional time-dependent incompressible NavierStokes equations (NSE) [45] and driven by the upper wall which moves from left to right, see Figure 1. The nondimensionalized NSE in primitive variables with Dirichlet boundary conditions over the domain Ω = ]0, 1[×]0, 1[ read ∂u 1 − ∆u + (u · ∇)u + ∇p = f in Ω, t > 0, ∂t Re ∇ · u = 0 in Ω, (1) u = ϕ on Γ = ∂Ω, u(x, 0) = u0 (x) in Ω, where u is the velocity; p, the pressure; Re, the Reynolds number; f, the external force. Here, f = 0 and the boundary conditions read ( u(x, ·) = (1, 0) if x ∈ upper wall, (2) u(x, ·) = 0 if x ∈ left, bottom, or right wall.
At first, an unexpected balance of viscous and pressure forces makes the fluid to turn into the unit square cavity. The properties of these forces depending upon the Reynolds number, a hierarchy of eddies develops. First, the large clockwise-rotating primary eddy (1) arises. Its location occurs toward the geometric center of the unit square cavity. And then several small eddies arise: the counterclockwise-rotating secondary eddies (2), the clockwise-rotating tertiary eddies (3), the counterclockwise-rotating quaternary eddies (4), the clockwise-rotating quinary eddies (5). Their locations 2
occur at the three relevant corners of the square cavity: bottom right (BR), bottom left (BL), and top left (TL), appearing hierarchically at the inclined ellipses in Figure 1. This is the stationary configuration. A combination of known methods is used to discretize and solve the NSE (see [24] for further details). But yet, to summarize, the linear Linθ∗ scheme [39] is the frame of work for the temporal discretization of the NSE, a scheme possessing the following advantages: second-order accuracy in time— detachability of the incompressibility and the nonlinearity—unconditional stability. In practice, at each temporal iteration, two generalized Stokes equations and two linear elliptic equations with variable coefficients must be solved, equations which then must be discretized in space. ∆t is the time step—and a staggered marker-and-cell (MAC) mesh with second-order finitedifferences [33] is used for the spatial discretization—h is the spatial mesh size. An orthogonal projection algorithm [36], the Conjugate Gradient (CG) method [12], the Fast Fourier Transform (FFT) method [43, 17] and SuperLU [14] are used to solve the generalized Stokes equation. The orthogonal projection algorithm transforms the generalized Stokes equation into a linear system whose coefficients matrix is symmetric, positive semidefinite and remains fixed throughout the temporal iterations. But yet, it includes a projector as a factor. This linear system is then solved by the CG method. The FFT method and SuperLU are used to compute the action of the projector. The variable coefficients in the linear elliptic equation come from the discretization of the nonlinear terms. The coefficients matrix is then nonsymmetric and moreover varies throughout the temporal iterations. Furthermore, the incremental unknowns method [44, 6, 7, 8, 9, 10, 11, 5, 18, 19, 20, 21, 35, 22, 23, 28, 40, 41] acts as a spatial preconditioner over the linear system. It is then solved by the preconditioned Bi-CGSTAB method [46], using in addition the block diagonal scaling preconditoner where the variables coefficients are discarded. This preconditioner remains then fixed throughout the temporal iterations, and furthermore the Fast Fourier Transform method [42, 17] is used to compute it. The drawback is that the discarding of the variable coefficients is not efficient at extreme Reynolds numbers—except that the time step diminishes. This is the reason why in the end ∆t = h. The combination of known methods before produces a Direct Numerical Simulation (DNS) which simulates the flow. It runs in double precision from t = 0 to a sufficiently long time t∞ . The spatial mesh size h = 1/256 is used throughout. As many as needed DNS’s are produced and then direct 3
observations of them are at the core of the temporal limit’s study. Thus far, Reynolds numbers from Re = 102 to Re = 109 have been considered [24, 25, 26, 27]. As the Reynolds number Re −→ ∞, the temporal limit [13] (at the Reynolds number Re fixed and as the time t −→ ∞) of the lid-driven square cavity flow evolves from stationary to periodic to aperiodic [32] and then to chaotic. The primary eddy switches from clockwiseto counterclockwise-rotating. The primary eddy is the eddy leading the dynamic. The primary eddy’s switching from clockwise- to counterclockwise-rotating does not mean that the original primary eddy switches direction of rotation. But it means that the original primary eddy reduces in size, becoming a thin strip attached to the top wall, extending from the top left corner to the top right corner, with a small clockwise-rotating eddy attached to the strip at the top right corner. And, meanwhile, the adjacent counterclockwise-rotating eddies coalesce occupying a larger and larger extension of the square cavity’s central region, making in addition the adjacent clockwise-rotating eddies to circulate counterclockwise around it. This wide counterclockwise-rotating eddy becomes the actual primary eddy, producing the primary eddy’s switching direction of rotation. Furthermore, the wide counterclockwise-rotating eddy may exhibit break-up and subsequent coalescence—in this temporally evolving flow. But yet, two ultimate questions deserve further attention. In the first place, when the primary eddy first switches from clockwise- to counterclockwise-rotating? The switching first is observed at Re = 307, 000 and, at Re = 310, 000 and beyond, unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. In the second place, how the temporal limit behaves at extreme Reynolds 1 numbers? Here, Re = O( ) = O(1015 ). And then, besides discretizaeps tion errors, round-off errors are also significant and may further distort the inherent flow in the lid-driven square cavity. All the same, an interesting dynamic still endures. Notwithstanding, it is not certain that the dynamic observed at high and extreme Reynolds numbers represents indeed true physics. But yet, it appears to agree with the theory (see [16, page 13]) where wide eddies of each kind form out as the fluid evolves. On the other hand, DNS’s at high and extreme Reynolds numbers appear scarce—if not absent. For instance, Reference [48] reaches out a DNS at 4
Re = 106 where many small counterclockwise-rotating eddies stand up at the square cavity’s central region, a fact disagreeing with ours (see [27]). Once and again, a benchmark for checking the accuracy and efficiency of NSE solvers [29, 37, 47, 3, 34, 31, 1, 15, 22, 4, 24, 25, 26], the lid-driven square cavity flow [38] becomes a benchmark for checking two-dimensional chaos in fluids [16]. This article is organized as follows. Section 2 highlights the facts from Re = 102 to Re = 109 . Section 3 considers the primary eddy’s switching from clockwise- to counterclockwise-rotating. Section 4 deals with flows at extreme Reynolds numbers: Re = 1014 , 1015 , 1016 , 1017 . And Section 5 provides the conclusion. 2. The facts: from Re = 102 to Re = 109 First of all, four principles appear [26], [27] at the genesis of small counterclockwiseor clockwise-rotating eddies: First principle A strong, still clockwise- or counterclockwise-rotating eddy attached to a rigid wall leaves it somewhere to turn away from it, arousing nearby, within itself, a tiny counterclockwise- or clockwiserotating eddy attached to it—perhaps stirring this up. Second principle The interaction of a clockwise- or counterclockwise-rotating eddy and another or itself on a counterclockwise- or clockwise-rotating eddy may split this out—whose parts may then merge owing to the same interaction. Third principle The currents nearby within a clockwise- or counterclockwiserotating eddy eventually may create a bend. It may then influence the globally clockwise- or counterclockwise-rotating fluid inside it to whirl locally counterclockwise or clockwise, a whirling which may further arouse deeper inside it a tiny, loose counterclockwise- or clockwiserotating eddy. Fourth principle A strong, still clockwise-rotating eddy attached to the moving (lid) top wall leaves it somewhere to turn away from it, arousing closely rightward, within itself, near attached to it, a tiny, loose counterclockwise-rotating eddy—perhaps stirring this up. This principle resembles the First principle but for a moving wall. 5
Figure 2 concretizes the first principle at Re = 5, 000. Here, exceptionally, h = 1/1024 is the spatial mesh size. The use of this fine mesh size allows to compute tiny quaternary eddies at the bottom right and bottom left corners, which otherwise—for coarser spatial mesh sizes—are missed. The First principle appears at the genesis of all the small eddies at the three relevant corners of the square cavity. Specifically, the primary eddy attaches to the right wall, the bottom wall, and the left wall. It leaves the right wall a little bit below the center of the right wall to turn leftwards, arousing below a small secondary eddy attached to both the right wall and the bottom wall. It leaves the bottom wall a little bit past the center of the bottom wall to turn upwards, arousing leftwards a small secondary eddy attached to both the bottom wall and the left wall. And it leaves the left wall a little bit above the center of the left wall to turn rightwards, arousing above a small secondary eddy attached to the left wall. Furthermore, at the bottom right corner the secondary eddy leaves the bottom wall a little bit to the left of the bottom right corner to turn upwards, arousing rightwards a small tertiary eddy attached to both the bottom wall and the right wall. In turn, this leaves the right wall just above the bottom right corner to turn leftwards, arousing below a tiny quaternary eddy attached to both the right wall and the bottom wall. On the same hand, at the bottom left corner the secondary eddy leaves the left wall a little bit above the bottom left corner to turn rightwards, arousing below a small tertiary eddy attached to the left wall but detached from the bottom wall—a feeble quaternary eddy stands up at the bottom left corner merging along the bottom wall with the secondary eddy. Figure 3, the second at Re = 12, 000. A small secondary eddy stands up at the bottom right corner attached to both the right wall and the bottom wall with two secondary eddies surfacing within it. The interaction of the primary eddy with the tertiary eddy at the bottom right corner on it is about to split it out. However, the splitting do not occur at Re = 12, 000 but certainly does at Re = 12, 500. Figure 4, the third at Re = 30, 000. At the bottom left corner clockwiserotating fluid within the primary eddy turns upwards and immediately falls downwards creating a clockwise-rotating bend. It influences the fluid within it to whirl locally counterclockwise, giving rise to a tiny, loose counterclockwiserotating eddy. Figure 5, the fourth at Re = 5 · 108 . A widespread clockwise-rotating eddy, a strip attached to the moving (lid) top wall, stretches out from the 6
left wall to past the center of the moving (lid) top wall with, horizontally, two clockwise-rotating eddies surfacing within it. The one on the left leaves the moving (lid) top wall a little bit to the right of the top left corner, arousing closely rightwards, within itself, near attached to the moving (lid) top wall, a tiny, loose counterclockwise-rotating eddy. Thus far, ranging from Re = 102 to Re = 109 , as many as needed DNS’s are produced to watch the nature of the temporal limit [24, 25, 26, 27]. Firstly, Reference [24] includes at low Reynolds numbers whenever possible a detailed comparison of quantities with results published by others, these sometimes agree—at other times disagree. Secondly, Reference [25] considers the Hopf bifurcations. Up to Re = 7, 307.75 the temporal limit is stationary. But at Re = 7, 308 the temporal limit is time periodic, not stationary. So the critical Reynolds number for the first Hopf bifurcation localizes between Re = 7, 307.75 and Re = 7, 308, a value smaller than the one reported elsewhere (see [2, 30] and the references therein). On the same hand, at Re = 13, 393.5 the temporal limit is time periodic. But at Re = 13, 393.75 the temporal limit is not time periodic anymore: losing unambiguously, abruptly time periodicity, it becomes aperiodic. So the critical Reynolds number for the second Hopf bifurcation localizes between Re = 13, 393.5 and Re = 13, 393.75. Also, Reference [25] considers the effect of the spatial-mesh-size refining. Several spatial mesh sizes further and further refined are considered: h = 1/128, 1/256, 1/512, 1/1024. For low Reynolds numbers, while not being reached yet, mesh convergence of temporal limits seems to take place. It should appear for h equal to or smaller than h = 1/2048, however. No sign of stagnation in the computations is observed. But yet, those DNS’s are out of reach for chaotic temporal limits because of the computational costs involved. Thirdly, Reference [26] extends the study up to Re = 500, 000. Up to Re = 200, 000 the temporal limit is aperiodic. Always the primary eddy rotates clockwise. But yet, at Re = 500, 000 the temporal limit is chaotic. At the upside of the total kinetic energy range, at first the primary eddy rotates clockwise, but then a competition for becoming the primary eddy sets up, the counterclockwise-rotating eddies widening enough, until a counterclockwiserotating primary eddy prevails. At the downside, the primary eddy rotates counterclockwise. At the mid side, at first the primary eddy rotates counterclockwise, but then a competition for becoming the primary eddy sets up, the clockwise-rotating eddies widening enough, until a clockwise-rotating pri7
mary eddy prevails. This alternating behavior persists as the time t −→ ∞. And then, Reference [27] further extends the study up to Re = 109 . As the Reynolds number Re −→ ∞, it persists up to Re = 5 · 108 . But yet, at Re = 109 the primary eddy—at its qualitative temporal limit’s pace— appears to rotate solely counterclockwise. Figure 6 displays the kinetic energy from t = 0 to t∞ = 150, 000. The time step ∆t = h Here, the fluid starts from data: the initial condition is the numerical flow computed at some time at a prior Reynolds number. Finally, this research concludes the study, reaching out a DNS at Re = 1017 . 3. The primary eddy’s switching from clockwise- to counterclockwiserotating From now on, the kinetic energy displays from t = 0 to t∞ = 50, 000. The time step ∆t = 2h. Here, the fluid starts from rest. At Re = 300, 000 Figure 7 displays the kinetic energy for several restartings. The primary eddy appears to rotate always clockwise on the total kinetic energy range. Besides, no lowest spike appears noticeably therein. Moreover, Figure 10 (top left) displays the flow which stands at t∞ = 50, 000: the primary eddy rotates clockwise. But then, a primary eddy’s switching from clockwise- to counterclockwiserotating first is observed at Re = 307, 000 at the upper side of the total kinetic energy range. Afterward, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range, and unexpectedly it relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. Figure 8 displays the kinetic energy, and Figure 10 (bottom left) displays the flow which stands at t∞ = 50, 000: the primary eddy rotates clockwise. At Re = 307, 000 the momentary counterclockwise-rotating primary eddy’s prevailing at the upper side of the total kinetic energy range displays in Figure 11. It is remarkable the limited extension of the counterclockwise-rotating primary eddy. At Re = 310, 000 the momentary counterclockwise-rotating primary eddy’s prevailing around the unique lowest spike in the total kinetic energy range dis-
8
plays in Figure 12. It is remarkable the wide extension of the counterclockwiserotating primary eddy. Furthermore, at Re = 305, 000 the flow behaves as at Re = 300, 000. And at Re = 400, 000 the flow reinforces the flow’s characteristics at Re = 310, 000. See Figure 9 and Figure 10. 4. Flows at extreme Reynolds numbers Here, the time step ∆t = h. The fluid starts from data. 1 For Reynolds numbers larger than Re = O( 2 ) = O(65, 536), discretizah tion errors are significant and may distort the inherent flow in the liddriven square cavity. Notwithstanding, an interesting dynamic endures at Re = 106 , 108 , 5 · 108 , 109 [27]. But yet, what happens at extreme Reynolds numbers where, besides discretization errors, round-off errors are also significant? 1 Let us consider the cases Re = 1014 , 1015 , 1016 , 1017 . Here, Re = O( )= eps O(1015). Figure 13 displays the kinetic energy, and Figure 14 displays the flow which stands at t∞ = 50, 000: the primary eddy rotates counterclockwise. At Re = 1014 , 1016 the primary eddy rotates mostly counterclockwise, but at the tall spikes it rotates clockwise. At Re = 1015 , 1017 the primary eddy appears to rotate solely counterclockwise. Figure 15 displays the flow at Re = 1017 past the time t∞ = 50, 000. Again, it is remarkable the wide extension of the counterclockwise-rotating primary eddy. In other words, in the end, the adjacent counterclockwise-rotating eddies coalesce into a wide, prevailing counterclockwise-rotating eddy in the square cavity’s central region—the actual primary eddy. The adjacent clockwiserotating eddies circulating counterclockwise around it are unable to coalesce and to break it up. It seems that the eddies of each kind coalesce to attain the maximum extension available. But there is not enough room in the unit square cavity for a wide eddy of each kind to stand up as the fluid evolves. Discretization and round-off errors seem to modify somehow the Reynolds number, and although they are significant—all the same—an interesting dynamic still endures.
9
5. Conclusion As the Reynolds number Re −→ ∞ the primary eddy switches from clockwise- to counterclockwise-rotating, a switching repeatedly reversing as the time t −→ ∞. The primary eddy is the eddy leading the dynamic. Marking the genesis of chaos, the switching first is observed at Re = 307, 000 at the upper side of the total kinetic energy range. But then, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range and unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. At extreme Reynolds numbers 1 Re = O( ) = O(1015 ), where eps is the machine epsilon in double precieps sion, as the time t −→ ∞, an interesting dynamic still endures where the switching repeatedly reverses or where it stays put, the primary eddy rotating solely counterclockwise. Acknowledgments This research was supported in part by the National Science Foundation through grants DMS-0906440 and DMS-1206438. References [1] E. Barragy, G. F. Carey, Stream function-vorticity driven cavity solution using p finite elements, Comput. & Fluids 26 (5) (1997) 453–468. [2] A. Brezillon, G. Girault, J. M. Cadou, A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics, Comput. & Fluids 39 (7) (2010) 1226–1240. [3] C.-H. Bruneau, C. Jouron, An efficient scheme for solving steady incompressible Navier-Stokes equations, J. Comput. Phys. 89 (2) (1990) 389–413. [4] C.-H. Bruneau, M. Saad, The 2D lid-driven cavity problem revisited, Comput. & Fluids 35 (3) (2006) 326–348.
10
[5] M. Chen, A. Miranville, R. Temam, Incremental unknowns in finite differences in three space dimensions, Comput. Appl. Math. 14 (3) (1995) 219–252. [6] M. Chen, R. Temam, The incremental unknown method I, Appl. Math. Lett. 4 (3) (1991) 73–76. [7] M. Chen, R. Temam, The incremental unknown method II, Appl. Math. Lett. 4 (3) (1991) 77–80. [8] M. Chen, R. Temam, Incremental unknowns for solving partial differential equations, Numer. Math. 59 (3) (1991) 255–271. [9] M. Chen, R. Temam, Incremental unknowns for convection-diffusion equations, Appl. Numer. Math. 11 (5) (1993) 365–383. [10] M. Chen, R. Temam, Incremental unknowns in finite differences: Condition number of the matrix, SIAM J. Matrix Anal. Appl. 14 (2) (1993) 432–455. [11] M. Chen, R. Temam, Nonlinear Galerkin method in the finite difference case and wavelet-like incremental unknowns, Numer. Math. 64 (3) (1993) 271–294. [12] P. Concus, G. H. Golub, D. P. O’Leary, A generalized conjugate gradient method for the numerical solution of elliptic partial differential equations, Sparse Matrix Computations, Academic Press, 1976, J. R. Bunch and D. J. Rose (Eds.). [13] P. Constantin, A few results and open problems regarding incompressible fluids, Notices Amer. Math. Soc. 42 (6) (1995) 658–663. [14] J. W. Demmel, S. C. Eisenstat, J. R. Gilbert, X. S. Li, J. W. H. Liu, A supernodal approach to sparse partial pivoting, SIAM J. Matrix Analysis and Applications 20 (3) (1999) 720–755. [15] E. Erturk, T. C. Corke, C. G¨ok¸c¨ol, Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, Internat. J. Numer. Methods Fluids 48 (7) (2005) 747–774. [16] C. Eugene Wayne, Vortices and two-dimensional fluid motion, Notices Amer. Math. Soc. 58 (1) (2011) 10–19. 11
[17] M. Frigo, S. G. Johnson, The design and implementation of FFTW3, Proc. IEEE 93 (2) (2005) 216–231. [18] S. Garcia, The matricial framework for the incremental unknowns method, Numer. Funct. Anal. Optim. 14 (1 & 2) (1993) 25–44. [19] S. Garcia, Numerical study of the incremental unknowns method, Numer. Methods Partial Differential Equations 10 (1) (1994) 103–127. [20] S. Garcia, Higher-order incremental unknowns, hierarchical basis, and nonlinear dissipative evolutionary equations, Appl. Numer. Math. 19 (4) (1996) 467–494. [21] S. Garcia, Algebraic conditioning analysis of the incremental unknowns preconditioner, Appl. Math. Modelling 22 (4–5) (1998) 351–366. [22] S. Garcia, Incremental unknowns for solving the incompressible NavierStokes equations, Math. Comput. Simulation 52 (5–6) (2000) 445–489. [23] S. Garcia, Incremental unknowns and graph techniques in three space dimensions, Appl. Numer. Math. 44 (3) (2003) 329–365. [24] S. Garcia, The lid-driven square cavity flow: From stationary to time periodic and chaotic, Commun. Comput. Phys. 2 (5) (2007) 900–932. [25] S. Garcia, Hopf bifurcations, drops in the lid-driven square cavity flow, Adv. Appl. Math. Mech. 1 (4) (2009) 546–572. [26] S. Garcia, Aperiodic, chaotic lid-driven square cavity flows, Math. Comput. Simulation 81 (9) (2011) 1741–1769. [27] S. Garcia, Chaotic lid-driven square cavity flows at extreme Reynolds numbers, Commun. Comput. Phys. 15 (3) (2014) 596–617. [28] S. Garcia, F. Tone, Incremental unknowns and graph techniques with in-depth refinement, Int. J. Numer. Anal. Model. 4 (2) (2007) 143–177. [29] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48 (3) (1982) 387–411.
12
[30] A. N. Gorban, D. J. Packwood, Enhancement of the stability of lattice Boltzmann methods by dissipation control, Phys. A 414 (2014) 285–299. [31] O. Goyon, High-Reynolds number solutions of Navier-Stokes equations using incremental unknowns, Comput. Methods Appl. Mech. Engrg. 130 (3–4) (1996) 319–335. [32] J. Guckenheimer, Strange attractors in fluids: Another view, Annu. Rev. Fluid Mech. 18 (1986) 15–31. [33] F. H. Harlow, J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (12) (1965) 2182–2189. [34] M. Li, T. Tang, B. Fornberg, A compact fourth-order finite difference scheme for the steady incompressible Navier-Stokes equations, Internat. J. Numer. Methods Fluids 20 (10) (1995) 1137–1151. [35] P. Poullet, Staggered incremental unknowns for solving Stokes and generalized Stokes problems, Appl. Numer. Math. 35 (1) (2000) 23–41. [36] V. Sarin, A. Sameh, An efficient iterative method for the generalized Stokes problem, SIAM J. Sci. Comput. 19 (1) (1998) 206–226. [37] R. Schreiber, H. B. Keller, Driven cavity flows by efficient numerical techniques, J. Comput. Phys. 49 (2) (1983) 310–333. [38] P. N. Shankar, M. D. Deshpande, Fluid mechanics in the driven cavity, Annu. Rev. Fluid Mech. 32 (2000) 93–136. [39] A. Smith, D. Silvester, Implicit algorithms and their linearisation for the transient incompressible Navier-Stokes equations, IMA J. Numer. Anal. 17 (4) (1997) 527–545. [40] L. Song, Y. Wu, Incremental unknowns in three-dimensional stationary problem, Numer. Algorithms 46 (2) (2007) 153–171. [41] L. Song, Y. Wu, Incremental unknowns method based on the θ-scheme for time-dependent convection-diffusion equations, Math. Comput. Simulation 79 (7) (2009) 2001–2012.
13
[42] P. N. Swarztrauber, The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson’s equation on a rectangle, SIAM Rev. 19 (3) (1977) 490–501. [43] R. A. Sweet, Direct methods for the solution of Poisson’s equation on a staggered grid, J. Comput. Phys. 12 (3) (1973) 422–428. [44] R. Temam, Inertial manifolds and multigrid methods, SIAM J. Math. Anal. 21 (1) (1990) 154–178. [45] R. Temam, Navier-Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing. American Mathematical Society, 2001. [46] H. A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 13 (2) (1992) 631–644. [47] S. P. Vanka, Block-implicit multigrid solution of Navier-Stokes equations in primitive variables, J. Comput. Phys. 65 (1) (1986) 138–158. [48] C. Zhen-Hua, S. Bao-Chang, Z. Lin, Simulating high Reynolds number flow in two-dimensional lid-driven cavity by multi-relaxation-time lattice Boltzmann method, Chinese Phys. 15 (8) (2006) 1855–1863.
14
u 3
2
1
3 Figure 1:
2
2
3
4
The lid-driven square cavity flow. The stationary configuration.
15
Figure 2:
Re = 5, 000. h = 1/1024. ∆t = 4h. Stationary temporal limit.
16
Figure 3:
Re = 12, 000. ∆t = 4h. Time periodic temporal limit.
17
Figure 4:
Re = 30, 000. ∆t = 3h. Aperiodic temporal limit.
18
Figure 5:
Re = 5 · 108 . ∆t = h. Chaotic temporal limit.
19
8
Re = 10
8
Re = 5 ⋅ 10
Re = 109
Figure 6:
The kinetic energy from t = 0 to t∞ = 150, 000. The time step ∆t = h.
20
Re = 300,000
Figure 7:
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = 2h.
21
Re = 307,000
Re = 310,000
Figure 8:
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = 2h.
22
Re = 305,000 Re = 310,000
Figure 9:
Re = 300,000
Re = 400,000
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = 2h.
23
Figure 10:
Re = 300,000
Re = 305,000
Re = 310,000
Re = 400,000
The flow at t∞ = 50, 000. A clockwise-rotating primary eddy’s standing.
24
Figure 11:
Re = 307, 000. A counterclockwise-rotating primary eddy’s prevailing.
25
Figure 12:
Re = 310, 000. A counterclockwise-rotating primary eddy’s prevailing.
26
16
Re = 10
14
Re = 10
17
Re = 10
Figure 13:
15
Re = 10
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = h.
27
Figure 14:
standing.
Re = 1014
Re = 1015
Re = 1016
Re = 1017
The flow at t∞ = 50, 000. A counterclockwise-rotating primary eddy’s
28
Figure 15:
Re = 1017 . A counterclockwise-rotating primary eddy’s prevailing.
29
PII: DOI: Reference:
S0378-4754(17)30186-6 http://dx.doi.org/10.1016/j.matcom.2017.04.010 MATCOM 4457
To appear in:
Mathematics and Computers in Simulation
Received date : 30 June 2015 Revised date : 10 April 2017 Accepted date : 26 April 2017 Please cite this article as: S. Garcia, Chaos in the lid-driven square cavity, Math. Comput. Simulation (2017), http://dx.doi.org/10.1016/j.matcom.2017.04.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Chaos in the lid-driven square cavity Salvador Garcia Instituto de Ciencias F´ısicas y Matem´ aticas Universidad Austral de Chile Casilla 567, Valdivia Chile
Abstract As the Reynolds number Re −→ ∞ the primary eddy switches from clockwiseto counterclockwise-rotating, a switching repeatedly reversing as the time t −→ ∞. The primary eddy is the eddy leading the dynamic. Marking the genesis of chaos, the switching first is observed at Re = 307, 000 at the upper side of the total kinetic energy range. But then, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range and unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. At extreme Reynolds numbers 1 Re = O( ) = O(1015 ), where eps is the machine epsilon in double precieps sion, as the time t −→ ∞, an interesting dynamic still endures where the switching repeatedly reverses or where it stays put, the primary eddy rotating solely counterclockwise. Keywords: Navier-Stokes equations, lid-driven square cavity flow, chaos 1. Introduction As the Reynolds number Re −→ ∞ the primary eddy switches from clockwise- to counterclockwise-rotating, a switching repeatedly reversing as the time t −→ ∞. The primary eddy is the eddy leading the dynamic. Marking the genesis of chaos, the switching first is observed at Re = 307, 000 at Email address: [email protected] (Salvador Garcia) Preprint submitted to Mathematics and Computers in Simulation
April 10, 2017
the upper side of the total kinetic energy range. But then, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range and unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. At extreme Reynolds numbers 1 Re = O( ) = O(1015 ), where eps is the machine epsilon in double precieps sion, as the time t −→ ∞, an interesting dynamic still endures where the switching repeatedly reverses or where it stays put, the primary eddy rotating solely counterclockwise. Chaos appears as a systematic deterioration of a stationary configuration, while always obeying a definite set of principles [26, 27] and while in the end reversing somehow such a stationary configuration. The domain is the unit square cavity. And the viscous incompressible flow is governed by the two-dimensional time-dependent incompressible NavierStokes equations (NSE) [45] and driven by the upper wall which moves from left to right, see Figure 1. The nondimensionalized NSE in primitive variables with Dirichlet boundary conditions over the domain Ω = ]0, 1[×]0, 1[ read ∂u 1 − ∆u + (u · ∇)u + ∇p = f in Ω, t > 0, ∂t Re ∇ · u = 0 in Ω, (1) u = ϕ on Γ = ∂Ω, u(x, 0) = u0 (x) in Ω, where u is the velocity; p, the pressure; Re, the Reynolds number; f, the external force. Here, f = 0 and the boundary conditions read ( u(x, ·) = (1, 0) if x ∈ upper wall, (2) u(x, ·) = 0 if x ∈ left, bottom, or right wall.
At first, an unexpected balance of viscous and pressure forces makes the fluid to turn into the unit square cavity. The properties of these forces depending upon the Reynolds number, a hierarchy of eddies develops. First, the large clockwise-rotating primary eddy (1) arises. Its location occurs toward the geometric center of the unit square cavity. And then several small eddies arise: the counterclockwise-rotating secondary eddies (2), the clockwise-rotating tertiary eddies (3), the counterclockwise-rotating quaternary eddies (4), the clockwise-rotating quinary eddies (5). Their locations 2
occur at the three relevant corners of the square cavity: bottom right (BR), bottom left (BL), and top left (TL), appearing hierarchically at the inclined ellipses in Figure 1. This is the stationary configuration. A combination of known methods is used to discretize and solve the NSE (see [24] for further details). But yet, to summarize, the linear Linθ∗ scheme [39] is the frame of work for the temporal discretization of the NSE, a scheme possessing the following advantages: second-order accuracy in time— detachability of the incompressibility and the nonlinearity—unconditional stability. In practice, at each temporal iteration, two generalized Stokes equations and two linear elliptic equations with variable coefficients must be solved, equations which then must be discretized in space. ∆t is the time step—and a staggered marker-and-cell (MAC) mesh with second-order finitedifferences [33] is used for the spatial discretization—h is the spatial mesh size. An orthogonal projection algorithm [36], the Conjugate Gradient (CG) method [12], the Fast Fourier Transform (FFT) method [43, 17] and SuperLU [14] are used to solve the generalized Stokes equation. The orthogonal projection algorithm transforms the generalized Stokes equation into a linear system whose coefficients matrix is symmetric, positive semidefinite and remains fixed throughout the temporal iterations. But yet, it includes a projector as a factor. This linear system is then solved by the CG method. The FFT method and SuperLU are used to compute the action of the projector. The variable coefficients in the linear elliptic equation come from the discretization of the nonlinear terms. The coefficients matrix is then nonsymmetric and moreover varies throughout the temporal iterations. Furthermore, the incremental unknowns method [44, 6, 7, 8, 9, 10, 11, 5, 18, 19, 20, 21, 35, 22, 23, 28, 40, 41] acts as a spatial preconditioner over the linear system. It is then solved by the preconditioned Bi-CGSTAB method [46], using in addition the block diagonal scaling preconditoner where the variables coefficients are discarded. This preconditioner remains then fixed throughout the temporal iterations, and furthermore the Fast Fourier Transform method [42, 17] is used to compute it. The drawback is that the discarding of the variable coefficients is not efficient at extreme Reynolds numbers—except that the time step diminishes. This is the reason why in the end ∆t = h. The combination of known methods before produces a Direct Numerical Simulation (DNS) which simulates the flow. It runs in double precision from t = 0 to a sufficiently long time t∞ . The spatial mesh size h = 1/256 is used throughout. As many as needed DNS’s are produced and then direct 3
observations of them are at the core of the temporal limit’s study. Thus far, Reynolds numbers from Re = 102 to Re = 109 have been considered [24, 25, 26, 27]. As the Reynolds number Re −→ ∞, the temporal limit [13] (at the Reynolds number Re fixed and as the time t −→ ∞) of the lid-driven square cavity flow evolves from stationary to periodic to aperiodic [32] and then to chaotic. The primary eddy switches from clockwiseto counterclockwise-rotating. The primary eddy is the eddy leading the dynamic. The primary eddy’s switching from clockwise- to counterclockwise-rotating does not mean that the original primary eddy switches direction of rotation. But it means that the original primary eddy reduces in size, becoming a thin strip attached to the top wall, extending from the top left corner to the top right corner, with a small clockwise-rotating eddy attached to the strip at the top right corner. And, meanwhile, the adjacent counterclockwise-rotating eddies coalesce occupying a larger and larger extension of the square cavity’s central region, making in addition the adjacent clockwise-rotating eddies to circulate counterclockwise around it. This wide counterclockwise-rotating eddy becomes the actual primary eddy, producing the primary eddy’s switching direction of rotation. Furthermore, the wide counterclockwise-rotating eddy may exhibit break-up and subsequent coalescence—in this temporally evolving flow. But yet, two ultimate questions deserve further attention. In the first place, when the primary eddy first switches from clockwise- to counterclockwise-rotating? The switching first is observed at Re = 307, 000 and, at Re = 310, 000 and beyond, unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. In the second place, how the temporal limit behaves at extreme Reynolds 1 numbers? Here, Re = O( ) = O(1015 ). And then, besides discretizaeps tion errors, round-off errors are also significant and may further distort the inherent flow in the lid-driven square cavity. All the same, an interesting dynamic still endures. Notwithstanding, it is not certain that the dynamic observed at high and extreme Reynolds numbers represents indeed true physics. But yet, it appears to agree with the theory (see [16, page 13]) where wide eddies of each kind form out as the fluid evolves. On the other hand, DNS’s at high and extreme Reynolds numbers appear scarce—if not absent. For instance, Reference [48] reaches out a DNS at 4
Re = 106 where many small counterclockwise-rotating eddies stand up at the square cavity’s central region, a fact disagreeing with ours (see [27]). Once and again, a benchmark for checking the accuracy and efficiency of NSE solvers [29, 37, 47, 3, 34, 31, 1, 15, 22, 4, 24, 25, 26], the lid-driven square cavity flow [38] becomes a benchmark for checking two-dimensional chaos in fluids [16]. This article is organized as follows. Section 2 highlights the facts from Re = 102 to Re = 109 . Section 3 considers the primary eddy’s switching from clockwise- to counterclockwise-rotating. Section 4 deals with flows at extreme Reynolds numbers: Re = 1014 , 1015 , 1016 , 1017 . And Section 5 provides the conclusion. 2. The facts: from Re = 102 to Re = 109 First of all, four principles appear [26], [27] at the genesis of small counterclockwiseor clockwise-rotating eddies: First principle A strong, still clockwise- or counterclockwise-rotating eddy attached to a rigid wall leaves it somewhere to turn away from it, arousing nearby, within itself, a tiny counterclockwise- or clockwiserotating eddy attached to it—perhaps stirring this up. Second principle The interaction of a clockwise- or counterclockwise-rotating eddy and another or itself on a counterclockwise- or clockwise-rotating eddy may split this out—whose parts may then merge owing to the same interaction. Third principle The currents nearby within a clockwise- or counterclockwiserotating eddy eventually may create a bend. It may then influence the globally clockwise- or counterclockwise-rotating fluid inside it to whirl locally counterclockwise or clockwise, a whirling which may further arouse deeper inside it a tiny, loose counterclockwise- or clockwiserotating eddy. Fourth principle A strong, still clockwise-rotating eddy attached to the moving (lid) top wall leaves it somewhere to turn away from it, arousing closely rightward, within itself, near attached to it, a tiny, loose counterclockwise-rotating eddy—perhaps stirring this up. This principle resembles the First principle but for a moving wall. 5
Figure 2 concretizes the first principle at Re = 5, 000. Here, exceptionally, h = 1/1024 is the spatial mesh size. The use of this fine mesh size allows to compute tiny quaternary eddies at the bottom right and bottom left corners, which otherwise—for coarser spatial mesh sizes—are missed. The First principle appears at the genesis of all the small eddies at the three relevant corners of the square cavity. Specifically, the primary eddy attaches to the right wall, the bottom wall, and the left wall. It leaves the right wall a little bit below the center of the right wall to turn leftwards, arousing below a small secondary eddy attached to both the right wall and the bottom wall. It leaves the bottom wall a little bit past the center of the bottom wall to turn upwards, arousing leftwards a small secondary eddy attached to both the bottom wall and the left wall. And it leaves the left wall a little bit above the center of the left wall to turn rightwards, arousing above a small secondary eddy attached to the left wall. Furthermore, at the bottom right corner the secondary eddy leaves the bottom wall a little bit to the left of the bottom right corner to turn upwards, arousing rightwards a small tertiary eddy attached to both the bottom wall and the right wall. In turn, this leaves the right wall just above the bottom right corner to turn leftwards, arousing below a tiny quaternary eddy attached to both the right wall and the bottom wall. On the same hand, at the bottom left corner the secondary eddy leaves the left wall a little bit above the bottom left corner to turn rightwards, arousing below a small tertiary eddy attached to the left wall but detached from the bottom wall—a feeble quaternary eddy stands up at the bottom left corner merging along the bottom wall with the secondary eddy. Figure 3, the second at Re = 12, 000. A small secondary eddy stands up at the bottom right corner attached to both the right wall and the bottom wall with two secondary eddies surfacing within it. The interaction of the primary eddy with the tertiary eddy at the bottom right corner on it is about to split it out. However, the splitting do not occur at Re = 12, 000 but certainly does at Re = 12, 500. Figure 4, the third at Re = 30, 000. At the bottom left corner clockwiserotating fluid within the primary eddy turns upwards and immediately falls downwards creating a clockwise-rotating bend. It influences the fluid within it to whirl locally counterclockwise, giving rise to a tiny, loose counterclockwiserotating eddy. Figure 5, the fourth at Re = 5 · 108 . A widespread clockwise-rotating eddy, a strip attached to the moving (lid) top wall, stretches out from the 6
left wall to past the center of the moving (lid) top wall with, horizontally, two clockwise-rotating eddies surfacing within it. The one on the left leaves the moving (lid) top wall a little bit to the right of the top left corner, arousing closely rightwards, within itself, near attached to the moving (lid) top wall, a tiny, loose counterclockwise-rotating eddy. Thus far, ranging from Re = 102 to Re = 109 , as many as needed DNS’s are produced to watch the nature of the temporal limit [24, 25, 26, 27]. Firstly, Reference [24] includes at low Reynolds numbers whenever possible a detailed comparison of quantities with results published by others, these sometimes agree—at other times disagree. Secondly, Reference [25] considers the Hopf bifurcations. Up to Re = 7, 307.75 the temporal limit is stationary. But at Re = 7, 308 the temporal limit is time periodic, not stationary. So the critical Reynolds number for the first Hopf bifurcation localizes between Re = 7, 307.75 and Re = 7, 308, a value smaller than the one reported elsewhere (see [2, 30] and the references therein). On the same hand, at Re = 13, 393.5 the temporal limit is time periodic. But at Re = 13, 393.75 the temporal limit is not time periodic anymore: losing unambiguously, abruptly time periodicity, it becomes aperiodic. So the critical Reynolds number for the second Hopf bifurcation localizes between Re = 13, 393.5 and Re = 13, 393.75. Also, Reference [25] considers the effect of the spatial-mesh-size refining. Several spatial mesh sizes further and further refined are considered: h = 1/128, 1/256, 1/512, 1/1024. For low Reynolds numbers, while not being reached yet, mesh convergence of temporal limits seems to take place. It should appear for h equal to or smaller than h = 1/2048, however. No sign of stagnation in the computations is observed. But yet, those DNS’s are out of reach for chaotic temporal limits because of the computational costs involved. Thirdly, Reference [26] extends the study up to Re = 500, 000. Up to Re = 200, 000 the temporal limit is aperiodic. Always the primary eddy rotates clockwise. But yet, at Re = 500, 000 the temporal limit is chaotic. At the upside of the total kinetic energy range, at first the primary eddy rotates clockwise, but then a competition for becoming the primary eddy sets up, the counterclockwise-rotating eddies widening enough, until a counterclockwiserotating primary eddy prevails. At the downside, the primary eddy rotates counterclockwise. At the mid side, at first the primary eddy rotates counterclockwise, but then a competition for becoming the primary eddy sets up, the clockwise-rotating eddies widening enough, until a clockwise-rotating pri7
mary eddy prevails. This alternating behavior persists as the time t −→ ∞. And then, Reference [27] further extends the study up to Re = 109 . As the Reynolds number Re −→ ∞, it persists up to Re = 5 · 108 . But yet, at Re = 109 the primary eddy—at its qualitative temporal limit’s pace— appears to rotate solely counterclockwise. Figure 6 displays the kinetic energy from t = 0 to t∞ = 150, 000. The time step ∆t = h Here, the fluid starts from data: the initial condition is the numerical flow computed at some time at a prior Reynolds number. Finally, this research concludes the study, reaching out a DNS at Re = 1017 . 3. The primary eddy’s switching from clockwise- to counterclockwiserotating From now on, the kinetic energy displays from t = 0 to t∞ = 50, 000. The time step ∆t = 2h. Here, the fluid starts from rest. At Re = 300, 000 Figure 7 displays the kinetic energy for several restartings. The primary eddy appears to rotate always clockwise on the total kinetic energy range. Besides, no lowest spike appears noticeably therein. Moreover, Figure 10 (top left) displays the flow which stands at t∞ = 50, 000: the primary eddy rotates clockwise. But then, a primary eddy’s switching from clockwise- to counterclockwiserotating first is observed at Re = 307, 000 at the upper side of the total kinetic energy range. Afterward, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range, and unexpectedly it relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. Figure 8 displays the kinetic energy, and Figure 10 (bottom left) displays the flow which stands at t∞ = 50, 000: the primary eddy rotates clockwise. At Re = 307, 000 the momentary counterclockwise-rotating primary eddy’s prevailing at the upper side of the total kinetic energy range displays in Figure 11. It is remarkable the limited extension of the counterclockwise-rotating primary eddy. At Re = 310, 000 the momentary counterclockwise-rotating primary eddy’s prevailing around the unique lowest spike in the total kinetic energy range dis-
8
plays in Figure 12. It is remarkable the wide extension of the counterclockwiserotating primary eddy. Furthermore, at Re = 305, 000 the flow behaves as at Re = 300, 000. And at Re = 400, 000 the flow reinforces the flow’s characteristics at Re = 310, 000. See Figure 9 and Figure 10. 4. Flows at extreme Reynolds numbers Here, the time step ∆t = h. The fluid starts from data. 1 For Reynolds numbers larger than Re = O( 2 ) = O(65, 536), discretizah tion errors are significant and may distort the inherent flow in the liddriven square cavity. Notwithstanding, an interesting dynamic endures at Re = 106 , 108 , 5 · 108 , 109 [27]. But yet, what happens at extreme Reynolds numbers where, besides discretization errors, round-off errors are also significant? 1 Let us consider the cases Re = 1014 , 1015 , 1016 , 1017 . Here, Re = O( )= eps O(1015). Figure 13 displays the kinetic energy, and Figure 14 displays the flow which stands at t∞ = 50, 000: the primary eddy rotates counterclockwise. At Re = 1014 , 1016 the primary eddy rotates mostly counterclockwise, but at the tall spikes it rotates clockwise. At Re = 1015 , 1017 the primary eddy appears to rotate solely counterclockwise. Figure 15 displays the flow at Re = 1017 past the time t∞ = 50, 000. Again, it is remarkable the wide extension of the counterclockwise-rotating primary eddy. In other words, in the end, the adjacent counterclockwise-rotating eddies coalesce into a wide, prevailing counterclockwise-rotating eddy in the square cavity’s central region—the actual primary eddy. The adjacent clockwiserotating eddies circulating counterclockwise around it are unable to coalesce and to break it up. It seems that the eddies of each kind coalesce to attain the maximum extension available. But there is not enough room in the unit square cavity for a wide eddy of each kind to stand up as the fluid evolves. Discretization and round-off errors seem to modify somehow the Reynolds number, and although they are significant—all the same—an interesting dynamic still endures.
9
5. Conclusion As the Reynolds number Re −→ ∞ the primary eddy switches from clockwise- to counterclockwise-rotating, a switching repeatedly reversing as the time t −→ ∞. The primary eddy is the eddy leading the dynamic. Marking the genesis of chaos, the switching first is observed at Re = 307, 000 at the upper side of the total kinetic energy range. But then, it is observed at Re = 310, 000 at the lower side of the total kinetic energy range and unexpectedly relies upon the occurrence of a unique lowest spike in the total kinetic energy range. However, this reliance is not observed at Re = 307, 000, even when the lowest spike is more pronounced herein, but continues occurring systematically beyond Re = 310, 000. At extreme Reynolds numbers 1 Re = O( ) = O(1015 ), where eps is the machine epsilon in double precieps sion, as the time t −→ ∞, an interesting dynamic still endures where the switching repeatedly reverses or where it stays put, the primary eddy rotating solely counterclockwise. Acknowledgments This research was supported in part by the National Science Foundation through grants DMS-0906440 and DMS-1206438. References [1] E. Barragy, G. F. Carey, Stream function-vorticity driven cavity solution using p finite elements, Comput. & Fluids 26 (5) (1997) 453–468. [2] A. Brezillon, G. Girault, J. M. Cadou, A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics, Comput. & Fluids 39 (7) (2010) 1226–1240. [3] C.-H. Bruneau, C. Jouron, An efficient scheme for solving steady incompressible Navier-Stokes equations, J. Comput. Phys. 89 (2) (1990) 389–413. [4] C.-H. Bruneau, M. Saad, The 2D lid-driven cavity problem revisited, Comput. & Fluids 35 (3) (2006) 326–348.
10
[5] M. Chen, A. Miranville, R. Temam, Incremental unknowns in finite differences in three space dimensions, Comput. Appl. Math. 14 (3) (1995) 219–252. [6] M. Chen, R. Temam, The incremental unknown method I, Appl. Math. Lett. 4 (3) (1991) 73–76. [7] M. Chen, R. Temam, The incremental unknown method II, Appl. Math. Lett. 4 (3) (1991) 77–80. [8] M. Chen, R. Temam, Incremental unknowns for solving partial differential equations, Numer. Math. 59 (3) (1991) 255–271. [9] M. Chen, R. Temam, Incremental unknowns for convection-diffusion equations, Appl. Numer. Math. 11 (5) (1993) 365–383. [10] M. Chen, R. Temam, Incremental unknowns in finite differences: Condition number of the matrix, SIAM J. Matrix Anal. Appl. 14 (2) (1993) 432–455. [11] M. Chen, R. Temam, Nonlinear Galerkin method in the finite difference case and wavelet-like incremental unknowns, Numer. Math. 64 (3) (1993) 271–294. [12] P. Concus, G. H. Golub, D. P. O’Leary, A generalized conjugate gradient method for the numerical solution of elliptic partial differential equations, Sparse Matrix Computations, Academic Press, 1976, J. R. Bunch and D. J. Rose (Eds.). [13] P. Constantin, A few results and open problems regarding incompressible fluids, Notices Amer. Math. Soc. 42 (6) (1995) 658–663. [14] J. W. Demmel, S. C. Eisenstat, J. R. Gilbert, X. S. Li, J. W. H. Liu, A supernodal approach to sparse partial pivoting, SIAM J. Matrix Analysis and Applications 20 (3) (1999) 720–755. [15] E. Erturk, T. C. Corke, C. G¨ok¸c¨ol, Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, Internat. J. Numer. Methods Fluids 48 (7) (2005) 747–774. [16] C. Eugene Wayne, Vortices and two-dimensional fluid motion, Notices Amer. Math. Soc. 58 (1) (2011) 10–19. 11
[17] M. Frigo, S. G. Johnson, The design and implementation of FFTW3, Proc. IEEE 93 (2) (2005) 216–231. [18] S. Garcia, The matricial framework for the incremental unknowns method, Numer. Funct. Anal. Optim. 14 (1 & 2) (1993) 25–44. [19] S. Garcia, Numerical study of the incremental unknowns method, Numer. Methods Partial Differential Equations 10 (1) (1994) 103–127. [20] S. Garcia, Higher-order incremental unknowns, hierarchical basis, and nonlinear dissipative evolutionary equations, Appl. Numer. Math. 19 (4) (1996) 467–494. [21] S. Garcia, Algebraic conditioning analysis of the incremental unknowns preconditioner, Appl. Math. Modelling 22 (4–5) (1998) 351–366. [22] S. Garcia, Incremental unknowns for solving the incompressible NavierStokes equations, Math. Comput. Simulation 52 (5–6) (2000) 445–489. [23] S. Garcia, Incremental unknowns and graph techniques in three space dimensions, Appl. Numer. Math. 44 (3) (2003) 329–365. [24] S. Garcia, The lid-driven square cavity flow: From stationary to time periodic and chaotic, Commun. Comput. Phys. 2 (5) (2007) 900–932. [25] S. Garcia, Hopf bifurcations, drops in the lid-driven square cavity flow, Adv. Appl. Math. Mech. 1 (4) (2009) 546–572. [26] S. Garcia, Aperiodic, chaotic lid-driven square cavity flows, Math. Comput. Simulation 81 (9) (2011) 1741–1769. [27] S. Garcia, Chaotic lid-driven square cavity flows at extreme Reynolds numbers, Commun. Comput. Phys. 15 (3) (2014) 596–617. [28] S. Garcia, F. Tone, Incremental unknowns and graph techniques with in-depth refinement, Int. J. Numer. Anal. Model. 4 (2) (2007) 143–177. [29] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48 (3) (1982) 387–411.
12
[30] A. N. Gorban, D. J. Packwood, Enhancement of the stability of lattice Boltzmann methods by dissipation control, Phys. A 414 (2014) 285–299. [31] O. Goyon, High-Reynolds number solutions of Navier-Stokes equations using incremental unknowns, Comput. Methods Appl. Mech. Engrg. 130 (3–4) (1996) 319–335. [32] J. Guckenheimer, Strange attractors in fluids: Another view, Annu. Rev. Fluid Mech. 18 (1986) 15–31. [33] F. H. Harlow, J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (12) (1965) 2182–2189. [34] M. Li, T. Tang, B. Fornberg, A compact fourth-order finite difference scheme for the steady incompressible Navier-Stokes equations, Internat. J. Numer. Methods Fluids 20 (10) (1995) 1137–1151. [35] P. Poullet, Staggered incremental unknowns for solving Stokes and generalized Stokes problems, Appl. Numer. Math. 35 (1) (2000) 23–41. [36] V. Sarin, A. Sameh, An efficient iterative method for the generalized Stokes problem, SIAM J. Sci. Comput. 19 (1) (1998) 206–226. [37] R. Schreiber, H. B. Keller, Driven cavity flows by efficient numerical techniques, J. Comput. Phys. 49 (2) (1983) 310–333. [38] P. N. Shankar, M. D. Deshpande, Fluid mechanics in the driven cavity, Annu. Rev. Fluid Mech. 32 (2000) 93–136. [39] A. Smith, D. Silvester, Implicit algorithms and their linearisation for the transient incompressible Navier-Stokes equations, IMA J. Numer. Anal. 17 (4) (1997) 527–545. [40] L. Song, Y. Wu, Incremental unknowns in three-dimensional stationary problem, Numer. Algorithms 46 (2) (2007) 153–171. [41] L. Song, Y. Wu, Incremental unknowns method based on the θ-scheme for time-dependent convection-diffusion equations, Math. Comput. Simulation 79 (7) (2009) 2001–2012.
13
[42] P. N. Swarztrauber, The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson’s equation on a rectangle, SIAM Rev. 19 (3) (1977) 490–501. [43] R. A. Sweet, Direct methods for the solution of Poisson’s equation on a staggered grid, J. Comput. Phys. 12 (3) (1973) 422–428. [44] R. Temam, Inertial manifolds and multigrid methods, SIAM J. Math. Anal. 21 (1) (1990) 154–178. [45] R. Temam, Navier-Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing. American Mathematical Society, 2001. [46] H. A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Statist. Comput. 13 (2) (1992) 631–644. [47] S. P. Vanka, Block-implicit multigrid solution of Navier-Stokes equations in primitive variables, J. Comput. Phys. 65 (1) (1986) 138–158. [48] C. Zhen-Hua, S. Bao-Chang, Z. Lin, Simulating high Reynolds number flow in two-dimensional lid-driven cavity by multi-relaxation-time lattice Boltzmann method, Chinese Phys. 15 (8) (2006) 1855–1863.
14
u 3
2
1
3 Figure 1:
2
2
3
4
The lid-driven square cavity flow. The stationary configuration.
15
Figure 2:
Re = 5, 000. h = 1/1024. ∆t = 4h. Stationary temporal limit.
16
Figure 3:
Re = 12, 000. ∆t = 4h. Time periodic temporal limit.
17
Figure 4:
Re = 30, 000. ∆t = 3h. Aperiodic temporal limit.
18
Figure 5:
Re = 5 · 108 . ∆t = h. Chaotic temporal limit.
19
8
Re = 10
8
Re = 5 ⋅ 10
Re = 109
Figure 6:
The kinetic energy from t = 0 to t∞ = 150, 000. The time step ∆t = h.
20
Re = 300,000
Figure 7:
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = 2h.
21
Re = 307,000
Re = 310,000
Figure 8:
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = 2h.
22
Re = 305,000 Re = 310,000
Figure 9:
Re = 300,000
Re = 400,000
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = 2h.
23
Figure 10:
Re = 300,000
Re = 305,000
Re = 310,000
Re = 400,000
The flow at t∞ = 50, 000. A clockwise-rotating primary eddy’s standing.
24
Figure 11:
Re = 307, 000. A counterclockwise-rotating primary eddy’s prevailing.
25
Figure 12:
Re = 310, 000. A counterclockwise-rotating primary eddy’s prevailing.
26
16
Re = 10
14
Re = 10
17
Re = 10
Figure 13:
15
Re = 10
The kinetic energy from t = 0 to t∞ = 50, 000. The time step ∆t = h.
27
Figure 14:
standing.
Re = 1014
Re = 1015
Re = 1016
Re = 1017
The flow at t∞ = 50, 000. A counterclockwise-rotating primary eddy’s
28
Figure 15:
Re = 1017 . A counterclockwise-rotating primary eddy’s prevailing.
29