Analysis and visualization tools in CFD, part II: A case study in grid adaptivity

FINITE ELEMENTS IN ANALYSIS A N D DESIGN

ELSEVIER

Finite Elements in Analysis and Design 19 (1995) 309-324

Analysis and visualization tools in CFD, Part II: A case study in grid adaptivity Francois Guibault "'b, Andr6 G a r o n a'b, Benoit Ozell a'b'*, Ricardo Camarero "'b a CERCA, 5160 Dbcarie, bureau 400, Montrbal, Canada H3X 2H9 b Ecole Polytechnique de Montrbal, C.P. 6079, Stn "Centre-Ville", Montrbal, Canada H3C 3A7

Abstract The objective of the present paper is to show, in the context of a configurable resolution and analysis environment, how grid adaptivity methodology can be made configurable and rapidly applied to various problems• In order to reach this objective, our error estimation and grid adaptation tools were configured to integrate Zienkiewicz and Zhu's local error estimator. The environment was then applied to the simulation of incompressible isothermal flow and natural convection problems.

1. Introduction

For engineers 1:o be able to use CFD both as an analysis and as a design tool, they need a methodology to assess and control the accuracy of computed solutions. This methodology should be, as far as possible, independent of the resolution technique used, and applicable to any of the mathematical and physical models involved in the description of the problem. Keeping this primary requirement in mind, resolution methodologies must then be looked at from the perspective of accuracy and error control. In this respect, finite element method (FEM) distinguishes itself, since it provides the necessary theoretical tools needed for the error analysis process. Indeed, for every class of equation or set of equations, be it elliptic, parabolic or hyperbolic, one can devise adequate tools to perform error analysis on those equations. These tools, which estimate the error of the numerical solution, and relate this error to the local grid size, are, in fact, generic in two ways: first, one error estimator can be applied to a wide variety of equations, and second, for one set of equations, numerous error estimators can be used. Furthermore, the overall efficiency of one estimator, when applied to a specific set of equations, can

* Corresponding address: CERCA, 5160 D6carie, bureau 400, Montr6al, Canada H3X2H9. 0168-874X/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 1 6 8 - 8 7 4 X ( 9 4 ) 0 0 0 7 0 - 0

F. Guibault et aL / Finite Elements in Analysis and Design 19 (1995) 309-324

310

be evaluated; and the most appropriate error estimator for a given set of equations can thus be found. Once the error estimate has been computed for a solution, the relationship of the local error to the grid size can be built, knowing the solution space and how the variables of the equations are interpolated. This provides the mechanism to adapt grid size in order to achieve error equidistribution through either global error reduction or by reaching a target relative error. Our goal, in this paper, is to show that the technology is mature enough that any elliptical equation can be described, solved and then analyzed for accuracy of the solution within a generic framework that is independent of the specific equation solved. For that purpose, we intend to describe the analysis technology for the specific case of the laminar Navier-Stokes equations. The same analysis tools will then be reapplied to laminar natural convection. The development of this generic framework has had a significant impact on the current work carried out within our research team. Most importantly, advancements in error estimation and grid adaptivity developed for one application have been carried over immediately to all the other fields of activity of the group; but this technology has also widened the range of applications our team could tackle efficiently. What we intend to demonstrate here, is how the physical model can grow in complexity, while the accuracy analysis tools remain the same. The actual error estimation process is based on a local projection scheme introduced by Zienkiewicz and Zhu I-1, 2-1, from which optimal size for mesh elements is computed and used in an adaptive mesh refinement scheme. This adaptive scheme tends to equidistribute a known relative error over the domain and thus provides a mechanism to control the accuracy of the solution.

2. Error estimation for elliptic equations In order to meet the accuracy requirements imposed on numerical solutions by practical applications, and this within reasonable computation costs, one must inevitably resort to some technique of local mesh concentration. Moreover, valid mesh concentration generally depends greatly not only on the geometry of the problem at hand, but also on the solution itself. This immediately points to the need of mesh regeneration or adaptation based on some analysis of the solution; the final goal being to equidistribute the overall error of the solution and thus achieve greatest accuracy at the lowest expense. The approach used in this study is twofold: first an error estimate is computed across the computational domain, and second, the mesh used for the computation is locally refined or coarsened according to the error estimate computed on each element. The error estimation process relies mainly on an appropriate norm of the solution based upon a priori knowledge of the Sobolev space H m in which the solution in sought. Indeed, any norm Ilullm of H m over a domain f2 becomes, when discretized with M dements: M

