slightly more complicated. One way is to define cells based on the vertices on the boundary, that then have volume on both sides. In this way, no fluxes along the periodic boundary need to be computed.
3.7. FINITE VOLUME METHODS OF HIGHER ORDER 53 the time discretization has to be considered as well. For the sake of simplicity, we will only consider explicit Euler time integration for the remainder of this chapter and examine different time integration methods in the next chapter.
Unfortunately, it turns out that even with this relaxed property, TVD schemes are not the final answer. First, Godunov proved that a linear TVD scheme can be at most of first order [70]. For this reason, higher order schemes have to be nonlinear. Furthermore, they have to reduce to first order at spatial local maxima, as was proven by Osher and Chakravarty [152]. In multiple dimensions, as was proven by Goodman and LeVeque [71], TVD schemes are at most first order. Finally, on unstructured grids, not even a plane wave can be transported without increasing the total variation. Nevertheless, nonlinear TVD finite volume schemes using higher order ideas have become the workhorse in academia and industry in the form of the MUSCL schemes using reconstruction and limiters.
An alternative development are positive schemes, as made popular by Spekreijse [185], respectively local extrema diminishing schemes (LED) as suggested by Jameson [99]. This were developments in parallel to the development of TVD schemes that looked particularly interesting in the context of unsteady flows. However, it turns out that these schemes have the same restrictions as TVD schemes. Nevertheless, this led to the development of interesting flux functions like CUSP and SLIP, which are widely used in aerodynamic codes.
3.7.2 Reconstruction
The reconstruction procedure is based on the primitive variablesq= (ρ, v1, v2, p)T, as this is numerically more stable than using the conservative variables [140]. At a given time t, the linear representation of a primitive variableq∈ {ρ, v1, ..., vd, p}in celliwith barycenter xi is given by
q(x) = qi+∇q·(x−xi), (3.18)
where qi is the primitive value corresponding to the conservative mean values in Ωi. The unknown components of∇q represent the slopes and are obtained by solving a least square problem []. In the case of the Navier-Stokes equations, these slopes can sometimes be reused for the computation of the viscous fluxes.
For a cell-centered method on cartesian grids, the computation of the gradients is straightforward. If unstructured grids are employed, the following procedure is suitable.
Let C be the closed polygonal curve that connects the barycenters of the neighbouring cells, see figure 3.6a. We then define the piecewise linear function qc onC by setting
qc(xj) =qj
for all barycenters xj of neighbouring cells. The least squares problem, which has to be solved for all primitive variables, is then:
Figure 3.6: Left: The curve C used in the construction of the least squares problem for an unstructured cell. Right: The intersection of the primary triangle with a dual cell.
min
∇q∈Rd
L(∇q) :=
Z
C
(q(x)−qc(x))2ds. (3.19) For a cell-vertex scheme, first the unique linear interpolants of the cell averages located at the nodes are computed in each primary cell Tj [142], as described in the section on viscous flows. Let these have the gradient ∇qj. The gradient on the dual box Ωi is then defined as
∇q= 1
|Ωi| X
j
∇qj|Tj∩Ωi|. (3.20)
Again, this procedure has to be done for each primitive variable. Note that|Tj∩Ωi|is zero except for the few primary cells that have a part in the creation of the box Ωi, see figure 3.6b. Thus, the reconstruction procedure is significantly easier to implement and apply for cell-vertex methods than for cell-centered methods in the case of unstructured grids.
3.7.3 Modification at the Boundaries
At a fixed wall, the boundary condition requires the flux to satisfyv·n= 0 in the evaluation point, thus the slopes have to satisfy the same condition. For the cell-vertex method, nothing needs to be done, since there the values at the boundary will be zero in the first place and thus, the linear interpolants computed on the primary cells have the required property.
However, for the cell-centered method, we have to modify the original velocity slopes (v1)x1,(v1)x2,(v2)x1 and (v2)x2. We first define the unit vector η = (η1, η2)T to be the one pointing from the cell barycenter to the evaluation point on the wall and the unit vector ϑ = (ϑ1, ϑ2)T = (−η2, η1)T to be perpendicular to η. To ease the following construction,
3.7. FINITE VOLUME METHODS OF HIGHER ORDER 55 we now restrict the set of possible slope modifications ( ¯v1)x1,( ¯v1)x2,( ¯v2)x1,( ¯v2)x2 satisfying
¯
v·n= 0 to those with a constant ϑ derivative. Thus we only modify the derivatives (v1)η and (v2)η. This does not define the new slopes uniquely, therefore we require the new slopes to deviate from the old slopes minimally in the euclidian norm. We obtain the least squares problem:
min(((¯v1)η −(v1)η)2 + ((¯v2)η −(v2)η)2|( ¯v1)x1,( ¯v1)x2,( ¯v2)x1,( ¯v2)x2 :v¯·n= 0).
With ∆ being the euclidian distance from the cell barycenter to the evaluation point on the edge, we can write the modified velocity on the evaluation point as v¯ = vi +
∆((v1)η,(v2)η)T, where vi is the mean velocity vector of the cell and we obtain:
¯
v·n = (vi+ ∆((v1)η,(v2)η)T)·n= 0
⇔(¯v1)ηn1+ (¯v2)ηn2 =−vi·n/∆.
Inserting this condition directly into the least squares functional allows us to solve the minimization problem to obtain the solution (¯v1)?η,(¯v2)?η. Then we can compute the modified slopes via
(¯v1)x1 =η1(¯v1)?η −η2(v1)ϑ, (¯v2)x1 =η1(¯v2)?η−η2(v2)ϑ, (¯v1)x2 =η2(¯v1)?η +η1(v1)ϑand (¯v2)x2 =η2(¯v2)?η +η1(v2)ϑ.
3.7.4 Limiters
As mentioned before, to obtain a TVD scheme, the scheme can be at most of first order at shocks and local extrema. Thus, on top of the reconstruction scheme, a so called limiter function is needed. Typically a slope limiter φ is employed for the switching between first and higher order spatial discretization:
q(x) = qi+φ∇q·(x−xi). (3.21)
If the limiter function is zero, the discretization is reduced to first order, while it can be of second order for other values. In one dimension, conditions can be found for a limiter to result in a TVD scheme. For these, the limiter is taken as a function of the left and right slope and should be between zero and two, with the more precise shape given in figure 3.7a).
Furthermore, the conditions for the resulting scheme to be second order are demonstrated in figure 3.7b). A large number of limiters has been suggested that fulfil these conditions, for example the superbee limiter or the van Albada limiter. On structured grids, these can be applied directionwise and it is possible to prove the TVD property.
Figure 3.7: Left: Admissible region for a TVD slope limiter. Right: Admissible region for a 2nd order TVD slope limiter.
If the grid is unstructured, heuristic approaches must be used. We suggest either the Barth-Jesperson limiter [8] or the limiter proposed by Venkatakrishnan [212]. The latter is defined in cell i as follows:
φ = min
j∈N(i)
(∆2ij +) + 2∆ij∆ij
∆2ij+ 2∆2ij+ ∆ij∆ij +
. (3.22)
Here, = 10−6 and
∆ij =
0, if ∆ij∆ij <0, qj−qi, else,
where the indicator ∆ij is given by
∆ij :=qx1(x1ij −x1i)−qx2(x2ij −x2i) with the edge center (x1ij, x2ij).
Regarding boundaries, in the ghost cells at the inlet boundaries, the slopes and limiters are computed in the usual way. As the ghost cells are the definite border of the computa-tional domain, no values beyond the ghost edges are interpolated or incorporated in any way, as is done for the fixed wall.