Hyperbolic Method to Explore Multiplicity Flow Solutions in a Four-Sided Lid-Driven Cavity ()
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
. 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
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,
(1)
for the velocity flow
, 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 being the time) Navier-Stokes equation follows,
(2)
The incompressibility property must be also imposed via,
(3)
For an homogeneous viscosity, an equivalent conservative form of Equation (2) and Equation (3) is,
(4)
with
being the velocity gradient tensor, and where
is the identity rank-2 tensor. The last equation makes appear two flux tensors, the advective flux modified by the pressure effect
which can be also called the inviscid flux below, and the viscous stress tensor
. The full flux is thus
.
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,
(5)
and,
(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.
) 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,
(7)
where
is the relaxation parameter defined by
, with
a length-scale parameter. The value of
value is determined in order to achieve an optimal fast convergence [6] [7] [8]. In this way, we have
as expected in the steady-state.
In two dimensions, the system 5 - 7 can be expressed in the following compact form,
(8)
with
being a source term containing the
tensor, and where
and
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
. Then, the set of equations becomes,
(9)
(10)
(11)
In this one dimensional version,
is a column vector containing three components, the pressure P, the velocity u, and g the gradient of u. The corresponding flux
has three components that are,
,
, and
. Finally, the source term
has three components, that are zero for the two first ones and
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
.
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
,
. Our finite-volume cell-centered discretization is thus,
(12)
where
denotes the solution value at the node j,
is an interface flux to be defined (see below), and
is a dual cell.
Steady state solutions of the previous system (i.e. with
) 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
. Moreover, in our case, the method is a relaxation procedure towards a steady-state, with an efficiency determined by the choice of the
parameter (or the associated
parameter). Steady state solutions of the previous system can be thus obtained by using the simply pseudo-time explicit iteration,
(13)
where n is the iteration counter,
is the pseudo-time step, and
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,
(14)
where the subscripts L and R stand for the left- and right-hand sides of the cell interface situated at
respectively. The first term is computed from an average value of the two fluxes
and
. 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
for the inviscid and viscous terms respectively. As a consequence, For the inviscid flux,
is taken as the maximum absolute value of the two eigenvalues
and represents a maximum inviscid wave speed. For the diffusive flux,
is simply
, 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,
(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,
(16)
where
is the gradient of
. A simple and economic method is to use
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
(see Figure 1). Using the numerical scheme described above in two dimensions with a spatial domain
, the resulting solution for the flow field lines is plotted in Figure 1 for a Reynolds number
(or a corresponding viscosity
). Indeed, we use a normalization where
, and
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
, 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
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
of the
velocity flow component as a function of x for two Reynolds number values,
and
. Note that an optimal value of
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
. Negative values in the two lower corners indicate reversal of the vortices rotation compared to the rest of the domain.
Figure 2.
profile obtained for
, in two runs using
and
. 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
, 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
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
. 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
and
. Indeed, contrary to the monotonic decrease of the residual obtained for the two lowest
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
. 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
. The direction imposed on each flow plate is indicated.
Figure 4. Residual term taken on one velocity flow component (that is
) as a function of the iteration counter, for a run of a four-sided cavity flow using
(see previous figure).
Figure 5. Residual term taken on one velocity flow component (that is
) as a function of the iteration counter, for runs of a four-sided cavity flow using
.
(a) (b)
Figure 6. The two stable steady state flow solutions (stream lines) obtained in two runs of a four-sided cavity flow using
. 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
for the bifurcation. We precisely get
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
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
, which is in the range of values
previously deduced in the literature.