PaperlayerView exact revision on arXiv

Source-linked AI summary

High order exactly divergence-free Hybrid Discontinuous Galerkin Methods for unsteady incompressible flows

Christoph Lehrenfeld, Joachim Schöberl

arXiv:1508.04245v2math.NAcs.CE

TL;DR

The paper addresses efficient discretization of unsteady incompressible Navier–Stokes equations. It combines operator splitting with Upwind DG for transport and H(div)-conforming HDG for Stokes problems, including projected jumps; benchmark cases demonstrate the method’s performance.

  • Problem

    Efficiently discretizing unsteady incompressible Navier–Stokes equations requires treating nonlinear transport and stiff linear Stokes terms with suitable temporal and spatial methods.

  • Method

    The method combines explicit Upwind DG treatment of hyperbolic transport with H(div)-conforming HDG for Stokes problems and a projected jumps operator that reduces global coupling.

  • Results

    Numerical test cases demonstrate the performance of the combined DG/HDG discretization for incompressible Navier–Stokes equations.

  • Takeaways & Limitations

    The combined formulation is tailored for efficiency by pairing standard DG transport discretization with an efficiently solvable, exactly divergence-free HDG Stokes discretization.

  • Takeaways & Limitations

    The explicit and implicit parts of the additive operator-splitting decomposition must use the same time-step size, so convection stability restricts implicit solution-step counts.

Abstract

from arXiv · show

In this paper we present an efficient discretization method for the solution of the unsteady incompressible Navier-Stokes equations based on a high order (Hybrid) Discontinuous Galerkin formulation. The crucial component for the efficiency of the discretization method is the disctinction between stiff linear parts and less stiff non-linear parts with respect to their temporal and spatial treatment. Exploiting the flexibility of operator-splitting time integration schemes we combine two spatial discretizations which are tailored for two simpler sub-problems: a corresponding hyperbolic transport problem and an unsteady Stokes problem. For the hyperbolic transport problem a spatial discretization with an Upwind Discontinuous Galerkin method and an explicit treatment in the time integration scheme is rather natural and allows for an efficient implementation. The treatment of the Stokes part involves the solution of linear systems. In this case a discretization with Hybrid Discontinuous Galerkin methods is better suited. We consider such a discretization for the Stokes part with two important features: H(div)-conforming finite elements to garantuee exactly divergence-free velocity solutions and a projection operator which reduces the number of globally coupled unknowns. We present the method, discuss implementational aspects and demonstrate the performance on two and three dimensional benchmark problems.

1. Introduction

The paper introduces an efficient high-order DG/HDG discretization for unsteady incompressible Navier–Stokes equations by pairing operator splitting with spatial methods tailored to transport and Stokes subproblems.

  • Problem and approach: Upwind DG provides a natural stable discretization for convection-dominated transport problems.The convection operator is treated explicitly through operator evaluations, avoiding repeated nonlinear linear-system setups.
  • Problem and approach: For the Stokes part, H(div)-conforming HDG yields exactly divergence-free velocity solutions, while hybridization reduces globally coupled unknowns and enables static condensation.These properties target efficient implicit solution of the linear Stokes systems.
  • Contributions: The paper introduces a projected jumps formulation, also called reduced stabilization or reduced-order HDG, for high-order H(div)-conforming Stokes discretization.The projection operator is presented as a general implementation of an approach previously limited to two dimensions with special integration rules.
  • Evaluation: The paper discusses efficiency and practicability through implementational aspects and numerical performance studies.The stated evaluations include two- and three-dimensional benchmark problems.
  • Problem and approach: The method decomposes Navier–Stokes into an unsteady Stokes problem and a hyperbolic transport problem using operator-splitting time integration.Different finite element spaces are used for the different spatial operators.
  • Contributions: The combined spatial discretization uses standard Upwind DG for hyperbolic transport and the new HDG method for Stokes-Brinkman problems.Transfer operations connect the different discretization spaces.

2. DG/HDG spatial discretization

The spatial discretization separates the Stokes, convection, and transfer operators, using HDG for the Stokes part and standard DG for convection.

  • Spatial discretization: The spatial formulation introduces discretizations for the Stokes operator, convection operator, and transfer operations between their spaces.These components are combined to formulate the semi-discrete problem.
  • Stokes discretization: The Stokes-Brinkman problem uses an H(div)-conforming Hybrid DG discretization improved by a projected-jumps modification.Projected jumps are also termed reduced stabilization or reduced-order HDG.
  • Convection and transfer: The convection part uses standard DG approaches, while transfer operations connect the different Stokes and convection discretization spaces.The resulting components support the final semi-discrete formulation.

2.1. H(div)-conforming HDG formulation of the Stokes-Brinkman problem

