Abstract

TECHNICAL PAPERS

Numerical simulation of two dimensional compressible and incompressible flows

G. K. Costa^I; P. R. M. Lyra^I; C. A. B. de Oliveira Lira^II

^IUniversidade Federal de Pernambuco; Rua Acadêmico Hélio Ramos S/Nº; 50740-530 Recife, PE. Brazil; gkc@demec.ufpe.br; prmlyra@demec.ufpe.br

^IIUniversidade Federal de Pernambuco; Av. Prof Luiz Freire nº 1000; 50740-540 Recife, PE. Brazil; cabol@ufpe.br

ABSTRACT

In this article, we make use of a stabilized Finite Element method to solve the complete set of Navier-Stokes equations. The methodology adopted is such that it allows for the use of different sets of variables, particularly the so called conservative and pressure variables. A space-time formulation using a simple augmented SUPG stabilizing term is proposed for the particular case of pressure variables. Comparison with data published in the available literature is done and a reasonably good agreement is obtained.

Keywords: Navier-Stokes equations, compressible flows and incompressible flows, SUPG

Introduction

Computation of fluid flows through the Finite Element Method has been the subject of much research over the last decades. The development of the methodology to be used throughout this paper can be traced back to the SUPG (Streamline Upwind Petrov-Galerkin) formulation for the simple advection-diffusion equation (a milestone reference on the theme is the work of Brooks & Hughes, 1982). Several attempts to extend the method to systems of equations can be reported since then, giving rise to different versions of SUPG-like formulations and discontinuity capturing operators, which are added to account for shocks in the case of compressible flows (see, for instance, Hughes & Mallet, 1986 and Shakib et al., 1991). In addition to that, is the fact that basically two different approaches are at sake: one for compressible and another for incompressible flows, this last one usually causing the decoupling of the energy balance equation and reducing the number of variables involved. Attempts have been made, however, towards a unified approach for both flow regimes (see, for instance: Moussaoui, 2003, Mittal & Tezduyar, 1998 and Hauke & Hughes, 1994. For references on works using techniques other than Finite Elements only, the reader can consult, for example: Gustafsson & Stoor, 1991 and Yoon et al., 1998). In this work, we follow the methodology introduced by Hauke & Hughes (1994), introducing a SUPG stabilizing matrix made out of two simple terms, one of them suitable for compressible flows and the other, for incompressible flows, in an attempt to deal with the two regimes with a single formulation. The resulting approach is such that all the element matrices and vectors can be analytically evaluated, whereas in typical formulations, the stabilizing term generally needs to be evaluated numerically (Hauke, 1995). Results show that this approach has given reasonably good results, when comparing to others found in the literature, though, due to the very simple stabilizing matrix, our formulation has revealed to be over diffusive in both regimes.

This paper is divided as follows: first a brief review of the governing equations is given, followed by the finite element formulation presentation. The drawbacks of the compressible scheme when used for incompressible flows, are discussed within the context of a formulation which makes use of a predictor-multicorrector algorithm, where an augmented SUPG matrix is, then, proposed. Finally, numerical results are given and compared with results published in the literature.

Nomenclature

The Navier-Stokes Equations

Fluid flows can be mathematically modeled by the Navier-Stokes equations which, in the absence of body forces and internal heat generation, can be written in vector conservation form as (Anderson, 1995):

where nd is the number of spatial dimensions, U is the vector of conservation variables, F_k is the advective flux vector and D_k is the diffusive flux vector.

Each of these vectors: U, F_k and D_k, may be written as a function of a general set of variables Y, that is: U = U(Y), F_k = F_k(Y) and D_k = D_k(Y). In this paper,Y will always refer to pressure primitive variables, so that:

where:

r is the density

u_k is the kth component of the velocity vector u

q is the absolute temperature

p is the absolute pressure

E is the total energy per unit mass:

In Eq. , e is the internal energy per unit mass. It is related with the temperature q by:

where c_v is the specific heat at constant volume.

Also,k_c is the kinetic energy per unit mass:

The advective and diffusive flux vectors are given by:

where are the components of the viscous tensor t^v, q_k is the kth component of the heat flux vector q and d_kh is the Kronecker delta. For a Newtonian fluid, we can write:

where µ is the dynamic viscosity and l is the molecular viscosity.

In Eq. , we assumed that the Stoke's hypothesis was valid, that is:

For a detailed description of equations to , the reader can consult Schlichting (1979).

