A method to numerically solve the Euler equations for fluids with general equations of state is presented. It is based on a formulation solving the conservation equations for either pressure primitive variables or entropy variables, instead of the commonly used conservation variables. We use a space-time discontinuous Galerkin finite-element discretization, which yields a highly local, potentially higher-order scheme. The algorithm is applied to test cases for compressible fluids to demonstrate its capabilities and the performance of the different variable sets

Pesch, Lars

Vegt, Jaap J.W. van der

Pesch, L.

van der Vegt, Jacobus J.W.

University of Twente Research Information

European Conference on Computational Fluid DynamicsECCOMAS CFD 2006P. Wesseling, E. Oñate and J. Périaux (Eds)c© TU Delft, The Netherlands, 2006A SPACE-TIME DISCONTINUOUS GALERKIN FINITE-ELEMENTDISCRETIZATION OF THE EULER EQUATIONS USING ENTROPYVARIABLESLars Pesch and Jaap J. W. van der VegtUniversity of Twente,Numerical Analysis and Computational Mechanics Group,Department of Applied Mathematics,P.O. Box 217, 7500 AE Enschede, The Netherlandse-mail: {L.Pesch, J.J.W.vanderVegt}@math.utwente.nlweb page: http://www.math.utwente.nl/nacm/Key words: Euler Equations, Discontinuous Galerkin, Compressible/Incompressible Fluids,Entropy Variables, Equation of StateAbstract. A method to numerically solve the Euler equations for fluids with general equa-tions of state is presented. It is based on a formulation solving the conservation equations foreither pressure primitive variables or entropy variables,instead of the commonly used con-servation variables. We use a space-time discontinuous Galerkin finite-element discretization,which yields a highly local, potentially higher-order schem . The algorithm is applied to testcases for compressible fluids to demonstrate its capabilities and the performance of the differentvariable sets.1 INTRODUCTIONMany numerical methods for fluid dynamics are suitable only for idealized types of fluids,like those for compressible flow being tailored to ideal gasses, and on the other hand methodsthat are based on formulations for (nearly) incompressiblemedia. The inclusion of complicatedequations of state, necessary to describe the physics of complex fluids, raises difficulties, as doesthe computation on domains where compressible and incompressibl flow conditions coexist.One way to derive numerical methods that are suitable both for compressible and incom-pressible flows is to solve for other variables than the set ofconserved variables, for which theequations are originally derived. Two sets of variables that yield well-posed formulations inthe incompressible limit are pressure primitive variablesand entropy variables1. The entropyvariables have the additional advantage that they symmetrize the equations and link to nonlin-ear stability2. Another key feature of the entropy variable formulation isits close relation tothe underlying thermodynamics, which can be incorporated by expressing the equations of statein terms of two material coefficients, the volume expansivity and isothermal compressibility,which are available in analytical or tabulated form for relevant substances.1Lars Pesch and Jaap J. W. van der VegtThe Euler equations are discretized with respect to the abovementioned sets of variables(entropy and pressure primitive) using a discontinuous Galerkin (DG) finite element method.It extends the space-time DG discretization discussed by van der Vegt and van der Ven3 togeneral equations of state and a formulation suitable for both c mpressible and incompressibleflows. The algorithm results in a highly local and potentially higher-order discretization. Sincespace-time finite-element basis functions are used, the method is well-suited for problems withmoving and deforming boundaries, and localhp-adaptation can be accommodated naturally.In this article, the algorithm will be demonstrated with test cases for media with ideal gasequations of state.2 GOVERNING EQUATIONS2.1 Euler equationsWe consider a fluid in the two-dimensional domainΩ during the time intervalT = [ts, te],described by the Euler equations of fluid dynamics written inco servation form, i.e.∂Ui∂t+∂Fir∂xr= 0 , i ∈ 1, . . . , 4 , ∀ (x, t) ∈ Ω × T , (1)combined with suitable boundary and initial conditions of the formU(x, t) = B(U,Uw) , ∀ (x, t) ∈ ∂Ω × T ,U(x, ts) = U0(x) , ∀ x ∈ Ω ,where the summation convention applies to the repeated index r for the space dimensions,r = 1, 2. The vector ofconservation variables Uand the flux matrixF are given byU =ρρv1ρv2ρetot, F =ρv1 ρv2ρv21 + p ρv1v2ρv2v1 ρv22 + pv1(ρetot + p) v2(ρetot + p). (2)Hereρ denotes density,v1, v2 are the velocity components with respect to a Cartesian coordinatesystem, andp is pressure. The total energyetot in this case is the sum of internal and kineticenergy,etot = e+ k = e+ 12v2r . (3)2.2 Equations of stateThe system (1) for the variables and flux (2) contains two moreunknowns than equations,hence it is not closed. One additional relationship has beengiven by Equation (3), but to com-plete the description, we have to specify how the pressurep d pends on the other variables.Remark: In the more general context of the Navier-Stokes equations,a additional quantity,the temperatureT, is needed to describe the diffusive effect of internal energy: heat conduction.2Lars Pesch and Jaap J. W. van der VegtTo remain within the framework of thermodynamics, which relat s state variables of divariantfluids (i.e. those whose state is determined by two independent variables, e.g.p andT), weinclude temperature in our formulation. Equations that diagnostically relate state variables of a fluid, for instanceρ, p, e,T, are calledequations of state. The system (1,2) can be closed by specifying, for example, the functionse= e(p,T) andp = p(ρ,T).For an ideal gas with constant specific heat (at constant volume),cv, and the gas constantR,these relations take the forme= e(T) = cv T , (4)p = ρR T . (5)A different idealization, the incompressible fluid, represents the limit in which there is only oneindependent thermodynamic variable as the densityρ is constant.To characterize the compressibility of a material we use twocoefficients, the volume expan-sivity αp and the isothermal compressibilityβT , which are defined as the relative changes of thespecific volumev = 1ρasαp ≔1v(∂v∂T)p, βT ≔ −1v(∂v∂p)T. (6)General equations of state can be formulated based on these two quantities as functions of thethermodynamic state4. For an ideal gas described by Equation (5) the coeffici nts evaluate toαp = 1/T andβT = 1/p; for incompressible fluids both coefficientsαp andβT are zero.2.3 Entropy variable formulationWe now interpret the conservation variablesU as dependent on some other set of variables. Inparticular we are interested in theprimitive variables including pressureandentropy variables,Y =pv1v2T, V =1Tµ − 12v2mv1v2−1, (7)respectively, with the chemical potentialµ.By rewriting (1) as quasi-linear system with respect to any suitable set of variables,W, weobtainAW0∂W∂t+ AWr∂W∂xr= 0 , (8)with AW0 =∂U∂W , AWj =∂F j∂W . Analysis shows that1:3Lars Pesch and Jaap J. W. van der Vegt• For conserved variablesU, the flux JacobiansAUj , ( j = 1, 2), are not well defined in theincompressible limitαp, βT → 0 (i.e. they contain entries that diverge or take the form00). This is connected to the fact that the equation for densityin (1) loses its prognosticcharacter forρ = const, and causes the breakdown of many numerical methods usingconservation variables when simulating (almost) incompressible fluids.• The JacobiansAYj , AVj , ( j = 0, 1, 2), are well defined in the incompressible limit. Thiscorresponds to the suitability of numerical methods based on these sets of variables tocompute both compressibleand incompressible flows (in particular many methods forincompressible flow are based on the pressure variables).• The matrixAV0 (for entropy variables) is symmetric and positive definite,in the incom-pressible case this weakens to semi-definiteness. The matricesAVj , ( j = 1, 2), are sym-metric. Discrete methods based on the entropy variables satisfy the entropy productioninequality without (additional) dissipative mechanisms.These properties led us to adapt a discontinuous Galerkin discret zation3 to use the pressureprimitive and in particular the entropy variables, in placeof conservation variables. Note thatthis method is not built on the quasi-linearized form (8), but discretizes the original conservationequations (1). Hence the method is conservative, even in thei compressible limit.3 DISCONTINUOUS GALERKIN DISCRETIZATIONTo apply a finite element discretization, the domainΩ ⊂ 2 is approximated by a suitablemeshT h with NE elements whose minimal element in-circle diameter ish. We extend thisspace tessellation into the time dimension and consider theproblem immediately inE ⊂ 3.The space-time elementsKe, e = 1, . . . ,NE, are related to the reference elementK̂ = [−1; 1]3by mappingsGe : K̂ → Ke.The basis functions are defined on the reference elementK̂ as tensor products of the one-dimensional monomial basisPn = {p : [−1; 1]→ , pk(ξ) = ξk, k = 0, . . . , n}, and we allow dif-ferent maximum ordersnt andns for the time and space dimensions. Hence the basis on thereference element isQnt ,ns(K̂) = {ϕ : K̂ →  | ϕ ∈ Pnt × (Pns)2}. Using these basis functions onthe reference element, the space of the test and trial functions for the finite element formulationis defined asVnt ,nsT h={W : T h × T→ 4 |WKe ◦Ge ∈ span{(Qnt ,ns(K̂))4} ∀e= 1, . . . ,NE}. (9)Functions in this space are representedp r elementasWKe(x, t) =(nt+1)(ns+1)2∑i=1Wi ϕi(G−1e (x, t)) , (10)with Wi ∈ 4 andϕi ∈ Qnt ,ns(K̂). Note that these functions may be discontinuous at the elementboundaries.4Lars Pesch and Jaap J. W. van der VegtAs we will shortly integrate Equation (1) over the space-time domainE, we rewrite it inspace-time divergence form as2∑r=0∂Fir∂xr= 0 , i ∈ 1, . . . , 4 , (11)with x0 = t and the flux in the time dimensionF· 0(U) = U.The DG discretization is obtained by the following steps: Wemultiply Equation (11) withan arbitrary test functionφ ∈ Vnt ,nsT hand chose the expansion of the variables,Wh, from the samespace, too. Integrating over each space-time elementKe and summing over all elements wereceive the weak formulation:Find aWh ∈ Vnt ,nsT h, such that for allφ ∈ Vnt ,nsT hthe following relation is satisfied:NE∑e=1−∫Ke∂φp∂xrFpr(U(Wh))dτ +∫∂Keφp F̂pr(U(W−h ),U(W+h )) nrds= 0 . (12)The boundary of each elementKe is denoted as∂Ke with ni the components of the exteriorpointing normal vector at this surface. Since the discrete functionWh is not necessarily contin-uous at the element faces we replace theF (U) in the integral over the element boundaries by aconsistent, symmetric numerical flux̂F (U−h ,U+h ). Herein, the limit of the discrete functionWhfrom the interior ofKe is calledW−h , whereas the exterior limit is denotedW+h ; the latter is eitherbased on theuh in a neighboring element or on the boundary data.The numerical fluxF̂ ensures that the discretization is conservative. We can expos theconservativity of the discretization and indicate the computational structure of the DG finite el-ement method by introducing the facesFi , i ∈ {1, . . . ,NF}, which are the objects of codimensionone that bound the elements. In Equation 12, the integrals over the boundaries of elements areevaluated twice for each internal faceFi. The two integrals differ only in the sign of the normalvectorn, so that we can combine the integrands asφLpnrF̂pr + φRp(−nr) · F̂pr = (φLp − φRp)nrF̂pr,where superscripts (·)L and (·)R denote values taken from the elements on the (arbitrarily de-fined) left and right side of the face, respectively. (Because the support of each basis function islimited to one element only one of the contributionsφL, φR will be nonzero.) Then we have−NE∑e=1∫Ke∂φp∂xrFpr(U(Wh))dτ +NF∑i=1∫Fi(φLp − φRp) nrF̂pr(U(WLh ),U(WRh ))ds= 0 . (13)This constitutes a nonlinear system of equations for the expansion coefficients ofWh in eachelement, which is solved using a pseudo-time method5.As the spatial numerical flux, we use the HLLC approximate Riemann solver for the idealgas test cases in Section 4; for other cases this may be replaced by the analytical flux of themean state on the two sides of the face (with an added stabiliztion term). On the faces whichseparate the time slabs, an upwind flux is used conforming to causality.5Lars Pesch and Jaap J. W. van der VegtThe treatment of boundary conditions through external state in the Riemann solver, whichis a bye-product of the DG discretization, considerably simplifies the specification of variousconditions, like sub- and supersonic in-/outflow or slip flow. In the spatially continuous dis-cretization of Hauke and Hughes1, the discrete equation system is solved by a Newton method,which requires not only a costly global system assembly but also the linearization of the bound-ary conditions.4 RESULTSTo demonstrate the possibilities arising from the formulation presented above, we give somecomputational results.4.1 Flow around a circular cylinderWe consider the subsonic flow of an ideal gas around a fixed circular cylinder with unitradius. The far-field flow is aligned with thex-axis, constant, and homogeneous:v∞ = u∞ex. Inthe limit of vanishing Mach number,M∞ → 0, the flow field is given by the classical potentialflow solution. To conform to the study by van der Vegt and van der Ven6, we chose the free-stream Mach number asM∞ = 0.38. At the inflow, we use far-field boundary conditions for themass flux, stagnation enthalpyhtot = h+ 12v2r , and pressure:ρRvR = ρ∞v∞, htotR = htot∞ , pR = pL,respectively, with (·)L denoting the flow state on the inside of the flow domain, while (·)R is theboundary data. Corresponding outflow conditions areρRvR = ρLvL, sR = sL, pR = p∞, with thespecific enthalpys. Due to the absence of viscosity in the Euler equations, a slip flow conditionis applied on the cylinder surface:v · n = 0. This is enforced by setting the momentum in theboundary data for the HLLC flux toρRvR = ρLvL − 2ρL(vL · n)n.We indicate the error of the numerical method using the totalpressure lossπ, defined asπ = 1−pp∞1+ 12(γ − 1)M21+ 12(γ − 1)M2∞γγ−1. (14)The pressure loss should be zero everywhere, deviations occur nly due to the discrete ap-proximation, which leads to entropy production at the surface of the cylinder. The resultingerror according to a computation with conservation variables on a mesh with 48× 32 elementsis plotted in Figure 1. Note that in all results we plot thediscontinuoussolution without anaveraging per element or per node. Figure 2 shows the total pressu e loss along the cylindersurface for entropy variables using two different orders of the spatial representation (since thisis a steady state case, we usent = 0 throughout). Results for different variable sets are hardlydistinguishable in this test, so we do not present a comparison.4.2 Oblique shockWe consider an inviscid supersonic flow atM = 2 on the unit square. The inflow is at an angleof 10◦ with the lower boundary of the domain, which is a wall (with slip condition). Conse-quently a shock forms at an angle of 29.3◦, cf. Hauke and Hughes1. The density field including6Lars Pesch and Jaap J. W. van der VegtFigure 1: Total pressure loss close to the cylinder.Computed on a 48× 32 mesh using conservation vari-ablesU.Figure 2: Comparison of the total pressure loss com-puted usingV-variables on a 48× 32. Dashed line:ns = 0 (constant representation ofVh per element),solid line:ns = 1 (bilinear representation ofVh).the shock is plotted in Figure 3. Note that we are not using anyform of discontinuity capturingor stabilization, hence the over- and undershoot close to the shock. In this way it is also rea-sonable to compare the results for different sets of variables, cf. Figure 4. Minor differencescan be observed. Presumably these arise due to amplificationby the different transformationsbetween the pressure primitive or entropy variables and theplotted variableρ. On the otherhand, the overall location of the shock is correctly represented, in contrast to formulations thatuse non-conservative formulations.5 CONCLUSIONS & OUTLOOKWe have described the extension of a discontinuous Galerkinfinite element discretizationof the Euler equations to allow using different sets of variables. In particular we have consid-ered entropy variables, because of their theoretically andpractically appealing properties. Thischoice allows to treat compressible and incompressible flows with the same numerical method.We have presented numerical examples for the compressible case.The extension of the numerical method to include the diffus ve terms of the Navier-Stokesequations is under way using the discretization described by Klaij et al.7. The goal is to combinethe presented scheme with an interface tracking method8 t compute multiphase flows withcompressible and incompressible media.References[1] G. Hauke and T. J. R. Hughes. A comparative study of different sets of variables for solvingcompressible and incompressible flows.Comput. Methods Appl. Mech. Eng., 153:1–44,7Lars Pesch and Jaap J. W. van der VegtFigure 3: Oblique shock test case; density field com-puted on a 20× 20 mesh using entropy variables.Figure 4: Comparison of the computed density fieldsalong x = 0.9. Solid line:U, dashed:Y, dotted:V-variables. All computations carried out on a 20× 20mesh withns = 1 and no stabilization added to exposethe differences between the different variable sets.1998.[2] G. Hauke, A. Landaberea, I. Garmendia, and J. Canales. Onthe thermodynamics, stabilityand hierarchy of entropy functions in fluid flow.Comput. Methods Appl. Mech. Eng., 195(33–36):4473–4489, 2006. doi: 10.1016/j.cma.2005.09.010.[3] J. J. W. van der Vegt and H. van der Ven. Space-time discontinuous Galerkin finite elementmethod with dynamic grid motion for inviscid compressible flows: I. General formulation.J. of Comput. Phys., 182(2):546–585, 2002.[4] R. L. Panton.Incompressible Flow. John Wiley & Sons, Inc., 1984.[5] C. M. Klaij, J. J. W. van der Vegt, and H. van der Ven. Pseudo-time stepping methodsfor space-time discontinuous Galerkin discretizations ofthe compressible Navier-Stokesequations.J. Comput. Phys. (accepted), 2006.[6] J. J. W. van der Vegt and H. van der Ven. Slip flow boundary conditions in discontinu-ous Galerkin discretizations of the Euler equations of gas dynamics. In H. A. Mang, F. G.Rammerstorfer, and J. Eberhardsteiner, editors,Proceedings Fifth World Congress on Com-putational Mechanics, pages 1–16, July 2002.[7] C. M. Klaij, J. J. W. van der Vegt, and H. van der Ven. Space-time discontinuous Galerkinmethod for the compressible Navier-Stokes equations.J. Comput. Phys., accepted, 2006.[8] W. E. H. Sollie, J. J. W. van der Vegt, and O. Bokhove. A space-time discontinuous Galerkinfinite element method for two-fluid problems.To be submitted to J. Comput. Phys., 2006.8

English

A space-time discontinuous Galerkin finite-element discretization of the Euler equations using entropy variables

http://doc.utwente.nl/65616/1/CFD06paper.submitted.pdf

A space-time discontinuous Galerkin finite-element discretization of the Euler equations using entropy variables

Abstract

Similar works

Full text

Available Versions

University of Twente Research Information