The Stokes-Brinkman discretization combines H(div)-conforming velocity-pressure spaces with tangential HDG hybridization. This structure yields exactly divergence-free velocities, stable viscous discretization, reduced coupling, and optimal-order error estimates.

  • 2.1. H(div)-conforming HDG formulation of the Stokes-Brinkman problem: The Stokes-Brinkman problem adds an inertia reaction term τ^-1 to the stationary incompressible Stokes equations.The resulting problem is treated as a Stokes-Brinkman problem within the HDG formulation.
  • 2.1.1. H(div)-conforming Finite Elements for Stokes.: H(div)-conforming velocity functions enforce continuity of normal components across element interfaces, while tangential continuity is imposed weakly.The velocity space uses element-wise polynomial functions with continuous normal traces; a DG formulation handles tangential components.
  • 2.1.1. H(div)-conforming Finite Elements for Stokes.: The chosen velocity-pressure pair converts weak incompressibility into strong incompressibility, so discrete velocities are exactly divergence-free.The equivalence is expressed by div(uT)q dx = 0 for all qh in Qh if and only if div(uT) = 0 in Ω.
  • 2.1.1. H(div)-conforming Finite Elements for Stokes.: Exact divergence-freeness supports energy stability and permits velocity error estimates independent of the pressure field.The paper identifies both consequences directly for the incompressible Navier-Stokes discretization.
  • 2.1.2. The HDG space for the velocity.: The HDG velocity space introduces facet unknowns for tangential traces, coupling neighboring elements only through these unknowns.Facet variables are introduced to improve linear-system handling; they contribute here only to viscous-force discretization.
  • 2.1.3. Viscous forces.: The viscous bilinear form assigns consistency, adjoint consistency, and stability roles to its terms, and is consistent, bounded, and coercive for sufficiently large α.The latter properties are stated with respect to a modified discrete norm.
  • 2.1.3. Viscous forces.: Static condensation preserves local element-or-aligned-facet sparsity while reducing the globally coupled velocity unknowns to skeleton variables.Element unknowns couple only within the same element or with aligned facets.
  • 2.1.5. Pressure force and incompressibility constraint: The discretization satisfies an h-independent LBB condition, and Brezzi’s theorem yields optimal-order a priori error estimates.The result follows from coercivity, the LBB condition, and consistency and continuity of the bilinear forms.

2.2. Projected jumps: An enhancement of the HDG Stokes discretization

The projected jumps formulation reduces facet polynomial degree while retaining element degree k, preserving consistency, boundedness, and coercivity. Local elimination implements the projection implicitly and reduces globally coupled unknowns.

  • Projected jumps: The method.: The projected jumps formulation overcomes the sub-optimality of using degree k for both primal and facet unknowns.A related hybrid mixed construction can postprocess degree-k unknowns to degree k + 1, whereas the projection improves the equal-degree HDG approach.
  • Projected jumps: The method.: Facet unknowns use degree k −1 while element unknowns retain degree k, reducing the facet-space polynomial degree.The reduction is introduced through a projection operator and projected jumps formulation.
  • Projected jumps: The method.: Only the final boundary integral is changed; the modification preserves consistency and boundedness and permits coercivity in a weaker norm.The first two boundary integrals are reformulated rather than substantively modified.
  • Projected jumps: The method.: The unchanged element space keeps the normal component of the element velocity polynomial at degree k.The reduction concerns the facet representation, not the element space Wh.
  • Interplay with operator-splitting: The reduction is appropriate within the paper’s operator splitting because facet variables are excluded from implicit convection treatment, even for convection-dominated problems.The stated premise is that convection does not enter implicit solution steps.
  • Implementation: The highest-order facet component is duplicated element-locally and eliminated, leaving only element unknowns and reduced facet unknowns in the global system.A Schur complement removes the eliminated degrees of freedom from the final stiffness matrix.

2.3. DG formulation for the convection

The convection component uses a standard Upwind DG discretization on elementwise polynomial spaces. An L2-orthogonal basis diagonalizes the mass matrix, while exact divergence-free convection velocity yields stability.

  • DG formulation for the convection: The convection discretization uses the standard DG space of vector-valued elementwise polynomials of degree k.The space is defined by Vh := {u : u ∈ [Pk(T)]d for every element T}.
  • DG formulation for the convection: An L2-orthogonal element basis makes the associated mass matrix diagonal.This supports efficient mass-operation evaluation.
  • DG formulation for the convection: The standard Upwind DG trilinear form provides a stable spatial discretization when the convection velocity is exactly divergence-free.The stability statement includes the condition div(uT) = 0.

2.4. Transfer operations and embeddings