It is generally admitted that the viscosity is constant in the case of incompressible flows. For compressible flows, specifically speaking of perfect gases in the common regimes of interest, it is possible to relate the viscosity with the absolute temperature through the empirical Sutherland's law (see, for instance, Anderson, 1995).

The Fourier law relates the heat flux with the absolute temperature:

where k is the thermal conductivity tensor. For an isotropic fluid, the tensor k can be substituted by a scalar value, k. This assumption will be used hereafter.

It is important to note that Eq. can also be written in its quasi-linear form:

where:

and K_kh is such that:

Expressions for the matrices A₀, A_k and K_kh can be found in the ^appendixappendix.

The two-dimensional Navier-Stokes equations can not be solved without a further state equation since there are four equations and five variables (r,p,q, u₁ and u₂). In this paper, we assume that the fluid can be either a liquid or a gas, and make use of the thermodynamic Maxwell relations (Wark, 1995) to obtain:

where C_p is the specific heat at constant pressure, a_p and b are the isochoric expansivity and the isothermal compressibility, respectively, and are given by:

where: v is the specific volume.

Equations to are enough to complement the system of Navier-Stokes equations. It can be shown that Eq. reduces to the perfect gas law: p = rRq, once the hypothesis of a perfect gas is made. Equations to are particularly useful for the determination of the Jacobian matrix A₀ in the case of liquids, when the necessary derivatives are not easily obtained from Eq.(20).

Space-Time Formulation

A solution Y = Y(x,t) for the Eq.(1) is sought in the space-time domain Q = W x [0,T]. From this moment on, we define that all the approximations to geometric entities, variables and functions will be marked with an over "hat". Having said this, let be an approximate partition of the domain Q such that is composed by space-time slabs (Shakib, 1991). Figure 1 illustrates the partitioned domain. The n^th space-time slab partition can be represented as:

where is the spatial approximate partition at time level n and l_n is an open time interval.

The time discontinuous Galerkin approach, which will be used hereafter, allows the approximate solution to be discontinuous at time level t_n. We then define l_n in Eq. as:

The whole approximate domain is written as:

Figure 1 also highlights a prismatic element at time level n. The elementary domain being:

Let us also denote the border of the spatial domain by and the space-time "border" by , the superscript e , when used, indicating that those entities are taken element-wise. We finally observe that nodes need not have the same spatial location at each time level. This particularity makes the time discontinuous Galerkin approach specially suitable for moving boundary problems (Aliabadi and Tezduyar, 1995).

A Constant-in-Time Approach

In order to write the finite element variational form of Eq.(1) , let S^h_n and V^h_n be, respectively, the set of trial functions and the space of weighting functions, taken at time level n and given by:

where represents the portion of the space-time border, at time level n, where Dirichlet boundary conditions are prescribed.

The Time Discontinuous Galerkin method can be stated as: given W Î and g(t) Î Â^m, we search for Î such that (Shakib, 1991):

where, and is the kth component of the vector N_k , normal to the spatial domain.

We approximate the solution vector in the usual manner, that is:

where N_n is the number of nodes in the spatial domain, j_j is the shape function corresponding to node j and is the nodal value of the solution vector. Let us also define the nodal solution vector as being m ´ N_n the hypervector containing the solution at each node.

As for the weighting function W, a Petrov-Galerkin definition is due, to avoid the well known oscillations and to stabilize the method:

where I is the identity matrix and P_i is the SUPG perturbation term, given by:

where t is the stabilizing SUPG matrix.

The shape function j_j may be a polynomial function of x and t, though it is not strictly necessary to be so (see, for instance, Bonhaus, 1998). In our work we choose it to be linear in x and constant in t, that is:

Such assumption, considerably simplifies the formulation given by Eq.(29) since it implies that:

Substitution of equations (30) through (34)in Eq.(29) finally gives:

where and d_c is the shock capturing operator to be defined later on. The sum in Eq. is necessary because the stabilizing terms are not defined in the element boundaries.

Equation (35) can be solved with the use of a predictor-multicorrector algorithm (Shakib, 1988) (the reader can consult that reference for details). By applying a change of variables U ® Y, and considering the counterpart of the nodal solution hipervector , the following linear system can ultimately be obtained:

where, for the case of a single correction step in the predictor-multicorrector algorith, one can write:

In Eq.(36) the residue R_i and the matrix M_i can be assembled from their corresponding elementary counterparts M_e and R_e which hold the following pattern:

