chapter. Principles of Solution of the Governing Equations 第3章 控制方程求解原理[cfd-0004]

Let us at the beginning turn our attention to the first step - the spatial discretisation of the Navier-Stokes equations, i.e., the numerical approximation of the convective and viscous fluxes, as well as of the source term. Many different methodologies were devised for this purpose in the past and the development still continues. In order to sort them, we can at first divide the spatial discretisation schemes into the following three main categories: finite difference, finite volume, and finite element. All these methods rely on some kind of grid in order to discretise the governing equations (2.19). Basically, there exist two different types of grids:

• Structured grids (Fig. 3.2) - each grid point (vertex, node) is uniquely identified by the indices \(i,j,k\) and the corresponding Cartesian coordinates \(x_{i,j,k}\), \(y_{i,j,k}\), and \(z_{i,j,k}\). The grid cells are quadrilaterals in 2-D and hexahedra in 3-D. If the grid is body-fitted like in Fig. 3.1a, we speak also of curvilinear grid.
• Unstructured grids (Fig. 3.3) - grid cells as well as grid points have no particular ordering, i.e., neighbouring cells or grid points cannot be directly identified by their indices (e.g., cell 6 adjacent to cell 119). In the past, the grid cells were triangles in 2-D and tetrahedra in 3-D. Today, unstructured grids usually consist of a mix of quadrilaterals and triangles in 2-D and of hexahedra, tetrahedra, prisms and pyramids in 3-D, in order to resolve the boundary layers properly. Therefore, we speak in this case of hybrid or of mixed grids.

The main advantage of structured grids follows from the property that the indices \(i,j,k\) represent a linear address space - also called the computational space, since it directly corresponds to how the flow variables are stored in the computer memory. This property allows it to access the neighbours of a grid point very quickly and easily, just by adding or subtracting an integer value to or from the corresponding index (e.g., like (i+1), (k-3), etc. - see Fig. 3.2). As one can imagine, the evaluation of gradients, fluxes, and also the treatment of boundary conditions is greatly simplified by this feature. The same holds for the implementation of an implicit scheme, because of the well-ordered, banded flux Jacobian matrix. But there is also a disadvantage. This is the generation of structured grids for complex geometries. As sketched in Fig. 3.4, one possibility is to divide the physical space into a number of topologically simpler parts - blocks - which can be more easily meshed. We therefore speak of multiblock approach [20]-[24]. Of course, the complexity of the flow solver is increased, since special logic is required to exchange physical quantities or fluxes between the blocks. Additional flexibility is added, if the grid points at both sides of an interface can be placed independently of each other, i.e., if the grid lines are allowed not to meet at a block boundary (like inside C or F in Fig. 3.4). Those grid points, which are located only on one side of a block interface are called hanging nodes. The advantage of this approach is quite obvious - the number of grid lines can be chosen separately for each block as required. The price paid for the enhanced flexibility is an increased overhead for the conservative treatment of the hanging nodes [25], [26]. The multiblock methodology also offers interesting possibilities with respect to the implementation of the flow solver on a parallel computer by means of domain decomposition. However, very long times (often weeks or months) are still required for the grid generation in the case of complex configurations.

Another methodology, related to block structured grids, represents the so-called Chimera technique [27]-[32]. The basic idea here is to generate first the grids separately around each geometrical entity in the domain. After that, the grids are combined together in such a way that they overlap each other where they meet. The situation is depicted in Fig. 3.5 for a simple configuration. The crucial operation is an accurate transfer of quantities between the different grids at the overlapping region. Therefore, the extension of the overlap is adjusted accordingly to the required interpolation order. The advantage of the Chimera technique over the multiblock approach is the possibility to generate the particular grids completely independent of each other, without having to take care of the interface between the grids. On the other hand, the problem of the Chimera technique is that the conservation properties of the governing equations are not satisfied through the overlapping region.