ilul12.,~ -- Y.IlutI,,~. 2

(1)

K

The norm of the error between the exact solution u and the numerical solution uh can thus be

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

311

Uhllm.

(2)

written as Ilell, = Ilu -

In practical applications, the exact solution is not known, and must be approximated using a projection operator ~ . The error can now be written as Ilellm = Ilu - uhllm I]~(Uh) --

uhl[m.

(3) (4)

On the other hand, the norm of the error on an element is linked to the size of the element by the relation Ilellm, x =

Ch~+1-'~,

(5)

where d is the del~:ree of the interpolation polynomials used in the finite element solution, C is a constant, h the size of the dement, and m the index of the Sobolev space. In this approach, the degree of the interpolation polynomials, the norm on the Sobolev space and the projector can all be changed, while the error estimation procedure remains the same. This generates two families of adaptive techniques, namely h and p. In h techniques, the grid size is modified, while the: degree of the polynomials is kept constant; while for p techniques, the degree of the polynomials is modified and the grid remains unchanged. Still, for any of these adaptive techniques, and for a given Sobolev space H ' , many projection operators can be considered l"31, that will produce so many ways of computing a more accurate solution from the numerical solution. The class of operators we will actually use in this work are based on the fact that a solution sought in H" will have its m - 1 first derivatives continuous, while its mth derivative will be discontinuous; one approach in building a "better" solution thus consists in projection this mth derivative at every node of the mesh to make it continuous. Of course, there is not only one way of performing this projection, and some projector might work better on some problems, while others are more appropriate in other circumstances. A first, obvious, way of projection a discontinuous function u~' at every node, is to compute that function at the node on every element, and assign to the node the average value found among all elements connected to that node. This average can either be direct or weighted by the area (2D) or the volume (3D) of each element. Another way, would consist in solving a global projection problem on the domain, which, in fact, would amount to solving a least square problem giving a single value of the function at every node. Of course, this is a rather expensive way of making the function continuous, and a local approximation of that procedure might prove far more cost effective. This is exactly the method proposed by Zienkiewicz and Zhu [41, that was further studied by Labb6 [5,1. The operator is defined by the construction of a polynomial of degree p built from the evaluation of the solution and its derivatives at certain points on all elements adjacent to a given node of the mesh. Typically the points of evaluation are the ones provided by the integration rule used in the finite dement computation. In this scheme, both the solution itself and its derivatives are approximated in the same function space, namely the one used by the FEM. Of course, at least

312

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

p + 1 points are needed for the construction, and when more are available, the polynomial is determined in the least-square sense. The fact that this new function space is exactly the one used for the finite element solution, ensures the following properties of the projected solution: • The projected and numerical solutions are identical, as well as the first m - 1 derivatives of the solution: ~(uh) - uh = 0 .... , ~(u~'-1) _ u~'-1 = 0 everywhere on the domain. • The projected mth derivatives of the solution are represented using a richer function space than the original numerical derivatives and should thus be more precise: i/gmUh~ ~m Uh ~{"Z"-gm/\OXF/--C~X------~#

O.

(6)

This implies that calculating the n o r m of the error in H " reduces to the evaluation of the semi-norm:

(~Uh) ~Uh2

~muh~ ~muh

(7)

Ilel12,,~ ~ II~(uh)- uhll2,~ + I1~ ~X~ --7X~ II°'° + "'" + II~k~xT' ] -- ~--~-~ 0.~.

:

I1~ ~mUh ~ll

(8)

oo

Once an error estimate becomes available on every element of the mesh, on optimal size can be determined for the elements. Indeed, for a given relative error level %pt, we have

Ileoptll,,,.o flop, = Ilullm.o

(9)

and from Eq. (1) ][eoptllm,K ~- ~]°Pt ][UhI]m'f~

M1/2

(10)

F r o m Eq. (5), we can link the optimal error on the mesh to an optimal size for each element

[leoptllm,r

=

Chdop+ 1 - m

(11)

and by assuming that the constant C does not vary from one mesh to the next, express the optimal size of the element as a function of its current size, current error and optimal error: hdp+l - m = h~c+x-,n IIeoptIIm,r

(12)

[[ellm.K Finally, these optimal sizes for every element of the domain can be fed to a mesh adaptation program that will compare sizes of elements with the provided targets and combine, split or otherwise regularize the mesh so that is best fits the solution requirements.

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

313

3. Computational scheme for elliptical equations The solver used for this study is based on a highly modular finite element environment coined as PIRATE (Programmable, Interactive, Resolution, Analysis and Treatment Environment) [6]. This environment enables the rapid prototyping of new projection or differential operators, new norms and even of completely new sets of equations as they are needed by the numerical application; it also promotes the seamless integration of the different components (or modules) needed by the solution cycle such as geometric modelers, grid generators, solvers and analysis packages. For the purpose of this study, a solution strategy has been developed that will be applied with little variation to all the problems considered. • Loop on increasing adimensional parameter values. - While the error is not equidistributed over the mesh 1. Solve the numerical problem by FEM. 2. Estimate the error over each element for selected unknowns 3. Take the minimum computed size of the element and modify the mesh accordingly. 4. Determine if the error is equidistributed. Examples of adimensional parameters are the Reynolds and Rayleigh numbers, usually chosen for isothermal or thermally driven flows. Elemental statistics are used to assess error equidistribution; for each selected unknown, minimal, maximal and average error are computed, as well as the variance. After a predetermined number of adaptive cycles, the statistics are used to determine whether the target error is satisfied, or further computations are required.

4. Laminar viscous flow - equations This section describes the initial set of equations and the discretization techniques on which the error control scheme is put to work. These equations form a basis from which other physical models will be built through the addition of supplementary equations. So, let's introduce the Navier-Stokes equations for incompressible laminar flow: p V . V V = 17,a,

(13)

I7. V = 0

(14)

with a=

- Pg + #z,

z=

VV + vTv,

which are discreti:zed using a Galerkin weighted residual method. Throughout these computations, the Crouzeiz-Raviart element was chosen. For this set of equations, the solution for the velocity components is sought in H a, and the norm on the error then becomes, for one of these components (for instance u): Ile, ll~,~ = II:'(uh) - uhllo2,a+ II 17(~'(uh)) - Vuhllg,o.

(15)

314

F. Guibault et al. / Finite Elements in Analysis and Design 19 (1995) 309-324

This enables us to compute, for each velocity component, an optimal element size; and in turn, the minimum of the two computed sizes, is used as the optimal size for the element.

5. L a m i n a r viscous flow - Results

The solution procedure and error control methodology has first been applied to laminar viscous flow, and for this set of equations, the only dimensional parameter to consider was the Reynolds number. The unknowns of the problem are the velocity components, which are computed in a fully coupled manner, and the pressure, which is computed as the Lagrange multiplier for the incompressibility constraint, and is imposed through an augmented Lagrangian formulation using Uzawa's algorithm. This section presents some laminar test cases and comparisons with available experimental data, as well as a parametric study of the flow inside a 2D linear cascade for various angles of attack. These three results show how grid adaptivity can be put to work to validate an error estimator and to compute accurately aerodynamic forces.

5.1. Backward facing step The first test case consisted to solve the laminar flow inside an expanding channel. Experimental data from Armaly [7] is available for that test case for 2D laminar flow up to a Reynolds number of 389. At that Reynolds number, the reattachment point is situated downstream of the expansion at X/S ~- 7.65 where S is the height of the expansion step. The Reynolds number is referenced to the average inlet velocity and the height of the step (S), which represents approximately 48% of the channel height. The experimental data comprises measurements of axial velocity across the channel height at four stations downstream of the step. The positions of the measurement stations are referenced in terms of the distance from the step, normalized by S. We will, however, show results only for stations X/S = 2.55 and X/S = 7.76, which are a representative sample of our results. Figs. 1 and 2 present the horizontal component of velocity along a vertical cut at two positions downstream of the expansion. Good agreement with the experimental results can be observed even in the region of the reattachment point for a relatively modest number of elements (2400). Fig. 3 shows the adapted grid resulting from the above procedure satisfying a target error of 10% on both velocity components; the pressure is not included in the error analysis.

5.2. Obstructed channel The next test case was an obstructed channel for which experimental data was also available. Simulations were conducted at Re = 82.5 in order to compare with the experimental data of Carvalho [8]. At that Reynolds number, the flow is bidimensional, laminar, stationary and it presents an important jet downstream of the obstacle, with two large recirculation zones above and under it. The Reynolds number is referenced to the average inlet velocity and the height of the obstruction (S), which represents 75% of the channel height. The experimental data comprises measurements of axial velocity across the channel height at six stations downstream of the obstacle.

D~ Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324 i

i

i

i

i

i

i

315

i

I

Exp. Local Prj. - -

0.13

O . dS-"

0.4

0. :2

0.001

0.002

0.003

0.004

0.005 Y (M)

0.006

0.007

0.008

0.009

.01

Fig. 1. Backward facing step: X component of velocity for X/S = 2.55.

1. :3

i

i

Exp. O Local Prj. - -

1.6

1.,4

I. 12

l

0.8

0.6

0.4

0.2

0 0

I

I

I

I

i

I

I

I

i

0.001

0.002

0.003

0.004

0.005 Y (M)

0.006

0.007

0.008

0.009

Fig. 2. Backward facing step: X component of velocity for X/S = 7.76.

0.01

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

316

Fig. 3. Backward facing step: adapted mesh.

i

i

i

Exp. o maillage f i n - -

0.8

0.6

<)

0.4

0.2

<)

<)

0

-0.2

0

I

I

I

I

I

I

I

I

I

0.001

0.002

0.003

0.004

0.005 Y (M)

0.006

0.007

0.008

0.009

0.01

Fig. 4. Obstructed channel: X component of velocity for X/S = 5.0, uniform grid.

The positions of the measured stations are referenced in terms of the distance from the obstacle, normalized by S. Even though this might not seem so at first glance, this second problem is far more chalenging than the first. Indeed, for the backward facing step problem, an adequate mesh can be build without any recourse to adaptive, grid technology. On the other hand, for the obstructed channel, the brute force approach, which consists of using a uniform but very fine grid throughout the computational domain, will general fail, as illustrated by Fig. 4. This solution was obtained on a uniform grid of 35 985 nodes, and as might be observed, agreement between experimental and numerical results is quite poor. Fig. 5 shows a remarkable improvement on the previous solution, using an adapted grid of only 4051 nodes which satisfies a target error of 10% on the velocity components. Still, numerical results remain quite distant from measurements at the reattachment point. Indeed, at that station, velocities are very low and the flow is mostly governed by pressure effects. This constitutes

317

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324 I

I

m~-p. Local Prj. - -

0.8

0.6

0.4

0.2

0

-0.2 I

I

I

I

I

I

I

I

I

0.001

0.002

0.003

0.004

0.005 Y (M)

0.006

0.007

0.008

0.009

0.01

Fig. 5. Obstructed channel: X component of velocity for X / S = 5.0, adapted grid.

a warning about the sensibility of the error estimator to detect error on pressure driven flow. In fact, this hypothesis was studied by H6tu [3], whose estimator includes implicitly the pressure. However, his best result only shows marginal improvement over our local projection. Outside this critical region, excellent agreement is observed with experimental data, as shown in Figs. 6 and 7.

5.3. Laminar 2D cascade flow Our goal is to compute the forces acting on the cascade blades as a function of the inlet flow angle for values of 0 °, 15 °, 30 ° and 45 °. The Reynolds number was set to 103 based upon the blade to blade distance and the inlet velocity magnitude. The simulations were realized for a target error of 10% and 5% to estimate the precision of these forces. As suggested by Szabo [9], forces on the blades were computed by summing up the reactions to the no-slip boundary condition to guarantee maximum accuracy. These values were compiled as a function of inlet flow angle and are presented in Table 1. The best results are the ones at 5% target error, as illu,~trated by the value of Fy on the blade for an inlet flow angle of 0 °. By comparing these two sets of walues, it can be estimated that the results for 5% target error are precise up to 2 significant digits.

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

318

i

Local

i

Exp. Prj.

o --

0.8

0.6

0.4

0.2

-0.2 I

I

I

I

I

I

I

I

I

0.001

0.002

0.003

0.004

0.005 Y (M)

0.006

0.007

0.008

0.009

0.01

Fig. 6. Obstructed channel: X component of velocity for X/S = 0.0, adapted grid. Obstacle !

: X/S !

= 4.0 !

i

i

Prj.

i

Local

0.8

0.6

O

0.4

0.2

-0.2 --Y-

0.001

I

0.002

,

I

I

I

I

I

I

I

0.003

0.004

0.005

0.006

0.007

0.008

0.009

Y

(M)

Fig. 7. Obstructed channel: X component of velocity for X/S = 4.0, adapted grid.

0.01

F. Guibault et al./Finite Elements in Analysis and Design 19 H995) 309-324

319

Table 1 Forces o n the blades at target e r r o r of 10% a n d 5 % Angles

. . . . . .

Fx

Fy

Target error of10% 45 ° 0.056852 30 ° 0.053559 15 ° 0.053796 0° 0.074369

0.385199 0.305316 0.199691 -0.000226

T a r g e t e r r o r of 5 % 45 ° 0.056895 30 ° 0.053664 15 ° 0.053754 0° 0.074326

0.386358 0.304962 0.198919 0.000002

\ ~

X~-TN......7~

Fig. 8. L a m i n a r cascade flow: grids for 0 °, 15 °, 300,45 °.

Fig. 8 shows the grid evolution as a function of the inlet flow angle for an error of 5%. Periodicity, the boundary layer on the pressure side, the leading edge separation and the wake are

320

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

f

13

r

?

f

----

Fig. 9. Laminar cascade flow: IIV II isovalues for 0°, 15°, 30°,45°. well captured. When the inlet flow angle is increased, the flow separates in two regions, one of low velocity, and one of high velocity, as illustrated in Fig. 9. The viscous stress at the interface of these two regions is incompatible with the outlet boundary condition, which explains the unusual concentration of grid points at the outlet for high angles of attack.

6. Natural convection - Equations

The second set of equations we have analyzed with our error control tool is obtained from the incompressible Navier-Stokes equations of the previous section by adding a new transport equation for temperature and a source term in the m o m e n t u m equations. Indeed, under the Boussinesq approximation, natural convection phenomenon can be described by adding a transport equation for temperature, which can be written as: p C p V . V T + V. (k V T ) = 0

(16)

and the source term in the m o m e n t u m equation, that is fl(T-

To)g

(17)

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

321

with g the gravitational acceleration vector. Again, this new equation and supplementary term are discretized with a Galerkin weighted residuals method, using for temperature the same quadratic interpolant that was used for velocity components. In this application, the physics has increased in complexity, hence, to accurately capture the thermal boundary layer as well as the m o m e n t u m one, error analysis on the temperature field must be performed. Since the physics has changed, the projector and the norms must be validated once again before being applied to a new type of application. This section presents a laminar test case to validate the error analysis tools in the context of natural convection, as well as parametric study of the thermally driven flow between two concentric cylinders. Here, the analysis issue is to accurately compute the average Nusselt number on the inner cylinder as a function of the Rayleigh number.

7. Natural convection- Results For this second set of equations, all the basic solution and error control algorithms remain exactly the same. The solution vector space for T is again H 1, and all the norms in the error analysis can thus be used directly. When computing the optimal size for each element, we then only need to consider the optimal size estimated by the norm on temperature together with the sizes calculated for each velocity component, and use the minimum of the three. 7.1. Natural convection in horizontal annuli

This test case consisted to solve the thermally driven flow between two concentric cylinders. The inner and outer cylinders are kept at constant temperatures; experimental data from Kuehn [10] is available for that test case at a Rayleigh number of 50 000 and Prandtl number of 0.7. Comparisons are based on the equivalent thermal conductivity (Keq), the ratio of the local Nusselt number at Rayleigh 50 000 to the local Nusselt number at Rayleigh 1, on the inner and outer cylinders. Results are presented in Fig. 10, where the angle 0 is a measure of the position on the cylinders, 0 = 0 ° being the topmost position, and 0 = 180 °, the bottommost. Simulations were carried out for target error of 2.5 and 15%. The numerical results are in very good agreement with the experiments for both error levels, and show that a target error of 15% is adequate and cost effective in terms of computing resources. 7.2. Natural convection in internally finned horizontal annuli

This parametric study aims at determining the influence of fins on the performance of a simple heat exchanger. For that purpose, simulations were carried out to compute the average Nusselt number on both the inner and outer cylinders of a finned annuli at various Rayleigh numbers. These computations span a range of Rayleigh numbers from 104 up to 10 7, for target errors of 15% and 10%. F r o m these computations, it is expected that optimal configuration of fins can be determined with respect to both their number and length. Table 2 presents the average Nusselt number as a function of Rayleigh number, computed by summing up the reactions on the inner and outer cylinders from the imposition of the temperature

322

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

12

i

Rayleigh 50000 i i

i

i

i

J

~x,

2.5-15 . . . . Into

8 ~.,..

ext +

10

E

r

6

\

',,+

4

"". .

.

.

.

.

.

.

.

.

.

.

~_

_o

o

140

160

2

o

20

40

60

80

0

100

120

180

Fig. 10. Natural convection in horizontal annuli: equivalent conductivity.

Table 2 Average Nusselt numbers at target error of 15% and 10% Rayleigh numbers Target error of 15% 104 105 106 107 Target error of 10% 104 105 106 107

NUinne r

NUoute r

2.0983 3.7106 6.7058 12.1275

1.1791 2.0896 3.8039 6.8225

2.0970 3.7077 6.7204 12.1618

1.1824 2.1006 3.8154 6.8954

as an essential b o u n d a r y condition. By comparing these two sets of values, it can be estimated that the results for 10% target error are precise up to 2 significant digits. Fig. 11 shows the adapted grid near the fins of a typical heat exchanger. Observe the clustering of nodes near fin tips, which can be attributed to the discontinuity in the gradient of the temperature field at these points. This is a side effect of the underlying hypothesis incorporated in the construction of the error estimator, which assumes that discontinuities in gradients will vanish as the d e m e n t size tends to zero. Of course, this assumption will always prove wrong in the neighborhood of a physical discontinuity. In that case, the error analysis will not converge to

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

323

Fig. 11. Natural convection in finned annuli: adapted grid.

a stable grid since the projector ~ constructs a continuous approximation of the gradient of the solution (Eq. (4)).

8. Conclusion Grid adaptation and equirepartition of the error was introduced a few years ago by many authors (I-1, 3] etc.). But, as we have shown, new applications of this methodology constantly arise. Hence, a framework for the coherent introduction and testing of these analysis tools was strongly advocated in our research group; in particular since we must fulfill a double mandate of both academic research and technological transfer to industrial partners. The test cases presented in this paper illustrate some possible applications in the field of engineering and the underlying difficulties which may come from the use of this methodology. Even though great benefit can be gained from its use in design and analysis, research must be pursued to find the error estinaator that best captures relevant physical phenomenons of the problem being investigated. In this identification process, complementary laboratory experiments will always play a crucial role.