The elements of the matrix M_e and of the vector V_e can be easily determined for the case of pressure variables. Although their expressions will not be given here, it is worth noticing that neither the SUPG stabilizing matrix t, nor the discontinuity capturing operator d_c are well defined for purely incompressible flows, when both a_p and b, given by equations and become zero (the reader may consult Hauke, 1995, for further details). The key for developing a methodology that embraces both compressible and incompressible flows, therefore, reside in handling these two terms of stabilization in such a way that both limits are well attained. In this work, we propose two things:

a) To make use of an augmented SUPG stabilizing matrix t, given by:

where the sub-indexes in and cp stand for incompressible and compressible and a is a function of the Mach number M that must be such that:

b) To annul the effect of the discontinuity capturing operator whenever the flow becomes incompressible. This can be done if we replace the operator d_c in Eq. by the augmented operator (1 -a)d_c.

In Eq.(40), t_in is suited for the incompressible limit whereas t_cp is valid for compressible flows, so that the both extremes of the spectrum are covered.

Hauke (1995) used a full stabilizing matrix for both t_in and t_cp. We propose to use a very simple design for t_cp and a diagonal pressure based d_in. As will be seen in the numerical examples, such an approach gives relatively good results for the usual testing problems.

Stabilization Matrices and DC Operator

Before proceeding, some important definitions are convenient. Let X be a (m.nd)x1 vector and M be a mxm matrix. We define the following norms:

Let us also define the gradient of the vector X with respect to the Cartesian coordinates:

Under the assumption that there is a mapping, between the Cartesian coordinates x_k® e_k and the natural coordinates x_k, we proceed to define the gradient of the vector e_k with respect to the natural coordinates as:

where J^-1 is given by:

Making use of these definitions, let us first write an expression for the discontinuity operator. We have chosen to follow Aliabadi & Tezduyar (1986) and use the following expression for d_c , slightly modified in order to accommodate only compressible flows:

where:

V being the set of entropy variables (see, for example, Shakib, 1988).

The factor (1-a) in Eq.(46) is very important, since becomes singular in the incompressible limit.

The SUPG stabilizing matrices are given by (see HAUKE, 1995 and Aliabadi & Tezduyar, 1995):

where:

In equations (51) to (54), above, h_e is a characteristic dimension of the element. It may be taken as:

where A_e corresponds to the element's area. In the previous equations, we also have that c is the speed of sound.

Numerical Examples

In what follows, a GMRES (General Minimum Residue) solver with reverse communication was used to obtain the solution of Eq. , in the predictor–multicorrector algorithm. The solver was obtained as a closed package from the public domain directory of CERFACS² (In French: "Centre Europeen de Recherche et de Formation Avancée en Calcul Scientifique"). The code used was copyrighted in 1998 and is a single precision package. Only a block diagonal preconditioner was used with the solver and the GMRES tolerance was set 0.01 in all the examples. Also, the dimension of the Krylov subspace was set to 50. Only one correction pass was used for the predictor-multicorrector algorithm (see Shakib, 1988), so, the prefix "multi" may be dropped and we have, in fact, a predictor-corrector algorithm.

Pressure variables were used in all the examples and, concerning the time marching algorithm, we observed the following equation for the global time step (Shakib, 1991):

where C_r stands for the Courant-Friedrichs-Lewy number (CFL) and the coefficients f_c and f_d are given by:

Leaky lid driven cavity. This is an example of a totally incompressible flow, where the density r is assumed to be constant. A unit square recipient with infinite transversal length contains a hypothetic fluid which is driven to the left by the shearing tensions caused by the upper lid uniform rightwards motion, as shown in Fig. 2. The moving lid will produce a whirl inside the recipient (cavity). Hauke (1995) solved a similar problem for a Reynolds number of 1 and of 400 and we intend to compare our numerical results with his and with those published by Ghia et al. (1982). The boundary conditions applied are as follows: at the upper lid (C), a unit velocity to the right plus a reference temperature were prescribed (u1 = 1.0m/s and q = 0K), at the side and bottom walls (B, D and A) a non-slip condition was imposed and we also prescribed the pressure to be equal to zero in the mid point of the bottom wall. The problem marches to a steady state from a set of initial conditions which we have chosen to be:_uk = 0.0m/s, q = 0k and p = 0.0Pa. The fluid properties were all assumed to be equal to the air's at 293K and 101325Pa (r = 1kg/m³, C_p = 1005.0J/(kgK), c_v = 718J/(kgK), b = 1.0x10^-6Pa^-1 and a =0.0036K^-1). The viscosity, however, was allowed to vary, to adjust to the Reynolds number in each example and the thermal conductivity was also modified under the assumption of a constant Prandtl number of 0.75. Then, for Re=1.0 we have that: µ = 2.0Ns/m² and k = 2700.0W/mK whereas for Re=400, µ = 0.005Ns/m² and k = 6.7W/mK.