The second type of grids are the unstructured grids. They offer the largest flexibility in the treatment of complex geometries [33]. The main advantage of the unstructured grids is based on the fact that triangular (in 2-D) or tetrahedral grids (in 3-D) can be in principle generated automatically, independent of the complexity of the domain. In practice, it is of course still necessary to set some parameters appropriately, in order to obtain a good quality grid. Furthermore, in order to resolve the boundary layers accurately, it is advisable to employ in 2-D rectangular and in 3-D prismatic or hexahedral elements near solid walls [34]-[41]. Another benefit of such mixed grids is the reduction of the number of grid cells, edges, faces and possibly also of grid points. One should keep in mind though that the generation of mixed grids is not trivial for geometrically demanding cases. Nevertheless, the time required to built an unstructured, mixed grid for a complex configuration is still significantly lower than what is necessary for a multiblock structured grid. Since the geometrical fidelity of the flow simulations is nowadays rapidly increasing, the ability to generate grids fast and with minimum user interaction becomes more and more important. This is particularly true in industrial environment. Further advantage of unstructured grids is that solution dependent grid refinement and coarsening can be handled in a relatively native and seamless manner. To mention also the disadvantages of unstructured methods, one of them is the necessity to employ sophisticated data structures within the flow solver. Such data structures work with indirect addressing, which, depending on the computer hardware, leads to more or less reduced computational efficiency. The memory requirements are in general higher as compared to the structured schemes, too. But despite all difficulties, the capability to handle geometrically complex problems in short turn-around times still weights much more. From this point, it is not surprising that for example nearly all vendors of commercial CFD software switched over to unstructured flow solvers. A detailed review of various methodologies for spatial and temporal discretisation on unstructured grids appeared recently in [42].

Having generated the grid, the next question is how to actually discretise the governing equations. As we already said, we can basically choose between three methodologies: finite differences, finite volumes, and finite elements. We want to discuss all of them briefly in the next subsections.

As already mentioned at the beginning of this chapter, the prevailing number of numerical schemes for the solution of the Euler and the Navier-Stokes equations applies the method of lines, i.e., a separate discretisation in space and in time. This approach offers the largest flexibility, since different levels of approximation can be easily selected for the convective and the viscous fluxes, as well as for the time integration - just as required by the problem solved. Therefore, we shall follow this methodology here. For the discussion of other methods, where time and space discretisations are coupled, the reader is referred to Ref. [43].

When the method of lines is applied to the governing equations (2.19), it leads, written down for each control volume, to a system of coupled ordinary differential equations in time

For clarity, we omitted any cell indices. In Eq. (3.3), \(\Omega\) denotes volume of the control volume and \(\vec{R}\) stands for the complete spatial (finite volume) discretisation including the source term - the so-called residual. The residual is a non-linear function of the conservative variables \(\vec{W}\). Finally, \(\overline{M}\) represents what is termed the mass matrix. For a cell-vertex scheme, it relates the average value of \(\vec{W}\) in the control volume to the point values at the associated interior node and the neighbouring nodes [110], [111]. In the case of a cell-centred scheme, the mass matrix can be substituted by an identity matrix, without compromising the temporal accuracy of the scheme. The same holds for a cell-vertex scheme applied on a uniform grid, since then the nodes coincide with the centroids of the control volumes. The mass matrix is a function of the grid only and couples the system of differential equations (3.3). For steady-state cases, where time accuracy is not a concern, the mass matrix can be "lumped", i.e., replaced by the identity matrix. In this way, the expensive inversion of \(\overline{M}\) can be avoided and the system (3.3) is decoupled. In this respect, it is important to realise that at the steady-state, the solution accuracy is determined solely by the approximation order of the residual. Thus, the mass matrix becomes important only for cell-vertex schemes applied to unsteady flows.

If we assume a static grid, we may take the volume \(\Omega\) and the mass matrix outside the time derivative. Then, we can approximate the time derivative by the following non-linear scheme [43]

with

being the solution correction. The superscripts \(n\) and \((n+1)\) denote the time levels (\(n\) means the current one). Furthermore, \(\Delta t\) represents the time step.

The scheme in Eq. (3.4) is 2nd-order accurate in time if the condition

is fulfilled, otherwise the time accuracy is reduced to 1st-order. Depending on the settings of the parameters \(\beta\) and \(\omega\), we can obtain either explicit (\(\beta=0\)) or implicit time-stepping schemes. We shall discuss this two main classes briefly in the following paragraphs, and in more detail later in Chapter 6.

