HIGH-ORDER SPECTRAL METHODS FOR COMPRESSIBLE FLOWS ON UNSTRUCTURED MESHES

The present work considers the use of high-order numerical methods for compressible aerodynamic solutions on unstructured meshes. The primary interest of the paper is to investigate the proper solution accuracy and efficiency of the Spectral Finite Volume (SFV) and Spectral Difference (SD) methods. The SFV method can easily be adapted to a finite volume solver while the SD method is designed for a finite difference framework. The mesh element support is also different for these schemes and will most likely affect the solution of a particular problem. The present work compares the order of accuracy and resolution capabilities of the two schemes for literature test cases.


INTRODUCTION
High order methods are necessary on the analysis of complex flows to reduce the number of mesh elements one would otherwise need if using traditional second-order schemes.Assessment of solution reliability in such cases may also require that a specified solution accuracy be met by the spatial discretization scheme.The spectral methods here considered were developed as an alternative to the discontinuous Galerkin, k-exact, ENO/WENO highorder schemes.Their common objective is to allow the implementation of a simpler and more efficient scheme, while still achieving high-order spatial accuracy.
The spectral finite volume (SFV) and spectral difference (SD) schemes are implemented in an unstructured, two dimensional Euler framework.The schemes allow for 2nd-, 3rd-and 4th-order accurate solutions, which use the Roe Riemann solver as the basic numerical flux scheme.The SFV method, as developed by Wang [1] circa 2002, was further extended to handle hyperbolic systems, two and three dimensions, smooth and discontinuous flow solutions.This method was originally incorporated by the authors into a 2-D finite volume solver due to its similar underlying structure.However, the SFV method, at its current state, supports only triangular mesh elements.This proves to be a drawback if one is concerned with viscous flows, mainly due to the boundary layer discretization, where quadrilateral elements are better suited for the gradients discretization.The SD method, on the other hand, shares many of the SFV scheme attributes, as previously described, only that its formulation is suited for quadrilateral unstructured meshes.The present work represents an effort to evaluate the SD method as an alternate high-order method for the applications of interest to the authors, namely external compressible aerodynamic flows.The CFD tools developed by the authors and used in the present work have been verified and validated previously [2,3,4].
The paper is organized as follows.First, the theoretical formulation of the SFV method for the Euler equations is presented.The Spectral Difference method formulation is shown next.The order of accuracy of both schemes are measured and compared and the final section presents numerical results for well-known literature test cases.

General Formulation
In the present work, the 2-D Euler equations are solved in their integral form as where P = Eî + F .The application of the divergence theorem to Eq. ( 1) yields The vector of conserved variables, Q, and the convective flux vectors, E and V , are given by The standard CFD nomenclature is being used here.Hence, ρ is the density, u and v are the Cartesian velocity components in the x and y directions, respectively, p is the pressure, and e t is the total energy per unit volume.The system is closed by the equation of state for a perfect gas where e i is the specific internal energy, and the ratio of specific heats, γ, is set as 1.4 for all computations in this work.In the finite volume context, for fixed meshes, Eq. ( 2) can be rewritten for the i-th mesh element as where Q i is the cell averaged value of Q at time t and V i is the volume, or area in 2-D, of the i-th mesh element.

