3.1 Spatial Discretisation 空间离散[cfd-3-1]

The finite difference method was among the first approaches applied to the numerical solution of differential equations. It was first utilised by Euler, probably in 1768. The finite difference method is directly applied to the differential form of the governing equations. The principle is to employ a Taylor series expansion for the discretisation of the derivatives of the flow variables. Let us for illustration consider the following example.

Suppose we would like to compute the first derivative of a scalar function \(U(x)\) at some point \(x_0\). If we develop now \(U(x_0+\Delta x)\) as a Taylor series in \(x\), we obtain

With this, the first derivative of \(U\) can be approximated as

The above approximation is of first order, since the truncation error (abbreviated as \(\mathcal{O}(\Delta x)\)), which is proportional to the largest term of the remainder, goes to zero with the first power of \(\Delta x\) (for a discussion on the order of accuracy see Chapter 10). The same procedure can be applied to derive more accurate finite difference formulae and to obtain approximations to higher-order derivatives.

An important advantage of the finite difference methodology is its simplicity. Another advantage is the possibility to easily obtain high-order approximations, and hence to achieve high-order accuracy of the spatial discretisation. On the other hand, because the method requires a structured grid, the range of application is clearly restricted. Furthermore, the finite difference method cannot be directly applied in body-fitted (curvilinear) coordinates, but the governing equations have to be first transformed into a Cartesian coordinate system - or in other words - from the physical to the computational space (Fig. 3.2). The problem herewith is that the Jacobian of coordinate transformation appears in the flow equations (see, e.g., Appendix A.1). This Jacobian has to be consistently discretised in order to avoid the introduction of additional numerical errors. Thus, the finite difference method can be applied only to rather simple geometries. Nowadays, it is sometimes utilised for the direct numerical simulation of turbulence (DNS), but it is only very rarely used for industrial applications. More details to the finite difference method can be found for example in [43], or in textbooks on the solution of partial differential equations.

The finite volume method directly utilises the conservation laws - the integral formulation of the Navier-Stokes/Euler equations. It was first employed by McDonald [44] for the simulation of 2-D inviscid flows. The finite volume method discretises the governing equations by first dividing the physical space into a number of arbitrary polyhedral control volumes. The surface integral on the right-hand side of Equation (2.19) is then approximated by the sum of the fluxes crossing the individual faces of the control volume. The accuracy of the spatial discretisation depends on the particular scheme with which the fluxes are evaluated.

There are several possibilities of defining the shape and position of the control volume with respect to the grid. Two basic approaches can be distinguished:

• Cell-centred scheme (Fig. 3.6a) - here the flow quantities are stored at the centroids of the grid cells. Thus, the control volumes are identical to the grid cells.
• Cell-vertex scheme (Fig. 3.6b) - here the flow variables are stored at the grid points. The control volume can then either be the union of all cells sharing the grid point, or some volume centred around the grid point. In the former case we speak of overlapping control volumes, in the second case of dual control volumes.

We shall discuss the advantages and disadvantages of cell-centred and cell-vertex formulations in both chapters on spatial discretisation.

The main advantage of the finite volume method is that the spatial discretisation is carried out directly in the physical space. Thus, there are no problems with any transformation between the physical and the computational coordinate system, like in the case of the finite difference method. Compared to the finite differences, one further advantage of the finite volume method is that it is very flexible - it can be rather easily implemented on structured as well as on unstructured grids. This renders the finite volume method particularly suitable for the treatment of flows in complex geometries.

Since the finite volume method is based on the direct discretisation of the conservation laws, mass, momentum and energy are also conserved by the numerical scheme. This leads to another important feature of the method, namely the ability to compute weak solutions of the governing equations correctly. However, one additional condition has to be fulfilled in the case of the Euler equations. This is known as the entropy condition. It is necessary because of the non-uniqueness of the weak solutions. The entropy condition prevents the occurrence of unphysical features like expansion shocks, which violate the second law of thermodynamics (decrease of the entropy). As a further consequence of the conservative discretisation, the Rankine-Hugoniot relations, which must hold across a solution discontinuity (such as a shockwave or a contact discontinuity), are satisfied directly.

It is interesting to note that under certain conditions, the finite volume method can be shown to be equivalent to the finite difference method, or to a low-order finite element method. Because of its attractive properties, the finite volume method is nowadays very popular and in wide use. It will be presented in the following chapters.

