Hyperbolic Method to Explore Multiplicity Flow Solutions in a Four-Sided Lid-Driven Cavity

Abstract

In this study, the hyperbolic method is adopted to explore the flow field states of incompressible flow in a four-sided lid-driven square cavity. In particular, we focus on the flow bifurcation obtained at the critical Reynolds number . In the hyperbolic method, the diffusive term is transformed into a hyperbolic one by introducing a diffusion flux term, which is the solution of an additional equation. A classical Riemann-like solver with a finite-volume discretization is thus employed for the full flux (split into advective and diffusive parts), in order to solve the steady-state incompressible Navier-Stokes equation. The incompressibility of the flow is treated via the artificial pseudo-compressibility method. It is shown that our numerical code is able to detect the bifurcation by the analysis of the residual term relaxation during the pseudo-time iteration procedure. Moreover, depending on the combination choice of slope limiters for the two spatial directions, our method is able to select the first or the second stable solution among the double flow field state obtained when the Reynolds number is higher than the critical value that is estimated to be 129.4 in our study.

Share and Cite:

Baty, H. (2022) Hyperbolic Method to Explore Multiplicity Flow Solutions in a Four-Sided Lid-Driven Cavity. American Journal of Computational Mathematics, 12, 314-323. doi: 10.4236/ajcm.2022.123021.

1. Introduction

The two-dimensional flow driven by boundaries moving inside a cavity is considered to be a classical problem of computational fluid mechanics. In this work, as in many studies, we consider a simple square cavity. This problem displays many interesting solutions, like vortex structures and multiple flow solutions [1] [2] [3] [4] [5].

The different studies have focused on cavity driven by a single plate, two parallel plates, two non-facing plates, or four-sided plates. In the latter configuration, the upper plate moves rightwards, the lower plate moves leftwards, the left plate moves downwards, and the right plate moves upwards. The four plates are also considered to move at the same velocity in the present work. A remarkable result is the existence of a bifurcation for a particular Reynolds number (Re) value close to R c = 130 . Indeed, the steady-state flow changes from a stable symmetric solution for a Reynolds number value lower than Rc to two stable asymmetric solutions for a higher Reynolds number. This is a pitchfork bifurcation as the symmetric solution becomes unstable beyond bifurcation. The value for Rc deduced from the different previous studies is not very well determined, as it varies between 129 and 130.4. A new evaluation is thus needed.

In the present study, we propose to revisit this problem by using the hyperbolic method. The hyperbolization idea is to transform the elliptic diffusive term into a hyperbolic one by introducing the diffusive flux as an additional variable. This new variable is thus the solution of an additional hyperbolic equation. The advantages of this method are: 1) the discretization used for the non-dissipative (inviscid) part of the equations can be directly applicable, 2) a speed-up factor of O ( 1 / h ) is obtained (h being the typical mesh spacing), 3) the diffusive terms are computed to the same order of accuracy as the main solution. We also consider the two-dimensional incompressible Navier-Stokes equation, with the incompressibility of the flow being ensured via the artificial pseudo-compressibility technique. The hyperbolic technique has been successfully applied to diffusion, advection-diffusion, Navier-Stokes, and magnetohydrodynamic equations [6] [7] [8] [9] [10].

The paper is organized as follows. Section 2 is devoted to the model and numerical method. In Section 3, the validation of the scheme and tests are presented for the single-sided cavity case. The results for the four-sided lid-driven cavity are presented in Section 4, with a particular emphasis on the bifurcation of the steady-state flow solution. Finally, conclusions are drawn in Section 5.

2. Model and Numerical Method

2.1. Navier-Stokes Model

The time dependent Navier-Stokes (NS) equation is generally written as,

v t + ( v ) v + P ν 2 v = 0 , (1)

for the velocity flow v , where P is the kinematic thermal pressure (i.e. thermal pressure divided by the fluid density), and ν is the kinematic viscosity coefficient.

Thus, the steady-state (i.e. t 0 , t being the time) Navier-Stokes equation follows,

( v ) v + P ν 2 v = 0. (2)

The incompressibility property must be also imposed via,

v = 0 . (3)

For an homogeneous viscosity, an equivalent conservative form of Equation (2) and Equation (3) is,

div ( v v + P I ¯ ¯ ) div ( ν g ¯ ¯ ) = 0, (4)

with g ¯ ¯ = grad v being the velocity gradient tensor, and where I ¯ ¯ is the identity rank-2 tensor. The last equation makes appear two flux tensors, the advective flux modified by the pressure effect F ¯ ¯ a = v v + P I ¯ ¯ which can be also called the inviscid flux below, and the viscous stress tensor F ¯ ¯ v = ν g ¯ ¯ . The full flux is thus F ¯ ¯ = F ¯ ¯ a + F ¯ ¯ v .