Spatial Discretization
The spatial discretization process determines a k-th order discrete approximation to the integral in the right-hand side of Eq. (5).In order to solve it numerically, the computational domain, Ω, with proper initial and boundary conditions, is discretized into N non-overlapping triangles, the spectral volumes (SVs), such that One should observe that the spectral volumes could be composed by any type of polygon, given that it is possible to decompose its bounding edges into a finite number of line segments Γ K , such that In the present paper, however, the authors assume that the computational mesh is always composed of triangular elements.Hence, although the theoretical formulation is presented for the general case, the actual SV partition schemes are only implemented for triangular grids.The boundary integral in Eq. ( 5) can be further discretized into the convective operator form where nf is the number of faces, or edges in 2-D, of SV i , and A r represents the area, or the length in 2-D, of the r-th edge of the SV.Given the fact that n is constant for each line segment, the integration on the right side of Eq. ( 8) can be performed numerically with a k-th order accurate Gaussian quadrature formula where (x rq , y rq ) and w rq are, respectively, the Gaussian quadrature points and the weights on the r-th edge of SV i , nq = integer((k + 1)/2) is the number of quadrature points required on the r-th edge, and h will be defined in the forthcoming discussion.For the second-order schemes, one Gaussian quadrature point is used in the integration.Given the coordinates of the end points of the element edge, z 1 and z 2 , one can obtain the quadrature point as the middle point of the segment connecting the two end points, For the third and fourth order schemes, two Gaussian points are necessary along each line segment.Their values are given by Using the method described above, one can compute values of Q i for instant t for each SV.Due to the discontinuity of the reconstructed values of the conserved variables over SV boundaries, one must use a numerical flux function to approximate the flux values along the cell boundaries.
The previously described procedure follows exactly the standard finite volume method.For a given order of spatial accuracy, k, for Eq. ( 8), using the SFV method, each SV i element must have at least degrees of freedom (DOFs).This corresponds to the number of control volumes (CVs) that SV i shall be partitioned into.If one denotes by CV i,j the j-th control volume of SV i , the cell-averaged conservative variables, q, at time t, for CV i,j is computed as where V i,j is the volume of CV i,j .Once the cell-averaged conservative variables, or DOFs, are available for all CVs within SV i , a polynomial, p i (x, y) ∈ P k−1 , with degree k − 1, can be reconstructed to approximate the q(x, y) properties inside SV i , i.e., where h represents the maximum edge length of all CVs within SV i .The polynomial reconstruction process is discussed in details in the following section.For now, it is sufficient to say that this high-order reconstruction is used to update the cell-averaged state variables at the next time step for all the CVs within the computational domain.Note that this polynomial approximation is valid within SV i and the use of numerical fluxes are necessary across SV boundaries.Integrating Eq. ( 5) in CV i,j , one can obtain the integral form for the CV mean state variable where f represents the E and F fluxes, and nf is the number of edges of CV i,j .The numerical integration can be performed with a k-th order accurate Gaussian quadrature formulation, similarly to the one used for the SV elements in Eq. ( 9).As previously stated, the flux integration across SV boundaries involves two discontinuous states, to the left and to the right of the edge.This flux computation can be carried out using an exact or approximate Riemann solver, or even a flux splitting procedure, which can be written in the form where q L is the conservative variable vector obtained by the p i polynomial applied at the (x rq , y rq ) coordinates and q R is the same vector obtained with the p nb polynomial in the same coordinates of the edge.Note that the nb subscript represents the element to the right of the edge, whereas the i subscript denotes the CV to its left.As the numerical flux integration in the present paper is based on one of the forms of a Riemann solver, this is the mechanism which introduces the upwind and artificial dissipation effects into the method, making it stable and accurate.In this work, the authors have used the Roe flux difference splitting method [5] to compute the numerical flux, i.e., Here, B is the Roe dissipation matrix in the direction normal to the edge.The B matrix, in the edge-normal direction, is defined as The B matrix has four real eigenvalues, namely, where v n is the velocity component normal to the edge and a is the speed of sound.Let T be the matrix composed of the right eigenvectors of B. Then, this matrix can be diagonalized as where Λ is the diagonal matrix composed of the eigenvalues of B, which can be written as where Λ and T are calculated as a function of the Roe averaged properties [5].Furthermore, Λ uses the magnitude of the eigenvalues.Finally, one ends up with the semi-discrete SFV scheme for updating the CVs, which can be written as where the right hand side of Eq. ( 21) is the equivalent convective operator, C(q i,j ), for the j-th control volume of SV i .It is important to emphasize that some edges of the CVs, resulting from the partition of the SVs, lie inside the SV element in the region where the polynomial is continuous.For such edges, there is no need to compute the numerical flux, as described above.Instead, one uses analytical formulas for the flux computation, i.e., without numerical dissipation.