The finite element method was originally employed for structural analysis only. It was first introduced by Turner et al. [45] in 1956. About ten years later, researchers started to use the finite element method also for the numerical solution of field equations in continuous media. However, only with the beginning of the 90's, did the finite element method gain popularity in the solution of the Euler and the Navier-Stokes equations. A good introduction into the classical finite element methodology can be found in [46]. Applications to flow problems are described in [47], [48], and more recently in [49].

The finite element method, as it is in general applied to the solution of the Euler/Navier-Stokes equations, starts with a subdivision of the physical space into triangular (in 2-D) or into tetrahedral (in 3-D) elements. Thus, an unstructured grid has to be generated. Depending on the element type and the required accuracy, a certain number of points at the boundaries and/or inside an element is specified, where the solution of the flow problem has to be found. The total number of points multiplied with the number of unknowns determines the number of degrees of freedom. Furthermore, the so-called shape functions have to be defined, which represent the variation of the solution inside an element. In practical implementations, linear elements are usually employed, which use the grid nodes exclusively. The shape functions are then linear distributions, whose value is zero outside the corresponding element. This results in a second-order accurate representation of the solution on smooth grids.

Within the finite element method, it is necessary to transform the governing equations from the differential into an equivalent integral form. This can be accomplished in two different ways. The first one is based on the variational principle, i.e., a physical solution is sought, for which a certain functional possesses an extremum. The second possibility is known as the method of weighted residuals or the weak formulation. Here, it is required that the weighted average of the residuals is identically zero over the physical domain. The residuals can be viewed as the errors of the approximation of the solution. The weak formulation has the same advantage as the finite volume discretisation of the conservation laws - it allows the treatment of discontinuous solutions such as shocks. Therefore, the weak formulation is preferred over the variational methodology.

The finite element method is attractive because of its integral formulation and the use of unstructured grids, which are both preferable for flows in or around complex geometries. The method is also particularly suitable for the treatment of non-Newtonian fluids. The finite element method has a very rigorous mathematical foundation, particularly for elliptic and parabolic problems. Although it can be shown in certain cases that the method is mathematically equivalent to the finite volume discretisation, the numerical effort is noticeably higher. This may explain why the finite volume method became more popular. However, both methods are sometimes combined - particularly on unstructured grids. So for example, the treatment of the boundaries and the discretisation of the viscous fluxes is usually "borrowed" from the finite element method.

There are few other numerical schemes which are only seldom used in practice, but which are, in certain situations, superior to the methods discussed above. Two particular approaches should be mentioned here briefly.

The first such example is the Spectral Element Method [50]-[53]. The spectral element method combines the geometrical flexibility of the finite element technique with the high-order spatial accuracy (e.g., 10th-order) and the rapid convergence rate of the spectral schemes [54]. The method is based on a high-order polynomial representation of the solution (usually Lagrangian interpolants), combined with a standard Galerkin finite element method, or the method of weighted residuals. The spectral element method is appropriate either for a particular problem in which high-order regularity is guaranteed, or for which high-order regularity is not the exception, like in incompressible fluid mechanics. It is especially suitable for vortical flows. The advantage of the spectral element method is primarily its non-diffusive, non-dispersive approximation of the convection operator, and its good approximation of convection-diffusion boundary layers. The method can treat geometrically and physically complex problems, supposed the condition of high-order regularity is fulfilled. Apart from the rather narrow range of applications, the principal disadvantage of the spectral element method is its very high numerical effort as compared for example to the finite volume method.

Another discretisation scheme, which gained recently some interest, is the so-called Gridless Method [55]-[57]. This method employs only clouds of points for the spatial discretisation. It does not require that the points are connected to form a grid as in conventional structured or unstructured grid schemes. The gridless method is based on the differential form of the governing equations, written in the Cartesian coordinate system. Gradients of the flow variables are determined by a least-squares reconstruction, using a specified number of neighbours surrounding the particular point. The gridless method is neither a finite difference nor a finite volume or a finite element approach since coordinate transformations, face areas or volumes do not have to be computed. It can be viewed as a mix between the finite difference and the finite element method. The principal advantages of the gridless method are its flexibility in solving flows about complex configurations (similar to unstructured methods), and the possibility to locate or cluster the points (or the clouds of points) where it is appropriate. For example, it would be easily possible to select only the neighbours in the characteristic directions when computing gradients. However, there is one unresolved problem. Although the gridless method solves the conservation law form of the Euler or the Navier-Stokes equations, it is not clear whether conservation of mass, momentum and energy is really ensured.

