Numerical simulation of two-dimensional laminar mixed-convection in a lid-driven cavity using the mixed finite element consistent splitting scheme

The Authors

Jeff C.-F. Wong, Department of Mathematics, The Chinese University of Hong Kong, Hong Kong, People's Republic of China

Acknowledgements

The author would like to thank the reviewers for their careful reading of the manuscript and valuable suggestions. This study was partly initiated while the author was a post-doctoral fellow at the Department of Physics, the Chinese University of Hong Kong (CUHK). The author would like to acknowledge the encouragement and helpful suggestions from Prof. Jun Zou and Prof. Emily S.C. Ching at CUHK. He would like to gladly acknowledge many stimulating discussions with Dr Michael Chan. During his visit at the National Singapore University, the author benefited from the advice and discussion of Prof. J-L. Guermond, who encouraged him to go public with his interest in the aspect of the consistent splitting scheme. He would like to express his appreciation to Frank Ng for help with the computing facility. The computations were performed on the IBM RS/6000 SP system at CUHK. Last but not least, the author is grateful beyond all exponential orders to Anita Spence who read his earlier manuscript.

Abstract

Purpose – The purpose of this paper is to propose an efficient/robust numerical algorithm for solving the two-dimensional laminar mixed-convection in a lid-driven cavity using the mixed finite element (FE) technique.

Design/methodology/approach – A numerical algorithm was based on the so-called consistent splitting scheme, which improved the numerical accuracy of the primitive variables. In order to obtain a stable solution, two choices of mixed FEs, the Taylor-Hood and Crouzeix-Raviart types, were used. Two mesh layouts were considered; uniform and non-uniform.

Findings – To verify that the proposed scheme had a second-order accuracy, some numerical results are presented and compared with the known solution. The answer was confirmative. Numerically accurate solutions were obtained for a fixed Prandtl number, Pr=0.71, for a range of the Reynolds number, Re from 100 to 3,000, and for a range of the Richardson number, Ri from 0.001 to 100. The results from these calculations, using the mixed FE consistent splitting scheme, agreed with the existing ones.

Research limitations/implications – Further extensions of this work could include the influence of various choices of Reynolds numbers, Prandtl numbers and Richardson numbers, and the effect of aspect ratio.

Originality/value – The present work was the first to apply a mixed FE in association with the consistent splitting scheme to the mixed convection problem.

Article Type:

Research paper

Keyword(s):

Finite element analysis; Convection; Programming and algorithm theory.

Journal:

International Journal of Numerical Methods for Heat & Fluid Flow

Volume:

17

Number:

1

Year:

2007

pp:

46-93

Copyright ©

Emerald Group Publishing Limited

ISSN:

0961-5539

1 Introduction

A situation where both the forced and natural convection effects are of comparable order is known as mixed convection. Well-known examples of convection applications include food processing and float glass production (Pilkington, 1969), thermal-hydraulics of nuclear reactors (Ideriah, 1980), dynamics of lakes and reservoirs (Imberger and Hamblin, 1982), heat transfer in solar collectors (Cha and Jaluria, 1984) and thermal convection of micropolar fluids (Hsu et al., 1995). Most of the work already done in a square domain considers thermal boundary conditions applied uniformly along two horizontal (or vertical) fluid layers with uniform temperatures and two vertical (or horizontal) adiabatic walls, and velocity boundary conditions established by shear and induced by a mechanically-driven sliding lid. Inspired by the work of Torrance et al. (1972), the crux of the paper is to study the value of the Richardson number (Gr/Re ²=Grashof number/Reynolds number²), which provides a measure of the importance of buoyancy-driven natural convection relative to lid-driven forced convection. Two competing convection mechanisms are also determined by the choice of Reynolds number and Prandtl number. Under the framework of experimental validation and numerical verification, Prasad and Koseff (1989) investigated the recirculation of the flow patterns induced by combined forced and natural convection heat transfer in a deep lid-driven cavity filled with water. Their results indicate that the overall heat transfer rate is a very weak function of the Grashof number for the examined range of the Reynolds number.

Numerical approximation of mixed convection problems in a lid-driven cavity offers a viable approach for analyzing the interaction of velocity and temperature variations. The following are some of the numerical studies on the mixed convection problems in a lid-driven cavity.

Moallemi and Jang (1992) used the Semi-Implicit Method for Pressure-Linked equations-revised (SIMPLER) algorithm with an inertial relaxation method to examine the effect of the Prandtl numbers on mixed convection problems in a lid-driven cavity. They considered the lower wall being heated and upper wall being cooled or vice-versa, which gives rise to mixed convective flows in a shallow lid-driven cavity, for a range of the Prandtl numbers. See also Amirthaganesh and Raghunandan (1998).

Iwatsu et al. (1993) extended the investigation of the mixed convection flows in a lid-driven cavity by including the unsteady term in the momentum and energy equations. They used an amended version of the Marker and Cell (MAC) method with the finite difference methods to study the effect of the Richardson numbers. Their results indicated that a decrease in Gr/Re ² produces flow features that are similar to those of a convectional driven cavity of non-stratified fluids. As Gr/Re ² increases, much of the fluid in both the middle and bottom positions of the cavity interior becomes stagnant. Later, Iwatsu and Hyun (1995) studied three-dimensional driven-cavity flows with a vertical temperature gradient. Tang and Tsang (1995) used the Least-Squares Finite Element Method (LSFEM) and Jacobi conjugate gradient technique for solving the transient mixed convection problem and provided comparisons with the earlier work of Iwatsu et al. (1993).

Buoyancy-aided and buoyancy-opposed cases of the mixed convection problem in a shear- and buoyancy-driven cavity were studied by Aydin (1999). His results indicated that the range of the Richardson numbers can be classified into three regimes; the pure forced, mixed and pure natural convection mechanisms. He employed a control volume-based finite difference method with the Alternative Direction Implicit (ADI) method for solving the steady vorticity transport and energy equation, and the Successive Over-Relaxation (SOR) method for solving the stream-function. Chamkha (2002), who did similar numerical works included the unsteady term in the system of governing equations.

Recently, Nicolás and Bermúdez (2005) used the stream-function/vorticity formulation for solving the transient mixed convection problem and provided comparisons with the earlier work of Iwatsu et al. (1993). Based on the operator splitting scheme, Bermúdez and Nicolás (1999) also considered the primitive variable formulation for solving a similar problem. Emphasis is placed on the choice of high Rayleigh and Reynolds numbers.

Despite the importance of combined buoyancy, thermal-driven and lid-driven convection, little prior work appears to have been done on the effect of pressure in a primitive variable formulation. The present investigation of combined buoyancy, thermal-driven and lid-driven convection in a square cavity is therefore undertaken to show the influence of pressure on a single moving-lid. It will be shown that, for some selected choices of Richardson numbers, the laminar mixed convection can have a significantly strong influence on the flow field, temperature and pressure in a cavity.

Our objective is to obtain accurate numerical solutions to the above stated problem for various choices of Gr/Re ², in order to characterize the nature of convective flows in a square cavity using the consistent splitting scheme. A scheme for solving the momentum and energy equations by using the mixed finite element (FE) methods is proposed. In this study, the pressure contributions to the boundary layer effect are reduced significantly, except at the corners. FE methods have been extensively developed for the mixed convection problems (Oosthuizen and Paul, 1985; Morzynski and Popiel, 1988) and other related problems and are now in wide use; to our knowledge, numerical results for detailed comparison by means of the choices of mixed FEs using the present approach are not yet readily available. Two classical types of mixed FEs that are known to be stable and of optimal convergence will be conducted in this study: Taylor-Hood and Crouzeix-Raviart elements with continuous pressure approximation.