A family of implicit time integration schemes is obtained from Eq. (3.4) by setting \(\beta\neq 0\). Very popular for the simulation of unsteady flows is the 3-point implicit backward-difference scheme with \(\beta=1\) and \(\omega=1/2\), which is 2nd-order accurate in time. In this case, the scheme is mostly employed within the so-called dual time-stepping approach [148]-[150], [110], [111], where a steady-state problem is solved in pseudo-time at each physical time step.

For the solution of stationary flow problems, a scheme with \(\omega=0\) is more suitable, since it requires less computer memory. Herewith, if we linearise the residual \(\vec{R}^{n+1}\) in Eq. (3.4) about the current time level, we obtain the scheme

The term \(\partial\vec{R}/\partial\vec{W}\) is denoted as the flux Jacobian. It constitutes a large sparse matrix. The expression enclosed in parenthesis on the left-hand side of Eq. (3.8) is also referred to as the implicit operator. As already discussed above, the mass matrix \(\overline{M}\) can be replaced by the identity matrix, without influencing the steady state solution. The parameter \(\beta\) in Eq. (3.8) is generally set to 1, which results in a 1st-order accurate temporal discretisation. A 2nd-order time accurate scheme is obtained for \(\beta=1/2\). However, this is not advised since the scheme with \(\beta=1\) is much more robust, and the time accuracy plays no role for steady problems anyway.

The principal advantage of implicit schemes as compared to explicit ones is that significantly larger time steps can be used, without hampering the stability of the time integration process. In fact, for \(\Delta t\rightarrow\infty\) the scheme (3.8) transforms into standard Newton's method, which exhibits quadratic convergence. However, the condition for quadratic convergence is that the flux Jacobian contains the complete linearisation of the residual. Another important advantage of implicit schemes is their superior robustness and convergence speed in the case of stiff equation systems and/or source terms, which are often encountered in real gas simulations, turbulence modelling, or in the case of highly stretched grids (high Reynolds number flows). On the other hand, the faster (in terms of time steps or iterations) and the more robust an implicit scheme is, the higher is usually the computational effort per time step or iteration. Therefore, an explicit scheme accelerated by multigrid can be equally or even more efficient. Furthermore, implicit schemes are significantly more difficult to vectorise or to parallelise than their explicit counterparts.

Written down for each control volume, the implicit scheme in Eq. (3.8) represents a large system of linear equations, which has to be solved for the update \(\Delta\vec{W}^{n}\) at each time step \(\Delta t\). This task can be accomplished using either a direct or an iterative method.

The direct methods are based on the exact inversion of the left-hand side of Eq. (3.8) using either the Gaussian elimination or some direct sparse matrix method [151], [152]. Although quadratic convergence was demonstrated on structured [153]-[156] as well as on unstructured grids [157], direct methods are not an option for 3-D problems because they require an excessively high computational effort and a huge amount of computer memory.

Thus, the only practical method for larger grids or 3-D problems are iterative methods. Here, the linear system is solved for \(\Delta\vec{W}^{n}\) at each time step using some iterative matrix inversion methodology. In order to reduce the memory requirements and also to increase the diagonal dominance, the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\) is mostly based on linearisation of a 1st-order accurate spatial discretisation of the right-hand side. The two main consequences of this approximation are that the quadratic convergence of Newton's scheme cannot be achieved and that the maximum time step becomes limited. On the other hand, the numerical effort of an iteration step is significantly reduced, which leads to a numerically highly efficient scheme.

In the case of structured grids, iterative methods like the Alternating Direction Implicit (ADI) scheme [158]-[161], the (line) Jacobi or the Gauss-Seidel relaxation scheme [162]-[166], and particularly the Lower-Upper Symmetric Gauss-Seidel (LU-SGS; also referenced to as LU-SSOR - Lower-Upper Symmetric Successive Overrelaxation) scheme [167]-[171] are mainly employed. All these methods are based on splitting of the implicit operator into a sum or product of parts, which can be each inverted more easily. Because of the associated factorisation error (the difference with respect to the original matrix) and also the simplification of the flux Jacobian, it does not pay off to solve the linear system very accurately. In fact, only one iteration is carried out at each time step of the ADI and the LU-SGS method.