324

F. Guibault et al./ Finite Elements in Analysis and Design 19 (1995) 309-324

The resolution framework mentioned here plays an important role in the development of new estimators by reducing the time needed for testing and for integration into real life engineering applications. One major asset of the resolution environment, lies in its configurable data extraction capabilities - as presented in the companion paper - which enable the high level description of the estimators to be tested. Indeed without resorting to any traditional code development (e.g. C or C + +) performance of new estimators can be graphically visualized and their efficiency quantitatively rated. Work is currently under way to apply this same methodology to investigate two equation turbulent flow models, such as k-e and k - o , and to identify the best estimator for each model.

References [1] O.C. Zienkiewicz and J.Z. Zhu, "The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique", Int. J. Numer. Methods Eng. 33(7), pp. 1331-1364, 1992. [2-10.C. Zienkiewicz and J.Z. Zhu, "The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity", Int. J. Numer. Methods Eng. 33(7), pp. 1365-1382, 1992. 1-3-1 J.F. Hrtu, Mrthodes d'61ements finis adaptatives pour les 6coulements visqueux incompressbiles, Ph.D. Thesis, Ecole Polytechnique de Montrral, 1991. [4] O.C. Zienkiewicz, J.Z. Zhu, and N.G. Gong, "Effective h-p-version adaptive analysis procedures for the finite element method", Int. J. Numer. Methods Eng. 28, pp. 879-891, 1989. 1-5] P. Labbr, Analyse per 616ments finis de la rrsistance 61ectrique dans les joints brasrs, Mrmoire de maitrise, Ecole Polytechnique de Montrral, 1993. [6-1 F. Guibault, R. Magnan, A. Garon and R. Camarero, "PIRATE: A CFD research environment", in: Proc. lOth Int. Conf. on Computing Methods in Applied Sciences and Engineering, Paris, France, Nova Science, 1992. [7] B.F. Armaly, F. Durst, J.C.F. Pereira and B. SchSnung, "Experimental and theoretical investigation of backward facing step", J. Fluid Mech. 127, pp. 473-496, 1983. [8] M.G. Carvalho, F. Durst and J.C.F. Pereira, "Predictions and measurements of laminar flow over two-dimensional obstacles", Appl. Math. Modelling 11, pp. 23-34, 1987. [9] B.A. Szabo and I. Babuska, Finite Element Analysis, Wiley, New York, 1991. [10-1 T.H. Kuehn and R.J. Goldstein, "Natural convection between horizontal concentric cylinders, Part 4", J. Fluid Mech. 74, pp. 675-719, 1976.