Journal of the Brazilian Society of Mechanical Sciences and Engineering
Print version ISSN 16785878
J. Braz. Soc. Mech. Sci. & Eng. vol.27 no.4 Rio de Janeiro Oct./Dec. 2005
http://dx.doi.org/10.1590/S167858782005000400005
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}
^{I}Universidade Federal de Pernambuco; Rua Acadêmico Hélio Ramos S/Nº; 50740530 Recife, PE. Brazil; gkc@demec.ufpe.br; prmlyra@demec.ufpe.br
^{II}Universidade Federal de Pernambuco; Av. Prof Luiz Freire nº 1000; 50740540 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 NavierStokes 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 spacetime 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: NavierStokes 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 PetrovGalerkin) formulation for the simple advectiondiffusion 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 SUPGlike 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 predictormulticorrector algorithm, where an augmented SUPG matrix is, then, proposed. Finally, numerical results are given and compared with results published in the literature.
Nomenclature
The NavierStokes Equations
Fluid flows can be mathematically modeled by the NavierStokes 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 quasilinear form:
where:
and K_{kh} is such that:
Expressions for the matrices A_{0}, A_{k} and K_{kh} can be found in the appendix.
The twodimensional NavierStokes equations can not be solved without a further state equation since there are four equations and five variables (r,p,q, u_{1} and u_{2}). 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 NavierStokes 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_{0} in the case of liquids, when the necessary derivatives are not easily obtained from Eq.(20).
SpaceTime Formulation
A solution Y = Y(x,t) for the Eq.(1) is sought in the spacetime 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 spacetime slabs (Shakib, 1991). Figure 1 illustrates the partitioned domain. The n^{th} spacetime 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 spacetime "border" by , the superscript e , when used, indicating that those entities are taken elementwise. 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 ConstantinTime 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 spacetime 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 PetrovGalerkin 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 predictormulticorrector 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 predictormulticorrector 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 subindexes 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 (1a) 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^{2} (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 predictormulticorrector algorithm (see Shakib, 1988), so, the prefix "multi" may be dropped and we have, in fact, a predictorcorrector 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 CourantFriedrichsLewy 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 nonslip 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:_{u}_{k} = 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^{3}, C_{p} = 1005.0J/(kgK), c_{v} = 718J/(kgK), b = 1.0x10^{6}Pa^{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^{2} and k = 2700.0W/mK whereas for Re=400, µ = 0.005Ns/m^{2} 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.
Figure 2 shows the unstructured mesh that was used. The mesh has 3688 triangular elements and 1936 nodes and is somewhat refined near the lid to smoothen oscillations in the pressure field (Hauke, 1995). A similar sized mesh was used by Hauke (1995) though, in his case, the mesh was structured with 40 x 40 square elements, slightly refined near the two upper corners. Ghia et al. (1982) used a much more refined structured mesh (121 x 121) and a multigrid technique to solve the same problem.
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_{1} and u_{2} 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} = 1.0, u_{2} = 0.0 q = 2.769 x 10^{4}K and r = 1.0kg/m^{3}

on C (0.2 < x < 1.2 and y = 0.8): u_{1} = 1.0, u_{2} = 0.0, q = 2.769x10^{4}K and r = 1.0kg/m^{3}

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_{1} = 0.0, u_{2} = 0.0, q = 7.753x10^{4}K
The initial conditions have been chosen to match the values on boundary B.
Results are plotted after 3500 interactions when steadystate 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.
Appendix: Matrices Used in This Paper
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:
References
Aliabadi, S. K. & Tezduyar, T. E., 1995, "Parallel fluid dynamics computations in aerospace applications", International Journal for Numerical Methods in Fluids, Vol 21, pp. 783805. [ Links ]
Almeida, R. C. C., 1993, "Uma formulação de PetrovGalerkin para a resolução das equações de Euler e NavierStokes compressível usando técnicas adaptativas", Ph.D Thesis (in portuguese), COPPE/UFRJ, Rio de Janeiro, Brazil, 105p. [ Links ]
Anderson, J. D., 1995, "Computational fluid dynamics: the basics with applications", McGrawHill, EUA, 550p. [ Links ]
Bohnaus, D. L., 1998, "A higher order accurate finite element method for viscous compressible flows", Ph.D Thesis, Virginia Polythecnic Institute and State University, Blacksburg, USA, 78p. [ Links ]
Brooks, N. & Hughes, T. J. R., 1982, "Streamline Upwind / PetrovGalerkin formulations for advection dominated flows with particular emphasis on the inccompressible NavierStokes equations" Computer Methods in Applied Mechanics and Engineering, Vol 32, pp.199259. [ Links ]
Codina, R., 1993, "A finite element formulation for the numerical solution of the advectiondiffusion equation", Monography, International Center for Numerical Methods in Engineering (CIMNE), Barcelona, Spain, 120p. [ Links ]
De Sampaio, P. A. B, Lyra, P. R. M, Morgan, K. & Weatherhill, N. P., 1993, "PetrovGalerkin solutions of incompressible NavierStokes equations in primitive variables with adaptive remeshing" Computer Methods in Applied Mechanics and Engineering, Vol 106, pp 143178. [ Links ]
Ghia, U., Ghia, K. N. & Shin, C., T., 1982, "Highre solutions for incompressible flow using NavierStokes equations and a multigrid method", Journal of Computational Physics, Vol 48, pp. 387411. [ Links ]
Gustafsson, B. & Stoor, H., 1991, "NavierStokes equations for almost incompressible flow" SIAM Journal on Numerical Analysis, Vol 28  6, pp 15231547. [ Links ]
Hauke, G. & Hughes, T. J. R., 1994, "A unified approach to compressible and incompressible flows" Computer Methods in Applied Mechanics and Engineering, Vol 113, pp 389395. [ Links ]
Hauke, G., 1995, "A unified approach to compressible and incompressible flows and a new entropyconsistent formulation of the Kepsilon model", Ph.D Thesis, Stanford University, Stanford, USA, 185p. [ Links ]
Hughes, T. J. R. & Mallet, M., 1986, "A new finite element formulation for computational fluid dynamics:III. The generalized streamline operator for multidimensional advectivediffusive systems" Computer Methods in Applied Mechanics and Engineering, Vol 58, pp 305328. [ Links ]
Lyra, P. R. M., 1994, "Unstructured grid adaptive algorithms for fluid dynamics and heat conduction", Ph.D. Thesis, University of Wales, Swansea, UK, 333p. [ Links ]
Mittal, S. & Tezduyar, T., 1998, "A unified finite element formulation for compressible and incompressible flows using augmented conservation variables" Computer Methods in Applied Mechanics and Engineering, Vol 161, pp 229243. [ Links ]
Moussaoui, F., 2003, "A unified approach for inviscid compressible and nearly incompressible flow by least squares finite element method" Applied Numerical Mathematics, Vol 44, pp.183199. [ Links ]
Shakib, F., 1988, "Finite Element analisys of the compressible Euler and NavierStokes equations", Ph.D. Thesis, Stanford University, Stanford, USA, 204p. [ Links ]
Shakib, F., Hughes, T. J. R. & Johan, Z., 1991, "A new finite element formulation for computational fluid dynamics: X. The compressible Euler and NavierStokes equations" Computer Methods in Applied Mechanics and Engineering, Vol 89, pp 141219. [ Links ]
Schlichting, H., 1979, "Boundary Layer Theory", McGrawHill, Antigua, Guatemala, 817p. [ Links ]
Wark, J. S., 1995, "Advanced thermodynamics for engineers", McGrawHill, EUA, 564p. [ Links ]
Yoon, K. T., Moon, S. Y, Garcia, S.A., Heard, G. W. & Chung, T. J., 1998, "Flowfield dependent mixed explicitimplicit (FDMEI) methods for high and low speed and compressible and incompressible flows" Computer Methods in Applied Mechanics and Engineering, Vol 151, pp 75104. [ Links ]
Zienkiewicz, O.C., Szmelter, J. & Peraire, J., 1990, "Compressible and incompressible flow; an algorithm for all seasons" Computer Methods in Applied Mechanics and Engineering, Vol 78, pp.105121. [ Links ]
Paper accepted July, 2005.
Technical Editor: Aristeu da Silveira Neto.