Discontinuous Galerkin

From CFD-Wiki

Revision as of 09:15, 3 January 2012 by Peter (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

The discontinuous Galerkin (DG) method is often referred to as a hybrid, or mixed, method since it combines features of both finite element and finite volume methods. The solution is represented within each element as a polynomial approximation (as in FEM), while the interelement convection terms are resolved with upwinded numerical flux formulas (as in FVM). Theoretically, solutions may be obtained to an arbitrarily high order-of-accuracy.

While the discontinuous Galerkin method was developed in the early 1970's, it was not used for CFD simulations until the early 1990's when it was first used to solve the Euler equations by Cockburn and Shu. The solution of the Navier-Stokes equations with the DG method was first accomplished by Bassi and Rebay in 1997. As the method gained more attention in the CFD research community, further advances have come fairly rapidly. Researchers are now using the DG method to perform simulations of a wide variety of flow regimes. The method has been adapted for use with compressible and incompressible, steady and unsteady, as well as laminar and turbulent conditions.

In addition to being arbitrarily higher-order accurate, the discontinuous Galerkin also permits the formulation of very compact numerical schemes. The is due to the fact that the solution representations in each element are purposefully kept independent of the solutions in other cells, with inter-element communication occurring only with adjacent cells (elements sharing a common face). This characteristic, along with other favorable numerical properties, makes this method extremely flexible (easily handling a wide variety of element types and mesh topologies) and also allows a number of adaptive techniques (both h- and p- refinement) and solver acceleration strategies to be implemented in a rather straightforward manner.

Derivation for Hyperbolic Conservation Laws

The discontinuous Galerkin method is derived from the finite element method, which is itself a variational method. To obtain the governing equations for the DG method, we begin with the strong form of the hyperbolic conservation laws:

$\frac{\partial \mathbf{Q}}{\partial t} + \vec{\nabla} \cdot \vec{\mathbf{F}}(\mathbf{Q}) = \mathbf{S}$

Here, $\mathbf{Q}$ is an array of conserved quantities, and $\vec{\mathbf{F}}(\mathbf{Q})$ is an array of flux vectors describing the local transport of $\mathbf{Q}$ . $\mathbf{S}$ is an array of source terms which represents the (non-conservative) production or destruction of the conserved quantities. Following the variational techniques employed in the finite element method, this expression then is multiplied by a test function $\psi$ and integrated over the domain of interest. This results in the so-called weak form.

$\int_{\Omega} \psi \frac{\partial \mathbf{Q}}{\partial t} \, d\Omega + \int_{\Omega} \psi \left( \vec{\nabla} \cdot \vec{\mathbf{F}}(\mathbf{Q}) \right) \, d\Omega = \int_{\Omega} \psi \mathbf{S} \, d\Omega$

Communication with the boundary conditions are established by performing an integration by parts on the second term.

$\int_{\Omega} \psi \frac{\partial \mathbf{Q}}{\partial t} \, d\Omega + \int_{\partial \Omega} \psi \left( \vec{n} \cdot \vec{\mathbf{F}}(\mathbf{Q}) \right) \, d\sigma - \int_{\Omega} \vec{\nabla} \psi \cdot \vec{\mathbf{F}}(\mathbf{Q}) \, d\Omega = \int_{\Omega} \psi \mathbf{S} \, d\Omega$

This form is typically the starting point for most integral based methods, including the finite volume and finite element methods. (To obtain the finite volume method, simply set $\psi=1$ .) It is the decisions made at this point that ultimately determine the type of numerical scheme used to solve this system of equations.

Nearly all numerical techniques will choose to discretize the domain of interest $\Omega$ into a set of non-overlapping elements $\Omega_e$ such that $\cup \Omega_e \approx \Omega$ . Thus the above integrals can be distributed over the elements.

$\sum_e \left( \int_{\Omega_e} \psi \frac{\partial \mathbf{Q}}{\partial t} \, d\Omega + \int_{\partial \Omega_e} \psi \left( \vec{n} \cdot \vec{\mathbf{F}}(\mathbf{Q}) \right) \, d\sigma - \int_{\Omega_e} \vec{\nabla} \psi \cdot \vec{\mathbf{F}}(\mathbf{Q}) \, d\Omega \right) = \sum_e \int_{\Omega_e} \psi \mathbf{S} \, d\Omega$

The next choice concerns the representation of the solution state within each element. The finite element method is able to achieve higher-orders of accuracy by representing the solution state as a collection of polynomial functions. Depending on how these polynomials are specified, the local description of the solution may possess certain desirable properties such as $C^0$ or $C^1$ continuity. These continuity constraints typically require the local solution state in each element to be dependent on a number of degrees of freedom which are often distributed over some neighborhood of the element.

In contrast to most other finite element techniques, the discontinuous Galerkin method specifies that there be no explicit continuity requirements for the the solution representation at the element boundaries. The solutions state is represented by a collection of piecewise discontinuous functions. Formally, this is accomplished by setting the basis (or shape) functions within each element to zero everywhere outside of their associated element. Thus, the representation of the current solution state within each element $\mathbf{Q}_e$ is completely determined by a set of coefficients $\mathbf{q}_i$ which are unique to that element and is independent of the solution state in all other elements.

$\mathbf{Q}_e = \sum_{i \in I_e} \mathbf{q}_i \phi_i$

For Galerkin methods, the test functions in the governing equations are taken from the same set of functions used to describe the solution. Since these functions are nonzero in at most one element, then all of the volume integrals in the resulting system of equations need only be evaluated in one and only one element. Due to the discontinuous nature of the solution, however, the surface integrals in the above equations will be double valued.

$\sum_e \left( \int_{\Omega_e} \phi \frac{\partial \mathbf{Q}}{\partial t} \, d\Omega - \int_{\Omega_e} \vec{\nabla} \phi \cdot \vec{\mathbf{F}}(\mathbf{Q}) \, d\Omega - \int_{\Omega_e} \phi \mathbf{S} \, d\Omega \right)_e + \sum_{f \in \Sigma} \left( \int_{\sigma_f} \phi^- \left( \vec{n}^- \cdot \vec{\mathbf{F}}(\mathbf{Q}^-) \right) \, d\sigma + \int_{\sigma_f} \phi^+ \left( \vec{n}^+ \cdot \vec{\mathbf{F}}(\mathbf{Q}^+) \right) \, d\sigma \right) + \sum_{f \in \Gamma} \left( \int_{\sigma_f} \phi^- \left( \vec{n}^- \cdot \vec{\mathbf{F}}(\mathbf{Q}^-) \right) \, d\sigma + \int_{\sigma_f} \vec{n} \cdot \vec{\mathbf{F}}^b \, d\sigma \right) = 0$

where $\Sigma$ is the set of interior faces, and $\Gamma$ is the set of boundary faces. This situation can be neatly resolved through the use of upwinded numerical flux functions, similar to those employed by the finite volume method.

$\sum_e \left( \int_{\Omega_e} \phi \frac{\partial \mathbf{Q}}{\partial t} \, d\Omega - \int_{\Omega_e} \vec{\nabla} \phi \cdot \vec{\mathbf{F}}(\mathbf{Q}) \, d\Omega - \int_{\Omega_e} \phi \mathbf{S} \, d\Omega \right)_e + \sum_{f \in \Sigma} \int_{\sigma_f} \left( \phi^- - \phi^+ \right) \mathbf{F}^*(\mathbf{Q}^-,\mathbf{Q}^+;\vec{n}) \, d\sigma + \sum_{f \in \Gamma} \int_{\sigma_f} \phi^- \mathbf{F}^*(\mathbf{Q}^-,\mathbf{Q}^b;\vec{n}) \, d\sigma = 0$

where $\vec{n} = \vec{n}^- = -\vec{n}^+$ . We take note that

$\mathbf{F}^*(\mathbf{Q}^-,\mathbf{Q}^+;\vec{n}) = \mathbf{F}^*(\mathbf{Q}^-,\mathbf{Q}^+;\vec{n}^-) = -\mathbf{F}^*(\mathbf{Q}^-,\mathbf{Q}^+;\vec{n}^+)$

and rewrite the governing equations as

$\sum_e \left( \int_{\Omega_e} \phi \frac{\partial \mathbf{Q}}{\partial t} \, d\Omega - \int_{\Omega_e} \vec{\nabla} \phi \cdot \vec{\mathbf{F}}(\mathbf{Q}) \, d\Omega - \int_{\Omega_e} \phi \mathbf{S} \, d\Omega + \sum_{f \in \Sigma_e} \int_{\sigma_f} \phi \mathbf{F}^*(\mathbf{Q},\mathbf{Q}^+;\vec{n}) \, d\sigma + \sum_{f \in \Gamma_e} \int_{\sigma_f} \phi \mathbf{F}^*(\mathbf{Q},\mathbf{Q}^b;\vec{n}) \, d\sigma \right) = 0$

where $\vec{n}$ is taken as the outward facing normal for $\Omega_e$ and $\Sigma_e = \Omega_e \cap \Sigma$ and $\Gamma_e = \Omega_e \cap \Gamma$ are, respectively, the interior and boundary faces of $\Omega_e$ . $\mathbf{Q}^+$ is the solution state in the element which shares face $\sigma_f$ and $\mathbf{Q}^b$ is the solution state imposed by a given boundary condition.

Diffusive Operators

The following are known methods for dealing with diffusive flux terms.

Interior Penalty (IP) -- Arnold
Lifting Operators (BR2) -- Bassi and Rebay
Baumann-Oden (BO) -- Baumann and Oden
Local Discontinuous Galerkin (LDG) -- Cockburn and Shu
Recovery -- van Leer and Nomura
Compact Discontinuous Galerkin (CDG) -- Peraire and Persson

Discontinuous Galerkin

From CFD-Wiki

Derivation for Hyperbolic Conservation Laws

Diffusive Operators

References

Views

My wiki

wiki navigation

wiki search

wiki toolbox