Transfer operations connect the distinct HDG and DG spaces used for Stokes and convection. For affine elements, the embedding uses mass matrices efficiently, but the allowed transfers restrict the time-integration scheme.

  • Transfer operations and embeddings: The convection and Stokes operators act on different spaces, so transfer operations based on embeddings make the discretizations compatible.The combined formulation uses these operations to apply convection and mass operators in the HDG space.
  • Embeddings: With Wh ⊂ Vh, the canonical embedding maps HDG functions to the DG convection space, while its coefficient matrix is not an identity.The dual embedding maps in the reverse functional direction.
  • Embeddings: The coefficient embedding is I = (MV)^−1MU,V, obtained from an equivalent L2 problem.Here MV is the DG mass matrix and MU,V is the mixed mass matrix.
  • Implementation: For affine linear transformations, MV is diagonal and transfer application costs essentially one sparse matrix multiplication.Thus applying the inverse mass matrix remains efficient.
  • Restriction on time integration scheme: The transfer restriction forbids applying Stokes operators to DG functions and using HDG functionals in convection solution steps.This imposes a restriction on the time-integration scheme.
  • Curved elements: On curved elements, the embedding becomes an L2 projection and MV has non-diagonal blocks only on curved elements, so inverse-mass applications remain cheap.The usual Piola transform means Wh is no longer a subspace of Vh.

2.5. The semidiscrete formulation

The semidiscrete problem combines the HDG Stokes bilinear form, DG convection operator, pressure coupling, forcing, incompressibility, and an initial condition. The formulation preserves kinetic-energy stability through the convection discretization.

  • The semidiscrete formulation: The spatially discrete problem seeks velocity u(t) in Uh and pressure p(t) in Qh.The formulation is posed as a differential-algebraic problem with incompressibility constraints.
  • The semidiscrete formulation: The momentum equation combines the mass derivative, viscous form Ah, convection form Ch, pressure coupling Dh, and forcing.The equation is enforced for every test velocity v and time t in [0,T].
  • The semidiscrete formulation: The pressure equation imposes Dh(u,q) = 0 for every pressure test function q, enforcing the discrete incompressibility condition.The initial velocity is specified through the mass bilinear form Mh.
  • The semidiscrete formulation: The convection stability property yields stability of the kinetic energy in the semidiscrete formulation.The convection form is defined on Uh using the embedding from the HDG space to the DG space.
  • The semidiscrete formulation: Operator-splitting time-integration methods are introduced next to solve the semidiscrete system efficiently.The semidiscrete formulation is therefore the starting point for the temporal treatment.

3. Operator-splitting time integration

The paper separates convection from Stokes terms in time, using operator-splitting schemes to treat convection explicitly while preserving incompressibility through stiffly accurate Stokes solves. Product and modified additive decompositions relax shared time-step restrictions, but product splitting introduces consistency error.

  • The semi-discrete system separates mass, Stokes, convection, and pressure operators while enforcing the discrete incompressibility constraint.
  • Additive decomposition using IMEX schemes: IMEX schemes treat convection explicitly and the remaining operators implicitly, but both parts must use the same time-step size.The explicit convection stability restriction can therefore increase the number of expensive implicit solution steps.
  • Product decomposition methods: Product decomposition solves separated Stokes and convection subproblems successively, allowing their time steps to be chosen independently.Its main disadvantage is an additional consistency error.
  • Product decomposition methods: The propagation operator advances a linear hyperbolic convection problem using extrapolated divergence-free velocities and a numerical time integrator.Extrapolation makes the convection problem linear hyperbolic and supports the stated stability property.
  • Modified fractional-step-θ scheme: A modified fractional-step-θ scheme combines implicit Stokes treatment and explicit or implicit convection treatment to avoid severe CFL restrictions and splitting errors simultaneously.The three-step scheme treats Stokes implicitly in the first and last stages, while the middle stage treats convection implicitly and viscosity explicitly.
  • Modified fractional-step-θ scheme: The modified scheme is second-order accurate, with stability restrictions reported empirically as much less restrictive than comparable IMEX schemes.Its formal stability analysis remains open.

4. Numerical examples