The outline of the paper is as follows: in Section 2, we review the governing equations of the unsteady incompressible viscous flows and the energy equation. We provide a brief description of the consistent splitting scheme in Section 3, with more detail on the concept of numerical implementation and discretization. The numerical performances of the scheme and the computational test cases are illustrated and discussed in Section 4. In Section 5, we present our numerical examples, and end with concluding remarks in Section 6.

Parts of this work have been reported by the author Wong (2005) at the International Conference on Scientific Computing, Nanjing, China.

2 Governing equation

Laminar mixed convection in a lid-driven cavity is governed by unsteady Navier-Stokes (NS) equations and an energy equation with the Boussinesq approximation for buoyancy. Let Ω be a bounded domain in Euclidean space R² with a piecewise smooth boundary ∂Ω. A fixed final time is denoted by t _f. Let us begin with the equation of motion as follows: Equation 1 with: Equation 2 where t is the time (s), x=(x, y) is a function of 2D position in the Cartesian coördinate system; u=u(x, t)=u(x, t)i+v(x, t)j, ρ, f=f(x, t) are the velocity vector, density, the body force vector acting on the system, respectively; ∇ is the gradient operator and ∇ ·  is divergence operator. The spatial and temporal domain of a solution is denoted by Ω×(0, t _f). For the Newtonian fluids, the constitutive relationship for the stress tensor σ is given by Newton's (or the NS) law as: Equation 3 Equation 4 where I is a unit tensor, τ is the deviator stress tensor, ɛ(u)=(1/2)(∇u+(∇u)^{T^}) is the rate-of-strain tensor, p=p(x, t) is the pressure, μ is the coefficient of the dynamic viscosity, and λ is the second coefficient of viscosity; T^ stands for the transpose. Substituting equations (2.3) and (2.4) into equation (2.1) yields: Equation 5 where ∇² is the Laplacian operator. If the fluid is incompressible and homogeneous, then both μ and λ are constants, and the fluid properties are assumed to be constant except the density in the thermal buoyancy term of the governing equations. Equation (2.5) can be further simplified as: Equation 6 where ν=μ/ρ is called a kinematic viscosity; g is the acceleration due to gravity, β is the volumetric coefficient of thermal expansion, T is the temperature and T _r is a reference temperature. By taking into consideration the Boussinesq approximation for buoyancy, we shall be concerned with the two-dimensional (2D) unsteady (transient or time-dependent) NS equation of incompressible viscous fluid in a dimensional form, cf. equations (2.6) and (2.2).

Dirichlet boundary conditions on the velocity vector are: Equation 7 with[1]: Equation 8 where b is a known velocity vector function on the entire boundary ∂Ω of Ω (or ∂Ω being the piecewise smooth boundary of Ω), and n is the outward-pointing normal vector on ∂Ω. Initial conditions are: Equation 9 with[2]: Equation 10 and: Equation 11 where u ₀ is a prescribed velocity vector function and Ω=Ω∪∂Ω will denote the closure of Ω. Note that if any of the three constraints on the input data, equations (2.8), (2.10) and (2.11), are violated, the problem is ill-posed and no solution exists (Gresho et al., 1998). In addition, the data are assumed to be regular enough and to satisfy all the compatibility conditions needed for a smooth solution to exist for all time (Quartapelle, 1994).

Let us mention the crucial idea that the pressure is uniquely determined up to an additive constant. Problems of the form: Equation 12 with p(x, t) + Const can be considered as a variety of pressure solutions to the NS equations. It is noted that the mathematical expression p(x, t) + Const can be traced back to the Ladyžhenskaya (1963, p. 47) book. It is worth mentioning that an arbitrary constant does not depend on spatial variables. To see this, take the gradient operator acting on the pressure scalar field, so we have: Equation 13 Note that the constant term ∇Const goes to zero. In the Ladyžhenskaya (1975) review article, an arbitrary constant may depend on temporal variables. We adopt her notations and define the pressure scalar field such that: Equation 14 Equation (2.14) states that the pressure plus any arbitrary constant varying with time is also a solution of the pressure in the NS system[3].

The governing equation of energy for unsteady flow is given by the following: Equation 15 where α is the thermal diffusivity given by: Equation 16 where k is thermal conductivity and C _p is the specific heat capacity. The entire boundary ∂Ω can be subdivided into two disjoint sets, ∂Ω₁ and ∂Ω₂. Thus, ∂Ω admits the following decomposition ∂Ω=∂Ω₁∪∂Ω₂ and =∂Ω₁∩∂Ω₂. Neumann and Dirichlet boundary conditions on the temperature are: Equation 17 and: Equation 18 where q is a specified normal heat flux on ∂Ω₁ and T _b is a known temperature scalar function on ∂Ω₂. Initial conditions are: Equation 19 where T ₀ is a prescribed temperature scalar function on Ω.

In the above equations the variables are reduced to dimensionless form by the introduction of the following scales: Equation 20 The symbol * denotes dimensional variables. The notations T _hot and T _cold denote the temperatures of hot and cold walls, respectively. The one side of the cavity wall with constant velocity for all times is U ₀. The variables u* and v* are the velocity components in the x* and y* directions, p* is the pressure, T* is the temperature and t* is the time.

With the introduction of the dimensionless form, a system of non-dimensional equations in the primitive variable formulation is obtained. In conservation form, this is written as:

Momentum equation: Equation 21 Incompressibility constraint: Equation 22 Energy equation: Equation 23 where we have denoted the dimensionless variables by the same symbols as the corresponding dimensional ones, and the term j that represents the buoyancy force effect on the flow field has ± signs; the plus sign indicates the buoyancy-upward (or buoyancy-assisted) flow while the negative stands for buoyancy-downward (or buoyancy-opposed) flow. The former case will be used for all the numerical calculations. The four governing parameters are: Equation 24 Equation 25 Equation 26 Equation 27 The cavity height and width is denoted by H and W(=H). The aspect of the cavity is A=H/W=1. The Grashof number can be defined as the ratio of the Rayleigh number to the Prandtl number.

The dimensionless initial and boundary conditions that correspond to Figure 1 are:

• The square cavity is initially produced by motionless fluid, the initial pressure is zero and the temperature is zero throughout the entire domain Ω=[0, 1]×[0, 1].
• The left- and right-hand walls are insulated (or adiabatic) and: Equation 28

We select T _cold=0, T _hot=1 and U ₀=1. For an adiabatic wall, no heat transfer is permitted through the lower and upper boundaries.

The total kinetic energy of the flow is defined by: Equation 29 which indicates important characteristic of the mixed convection.

3 Consistent splitting scheme

One of the easiest ways for solving the unsteady NS equation is the fractional-step projection method (or the operator-splitting method), on which the present algorithm is partially based, that was independently introduced by Chorin (1968, 1969) and Témam (1969a, b) in the late 1960s. The method is based on the calculation of an intermediate velocity from the momentum equation where the pressure gradients are neglected; under the umbrella of the Helmholtz decomposition theorem, the pressure link between the incompressibility constraint and the momentum equations is established by transforming the incompressibility constraint into a Poisson equation for pressure, with particular emphasis on the homogeneous Neumann boundary condition. Finally the velocity is corrected using the computed pressure.

Based on the operator-splitting technique for solving the unsteady-state Stokes and hyperbolic systems, Timmermans et al. (1996) have investigated the palliations of the artificial Neumann boundary condition for the pressure, in which the correction term for the pressure should be added. Later, the same investigators extended their numerical methods to the thermally-driven flow problems in Minev et al. (1996). By making use of Goda's (1979) pressure-correction technique, high-order accuracies for the pressure can be achieved. More significantly, the pressure approximation is no longer spoiled by a numerical boundary layer produced by an artificial Neumann boundary condition (equation (3.36)). The mathematical proof of stability analysis on this development was put forth by Guermond and Shen (2003). The latest developments on the consistent splitting scheme can be found in Guermond et al. (2006).

The main ideas of the consistent splitting scheme can be summarized as follows:

• multiplying equation (2.21) by a test function[4] ∇q, ∀q∈H ¹ and isolating the pressure gradient as an unknown variable, and then obtaining the weak form of the Poisson equation for the pressure;
• using the div-curl identity ∇² u=∇(∇ · u)−∇×∇×u=−∇×∇×u for replacing ∇² u;
• applying equation (2.21) for simplifying the non-linear inertial term and source term, and putting the viscous term and time derivative term into the final form;
• using the gradient operator property for grouping the divergence-free term times the Reynolds number by selecting the arbitrary constant 1/Re,, i.e. (1/Re)(∇² u−∇×∇×u)=(1/Re)∇(∇ · u);
• introducing the auxiliary pressure variable φ, and extracting the pressure correction term plus the new term (1/Re)∇ · u;
• solving φ using equation (3.35) first and then the pressure by equation (3.36); and
• keeping the fractional-step projection methods in place, i.e. equations (3.34) and (3.35).

Three main characteristics are present in the scheme. One is the self-maintenance concept of the usual fractional-step projection method, such that the incompressibility constraint is directly satisfied by solving the Poisson problem for the pressure. Instead of determining the pressure, the auxiliary pressure is calculated. Another is the diminishment of the boundary-layer effect for the pressure, such that an extra divergence-free term derived from the viscosity term is included in the pressure Poisson equation in a weak sense. The last is that there is no need to determine an intermediate velocity.

A second-order decoupled approximation to the non-dimensional NS and energy equations in the mixed convection problem is defined as follows: let u ⁰=u ⁻¹=u(x, 0), p ⁰=p ⁻¹=p(x, 0), T ⁰=T ⁻¹=T(x, 0). Let (u ⁿ , p ⁿ , and T ⁿ ) be the nth-time step to (u(x, nΔt), p(x, nΔt), and T(x, nΔt)). For n=1, find u ¹, p ¹ and T ¹ such that: Equation 30 Equation 31 Equation 32 Equation 33 Then, for n≥1, find u ⁿ⁺¹, p ⁿ⁺¹ and T ⁿ⁺¹ such that: Equation 34 Equation 35 Equation 36 Equation 37 The boundary conditions on the velocity u ⁿ⁺¹|_∂Ω=b ⁿ⁺¹ are essential; the boundary conditions on the temperature T ⁿ⁺¹|_∂Ω1 =T _hot or T _cold are essential and (∂T ⁿ⁺¹/∂n)|_∂Ω2 =0 are natural, i.e. no-permeability.

3.1 Stability results

Since, the present stability analysis closely follows that of Guermond and Shen (2003), it will only be grouped the energy equation into the NS system in the consistent splitting scheme framework. Let us consider the system of equations (3.34)-(3.37): Equation 38 Equation 39 Equation 40 Equation 41 where g ⁿ⁺¹=(Gr/Re ²)T ⁿ j−(u ⁿ  · ∇)u ⁿ⁺¹−(1/2)(∇ · u ⁿ )u ⁿ⁺¹ and h ⁿ⁺¹=−(u ⁿ⁺¹ · ∇)T ⁿ⁺¹. For simplicity, we restrict ourselves to considering the homogeneous Dirichlet data u ⁿ⁺¹|_∂Ω=0 and T ⁿ⁺¹|_∂Ω=0. The backward difference operator δ can be defined as follows: δ u ⁿ⁺¹=u ⁿ⁺¹−u ⁿ , δ u ⁿ =u ⁿ −u ⁿ⁻¹, and δ ² u ⁿ⁺¹=δ u ⁿ⁺¹−δ u ⁿ , for instance. What we need to provide is an unconditionally bound estimation of the system (equations (3.38)-(3.41)), where the system only depends on the initial conditions u ⁰, u ¹, T ⁰ and T ¹.

Lemma 3.1

There exists a positive constant C such that the solution of equations (3.38)-(3.41) is bounded in the following sense: Equation 42

Proof

Taking the difference between two incremental time-steps on equation (3.38), applying the operator δ to corresponding equations, using the div-curl identity −∇² u=∇(∇ · u)−∇×∇×u, and grouping and replacing p ⁿ −p ⁿ⁻¹+(1/Re)∇ · u ⁿ by φ ⁿ , we have: Equation 43 Multiplying 2Δtδ u ⁿ⁺¹ by equation (3.43), integrating the corresponding equation over Ω taking inner product of equation (3.43) and using integration by parts w.r.t. the second and the third terms, we have: Equation 44 Using the identity (a−b, 2a)=|a|²−|b|²+|a−b|², we express the first and the third terms as: Equation 45 and: Equation 46 To introduce ∇φ ⁿ⁺¹, we use the operator δ in equation (3.39), i.e.: Equation 47 Substituting q=2Δt ² φ ⁿ into equation (3.47), and using the identity (b−a, 2a)=|a|²−|b|²−|a−b|², we have: Equation 48 Putting q=2Δt ² φ ⁿ⁺¹ into equation (3.47), and changing the index n+1 to n, we obtain: Equation 49 Likewise, putting q=Δt ² δφ ⁿ⁺¹ into equation (3.47), using the Cauchy-Schwarz inequality on the right-hand side, multiplying by 1/‖∇δφ ⁿ⁺¹‖² and taking the square of both sides, we get: Equation 50 Taking the difference between two consecutive times on equation (3.41), and applying the operator δ to the corresponding equations, we have: Equation 51 Multiplying 2ΔtδT ⁿ⁺¹ by equation (3.51), integrating the corresponding equation over Ω taking the inner product of the resultant equations, and using the above algebraic identity, we have: Equation 52 Likewise, multiplying 2Δt(δT ⁿ⁺¹−δT ⁿ ) by equation (3.51), we have: Equation 53 Summing up equations (3.44), (3.45), (3.46), (3.48), (3.49), (3.50), (3.52) and (3.53), and using the Cauchy-Schwarz and Poincaré inequalities, and the results ‖δ u ⁿ⁺¹‖≤c ₁‖∇δ u ⁿ⁺¹‖ ‖δT ⁿ⁺¹‖≤c ₂‖∇δT ⁿ⁺¹‖ and ‖δ ² T ⁿ ‖≤c ₃‖∇δ ² T ⁿ⁺¹‖ and using the arithmetic-geometric mean inequality, we have: Equation 54 Using the identity ‖∇u‖²=‖∇×u‖²+‖∇ · u‖², which holds for all u∈H ₀ ¹(Ω), and taking the sum of equation (3.54) from n=1 to k, we finally get: Equation 55 where c _*=4c ₂ ²−(1/2)c ₃ ²≥0. This completes the proof.

3.2 Temporal discretization

A fully implicit second-order backward differentiation formula (BDF), that is of accuracy O(Δt ²), is used for the time derivative that is of second-order accuracy, i.e.: Equation 56 where w replaces the velocity vector field u and the temperature scale field T. This scheme is usually recommended for the time discretization of the time derivative term, because the scheme is stable and second-order accurate.

3.3 Treatment of non-linear term

The treatment of the non-linear advection term transformed into a linear advection one using a linear extrapolation in time was discussed by Turek (1996). The non-linear term (u ⁿ⁺¹ · ∇)u ⁿ⁺¹ is replaced in a semi-implicit way by (u ⁿ  · ∇)u ⁿ⁺¹ or ((2u ⁿ −u ⁿ⁻¹) · ∇)u ⁿ⁺¹. The two corresponding terms have different numerical accuracies in each time step. The former is of first-order only, while the latter is of second-order. As suggested by Témam (1969a, b) and considered by many others (Quarteroni and Valli, 1994, p. 444), adding the kinetic term that is regarded somewhat as the skew-symmetric form can have no restriction on the choice of time step Δt, which guarantees an unconditional stability. To enhance the accuracy, we use (1/2)(∇ · (2u ⁿ −u ⁿ⁻¹))u ⁿ⁺¹. More importantly, when the approximate velocity fields do not exactly satisfy the incompressibility condition, this term plays an essential role in preserving the dissipativity of the discrete system (Shen, 1992).

