KozCG

The equations of compressible flow

The governing equations solved are the 3D unsteady Euler equations, governing inviscid compressible flow, written here in the flux-conservative form

where the summation convention on repeated indices has been applied, and is the density, the velocity vector, is the specific total energy, is the specific internal energy, and is the pressure. , , and are source terms that arise from the application of the method of manufactured solutions, used for verification; these source terms are zero when computing practical engineering problems. The system is closed with the ideal gas law equation of state

where is the ratio of specific heats. The exact form of the equation of state is of secondary importance, any arbitrary (i.e. analytic or tabulated) equation of state can be used to relate density and internal energy to pressure.

Discretization in space

To arrive at a continuous Galerkin finite element method we start from the weak form of the governing equations above,

which requires that the error in the numerical solution, sampled at discrete points, , using the weighting functions , vanish on the whole domain, , in an integral sense. Numerically approximating the solution as , where denotes the unknowns at the discrete node , leads to the Galerkin weak statement

Integrating the flux term by parts, neglecting the resulting boundary integral (assuming zero flux on the problem boundary), and applying and , yields the final weak form for the whole domain, ,

All integrals above are evaluated by breaking up the domain, , into sub-domains as a sum over integrals over discrete elements, ,

where the inner summation is over points forming (gather) and the outer summation is over tetrahedra connected to point (scatter). The above equations result in the following semi-discrete system of equations

where the comma denotes a derivative. The size of the consistent mass matrix is , where is the number of nodes of the computational mesh. According to the sum above, only those elements contribute to a given row of which contain node (scatter), thus is sparse.

Discretization in time

The equation above is discretized in time using a Lax-Wendroff (Taylor-Galerkin) scheme, implemented using two stages:

The above scheme combined with linear shape functions is identical to a two-stage Runge-Kutta Galerkin finite element method, and thus central differencing without damping. Stabilization is obtained by using constant shape functions for the half step solution where the gather results in element quantities, followed by a scatter step using linear shape functions resulting in nodal quantities. Assuming linear tetrahedra, the combined spatial and temporal discretization that achieves this yields the following staggered scheme:

where is the vector of element-centered solutions at the half step. Note that the first step discretizes the flux integral before integration by parts, while the second one after integration by parts, hence the difference in sign.

Flux-corrected transport

Flux-corrected transport (FCT) is a solution to circumvent the consequence of Godunov's theorem, which states that no linear scheme of order greater than 1 will yield monotonic solutions. Accordingly, FCT is a nonlinear scheme that combines a high-, and a low-order scheme using limiting.

In the FCT scheme used, the high-order solution at the new time step can be written as

where and denote the solution increments of the high-, and low-order schemes, respectively. In the equations above it is the last term that is limited in a way to avoid spurious oscillations as

The high-order scheme scheme, given above is symbolically written as

From the equation above we construct a low order scheme by lumping the mass matrix and adding mass diffusion

where is the lumped mass matrix and is a diffusion coefficient. Using guarantees a monotone low order solution. If we rewrite the above equation as

the difference between the right hand sides of the high and low order schemes can be recognized as

also called as the anti-diffusive element contributions (AEC). The AEC is then limited, , and applied to advance the solution to the next time step using the equation above, where is the limit coefficient for a given element .

The limiting procedure

The weak sum above shows that the left hand side is the consistent mass matrix. is lumped by summing the rows to the diagonals. Inverting distributed across computers would be costly. Using a lumped (diagonal) matrix instead reduces computational cost and software complexity at the cost of some additional numerical error but does not reduce the order of accuracy of the method. If the mesh does not move and its topology does not change, the left hand side needs no update between time steps.

The high-order right hand side is computed by the two-step procedure given above. Since the two steps are staggered, the gather takes information from nodes to cell centers, followed by a scatter, moving information from cells to nodes, both steps are contained within a single right hand side calculation within a time step. Within the two steps there is no need for parallel communication as an element always resides on a given mesh partition and only mesh nodes are shared between processing elements (PEs).

A step of the limiting procedure is to compute the mass diffusion term. Using linear tetrahedra, this is given for each element by

where is the element Jacobian, computed from the triple product of the edge vectors of the tetrahedron given by vertices A,B,C, and D. The above equation is also used to compute the anti-diffusive element contributions. Once the AECs are computed for each element, the next step is to sum all positive (negative) anti-diffusive element contributions to node i

The limiting procedure requires the maximum and minimum nodal values of the low-order solution and the previous solution,

Another alternative is to only consider the low order solution, when computing the allowed solution bounds, which leads to the so-called 'clipping limiter' in place of the equation above. This is followed by computing the maximum and minimum nodal values of all elements,

then computing the maximum and minimum unknowns of the elements surrounding each node,

The limit coefficients will be computed (see below) based on and the maximum and minimum increments and decrements the nodal solution values are allowed to achieve,

Defining the ratios of positive and negative element contributions for each node that ensure monotonicity as

the limit coefficient for each element is taken as the most conservative ratio

The limited AEC is then scatter-added to nodes of the low-order solution:

The above procedure is general, works on the numerical solution (instead on fluxes or slopes), and does not require a Riemann solver. The same procedure can be used for each scalar in a system of equations. This works well for independent scalars, but for coupled system of equations additional techniques have been developed to reflect the coupled nature of the equations in the limiting procedure.