CC BY-NC-ND
0Available online at www.sciencedirect.com
ScienceDirect
Procedía Engineering
ELSEVIER
Procedía Engineering 61 (2013) 100 - 107
www.elsevier.com/locate/proeedia
Parallel Computational Fluid Dynamics Conference (ParCFD2013)
MRT-LBM simulation of four-lid-driven cavity flow bifurcation
Congshan ZHUOa, Chengwen ZHONGa, Xixiong GUOa, Jun CAOb^
aNational Key Laboratory ofScience and Technology on Aerodynamic Design and Research, Northwestern Polytechnical University, Xi'an 710072, PR China bDepartment of Mechanical and Industrial Engineering, Ryerson University, Toronto, Ontario, Canada, M5B 2K3
Abstract
This paper seeks to make a systematic study over the complex four-lid-driven cavity flows using the multiple-relaxation-time (MRT) lattice Boltzmann method (LBM). The flow is generated by moving the top wall to the right and the bottom wall to the left, while moving the left wall downwards and the right wall upwards, with an identical moving speed. The present MRT-LBM results reveal a lot of important features of bifurcated flow, such as the multiplicity of stable asymmetric and unstable symmetric cavity flow patterns when the Reynolds number exceeds its first critical value (corresponding to the first steady bifurcation), and the second steady bifurcation phenomena on the first unstable solution at the second critical Reynolds number (corresponding to the second steady bifurcation), as well as the flow periodicity after the third critical Reynolds number is reached (referred to as Hopf bifurcation point). The present MRT simulations have predicted that the critical Reynolds numbers are at 359±1 and 721±6 for the second steady bifurcation and the Hopf bifurcation, respectively. For the study of periodic four-lid-driven flows, the stream function and the phase-space trajectory are investigated in detail. Through comparison against the stability analysis and numerical results reported elsewhere, not only does the MRT-LBM approach exhibit its fairly satisfactory accuracy, but also its remarkable capability for investigating the multiplicity of complex flow patterns.
© 2013 The Authors. Published by Elsevier Ltd.
Selection and peer-review under responsibility of the Hunan University and National Supercomputing Center in Changsha (NSCC)
Keywords', lattice Boltzmann method; multiple-relaxation-time model; four-lid-driven cavity flow; bifurcation.
1. Introduction
The lid-driven cavity flows have been considered as a fundamental problem in fluid mechanics and also relevant to a large number of industrial applications [1]. When building a computational model of cavity flow, its geometry is very simple but the flow behavior could be fairly complex that may involve the vortex dynamics [2], flow stability [3], multiplicity [4], and bifurcation [5-6], etc. Hence, the cavity flow has been routinely attracting tremendous attention from the Computational Fluid Dynamics (CFD) community. For instance, the cavity flows driven by two facing or non-facing lids have been extensively investigated both numerically and experimentally in the literatures [7-8], and the results showed that multiple solutions can be generally observed when the Reynolds number is particularly chosen.
Through the finite difference investigation of a series of four-lid-driven cavity flow test cases for Reynolds numbers up to 350, Wahba [5] demonstrating that, when the Reynolds number increases, flow instability becomes increasingly remarkable that typically appears in the form of steady bifurcation provided the Reynolds number does not exceed 129; beyond this critical Reynolds number, multiple solutions of this four-lid-driven cavity flow could arise. More recently, Wahba [9] performed a new series of four-lid-driven cavity flow simulations for the Reynolds number up to 1000, and a Hopf bifurcation was detected at the critical Reynolds number of 735±4, right before the flow becomes periodic. On the other hand, the complex four-lid-driven cavity flow cases were also studied by Cadou et al [10] with focus particularly
* Corresponding author. Tel.: +1-416-979-5000 ext. 7694; fax: +1-416-979-5265. E-mail address: j cao@ryerson. ca
1877-7058 © 2013 The Authors. Published by Elsevier Ltd.
Selection and peer-review under responsibility of the Hunan University and National Supercomputing Center in Changsha (NSCC) doi: 10. 1016/j .proeng.2013.07. 100
placed on flow stability analysis. In [10], the steady bifurcation at a critical Reynolds number of 130 was confirmed; however, from the bifurcated solutions, no Hopf bifurcation was found for a Reynolds number between 130 and 1000; instead, a second steady bifurcation was found on the unstable and stable solutions for a critical Reynolds number close to 360 and 860, respectively. Obviously, these two sets of results exhibit some disagreements on the second steady bifurcation and the Hopf bifurcation.
Compared to conventional numerical methods used in the above cited cavity flow simulations, the lattice Boltzmann method (LBM) has been also frequently applied to the study of flow solution multiplicity. Arumuga and Anoop Dass [6] investigated the multiplicity of steady solutions in two- and four-lid-driven cavity flows using the lattice Bhatnagar-Gross-Krook (LBGK) model. In that study, beyond the critical Reynolds number, the multiple solutions can be obtained by adjusting the lid-driven velocity. However, so far no investigation is found that uses LBM to investigate the second steady bifurcation or the Hopf bifurcation for the four-lid-driven cavity flows. Besides, the multi-relaxation-time (MRT) model, which has shown to be superior to the LBGK model in terms of numerical accuracy and stability [11-12], can render the LBM with free parameters for adjustment in order to achieve the desired flow solution multiplicity; promisingly, the MRT-LBM model should be able to play an important role in complex four-lid-driven cavity flow simulations although it has not been so employed as of today. Also, the existing literature regarding the investigation of flow bifurcation in the four-lid-driven cavity showed somehow significant discrepancy in terms of the second steady bifurcation and the Hopf bifurcation, demanding further study in this direction through other novel numerical approaches.
In the present study, the MRT-LBM model is applied to the examination of complex four-lid-driven cavity flows involving different flow bifurcation phenomena. First, the multiple solutions with asymmetric flow patterns at Re=300 are numerically investigated and then verified through the comparison against the available benchmark results, so that the capabilities of the MRT model can be validated. Then, as the main objective of the present study, the MRT simulations for a range of higher Reynolds numbers (up to 1000) are performed, with the aid of the plots of convergence, stream function, and phase-space trajectory, to investigate the second steady bifurcation as well as the Hopf bifurcation arising in four-lid-driven cavity flows.
The rest of this paper is organized as follows. The MRT-LBM model used in the present simulation work is described in Section 2. Section 3 demonstrates and discusses in detail the present simulation results of multiple solutions and flow bifurcations in the four-lid-driven cavity. Concluding remarks drawn from this study are made in Section 4.
2. Numerical methods
2.1. MRT-LBM model
The Lattice-Boltzmann equation (LBE) with multiple-relaxation-time collision operator is given by [11-12]:
f;.(x+c,S„t+8,)-f,(x,t) = -(M-1 SM)..[fy(x,t)-jy{x,/)] ^ \* MERGEFORMAT (1)
where ft and JJ9 are the discrete distribution function and its equilibrium state, respectively, jris the spatial vector, ct is the discrete velocity in the Ah direction; M is an orthogonal transformation matrix projecting the discrete distribution functions fi and f^ onto the moment spaces m = M/" and m^ = tAf^, and S is a non-negative diagonal matrix with j\ denoting the reciprocal of relaxation time for the Ah moment.
In the present two-dimensional applications, the nine-velocity (D2Q9) model is used, and the discrete velocities et are defined by c0 = (0,0)c, c\ = -c, = (1,0)c, c2 = -c4 = (0,1 )c, cj5 = -c7 = (1,1 )c, c6 = -cg = (-1,1 )c, where c = 8x/8, is the lattice speed with 8X and 8t respectively representing the lattice and time step sizes, which are generally set at 1 in the simulations. In the MRT-D2Q9 model, without presence of external forces, the reduced diagonal matrix appears in the form of S = diag(0,se,se,0,^,0,^,sv,sv) , where the relaxation parameter sv is determined by the dynamic viscosity
fj. = (j-"1 - 0.5)pt?s8,, and all others free relaxation parameters can be adjusted within the interval of (0,2). In this study,
unless otherwise notified, se, se and sq are generally set at 1.2, 1.2 and 1, respectively. More details on M, m and ma? can
be found in Refs. [11-12].
The equilibrium distribution function in Eq. (1) can be expressed as:
Л" = щр
с,и U-и)1 и1 cl 2с.\ 2d
, /=0,1,-,8
\* MERGEFORMAT (2)
where wQ = 4/9 , wx_A = 1/9 , w5_H = 1/36, if = u- u, and cs=c/\f3 representing the lattice sound speed.
Similar to the LBGK model, the MRT evolution process can also be split up into two steps, namely collision and streaming, which can be respectively described by
m' = m-M_1S(m-m^),
\* MERGEFORMAT (3)
fi(x+ci8t,t+8t) = f(x,t),
\* MERGEFORMAT (4)
where f = M m' is the post-collision distribution functions. The macroscopic density and momentum are given by
\* MERGEFORMAT (5)
2.2. Boundary conditions
Fig. 1 Geometry and boundary conditions for a four-lid-driven cavity.
The geometry and boundary conditions of the four-lid-driven cavity flow are shown in Fig. 1. The top wall moves from left to the right and the left wall moves downwards, while the other two parallel faces move in opposite directions. All four driven walls take the same speed ¿70 . In this study, the non-equilibrium extrapolation boundary scheme, which was proposed by Guo et al [13], is adopted for the four moving walls.
3. Numerical results and discussions
All numerical tests performed in this study use Cartesian coordinates with a uniform gird. Throughout this section, the initial fluid density is set at p0 = 1, and the characteristic length L = N8X is determined by the width (equal to the height) of the cavity, which is equal to the numbers of lattices Жоп one side since 8X= 1. Here, the dynamic viscosity is defined by ¡л = pC/0Z/Re, where ¿70 is the aforementioned speed of the moving walls which is set at 0.1 unless otherwise notified.
In the present simulations, a convergence parameter is defined as follows,
[ч 1/2
Y\u(x,f)-u(x,f , \* MERGEFORMAT (6)
which can be used to quantify the convergence and can also help to detect the onset of flow bifurcation.
Also, note that, in this simulation work, the stream function y/ is defined as i// = (J udy+ J - vdxj j2 , and is calculated
by using the Simpson's rule. Here, the two integrations in the x and y direction are computed respectively with the left and bottom walls set as a "datum" (y/ = 0) of the stream function.
All test cases reported in this paper were run on an NVIDIA Tesla CI060 Graphic Processing Unit (GPU) via the Compute Unified Device Architecture (CUDA) interface. By utilizing the explicit parallelism provided by the graphics hardware,computing time can be greatly saved. Compared to the time consumption by using a CPU (Intel(R) Xeon(R) CPU E5520 @ 2.27 GHz), test cases of this present study employing GPU may generally reduce the computational time by 90%.
Fig. 2 The four-lid-driven cavity flows at Re=300: (a) Convergence history; (b) Intermediate symmetry; (c) Final asymmetry A; (d) Final asymmetry B
To demonstrate the existence of multiple steady solutions, a set of test cases at a post-critical Reynolds number, Re = 300, are performed using the present MRT model with the grid resolution of A^160. Here, the wall moving velocity ¿70 is set at 0.01. The convergence history represented by s , defined by Eq. , and three representative flow states obtained in the present simulations are plotted in Fig. 2. As indicated by the whole convergence history in Fig. 2(a), the flow changes gradually from the initial state "1" to an unstable symmetric state "2", as shown in Fig. 2(b), then experiences an intermediate strongly unstable state "3", and eventually tends to the stable asymmetric state "4", as shown in Fig. 2(c) or Fig. 2(d) that is referred to as final asymmetries A or B. Note that the asymmetric A and B states particularly depend on the choice of the MRT initialization parameters, including the lid velocity, free relaxation rates, and boundary condition. The flow patterns of one unstable symmetric and two stable asymmetric states obtained here agree well with those reported in previous literatures [5-6,10].
3.2. Second steady bifurcation
—•— Stable symmetric solution —<>— Stable asymmetric solution
Fig. 3 Bifurcation diagram for the four-lid driven cavity flow (N=160)
Cadou et al. [10] reported that, in addition to the above observed steady bifurcated branches, there exists a second steady bifurcation leading to possibly solutions of both unstable and stable nature provided the Reynolds number is further increased. In order to investigate the second steady bifurcation phenomena, a new series of test cases using an identical mesh with resolution of 7VM60 are examined by gradually increasing the Reynolds number that starts at 300. The entire flow bifurcation diagram over the tested Reynolds numbers ranging from near zero up to 700 is depicted in Fig. 3,
employing the stream function value monitored at the cavity center as the representative quantity of the flow behavior, as suggested by Wahba [5]. The first portion of the diagram with the Reynolds number increasing until an intermediate value, i.e., about Re = 350, reflects the flow behavior in the vicinity of the first critical Reynolds number, which is nearly 150 in this case; accordingly, this first flow bifurcation stage is referred to as "first steady bifurcation Remark that the bifurcated flow at Re = 300 examined in the previous subsection falls into the first steady bifurcation stage that appears with three possible solutions: one unstable symmetric solution (within the dash line in the middle) and two stable asymmetric solutions (the upper and lower lines with circles). Now, the focus is placed on the second critical Reynolds number, which is recognized within a fairly small interval including Re = 350. As shown in Fig. 3, the stable solution appears briefly in the stable form (the short middle solid segment) whereas both upper and lower asymmetric solutions are discontinued accordingly. Right after this brief "pause", the solution multiplicity reoccurs that includes again the three possible forms as aforementioned, which can be referred to as "second steady bifurcation. The numerical experiments performed in this study indicat that the second steady bifurcation phenomena can no longer be observed once the Reynolds number goes up beyond 1000.
U " --Re=350
--Re=352
3 --Re=354
0 "----Re=356
----------------Re=358
n Re=360
0 500 1000 1500 2000
timestep/104
Fig. 4 History of convergence for Reynolds numbers associated with the second steady bifurcation.
For a more detailed study of the second steady bifurcation phenomena, the history of convergence index, s , for the four-lid-driven cavity flow at Re=350, 352, 354, 356, 358 and 360 are plotted in Fig. 4. It can be found in Fig. 4 that, for Reynolds numbers Re=350, 352 and 354, the convergence index s decreases rapidly to a very small value that remains below 10"14, and then appears unceasingly fluctuating with a small amplitude. However, given a further slight increase of Reynolds number (Re=356, 358 and 360), the convergence index s escapes the line of 10"14, and rise up steeply; the more increased Reynolds number, the bigger uprising rate. The present simulation results have shown that all the three possible flow solutions can be captured for Re=356, 358 and 360, suggesting the second steady bifurcation on the unstable symmetric solution of the first steady bifurcated branches should stay at 355±1.
Fig. 5 Streamline patterns for Re=358 at (a) 1000000th step; (b) 16000000th step; and (c) 20000000th step
Representatively, the flow pattern variation for Re=358 is demonstrated in Fig. 5. Note that Fig. 5(a), 5(b), and 5(c) correspond to the states "a", "c", and "d" in Fig. 4, respectively. Also, it can be noticed that the stable asymmetric state (Fig. 5(c)) looks in a similar fashion as in the asymmetric "A" state as shown in Fig. 2(c). Consistently, the asymmetric B state, as shown in Fig. 2(d), is also achievable if the MRT model is re-initialized by using another set of adjusting parameters.
Next, in order to investigate the influence of the lid velocity ¿/0 and the grid resolution N, a new set of MRT tests are performed using three different lid velocities, ¿/0 = 0.02, 0.04, and 0.1, and three different grid resolutions N= 160, 240, and 320. Here, all predicted values of the second critical Reynolds number are grouped in the Table 1. A reasonable conclusion is that the second critical Reynolds number converges towards 359±1. Note the second critical Reynolds number was predicted at Re = 360 by Cadou et al. [10] in the same case, which satisfactorily confirms the effectiveness of the present MRT model for the study of solution multiplicity in complex flow scenarios.
Table 1. The critical Reynolds number for the second steady bifurcation in the four-lid-driven cavity flows at different ¿/0 and N
¿7o=0.02 ¿/o=0.04 ¿/0=0.1
357+1 355 ±1 355 + 1
71^240 359+1 357+1 357+1
7V^320 359+1 359+1 359+1
3.3. Hopf bifurcation
For the four-lid-driven cavity flow, the Hopf bifurcation phenomena, which signify the change from the steady flow regime to the unsteady one, may arise once the Reynolds number goes up beyond its third critical value, which is often referred to as "Hopfbifurcation point". First, in order to investigate the existence of the Hopf bifurcation, a series of tests with grid resolution of 71M60 are performed by gradually increasing the Reynolds number from 400 to 1000. To help observe the Hopf bifurcation phenomena, a monitor point is placed at the cavity center to record the stream function, resulting in Figs. 6(a)-(d) for Reynolds number Re=400, 600, 750 and 1000, respectively. At Re=400 and 600, neither of the plotted stream function behaviors changes at the end, indicating the steady-state of the flow is maintained. However, for the tests at Re=750 and 1000, the periodic behavior looks evident for both, revealing the flows have turned to be unsteady.
Fig. 6 Histories of stream function at the cavity center for different Reynolds numbers: (a) Re=400; (b) Re=600; (c) Re=750; (d) Re=1000
Next, a series of detailed numerical simulations are performed on three grids with resolutions N= 160, 240, and 320 under different lid velocities ¿70 = 0.01, 0.02, 0.03, 0.04, and 0.1. Each test unfolds by gradually increasing the Reynolds number from 696 to 750. All resulting critical Reynolds numbers corresponding to the Hopf bifurcation are grouped in Table 2. Since numerical computation is usually unable to capture the exact single critical Reynolds number for Hopf bifurcation, a small interval that includes the critical Reynolds number in question becomes the hunting target in practice. As shown in Table 2, the Hopf bifurcation points predicted by the MRT model are located at 706±7, 717±6 and 721±6, using the different grids of TIM 60, 240 and 320, respectively. Also, in general, the finer mesh, the better convergence for each of the five lid velocities selected for this simulation series. As a result, the Hopf bifurcation at 721±6 can be identified in the present MRT simulations, which practically agrees pretty well with the Hopf bifurcation point (735±4) reported by Wahba [9].
Table 2. The critical Reynolds numbers for Hopf bifurcation in the four-lid-driven cavity flows at different ¿/0 and N
¿/5=0.01 ¿>6=0.02 ¿yo=0.03 ¿/o=0.04 ¿/5=0.1
N=160 703 705-713 705-711 707 699-709
N=240 715-719 713 715-723 711-715 707
N=320 723 727 719 719 715
Particularly, Fig. 7 resulting from 240 and ¿>&=0.04 is used to illustrate the way of narrowing down towards the interval of the Hopf bifurcation point for the four-lid-driven cavity flows. The simulation results show that the stable states are obtained at Re= 704, 706, 708, 710 and 714, while the unstable (periodic) states are observed at Re= 712, 716, 718, 720 and 722. More interestingly, when Re increases from 710 to 716, the stable and unstable flow solutions appear alternately; moreover, for Re= 704, 706, 708, 710, the flow remains stable whereas, for Re= 716, 718, 720, 722, the flow turns to be unstable. Thus, the critical Reynolds number can be approximately determined within the interval [711,715].
"°'20 100 200 300 400 500 600 "°'20 100 200 300 400
timestep/104 (b) timestep/104
Fig. 7 Histories of stream function at cavity center for (a) Re=710, 714; and (b) Re=712, 716
Alternatively, the phase-space trajectory is another means for identifying the interval of the Hopf bifurcation point. Here, the phase-space trajectories for Re=710 and 714 are respectively plotted in Figs. 8(a) and (b), both featuring the fixed point attractor. On the other hand, the two trajectories for Re=712 and 716 are plotted in Fig. 8(c); the obvious overlapping of the two trajectories clearly suggests a limit cycle attractor and, thus, a periodic flow state.
Fig. 8 Phase-space trajectories for (a) Re=710; (b) Re=714; (c) Re=712 and 716
4. Conclusions
In this paper, the multiple-relaxation-time (MRT) Lattice-Boltzmann method (LBM) is applied to the investigation of the complex four-lid-driven cavity flows for a wide range of Reynolds numbers (up to 1000). The present MRT code is first validated through the simulation of flow solution multiplicity at Re=300. The multiple flow patterns are all captured by the present MRT approach. Then, using further increased Reynolds numbers, the second steady flow bifurcated only from the unstable solution within the first bifurcation stage is also successfully observed in the present MRT simulations, and the second critical Reynolds number is predicted at 359±1, which satisfactorily agrees with the prediction achieved by the stability analysis conducted by others. Above the second critical Reynolds number, three multiple solutions including one unstable symmetric state and two possible asymmetric states can also be obtained. Furthermore, in the present MRT simulations, the four-lid-driven cavity flows are shown to undergo a Hopf bifurcation at the critical Reynolds number of 721 ±6, which also echoes other relevant available numerical results.
The numerical experiments reported here have addressed different sorts of multiple flow solutions, involving the first steady flow bifurcation, the second steady flow bifurcation, and the Hopf bifurcation, and reached good agreement with other available simulation results obtained using other numerical approaches. The effectiveness and robustness of the present MRT model may suggest its extensive applicability towards simulations of flow scenarios with the nature of more complex solution multiplicity.
Acknowledgements
This project was financially supported mainly by the funds from the National Pre-Research Foundation of China (Grant No. 9140A13040111HK0329), the Aeronautical Science Fund of China (Grant No. 20111453012), the opening project of BIT (Beijing Institute of Technology) under grant No.KFJJ10-4M, as well as the Discovery Grant of the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors also wish to thank the Supercomputing Center of the University of Science and Technology of China for providing the computational facilities to conduct the major part of the computations reported in this paper.
References
[1] Aidun, C.K., Triantafillopoulos, N.G., Benson, J.D., 1991. Global stability of a lid-driven cavity with through flow: flow visualization studies, Physics
of Fluids 3, p. 2081.
[2] Cheng, M., Hung, K.C., 2006. Vortex structure of steady flow in a rectangular cavity, Computers & Fluids 35, p. 1046.
[3] Sahin, M, Owens, RG., 2003. A novel fully-implicit finite volume method applied to the lid-driven cavity problem—Part II: Linear stability analysis,
International journal for numerical methods in fluids 42, p. 79.
[4] Shankar, P.N., Deshpande., M.D., 2000. Fluid mechanics in the driven cavity, Annual Review of Fluid Mechanics 32, p. 93.
[5] Wahba, E.M., 2009. Multiplicity of states for two-sided and four-sided lid-driven cavity flows, Computers & Fluids 38, p. 247.
[6] Arumuga Perumal, D., Anoop Dass, K., 2011. Multiplicity of steady solutions in two-dimensional lid-driven cavity flows by Lattice Boltzmann Method,
Computers and Mathematics with Applications 61, p. 3711.
[7] Kuhlmann, H.C., Wanschura, M., Rath, H.J., 1997. Flow in two-sided lid-driven cavities: non-uniqueness, instabilities, and cellular structures, Journal
of Fluid Mechanics 336, p. 267.
[8] Albensoeder, S., Kuhlmann, H.C., Rath, H.J., 2001. Multiplicity of steady two-dimensional flows in two-sided lid-driven cavities, Theoretical
Computational Fluid Dynamics, 14, p. 223.
[9] Wahba, E.M., 2012. Numerical Simulations of Flow Bifurcations Inside a Driven Cavity, CFD Letters 3, p. 100.
[10] Cadou, J.M., Guevel, Y., Girault, G., 2012. Numerical tools for the stability analysis of 2D flows: application to the two- and four-sided lid-driven cavity, Fluid Dynamic Research 44,031403.
[11] Lallemand, P., Luo, L.S., 2000. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance and stability, Physical Review E 61, p. 6546.
[12] Luo, L.S., Liao, W., Chen, X., Peng, Y., Zhang, W., 2011. Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations, Physical Review E 83,056710.
[13] Guo, Z., Zheng, C., Shi, B., 2002 Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Physics 11, p. 366.