At the two upper corner nodes, the velocity field is discontinuous, and that adds an extra difficulty to the problem. We have chosen to set a linear gradual variation of the speed at the two corners, from 0,0m/s to 1,0m/s (Lyra, 1994 and De Sampaio et al., 1993). For linear shape functions, it suffices to set the speed at the two upper corner nodes to 1,0m/s whereas at their adjacent nodes belonging to the side walls, the speed is set to zero.

Steady state was reached after 20000 interactions for the case or Re=1 and 30000 interactions for the case of Re=400. The average global time step could not be chosen arbitrarily. We used C_r = 1 (see Eq. (57)) for Re = 1 and C_r = 500for Re = 400. No convergence was achieved when we tried to use C_r = 10 for the case of Re=1.0, indicating that there might be a maximum CFL value, over which there is no convergence. Figure 3 shows the velocity vectors and the pressure isolines obtained. A relatively good agreement is found between our results and those published in the literature. The oscillations in the pressure isolines near the top lid, which become more evident in the case of Re=400 can also be found in other authors' results (see, for instance, Zienkiewicz et al., 1990). Hauke (1995) seems to have reduced these oscillations drastically by making a local mesh refinement near the upper left and right corner nodes.

The behaviour of the velocity components u₁ and u₂ can be seen in Fig. 4 where some local approximate results obtained from a visual observation of the charts published in the work of Hauke (1995) are plotted for comparison for Re=1. In the case of Re=400, our results are compared with the results found in Ghia et al. (1982), usually considered as a benchmark for this particular example. We observe that there is a somewhat over diffusive behaviour of our solution when a comparison is made. This is more clearly seem in the case of Re=400 where advective effects become more pronounced and the stabilizing terms have a greater impact on the solution. Such over diffusive behaviour can be attributed to the simpler SUPG matrix used. Nevertheless, our results show a better agreement when compared to those obtained by Zienkiewicz et al. (1990), who used a similar mesh.

Supersonic flow past a flat plate. Next, we give a classic example of a compressible flow known as Carter's problem (Shakib, 1988). A Mach 3.0 flow of air (considered here as a perfect gas) past a fixed flat plate will be studied as shown in Fig. 5. The unstructured mesh used consists of 3707 elements and 1940 nodes and it is dimensionally equivalent to the mesh found in Almeida (1993), who also used an unstructured mesh.

For the dynamic viscosity, we take the expression for the Sutherland's law given in Almeida (1993):

The Prandtl number was assumed to be constant and equal to 0.75.

The boundary conditions for this problem are, as follows:

• on A (-0.2 < x < 0.0 and y = 0.0): u.n = 0.0, = 0.0 and q.n = 0.0

• on B (x = -0.2 and 0.0 < y < 0.8): u₁ = 1.0, u₂ = 0.0 q = 2.769 x 10^-4K and r = 1.0kg/m³

• on C (-0.2 < x < 1.2 and y = 0.8): u₁ = 1.0, u₂ = 0.0, q = 2.769x10^-4K and r = 1.0kg/m³

• on D (x = 1.2 and 0.0 < y < 0.8): free flow

• on E (0.0 < x 1.2 and y = 0.0): u₁ = 0.0, u₂ = 0.0, q = 7.753x10^-4K

The initial conditions have been chosen to match the values on boundary B.

Results are plotted after 3500 interactions when steady-state was reached. Figure 6 shows the Mach isolines and figures 7 and 8 show the density and temperature variation along the vertical axis which passes through x = 1.0m. Again, we notice that our method is more diffusive. The discrete points were obtained from Almeida (1993) in which the CAU (Consistent Approximate Upwind) method was used along with SUPG.

Conclusions

The procedure outlined above has proved to work for solving both compressible and incompressible flows with reasonable accuracy. The results show that in both regimes, reasonably good agreement is found between our results and those published in the literature, even with the use of simpler matrices though our results have shown to be over diffusive in both cases. There is much place for further work, mainly in what concerns to the SUPG stabilizing matrix and discontinuity capturing operator. Extension to three dimensional problems, which is a natural path from here, is already at course.

Acknowledgements

The authors would like to thank the Brazilian Research Council (CNPq) for the financial support granted.

All the matrices below can be found in Hauke (1995). Some parameters which are used in the matrices are given below:

Proceeding with the matrices, we have:

Paper accepted July, 2005.

Technical Editor: Aristeu da Silveira Neto.

appendix

Publication Dates

History