Reconstruction Procedure
The evaluation of the conserved variables at the quadrature points is necessary in order to perform the flux integration over the mesh element edges.These evaluations can be achieved by reconstructing conserved variables in terms of some base functions using the DOFs within a mesh element, denoted here as a Spectral Volume, or SV for short [1,6].The present work has carried out such reconstructions using polynomial functions.Let P m denote the space of m-th degree polynomials in two dimensions.Then, the minimum dimension of the approximation space that allows P m to be complete is In order to reconstruct q in P m , it is necessary to partition the SV into N m non-overlapping control volumes, or CVs, such that The reconstruction problem, for a given continuous function in S i and a suitable partition, can be stated as finding p m ∈ P m such that With a complete polynomial basis, e l (x, y) ∈ P m , it is possible to satisfy Eq. ( 24).Hence, p m can be expressed as where e is the base function vector, and b is the reconstruction coefficient vector, The substitution of Eq. ( 25) into Eq.( 24) yields e l (x, y)dS = q i,j .
If q denotes the [q i,1 , • • • , q i,N m ] T column vector, Eq. ( 26) can be rewritten in matrix form as where the S reconstruction matrix is given by and, then, the reconstruction coefficients b can be obtained as provided that S is non-singular.With the substitution of Eq. (29) into Eq.( 24), p m is, then, expressed in terms of shape functions Equation ( 30) gives the value of the conserved state variable, q, at any point within the SV and its boundaries, including the quadrature points, (x rq , y rq ).The major advantage of the SFV method is that the reconstruction process does not need to be carried out for every mesh element S i .One can compute these coefficients as a pre-processing step and they do not change along the simulation.This is a major difference when compared to methods such as ENO and WENO, for which every mesh element has a different reconstruction process at each time step.The polynomial base functions for the linear, quadratic and cubic reconstructions are listed in Table 1.Clearly, the linear, quadratic and cubic polynomial reconstructions will yield, respectively, 2nd-, 3rd-and 4th-order spatial discretization numerical schemes.
For most of the proposed applications in compressible aerodynamics, the use of a limiter technique is necessary to render the scheme stable.The authors conducted several studies regarding different limiter formulations and their performance.The interested reader is referred to Ref. [3] for a throughout discussion of this subject.The cited reference also presents the curved boundary formulation.This technique is responsible to maintain the methods design order of accuracy near the domain boundaries and also helps to reduce the overall number of mesh elements.The extension of such techniques for the Spectral Difference method have already been carried out elsewhere in the literature and is currently under development by the authors on their implementation.

Linear Reconstruction
For the linear SFV method reconstruction, m = 1, one needs to partition a SV in three sub-elements as in Eq. ( 22) and use the base vector as defined in Table 1.The partition scheme is given for a standard element.The partition for this case is uniquely defined.The structured aspect of the CVs within the SVs is used to determine neighborhood information and generate the global connectivity data considering a hash table search algorithm [7].
The linear partition is presented in Fig. 1(a).It yields a total of 7 points, 9 edges (6 are external edges and 3 are internal ones) and 9 quadrature points.The linear polynomial for the SFV method depends only on the base functions and on the partition shape.The integrals of the reconstruction matrix in Eq. ( 28) are obtained analytically [8] for fundamental shapes.The shape functions, in the sense of Eq. (30), are calculated and stored in memory for the quadrature points, (x rq , y rq ), of the standard element.Such shape functions have the exact same value for the quadratures points of any other SV of the mesh, provided they all have the same partition.There is one quadrature point located at the middle of the every CV edge.

Quadratic Reconstruction
For the quadratic reconstruction, m = 2, one needs to partition a SV in six subelements and use the base vector as defined in Table 1.The partition scheme is also given in this work for a right triangle.The nodes of the partition are obtained in terms of barycentric coordinates of the SV element nodes in the same manner as the linear partition.The structured aspect of the CVs within the SVs is used to determine neighborhood information and generate the connectivity table.The ghost creation process and edge-based data structure is the same as for the linear reconstruction case.The partition used in this work [9] is shown in Fig. 1(b).It has a total of 13 points, 18 edges (9 external edges and 9 internal ones) and 36 quadrature points.The shape functions, in the sense of Eq. (30), are obtained as in the linear partition.The reader should note that, in this case, the base functions have six terms that shall be integrated.Again, these terms are obtained exactly and kept in memory [8].In this case, two quadrature points are required per CV edge.

Cubic Reconstruction
For the cubic reconstruction, m = 3, one needs to partition the SV in ten sub-elements and to use the base vector as defined in Table 1.The ghost creation process and edge-based data structure is the same as for the linear and quadratic reconstruction cases.As a matter of fact, the same algorithm utilized to perform these tasks can be applied to higher order reconstructions.The partition used in this work is the improved cubic partition [9], presented in Fig. 1(c) and it has a total of 21 points, 30 edges (12 external edges and 18 internal ones) and 60 quadrature points.The shape functions, in the sense of Eq. (30), are obtained as in the linear partition in a pre-processing step.As with the quadratic reconstruction, each CV edge has two quadrature points [2].