Implicit iterative methods for unstructured grids are in the most cases based on the Gauss-Seidel relaxation scheme [172]-[175]. In order to improve the convergence, it is possible to use the red-black Gauss-Seidel methodology. Its extension to unstructured grids was demonstrated in Refs. [176]-[178]. A particularly interesting possibility is also offered by an implementation of the LU-SGS scheme on unstructured grids [18], [179], [180], because of its very low memory requirements and numerical effort.

Because of the success of the line-implicit methods on structured grids, a few attempts were made to adopt this methodology on unstructured grids [181], [182]. The approach was to construct continuous lines such that each grid point or each grid cell (in the case of a cell-centred scheme) is visited only once - the so-called Hamiltonian tour [183]. The lines were oriented primarily in coordinate directions, but they were folded at the boundaries and where necessary (therefore they were nicknamed "snakes"). A tri-diagonal solver was then employed to invert the left-hand side of Eq. (3.8). Later on, it was recognised that folding the lines can slow down the convergence. To overcome this, each line was broken up into multiple linelets [184]. However, the performance on a vector computer was rather poor. The idea of linelets was also employed to improve the convergence of an explicit scheme on highly stretched viscous unstructured grids using an implicit solver in the direction across the boundary layer [145].

More sophisticated iterative techniques, which treat the linear equation system in a more global way, are the so-called Krylov subspace methods. Their development was triggered by the introduction of an efficient iterative scheme for solving large, sparse linear systems - namely the conjugate gradient method by Hestenes and Stiefel [185]. The original conjugate gradient method is restricted to Hermitian positive definite matrices only, but for an \(n\times n\) matrix it converges in at most \(n\) iterations. Since then, a variety of Krylov subspace methods was proposed for the solution of arbitrary non-singular matrices, as they occur in CFD applications. For example, there are methods like the Conjugate Gradient Squared (CGS) [186], the Bi-Conjugate Gradient Stabilised (Bi-CGSTAB) [187], or the Transpose-Free Quasi-Minimum Residual (TFQMR) [188] scheme.

However, the most widely employed method is the Generalised Minimal Residual (GMRES) scheme developed by Saad and Schultz [189]. If we rewrite the implicit scheme in Eq. (3.8) as

then \(\bar{J}\) represents a large, sparse, and non-symmetric matrix (the left-hand side). Starting from an initial guess \(\Delta\vec{W}_0^{n}\), the GMRES(\(m\)) method seeks a solution \(\Delta\vec{W}^{n}\) in the form \(\Delta\vec{W}^{n}=\Delta\vec{W}_0^{n}+\vec{y}_m\), where \(\vec{y}_m\) belongs to the Krylov subspace

such that the residual \(\Vert\bar{J}\,\Delta\vec{W}^{n}+\vec{R}^{n}\Vert\) becomes a minimum. The parameter \(m\) specifies the dimension of the Krylov subspace, or in other words the number of search directions (\(\bar{J}^{i}\vec{r}_0\)). Since all directions have to be stored, \(m\) is usually chosen between 10 and 40, the higher number being necessary for poorly conditioned matrices (which arise in the simulation of turbulent flows, real gas, etc.). GMRES has to be restarted, if no convergence is achieved within \(m\) sub-iterations. The GMRES method requires significantly more memory than, e.g., Bi-CGSTAB or TFQMR, but it is more robust, smoothly converging and usually also faster. A very detailed comparison of the various methodologies can be found in Ref. [190].

Nevertheless, as with other conjugate gradient methods, preconditioning is absolutely essential for CFD problems. Here, we solve

instead of the system in Eq. (3.9). The matrices \(\bar{P}_L\) and \(\bar{P}_R\) denote left and right preconditioners, respectively. The preconditioner should approximate \(\bar{J}^{-1}\) as close as possible, in order to cluster the eigenvalues near unity. On the other hand, it should be of course easy to invert. One particularly efficient preconditioner is the Incomplete Lower Upper factorisation method [191] with zero fill-in (ILU(0)). For the discussion of different preconditioning techniques in connection with GMRES the reader is referred to [97], [192]-[195].

