Introduction

  In spectral methods the question as to what basis functions $\phi_n(x)$ to use to approximate a smooth function $f(x) \approx \sum_{n=0}^N \hat{f}_n \phi_n(x)$ in a separable domain is relatively easily answered. Typically, boundary conditions, fast transforms, and numerical quadrature dictate this choice depending on the particular partial differential equation considered. In multi-domain spectral methods which are used in complex-geometry computational domains the choice of the subdomain shape is an important factor as well. One can choose arbitrary polygons, for example, to accommodate geometric complexity but such a choice may prevent the use of polynomial basis functions and lead to prohibitive computational difficulty in maintaining interfacial continuity among subdomains.

For the version of the finite element method the choice of subdomain shape and associated functions bases dictates the efficiency of the method. Starting with the pioneering work of B.A. Szabo in the mid-seventies several versions of this approach have been formulated for both solid mechanics as well as fluid dynamics in the founding papers of Babuska et al [5,6] and by Patera et al Patera84, KBP85. Similar formulations were developed to allow greater flexibility in handling geometric complexities and local refinement requirements in the works of [9], [10], [11] and [12]. In figure 1.1 we demonstrate how a complex three-dimensional volume can be formed using prisms. We solved the Helmholtz problem within this volume, and show that as we increased the polynomial order of our representation we obtained an exponential decrease in the approximation error.

  

Figure 1.1: Top: Unstructured prism mesh, Bottom: p-type convergence for the Helmholtz problem with $\lambda=1$ and solution sin(y)sin(z). The face geometry was downloaded from the web-site: http://www.3dcafe.com/

The basis or shape functions employed in the aforementioned formulations typically involve either Legendre or Chebyshev one-dimensional polynomials, and multi-dimensional expansions are constructed using tensor products for quadrilateral or hexahedral elements [13], [14]. A basis directly associated with boundary and interior elemental nodes usually lacks hierarchy and in addition a variable p-order per element is not readily achievable [7]. However, a modal basis can be constructed that consists of vertex modes, edge modes, and interior modes in two-dimensions that are hierarchical and allow for a variable order within as well as along the boundary of the element [13]. In addition to the flexibility in handling more easily nonuniform resolution requirements, hierarchical bases can lead to better conditioning of mass and stiffness matrices [14], and thus fast iterative algorithms can be effectively employed in the solution algorithms.

The need for constructing p-type bases in triangular domains became evident from the beginning. A. Peano [15] constructed a hierarchical triangular basis based on area (barycentric) coordinates by selecting as nodal variables high-order tangential derivatives along each side evaluated at the midside and imposing appropriate continuity constraints for C⁰ and C¹ elements. A variation of this construction was later developed [13] that introduces Legendre polynomials to avoid round-off error for high-order p-expansions. Both approaches however require special integration rules which are quite complicated at high polynomial order. Dubiner [16] first developed an alternative hierarchical basis for triangular domains, which unlike Peano's basis is based on cartesian coordinates and preserves the tensor product property that leads to sum factorization and low operation count, a crucial property for p-type finite elements and spectral methods.

Dubiner's basis was implemented in [20] using a Galerkin formulation of the Navier-Stokes equations, and it was found to be competitive in cost with the nodal basis on quadrilaterals employed in the spectral element method [7]. This new modal basis employs Jacobi polynomials of mixed weight to automatically accommodate exact numerical integration using standard Gauss-Jacobi one-dimensional quadrature rules. In particular, exploiting the tensor product property of the basis, multi-dimensional integrals can be evaluated efficiently as a series of one-dimensional integrals; similarly the cost of evaluating derivatives or squares of a function is maintained at operation count ${\cal O}(N^{d+1})$ where N is the number of modes per direction as in the quadrilateral or hexahedral elements and d is the dimension of the element.

Recent trends in mesh generation have emphasized the use of hybrid discretizations consisting of triangular, quadrilateral, tetrahedral, prismatic, pyramidal and hexahedral subdomains in order to provide greater flexibility in refinement and derefinement procedures. For high-order discretizations this implies that suitable basis functions have to be developed to conform to the subdomain shapes while at the same time preserve the good approximation properties and efficiency of the aforementioned bases. In this thesis, we develop a unified description of such hybrid basis functions following earlier developments in [16] and [21]. We also develop five types of basis functions which are either modal or nodal or mixed and which may or may not be hierarchical.

There has been recently an interesting debate as to what are the relative advantages of structured and unstructured grids, with a renewed interest in cartesian grids for complex-geometry aerodynamic flows [22]. The current confusion stems from the fact that most finite element and finite volume formulations in use today produce solutions, which depend strongly on the quality of the mesh. Specifically, for highly distorted grids convergence is questionable, and in most cases convergence rates are typically less than second-order. To this end, some - relatively few - efforts have addressed the development of hybrid grids, i.e. a mixture of structured and unstructured grids in order to combine the merits of both discretizations in the context of complex-geometry aerodynamic flows [24], [25], [26]. Such methods lead to more flexible geometric discretization and resolution placement, but they are still of low-order accuracy, i.e. at most second-order.

We also develop a new formulation for compressible Navier-Stokes and compressible magnetohydrodynamics solutions employing high-order spectral/hp element discretization on mixed element grids. A discontinuous Galerkin formulation is developed both for the convective as well as the diffusion contributions that allows multi-domain representation with a discontinuous (i.e. globally L²) trial basis. This discontinuous basis is orthogonal, hierarchical, and maintains a tensor-product property (even for non-orthogonal elements), a key property for the efficient implementation of high-order methods. Due to this basis orthogonality the resulting mass matrix is diagonal, and thus the proposed method is matrix-free given that an explicit time-stepping is involved. The conservativity property is maintained automatically by the discontinuous Galerkin formulation, and monotonicity is controlled by varying the order of the spectral expansion around discontinuities.

  

Figure 1.2: Hybrid discretization around a multi-element airfoil. Only part of the domain is shown.

An example of a hybrid discretization is shown in figure 1.2 for simulation of flow past a multi-element airfoil. This problem was studied by Barth [4] who used unstructured disretizations, i.e. only triangles. There are two regions of structured discretizations as shown in these two plots: First, around the four airfoils where boundary layers are present and where all the vorticity is generated. Second, in the wake region where one single family of square elements is embedded within the unstructured discretization. This uniform meshing minimizes dispersion errors of the traveling vortex street and facilitates maximum storage efficiency. We will revisit this example in section 7.2.1 in order to present simulation results at subsonic speeds.