3.4 Treatment of pressure gradient and temperature in the NS equation

One of the most important features of the consistent splitting scheme in the NS equation (3.34) is the linear extrapolation in time of the pressure gradient term ∇(2p ⁿ −p ⁿ⁻¹) and the temperature term 2T ⁿ −T ⁿ⁻¹, for n≥1 in order to enhance the accuracy of the second-order rate.

3.5 Homogenous Neumann boundary condition

When working on the homogenous Neumann boundary condition, one has to pay special attention; it is well known that the solution of equation (3.35) is known up to an additional constant when subjected to the homogeneous Neumann boundary condition (Kellogg, 1929, Chapter xi). As noted in literature, there are two common ways of obtaining a pressure scalar field solution unique to the NS system. It can be done by:

1. Option I. Imposing an additional condition on φ(x, t; Const(t)) such as specifying one point φ(x ₀, y ₀, t; Const(t)) on the boundary ∂Ω (a point of interest (x ₀, y ₀) must be adjusted).
2. Option II. Requiring that the integral of φ(x, t; Const(t)) over Ω is equal to zero or can be interpreted as the average value, i.e.: Equation 57

In the course of the numerical calculation for solving a singular (elliptic) Poisson problem, one may run into these two options. For the study of the first option in hydrodynamic stability, see, e.g. Ladyžhenskaya (1963, Chapter 1, p. 24) and especially the surveys by Rosenhead (1963, p. 55) and Chu (1978, p. 152). To remedy this difficulty more handily in terms of handling the homogeneous Neumann boundary condition for pressure, a fixed Dirichlet pressure boundary value on the boundary that does not vary with time is commonly found in applications of flow physics problems such as the 2D lid-driven cavity (Gresho et al., 1980, p. 33; Donea et al., 1981, p. 121; Comini and Del Giudice, 1982, p. 468). Assigning a fixed value is also directly suited for the numerical verifications of some existing numerical schemes (Timmermans et al., 1996). For further discussion, the reader is referred to Wong and Chan (2006). In this paper, for all numerical calculations, we shall use the augmented system (Henshaw and Petersson, 2003) for handling the uniqueness of pressure.

3.6 Spatial discretization

The FE formulations based on the Galerkin weighted residual approach are used for the discretization of the governing equations. Velocities, auxiliary pressures, pressures, and temperatures are taken as the unknown variables. The FE of choice are the Taylor-Hood and the Crouzeix-Raviart types, {P₂,P₁,P₁,P₂} and {P₂ ⁺,P₁,P₁,P₂ ⁺} elements, respectively, and their respective shape functions are used as shown in Figure 2; h stands for the maximum length of element measured from the triangulation. These choices verify the inf-sup condition which insures the uniqueness and existence of the solution of the discrete Stokes problems. For a comprehensive discussion, we refer the reader to Girault and Raviart (1981) or Brezzi and Fortin (1991). Use of the Crouzeix-Raviart FE types with discontinuous pressure fields for solving the incompressible flow and heat flow problems was carried out by Cuvelier et al. (1986), while the choice of continuous pressure fields, for instance, was employed by Soulaïmani et al. (1987).

In order to solve the corresponding equations by FE techniques, a variational form of the governing equations in a suitable function space is constructed, the domain is discretised into FEs, and local element matrices are calculated and assembled to form global matrices. The methods for obtaining the variational form of the governing equations are well-known and will not be addressed here (Gresho et al., 1998). The matrix forms of equations (3.34)-(3.36) are:

1. Given approximate velocities, {U _h ⁿ ,U _h ⁿ⁻¹,V _h ⁿ ,V _h ⁿ⁻¹}, corresponding pressures {P _h ⁿ ,P _h ⁿ⁻¹} and corresponding temperature {T _h ⁿ ,T _h ⁿ⁻¹}, solve: Equation 58 and: Equation 59 for approximate velocities at time t ⁿ⁺¹.
2. Using approximate velocities, solve the auxiliary pressure Φ _h ⁿ⁺¹: Equation 60
3. Using approximate velocities, solve the transfer pressure Π _h ⁿ⁺¹: Equation 61
4. Update the pressure at time t ⁿ⁺¹ viz.: Equation 62
5. Using approximate velocities, solve the temperature: Equation 63
6. Repeat steps 1-5 until a final number of time steps is reached.

For the sake of clarity, U _h ⁿ⁺¹,V _h ⁿ⁺¹,Φ _h ⁿ⁺¹,Π _h ⁿ⁺¹,P _h ⁿ⁺¹ and T _h ⁿ⁺¹are the vectors containing unknowns at the (n+1)th time step and: Equation 64 represent the consistent mass matrix, the convective matrix w.r.t (with respect to) the u-comp. (component), the convective matrix w.r.t the v-comp., the kinetic matrix, the diffusion matrix, and the divergence matrix w.r.t the x-comp., the divergence matrix w.r.t the y-comp., respectively; S^,DX,DY and M^ denote the diffusion matrix, the divergence matrix w.r.t the x-comp., the divergence matrix w.r.t the y-comp. and the consistent mass matrix, respectively; C^ _{U h ⁿ⁺¹} and C^ _{V h ⁿ⁺¹} represent the convective matrix w.r.t the u-comp. and the convective matrix w.r.t the v-comp. in the energy equation, respectively.

3.6.1 Remarks

The following numerical implementations were made within an in-house code that was based on the above algorithm:

• To compute all the element integrals we use classical quadrature formulae of Gaussian type. The integration over the triangles is performed in terms of the Gaussian quadrature rules using a four-point formula for the P₁ interpolation and a seven-point formula for the P₂ and P₂ ⁺ interpolations.
• Except the convective and kinetic matrices, all the matrices are computed once and stored in sparse matrix formats.
• Operated with the SPARSKIT package borrowed from Saad, preconditioners for sparse GMRES iterative solvers derived from threshold-based ILUT factorizations were used. The following selective parameters were used for all performance tests unless otherwise stated: the number of fill-in elements per row is 50; calculation is terminated when the relative precision is below ɛ=10⁻⁸; and convergence of the iterative process was fixed by a specific number of iteration. We shall provide more details in Section 4.
• To match the same degree of freedom of {P _h ⁿ ,P _h ⁿ⁻¹} in step 1, in the case of the Taylor-Hood FE type, we take the average value between two functional values in step 4. For the choice of the Crouzeix-Raviart FE type, we take the average value between two functional values first in step 4 and then calculate the centroid value using the linear pressure.
• To join two different types of FE interpolations and match the same degree of freedom of {U _h ⁿ⁺¹,U _h ⁿ ,U _h ⁿ⁻¹} in step 2 and {V _h ⁿ⁺¹,V _h ⁿ ,V _h ⁿ⁻¹} in step 3, two rectangular divergence matrices, ^DX _^ and ^DY _^, are formed. The integration over the triangles is performed in terms of the Gaussian quadrature rules using a seven-point formula for a combination of the P₁ and P₂, and the P₁ and P₂ ⁺ interpolations.
• A unique solution to the matrix system (equation (3.60)) will be found by solving the augmented system (Henshaw and Petersson, 2003, p. 115) such that: Equation 65

where r={1, … ,1} is the right null vector of S^,r ^{T^} denotes the transpose of r and G ⁿ⁺¹ denotes the right-hand side of equation (3.60). The last equation will set the mean value of φ to zero.

4 Numerical results for the mixed convection problems