Since the GMRES method requires a considerable amount of computer memory for storing the search directions and possibly also the preconditioning matrix, it is a good idea to circumvent an explicit formation and storage of the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\). This is offered by the so-called matrix-free approach. The idea is based on the observation that GMRES (and some other Krylov subspace methods) employs only matrix vector products of the form

which can be simply approximated by finite-differences as

thus requiring only residual evaluations. The parameter \(\epsilon\) has to be chosen with some care, in order to minimise the numerical error (see, e.g., [196] or [197]). Another, and even more important, advantage of the matrix-free approach is that (numerically) accurate linearisation of a high-order residual \(\vec{R}^{n}\) can be easily utilised in the implicit scheme. Hence, the quadratic convergence of Newton's scheme can be achieved at moderate costs. In this case we speak of Newton-Krylov approach [197]-[201], [110]. Practical experience indicates that from all Krylov subspace methods, GMRES is best suited for the matrix-free implementation [202]. An interesting possibility is to utilise the LU-SGS scheme as a preconditioner for the matrix-free GMRES method. Since the LU-SGS scheme also does not require an explicit storage of the flux Jacobian, the memory requirements can be even further reduced. The computational efficiency of this approach was recently demonstrated for 3-D inviscid and laminar flows on unstructured grids [203].

The convergence of an implicit scheme can also be enhanced by using multigrid. There are basically two possible ways. First, we can employ multigrid inside an implicit scheme - as a solver for the linear equation system (3.9) arising at each time step, or as a preconditioner for one of the conjugate gradient methods [204], [205]. Second, the implicit scheme itself can serve as a smoother within the FAS multigrid method, which is applied directly to the governing equations [206]-[209], [178]. Some investigations show that at least for purely aerodynamic problems, rather "simple" implicit schemes (like Gauss-Seidel) combined with multigrid result in computationally more efficient solvers (in terms of the CPU time) than, e.g., GMRES [178], [198].

The solution of the governing equations (2.19) does not raise any fundamental difficulties in the case of inviscid or laminar flows. The simulation of turbulent flows, however, presents a significant problem. Despite the performance of modern supercomputers, a direct simulation of turbulence by the time-dependent Navier-Stokes equations (2.19), called the Direct Numerical Simulation (DNS), is still possible only for rather simple flow cases at low Reynolds numbers (\(Re\)). The restrictions of the DNS become quite obvious when recalling that the number of grid points needed for sufficient spatial resolution scales as \(Re^{9/4}\) and the CPU-time as \(Re^{3}\). This does not mean that DNS is completely useless. It is an important tool for understanding the turbulent structures and the laminar-turbulent transition. DNS also plays a vital role in the development and calibration of new or improved turbulence models. However, in engineering applications, the effects of turbulence can be taken into account only approximately, using models of various complexities.

The first level of approximation is reached for the Large-Eddy Simulation (LES) approach. The development of LES is founded on the observation that the small scales of turbulent motion posses a more universal character than the large scales, which transport the turbulent energy. Thus, the idea is to resolve only the large eddies accurately and to approximate the effects of the small scales by relatively simple subgrid-scale models. Since LES requires significantly less grid points than DNS, the investigation of turbulent flows at much higher Reynolds numbers becomes feasible. But because LES is inherently three-dimensional and unsteady, it still remains computationally very demanding. Thus, LES is still far away from becoming an engineering tool. However, LES is well suited for detailed studies of complex flow physics including massively separated unsteady flows, large scale mixing (e.g., fuel and oxidiser), aerodynamic noise, or for the investigation of flow control strategies. LES is also very promising for more accurate computations of flows in combustion chambers or engines, heat transfer and of rotating flows. An overview of research activities in LES was recently published in [210].

The next level of approximation is represented by the so-called Reynolds-Averaged Navier-Stokes equations (RANS). This approach, which was presented by Reynolds in 1895, is based on the decomposition of the flow variables into mean and fluctuating parts, followed by time or ensemble averaging [211] (see also [212], [213]). In cases where the density is not constant, it is advisable to apply the density (mass) weighted or Favre decomposition [214], [215] to the velocity components. Otherwise, the averaged governing equations would become considerably more complicated due to additional correlations involving density fluctuations. It is common to assume that Morkovin's hypothesis [216] is valid, which states that the turbulence structure of boundary layers and wakes is not notably influenced by density fluctuations for Mach numbers below 5.