2.2. Artificial Compressibility Method

The previous incompressible Navier-Stokes equation can be solved by introducing a pseudo-time parameter τ in the following compressible equations,

P τ + div ( c 2 v ) = 0, (5)

and,

v τ + div ( v v + P I ¯ ¯ ) div ( ν g ¯ ¯ ) = 0, (6)

with c2 the square value of c, an artificial sound speed. The value of c is arbitrarily chosen in order ensure that the steady-state solution of Equation (5) and

Equation (6) (i.e. τ 0 ) is equivalent to the originel incompressible

Navier-Stokes equation. In practical, the value of c is chosen to be one in this work. This pseudo-incompressible form is said to be equivalent to the incompressible NS equations in the steady-state [11].

2.3. The Hyperbolic Method

It is possible to introduce another differential equation for the stress tensor as,

T r g ¯ ¯ τ grad v = g ¯ ¯ , (7)

where T r is the relaxation parameter defined by T r = L r 2 / ν , with L r a length-scale parameter. The value of L r value is determined in order to achieve an optimal fast convergence [6] [7] [8]. In this way, we have g ¯ ¯ = grad v as expected in the steady-state.

In two dimensions, the system 5 - 7 can be expressed in the following compact form,

U τ + F x + G y = S , (8)

with S being a source term containing the g ¯ ¯ tensor, and where F and G represent x-directed and y-directed full flux respectively.

In order to fix the ideas, we can write a simpler one-dimensional 1D (x dependent) version below, where the fluid velocity flow is given by u ( x ) . Then, the set of equations becomes,

P τ + ( c 2 u ) x = 0, (9)

u τ + ( u 2 + P ν g ) x = 0, (10)

g τ 1 T r u x = g T r . (11)

In this one dimensional version, U is a column vector containing three components, the pressure P, the velocity u, and g the gradient of u. The corresponding flux F has three components that are, c 2 u , u 2 + P ν g , and u / T r . Finally, the source term S has three components, that are zero for the two first ones and g / T r for the third one. It is important for the following numerical implementation to determine the eigenvalues of the inviscid flux Jacobian (i.e. ignoring the viscous terms), that are λ ± = u ± u 2 + c 2 .

2.4. Discretization and Implementation in One Dimension

We consider the 1D version below, as it is simpler to write. For a one-dimensional grid of N nodes with a uniform spacing h, the solution can be evaluated at the nodes denoted by x j , j = 1 , 2 , 3 , , N . Our finite-volume cell-centered discretization is thus,

U j τ = 1 h [ F j + 1 / 2 F j 1 / 2 ] + 1 h I j S d x , (12)

where U j denotes the solution value at the node j, F j + 1 / 2 is an interface flux to be defined (see below), and I j = [ x j 1 / 2 , x j + 1 / 2 ] is a dual cell.

Steady state solutions of the previous system (i.e. with τ 0 ) can be

obtained in different ways. The simplest way consists in marching towards the steady-state using a pseudo-time iteration with a small increment Δ τ . This is similar to a standard scheme (e.g. Euler or Runge-Kutta method) used to integrate in time a general differential equation using a marching procedure with a time-step Δ t . Moreover, in our case, the method is a relaxation procedure towards a steady-state, with an efficiency determined by the choice of the T r parameter (or the associated L r parameter). Steady state solutions of the previous system can be thus obtained by using the simply pseudo-time explicit iteration,

U j n + 1 = U j n Δ τ R e s j n , (13)

where n is the iteration counter, Δ τ is the pseudo-time step, and R e s j n is the residual right-hand side of Equation (12). The latter residual term is required to vanish (up to a pre-defined given accuracy) when a steady state solution is obtained.

The numerical fluxes can be estimated using the standard upwind formula,

F j + 1 / 2 = 1 2 ( F L + F R ) 1 2 c M ( U R U L ) , (14)

where the subscripts L and R stand for the left- and right-hand sides of the cell interface situated at x j + 1 / 2 respectively. The first term is computed from an average value of the two fluxes F L = F ( U L ) and F R = F ( U R ) . Instead of using a full flux directly, the full flux can be separated into two parts (inviscid and viscous parts), for which the above numerical formula can be used separately (see paper [9] ). The two corresponding flux Jacobians have the following eigenvalues λ ± and ± ν / T r for the inviscid and viscous terms respectively. As a consequence, For the inviscid flux, c M is taken as the maximum absolute value of the two eigenvalues λ ± = u ± u 2 + c 2 and represents a maximum inviscid wave speed. For the diffusive flux, c M is simply ν / T r = ν / L r , and can be interpreted as an artificial viscous wave speed. The source term in Equation (12) is computed using a simple point integration approximation, that is shown to not degrade the accuracy of the method.

In the explicit pseudo-time iteration, the value of the increment in the pseudo-time variable τ is limited by a criterion,