Whichever spatial discretisation scheme we might select, it is important to ensure that the scheme is consistent, i.e., that it converges to the solution of the discretised equations, when the grid is sufficiently refined. It is therefore very important to check how much the solution changes, if the grid is refined (e.g., if we would double the number of grid points). If the solution improves only marginally, we speak of grid converged solution. Another rather self-evident requirement is that the discretisation scheme possesses the order of accuracy, which is appropriate for the flow problem being solved. This rule is sometimes given up in favour of faster convergence, particularly in industrial environment (bad solution is better than no solution). This is of course a very dangerous practice. We shall return to the question of accuracy, stability and consistency later in Chapter 10.

So far, we discussed only the basic choices which exist for the spatial discretisation. But within each of the above three main methods - finite difference, finite volume and finite element - various numerical schemes exist to perform the spatial discretisation. In this context, it is convenient to differentiate between the discretisation of the convective and the viscous fluxes (\(\vec{F}_c\) and \(\vec{F}_v\) in Eq. (2.19), respectively). Because of the physical nature of the viscous fluxes, the only reasonable way is to employ central differences (central averaging) for their discretisation. Thus, their discretisation on structured grids is straightforward. On unstructured triangular or tetrahedral grids, the viscous fluxes are best approximated using the Galerkin finite element methodology, even in the case of a finite volume scheme [58]. The situation becomes more complicated for unstructured mixed grids, where a modified averaging of gradients is more appropriate [59]-[63].

However, the real variety is found in the discretisation of the convective fluxes. In order to classify the individual methodologies, we will restrict our attention to schemes developed for the finite volume method, although most of the concepts are also directly applicable to the finite difference or the finite element method.

To the first category we may count schemes, which are based solely on central difference formulae or on central averaging, respectively. These are denoted as central schemes. The principle is to average the conservative variables to the left and to the right in order to evaluate the flux at a side of the control volume. Since the central schemes cannot recognise and suppress an odd-even decoupling of the solution (i.e., the generation of two independent solutions of the discretised equations), the so-called artificial dissipation (because of its similarity to the viscous terms) has to be added for stabilisation. The most widely known implementation is due to Jameson et al. [64]. On structured grids, it is based on a blend of 2nd- and 4th-differences scaled by the maximum eigenvalue of the convective flux Jacobian. A combination of an undivided Laplacian and biharmonic operator is employed on unstructured grids [65]. The scheme can be improved remarkably using different scaling factors for each equation. This approach is known as the matrix dissipation scheme [66]. It should be mentioned that on unstructured, mixed element grids the explicit Runge-Kutta time-stepping scheme can become unstable, when combined with the conventional central scheme [67].

On the other hand, there are more advanced spatial discretisation schemes, which are constructed by considering the physical properties of the Euler equations. Because they distinguish between upstream and downstream influences (wave propagation directions), they are termed upwind schemes. They can be roughly divided into four main groups:

• flux-vector splitting,
• flux-difference splitting,
• total variation diminishing (TVD), and
• fluctuation-splitting schemes.

Each of these is described briefly in the following subsections.

One class of the flux-vector splitting schemes decomposes the vector of the convective fluxes into two parts according to the sign of certain characteristic variables, which are in general similar to but not identical to the eigenvalues of the convective flux Jacobian. The two parts of the flux vector are then discretised by upwind biased differences. The very first flux-vector splitting schemes of this type were developed in the beginning of the 1980's by Steger and Warming [68] and by Van Leer [69], respectively. A second class of flux-vector splitting schemes decompose the flux vector into a convective and a pressure (an acoustic) part. This idea is utilised by schemes like AUSM (Advection Upstream Splitting Method) of Liou et al. [70], [71], or the CUSP scheme (Convective Upwind Split Pressure) of Jameson [72], [73], respectively. Further similar approaches are the Low-Diffusion Flux-Splitting Scheme (LDFSS) introduced by Edwards [74], or the Mach number-based Advection Pressure Splitting (MAPS) scheme of Rossow [75], [76]. The second group of flux-vector splitting schemes gained recently larger popularity particularly because of their improved resolution of shear layers, but only a moderate computational effort. An advantage of the flux-vector splitting schemes is also that they can be quite easily extended to real gas flows, as opposed to flux-difference splitting or TVD schemes. We shall return to real gas simulations further below.