By inserting the decomposed variables into the Navier-Stokes equations (2.19) and averaging, we obtain formally the same equations for the mean variables with the exception of two additional terms. The tensor of the viscous stresses is extended by one term - the Reynolds-stress tensor [211]

where \(v_i''\), \(v_j''\) denote the density-weighted fluctuating parts of the velocity components \(u,v,w\); \(\overline{\phantom{x}}\) and \(\widetilde{\phantom{x}}\) stand for ensemble and density weighted averaging, respectively. The Reynolds-stress tensor represents the transport of mean momentum due to turbulent fluctuations. Furthermore, the diffusive heat flux \(k\nabla T\) in the energy equation (cf. Eq. (2.8)) is enhanced by the so-called turbulent heat-flux vector [43]

Thus, we can see that the solution of the Reynolds-averaged Navier-Stokes equations requires the modelling of the Reynolds stresses (3.13) and of the turbulent heat flux (3.14). The advantages of this approach are that considerably coarser grids can be used as compared to LES, and that stationary mean solution can be assumed (at least for attached or moderately separated flows). Clearly, both features significantly reduce the computational effort in comparison to LES or even DNS. Therefore, the RANS approach is very popular in engineering applications. Of course, because of the averaging procedure, no detailed information can be obtained about the turbulent structures.

A large variety of turbulence models was devised to close the RANS equations and the research still continues. The models can be divided into first- and second-order closures, respectively.

The most complex, but also the most flexible, are second-order closure models. The Reynolds-Stress Transport (RST) model, which was first proposed by Rotta [217], solves modelled transport equations for the Reynolds-stress tensor. The partial differential equations for the six stress components have to be closed by one additional relation. Usually, an equation for the turbulent dissipation rate is employed. The RST models are able to account for strong nonlocal and history effects. Furthermore, they are able to capture the influence of streamline curvature or system rotation on the turbulent flow.

Closely related to the RST approach are the Algebraic Reynolds-Stress (ARS) models. They can be viewed as a combination of lower level models and the RST approach. The ARS models employ only two transport equations, mostly for the turbulent kinetic energy and the dissipation rate. The components of the Reynolds-stress tensor are related to the transport quantities by non-linear algebraic equations [218]. The ARS approach is capable of predicting rotational turbulent flows and secondary flows in channels with accuracy similar to the RST models. Detailed overviews of the RST and ARS models can be found in [219], [220].

Because of numerical problems with the RST and ARS models, which are primarily caused by the stiffness of the RST and the non-linearity of the ARS equations, first-order closures are more widely used in practice. In these models, the Reynolds stresses are expressed by means of a single scalar value, the so-called turbulent eddy viscosity. This approach is based on the eddy viscosity hypothesis of Boussinesq [221], [222], which assumes a linear relationship between the turbulent shear stress and the mean strain rate, similar to laminar flow. Herewith, the dynamic viscosity \(\mu\) in the viscous stress tensor (2.15) or in the governing equations (2.19) is replaced by the sum of a laminar and a turbulent component

As described earlier, the laminar viscosity is computed, for example, with the aid of the Sutherland formula (2.30). In analogy, the turbulent heat-flux vector (3.14) is modelled as

where \(k_T\) denotes the turbulent thermal conductivity coefficient. Hence, the thermal conductivity coefficient in Eq. (2.24) is evaluated as

The turbulent Prandtl number is in general assumed to be constant in the flow field (\(Pr_T=0.9\) for air). The coefficient of the turbulent eddy viscosity \(\mu_T\) has to be determined with the aid of a turbulence model. The limitations of the eddy viscosity approach are given by the assumption of equilibrium between the turbulence and the mean strain field, and by the independence on system rotation. The accuracy of the eddy-viscosity based models can be significantly improved either by using correction terms [223], [224], or by employing non-linear eddy viscosity approaches [225]-[227].

The first-order closures can be categorised into zero-, one-, and multiple-equation models, corresponding to the number of transport equations they utilise. Within the zero-equation or, as they are also denoted, algebraic models, the turbulent eddy viscosity is calculated from empirical relations, which employ only local mean flow variables. Therefore, no history effects can be simulated, which prevents a reliable prediction of separated flows. The most popular algebraic model, which is still in use for some applications, was developed by Baldwin and Lomax [228].

