An HDG Method for unsteady compressible flows Alexander Jaust, Jochen Sch¨ utz and Michael Woopen Preprint No. 410 Key words: October 2014 CFD, hybridized discontinuous Galerkin, shock-capturing AMS Subject Classifications: 65M60, 35Q31, 35L67 Institut f¨ ur Geometrie und Praktische Mathematik RWTH Aachen Templergraben 55, D-52056 Aachen, (Germany) Alexander Jaust MathCCES, RWTH Aachen University, Schinkelstrae 2, 52062 Aachen, e-mail: [email protected] Jochen Sch¨ utz IGPM, RWTH Aachen University, Schinkelstrae 2, 52062 Aachen, e-mail: [email protected] Michael Woopen AICES, RWTH Aachen University, Schinkelstrae 2, 52062 Aachen, e-mail: [email protected] An HDG Method for unsteady compressible flows Alexander Jaust and Jochen Sch¨utz and Michael Woopen Abstract The recent gain of interest in Discontinuous Galerkin (DG) methods shows their success in computational fluid dynamics. However, in the case of an implicit discretization, one potential drawback is the high number of globally coupled unknowns. By means of hybridization, this number can be significantly reduced. The hybridized DG (HDG) method has proven to be beneficial especially for steady flows. In this work we apply this method to time-dependent flow problem with shocks. Due to its inherently implicit structure, time integration methods such as diagonally implicit Runge-Kutta (DIRK) methods present themselves as natural candidates. Furthermore, as the application of flux limiting to HDG is not straightforward, an artificial viscosity model is applied to stabilize the method. Key words: CFD, hybridized discontinuous Galerkin, shock-capturing 1 Introduction A prominent class of high-order methods for computational fluid dynamics are socalled discontinuous Galerkin methods [3, 4, 11, 10, 9, 7]. Based on a partitioning of the domain into a set of N elements, the solution is approximated by piecewise polynomials on each of the elements. This allows local (non-conforming) refineAlexander Jaust MathCCES, RWTH Aachen [email protected] University, Schinkelstraße 2, 52062 Aachen, e-mail: Jochen Sch¨utz IGPM, RWTH Aachen University, Templergraben 55, 52056 Aachen, e-mail: [email protected] Michael Woopen AICES, RWTH Aachen University, Schinkelstraße 2, 52062 Aachen, e-mail: [email protected] 1 2 Alexander Jaust and Jochen Sch¨utz and Michael Woopen ment by varying the number of elements or the degree p of the polynomials used on each of the elements. However, these methods suffer from a large number of globally coupled unknowns when being used with implicit time-stepping methods or in the context of stationary problems. An approach to reduce this number is to use hybridization resulting in the hybridized discontinuous Galerkin (HDG) method [14, 17, 15, 16, 21, 18]. By introducing an additional unknown λh living on the element interfaces, i.e. the mesh skeleton, the system of equations can be formulated such that it is only globally coupled in this hybrid variable. Therefore, the num ber of globally coupled unknowns can be changed asymptotically from O pd N to O pd−1 Nˆ where Nˆ is the number of element interfaces in the mesh. In this paper we focus on supersonic inviscid flows. These flows can be described by the compressible Euler equations, and it is known that those solutions can develop discontinuities. These discontinuities are a severe issue for high-order methods, as they usually show oscillatory behavior that leads stability issues. In order to capture discontinuities and stabilize computations we use an artificial shockcapturing method that has been introduced by Persson and Peraire [19]. 2 Numerical method In this section, we give a brief introduction to the hybridized discontinuous Galerkin method. For more details, we refer to, e.g., [14, 12]. We shortly describe the applied time integration and shock-capturing schemes. Please note, that we focus on the two dimensional case in this work, i.e., d = 2. Therefore, we will refer to element interfaces as edges from here on. 2.1 Governing Equations We consider supersonic inviscid flows that can be described using the compressible Euler equations, given by ∂w + ∇ · f (w) = 0 ∂t ∀ (x,t) ∈ Ω × [0, ∞) (1) given on a domain Ω , equipped with appropriate boundary conditions. The vector of conserved variables is w = (ρ, ρu1 , ρu2 , E)T and involves the density ρ, velocities u1 and u2 and total energy E. Convective fluxes f = ( f1 , f2 ) are given by f1 = (ρu1 , P + ρu21 , ρu1 u2 , u1 (E + P))T , f2 = (ρu2 , ρu1 u2 , P + ρu22 , u2 (E + P))T . (2) An HDG Method for unsteady compressible flows 3 Pressure P is determined using the ideal gas law with the adiabatic constant γ = 1.4 for air. 2.2 The Hybridized Discontinuous Galerkin Method In order to discretize equation (1), we assume a partitioning of Ω as Ω = Additionally, we define the following ansatz spaces k=1 Ω k . SN Vh := { f ∈ L2 (Ω ) | f|Ωk ∈ Π p (Ωk ) ∀k = 1, . . . , N}4 (3) b ek ∈ Γ }4 , Mh := { f ∈ L2 (Γ ) | f|ek ∈ Π p (ek ) ∀k = 1, . . . , N, (4) where Γ is the set of all edges. Given these definitions, the semi-discrete formulation of the HDG method can be written as N ∑ k=1 ((wh )t , ϕh )Ωk − ( f (wh ), ∇ϕh )Ωk + h fˆ · n, ϕh− i∂ Ωk = 0 ∀ϕh ∈ Vh (5) ∀µh ∈ Mh . (6) b N ∑ hJ fˆ · nK, µh iΓ = 0 k=1 (·, ·) and h·, ·i denote element and edge scalar product, respectively. Equation (5) is very similar to the weak formulation one obtains for standard DG methods. However, the coupling between cells is done via the hybrid variable λh . The third equation is derived from the original problem and is necessary to both determine λh and ensure that the total flux has a weak divergence. Note that the fluxes over edges have been substituted by the numerical fluxes fˆ := f (λh ) − α(λh − w− h )n, (7) with positive real parameter α. This is a Rusanov-/Lax-Friedrichs inspired flux that retains the locality of the scheme. In this setting, locality means that the system of equations on each element only depends on the information given on the element and on its edges. This allows the application of static condensation [6]. Therefore, the linearized system of equations can be written such that it is globally coupled only in λh . This may lead to an extensive reduction of globally coupled unknowns. 2.3 Time Integration In equations (5)-(6), only the temporal derivative of wh occurs. Thus, the discretization leads to a system of differential-algebraic equations (DAEs) of index 1. This poses a severe restriction to the time integration methods that can be applied [8]. Therefore, only implicit methods with adequate stability properties can be used with 4 Alexander Jaust and Jochen Sch¨utz and Michael Woopen HDG [15, 16, 17]. In our case, we found backward differentiation (BDF) methods and diagonally implicit Runge-Kutta (DIRK) methods [2, 1, 5, 8] to be reliable time integrators for the HDG method [12, 20, 13]. 2.4 Shock-Capturing High-order methods such as the DG or the HDG method tend to show oscillatory behavior near discontinuities or steep gradients if they are not suitably stabilized. For DG, flux limiters and artificial viscosity models are a generic choice. However, the HDG method relies on the locality of the method, which makes it hard or even impossible to employ flux limiters. The use of artificial viscosity models, however, is possible. In this work, we use a shock-capturing method introduced by Persson and Peraire [19], which relies on adding an additional diffusive term to the equations. In our case, we add a (inconsistent) discretization of the Laplacian, (εk ∇wh , ∇φh ), to equation (5). In comparison to the viscous terms of the Navier-Stokes equations, this is much simpler to evaluate. The cell-wise constant viscosity εk is determined using a shock sensor ! (w¯ − w, w¯ − w)L2 (Ωk ) (8) sk := log10 (w, w)L2 (Ωk ) where w¯ represents the L2 -projection of w from Π p (Ωk ) to Π p−1 (Ωk ). Please note, that this projection step is very cheap as we employ orthogonal basis functions. Hence, sk measures the discrete smoothness of a solution by comparing its highest order terms to the complete solution. Then, the actual amount of viscosity ek for each element is computed from 0 , sk < s0 − κav ε0 π(sk − s0 ) εk = 1 + sin , s0 − κav ≤ sk ≤ s0 + κav . (9) 2 2κav ε0 , sk > s0 + κav where ε0 ∼ hp , s0 ∼ log(p) and κav are problem-dependent parameters. They have to be chosen such that the method is sufficiently stabilized while keeping discontinuities sharp. 3 Numerical Results In this section we present results obtained from the described method. At each time step, we employ a Newton-Krylov solver based on a restarted GMRES for the re- An HDG Method for unsteady compressible flows 5 sulting linear system. As a preconditioner we use an incomplete LU factorization without additional levels of fill. 3.1 Double Mach Reflection at a Wedge This famous test case is taken from the paper of Woodward and Colella [22]. Supersonic flow enters the domain from the left and hits a wedge of 30 degree angle (see Fig. 1(a)). Instead of a rectangular domain as in the original paper, we choose the domain such that the flow enters the domain on the left and is parallel to the x1 -direction. Therefore, we do not have to prescribe the shock on the upper boundary conditions. We have inflow boundary conditions on the left, slip-wall boundary conditions at the wedge and symmetry boundary conditions everywhere else. The domain is initialized with pre-shock values, (ρL = 8.0, u1,L = 8.25, u2,R = 0, PL = 116.5)T , in front of the wedge, denoted by L, and post-shock values, (ρR = 1.4, u1,R = 0, u2,R = 0, PR = 1.0)T , everywhere else (see Fig. 1(b)). We run the simulation up to t = 0.2 and use a DIRK method of second order with 2 stages described by Alexander [2]. Due to limited deflection angles for the given flow conditions, the shock is reflected such that a complex structure occurs. A convex shock beginning at the tip of the wedge is created and a region with interacting shocks develops close to where the shock hits the wedge. We run two simulations with N1 = 3167 and N2 = 8395 elements. In both cases polynomials of degree p = 3 are employed. For the coarse mesh, see Fig. 2, and the fine mesh, see Fig. 3, we show the mesh, the artificial viscosity, the density and isolines of the density. For both meshes, the general structure of the solution is captured while on the finer mesh the shocks are sharper. The shocks are detected by the smoothness sensor such that artificial viscosity is only applied in these regions. In the region where the shock interacts with an occurring jet much less viscosity is added than at the outer shocks. However, on neither of the meshes Kelvin-Helmholtz instabilities in the jet can be seen. There is possibly too much dissipation due to the applied viscosity, the mesh resolution or the applied time integration method. Note that there are small oscillations in the area in front of the shock that cannot be seen in this figures. These may be reduced by using other parameters for the shock-capturing. However, these are actually that small, such that the shock sensor hardly recognizes them. 4 Conclusion and Outlook We have presented a hybridized DG method for an unsteady compressible flow problem. By applying an artificial viscosity model we can stabilize the method successfully to approximate flows with shocks. 6 Alexander Jaust and Jochen Sch¨utz and Michael Woopen L 30◦ (a) Sketch of the physical domain. R 30◦ (b) Partitioning of the domain for setting the initial conditions. Fig. 1 Sketches of computational domain and how initial data is distributed. (a) Mesh with 3167 elements. (b) Cellwise artificial viscosity using 14 levels from εk,min = 0.0 and εk,max = 0.014 at t = 0.2. (c) Density distribution for 20 levels with (d) Isolines of the density distribution for 20 ρmin = 1 and ρmax = 19 at t = 0.2. levels with ρmin = 1 and ρmax = 19 at t = 0.2. Fig. 2 Mesh and solution at t = 0.2. Future work will include more detailed studies regarding the optimal parameters for the shock-capturing scheme as well as the behavior on refined grids. The meshes used here rely on uniformly large elements which is not suitable for flows with sharp flow features such as shocks. Therefore, local adaptation of the polynomial degree p and the mesh is work in progress. The latter is a challenging task for unsteady flows since the flow features are likely to move. This introduces the need for both mesh refinement and coarsening to keep the number of elements low. An HDG Method for unsteady compressible flows (a) Mesh with 8395 elements. 7 (b) Cellwise artificial viscosity using 15 levels from εk,min = 0.0 and εk,max = 0.007 at t = 0.2. (c) Density distribution for 20 levels with (d) Isolines of the density distribution for 20 ρmin = 1 and ρmax = 19 at t = 0.2. levels with ρmin = 1 and ρmax = 19 at t = 0.2. Fig. 3 Density distribution for 20 levels with ρmin = 1 and ρmax = 19 at t = 0.2. In addition to adaptation, we also plan to try further shock-capturing strategies to see whether there are methods in particular well-suited for the HDG method. This also includes different ways of applying the artificial viscosity. As suggested in the work by Persson and Peraire [19] it may be beneficial to use the physical diffusive terms present in the Navier-Stokes equations instead of a Laplacian. Another point that to address are time integration schemes. So far, classical one step and multistep methods have been applied, but it may be worthwhile to also consider other concepts such as general linear methods. At the end, we expect to have a versatile numerical method for unsteady flow applications once all these points have been worked out. References 1. A. H. Al-Rabeh. Embedded DIRK methods for the numerical integration of stiff systems of ODEs. International Journal for Computer Mathematicss, 21:1:65–84, 1987. 2. R. Alexander. Diagonally Implicit Runge-Kutta Methods for Stiff O.D.E.’s. SIAM Journal on Numerical Analysis, 14(14):1006–1021, 1977. 3. D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of Discontinuous Galerkin methods for elliptic problems. SIAM J. Num. Anal., 39(5):1749–1779, 2002. 8 Alexander Jaust and Jochen Sch¨utz and Michael Woopen 4. F. Bassi and S. Rebay. A high-order accurate discontinuous finite-element method for the numerical solution of the compressible navier-stokes equations. J. Comp. Phys., 131:267– 279, 1997. 5. J. R. Cash. Diagonally Implicit Runge-Kutta Formulae with Error Estimates. Journal of the Institute of Mathematics and its Applications, 24:293–301, 1979. 6. Bernardo Cockburn, Jayadeep Gopalakrishnan, and Raytcho Lazarov. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 01 2009. 7. Krzysztof J. Fidkowski, Todd A. Oliver, James Lu, and David L. Darmofal. p–Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier–Stokes equations. J. Comp. Phys., 207:92–113, 2005. 8. E. Hairer and G. Wanner. Solving Ordinary Differential Equations II, Stiff and DifferentialAlgebraic Problems. Springer Series in Computational Mathematics. Springer-Verlag, Berlin Heidelberg, 1991. 9. R. Hartmann and P. Houston. Symmetric interior penalty DG methods for the compressible Navier–Stokes equations II: Goal–oriented a posteriori error estimation. International Journal of Numerical Analysis and Modeling, 3(2):141–162, 2006. 10. Ralf Hartmann and Paul Houston. Symmetric interior penalty DG methods for the compressible Navier–Stokes equations I: Method formulation. International Journal of Numerical Analysis and Modeling, 3(1):1–20, 2006. 11. Paul Houston and Endre S¨uli. hp-adaptive discontinuous galerkin finite element methods for first order hyperbolic problems. SIAM J. Sci. Comput., 23(4):1226–1252, 2001. 12. A. Jaust and J. Sch¨utz. A temporally adaptive hybridized discontinuous Galerkin method for instationary compressible flows. Computers and Fluids, 98:177–185, 2014. 13. Alexander Jaust, Jochen Sch¨utz, and Michael Woopen. A hybridized discontinuous galerkin method for unsteady flows with shock-capturing. AIAA Paper 2014-2781, 44th AIAA Fluid Dynamics Conference, 2014. 14. N. C. Nguyen and J. Peraire. Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics. Journal of Computational Physics, 231:5955–5988, 2012. 15. N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations. Journal of Computational Physics, 228(9):3232–3254, 5 2009. 16. N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit high-order hybridizable discontinuous Galerkin method for nonlinear convection-diffusion equations. Journal of Computational Physics, 228(23):8841–8855, 12 2009. 17. N. C. Nguyen, J. Peraire, and B. Cockburn. High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics. Journal of Computational Physics, 230:36953718, 2011. 18. J. Peraire, N. C. Nguyen, and Bernardo Cockburn. A Hybridizable Discontinuous Galerkin Method for the Compressible Euler and Navier-Stokes Equations. AIAA Paper 2010-362, 48th AIAA Aerospace Sciences Meeting and Exhibit, 2010. 19. P.-O. Persson and J. Peraire. Sub-cell shock capturing for Discontinuous Galerkin methods. AIAA Paper 06-0112, American Institute of Aeronautics and Astronautics, 2006. 20. J. Sch¨utz, M. Woopen, and G. May. A combined hybridized discontinuous Galerkin / hybrid mixed method for viscous conservation laws. In F. Ancona, A. Bressan, P. Marcati, and A. Marson, editors, Hyperbolic Problems: Theory, Numerics, Applications, pages 915–922. American Institute of Mathematical Sciences, 2012. 21. Jochen Sch¨utz and Georg May. A Hybrid Mixed Method for the Compressible Navier-Stokes Equations. Journal of Computational Physics, 240:58–75, 2013. 22. P. Woodward and P. Colella. The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics, 54:115–173, 1984.

© Copyright 2017 ExploreDoc