SPECTRAL DIFFERENCE METHOD
The Spectral Difference method (SD) formulation in the present work considers the unsteady compressible 2-D Navier Stoker equations in conservative form where Q is the vector of conserved variables, E and F are the fluxes as described in the SFV method section.The implementation follows the formulation presented in Refs.[10,11].The SD method employes a Finite Difference like scheme rather than a Finite Volume formulation.In this sense, it represents a different approach to the solver programming and algorithm considerations.For instance, to achieve an efficient implementation, all elements in the physical domain (x, y) are transformed into a unit square element in the computational domain.Such transformation can be written as where K is the number of points used to define the physical element, (x i , y i ) are the Cartesian coordinates of those points, and M i (ξ, η) are the shape functions of the geometric transformation.The metrics and the Jacobian of the transformation can be computed in a pre-processing step and kept in memory given the stationary aspect of the problems here considered.The governing equations in the physical domain are then transferred into the computation domain, and can be rewritten as where In the standard element, two sets of points are defined, namely the solution points and the flux points, illustrated in Fig. 2. In order to construct a degree (N − 1) polynomial in each coordinate direction, solution at N points are required.The solution points in 1-D are chosen to be the Gauss points defined by The flux points are selected to be the Gauss-Lobatto points given by Using the solution at N solution points, a degree (N − 1) polynomial can be built using the following Lagrange basis defined as similarly, using the fluxes at (N + 1) flux points, a degree N polynomial can be built for the flux using a similar Lagrange basis defined as The reconstructed solution for the conserved variables in the standard element is just the tensor products of the two one-dimensional polynomial, Qi,j Similarly, the reconstructed flux polynomials take the following form The reconstructed fluxes are only element-wise continuous, but discontinuous across cell interfaces.For the inviscid flux, a Riemann solver is employed to compute a common flux at interfaces to ensure conservation and stability.In the present work, the Roe approximate Riemann solver is considered, as described in the previous sections.The numerical flux includes the artificial dissipation and upwind characteristic of the method.In order to compute the inviscid flux derivatives the method computes the conservative variables at the flux points.The inviscid fluxes are then updated followed by the numerical flux calculations.
The derivatives of the fluxes are computed at the solution points using the derivatives of Lagrange operators l as At the moment of writing, the SD method implementation is currently under way on the available numerical framework.A second, third and fourth-order reconstruction procedure will be considered for benchmark simulations against the SFV method.

Temporal Discretization
The temporal discretization is concerned with solving a system of ordinary differential equations.In the present work, the authors use a third-order, TVD Runge-Kutta scheme [12].Rewriting Eq. ( 21) in a concise ODE form, one obtains Figure 2. Possible flux points (blue squares) and solution points (red circles) distribution for the SD method. where and w rq f Roe (q L (x rq , y rq ), q R (x rq , y rq ), n r )A r . (45) Hence, the time marching scheme can be written as where the n and n + 1 superscripts denote, respectively, the values of the properties at the beginning and at the end of the n-th time step.The α coefficients are α 1 = 3/4, α 2 = 1/4, α 3 = 1/3 and α 4 = 2/3.

Accuracy Assessment
The accuracy of the SFV and SD methods is tested for the linear wave equation, for −1 ≤ x ≤ 1 and u(x, 0) = u 0 (x) with periodic boundary condition at the domain extremes.For this setup, the initial condition is u 0 (x) = sin(πx).No limiters were used in this study.A third-order TVD Runge-Kutta method is employed for time integration with a ∆t value of 10 −4 , in order to make the discretization error time-step independent.shows the L 1 and L ∞ error norms produced using the SFV method with the Gauss-Lobatto point distribution.One can note that all schemes are able to achieve their expected order of accuracy.The same study is carried out for the SD method, presented in Table 3 and Fig. 3.In these tables, N DOF represents the number of degrees of freedom for the problem.One can observe that the expected order of accuracy, even for coarse grids, is achieved for both methods.However, the performance of the fifth and sixth-order SD scheme is questionable, specially for the finer meshes.This is mostly likely related to the oscillatory behavior of the polynomial interpolation, due to the internal points distribution chosen.As the grid is refined, the errors actually increase in both norms, which give a negative order of accuracy.Other test cases are considered for the 1-D SFV and SD methods.A Gaussian pulse with a half width equal to 0.05 is used as an initial condition, again with periodic boundary conditions.This simulation was carried out for various orders of accuracy, k = 1, 2, 3, 4 and 6 for t = 2 with 100 SVs.The Gauss-Lobatto points distribution was used for CV partitioning.The results are presented in Fig. 4 and 5.The first-order  simulation smeared the pulse to the point that it is hardly recognizable.This is due to the extra amount of dissipation associated with such low-order schemes.The second-order simulation retained the shape of the initial condition but also smeared the pulse and produced an oscillatory behavior, which, in turn, is associated with the dispersion properties of the method.The third, fourth and sixth-order simulation yielded nice results just like the analytical solution, also shown in the figure.