History effects are taken into account by the one- and two-equation models, where the convection and the diffusion of turbulence is modelled by transport equations. The most widely used one-equation turbulence model is due to Spalart and Allmaras [229], which is based on an eddy-viscosity like variable. The model is numerically very stable and easy to implement on structured as well as on unstructured grids.

In the case of the two-equation models, practically all approaches employ the transport equation for the turbulent kinetic energy. Among a large number of two-equation models, the \(K-\varepsilon\) model of Launder and Spalding [230] and the \(K-\omega\) model of Wilcox [231] are most often used in engineering applications. They offer a reasonable compromise between computational effort and accuracy. An interesting comparison between the Spalart-Allmaras model and various two-equation turbulence models was recently published in [232].

Regardless of the numerical methodology chosen to solve the governing equations (2.19), we have to specify suitable initial and boundary conditions. The initial conditions determine the state of the fluid at the time \(t=0\), or at the first step of an iterative scheme. Clearly, the better (the closer to the solution) the initial guess will be, the faster the final solution will be obtained. Moreover, the probability of breakdown of the numerical solution process will be reduced correspondingly. Therefore, it is important that the initial solution satisfies at least the governing equations and the additional thermodynamic relations. A common practice in external aerodynamics consists of prescribing freestream values of pressure, density and velocity components (given as Mach number, angle of attack and sideslip angle) in the whole flow field. In turbomachinery, it is important to specify the flow directions in the complete domain to one's best knowledge. The same holds also for the pressure field. It is therefore quite worthwhile to employ lower-order approximations (like potential methods) to generate a physically meaningful initial guess.

Any numerical flow simulation considers only a certain part of the physical domain. The truncation of the computational domain creates artificial boundaries, where values of the physical quantities have to be specified. Examples are the farfield boundary in external aerodynamics; the inlet, outlet and the periodic boundary in the case of internal flows; and finally the symmetry plane. The main problem when constructing such boundary conditions is of course that the solution on the truncated domain should stay as close as possible to a solution which would be obtained for the whole physical domain. In the case of the farfield, inlet and outlet boundaries, characteristic boundary conditions [233]-[235] are often used in order to suppress the generation of non-physical disturbances in the flow field. But despite this, the farfield or the inlet and outlet boundaries may still not be placed too close to the object under consideration (wing, blade, etc.). Otherwise, the accuracy of the solution would be reduced. For external flows, when a lifting body is considered, it is possible to correct the flow variables at the farfield boundary using a single vortex centred at the airfoil or the wing [235]-[237]. In this way, the distance between the body and the farfield boundary can be significantly reduced without impairing the solution accuracy, or improving the accuracy for a given outer boundary position [160], [237], [238]. For internal flow problems, formulations for the inlet and outlet boundaries based on linearised Euler equations and Fourier series expansion of the perturbations were developed [239]-[241]. These formulations allow for a very close placement of the inlet and outlet boundaries to a blade without influencing the solution.

A different type of boundary condition is found when the surface of a body is exposed to the fluid. In the case of inviscid flow governed by the Euler equations (2.45), the appropriate boundary condition is to require the flow to be tangential to the surface, i.e.,

By contrast, for the Navier-Stokes equations no relative velocity between the surface and the fluid immediately at the surface is assumed - the so-called no-slip boundary condition

The treatment of walls becomes more involved in cases, where, e.g., a specified wall temperature distribution has to be met, or when the heat radiation has to be taken into account (see, e.g., [242], [243]).

Furthermore, boundary conditions have to be defined for surfaces where different fluids (e.g., air and water) meet together [244]-[247]. But apart from the physical boundary conditions and those imposed by truncating the flow domain, there can be boundaries generated by the numerical solution method itself. These are for example coordinate cuts and block or zonal boundaries [20]-[26].

The correct implementation of boundary condition is the crucial point of every flow solver. Not only the accuracy of the solution depends strongly on a proper physical and numerical treatment of the boundaries, but also the robustness and the convergence speed are considerably influenced. More details of various important boundary conditions are presented in Chapter 8.