The second group - flux-difference splitting schemes - is based on the solution of the locally one-dimensional Euler equations for discontinuous states at an interface. This corresponds to the Riemann (shock tube) problem. The values on either side of the interface are generally termed as the left and right state. The idea to solve the Riemann problem at the interface between two control volumes was first introduced by Godunov [77] back in 1959. In order to reduce the numerical effort required for an exact solution of the Riemann problem, approximate Riemann solvers were developed, e.g., by Osher et al. [78] and Roe [79]. Roe's solver is often used today because of its excellent resolution of boundary layers and a crisp representation of shocks. It can be easily implemented on structured as well as on unstructured grids [80].

The idea of TVD schemes was first introduced by Harten [81] in 1983. The TVD schemes are based on a concept aimed at preventing the generation of new extrema in the flow solution. The principal conditions for a TVD scheme are that maxima must be non-increasing, minima non-decreasing, and no new local extrema may be created. Such a scheme is called monotonicity preserving. Thus, a discretisation methodology with TVD properties allows it to resolve a shock wave without any spurious oscillations of the solution. The TVD schemes are in general implemented as an average of the convective fluxes combined with an additional dissipation term. The dissipation term can either depend on the sign of the characteristic speeds or not. In the first case, we speak of an upwind TVD scheme [82], in the second case of a symmetric TVD scheme [83]. The experience shows that the upwind TVD scheme should be preferred since it offers a better shock and boundary layer resolution than the symmetric TVD scheme. The disadvantage of the TVD schemes is that they cannot be easily extended to higher than second-order spatial accuracy. This limitation can be overcome using the ENO (Essentially Non-Oscillatory) discretisation schemes [84]-[89].

The last group - the fluctuation-splitting schemes - provides for true multidimensional upwinding. The aim is to resolve accurately also those flow features which are not aligned with the grid. This is a significant advantage over all above upwind schemes, which split the equations according only to the orientation of the grid cells. Within the fluctuation-splitting methodology, the flow variables are associated with the grid nodes. Intermediate residuals are computed as flux balances over the grid cells, which consists of triangles in 2-D and of tetrahedra in 3-D. The cell-based residuals are then distributed in an upwind-biased manner to the nodes. After that, the solution is updated using the nodal values. In the case of systems of equations (Euler or Navier-Stokes), the cell-based residuals have to be decomposed into scalar waves. Since the decomposition is not unique in 2-D and in 3-D, several approaches were developed in the past. The variety reaches from the wave model of Roe [90], [91] over the algebraic scheme of Sidilkover [92] to the most advanced characteristic decomposition method [93]-[96]. Despite the above mentioned advantage over the dimensionally split Riemann, TVD, etc. solvers, the fluctuation-splitting approaches are so far used only in research codes. This can be attributed to the complexity and the high numerical effort, as well as to convergence problems.

You may ask now, what are the benefits and the drawbacks of the individual spatial discretisation methods. Generally speaking, central schemes require lower numerical effort, and hence less CPU time per evaluation, as compared to upwind schemes. On the other hand, upwind schemes are able to capture discontinuities much more accurately than central schemes. Furthermore, because of their lower numerical diffusion, the upwind schemes can resolve boundary layers using less grid points. Particularly, Roe's flux-difference splitting scheme and the CUSP flux-vector splitting scheme allow a very accurate computation of boundary layers. The negative side of the upwind schemes emerges for second- or higher-order spatial accuracy. The problem is that the so-called limiter functions (or simply limiters) have to be employed in order to prevent the generation of spurious oscillations near strong discontinuities. Limiters are known to stall the convergence of an iteration scheme, because of their accidental switching in smooth flow regions. A remedy was suggested by Venkatakrishnan [97]-[99], which works satisfactorily for most practical cases. However, small wiggles in the solution must be taken into account. Another disadvantage of the limiter functions is that they require high computational effort, particularly on unstructured grids.

With respect to real gas simulations, and in particular to chemically reacting flows, several extensions of the upwind discretisation schemes were presented. For the case of fluids in thermodynamic and chemical equilibrium, modifications of the Van Leer flux-vector splitting [69] and of Roe's approximate Riemann solver [79] were described in Refs. [100]-[102]. Formulations of the both upwind methods for the more complex case of flows with non-equilibrium chemistry and thermodynamics were provided in [103]-[108] and in the references cited therein. In particular, the articles [102], [106] and [107], respectively, give a good overview of the methodologies employed for the upwind discretisation schemes. Recently, a summary of the governing equations together with Jacobian and transformation matrices, which may be required by an upwind scheme, was presented in Ref. [109] for the case of chemically reacting flows.