The numerical examples validate convergence, examine linear-system sparsity, and evaluate the discretization on benchmark problems in two and three dimensions. Projected jumps preserve optimal convergence while reducing facet unknowns and can improve linear-system sparsity.

  • The experiments target convergence, dependence of linear-system sparsity on velocity spaces, and performance on benchmark problems.The convergence tests use two-dimensional Stokes and Navier–Stokes problems with known exact solutions; later tests use three-dimensional settings.
  • The computations use the ngsflow add-on for NGSolve, direct solvers, and the reduced H(div)-conforming space with projected jumps for benchmark problems.
  • 4.1. Stokes flow around obstacle: The obstacle-flow test obtains optimal orders O(h^k+1) for velocity L2 error and O(h^k) for velocity-gradient and pressure errors.Here k is the polynomial order of the element-interior velocity field.
  • 4.1. Stokes flow around obstacle: In the obstacle-flow test, projected jumps reduce tangential facet unknowns by one order while changing velocity errors only marginally.Pressure errors differ between formulations, but both converge optimally and the difference decreases.
  • 4.2. Kovasznay flow: The Kovasznay-flow results with and without projected jumps are almost identical across the reported convergence comparisons.The reported difference is very small and appears only for k = 1.

4.3. Linear systems - A comparison between DG and HDG methods

The comparison studies four DG or HDG discretizations for a three-dimensional vector-valued reaction Poisson problem, focusing on sparsity before and after static condensation. Projected-jump HDG has the most favorable reported sparsity pattern across polynomial degrees.

  • The comparison uses symmetric positive definite systems from a three-dimensional reaction Poisson problem and evaluates sparsity with sparse Cholesky factorization.Accuracy is not compared.
  • The reported quantities are total unknowns, condensed unknowns, nonzero entries in the system matrix, and nonzero entries in the Cholesky factor.
  • Four methods—HDG, projected-jump HDG, standard DG, and nodal DG—are tested on polynomial degrees k = 1 through k = 6.
  • For k ≥4, HDG has fewer nonzero entries in both the system matrix and Cholesky factor than the compared alternatives.At low degree, HDG requires many additional unknowns and system-matrix nonzeros.
  • Projected-jump HDG substantially reduces the HDG overhead at small k and outperforms all alternatives for every tested polynomial degree.

4.4. A two-dimensional benchmark problem

The two-dimensional benchmark evaluates spatial and temporal discretization for flow around a circular obstacle at Re = 100. Increasing polynomial degree yields rapid spatial convergence, while alternative operator-splitting methods permit larger time steps with second-order convergence.

  • Problem setup: The benchmark models laminar flow around a circular obstacle at Re = 100, producing a periodic vortex street.The quantities of interest are maximal and minimal drag and lift forces.
  • Numerical setup: The mesh uses unstructured triangles with one anisotropically refined quadrilateral layer around the disk.Figure 4 depicts the geometry, mesh, and a typical solution.
  • Spatial discretization: Spatial results show rapid convergence as the polynomial degree k increases.Table 3 reports accuracy for varying polynomial degrees.
  • Spatial discretization: The method achieves the same order of accuracy as literature results with substantially fewer degrees of freedom.The comparison uses FEATFLOW with continuous quadratic velocity and discontinuous piecewise-linear pressure.
  • Temporal discretization: Alternative operator-splitting methods allow much larger time steps and exhibit second-order convergence than the second-order IMEX scheme permits.The product decomposition method is more accurate than the modified fractional step method by one time level.

4.5. A three-dimensional benchmark problem

The three-dimensional benchmark studies computational effort and accuracy for unsteady flow around a cylindrical obstacle. High-order discretizations achieve literature-level accuracy with dramatically smaller computing times, although direct solvers limit scalability as problem complexity grows.

  • Numerical setup: The geometry generalizes the two-dimensional case to a cuboid channel containing a cylindrical obstacle.Figure 5 shows the used mesh and solution at t = 0.4.
  • Problem setup: The three-dimensional benchmark uses time-dependent inflow over [0, 8s], with Reynolds number varying from Re = 0 to Re = 100.The mesh contains 5922 tetrahedral elements and the study focuses on computational effort over accuracy.
  • Numerical results: The method achieves the same accuracy as literature results with dramatically smaller computing time.The comparison includes results from and.
  • Numerical results: Using polynomial degrees k > 2 increases efficiency further while retaining high-order implementation efficiency.The authors extend the conclusion that third-order methods are more efficient than lower-order methods.
  • Implementation and limitations: For increasing problem complexity, efficient linear solvers are mandatory because the reported tests used direct solvers only.The paper plans suitable Stokes preconditioners based on scalar p-version domain-decomposition results.

5. Conclusion

The paper concludes that a combined DG/HDG discretization can efficiently treat unsteady incompressible Navier–Stokes problems by matching spatial and temporal treatments to Stokes and transport operators.

  • Conclusion: Operator splitting separates the Navier–Stokes equations into linear Stokes-type and hyperbolic transport problems.The Stokes part uses projected-jumps H(div)-conforming HDG, while transport uses standard DG.
  • Conclusion: The projected-jumps formulation further improves the efficiency of the Hybrid DG treatment for linear systems.Numerical test cases demonstrate the method’s performance.
Loading 1508.04245v2…