Δ τ = C F L × M i n [ h c M + ν / L r ] , (15)

where CFL is the Courant-Friedrichs-Lewy number less than or equal to one, and where we have simply added the maximum wave speed to the dissipative wave speed in the denominator.

2.5. Left/Right Interfaces Interpolation

The left and right values can be evaluated by a linear extrapolation from the cell centre, as for example for the u solution of the 1D discretization,

u L = u j + 1 2 h u j , u R = u j 1 2 h u j , (16)

where u j is the gradient of u j . A simple and economic method is to use g j to evaluate the gradient term, as g is the spatial derivative of u in the steady state. This is a second order scheme that has been previously shown to converge rapidly towards the steady state. Moreover, some classical limited reconstructions options can be easily added to the scheme, that are the superbee, Koren, monotized-central, or minmod limiters.

2.6. Discretization and Implementation in Two Dimensions

The previous one dimensional model can be also easily discretized in two dimensions. The resulting set of hyperbolic equations can be obtained from the work done by Nishikawa (2011) (see [9] ), and are not detailed in this paper. Our algorithm only differs by the splitting procedure of the full Flux (and the corresponding Jacobian) into inviscid and viscous parts.

3. Application to the Lid-Driven Cavity—One Sided Case

As a test case, we first consider the one sided square cavity, where only the upper plate is moving with a velocity V x = V u = 1 (see Figure 1). Using the numerical scheme described above in two dimensions with a spatial domain [ 0 : 1 ] 2 , the resulting solution for the flow field lines is plotted in Figure 1 for a Reynolds number R e = L V u / ν = 1 / ν = 1000 (or a corresponding viscosity ν = 10 3 ). Indeed, we use a normalization where V u = 1 , and L = 1 for the characteristic length scale.

Note that a minmod limiter is chosen for this run. A grid resolution of 130 × 130 uniformly spaced cells is also employed, that is shown to lead to well converged results (see below). Moreover, we have checked that using a higher grid resolution of 180 × 180 does not modify the results. We have also found that

taking L r = 1 20 π , allows an optimal and fast convergence of our results towards

steady-state solutions. The stream lines solution plotted in Figure 1 is similar to the solution expected from previous studies [1] [2] [3] [4]. Indeed, for this relatively high R e value, the two counter rotating vortices are barely visible at the two lower cavity corners.

In order to check the validity and convergence of our results, we have plotted in Figure 2 the cross cut along an horizonthal line at y = 1 / 2 of the V y velocity flow component as a function of x for two Reynolds number values, R e = 1000

and R e = 100 . Note that an optimal value of L r = 1 20 π is chosen in order to

Figure 1. Colored iso-contours of a stream function (deduced from the velocity field) showing the flow stream lines for a one-sided cavity flow, obtained for a run using R e = 1000 . Negative values in the two lower corners indicate reversal of the vortices rotation compared to the rest of the domain.

Figure 2. V y ( x ) profile obtained for y = 1 / 2 , in two runs using R e = 100 and R e = 1000 . A few selected values (see points) extracted from Ghia et al. (1982) are also shown for these two Reynolds number values.

ensure a fast convergence of the results. Our results are compared with the results obtained by Ghia et al. (1982) that is the usual numerical solution taken as the reference in the litterature [12]. As one can see in Figure 2, the agreement is very good.

4. Application to the Lid-Driven Cavity—Four Sided Case

We focus now on the four-sided square cavity case, where the velocity value imposed at each plate is taken to be one. In Figure 3, one can see the solution obtained for R e = 100 , that is in agreement with the expected solution shown in the literature. In order to examine the convergence of the solution towards the steady state, the residual value taken on V x is plotted as a function of the iteration step in Figure 4. Note that, even if a very good approximate solution for the true numerical steady state is already obtained after only 10,000 iterations (the residual being 10−9), it can be even further ameliorated until 40,000 iterations with a residual value of 10−12.

Now, we explore different Reynolds number values by increasing it over the previous value of R e = 100 . The results of the residual term evolution as a function of the iteration counter are plotted in Figure 5 for 4 values of the Reynolds number ranging between 100 and 160. A change in the convergence behavior is clearly visible between R e = 120 and R e = 140 . Indeed, contrary to the monotonic decrease of the residual obtained for the two lowest R e values, a clearly non monotonic variation is seen to precede the final convergence phase for the two highest Reynolds numbers. We have checked that it is in correspondence to the existence of two possible stable steady-state flow solutions, as shown in Figure 6 for R e = 140 . The two previous solutions only differ by the use of different combination of slope limiters in x and y directions. Consequently, we take this change of convergence behavior as a good diagnostic in order to determine

Figure 3. Colored is ocontours of a stream function (deduced from the velocity field) showing the flow stream lines for a four-sided cavity case, obtained for a run using R e = 100 . The direction imposed on each flow plate is indicated.