As a numerical verification for the order-accuracy in time and in space of the above schemes combined with the Reynolds number Re, Prandtl number Pr and Richardson number Gr/Re ², an analytical example used by Auteri and Parolini (2002), which belongs to a special sub-branch of the Beltrami flows (Ethier and Steinman, 1994), is considered. The exact solution of the problem equations (2.1) and (2.2) on the square domain Ω=[−1, 1]×[−1, 1] is chosen as follows: Equation 66 Equation 67 Equation 68 Equation 69 with the corresponding force terms: Equation 70 Equation 71 Equation 72 It shall be assumed that at time t=0, the velocity field, pressure and temperature are all zero, i.e. u(x, y, 0)=v(x, y, 0)=p(x, y, 0)=T(x, y, 0)=0. The horizontal velocity boundary conditions are: Equation 73 and the vertical velocity conditions are: Equation 74 Thus, the incompressibility constraint holds, i.e. n · u=0 on ∂Ω for t≥0.

The horizontal thermal boundary conditions are: Equation 75 and the vertical thermal conditions are: Equation 76 In order to reveal more hidden information both from the velocity and pressure boundary layers and from the corner vortices, a thin layer within the domain is made, which not only maintains the uniform mesh structure but also has a concentration of triangles along the boundaries. This so-called boundary-refined mesh was used in the Gervais et al. (1997) for studying a stability analysis of the 2D lid-driven cavity. The boundary-refined mesh is defined as the conformal image of the uniform mesh on [0, 1]² by the 2D mapping: Equation 77 where i denotes the number of nodes, and x(i), y(i) are sets of 2D uniform mesh cöordinates, and X(i), Y(i) are sets of 2D boundary-refined mesh cöordinates. The finer the mesh refinement, the thinner the layers are. For the purpose of illustration, the 2D boundary-refined Taylor-Hood FE meshes and the 2D uniform Crouzeix-Raviart FE meshes are shown in Figure 3. To verify the convergence rate with respect to the spatial discretization, we select the grid range from 129 × 129 to 513 × 513 for the Taylor-Hood element type, while we choose the uniform mesh layout for the Crouzeix-Raviart element type. Table I shows the convergence of the iterative process that was obtained in the optimal number of iterations.

By a series of numerical experiments we obtained the rate of convergence for Re=Gr=Pr=1. Figures 4 and 5 give, for different choices of the mesh size, the error norms against the time step size. Scales on both axes are logarithmic. The slopes of the lines allow an estimate of the rate of convergence.

In Figures 4 and 5, the observed rates of convergence are compared with the expected rates. Though the theoretical justifications of the present scheme are not developed, one may wonder that what is order of convergence? We shall measure all the unknown quantities in the L ²-norm sense. As for the ‖ · ‖ _{l ^∞(0,t f; L ²(Ω))}- and ‖ · ‖ _{l ²(0,t f; L ²(Ω))}-errors of the approximate velocity fields and the approximate temperature, the scheme appears to have expected convergence rates on O(Δt ²). As for the ‖ · ‖ _{l ^∞(0,t f; L ²(Ω))}- and ‖ · ‖ _{l ²(0,t f; L ²(Ω))}-errors of the approximate pressure, one infers from the figures, that the convergence rates are in O(Δt ²). As suggested, when the grid size is not fine enough, a reduction of the parameter Δt does not enhance the accuracy.

The error l ^∞(0, t _f; L ²(Ω)) and l ²(0, t _f; L ²(Ω))-norm calculations with respect to the FE mesh refinement of the spatial grids and the temporal refinement of the time-steps are shown in Figures 4 and 5. In all cases the convergent rate of the approximate primitive variables is of O(h ²+Δt ²). From these observations and other ones from various time-step sizes and mesh sizes, the scheme undoubtedly appears to be unconditionally stable.

5 Computational experiment

5.1 Calculation of stream-function

It is known that in the vorticity/stream-function formulation the stream-function ψ satisfies the Poisson equation: Equation 78 Here ω is the vorticity, and is defined by: Equation 79 The boundary condition for the stream-function is given by: Equation 80 which implies that there is no mass transfer passing through the wall of the cavity and that the given boundaries themselves represent one of the stream lines. We introduce the test function φ. Multiplying equation (5.77) by φ and using the Green-Gauss formula for the Laplacian operator, we obtain: Equation 81 Recall that the known boundary data ∂ψ/∂n=u _s , where s is the tangential direction, do not automatically vanish because the constant velocity U ₀ along the tangential direction is imposed on one of the walls. The above-mentioned formulations were used for all 2D graphical representation for the stream-function in this work. Unless otherwise stated, for the purpose of drawing figures, no attempt was made to smooth the obtained results, so the exact temperature and stream-function were used.

5.2 Solution procedure

In order to validate the precision and accuracy of the numerical scheme of the present study, a series of computations was run for aiding and opposing mechanisms of mixed convection in a shear- and buoyancy-driven cavity. When considering the mixed convection problems, we designate the flow as an aiding flow when the buoyancy force has a component in the direction of the shear velocity, and likewise, as an opposing flow when the buoyancy component is opposite to the shear velocity. The Reynolds number is kept fixed at Re=100. The Prandtl number of the air Pr=0.71 is used. Parametric studies of the effect of the mixed convection parameters Gr/Re ² that can be divided into three different heat-transfer regimes are reduced to the following three cases due to space restrictions:

1. Gr/Re ²=0.001 would be belong to a class of pure forced convection, i.e. Gr/Re ²≪1.
2. Gr/Re ²={0.1, 0.5, 1, 2, 5, 10} would lie within the between a range of mixed convection regimes, i.e. 0.1≤Gr/Re ²≤10.
3. Gr/Re ²=100 would fall into a class of pure natural convection, i.e. Gr/Re ²≫1.

For a more complete account of the classification of the three heat-transfer regimes, the reader is referred to Aydin (1999).

All the numerical results are completed at the dimensionless time t _f=90. The time step Δt=0.0005 was used for all the numerical testings. Boundary-refined and uniform spaced grid systems are used throughout the solution domain. Fine grids are used; 257 × 257 and 129 × 129 mesh layouts for the Taylor-Hood elements, and 257 × 257 + 32,768 and 129 × 129 mesh layouts for the Crouzeix-Raviart elements. The relative precision is set to 10⁻⁸; the number of iteration for the all the computations are fixed at 200.

5.2.1 Steady-state solution

Here, in order to gain some insight into a quasi-steady (or a nearly steady) situation, the time history of the total kinetic energy for certain cases, is shown in Figure 6. It can been seen that the increase of Gr/Re ² contributes to the large magnitude of the kinetic energy. The convergence criterion used to obtain the steady-state solution was the standard relative error, which is based on the L ²-norm and max-norm estimation given by: Equation 82 and: Equation 83 where ‖ · ‖₂ indicates the product of the variable over all the grid points in the computational domain and | · |_max refers to a maximum value over the entire field of grid point values. The author is fully aware of the fact that all the results presented here slowly vary with time. A longer time is needed before a steady-state solution is finally reached.

5.2.2 Incompressibility constraint

The time history of the divergence-free velocity field at the final dimensionless time t _f=90 is shown in Figure 7. As can be seen, the divergence-free velocity fields that are measured by the L ²-norm estimation have a constant rate.

5.2.3 Results and discussion

Figures 8-11 show the isotherms and stream lines for the aiding- and opposing-buoyancy flows for a range of the Richardson number, Ri from 0.001 to 100. Several interesting features are seen in Figures 8-11. When Gr/Re ²=0.001, the flow is strongly influenced by the force convection. However, as Gr/Re ² is set to 1, the flow displays a considerable change. Two convective cells, clearly seen in Figure 11(b), (d), (f) and (g) are the results of the interaction between the force convection and natural convection mechanisms. When Gr/Re ²=100, as shown in Figures 8(h) and 9(h), the movement of the flow reversal is obtained, resulting in a clockwise rotating cell (the buoyancy-aiding case) to a counter clockwise rotating one (the buoyancy-opposing case). The flow is mainly dominated by the natural convection. As can be seen, the results of the proposed numerical scheme are in good agreement with those done by Aydin (1999). This lends confidence in the results of the proposed scheme to be reported subsequently.