Ringleb Flow
The Ringleb flow simulation consists of an internal flow, which has an analytical solution for the Euler equations derived with the hodograph transformation [13].The analytical solution is used as initial condition for the SFV simulations here discussed.The flow depends on the inverse of the stream function, k, and the velocity magnitude, v t .In the present simulations, these parameters are chosen as k = 0.4 and k = 0.6, in order to define the bounding walls, and v t = 0.35 to define the inlet and outlet boundaries.For such configuration, the test case represents an irrotational and isentropic flow around a symmetric blunt obstacle.An interesting property of the Ringleb test case is that transition of flow regime, from subsonic to supersonic, for example, is shockless [14].In order to measure the order of the implemented SFV method, four meshes are considered for the mesh refinement study, corresponding to 128, 512, 2048 and 8192 spectral volume elements.The analytical solution is computed for all meshes in order to measure how close the numerical results are to the exact solution.The error with respect to the analytical solution is computed using the L 1 and L ∞ norms of the density.Figure 6 shows the 2048-element grid and the Mach number contours computed in this grid with the fourth order SFV method, using the corresponding high-order boundary representation.Table 4 presents the L 1 and L ∞ error norms of the density for the present calculations with the high-order boundary representation for the SFV method.

Airfoil Simulation
A simulation with the NACA 0012 airfoil is considered for the SFV method.The Euler equations are solved on a coarse mesh with 716 cells and 358 nodes, from which 40 define the airfoil wall.This O-grid mesh is presented in Fig. 7.The airfoil profile itself is collapsed on the trailing edge.The far field boundary radius is 50 chord units.
The free stream flow replicate the conditions of the experimental data [15], that is, free stream Mach number value of M ∞ = 0.8 and 0 deg angle-of-attack.Simulations with the second, third and fourth order SFV schemes are performed, along with the first-order Roe scheme. Figure 8 shows the computed Cp values obtained with with the first-order Roe scheme and third-order SFV method, considering quadratic curved boundaries.It is clear from the Cp distribution that the first-order scheme introduces too much dissipation and essentially smears out the shock wave.The high-order Cp distribution, on the other hand, is remarkably close to the experimental results in Fig. 8, particularly for the shock position, considering the crude mesh discretization.These show the potential for high-order methods, on such applications, to ease the mesh generation process.One should observe, however, that the experimental results consider the presence of the boundary layer and the consequent shock-boundary layer interaction that necessarily occurs in the experiment.For the numerical solution, the shock presents a sharper resolution, as one can expect for an Euler simulation.
Another relevant simulation is performed to assess the benefits of the curved boundary implementation, namely, the measure of entropy error s levels at the airfoil boundary.Because the diffusive flux vectors are zero, there is no physical dissipation mechanism that produces heat in regions of smooth flow, away from shocks.If no external heat is added into the flow, then it is called adiabatic and, from the first law of thermodynamics, it follows that entropy, give by is constant throughout the field if no shocks are present.Therefore, the entropy error s , defined as is a good measure of the accuracy of a numerical solution obtained with a method to approximately solve the Euler equations.Figure 9 presents the entropy error generate by the third-order SFV method with linear and curved boundary cells for the coarse mesh, which has only 40 cells to represent the whole airfoil geometry.The curved boundary approach is able to produce the smaller error levels than the linear approach.One can even note, form the figure, that at position x = 0 there is an increase of entropy error, due to the presence of the shock wave in this region.This, again, demonstrates that such extension indeed improves the overall accuracy of the SFV method.

CONCLUSIONS
The present work compares two high-order methods for several test cases.The formal order of accuracy of the SFV and SD methods is achieved for the model convection equation.Moreover, solutions of the Euler equations for the SFV method are in good agreement with the available analytical and experimental data.Even though additional numerical techniques are required to handle the desired flow regime of compressible aerodynamics, such as the limiter and high-order boundary representation, the SFV method is only capable to work with triangular mesh elements.The ability to handle quadrilateral meshes is a desirable feature of the SD method, specially for RANS simulations.Current work is in place to extend the SD method for the 2-D formulation and proper results and benchmark simulations against the SFV method will be presented at the conference.

Figure 3 .
Figure 3. Spectral Difference solution error versus mesh spacing plot order of accuracy measurement of a sine wave convection.

Figure 4 .
Figure 4. Simulation of a traveling Gaussian pulse with SFV schemes of various orders at t = 2.

Figure 5 .
Figure 5. Spectral Difference error versus mesh spacing plot order of accuracy measurement of a Gaussian pulse convection.

Table 2 Table 2 .
Accuracy assessment of the 1-D SFV method for the wave equation.Gauss-Lobatto CV partition.

Table 3 .
Accuracy assessment of the 1-D SD method for the wave equation.

Table 4 .
Figure 6.Ringleb flow mesh and Mach number contour results for fourth-order SFV method.Accuracy assessment of SFV method for the Ringleb flow test case.