Figure 4. Residual term taken on one velocity flow component (that is V x ) as a function of the iteration counter, for a run of a four-sided cavity flow using R e = 100 (see previous figure).

Figure 5. Residual term taken on one velocity flow component (that is V x ) as a function of the iteration counter, for runs of a four-sided cavity flow using R e = 100 , 120 , 140 , 160 .

(a) (b)

Figure 6. The two stable steady state flow solutions (stream lines) obtained in two runs of a four-sided cavity flow using R e = 140 . The left and right solutions correspond to the asymmetric-1 and asymmetric-2 solutions respectively, as quoted by Whaba (2011) in [1]. The two simulations only differ by the use of different slope limiters combination in x and y directions.

the critical value of the Reynolds number R c for the bifurcation. We precisely get R c = 129.4 with our procedure.

5. Conclusion

In this study, we use the hyperbolic method in order to study steady-state flow solutions of the four-sided lid-driven square cavity in two dimensions. A classical Riemann solver is used to solve the incompressible Navier-Stokes equation with the artificial compressibility method. A simple second-order finite-volume discretization on a rectangular grid using upwind advective and dissipative fluxes is employed. Our results show a remarkable efficiency to determine the critical Reynolds number R c corresponding to the bifurcation between a single stable solution at a low Reynolds number to a doubly stable one when the Reynolds number is increased. This is precisely obtained by the examination of the convergence behavior of the residual terms during the explicit pseudo-time iteration. We get the critical value R c = 129.4 , which is in the range of values [ 129 : 130.4 ] previously deduced in the literature.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] [1]Wahba, E.M. (2011) Multiplicity of States for Two-Sided and Four-Sided Lid-Driven Cavity Flows. Computers & Fluids, 38, 247-253.
https://doi.org/10.1016/j.compfluid.2008.02.001
[2] Perumal, D.A. and Dass, A.K. (2011) Multiplicity of Steady Solutions in Two-Dimensional Lid-Driven Cavity Flows by Lattice Boltzmann Method. Computers and Mathematics with Applications, 61, 3711-3721.
https://doi.org/10.1016/j.camwa.2010.03.053
[3] Cadou, J.M. and Guevel, Y. and Girault, G. (2012) Numerical Tools for the Stability Analysis of 2D Flows: Application to the Two-and Four-Sided Lid-Driven Cavity. Fluid Dynamics Research, 44, Article ID: 031403.
https://doi.org/10.1088/0169-5983/44/3/031403
[4] Chen, K.T. and Tsai, C.C. and Luo, W.J. and Chen, C.N. (2013) Multiplicity of Steady Solutions in a Two-Sided Lid-Driven Cavity with Different Aspect Ratios. Theoretical and Computational Fluid Dynamics, 27, 767-776.
https://doi.org/10.1007/s00162-013-0296-z
[5] Mohapatra, R.C. (2016) Study on Laminar Two-Dimensional Lid-Driven Cavity Flow with Inclined Side Wall. Open Access Library Journal, 3, Article ID: e2430.
http://dx.doi.org/10.4236/oalib.1102430
[6] Nishikawa, H. (2007) A First-Order System Approach for Diffusion Equation. I: Second-Order Residual-Distribution Schemes. Journal of Computational Physics, 227, 315-352.
https://doi.org/10.1016/j.jcp.2007.07.029
[7] Nishikawa, H. (2014) First-, Second-, and Third-Order Finite-Volume Schemes for Diffusion. Journal of Computational Physics, 256, 791-805.
https://doi.org/10.1016/j.jcp.2013.09.024
[8] Nishikawa, H. (2010) A First-Order System Approach for Diffusion Equation. II: Unification of Advection and Diffusion. Journal of Computational Physics, 229, 3989-4016.
https://doi.org/10.1016/j.jcp.2009.10.040
[9] Nishikawa, H. (2011) First-, Second-, and Third-Order Finite-Volume Schemes for Navier-Stokes Equations. Proceedings of the 20th AIAA Computational Fluid Dynamics Conference, Honolulu, AIAA Paper 2011-3044.
[10] Baty, H. and Nishikawa, H. (2016) Hyperbolic Method for Magnetic Reconnection Process in Steady State Magnetohydrodynamics. Monthly Notices of the Royal Astronomical Society, 459, 624-637.
https://doi.org/10.1093/mnras/stw654
[11] Chorin, H. (1997) A Numerical Method for Solving Incompressible Viscous Flow Problems. Journal of Computational Physics, 135, 118-125.
https://doi.org/10.1006/jcph.1997.5716
[12] Ghia, U., Ghia, K.N. and Shin, C.T. (1982) High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method. Journal of Computational Physics, 48, 387-637.
https://doi.org/10.1016/0021-9991(82)90058-4

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.