No attempt was considered to smooth out the corner singularities in this study. The Gibbs-like singularities of the pressure appeared at the sliding-lid velocity corners, as shown in Figure 12. Judging from the given figure, no boundary-layer effect for the pressure along the given boundaries was made. To demonstrate the negligible effect on the internal feature of the pressure, except in the neighborhood of the corner singularities, the following contour plots were cut down from the actual size of the computational grid, i.e. from 129 × 129 to 124 × 124. No node-to-node oscillations are made in Figure 12 by using the mixed FE approximation. Until otherwise stated, all the pressure contour plots followed the same procedure. As shown in Figures 13 and 14, the pressure is influenced by the cavity corners. When Gr/Re ² increases, the profile of the pressure distribution becomes more flat and uniform across the cavity, except the near corners. In the mixed convection region, the pressure distribution exhibit slight differences for the buoyancy-aiding and buoyancy-opposing cases.

Figures 15 and 16 show the velocity distributions at Ri={0.001, 0.01, 0.1, 0.5, 1, 2, 5, 10, 100}. The distributions of u along a vertical line and v along a horizontal line passing through the center of the cavity are shown in Figures 15 and 16. The results of eight cases with regard to the buoyancy-aiding and buoyancy-opposing cases are superimposed. Different Richardson numbers Gr/Re ² have different effects on the distributions near the walls. The shift of the heat transport mechanism from forced convection to natural convection has a strong effect on the velocity distribution when the Richardson number is large. Velocity profiles of the buoyancy-aiding and buoyancy-opposing cases mirror each other.

Profiles of the pressure distributions along the vertical line of the cavity at the mid-plane are shown in Figure 17 for various Richardson numbers Gr/Re ². We set a pressure datum point at the center of the cavity. It can be seen that both buoyancy-aiding and buoyancy-opposing cases have no major difference. For instance, no mirror profiles are observed. When Gr/Re ²=100, the magnitude of the pressure solution is significantly magnified, as compared with those of velocity solutions.

Representative profiles for the horizontal component of the temperature T, at the square cavity mid-plane perpendicular to the sliding lid for various Gr/Re ² are shown in Figure 18. As Gr/Re ² increases, the temperature distribution becomes visible and has greater effect near the wall.

5.2.4 Numerical performance of Crouzeix-Raviart FE

In this part, we investigate the computational performance of two classical types of mixed FEs. Computed pressure results using different mesh layouts in the mixed FE approximations are shown in Figures 8(g), 9(b), 13(g), 14(b), and 19. Judging from these figures, the computed results are almost identical to the characteristic of the pressure distribution, except the tiny difference near the wall corners. The present calculations are compared to the earlier ones with respect to the total kinetic energy time history as shown in Figure 20, which confirms the applicability of the consistent splitting scheme using the Crouzeix-Raviart element types. Overall, the comparison is quite satisfactory. But when solving larger problems by making use of the Crouzeix-Raviart elements, the cost of solving several global matrix systems may dominate.

5.2.5 More validation results

Our goal here is to validate the present numerical method against the works of Iwatsu et al. (1993) and Nicolás and Bermúdez (2005). The associated boundary conditions for the velocity and the temperature are given by: Equation 84 The Prandtl number is kept fixed at Pr=0.71. Two cases are considered:

1. The final time is set to t _f=200. We conclude that:
   • It can be seen from Figure 21(a) and (b) that the solution of the present numerical code agree with the numerical results of Iwatsu et al. (1993) for Gr=10² and Re=1,000. Good agreement was found between the minimum and maximum values of the velocities of the present results and the numerical results of Iwatsu et al. (1993) as shown in Table II.
   • Various choices of time-step sizes are used: Δt={0.005, 0.0025, 0.00125}. Judging from Table II, one infers that these results are consistent.
2. The final time is set to t _f=400. The Grashof number is Gr=10⁶ for Re=400, 1,000 and 3,000. The present numerical code was validated against the results of Iwatsu et al. (1993) and Nicolás and Bermúdez (2005). We conclude that:
   • It can be seen from the comparison that these solutions are in good agreement. We used different time-step sizes to check the results with the long time simulation. Table III confirms this agreement among the solutions. Owing to the space limitation, we only report the stream-function and temperature solutions for Gr=10⁶ and Re=1,000. As shown in Figure 21(c) and (d), these results agree with those produced by Iwatsu et al. (1993) and Nicolás and Bermúdez (2005).
   • Time histories of the total kinetic energy and the divergence-free velocity field at the final dimensionless time t _f=400 are shown in Figure 22. As can be seen, the increase of Re increases the magnitude of energy and the divergence-free velocity fields have a constant rate.
   • Figure 23 shows horizontal and vertical velocity profiles at the midsections of the cavity at X=0.5 and Y=0.5, respectively. When compared with the results of Iwatsu et al. (1993) and Nicolás and Bermúdez (2005), these results were found to be in good agreement.
   • Figure 24 shows velocity vector field and pressure scalar field plots for Gr=10⁶ and Δt=0.0025. When Re increases, the flow is dominated by the forced convection.

6 Conclusion

The purpose of this work was to propose a numerical algorithm that was based on the so-called consistent splitting scheme for solving the equations of the unsteady NS and energy equations so as to study the effect of Richardson numbers on the mixed convection in a 2D cavity. The right and left sides of the cavity are, respectively, heated and cooled isothermally and vice-versa, and the top and bottom sides are adiabatic. Both cases of the cavity with three rigid boundaries and a shear velocity lid have been considered. The convergence rate of the numerical scheme is of second-order in O(h ²+Δt ²) in the sense of the L ²-norm estimation. Two types of mixed FEs, Taylor-Hood and Crouzeix-Raviart, have been used and verified the expected order rate. One of the main features of the aiding- and opposing-buoyancy cases is that the pressure shows uniform patterns across the cavity, except near the corners as Gr/Re ² increases.

Further extensions of this work could include the influence of various choices of Reynolds numbers, Prandtl numbers and Richardson numbers, and the effect of aspect ratio; a study of the stability of this solution in order to find the threshold of interaction between the viscous dissipation and buoyancy driven convection would be of particular interest. These topics will be presented in a separate study.

Fortran 90 code that implements the presented algorithm is available over the internet from the author (contact: jwong@math.cuhk.edu.hk for details).

Equation 1

Equation 2

Equation 3

Equation 4

Equation 5

Equation 6

Equation 7

Equation 8

Equation 9

Equation 10

Equation 11

Equation 12

Equation 13

Equation 14

Equation 15

Equation 16

Equation 17

Equation 18

Equation 19

Equation 20

Equation 21

Equation 22

Equation 23

Equation 24

Equation 25

Equation 26

Equation 27

Equation 28

Equation 29

Equation 30

Equation 31

Equation 32

Equation 33

Equation 34

Equation 35

Equation 36

Equation 37

Equation 38

Equation 39

Equation 40

Equation 41

Equation 42

Equation 43

Equation 44

Equation 45

Equation 46

Equation 47

Equation 48

Equation 49

Equation 50

Equation 51

Equation 52

Equation 53

Equation 54

Equation 55

Equation 56

Equation 57

Equation 58

Equation 59

Equation 60

Equation 61

Equation 62

Equation 63

Equation 64

Equation 65

Equation 66

Equation 67

Equation 68

Equation 69

Equation 70

Equation 71

Equation 72

Equation 73

Equation 74

Equation 75

Equation 76

Equation 77

Equation 78

Equation 79

Equation 80

Equation 81

Equation 82

Equation 83

Equation 84

Equation 85

Equation 86

Equation 87

Figure 1Geometry of the 2D mixed convection in a lid-driven cavity problem

Figure 2{P₂,P₁,P₁,P₂} and {P₂ ⁺,P₁,P₁,P₂ ⁺} for {u _h , φ _h , p _h , T _h }

