Print version ISSN 1678-5878
J. Braz. Soc. Mech. Sci. & Eng. vol.27 no.4 Rio de Janeiro Oct./Dec. 2005
Numerical simulation of two dimensional compressible and incompressible flows
G. K. CostaI; P. R. M. LyraI; C. A. B. de Oliveira LiraII
IUniversidade Federal de Pernambuco; Rua Acadêmico Hélio Ramos S/Nº; 50740-530 Recife, PE. Brazil; email@example.com; firstname.lastname@example.org
IIUniversidade Federal de Pernambuco; Av. Prof Luiz Freire nº 1000; 50740-540 Recife, PE. Brazil; email@example.com
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
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.
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, Fk is the advective flux vector and Dk is the diffusive flux vector.
Each of these vectors: U, Fk and Dk, may be written as a function of a general set of variables Y, that is: U = U(Y), Fk = Fk(Y) and Dk = Dk(Y). In this paper,Y will always refer to pressure primitive variables, so that:
r is the density
uk 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 cv is the specific heat at constant volume.
Also,kc is the kinetic energy per unit mass:
The advective and diffusive flux vectors are given by:
where are the components of the viscous tensor tv, qk is the kth component of the heat flux vector q and dkh 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:
and Kkh is such that:
Expressions for the matrices A0, Ak and Kkh can be found in the appendix.
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, u1 and u2). 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 Cp is the specific heat at constant pressure, ap 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 A0 in the case of liquids, when the necessary derivatives are not easily obtained from Eq.(20).
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 nth space-time slab partition can be represented as:
where is the spatial approximate partition at time level n and ln 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 tn. We then define ln 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 Shn and Vhn 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 Nk , normal to the spatial domain.
We approximate the solution vector in the usual manner, that is:
where Nn is the number of nodes in the spatial domain, jj 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 ´ Nn 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 Pi is the SUPG perturbation term, given by:
where t is the stabilizing SUPG matrix.
The shape function jj 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 dc 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 Ri and the matrix Mi can be assembled from their corresponding elementary counterparts Me and Re which hold the following pattern:
The elements of the matrix Me and of the vector Ve 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 dc are well defined for purely incompressible flows, when both ap 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 dc in Eq. by the augmented operator (1 -a)dc.
In Eq.(40), tin is suited for the incompressible limit whereas tcp is valid for compressible flows, so that the both extremes of the spectrum are covered.
Hauke (1995) used a full stabilizing matrix for both tin and tcp. We propose to use a very simple design for tcp and a diagonal pressure based din. 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 xk ® ek and the natural coordinates xk, we proceed to define the gradient of the vector ek 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 dc , slightly modified in order to accommodate only compressible flows:
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):
In equations (51) to (54), above, he is a characteristic dimension of the element. It may be taken as:
where Ae corresponds to the element's area. In the previous equations, we also have that c is the speed of sound.
In what follows, a GMRES (General Minimum Residue) solver with reverse communication was used to obtain the solution of Eq. , in the predictormulticorrector algorithm. The solver was obtained as a closed package from the public domain directory of CERFACS2 (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 Cr stands for the Courant-Friedrichs-Lewy number (CFL) and the coefficients fc and fd 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/m3, Cp = 1005.0J/(kgK), cv = 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/m2 and k = 2700.0W/mK whereas for Re=400, µ = 0.005Ns/m2 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 Cr = 1 (see Eq. (57)) for Re = 1 and Cr = 500for Re = 400. No convergence was achieved when we tried to use Cr = 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 u1 and u2 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): u1 = 1.0, u2 = 0.0 q = 2.769 x 10-4K and r = 1.0kg/m3
on C (-0.2 < x < 1.2 and y = 0.8): u1 = 1.0, u2 = 0.0, q = 2.769x10-4K and r = 1.0kg/m3
on D (x = 1.2 and 0.0 < y < 0.8): free flow
on E (0.0 < x 1.2 and y = 0.0): u1 = 0.0, u2 = 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.
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.
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:
Aliabadi, S. K. & Tezduyar, T. E., 1995, "Parallel fluid dynamics computations in aerospace applications", International Journal for Numerical Methods in Fluids, Vol 21, pp. 783-805. [ Links ]
Almeida, R. C. C., 1993, "Uma formulação de Petrov-Galerkin para a resolução das equações de Euler e Navier-Stokes 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", McGraw-Hill, 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 / Petrov-Galerkin formulations for advection dominated flows with particular emphasis on the inccompressible Navier-Stokes equations" Computer Methods in Applied Mechanics and Engineering, Vol 32, pp.199-259. [ Links ]
Codina, R., 1993, "A finite element formulation for the numerical solution of the advection-diffusion 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, "Petrov-Galerkin solutions of incompressible Navier-Stokes equations in primitive variables with adaptive remeshing" Computer Methods in Applied Mechanics and Engineering, Vol 106, pp 143-178. [ Links ]
Ghia, U., Ghia, K. N. & Shin, C., T., 1982, "High-re solutions for incompressible flow using Navier-Stokes equations and a multigrid method", Journal of Computational Physics, Vol 48, pp. 387-411. [ Links ]
Gustafsson, B. & Stoor, H., 1991, "Navier-Stokes equations for almost incompressible flow" SIAM Journal on Numerical Analysis, Vol 28 - 6, pp 1523-1547. [ 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 389-395. [ Links ]
Hauke, G., 1995, "A unified approach to compressible and incompressible flows and a new entropy-consistent formulation of the K-epsilon 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 advective-diffusive systems" Computer Methods in Applied Mechanics and Engineering, Vol 58, pp 305-328. [ 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 229-243. [ 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.183-199. [ Links ]
Shakib, F., 1988, "Finite Element analisys of the compressible Euler and Navier-Stokes 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 Navier-Stokes equations" Computer Methods in Applied Mechanics and Engineering, Vol 89, pp 141-219. [ Links ]
Schlichting, H., 1979, "Boundary Layer Theory", McGraw-Hill, Antigua, Guatemala, 817p. [ Links ]
Wark, J. S., 1995, "Advanced thermodynamics for engineers", McGraw-Hill, EUA, 564p. [ Links ]
Yoon, K. T., Moon, S. Y, Garcia, S.A., Heard, G. W. & Chung, T. J., 1998, "Flowfield dependent mixed explicit-implicit (FDMEI) methods for high and low speed and compressible and incompressible flows" Computer Methods in Applied Mechanics and Engineering, Vol 151, pp 75-104. [ 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.105-121. [ Links ]
Paper accepted July, 2005.
Technical Editor: Aristeu da Silveira Neto.