Figure 3Figure (a) and (b) are 2D boundary-refined meshes in two FE (P₁/P₂) spaces. Figure (c) and (d) are 2D uniform meshes in two FE (P₁/P₂ ⁺) spaces

Figure 4Convergence results for the velocity, pressure and temperature at Re=Gr=Pr=1 and t _f=1 using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 5Convergence results for the velocity, pressure and temperature at Re=Gr=Pr=1 and t _f=1 using the Crouzeix-Raviart element type with the uniform mesh layout

Figure 6Time history of the kinetic energy is obtained by the buoyancy-aiding case in the left panel and buoyancy-opposing case in the right panel

Figure 7Time history of the divergence-free term is obtained by the buoyancy-aiding case in the left panel and buoyancy-opposing case in the right panel

Figure 8Isotherms for the buoyancy-aiding solutions

Figure 9Isotherms for the buoyancy-opposing solutions

Figure 10Stream lines for the buoyancy-aiding solutions

Figure 11Stream lines for the buoyancy-opposing solutions

Figure 12Mesh plots of pressure

Figure 13Pressure scalar field plots for the buoyancy-aiding solutions

Figure 14Pressure scalar field plots for the buoyancy-opposing solutions

Figure 15Non-dimensional horizontal velocity profiles at Y=0.5, Re=100, Pr=0.71 and R _i ={0.001, 0.01, 0.1, 0.5, 1, 2, 5, 10, 100} using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 16Non-dimensional vertical velocity profiles at X=0.5, Re=100, Pr=0.71 and R _i ={0.001, 0.01, 0.1, 0.5, 1, 2, 5, 10, 100} using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 17Non-dimensional pressure profiles at X=0.5, Re=100, Pr=0.71 and Ri={0.001, 0.01, 0.1, 0.5, 1, 2, 5, 10, 100} using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 18Non-dimensional temperature profiles at Y=0.5, Re=100, Pr=0.71 and Ri={0.001, 0.01, 0.1, 0.5, 1, 2, 5, 10, 100} using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 19Temperature and pressure scalar field plots for the buoyancy-aiding solutions in the left panel and the buoyancy-opposing solutions in the right panel

Figure 20Time history of the kinetic energy: obtained by the buoyancy-aiding case in the left panel and using the buoyancy-opposing case in the right panel

Figure 21Stream-function and temperature scalar field plots for Gr=10² in the left panel and Gr=10⁶ in the right panel

Figure 22Time histories of the kinetic energy and the divergence-free term: Re={400, 1,000, 3,000}, Pr=0.71, Gr=10⁶ and Δt=0.0025 using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 23Non-dimensional vertical velocity profiles at X=0.5 and horizontal velocity profiles at Y=0.5, Re={400, 1,000, 3,000}, Pr=0.71, Gr=10⁶ and Δt=0.0025 using the Taylor-Hood element type with the boundary-refined mesh layout

Figure 24Velocity vector field and pressure scalar field plots for Gr=10⁶ and Δt=0.0025

Table I

Table II Gr=10²: (1) Re=1,000

Table III Gr=10⁶: (1) Re=400, (2) Re=1,000, (3) Re=3,000

References

Amirthaganesh, S., Raghunandan, B.N. (1998), "Mixed convection in a square cavity – some new results", Computational Fluid Dynamics Journal, Vol. 6 No.4, pp.43-51.

[Manual request] [Infotrieve]

Auteri, F., Parolini, N. (2002), "Numerical investigation of the first instabilities in the differentially heated 8:1 cavity", Int. J. Numer. Methods Fluids, Vol. 40 pp.1121-32.

[Manual request] [Infotrieve]

Aydin, O. (1999), "Aiding and opposing mechanisms of the mixed convection problems in a shear- and buoyancy-driven cavity", Int. Commun. Heat Mass Transfer, Vol. 26 No.7, pp.1019-28.

[Manual request] [Infotrieve]

Bermúdez, B., Nicolás, A. (1999), "An operator splitting numerical scheme for thermal/isothermal incompressible viscous flows", Int. J. Numer. Methods Fluids, Vol. 29 pp.397-410.

[Manual request] [Infotrieve]

Brezzi, F., Fortin, M. (1991), Mixed and Hybrid Finite Element Methods, Springer-Verlag, Berlin, .

[Manual request] [Infotrieve]

Cha, C.K., Jaluria, Y. (1984), "Recirculating mixed convection flow for energy extraction", Int. J. Heat Mass Transfer, Vol. 27 No.10, pp.1801-12.

[Manual request] [Infotrieve]

Chamkha, A.J. (2002), "Hydromagnetic combined convection flow in a vertical lid-driven cavity with internal heat generation or absorption", Numer. Heat Transfer – Part A: Applications, Vol. 41 pp.529-46.

[Manual request] [Infotrieve]

Chorin, A.J. (1968), "Numerical solution of the Navier-Stokes equations", Math. Comput., Vol. 22 pp.745-62.

[Manual request] [Infotrieve]

Chorin, A.J. (1969), "On the convergence of discrete approximations to the Navier-Stokes equations", Math. Comput., Vol. 23 pp.341-53.

[Manual request] [Infotrieve]

Chu, C.K. (1978), "Numerical methods in fluid dynamics", Advances in Applied Mechanics, Vol. 18 pp.285-331.

[Manual request] [Infotrieve]

Comini, G., Del Giudice, S. (1982), "Finite-element solution of the incompressible Navier-Stokes equations", Numer. Heat Transfer, Vol. 5 pp.463-78.

[Manual request] [Infotrieve]

Cuvelier, C., Segal, A., Van Steenhoven, A.A. (1986), Finite Element Methods and Navier-Stokes Equations, D. Reidel Publishing, Dordrecht, .

[Manual request] [Infotrieve]

Donea, J., Giuliani, S., Laval, H., Quartapelle, L. (1981), "Solution of the unsteady Navier-Stokes equations by a finite-element projection method", in Taylor, C., Morgan, K. (Eds),Recent Advances in Numerical methods in Fluids, Prineridge, Swansea, Vol. Vol. 2 pp.97-132.

[Manual request] [Infotrieve]

Ethier, C.R., Steinman, D.A. (1994), "Exact fully 3D Navier-Stokes solutions for benchmarking", Int. J. Numer. Methods Fluids, Vol. 19 pp.369-75.

[Manual request] [Infotrieve]

Gervais, J.J., Lemelin, D., Pierre, R. (1997), "Some experiments with stability analysis of discrete incompressible flows in the lid-driven cavity", Int. J. Numer. Methods Fluids, Vol. 24 pp.477-92.

[Manual request] [Infotrieve]

Girault, V., Raviart, P. (1981), Finite Element Approximation of the Navier-Stokes Equations, Springer-Verlag, New York, NY, .

[Manual request] [Infotrieve]

Goda, K. (1979), "A multistep technique for the with implicit difference schemes for calculating two- or three-dimensional cavity flows", J. Comput. Physics, Vol. 30 pp.76-95.

[Manual request] [Infotrieve]

Gresho, P., Lee, R.L., Sani, R. (1980), "On the time-dependent solution of the incompressible Navier-Stokes equations in two and three dimensions", in Taylor, C., Morgan, K. (Eds),Recent Advances in Numerical methods in Fluids, Prineridge, Swansea, Vol. Vol. 1 pp.27-79.

[Manual request] [Infotrieve]

Gresho, P., Sani, R., Engelman, M.S. (1998), Incompressible Flow and the Finite Element Method Advection-Diffusion and Isothermal Laminar Flow, Wiley, Chichester, .

[Manual request] [Infotrieve]

Guermond, J.L., Shen, J. (2003), "A new class of truly consistent splitting schemes for incompressible flows", J. Comp. Physics, Vol. 192 pp.262-76.

[Manual request] [Infotrieve]

Guermond, J.L., Minev, P., Shen, J. (2006), "An overview of projection methods for incompressible flows", Comput. Methods Appl. Mech. Eng., Vol. 195 pp.6011-45.

[Manual request] [Infotrieve]

Henshaw, W.D., Petersson, N.A. (2003), "A split-step scheme for the incompressible Navier-Stokes equations", in Hafez, M.M. (Eds),Numerical Solutions of Incompressible Flows, World Scientific Publishing Co., Singapore, pp.108-25.

[Manual request] [Infotrieve]

Hsu, T.H., Hsu, P.T., Chen, C.K. (1995), "Thermal convection of micropolar fluids in a lid-driven cavity", Int. Commun. Heat Mass Transfer, Vol. 22 No.2, pp.189-200.

[Manual request] [Infotrieve]

Ideriah, F.J. (1980), "Prediction of turbulent cavity flow driven by buoyancy and shear", J. Mech. Engng Sci., Vol. 22 No.6, pp.287-95.

[Manual request] [Infotrieve]

Imberger, J., Hamblin, P.F. (1982), "Dynamics of lakes, reservoirs, and cooling ponds", Ann. Rev. Fluid Mech., Vol. 14 pp.153-87.

[Manual request] [Infotrieve]

Iwatsu, R., Hyun, J.M. (1995), "Three-dimensional driven-cavity flows with a vertical temperature gradient", Int. J. Heat Mass Transfer, Vol. 38 No.18, pp.3319-28.

[Manual request] [Infotrieve]

Iwatsu, R., Hyun, J.M., Kuwahara, K. (1993), "Mixed convection in a driven cavity with a stable vertical temperature-gradient", Int. J. Heat Mass Transfer, Vol. 36 No.6, pp.1601-8.

[Manual request] [Infotrieve]

Kellogg, O.D. (1929), Foundations of Potential Theory, Springer, Berlin, .

[Manual request] [Infotrieve]

Ladyžhenskaya, O.A. (1963), Mathematical Problems in the Dynamics of a Viscous Incompressible Flow, Gordon & Breach, New York, NY, .

[Manual request] [Infotrieve]

Ladyžhenskaya, O.A. (1975), "Mathematical analysis of Navier-Stokes equations for incompressible liquids", Ann. Rev. Fluid Mech., Vol. 7 pp.249-72.

[Manual request] [Infotrieve]

Minev, P.D., Van De Vosse, F.N., Timmermans, L.J.P., Van Steenhoven, A.A. (1996), "A second order splitting algorithm for thermally-driven flow problems", Int. J. Num. Meth. Heat Fluid Flow, Vol. 6 No.2, pp.51-60.

[Manual request] [Infotrieve]

Moallemi, M.K., Jang, K.S. (1992), "Prandtl number effects on laminar mixed convection heat transfer in a lid-driven cavity", Int. J. Heat Mass Transfer, Vol. 35 No.8, pp.1881-92.

[Manual request] [Infotrieve]

Morzynski, M., Popiel, C.O. (1988), "Laminar heat transfer in a two-dimensional cavity covered by a moving wall", Numer. Heat Transfer, Vol. 13 pp.265-73.

[Manual request] [Infotrieve]

Nicolás, A., Bermúdez, B. (2005), "2D thermal/isothermal incompressible viscous flows", Int. J. Numer. Methods Fluids, Vol. 48 pp.349-66.

[Manual request] [Infotrieve]

Oosthuizen, P.H., Paul, J.T. (1985), "Mixed convective heat transfer in a cavity", Fundamental Forced and Mixed Convection, Proceedings of the 23rd National Heat Transfer Conference, August 4-7, Denver, Colorado, HTD-vol., pp.159-69.

[Manual request] [Infotrieve]

Pilkington, L.A.B. (1969), "Review lecture: the float glass process", Proc. Roy. Soc. Lond. A., Vol. 314 pp.1-25.

[Manual request] [Infotrieve]

Prasad, A.K., Koseff, J.R. (1989), "Combined forced and natural convection heat transfer in a deep lid-driven cavity flow", Heat Transfer in Convective Flows, ASME, New York, NY, HTD-107, pp.155-62.

[Manual request] [Infotrieve]

Quartapelle, L. (1994), Numerical Solutions of the Incompressible Navier-Stokes Equations, Birkhäuser, Basel, .

[Manual request] [Infotrieve]

Quarteroni, A., Valli, A. (1994), Numerical Approximation of Partial Differential Equations, Springer, Berlin, Springer Series in Computational Mathematics, Vol. Vol. 23.

[Manual request] [Infotrieve]

Rosenhead, L. (1963), "Introduction – boundary layer theory", in Rosenhead, L. (Eds),Laminar Boundary Layers: An Account of the Development, Structure, and Stability of Laminar Boundary Layers in Incompressible Fluids, Clarendon Press, Oxford, pp.46-109.

[Manual request] [Infotrieve]

Shen, J. (1992), "On error estimates of some higher order projection and penalty-projection methods for Navier-Stokes equations", Numer. Math., Vol. 62 No.1, pp.49-73.

[Manual request] [Infotrieve]

Soulaïmani, A., Fortin, M., Ouellet, Y., Dhatt, G., Bertrand, F. (1987), "Simple continuous pressure elements for two- and three-dimensional incompressible flows", Comput. Methods Appl. Mech. Eng., Vol. 62 pp.47-69.

[Manual request] [Infotrieve]

Tang, L.Q., Tsang, T.H. (1995), "Transient solutions by a least-squares finite-element method and Jacobi-conjugate gradient technique", Numer. Heat Transfer, Part B: Fundamentals, Vol. 28 pp.183-98.

[Manual request] [Infotrieve]

Témam, R. (1969a), "Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires I", Arch. Rational Mech. Anal., (in French), Vol. 32 pp.135-53.

[Manual request] [Infotrieve]

Témam, R. (1969b), "Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II", Arch. Rational Mech. Anal., (in French), Vol. 33 pp.377-85.

[Manual request] [Infotrieve]

Timmermans, L.J.P., Minev, P.D., Van De Vosse, F.N. (1996), "An approximate projection scheme for incompressible flow using spectral elements", Int. J. Numer. Methods Fluids, Vol. 22 pp.673-88.

[Manual request] [Infotrieve]

Torrance, K., Davis, R.W., Eike, K., Gill, P., Gutman, D., Hsui, A., Lyons, S., Zein, H. (1972), "Cavity flows driven by buoyancy and shear", J. Fluid Mech., Vol. 51 pp.221-31.

[Manual request] [Infotrieve]

Turek, S. (1996), "A comparative study of time-stepping techniques for the incompressible Navier-Stokes equations: from fully implicit non-linear schemes to semi-implicit projection methods", Int. J. Numer. Methods Fluids, Vol. 22 pp.987-1011.

[Manual request] [Infotrieve]

Wong, J.C.F. (2005), "Numerical simulation of two-dimensional laminar mixed-convection in a lid-driven cavity using the mixed finite element consistent splitting scheme: influence of Richardson numbers", paper presented at the 20 min contributed talk in International Conference on Scientific Computing (ICSC05), Nanjing, June 4-8, .

[Manual request] [Infotrieve]

Wong, J.C.F., Chan, M.K.H. (2006), "A consistent splitting scheme for unsteady incompressible viscous flows I. Dirichlet boundary condition and applications", Int. J. Numer. Methods Fluids, Vol. 51 pp.385-424.

[Manual request] [Infotrieve]

Further Reading

Saad, Y. (2005), "SPARSKIT: a basic tool-kit for sparse matrix computations (Version 2)", available at: www-users.cs.umn.edu/ ∼ saad/software/SPARSKIT/sparskit.html, .

[Manual request] [Infotrieve]

A Corresponding author

Jeff C.-F. Wong can be contacted at: jwong@math.cuhk.edu.hk