Computational tools for supporting the testing of civil aircraft configurations in wind tunnels

ARTICLE IN PRESS

Progress in Aerospace Sciences 44 (2008) 67–120 www.elsevier.com/locate/paerosci

Computational tools for supporting the testing of civil aircraft configurations in wind tunnels S. Bosnyakova,, I. Kursakova, A. Lysenkova, S. Matyasha, S. Mikhailova, V. Vlasenkoa, J. Questb a

TsAGI, 1 Zhukovsky Street, TsAGI, 140180 Zhukovsky, Moscow Region, Russian Federation b ETW, Germany Available online 18 December 2007

Abstract Following a review of well-known CFD packages, this paper presents the mathematical and methodical ideas, which became the basis of the numerical method for a solution of the Reynolds (time)-averaged Navier–Stokes equations (RANS), relevant to supporting testing an aircraft configuration in wind tunnels (WTs). The estimation of finite speed of information propagation due to diffusion is given. Explicit and implicit methods for numerical solution of gasdynamic equations are compared. Stability conditions, taking into account both interaction of convection with diffusion and presence of exponentially varying modes of solution, are considered. Three methods of time stepping organization are described—global, local and fractional time stepping. Numerical boundary conditions for simulation of flow in the WT are proposed. Examples of method implementation for the European Transonic Windtunnel (ETW) case are presented. r 2007 Elsevier Ltd. All rights reserved.

Contents 0. 1.

2.

3.

4.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 1.1. General formulation of conservation laws and numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 1.2. Determination of grid characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 1.3. Favre-averaged Navier–Stokes equations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Selection of the numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 2.1. Selection of an approach to an approximation of diffusive terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 2.2. Methods for solving the small time step problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 2.3. Generalization of the Courant–Friedrichs–Levi condition to tasks with interaction of convection and diffusion. . . . . . 89 2.4. Selection of the method for an approximation of source terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Approximation of Favre averaged Navier–Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.1. Special features of an approximation of convective fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.2. Approximation of diffusive fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.3. Approximation of source terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3.4. Stability condition of the scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Numerical boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.1. General principles of the numerical boundary condition formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.2. Boundary condition on the boundary between adjacent blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3. The boundary condition ‘‘symmetry’’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.4. The boundary condition ‘‘TW_const’’ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

Corresponding author. Tel.: +7 495 556 4323; fax: +7 495 777 6332.

E-mail address: [email protected] (S. Bosnyakov). 0376-0421/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.paerosci.2007.10.003

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

68

5. 6. 7. 8. 9.

4.5. The boundary condition ‘‘Riemann’’. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6. The boundary conditions at the inlet and outlet sections of ETW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical modelling of the ETW test section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The flow behaviour within a single slot channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The test section flow with the installed short axial probe (SAP) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Flow calculations on a wing-body configuration in the test section of ETW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0. Introduction Simulating testing in a wind tunnel (WT) is a complex problem for any mathematical simulation. It consists of elements that essentially influence the test results under investigation, for example, the supporting sting, inclined side walls, diffuser, duct in front of the model, etc. If these parts are neglected, the mathematical model will not be representative. In practice, WTs of different types might be used to test a particular configuration. This may include WTs with open test section (T-101, TsAGI), closed test section (T-16, Arnold Engineering Center), Eiffel’s chamber (SVS-2, TsAGI), perforated walls (T-128, TsAGI) and slotted walls (ETW, Germany). In each case, a specific approach needs to take into account those peculiarities of test equipment and technique used in the experiments. This paper presents a description of calculation algorithms and special boundary conditions that simulate the above peculiarities. Preparing for investigations in a WT is a rather complex and many-sided problem. All stages of experimental research are typically accompanied by complementary numerical work. One may point out, for example, that D’Alambert’s paradox appeared as an outcome of theoretical investigations and has been validated afterwards experimentally. Other notable examples may be mentioned too, but it is beyond the frame of this review. It should be emphasized that a technological revolution in computer capabilities and in calculation methods has enabled the formulation of a problem, which is principally new in complexity. It is the problem of creating the ‘‘electronic wind tunnel (EWT)’’ [1–5]. It is well known that codes tackling 3D problems often operate at the limits of the available computing capacity. This usually results in a computer methodology, being closely linked to an available computer system for an optimization of its resources. Four groups of parameters probably affecting the results are: (1) the calculation algorithm, (2) professional programming, (3) code optimization and (4) parallel computing. All these factors are of equivalent importance, but the last one requires particular attention. Some approaches to parallel computing are well known, for example, the organization of shared memory. In this case, all processes use common RAM simultaneously. They asynchronously ‘‘call for’’ common memory both with reading and writing requests. This model is named as Parallel Random Access Machine

103 103 104 107 109 113 117 119 119

(PRAM). Any processor in an ‘‘ideal’’ PRAM may address to any cell of RAM during the same time. In practice, the access time to RAM for a real computer is non-uniform. So, a uniform memory access (UMA) system is used in practice. Usually this system is realized for 2- or 4-processor computer architectures. The number of computers increases in the case of not uniform memory access (NUMA) ideology or in the case of computers with Common Distributed Memory. In such cases, the memory belongs to different processors but it is connected with each other using a local unibus. More recently, the computer technique underwent sweeping changes. The power of personal computers has achieved that of work stations. Modern clusters with hundreds and thousands of processors have appeared. An essential factor is a huge diminution in system cost. It permits the use of clusters in individual WT laboratories and the integration of computational investigations in the technological loop of experiments. Ideally, each WT may be supplemented by its ‘‘Electronic analogue’’. One may refer to at least three areas of applying such ‘‘analogue’’ work in practice: (1) at the stage of the preparation of the experiment; (2) during the primary processing of experimental data; (3) at the stage of post-processing and creating mathematical models of the investigated phenomena. In many cases, it is useful to simulate the flow in an empty WT taking into account its specific geometry. Mathematical models of describing the phenomena in WTs are developed using different physical approaches. Presently, the nonlinear model of viscous flow is prevalent. For its description, a Reynolds-averaged Navier-Stokes equation system (RANS) [6] is used. A known turbulence model is implemented to close this equation system. It is common understanding that the scientific development of methods using RANS is more or less completed nowadays. This does not mean that further work on this subject is no longer progressing. Nevertheless, the main efforts of scientific schools are now turned towards developments of LES because the trend in physical realism [7] is dependent on these methods. Indeed, published differential turbulence models [8], as a rule, are adjusted to solve a narrow class of problems and may be used only in the framework of a zone-wise approach. For example, the wellknown model (ko) [9] provides acceptable results in nearwall zones, where a fully developed turbulent boundary layer exists. The other popular turbulence model (ke)

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Nomenclature t time xk (x1 ¼ x, x2 ¼ y, x3 ¼ z) Cartesian coordinates ~ r ¼ ðx; y; zÞT radius-vector of a point r density uk (u1 ¼ u, u2 ¼ v, u3 ¼ w) velocity components ~ ¼ ðu; v; wÞ velocity vector V p static pressure T static temperature E total energy per unit mass of gas M Mach number c speed of sound q parameter proportional to characteristic value of velocity turbulent fluctuations o parameter proportional to characteristic frequency of turbulent fluctuations m dynamic coefficient of molecular viscosity tik element of viscous stress tensor qk heat flux along xk-axis due to heat conductivity Cp specific heat per unit mass of gas at constant pressure mT dynamic coefficient of turbulent viscosity RikRik element of Reynolds stress tensor sk turbulent heat flux along xk axis ~ ¼ ðT; u; v; w; p; q; oÞ full vector of parameters P f~ ¼ ðT; u; v; w; q; oÞ vector of gasdynamic parameters without pressure ~ ¼ ðr; ru; rv; rw; rE; rq; roÞ vector of conservative U variables ~ ~ through the cell border F vector of fluxes of U L length t value of time step h grid step in space (in 1D tasks) ~ s ¼ ðsx ; sy ; sz Þ vector of cell border area V volume or absolute value of velocity vector ~ x ¼ ðxi ; xj ; xk Þ coordinates along grid lines ~ n ¼ ðnx ; ny ; nz Þ unit normal vector

z4 ¼ q, z5 ¼R o turbulence parameters R z6 ¼ V n þ dp z7 ¼ V n  dp rc ; rc acoustic Riemann invariants Re Reynolds number Pr Prandtl number CFL Courant number Cdiff analog of Courant number for diffusion a convection velocity in model equation D diffusion coefficient in model equation g ¼ 1.4 specific heat ratio for air RE287 J/(kg K) gas constant for air Upper indices n

[10–12] is adjusted for use in mixing layers and its implementation outside these regions may cause errors. A successful solution is by a combined approach, for example, SST [13], where a smooth transition from one turbulence model to another, depending on the ‘‘remoteness’’ from solid walls, is used. A special problem is the calculation of the position of the laminar–turbulent transition zone. Recently, a new SST model has been proposed that contains additional equations to solve this problem [14].

number of time step

Lower indices parameters along vector ~ n or along n-axis parameters along directions, which are perpendicular to vector ~ n or to n-axis (i, j, k) cell centre numeration (i71/2, j, k), etc. numeration of cell borders n t 1, t 2

Functions and operators used in the text 1.

2.

3.

4.

Riemann invariants for 1D flow along n-axis z1 ¼ Cp ln TR ln p entropy z2 ¼ Vt1, z3 ¼ Vt2 velocity components, which are perpendicular to n-axis

69

5.

Linear interpolation between the values aL ¼ a(xL) and aR ¼ a(xR) for the point xA[XL; XR], xxL ¼ DL, xRx ¼ DR: interpolation D a þD a a a ðaL ; aR ; DL ; DR Þ ¼ RDL þDL R ¼ aL þ DL DRþDL : L R L R Linear extrapolation through the values a1 ¼ a(x1) and a2 ¼ a(x2) for the point x, if xx1 ¼ D1, x1x2 ¼ D2: extrapolation ða1 ; a2 ; a a D1 ; D2 Þ ¼ a1 þ D1 1D 2 : 2 Increment of parameter a between the border centre of (i+1/2, j, k) cell and centre of (i, j, k) cell (see Fig. 11): DLiþ1=2 a  aiþ1=2;j;k  ai;j;k : Increment of parameter a between the centre of (i, j, k) cell and border centre of (i1/2, j, k) cell (see Fig. 11): DR i1=2 a  ai;j;k  ai1=2;j;k : Increment of parameter a between centres of (i+1, j, k) and (i, j, k) cells:Diþ1=2 a  aiþ1;j;k  ai;j;k ¼ DLiþ1=2 a þ DR iþ1=2 a:

A common point of the above-mentioned methods is their low accuracy in regions of separated flow. Many papers are devoted to this problem—for example [8,15–17]. These reveal that specialized turbulence models with an enlarged number of parameters permit an improvement of convergence of calculated and experimental data. In a number of publications, RSM models [18] are recommended for application. Those models give for example good results in two special cases: (1) an oblique shock caused by a wedge and (2) an oblique shock reflection from

ARTICLE IN PRESS 70

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

a channel wall. A serious problem arises during investigations of ‘‘hot’’ jets from nozzles in WTs. Already in the simplest case of an axis symmetric jet, ‘‘cold’’ turbulence model application results in an essential disagreement of calculated and experimental data [19]. This does not enable a correct estimation of propulsion performance and acoustic characteristics of real nozzles [20]. The list of possible problems may be expanded, but this paper is devoted to define, which practical tasks may be successfully solved in a framework of RANS applications. However, this question does not have a simple answer. The application of different industrial codes (for example [21–30]), does not reveal answers to many simple questions regarding the accuracy of obtained results. It is due to uncertainties that arise in the determination of computational domain boundaries and from the necessity to formulate additional boundary conditions for different turbulence models. There is a paradox that an investigator can obtain widely differing results by varying boundary conditions or adjusting turbulence model coefficients. In some cases such results are in good agreement with experimental data due to procedures named ‘‘tuning’’. In this case, the experimental study has to be performed in advance. Experimentalist may consider that such code ‘‘adjustment’’ is an elementary garbling of the results. The only possible solution of this problem is to perform a direct simulation of boundary conditions and to achieve the correspondence of turbulence levels and scales in zones by mounting the models indirectly. When an experiment in a WT is simulated, the influence of both WT walls and supporting systems need to be taken into account. Turbulence levels and scales at the entrance to a WT test section have to be determined. It is obvious that the complexity of the task in such a formulation essentially increases because the user spends considerable time on grid construction and because additional computer resources are spent for calculating the flow upstream and downstream of the model. But it provides a possibility of predicting investigated effects adequately at the stage of the preliminary preparation of an experiment. Aircraft development is a many-sided process, where specialists of different disciplines are typically involved. Usually, it includes several levels. At the stage of the preliminary development, ‘‘fast’’ codes based on simplified (aerodynamic) methods are used. An advantage of integrated approaches is the quickness of their solution but its drawback is the low accuracy of the result. A solution is to use a EWT-like code and to replace these simplified calculations by codes realized by RANS methods. On the one hand, it results in time-consuming calculation that is unacceptable for a practical application. On the other hand, it provides means of specific preliminary data. Let us consider that an aircraft project is ready and it is time to verify ideas, which are formulated in the project. The preparation of experiments needs time. Usually, model development and manufacturing typically will require 6 months. During this period, a detailed

Fig. 1. Applicability domain of EWT for practical purposes.

calculation of the flow around the model both in WT and in free flow is possible. The applicable domain of such a calculation results is sufficiently wide enough. Fig. 1 shows a block scheme, where links arising during the organization of collaboration in the framework of new aircraft development project are presented. An item ‘‘EWT’’ is placed into the centre of the scheme. It is connected with the following items: (1) ‘‘aerodynamics’’, (2) ‘‘dynamics’’, (3) ‘‘strength’’, (4) ‘‘gas dynamics of engine’’ and (5) ‘‘avionics’’. Some links in Fig. 1 are obvious but others may require some explanation, for example, the connection ‘‘EWT’’‘‘avionics’’. It is known that one practical problem is how to optimize the location of static and total pressure probes, vane sensors and other devices that are necessary to control the aircraft. This problem may be solved in the framework of an EWT implementation. A detailed calculation of the flow field around the aircraft nose at all main regimes provides the necessary information to make the best decision. Another important example is an investigation of an engine thrust reversal system. This needs the collaboration of different experts. For example, to avoid the penetration of hot jets into the inlet it is necessary to close the reverse system by flaps in a timely fashion. To predict this time it is possible to use EWT and to select suitable positions for sensors which register the presence of the jet being sucked into the engine. This job is for the experts of the engine department. Aerodynamicists also take part in this job and propose wing–pylon configurations to provide the maximal deflection of jets from the engine. Structural experts calculate loads to avoid destructions of the structure. Flight dynamics experts try to avoid an aircraft deviation from the runway as well as an aircraft turnover. It is obvious that all these groups of specialists need reliable data that may be provided using EWT. There are other possible areas for applying the codes prepared in this framework of RANS methodology. Fig. 2 presents a block scheme that demonstrates how EWT may be used for scientific research purposes. To provide high levels of flexibility the code needs to be easily to be modified and to contain different mathematical models, boundary conditions and methods. Usually it is necessary to include internal code templates providing important options as: (1) aircraft shape optimization, (2) wing–pylon–engine integration, (3) engine–fuselage integration; (4) inlet–engine–nozzle integration; (5) preliminary fast estimation of aerodynamic forces and momentums; (6) multi-disciplinary phenomena. Special

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

attention is usually paid to the accuracy of the static pressure calculation on surfaces because it permits both the determination of integral values of aerodynamic coefficients and the solution of aero-elastic problems. In this case, the RANS equation may be used because of many factors for example, the necessity to take into account the influence of separated flow zones on local peaks of pressure and temperature. In addition, the flow in separated zones is unsteady and so needs to be taken into account for problems of dynamic strength and flight dynamics. EWT may also be applied as an approach to fill in missing experimental data. It is useful, when the experiment has been performed using newly emerging methods, such as pressure sensitive paints, PSP [31]. As a positive experience, let us consider a case when complete surface distribution of the static pressure obtained in PSP is used to formulate boundary conditions in the EWT. The calculation performed in this formulation may be compared with the calculation performed using the traditional method. Displacement of shock waves and rarefaction zones in flow field may be treated as a real contribution to information in 3D space that is absent in the PSP method. Another example is linked to particular methodologies of experiments in WTs. For example the unique WT T-128 (TsAGI) permits the adaptation of the test section wall porosity during a run in order to minimize the influence of the WT walls on the experimental data. Let us consider a possible method of the organization of experimental data.

Fig. 2. Applicability domain of EWT for scientific purposes.

Fig. 3. Codes and used computers.

71

At the first stage, a calculation of the flow around the aircraft model in free flow is performed and CFD data are interpolated to the imaginary walls of the WT test section. Using such information and technological sheets, the operator provides a coincidence of the WT wall static pressure distributions with available CFD data (it may be performed using a special code developed for T-128). In this approach the influence of WT walls on experimental data are essentially diminished. Another example shows how the influence of a model sting is taken into account using a EWT code. It this case, two CFD ‘‘runs’’ are required: one calculating the total configuration (WT+model+sting) and another one calculating it without the sting. The difference between the obtained solutions represents the sting effect on test results. The numerical experiment implies that numerous CFD ‘‘runs’’ are performed to obtain parametrical data. A special attention, in this case, is paid to a computation system that may consist of one or several computers. Our experience shows that real work is only possible using supercomputers (clusters). An analysis of information presented in [32] shows that supercomputers of two types (CRAY and NEC) are used across the EU. Presently, the CRAY is widely going to be replaced by the sixth series of NEC. For example, code ElsA, which may be used to simulate flow in WTs, is prepared for an installation using CRAY or NEC (see Fig. 3). Versions of this code have also been developed to work both on clusters and personal computers. The FLOWer and TAU codes have similar properties. On the other hand, the FASTFLO code is mainly oriented to computers of the NEC series. Obviously, the computer selection is often dominated by the financial resources of the individual applier. Typically, universities and small companies perform the calculations using moderate-size clusters. Usually, EWT-like codes consist of four main tasks (see Fig. 4). The geometry is closely linked to the key task of numerical grid generation. Our experience shows [33] that all the main numerical methods are grid-dependent ones. For example, a scheme of low approximation order may result in more high-quality solutions for an adapted grid than a scheme of higher order for a grid of lower quality (for a similar number of nodes). Three types of grids are typically used: (1) structured; (2) unstructured; and (3) combined (hybrid). Fig. 5 presents information about the grids used in the considered codes. It is obvious that users prefer structured and combined grids because they may provide a higher quality of results. For example in the boundary layer region structured or

Fig. 4. Typical task structure of a CFD code.

ARTICLE IN PRESS 72

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 5. Typically used grid types.

Fig. 6. Methodologies used in grid development.

combined grids offer essential advantages due to the adapted direction of the gridline. Combined grids are constructed as follows: structured grids near the surface and unstructured grids in the free flow (code CANARI and SAUNA). A lot of codes use also unstructured grids (Airplane, EDGE, TAU1). This is due to ease a 3D grid generation for complicated applications. It is known that semi-automatic generators for unstructured grids are included into many engineering systems and work using a principle of ‘‘a mouse click’’. On the other hand, the structured grid development is a rather complicated task requiring intensive manual work. Fig. 6 schematically presents techniques and methods that are used by code developers to generate structured grids. First of all, the technology using chimera characteristic grids is of interest. It permits describing peculiarities of individual elements of the aircraft geometry in a most favourable way. For example, a grid may essentially be simplified when it is necessary to calculate flow around a wing with deflected flaps. A method of locally compressing the grid near body edges or around points of high flow gradients is widely used resulting in an essential increase of the quality of the numerical solution. Fig. 7 presents an example of a polar for a typical civil aircraft fuselage for two calculations: (1) with a uniform grid and (2) with a grid compressing nodes in the shock region. It is obvious that the polar becomes more accurate in the case of node compression. It should be mentioned however that, in the case of applying explicit schemes, a local compression 1 TAU is used with prismatic grids covering the boundary layer plus unstructured grids for the outer flow (low speed high lift applications).

Fig. 7. Influence of local grid compression near shock.

results in negative ‘‘follow-on’’ effects. There are losses in the calculation accuracy at large cells that remain after realization of local compression procedures. It is known [34] that an explicit scheme has its best accuracy if the Courant numbers are about CFL ¼ 1. In the case of a grid compression, the Courant number is defined by the minimum cell. At that, large cells are calculated using Courant numbers that are essentially lower than 1 and, as a result, the accuracy at these cells decreases. Fig. 6 shows the local compression technology being used in 7 codes among the 20 ones considered. In many cases, rotating turbine blades, turboprop engine propellers and rotors of helicopters have to be simulated. A methodology of ‘‘rotated matrices’’ is used for those purposes. It is used in 5 codes among those presented in

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 6. The extremely complicated problem of merging numerical solutions obtained from different grids is solved mathematically. Then in the merging regions the approximation order of the resulting solution may be lost. In addition, any convergence is not guaranteed, when the velocity of shifting the boundary is comparable to the velocity of the propagation of the perturbation. One of the most important problems in the calculation of wings and engines is the application of ‘‘C’’-type grids. It has been shown by many authors [35] that grids of other types essentially decrease the quality of the solution. In the simple cases of isolated wings or engines, a ‘‘C’’-type grid is easily developed and appropriate. The main problem arises in the case of an aircraft configuration in which it is necessary to merge grids of different types, for example ‘‘H’’ and ‘‘C’’ grids. Such joining may be performed in an automatic or manual way. Fig. 6 identifies codes that realize automatic algorithms capable of solving this problem with a good quality. Grid adaptability means that nodes may be rebuilt at each calculation step for the purpose of automatic compression, where it is necessary to solve the problem. An advantage of such an approach is an increase in accuracy but the drawback is the simultaneous increase in calculation time and computer resources. Adaptation does need an adjustment for each definite case. For simple cases (for example, shock-induced separation on an airfoil), they provide excellent results. In complicated cases, however, failures are possible. For example, simultaneous adaptation for shock and rarefaction waves due to a separated boundary layer causes complications. Hence, calculations using non-adaptive grids that are compressed in advance are presently more popular. The Lagrange–Euler (ALE condition) is the mostly used approach. Practically, it permits the calculation of moving objects in disturbed flows. As an example, a problem of solid object ingestion into an engine during take-off acceleration may be amenable to such an approach. The computational domain is divided into blocks. Some of

73

these blocks are movable and attached to the interested solid object. Special attention in ALE realization is paid to boundary conditions to avoid order of approximation losses near the joint of mobile and fixed blocks. Than a technology of Euler and Lagrange grid intersection followed by interpolation of the calculation would partly be used. An analysis shows that purely block-structured grids are rarely used now, see Fig. 6. In part, it is due to the highly laborious nature of their construction. The problem is simplified, if the condition of joining the adjacent blocks is not satisfied (one permits intersection, as it is shown above). But this results in an essential loss of calculation accuracy. On the other hand, a positive effect of using the blocks is obvious in the calculation capabilities of complicated geometries. Blocks enable separation of all the main peculiarities resulting in high accuracy. In addition, a block-structured grid enables a zonal approach to be realized as mentioned in the introduction. Indeed, different turbulence models or different numerical methods may be used in different blocks (for example, the Godunov method [34] in one zone and the Jameson’s potential method [36] in the other one). A block-structured grid provides essential advantages in the parallelization of problems using a multi-processor computer or cluster. Then, blocks are natural parts of the problem, which are easily stored on different processors. The process of information exchanges comes from the formulation of boundary conditions at each calculation step and does not need data transfer over the network. As a rule, information exchange takes place in the near-boundary zones of appropriate blocks. An analysis shows that explicit numerical schemes are used for the solution of steady problems in many cases (see Fig. 8). This may be explained, first of all, by the welldeveloped and investigated methodology for obtaining high-quality monotonic solutions in complicated situations, where the calculation accuracy depends on many geometrical and physical factors. Special procedures are of

Fig. 8. Types of used numerical schemes in steady problems.

Fig. 9. Additional methods used in explicit schemes.

ARTICLE IN PRESS 74

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

great importance (for example, implicit smoothing algorithms, see Fig. 9). They permit the elimination of oscillations due to the use of non-monotonic schemes, for example, schemes based on Runge–Kutta integration with artificial viscosity. The main problem of explicit schemes is the lengthy calculation necessary for realistic test cases. A number of methods are used to accelerate the process, but algorithms, which are based on the ideas of multi-grid, are of widest application. These approaches are based on the property of waves (that arise in the discrepancy of current flow fields and of given boundary conditions) to propagate with speeds depending on cell sizes. Long waves are quickly propagated on a large-cell grid; small cells propagate on the short-wave part of the spectrum. Our experience shows that the problem may be accelerated by orders of magnitudes (with comparable accuracy). It may be achieved in the case of an optimal organization of algorithms. It should be noted however that practical application of multi-grid algorithms is connected with serious problems. This method demands a higher requirement for RAM because, in practice, it is necessary to store data for all grids used in this multi-grid approach. Sometimes, multi-grid algorithms operate unstable, for example, in the case of using the function ‘‘minmod’’. The main problem is the 3D reconstruction [37] of algorithms sometimes violating the uniformity of the used method. Fig. 10 shows codes that perform simulations of laminar–turbulent transition, a difficult task which is not finally resolved even today. It has been mentioned above that, in a number of cases, it is necessary to specify a priori the laminar–turbulent transition of the boundary layer because it is defined by artificial factors that cannot be calculated. The drawback of this approach is the need for experimental data obtained ‘‘a priori’’. So, it may be used in the case of simple and known geometries but, in the case of a generally new design, there is no information and thus it will be difficult to predict the laminar–turbulent transition. Approaches relevant to the calculation of engine influence on aircraft aerodynamics, for example ‘‘activator disks’’, are extremely useful. There exist various algorithms characterized by the complexity of mathematical models. Simple variants account for energy and impulse growth due

to engine operations; complicated ones simulate swirl and the presence of unburned fuel components in the exhaust. Any approach may be used, depending on the complexity of the solved task. The results of flow calculations in a WT are presented in various articles. For example, the flow around a half-model of a civil aircraft consideration is considered in European Transonic Windtunnel (ETW) [38]. A peculiarity of this calculation is the deployed high-lift devices. The model is attached to the ceiling of the WT test section (all slots in the walls were closed). A spacer called ‘‘peniche’’ is located between the fuselage and the ceiling. The peniche shifts the model out of the wall boundary layer to reduce its effect on the flow field development. The calculation is performed using an unstructured prismatic grid within the boundary layer; the algorithms of the grid construction inside the boundary layer and in the inviscid main flow are different. The grid generator TRITET [39] is used. Calculated results and experimental data with and without applied wall interference corrections are compared. The conclusion is that the calculation provides a good description of the flow field around model in the WT. Article [40] contains an experience of flow simulation in two WTs-DNW-NWB and ETW. The calculations were performed using the TAU code [41] featuring unstructured grids. Solutions were obtained using an implicit numerical scheme. To close the RANS equation system, the SST turbulence model [13] is used. It is partly adjusted to solve problems including boundary layer separation. Effects of WT walls and model supports (stings) on test results have been investigated. Algorithms for experimental result correction using numerical solutions have been proposed. A model of a transport aircraft ALVAST (similar to an Airbus A320) in landing configuration is considered. In addition, special attention is paid to methodical calculations of the flow around an axis-symmetrical body named SAP in ETW (see chapter 7) to permit investigations of the influence of WT elements on the flow around the model [42]. A general problem in that context is the achievement of a pre-defined Mach number at the model position. In WTs it is an outcome of comprehensive calibrations to reference it to a measurement far upstream of the test article. A Mach number correction method taking into account experimental results is also proposed in [42]. Collercandy [43] has made an important contribution to the creation of

Fig. 10. Simulated physical processes.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

a mathematical model of ETW. He performed the calculations of flow around the HIRETT model in ETW using a RANS code with a structured grid. Unfortunately, these papers are restricted. In conclusion, it should be noted that considerable efforts for creating codes to simulate flows in WTs have been spent in Sweden. For example, a special attention has been paid in Ref. [44] to physical aspects of flow in an axis-symmetrical WT with slotted walls. At the same time, a mathematical model is developed using an inviscid gas model [45]. The results obtained are in fair correspondence with the experiments; it confirms an assumption formulated above that, in a number of cases, simplified methods may be used to create mathematical models of WTs. The following four sections provide: (i) the basic equations, (ii) the choice of the numerical scheme, (iii) the approximation of Favre-averaged Navier–Stokes equations and (iv) the numerical boundary conditions developed in the TsAGI method for the calculations of 3D flows. The remaining sections deal with how the method was used to tackle the problem of flow simulation for a civil aircraft configuration mounted in the test section of a slotted wall transonic WT. 1. Basic equations 1.1. General formulation of conservation laws and numerical scheme A numerical scheme is constructed in the framework of the finite-volume approach. Considered is an arbitrary flow of a gas with a fixed volume V inside this flow. Surface S of this volume is stationary in time, and the gas flows through it freely. The conservation laws of mass, momentum and energy and turbulence parameters for this volume are then written. Introducing a quantity of an arbitrary physical value in unit volume of gas, the integral conservation law of this physical value may be represented as follows: Z I Z q ~a d~ a dV ¼  F sþ W a dV . (1) qt V V S

On the left-hand side of this equation we see the variation rate of the parameter a and the total quantity in the volume V. The right-hand side describes the reasons of this variation: a (1) fluxes of a through the surface S (F~ is flux of a through unit area per unit time, d~ s is vector of area element, which is perpendicular to the surface); (2) local sources of a (Wa is source term that describes the rate of production or consumption of a value due to local sources).

By applying the Gauss–Ostrogradsky theorem to the integral and moving the time derivative under the integral

75

(because V ¼ const), in the limit V-0 one may obtain a differential equation that expresses this conservation law: a

qa qF ax qF y qF az þ þ þ ¼ W a. qt qx qy qz

(2)

Here F ax ; F ay ; F az are fluxes of a through the unit area, oriented perpendicular to the x, y, z axes. The applicability of the differential Eq. (2) is tighter than the applicability of the integral law (1), because (2) cannot be applied to discontinuous flows. We shall calculate non-stationary flow using the time marching approach. In this approach, the initial flow field (in the time moment t0) is known, and solution for t4t0 is obtained by sequential transitions from a known time layer tn, n ¼ 0,1,2,y, to the next time layer tn+1 ¼ tn+tn. To obtain a stationary solution, we shall use the stabilization method, when an initial approximation is given, and a stationary solution is obtained as a limit of non-stationary adaptation of the flow to the given stationary boundary conditions. We use regular grids with hexahedral cells. During computation, flow parameters are stored in the centres of grid cells. Let us consider grid cell (i, j, k) (Fig. 11). We introduce the following definitions: Z 1 ani;j;k ¼ að~ r; tn Þ dV V i;j;k V i;j;k average value of a over the cell in time moment tn, Z tnþ1 Z 1 a W i;j;k ¼ n W a ð~ r; tÞ dV dt t V i;j;k tn V i;j;k average value of Wa over the cell and over the time step tn, Z n nZ 1 t þt ~a d~ F s dt, ðF a Þiþ1=2 ¼ n t tn S iþ1=2;j;k a

~ ¼ ðF a ; F a ; F a Þ—average value of the a flux where F x y z through the cell border (i+1/2, j, k) over the time step tn. It is worth noting that fluxes F ax ; F ay ; F az are related to unit area; on the contrary, ðF a Þiþ1=2 is the flux through the whole cell border with area Si+1/2, j, k. By integration of Eq. (1) over time from the moment tn till the moment tn+1, one may obtain the following approximation of the conservation law (1): n anþ1 i;j;k ¼ ai;j;k 

tn ½ðF aiþ1=2  F ai1=2 Þ þ ðF ajþ1=2  F aj1=2 Þ V i;j;k

þ ðF akþ1=2  F ak1=2 Þ þ tn W ai;j;k .

ð3Þ

We shall consider the following conservation laws: (1) conservation law of mass (a ¼ r); (2) conservation law for each of the three velocity components (a ¼ ru; rv; rw); (3) conservation law of energy (a ¼ rE); (4) conservation laws for the parameters of the turbulence model (in this paper turbulence parameters are q and o, see Section 1.3). Combining approximations (3) for all conservation laws in

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

76

Fig. 11. A cell of a computational grid and a map of indices.

one system and introducing the vector of conservative ~ ¼ ðr; ru; rv; rw; rE; rq; roÞ, vectors of fluxes variables U ~ ~j1=2 ; F ~k1=2 , and of U through the cell borders F~i1=2 ; F ~ vector of source terms W , we obtain a general formulation of our numerical scheme:    tn h ~ ~i1=2 þ F ~jþ1=2  F~j1=2 F iþ1=2  F V i;j;k  i ~ i;j;k . ð4Þ þ F~kþ1=2  F~k1=2 þ tn W n

~ ~nþ1 ¼ U U i;j;k i;j;k



n

~ is known (it may be calculated In Eq. (4) the vector U i;j;k through the flow parameters at the known time layer n). nþ1 ~ Vector U i;j;k needs to be found. For this purpose, one should define the methods for calculation of fluxes ~i1=2 ; F~j1=2 ; F~k1=2 and of source terms W ~ . These F methods will be defined in the following sections.

Increments of coordinates along grid lines are defined as follows:   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   DLiþ1=2 xi ¼ DLiþ1=2~ r ¼ ðDLiþ1=2 xÞ2 þ ðDLiþ1=2 yÞ2 þ ðDLiþ1=2 zÞ2 ,   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  R  2 2 2 R R DR r ¼ ðDR iþ1=2 xÞ þ ðDiþ1=2 yÞ þ ðDiþ1=2 zÞ , iþ1=2 xi ¼ Diþ1=2~ Diþ1=2 xi ¼ DLiþ1=2 xi þ DR iþ1=2 xi .

(As an example, all formulas in this paper will be given for grid line i and for cell border (i+1/2, j, k). Along other grid lines and for other cell borders, analogous formulas are used.) To relate the derivatives of flow parameters along grid lines with the derivatives along Cartesian axes, one should know the Jacobian matrix of coordinate transformation: 31 2 3 2 J¼

1.2. Determination of grid characteristics Let us also describe how we calculate characteristics of the computational grid, which are used in our numerical method. Centres of grid cells are determined as the arithmetic mean of the cell vertices, and centres of cell borders, as the arithmetic mean of the border vertices. Advantages of such an approach are small quantities of calculations, stability of algorithm on poor quality grids (with strongly elongated, skewed cells or cells with degenerated borders), and also the fact that the centre of each cell is placed at the middle point of segments, which connect centres of opposite borders. In the formulation of numerical method, the coordinate frame, which is connected with grid lines, is introduced (Fig. 11).

ð5Þ

qxi qx 6 qx 6 j 6 qx 4 qxk qx

qxi qy qxj qy qxk qy

qxi qz qxj qz qxk qz

qx qx

i 7 6 qy 7 6 7¼6 qxi 5 6 4 qz

qxi

qx qxj

qx qxk

qy qxj

qy qxk

qz qxj

qz qxk

7 7 7 . 7 5

(6)

This matrix should be known both in the cell centres and in the centres of cell borders. In the cell centres the metric coefficients from (6) are calculated by the formula: !



q~ r qxi



¼ interpolation i;j;k

¼

DLiþ1=2~ r DLiþ1=2 xi

L DR r DLiþ1=2~ r DR i1=2~ i1=2 xi Diþ1=2 xi ; ; ; L 2 2 DR i1=2 xi Diþ1=2 xi

.

ð7Þ

(Here it is taken into account that DR DLiþ1=2~ r r i1=2~ ¼ , R L Di1=2 xi Diþ1=2 xi because the cell centre is placed at the middle point of the segment that connects centres of opposite borders.)

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

In the centres of cell borders, which are placed inside the computational domain (not at its boundary), the metric coefficients are calculated as follows (see Fig. 11):

77

where averaging interval T should be much larger than characteristic time of turbulent fluctuations, but much smaller than the characteristic time of the mean flow

!  DLiþ1=2~ r DR r DLiþ1=2 xi DR q~ r iþ1=2~ iþ1=2 xi ; ¼ interpolation L ; ; , qxi iþ1=2;j;k 2 2 Diþ1=2 xi DR iþ1=2 xi     ~ riþ1=2;j1=2;k1=2 þ ~ riþ1=2;j1=2;kþ1=2 riþ1=2;jþ1=2;k1=2  ~ riþ1=2;jþ1=2;kþ1=2  ~ q~ r   , ¼  qxj iþ1=2;j;k ~ riþ1=2;jþ1=2;kþ1=2  ~ riþ1=2;j1=2;k1=2 þ ~ riþ1=2;j1=2;kþ1=2  riþ1=2;jþ1=2;k1=2  ~     ~ riþ1=2;j1=2;k1=2 þ ~ riþ1=2;jþ1=2;k1=2 riþ1=2;j1=2;kþ1=2  ~ riþ1=2;jþ1=2;kþ1=2  ~ q~ r   . ¼  qxk iþ1=2;j;k ~ riþ1=2;jþ1=2;kþ1=2  ~ riþ1=2;j1=2;k1=2 þ ~ riþ1=2;jþ1=2;k1=2  riþ1=2;j1=2;kþ1=2  ~ 

In the centres of cell borders, which are placed at the boundary (these borders have numbers i ¼ 1/2 and Ni+1/ 2, where Ni is quantity of inner cells along i-axis in the given block of computational domain), the following formulas are used: 

ð8Þ

variation. For convenient description of compressible turbulent flow, in expansion of terms of time-averaged equations the Favre averaging is used for all gas parameters, except for r and p. Favre averaging is defined by the relation a~ ¼ ra=r¯ [46]. As a result of averaging,

!  R L DR r DL3=2~ r DR q~ r 1=2~ 1=2 xi D1=2 xi þ D3=2 xi ; ¼ extrapolation R ; ; qxi 1=2;j;k 2 2 D1=2 xi DL3=2 xi ¼



DR r 1=2~ DR 1=2 xi

,

 DLN i þ1=2~ r q~ r ¼ L , qxi N i þ1=2;j;k DN i þ1=2 xi

where it is taken into account that DR r 1=2~ R D1=2 xi

¼

DL3=2~ r , L D3=2 xi

because the cell centre is placed at the middle point of segment that connects centres of opposite borders. The exception is boundary of type ‘‘Joint’’ (see Section 4.2), where the same formulas (8) are used, as for inner borders of cells. For derivatives along directions, which are tangential to the boundary of computational domain, on the borders i ¼ 1/2 and Ni+1/2 for each type of boundary the same formulas (8) are used, as for inner borders of cells. It is worth noting that (7)–(9) gives values of metric coefficients with second order of approximation.

ð9Þ

additional terms arise, which are connected with fluctuations ~ For closure of the equation system for timea00 ¼ a  a. averaged flow, semi-empirical turbulence models are used; they should relate in some way to the additional terms with parameters of mean flow. In this paper, as an example only one turbulence model is used, namely the (qo)-model of Coakley [47,48]. This model uses two differential equations for the following parameters of turbulence: pffiffiffi g 00 2 þ vf 00 2 þ w 00 2 Þ=2 is kinetic (1) q ¼ k~, where k~ ¼ ðuf energy of turbulence. Parameter q is proportional to characteristic value of velocity fluctuations pffiffiffi ~ hV 0 i ¼ 2k; ~ where  ¼ ðP tij qu00 =qxj Þ=r¯ is the dissipation (2) o ¼ =k, i i;j rate of turbulence kinetic energy. Parameter o has dimension [Hz]; it should be proportional to the characteristic frequency of turbulent fluctuations

1.3. Favre-averaged Navier–Stokes equations Unlike in the previous paper [4], a time-averaged Navier–Stokes equation system is considered. Averaging is performed in accordance with the formula Z 1 tþT=2 aðxi ; tÞ ¼ aðxi ; tÞ dt, T tT=2

Z

Z

1

fEðf Þ df

hf i ¼ 0

1

1

Eðf Þ df

,

0

where E(f) is the spectral density of the turbulent fluctuation energy. The coefficient between o and /fS, generally speaking, depends on the considered class of flows.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

78

Below we shall omit the signs of averaging. All flow parameters in formulas will be averaged in accordance with the rules of Favre averaging. The (qo)-turbulence model belongs to the class of semi-empirical models based on Boussinesq hypothesis [49] about similarity between the processes of turbulent transport and processes of molecular transport. This hypothesis enables to use for turbulent fluxes the same expressions as for molecular fluxes, but with a replacement of the molecular viscosity coefficient m by the turbulent viscosity coefficient mT. The coefficient mT in (qo)-model is expressed through q and o using a modification of the Prandtl–Kolmogorov formula [8]: mT ¼ C m F wall F Dash r

q2 . o

(10)

Here Cm ¼ 0.09; Fwall ¼ 1exp(0.02rqyW/m) is a wall function that takes into account the influence of molecular viscosity on turbulent fluctuations in the viscous sublayer and in the buffer region of the boundary layer (yW is distance from the given point to the nearest solid wall). Many methods are available to take into account effects of compressibility (see e.g. [8]). In this work the simplest model will be considered, which is based on the model [50] for description of compressibility effects in the mixing layers. In this model, an additional function is introduced in formula (10) for turbulent viscosity: F Dash ¼ 0:25 þ 0:75½1 þ expð24:73ðM T  0:2ÞÞ1 , (11) pffiffiffi where M T ¼ 2q=c is the turbulent Mach number. Eq. (11) is a crude model that, however, has one important advantage in comparison with more modern models—it is computationally stable. For time-averaged Navier–Stokes equations, components of Eq. (4) have the following form: 3 2 3 2 0 r 7 6 7 6 6 0 7 6 ru 7 7 6 7 6 7 6 7 6 6 0 7 6 rv 7 7 6 7 6 7 6 7 6 ~ ¼ F~conv þ F ~ ¼ 6 0 7; F ~ ¼ 6 rw 7; W ~diff , U 7 6 7 6 7 6 7 6 6 0 7 6 rE 7 7 6 7 6 7 6 7 6 6 SðqÞ 7 6 rq 7 5 4 5 4 SðoÞ ro 3 2 Q 7 6 6 Qu þ psx 7 7 6 7 6 6 Qv þ psy 7 7 6 7 6 ~conv ¼ 6 Qw þ psz 7, F 7 6 7 6 6 QðE þ p=rÞ 7 7 6 7 6 7 6 Qq 5 4 Qo

2

~diff F

0

3

7 6 7 6 I xn 7 6 7 6 7 6 I xn 7 6 7 6 7. 6 I xn ¼6 7 7 6 6 ðI xn u þ I xn v þ I xn wÞ þ yn þ T kn 7 7 6 7 6 q 7 6 T n 5 4 o Tn

ð12Þ

conv

Vector F~ includes convective fluxes of mass, of three momentum components, of energy and of turbulence parameters. Convective fluxes are connected with the flow motion as continuous medium. First five of these fluxes are also present in Euler equation system for flows of inviscid gas. ~diff contains fluxes, connected with molecular Vector F diffusion (averaged effects of thermal motion of molecules) and with turbulent diffusion (averaged effects of fluctuation motion of flow volumes): (1) I in ¼ tin þ rRin ¼ ðtix þ rRix Þsx þ ðtiy þ rRiy Þ sy þ ðtiz þ rRiz Þsz —sum of molecular fluxes of momentum (viscous stresses tij ¼ mS ij , where S nij ¼ ðqui =qxj Þþ ~ ¼ ðqu=qxÞ þ ðqv=qyÞþ ðquj =qxi Þ  2ddij , d ¼ div V 3

ðqw=qzÞ) and of turbulent fluxes of momentum rRij . It is useful to represent the model for Reynolds stresses 00 u00 as follows: Rij ¼ ug i j 2 m Rij  q2  T Sij ¼ kaij , 3 r

(13)

where k ¼ q2 is the turbulence kinetic energy, aij ¼ ð2=3Þdij  M m S nij , M m ¼ mT =rk ¼ 0:09F wall F Dash ð1=oÞ. (2) yn ¼ qn þ rsn ¼ ðqx þ rsx Þsx þ ðqy þ rsy Þ sy þ ðqz þ rsz Þsz —sum of heat fluxes due to heat conductivity qm ¼ ðmC p =PrC p T Þð@T=@xm Þ and of turbulent 00 00 u  fluxes of thermal energy rsm ¼ rC p Tg

(3)

m Cp T ðmC p =PrT Þð@T=@xm Þ. Values of Prandtl numbers2: C T PrC p T ¼ 0:72 (for air), PrT p ¼ 0:9. T kn ¼ kn þ rK n ¼ ðkx þ rK x Þsx þ ðky þ rK y Þsy þ ðkz þ

rK z Þsz —sum of turbulent fluxes of turbulence kinetic energy- rK m ¼ ru00m ðu00 2 þ v00 2 þ w00 2 =2Þ and of molecular fluxes of k  km ¼ txm u00 þ tym v00 þ tzm w00 . (4) T qn ¼ T qx sx þ T qy sy þ T qz sz —summary diffusive flux of parameter q. o o o (5) T o n ¼ T x sx þ T y sy þ T z sz —summary diffusive flux of parameter o. 2 We shall use the term ‘‘Prandtl number for parameter f’’ for nondimensional value that is equal to the ratio of momentum diffusion coefficient (i.e. of the viscosity m) to diffusion coefficient of parameter f. As an example, for molecular heat transport f ¼ CpT, and classical Prandtl number we have PrC p T ¼ m=DC p T , where DC p T ¼ l=C p is coefficient of molecular diffusion of thermal energy (l is coefficient of heat conductivity). Therefore, diffusive flux of parameter f in direction of xk axis is represented as ðDf ðqf =qxk ÞÞ ¼ ðm=Prf Þðqf =qxk Þ.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

79

In the original version of the (qo)-model [47,48], the term T kn in the energy equation is omitted. But it may be calculated as follows. Fluxes of k should be derived from fluxes of Reynolds stresses in accordance with formula

The source terms for the variants [47,48] of the (qo)model has the following form: 8 1 ðrP  rÞ; < SðqÞ ¼ 2q   o : SðoÞ ¼ q2 ðC o11 þ C o12 F wall ÞrP  C o2 r :

k ¼ ðRxx þ Ryy þ Rzz Þ=2.

P Here P= i,jRij(qui/qxj)—production rate of turbulence kinetic energy, and e=q2o—its dissipation rate. For our realization of a numerical method, it is useful to rewrite these formulas as follows: (   SðqÞ ¼ r Aq o1 þ Bq þ C q o q;   (18) SðoÞ ¼ r Ao þ Bo o þ C o o2 ;

(14)

For fluxes of Reynolds stresses one may use the model by Donaldson and Rosenbaum [51]: !

 qRij qRim qRjm m Rij T T m  ru00i u00j u00m   m þ R þ þ , qxj qxi Pr ij qxm T

(15) where where the turbulent Prandtl number for Rij is equal to R PrT ij  0:81 [52]. Using (14), one may obtain from (15) the following model: ! Ryy Rxx Rzz T þ T þ T m m m   m þ TR T km  m 2 PrT ij

 qk qRxm qRym qRzm  þ þ þ . qxm qx qy qz Finally, using relation (13) and assuming that the structural coefficient aij varies slowly (that corresponds to assumption about equilibrium character of turbulence, which is also implied in Prandtl–Kolmogorov formula (10)), we obtain the model: !

mT qq qq qq qq k T m  2q m þ R þ axm þ aym þ azm . qxm qx qy qz Pr ij T

(16) For summary diffusive fluxes of parameters q and o, the following models are used:     mT qq mT qo q o Tm   m þ q ; Tm   m þ o . Prturb qxm Prturb qxm In the variant [47] of (qo)-model, the turbulent Prandtl numbers for parameters q and o have values Prqturb ¼ 1; Pro and in the variant [48]— turb ¼ 1:3, Prqturb ¼ Pro turb ¼ 2. Molecular viscosity of air is calculated using the Sutherland formula [53]:   T 3=2 273 þ 122 kg 5 , (17) m ¼ 1:72  10 273 T þ 122 ms where the temperature is measured in degrees Kelvin. In addition, equations for turbulent parameters contain source terms S(q) and S(o) that describe local effects of turbulence production (energy transmission from the mean motion to the largest turbulent eddies) and of turbulence dissipation (conversion of the energy of the smallest turbulent eddies into the energy of molecular thermal motion).

Aq ¼ 12C m F Dash F wall G;

Bq ¼ 13d;

C q ¼ 12,

Ao ¼ ðC o11 þ C o12 F wall ÞC m F Dash G, Bo ¼  23dðC o11 þ C o12 F wall Þ; C o ¼ C o2 . P Here G ¼ i;j S ij qui =qxj . In a variant [47] of the (qo)model the following set of constants is used: Co11 ¼ 0.045, Co12 ¼ 0.395, Co2 ¼ 0.92, and in another variant [48] Co11 ¼ 0.055, Co12 ¼ 0.5, Co2 ¼ 0.833. 2. Selection of the numerical scheme For the selection of a numerical method we used two key principles: 1. Methods ensuring high quality and reliability of results are more preferable than simplified fast numerical algorithms. As a consequence, one should not use for a construction of the numerical scheme only considerations of approximation, stability and minimization of computational efforts. The scheme should also take into account mathematical properties of the equation system to be solved. 2. In the transition from Euler equation system to Navier–Stokes equations and further to time-averaged Navier–Stokes equation system, the principle of inheritance should work: for an approximation of the same physical processes (e.g. convection) the same numerical methods should be used. The numerical method for the solution of Euler equations should be a ‘‘subset’’ of the numerical method for the solution of the Navier–Stokes equations, and so on. In accordance with the second principle, we use for an approximation of the convective terms in the timeaveraged Navier–Stokes Eq. (12) the same numerical scheme as in the method that was used in our previous paper [4] for the solution of the Euler equations. This is the explicit Godunov–Kolgan–Rodionov (GKR) scheme that nominally has a second order of approximation both in time and in space. It is a TVD method, and it uses an exact

ARTICLE IN PRESS 80

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

solution of Riemann problem. A brief description of the GKR scheme will be given in Section 3.1. For the following considerations it is important to point out that in our numerical method for the time-averaged Navier–Stokes equations, as well as in the GKR scheme, the time step is performed using a two-step procedure (predictor–corrector) that is equivalent to the two-step Runge–Kutta method (Euler method with central point). The predictor step is used for a preliminarily estimate of the ~n (with the first parameters at the time layer (n+1)-P ijk approximation order in time). This step is performed using scheme (4), where the values of fluxes and source terms are calculated at the time layer n. The Corrector step is used to obtain final values of parameters at the time layer (n+1)~n . Calculations are performed using the scheme (4) P ijk again, but fluxes and source terms are calculated at the time layer (n+1/2). Flow parameters in the cell centres at the time layer (n+1/2) are calculated using formula ~nþ1=2 ¼ ð1=2ÞðP ~n þ P ~n Þ. P i;j;k i;j;k i;j;k In the following section we shall present ideas, which were used in the selection of an approximation for diffusive terms (molecular and turbulent fluxes), and also for source terms that appear in the equations for turbulence parameters. 2.1. Selection of an approach to an approximation of diffusive terms The most popular schemes for a solution of the Navier–Stokes equations are implicit schemes. In particular, in widely distributed commercial codes (Fluent, ANSYS CFX, STAR CCM+, etc.) implicit schemes are used as the basic way for the solution of the Navier–Stokes equations, although some codes also allow the use of explicit methods. To justify the choice of an implicit scheme in the solution of Navier–Stokes equations, two main arguments are typically addressed: The first argument: Implicit schemes allow solving the problem of small time steps that arises in solution of tasks with strongly different time scales (e.g. in simulation of flows with boundary layers). But this problem may be solved also in the framework of explicit schemes. It will be shown in the next section of our paper. The second argument: Due to the presence of diffusive terms with second derivatives, Navier–Stokes equations describe the propagation of information with infinite speed. When an explicit scheme is used, during one time step the information passes a finite distance (not more than the scheme’s stencil). For example, let us consider the case, when the scheme’s stencil has five cells in each spatial direction (two cells in both directions from the current cell). If we introduce a perturbation in the cell j0, then after each time step the perturbations will pass two cells in both directions from the cell j0 along each index axis. However in implicit schemes all grid cells are connected with each other; as a result, the perturbation, which was introduced at the current time layer, produces the response in all cells

at the next time layer. Therefore, in implicit schemes the information propagates with infinite speed. Consequently, implicit schemes correspond better to the mathematical nature of Navier–Stokes equations. But a more detailed analysis of the propagation of information due to diffusion reveals that this argument is not completely satisfactory. Let us consider the model equation, which contains only diffusive terms: qu q2 u ¼D 2; qt qx

D ¼ const40.

(19)

This equation is parabolic, similar to all equations in the Navier–Stokes system, with the exception of the first equation (continuity equation). From a mathematical point of view, it describes the propagation of information in space with infinite speed. For example, remembering the classical task for Eq. (19), in which the delta-function is considered as the distribution of u(x) at the initial time moment: uðx; t ¼ 0Þ ¼ dðxÞ. The exact solution of this problem (see e.g. [54]) is

 1 x2 uðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi exp  . (20) 4Dt 4pDt Thus, at t40 the delta-function is converted into a Gaussian curve; with the growth of time, the amplitude of Gaussian curve diminishes, the width of its cupola increases; R þ1 at that, 1 uðx; tÞ dx remains constant (Fig. 12). This is the typical work of diffusion. Now we see that at t ¼ 0 u(x) ¼ 0 everywhere, except x ¼ 0; but at each indefinitely small t40 u(x)6¼0 everywhere. This fact just means that the speed of information propagation is infinite. But in nature infinite speeds are impossible, and not only because of relativistic restrictions, but also because of the

Fig. 12. Solution for a smoothing of the delta-function by diffusion.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

fact that molecular diffusion originates in chaotic motion of molecules, and consequently it is a slow process. The phenomenon of infinite speed is connected with the limited accuracy of assumptions that are used in derivation of the diffusion equation (19). The same statement may be made for Navier–Stokes equations: they include only a simplified model of the molecular transport (and may be derived through a simplification of the Boltzmann equations [55]). In [56], giving an example of the heat conductivity equation, it was shown that if we use more accurate assumptions in the derivation of the Navier– Stokes equations from Boltzmann equations, we may obtain equations with additional terms that describe the propagation of information with finite speed. Naturally, very small deviations of u(x) from a zero value at far distances from the point x ¼ 0 at t40 (see (20)) are physically unrealistic. They arise due to the inaccuracy of Eq. (19), and from a physical viewpoint we may neglect them. We may estimate the correct speed of information propagation, as the speed of widening the Gaussian curve cupola. It is easy to see from (20) that u(x) diminishes from the value umax (x ¼ 0) to the value (umax/N) at the distance pffiffiffiffiffiffiffiffiffi 2 KDt, where K ¼ ln N. Consequently, the speed of information propagation due to diffusion is equal to rffiffiffiffiffiffiffiffiffi d  pffiffiffiffiffiffiffiffiffi D 2 KDt ¼ K . V diff ðtÞ ¼ (21) dt t This is a finite speed. The only exception is at time moment t ¼ 0. But an infinite speed at t ¼ 0 is connected to the fact that we had given a delta function as the initial perturbation, and the delta function produces an infinite diffusive flux D(qu/qx). When the profile of u(x) becomes more smooth, the speed becomes finite; moreover, with the diminishing of the perturbation amplitude the speed of information propagation tends to zero! Naturally, this is merely estimation. The speed of information propagation due to diffusion cannot be determined exactly, because molecules are distributed over speeds in a random manner; for each K onep may find both ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi molecules with a speed higherpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi than KðD=tÞ and molecules with a speed less than KðD=tÞ. Formula (21) gives only characteristic values of the average speed of information propagation due to diffusion. The parameter K depends on our choice of the bounds for the Gaussian curve (20) cupola, i.e. of maximal value of u(x,t) that will be treated as non-physical and, hence, neglected. Therefore, from a physical viewpoint the information propagates in a viscous flow with finite speed. In this sense, we should make an inverse conclusion: explicit schemes are more preferable than implicit schemes. However, it is necessary to note that this conclusion may be made only about the physical speed of information propagation. From a mathematical viewpoint, nothing has changed: we approximate the Navier–Stokes equations that describe the propagation of information with infinite speed. Because explicit schemes give finite speed of

81

information propagation, numerical solutions of the Navier–Stokes equations with the use of an explicit scheme will inevitably be different from an exact solution of these equations. In the following section we shall show that this difference is small, and the quality of numerical solutions by an explicit scheme is quite acceptable and even better than the quality of solutions by implicit schemes. This difference may also be considered as small perturbations that are introduced at each time step. An explicit numerical scheme should be stable to such perturbations. In the Section 2.3 we shall show that the stability condition for explicit scheme enables the determination of the concrete value of the parameter K. This value represents the level of perturbations, for which the explicit scheme is still stable. Now we are going to compare quantitatively the errors in the description of diffusion, which are produced by explicit and implicit schemes. We shall use the model Eq. (19) again. This equation may also be rewritten in standard conservative form: qu qF þ ¼ 0, qt qx where F ¼ D(qu/qx) is the diffusive flux of the parameter u. An approximation of this flux, analogous to (4), may be represented as follows: unþ1 ¼ unj  j

 t F jþ1=2  F j1=2 . h

(22)

The numeration, which is used in (22), is shown in Fig. 13. In a 1D case the cell j is a segment of the grid line t ¼ const between the nodes (j1/2) and (j+1/2); unj and R xjþ1=2 Fj+1/2 are approximations of the values ð1=hÞ xj1=2 Rt uðx; tn Þ dx and ð1=tÞ tnnþ1 uðxjþ1=2 ; tÞ dt, accordingly.

Fig. 13. Grid and map of indices for model task.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

82

Scheme (22) for an approximation of Eq. (19) should satisfy at least the following requirements: 1. It should have at least a second-order approximation. First-order schemes are not good for approximation of (19). Indeed, errors produced by schemes of first order in time (see e.g. below schemes (32) and (33)) have a character of numerical viscosity that would compete with physical viscosity (and too detailed grids would be necessary to exclude this competition). 2. Because Eq. (19) is not modified after a change of x by (x), the information propagates due to diffusion symmetrically in all directions. Consequently, the stencil of the flux Fj+1/2 approximation should be symmetric with respect to the cell border (j+1/2). In the case of a uniform grid, these requirements are satisfied for a central-difference approximation F jþ1=2 ¼ D

Djþ1=2 u ujþ1  uj ¼ D . Djþ1=2 x h

For this approximation the scheme (22) may be rewritten in the form that more resembles the approximation of the differential equation (19): ujnþ1  unj uj1  2uj þ ujþ1 ¼D . t h2

Let us assume we are going to solve the following Cauchy problem for the diffusion equation (see Fig. 13): 8 qu q2 u > > qt ¼ D qx2 ; > > < uðx; 0Þ ¼ jðxÞ; x 2 ½0; L; > uð0; tÞ ¼ c1 ðtÞ; t40; > > > : uðL; tÞ ¼ c ðtÞ; t40: 2

Consider the time layer t ¼ const and expand the exact solution of this Cauchy problem in Fourier series:

Substituting this expansion into Eq. (19), we find the dependence of the amplitude factor Uk(t) on time:

uj  t

¼D

unj1



2unj 2

þ

unjþ1

h

nþ1=2

2 2 2 U k ð0Þeð4p =L Þk Dt eið2p=LÞkx .

(25)

k¼1

Now consider the numerical solution of the same Cauchy problem on a grid containing N cells of equal size h at the segment xA[0; L] (Fig. 13). We may expand the numerical solution in a finite Fourier series: unj ¼

þðN1Þ=2 X

U nk eið2p=LÞkjh .

(26)

(For brevity, we shall consider only the case, when N is odd.) It is worth noting that as a result of the task discretization, only harmonics with wavelength lXlmin ¼ L/((N1)/2)E2h can be are considered. Shorter wavelengths are not resolved by the grid. Substituting expansion (26) in the scheme equation (24), we find the dependence of the amplitude factor U nk on time. Finally, we represent the numerical solution in the following form: unj ¼

(23)

.

ðN1Þ=2 X

U k ð0Þrn ðjk ; C diff Þeijk j ;

k¼ðN1Þ=2

jk ¼

2pkh 2pk ¼ . L N (27)

2. Corrector step: uj1 unþ1  unj j ¼D t

þ1 X

uðx; tÞ ¼

1. Predictor step: unj

U k ðtÞeið2p=LÞkx .

k¼1

k¼ðN1Þ=2

If on the right-hand side only the parameters from the known time layer are used, then we obtain an explicit scheme; if any parameters from the unknown time layer (n+1) are also used, then we obtain an implicit scheme. Initially, let us generate an explicit scheme. To obtain a second-order approximation in time, we may use the ‘‘predictor–corrector’’ procedure (as in the GKR scheme for convective fluxes):

n

þ1 X

uðx; tÞ ¼

nþ1=2

 2uj

2

h

nþ1=2

þ ujþ1

;

nþ1=2

uj

¼

unj þ unj . 2

Introducing the diffusive analogue of the Courant number, a non-dimensional parameter Cdiff ¼ D(t/h2), and using (23), one finally obtain

Because (N1)/2pkp(N1)/2, jk ¼ (2pk/N)A(p;p). For scheme (24)   rðjk ; C diff Þ ¼ 1 þ 4C diff sin2 ðjk =2Þ 2C diff sin2 ðjk =2Þ  1 ¼ 1  2Z þ 2Z2 ;

Z ¼ 2C diff sin2 ðjk =2Þ.

ð28Þ

But the exact solution (25) may be represented in each point (xj, tn) in the form (analogous to (27)) uðxj ; tn Þ ¼

þ1 X

U k ð0Þrnex ðjk ; C diff Þeijk j ,

k¼1

ujnþ1 ¼ ð1  2C diff þ 3C 2diff Þunj þ ðC diff  2C 2diff Þðunj1 þ unjþ1 Þ þ 12C 2diff ðunj2 þ unjþ2 Þ.

ð24Þ

For an analysis of the properties of this scheme, we shall use the Neumann’ spectral method.

2

rex ðjk ; C diff Þ ¼ eC diff jk . Then let us compare expansions of exact and numerical solutions. Harmonics of the exact solution are positive and damp with the growth of time step number n: 0orexp1.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Harmonics of the numerical solution are also always positive, but at Cdiff40.5 the values r41 arise, and the amplitude of harmonics grows with no limit. Therefore, our scheme is stable at Cdiffp0.5, i.e. at tp(h2/2D). Now we consider an implicit scheme for Eq. (19), which uses the same method as for the approximation of diffusive fluxes and also has a second-order approximation (Crank– Nicolson scheme [57]): ujnþ1  unj t

" # nþ1 nþ1  2unþ1 þ ujþ1 unj1  2unj þ unjþ1 uj1 1 j ¼ D þ . 2 h2 h2 (29)

Analysis of this scheme by Neumann method gives rðjk ; C diff Þ ¼

1  2C diff sin2 ðjk =2Þ . 1 þ 2C diff sin2 ðjk =2Þ

(30)

It is easy to see that for each Cdiff we have |r(jk, Cdiff)|p1; consequently, this scheme is absolutely stable. Therefore, it is possible to perform calculations with a time step tbh2/2D and, accordingly, to move along the physical time axis much more quickly. But aside from this excellent property, implicit schemes have also essential drawbacks. First of all, let us compare the computational cost per time step. In an explicit scheme the value of unþ1 may be directly (explicitly) expressed from j the finite-difference equation (24). In an implicit scheme (see (29)) finite-difference equations for all cells are related with each other, and therefore it is necessary to solve a system of N equations (for all cells) to find unþ1 . As a result, j the computational cost per time step and the required volume of RAM become much larger than in the case of an explicit approach. Moreover, computational cost grows very quickly together with the quantity of spatial coordinates (the stencil of the scheme and the total quantity of cells increase quickly with the growth of the task dimension). For full Navier–Stokes equations this problem is aggravated, because the system of implicit equations for ujnþ1 becomes nonlinear and may be solved only by an iterative method. Implicit schemes may be more effective than explicit schemes only in the case when a calculation is performed with very large time steps, much larger than the stability limit of explicit schemes. The higher the dimension of task, the larger time steps are necessary to obtain an advantage in comparison with explicit schemes. But a growth of the time step leads to an essential decrease of the quality of the implicit scheme. The quality of explicit and implicit schemes may be compared on the basis of a Neumann spectral analysis. During some time T, in the exact solution the initial modulus of amplitude of harmonics is multiplied by |rex(jk, Cdiff)|T/t, in the numerical solution—by |r(jk, Cdiff)|T/t. Therefore, the relative amplitude error for the given

83

harmonic will be equal to   T =t    T=t     U k ð0Þ rex ðjk ; C diff Þ  U k ð0Þ rðjk ; C diff Þ   Eðjk ; C diff Þ ¼   T=t   U k ð0Þ rex ðjk ; C diff Þ 

 C =C ¼ rðjk ; C diff Þ=rex ðjk ; C diff Þ T diff  1,

where CT ¼ D(T/h2). The value of the parameter CT may be chosen arbitrarily. We shall take CT ¼ 0.5. It is important that due to such a definition we may compare the errors that have accumulated during one and the same physical time T, and the result does not depend on the value of time step in calculations with different schemes. Mean amplitude errors over some range of wavelengths lA[l0;+N) (i.e. jA[0; 2ph/l0])) may be defined as follows: Z 2ph=l0 ¯ diff Þ ¼ l0 ðEðj; C diff Þ  Eðj; 0ÞÞ dj. (31) EðC 2ph 0 Mean R p errors over the whole range of wavelengths ð1=pÞ 0 Eðj; C diff Þ dj principally cannot be equal to zero at all values of Cdiff, because the errors inevitably grow, when the wavelength of harmonics l approaches the limit of the scheme resolution lmin ¼ 2h (j-p). That is why for a comparison of errors of different schemes we subtract the value E(j,0) from the integrand in (31). In the case of diffusion the equation (contrary to tasks with convection) schemes usually give the best prediction of the damping rate of harmonics, if the non-dimensional time step Cdiff0. Indeed, let us consider a decrease in time step t at fixed grid step in space h. Then Cdiff will also diminish. Using n ¼ t/t in (27), we obtain for the longest wavelengths (|jk|51 and lk/L ¼ 1/|k|_b1), which are the best of all the resolved by the grid; harmonics of the numerical solution (27) becomes more and more close to harmonics of the exact solution (25) as Cdiff tends to zero. ¯ In Fig. 14 the curves EðCÞ are plotted for the secondorder approximation explicit scheme (24) (solid line), for the second-order implicit scheme (29) (dotted line) and also for two first-order schemes: (1) For the explicit scheme of the first-order approximation in time unþ1  unj j t

¼D

unj1  2unj þ unjþ1 h2

.

(32)

(2) For the implicit scheme of first-order approximation in time unþ1  unj j t

¼D

nþ1 unþ1 þ unþ1 j1  2uj jþ1

h2

.

(33)

In Figs. 14(a) and (b), data for CdiffA[0; 0.25] are presented. This is the range of Cdiff, where explicit schemes give good results. Near the stability limit Cdiff ¼ 0.5 the errors of explicit

ARTICLE IN PRESS 84

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 14. Behaviour of the error function for different schemes: (a) CdiffA[0; 0.25], wavelengths lX2 h; (b) CdiffA[0;0.25], wavelengths lX6 h; (c) CdiffA[0;10], wavelengths lX2 h; (d) CdiffA[0;10], wavelengths lX6 h.

schemes increase sharply; therefore, we may recommend using explicit schemes (24) and (32) with Cdiffp0.25. In Figs. 14(c and d) results for CdiffA[0; 10] are shown. Figs. 14(a and c), are obtained for the whole range of wavelengths (jA[0; p], lX2h), and Figs. 14(b and d), for the range jA[0; p/3], i.e. for the wavelengths lX6h. (For jA[0; p] the mean error is mostly based on errors in the representation of the shortest wavelengths, which are badly resolved by any scheme.) In Table 1 the maximum values of the mean error ¯ diff Þ are presented (for the explicit schemes, in the range EðC of their applicability, CdiffA[0; 0.25], and for implicit schemes, at Cdiff ¼ 10).

An analysis of Fig. 14 and Table 1 leads to following conclusions: 1. At small values of the non-dimensional time step Cdiff the errors of schemes with first-order approximation in time are higher than the errors of schemes with second order in time. But with the growth of Cdiff the errors of second-order schemes increase faster and become higher than the errors of first-order schemes. 2. The maximum error of an explicit scheme (24) in the range of its applicability is much less than the errors of an implicit scheme (29) in its working range (at Cdiffb1). For wavelengths lX2h the maximal amplitude error of

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

85

Table 1 Maximal values of mean relative error for different schemes (CT=0.5)

First order of approximation Second order of approximation

Explicit scheme, lX2h

Implicit scheme, lX2h

Explicit scheme, lX6h

Implicit scheme, lX6h

2.25 1.58

10.7 13.3

0.0151 0.0009

0.129 0.150

Table 2 Maximal values of mean relative error for different schemes (CT=10)

First order of approximation Second order of approximation

Explicit scheme, lX2h

Implicit scheme, lX2h

Explicit scheme, lX6h

Implicit scheme, lX6h

1.91  1023 4.4  1028

1.18  1039 4.38  1040

0.352 0.036

278.3 1777

the implicit scheme (29) is 8.5 times higher than maximum error of the explicit scheme (24), and for wavelengths lX6h it is 167 times higher. If we take a longer physical time interval T (we consider an accumulation of errors during the interval T), then the difference between explicit and implicit schemes becomes larger. Let us take, e.g. CT ¼ 10. Analogous results for this physical time are given in Table 2. Conclusions 1 and 2 are confirmed. But in this case for wavelengths lX2h the maximum amplitude error of the implicit scheme (29) is 1011 times higher than maximal error of the explicit scheme (24), and for wavelengths lX6h, it is 49,400 times higher! The relative errors for l 2h appear to be excessively high. Therefore, at CT ¼ 10 all schemes are inapplicable for l 2h. For harmonics with wavelength lX6h only explicit schemes give errors less than 1 (i.e. less than exact value of amplitude). For CT ¼ 10 the errors that are less than the exact value of the amplitude are achieved by the implicit scheme of the second-order approximation only at l\10h, and by the first-order implicit scheme only at l\11h. Therefore, the quality of an implicit scheme within its working range is much worse than the quality of an explicit scheme in its working range. In reality this problem is aggravated, because an implicit scheme is usually simplified to improve its efficiency. There are various ways to simplify schemes: by a linearization of implicit difference equations, by splitting of the difference operator or other approaches. But all these simplifications will reduce the scheme’s quality. It is necessary to note that all these arguments consider explicit and implicit schemes for the non-stationary diffusion Eq. (19). If we want to obtain stationary solutions by means of a stabilization method, than in the stationary limit the dependence of the numerical solution on time vanishes. In this situation, all four considered schemes (24), (29), (32) and (33) coincide. Thus, if the stationary solution is unique, it should be the same for all four schemes. Naturally, we cannot be sure in the uniqueness of a stationary solution. Existence and uniqueness theorem are

not proven for general cases even for the Navier–Stokes system of differential equations [58]. All the more, it is not proved for a finite-difference approximation of this system. But even if we would know exactly that a stationary solution is unique, we should remember that this stationary solution is obtained as the limit of a non-stationary (or pseudo- non-stationary) process. And if we use an implicit scheme, this non-stationary process will be described with a very poor quality. As a consequence, it is possible that the convergence of the numerical solution to the stationary state will become slow. At worst, we may have no convergence or we may obtain a convergence to a nonphysical stationary state (if the stationary state is still not unique). The risk of such an unsuccessful outcome increases, when a flow with a complex physical structure is simulated. That is why in the use of implicit schemes the authors often start calculations with small values of the Courant number (in spite of the fact that an implicit scheme is not efficient in this case), and only when approaching the stationary state the Courant number (and, accordingly, the speed of motion along the physical time axis) increases. Therefore, if it is necessary to obtain a high quality and reliability of a numerical solution, the use of explicit schemes is preferable. This conclusion is beyond doubt if we are talking about descriptions of non-stationary flows. If we look for a stationary solution, then the use of an implicit scheme is at least a risky procedure, which does not guarantee that a physically correct stationary solution will be obtained. That is why for our numerical method we have chosen an explicit approximation of the Navier– Stokes equations. 2.2. Methods for solving the small time step problem Despite of all advantages of explicit schemes described in the previous section, the most popular approach for a solution of the Navier–Stokes equations is the use of implicit schemes. The reason is that in the framework of implicit schemes there is no problem of small time steps.

ARTICLE IN PRESS 86

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

In the previous section we have demonstrated (see also Fig. 14(c)) that an explicit scheme is stable and has acceptable quality under the following condition: h2 . (34) 4D Let us compare this time step restriction with the wellknown Courant–Friedrichs–Levi condition that arises in the approximation of convective fluxes. A good model equation for an analysis of methods for the simulation of convection is the advection equation

tp

qu qu þ a ¼ 0; a ¼ const40. qt qx For this equation, the CFL restriction has the following form: h (35) tp . a Not infrequently the opinion occurs that condition (34) is stronger than (35), because the time step is proportional to h2. Nevertheless, let us compare restrictions (34) and (35). We see that the convective restriction (35) is stronger than (34), if ah 44. (36) D Our experience shows that condition (36) is satisfied in the majority of cells, except for the case of very small cells (e.g. within the viscous sublayer of a turbulent boundary layer) and the case of very high turbulent viscosity. Nevertheless, the problem of small time steps in the calculations of viscous gas flows with explicit schemes remains. For flows with boundary layers we should use grids with very small cells, for which even the CFLrestriction (35) becomes too strong. In the solution of a stationary task by the stabilization method, the stabilization of the flow is achieved, when information from the boundaries has passed the computational domain several times. If the characteristic linear size of the computational domain is equal to L, the characteristic value of the flow velocity inside the domain is V, and the characteristic speed of sound is c, then the flow will stabilize (i.e. will come to a steady state) after the physical time T 10(L/V+c). Initially, let us consider the calculation on a quasiuniform grid, containing about 100 cells in each index direction (the size of the grid step is h L/100). If we use an explicit scheme, then the calculation is performed with a time step that is determined by the CFL restriction, and it may be estimated as t (h/V+c). For this case, we may estimate the number of time steps until flow stabilization by N ¼ (T/t) 1000. That is why in calculations on a quasiuniform grid about several thousands of time steps are typically sufficient to achieve a steady solution.3

Reh ¼

3 This estimation is correct, if the flow velocity is comparable to the speed of sound. But if we consider the case M ¼ V/c-0, then the duration

Now let us consider the calculation of a flow with boundary layers. In this case we should compress the grid strongly towards the solid surface—with the objective to describe in detail the no-slip region, which determines friction and heat fluxes at the wall. Usually it is recommended to use Rehmin 1, where hmin is the size of the near-wall cells of the grid. As for the largest cells (in the inviscid core of the flow), their size may be estimated as hmax L/100, as before. Then the ratio of cell sizes (hmax/ hmin) (ReL/100). As a result, in the case of a fully developed turbulent boundary layer, at ReL 107, the value of hmin may be 105 times smaller than the size of largest cell hmax. Accordingly, the time step t (hmin/V+c) becomes 105 times smaller than in a calculation on a quasi-uniform grid, and the quantity of time steps until flow stabilization N T/t increases 105 times, i.e. now it is counted at hundreds of millions! But the corresponding time of calculation would be unacceptable large. However, an implicit scheme allows in principle a calculation with arbitrary large time steps and, consequently, to get a result immeasurably faster. It is this factor that makes implicit schemes so popular. But, as it was shown in the previous section, the penalty for speed of getting the results is loss of quality in the description of non-stationary processes and diminishing the reliability in the solution of stationary tasks. Fortunately, there are ways to solve the problem of small time steps within the explicit approach to approximation of the equations. Our numerical scheme includes two such methods. Let us consider for example a 1D non-stationary task for the model advection equation qu qu þ aðxÞ ¼ 0. qt qx Fig. 15(a) shows the distribution of local time step restrictions tmax , which are determined by imposing the j CFL-condition (35). Because the convective term in the advection equation describes the transmission of information with the speed dx/dt ¼ a(x), the information passes during one time step the distance CFLj hj, where CFLj ¼ aj(tj/hj) is the local Courant number. The standard approach to the organization of time stepping (global time stepping) may be described as follows: we choose the strongest time step restriction-tmin ¼ minj tmax and make in all cells the same j time step tmin. If tmin 5tmax ¼ max tmax (this is characterj istic for strongly non-uniform grids, which are used in the simulation of flows with boundary layers), then in the

(footnote continued) of flow stabilization is determined by time, which is necessary for flow to pass the computational domain, i.e. T L/V. As a result at M-0 we obtain N (100/M)-N. That is why explicit schemes are inapplicable in simulation of flows that are close to the limit of an incompressible liquid.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

87

Fig. 15. To explain global, local and fractional time stepping.

majority of cells the Courant number is

modified equation

tmin tmin CF Lj ¼ aj ¼ max 51. hj tj

CFLðxÞ

Consequently, in the majority of cells the information passes during one time step only a small part of the grid step hj. Thus, information propagates very slowly over the computational domain, and we need a very large quantity of time steps, until the flow will adjust to the boundary conditions. This problem is called the problem of small time steps. The first method of overcoming this problem in the framework of explicit scheme is applied, if the stationary flow is calculated by the stabilization method. In this case the quality of the description of transitional non-stationary process is not of interest; it is only important to obtain convergence of this process to a stationary state. In this situation we use the method of local time stepping well known since a long time (see e.g. [5]). In this approach, the calculation in each cell is performed with a time step that is determined by a local time step restriction for this separate cell (Fig. 15(a)). As a result, the value of the time step varies from one cell to the other. But when all parameters unþ1 are determined, these j parameters are attributed formally to the same time layer. After that, this procedure is repeated, until a stationary solution is reached. When such a procedure is used, convergence to stationary state accelerates dramatically, because in each cell we use the maximum possible value of the time step, CFLj 1, and information passes in one time step the maximum possible distance ( hj). As a result, the flow adjusts more quickly to stationary boundary conditions. It is easy to understand that the number of time steps should diminish O(tmax/tmin) times in comparison with global time stepping. Naturally, using such procedure leads to an incorrect description of the non-stationary process as a whole. Obviously, local time stepping is equivalent to an increase of the local speed of information propagation (dx/dt) ¼ a(x)(CFLj)1 times, i.e. to the solution of

with global time stepping. Because at q/qt ¼ 0 this equation coincides with the initial advection equation, the stationary solution should be the same, as in a solution of the initial equation with global time stepping—if the stationary solution is unique. As in local time stepping the non-stationary process as a whole is described non-physically, we meet the same problems as in the case of implicit schemes: we cannot guarantee the convergence and the uniqueness of the stationary solution. Nevertheless, local time stepping has advantages in comparison with very-high-Courant-number global time stepping using implicit schemes. As it was shown above, at very high values of time step the quality of implicit schemes becomes too poor. Therefore, in a calculation with the use of an implicit scheme the nonstationary process is described non-physically in each cell. However, in the case of local time stepping the nonstationary process is described physically in each separate cell. Indeed, in such an approach each cell may be considered as separate computational domain with nonstationary boundary conditions. The quality of the numerical scheme is kept in such case. In a simulation of tasks with turbulent boundary layers, which were considered above, local time stepping enables a diminishing number of time steps 100–1000 times in comparison with global time stepping and to obtain stationary solution within an appropriate CPU time. The second method of overcoming the problem of small time steps in the framework of an explicit scheme is applied in case of a calculation of an unsteady flow combined with a correct description of each step of its development. In such cases it is possible to use the method of fractional time stepping. The idea of this approach is obvious enough, but we have only a very small amount of information about its use by various authors. We may only point to paper [59], where this approach is used in conjunction with another

qu qu þ aðxÞ ¼ 0 qt qx

ARTICLE IN PRESS 88

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

numerical method for a simulation of non-stationary flows with finite-rate combustion. In fractional time stepping, the calculation is performed with different time steps in various cells. In each cell, a time step satisfies the local time step restriction. But the quantities of local time steps in different cells are different and are adjusted in such a way that periodically all cells reach the same level of physical time. Due to this fact, it is possible to describe the non-stationary process correctly. The moment, when all cells have reached the same level of physical time, will be called the ‘‘finishing of global time step’’. For example, if in the cell A the local time step is equal to tmax, and also in the cell B, tmax/2, and in the cell C, tmax/4, then during one global time step it is necessary to make one local time step in the cell A, in the cell B, two local time steps, and in the cell C, four local time steps (see Fig. 15(b)). Thus, in each cell we divide the global time step into smaller fractions, local time steps, in such a way that local time steps satisfy to local time step restriction for the current cell. This is origin of the term—fractional time stepping. Let us consider in detail the organization of fractional time stepping. Before the beginning of a global time step, we determine in each cell of the computational domain the maximum possible value of time step tstab ijk , which is determined by the local stability condition. After that, the maximum time step in the whole computational domain is determined— tmax ¼ maxi;j;k tstab ijk . Then in each cell we determine the value of the local time step. In the cell with number (i, j, k) we take the local value of time step to be equal to tijk ¼ tmax/2I. The integer parameter l is chosen to satisfy stab the local stability condition: tstab ijk =2otijk ptijk . As a result, all sets of cells of the computational domain are subdivided into several subsets, which will be called levels. In all cells of level m the value of local time step is equal to tijk ¼ tmax/2Mm, where M is the quantity of levels. In the cells of the 1st level the local time step is minimal, in the cells of the level M we have tijk ¼ tmax. Each local time step in our numerical method (which is based on the GKR scheme) consists of two inner steps— predictor and corrector. If at the beginning of a time step parameters in the cell correspond to physical time t0 and the value of local time step in this cell is equal to t, then during this local time step 2 inner steps will be performed: after the first inner time step (predictor) we shall obtain parameters at the time layer (t0+t/2), after the second inner time step (corrector), at the time layer (t0+t). In the organization of the fractional time stepping procedure it is very important to synchronize time steps in different cells. Let us consider for example the case when we have three groups (levels) of cells. At the first level of cells the value of local time step is tmax/4, at the second tmax/2, at the third tmax (Fig. 15(b)). The algorithm of calculations during one global time step is represented in Table 3.

Table 3 Algorithm of calculations during one global time step Number of passage

Time moment

1st level t=tmax/4

2nd level t=tmax/2

3rd level t=tmax

1 2 3 4 5 6 7 8

t0 t0=tmax/8 t0=tmax/4 t0=3tmax/8 t0=tmax/2 t0=5tmax/8 t0=3tmax/4 t0=7tmax/8

1st inner step 2nd inner step 1st inner step 2nd inner step 1st inner step 2nd inner step 1st inner step 2nd inner step

1st inner step

1st inner step

2nd inner step 1st inner step

2nd inner step

2nd inner step

If total quantity of cell levels is equal to M, then during one global time step we perform 2  2M1 ¼ 2M passages over the whole computational domain. During the passage, the structure of calculations in each cell depends upon the number of the level m, to which the current cell belongs. If r is number of passage and (r1) is divided by 2m without a remainder, then in this cell the predictor step is performed; if (r12m1) is divided by 2m without a remainder, then in this cell the corrector step is performed. In addition, in the adjacent cells we should know the flow parameters at the same time as in the given cell (e.g. for the sixth passage at the time moment t0+5tmax/8, see Table 3). These parameters are determined by interpolation. In each cell we store the values of flow parameters at two time layers: at the initial time layer, which corresponds to beginning of this local time step in this cell; and at the time layer, which was reached in the previous inner time step. For example, during the sixth passage we make neither predictor nor corrector steps in the cell of the second level; but in this cell we store values of the parameters at time layer t0+tmax/2 (beginning of the current local time step) and at the time layer t0+3tmax/4 (this is the time layer, which was reached after a predictor step in this cell). Interpolation between these two layers allows us to obtain parameters at the layer t0+5tmax/8 that is of interest for us at the current passage. The obvious disadvantage of the described algorithm of the fractional time stepping organization is that it is not conservative: fluxes through the common border of the cells (i, j, k) and (i+1, j, k) are calculated separately for each of these cells using independent formulas and in different time moments. A conservative algorithm of fractional time stepping is possible, but it requires much more calculations and much more RAM. But in tests, which were performed for the model equation, transition to conservative algorithm resulted in a very small variation of the numerical solution. So, we have decided to use a non-conservative version of the fractional time stepping. In the case of rational programming, the methodology of fractional time stepping diminishes the total time of computation in comparison with global time stepping (because in large cells we make only a small number of time steps). In addition, convection is described in the case of

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

fractional time stepping with higher quality than in the case of global time stepping. Indeed, with the use of the differential approximation method [60] one may show that approximation errors in the description of convection are proportional to (1CFL). If a calculation is performed at strongly non-uniform grids, then for global time stepping in the majority of cells we have CFL51 (as it was shown above), and the approximation errors are maximal. But for fractional time stepping in all cells we have CFLA[0.5; 1]. Naturally, in the calculation of stationary flows by the stabilization method the fractional time stepping decelerates convergence to a stationary state in comparison with local time stepping. Especially this is true for flows with boundary layers, where the values of maximal and minimal time steps differ by several orders of magnitude. In such tasks calculation with fractional time stepping will be too prolonged. This is an inevitable penalty for a correct description of non-stationary processes. If we are interested in a stationary state of a flow with boundary layers, we may choose local time stepping. But if we are interested in nonstationary processes in flows with boundary layers, local time stepping is inapplicable, and fractional time stepping does not solve the problem of small time steps. A possible solution of this problem is the use of the zonal approach, in which near-wall regions of flow (viscous sublayers and buffer zones of boundary layers, where grid cells should be too small) are calculated by an implicit scheme, and the remaining part of the flow—by an explicit scheme with fractional time stepping. We may assume that information propagates very quickly across the thin layer, where the implicit scheme is used. Then the flow in this thin layer in the most part of the computational domain (possibly except for the regions of flow separation and reattachment) should adjust quickly to non-stationary processes in the remaining part of the flow. Then we may consider the processes in this thin layer as quasi-stationary, and application of an implicit scheme in such layer will be quite justified, because for quasi-stationary processes quality of implicit schemes approaches the quality of explicit schemes. 2.3. Generalization of the Courant–Friedrichs–Levi condition to tasks with interaction of convection and diffusion Now we shall consider problems that arise in the solution of equations, containing both convective and diffusive terms. Let us take a scalar model equation with convective and diffusive terms: qu qu q2 u þ a ¼ D 2 ; a; D ¼ const40. (37) qt qx qx For the approximation of the convective term a(qu/qx) we may use the GKR scheme. The stability condition for this scheme is given by Eq. (35); this is the classical Courant–Friedrichs–Levi condition (CFL). For the diffusive term D(q2u/qx2) we shall use approximations (23) and

89

(24). For stability and acceptable quality of this scheme we have obtained the time step restriction (34). As for the combined numerical method for Eq. (37), it seems natural to obtain the stability condition by taking the stability conditions for each of the above two schemes and to choose the stronger of them. This leads to the following time step restriction:   h h2 . (38) tpmin ; a 4D Sometimes this condition provides stability of the scheme (e.g. in the simulation of mixing layers). But in more stiff problems (e.g. flows with boundary layers) the scheme becomes unstable for some reason. As the simplest example of such an instability, we may consider an attempt to solve numerically the following boundary-value problem for Eq. (37): 8 qu qu q2 u > a; D ¼ 1; > < qt þ a qx ¼ D qx2 ; (39) t ¼ 0; xe½4:5; 5:5 : uðxÞ ¼ 0; > > : t ¼ 0; x 2 ½4:5; 5:5 : uðxÞ ¼ 1: We shall solve this task at the interval xA[0;10] on a uniform grid of 100 cells. At boundaries x ¼ 0 and 10, let us impose condition u(t) ¼ 0. In the exact solution of this problem we should see, how the initial rectangular perturbation is gradually propagating downstream with the velocity a due to convection and simultaneously smoothed (spread in space) due to diffusion. In Fig. 16(a) we present distributions u(x), obtained at sequential stages of this task calculation by our scheme with stability condition (38). One may see that approximately after 150 time steps the rapid destruction of the numerical solution arises. It is important that oscillations occur along the whole length of the computational domain. Consequently, instability is produced not by perturbations from the boundaries, but by an inner development of the numerical solution. Therefore, condition (38) does not guarantee the stability of our scheme. So, we have undertaken investigations to improve this condition. Unfortunately, it is impossible to obtain exact time step restrictions using a Neumann spectral analysis, because our approximation of the convective term is nonlinear and depends on the local behaviour of the numerical solution. But even if we would succeed in the Neumann analysis, we would obtain some abstract time step restriction for the approximation of the model equation (37); correct generalization of such restrictions to the case of Navier–Stokes equations will be, most probably, impossible. It is desirable to obtain a reliable stability condition that would have clear physical sense and due to this fact would be generalized easily to the case of full Navier–Stokes equations. Because the convective term a(qu/qx) describes the transport of information in space with the speed dx/dt ¼ a, the CFL-condition (35) may be interpreted as follows: the

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

90

Fig. 16. Two numerical solutions of task (39): (a) with stability condition (38) and (b) with condition (42).

distance of information propagation in one time step due to convection should not be more than one grid step in space. But in Section 2.1 we have derived the estimation for finite speed of information propagation due to diffusion— (21). Then for the scheme (23) and (24) it would be possible to formulate a stability condition, analogous to CFLcondition: the distance of information propagation in one time step due to diffusion should not be more than one grid step in space. Because in one time step the ffi perturbation pffiffiffiffiffiffiffiffiffi passes due to diffusion the distance 2 KDt (see Section 2.1), then we obtain the following condition: pffiffiffiffiffiffiffiffiffiffi h2 2 KDtph or tp . 4KD

(40)

But for K ¼ 1 this condition coincides with condition (34)! Thus, we have managed to find a physical sense of condition (34); moreover, we have found the value of K factor.

Naturally, for different schemes we shall obtain different values of K. A coincidence of formulas (40) and (34) at K=1 means that in the model of this numerical scheme the cupola of Gaussian curve is bounded by the points, where uðxÞ ¼

umax umax umax ¼ ¼ . N expðKÞ e

If the amplitude of perturbation is more than e times smaller than the maximal amplitude in the perturbed region, then this perturbation has no influence on the stability of the given scheme. Now it is easy to formulate the stability condition for our scheme for combined equation (37). The total speed of the propagation of information in the flow direction is equal to sum of convective speed a and speed of information pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi propagation due to diffusion (21), i.e. to a þ KðD=tÞ. Thus the distance of information pffiffiffiffiffiffiffiffiffiffipropagation in one time step will be equal to at þ 2 KDt. According to the CFL principle, this distance should not be larger than the grid

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

step in space: pffiffiffiffiffiffiffiffiffiffi at þ 2 KDtph.

(41)

The solution of inequality (41) with respect to t at K ¼ 1 gives us the following time step restriction for our scheme: tptmax

4tdiff ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 , 1 þ 1 þ 4tdiff =tconv

(42)

where tconv ¼ (h/a) is a ‘‘convective’’ time step restriction (35), and tdiff ¼ (h2/4D) is a ‘‘diffusive’’ time step restriction (34). If tconv5tdiff, then (42) is converted into condition (35), and if tconvbtdiff, into condition (34). But in the flow regions, where tconv tdiff, we obtain more complex restriction, which takes into account the interaction of convection with diffusion. One may note that condition (41) at pKffiffiffiffiffiffiffiffiffi ¼ 1ffi may also be rewritten in the form CFL þ 2 C diff p1. From this formulation one may easily see that condition (41) appears to be stronger than each of individual conditions (34) and (35). For example, in the case, when tconv=tdiff, we shall obtain from (42) tmaxE(tdiff/2.6), that in this case gives CFLmaxE0.38 and (CFLmax)diffE0.096. For comparison: condition (35) always gives CFLmax ¼ 1, and condition (34), (Cdiff)max ¼ 0.25. Therefore, condition (42) is stronger than our initial restriction (38). If we apply (42) to a numerical solution of the task (39), then the calculation appears to be stable—see Fig. 16(b). In fact, condition (41) is a generalized variant of the Courant–Friedrichs–Levi stability condition, which takes into account the interaction of convection and diffusion. Due to clear physical sense, this condition generalizes easily to the full Favre averaged Navier–Stokes equation system. 2.4. Selection of the method for an approximation of source terms In differential equations for turbulence parameters, source terms that describe production and dissipation of turbulence appear. To substantiate our approach to an approximation of these terms, let us consider a model system of equations with source terms. To derive this model system, let us take the last two equations (for turbulence parameters) from the equation system (4)+(12), subtract from each equation the continuity equation, multiplied by the corresponding parameter of turbulence (q or o), and finally exclude the spatial terms (fluxes). We shall obtain the following model system of equations: " #

 SðqÞ=r q d~ p ~; ~ ~ ¼ ¼W p¼ ; W . (43) SðoÞ=r dt o Now let us linearize the source terms in the vicinity of some point of space–time. For the following analysis it will be useful to take for this point the state in the current grid

91

cell (i, j, k) at the known time layer n: !n ~ qW n ~ ~ ð~ p ~ pni;j;k Þ. W  W i;j;k þ q~ p

(44)

i;j;k

~ on other Here we have neglected the dependence of W parameters (mainly, on velocity gradients, see (18)). We have used the approximation that these other parameters are equal to their values at the time layer n. Introducing a new variable 2 !n 31 ~ q W 5 W ~n ¼ ~ ~ p¼~ p ~ pni;j;k þ 4 p þ const, i;j;k q~ p i;j;k

we convert system (43) to the following form: d~ p ¼ R ~ p; dt

R ¼ const,

where the matrix !n ~ qW R¼ ¼ q~ p i;j;k

~ qW q~ p

(45)

!n . i;j;k

Let us solve (45) as a system of ordinary differential equations. We may represent R as R ¼ TKT1, where K ¼ diag(li) is a diagonal matrix consisting of eigenvalues li of matrix R, and T ¼ colð~ ei Þ is the matrix with columns coinciding with eigenvectors ~ ei of matrix R. Multiplying (45) by T1 from the left and introducing ~ z ¼ T1~ p, we shall obtain the system d~ z=dt ¼ K~ z. But this system consists of independent scalar differential equations (dzi/dt) ¼ li zi with the solution zi ¼ zi0 eli t . Finally, we should return to the variable ~ p using formula ~ p ¼ T~ z. Consequently, the exact solution of (45) has the form: X ~ ~ ei zi0 eli t , pðtÞ ¼ (46) i

z0 that is defined by the where zi0 is ith component of vector ~ relation: ~ z0 ¼ T1~ pð0Þ. Consequently, due to the presence of source terms the solution of the time-averaged Navier–Stokes equation system contains modes, which vary exponentially in time. If lio0, the corresponding mode of solution (46) damps. If li40, the corresponding mode of solution grows. An approximation of source terms should have at least second order in time. This requirement is connected not only with the need to keep the approximation order of our numerical method. A more essential problem is that due to the presence of exponentially varying modes the errors of first order, growing exponentially, will strongly deteriorate the solution. As in Section 2.1, let us consider two methods for an approximation of source terms—explicit and implicit schemes of second-order accuracy in time.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

92

Initially let us consider an explicit approximation of the predictor–corrector type: 8 n n < ~p ~p ¼ R ~ pn ; t (47) : ~pnþ1 ~pn ¼ R ~pn þ~pn : t 2 System (47) may be rewritten in the form   t2 2 n nþ1 ~ p ¼ E þ tR þ R ~ p , 2

i

where qi ¼ 1 þ tli þ ð1=2Þt2 l2i . Designation ‘‘(qi)n’’ means qi raised to a power n (not the value at the time layer n). If t-0, then ½1 þ tli þ ð1=2Þt2 l2i n ! eli nt , and (48) converts to (46). But we work with finite time steps. For the scheme (47) at each value of t the geometric ratio qi40, and therefore the modes of numerical solution are always positive, and so in eli t . But we should also require that the modes of the numerical solution, as well as modes of the exact solution, should damp for lio0 and should grow for li40. If li40, then qi41, and the mode of the numerical solution grows, similar to the exact solution. But if lio0, we should require qio0. It is easy to see that to satisfy this condition one should take to2/|li|. If this condition will be violated even for one mode of solution, then we shall obtain an exponentially growing mode instead of an exponentially damping one. This would mean instability of the numerical scheme. Consequently, the explicit scheme (47) is stable under the condition 2 , max jl i j

qi ¼

(50)

where are negative eigenvalues of the matrix R. Now let us consider the possible situation, when max jl i jb i mini jl i j (by two orders of magnitude or more). Then the main part of the exact solution (46) will belong to the modes with a minimal value of jl i j, such modes will damp the most slowly. On the contrary, the modes with maximal value of jl i j damp very quickly and disappear from the sum (46). The characteristic time of the whole process may be estimated as tchar 1=ðmini jl i jÞbtstab . From a practical viewpoint, it would be enough for the description of the whole process to select a time step so that within the time tchar we would have 5–10 time steps. But the stability condition requires to move in time much more slowly, with the step tstab, which guarantees the stability of the most insufficient of the modes of solution. As a result, computational cost may be too large. This problem is

1 þ ðt=2Þli . 1  ðt=2Þli

It is easy to see that for lio0 follows |qi|o1. Consequently, the problem of stiffness does not appear: if some modes damp in reality, then they damp in the numerical solution, too. But it is possible that qio0, and then the mode of the numerical solution will behave non-physically (i.e. different from the corresponding mode of the exact solution, eli t ): it will change sign at each time step (however, modes of the exact solution are always positive). In addition, at lio0 it is possible that qio1, and in this situation the mode will oscillate and will grow to infinity. This behaviour may also be considered as instability of a numerical scheme. To formulate an appropriate time step restriction for an implicit scheme, we shall consider two possible variants. Variant 1: All lio0, i.e. all modes of the exact solution will damp. Then the main part of the exact solution belongs to modes, which correspond to eigenvalues liA0(10lmin,lmin), where lmin is the negative eigenvalue with the minimal modulus. Then we shall impose the following requirements on the modes of numerical solution: (1) The main modes (liA(10lmin,lmin)) should be positive (as well as corresponding modes of the exact solution—eli t ) and should damp, i.e.

0o

(49)

i

l i

~ ~ pn pnþ1 pnþ1  ~ pn þ ~ ¼R . t 2

The solution of this system may also be represented in the form (48), but for this scheme we have

where E is the unit matrix. After that, we may perform the same transformations as in the derivation of the exact solution of system (45); during these transformations, we should take into account that R2 ¼ TK2T1. As a result, we shall obtain the exact solution of the system of recurrent Eq. (47): X ~ ~ ei zi0 ðqi Þn , pn  ~ pðtn Þ ¼ (48)

totstab ¼

called the problem of stiffness of the differential equation system [61]. Now let us consider an implicit approximation of the second order in time:

1 þ ðt=2Þli o1. 1  ðt=2Þli

(51)

(2) Insignificant modes (lip10lmin) should only damp; but we may allow ourselves to describe them incorrectly, i.e. we may allow them to change sign in time. (In any case, these modes damp very quickly.) So, we shall impose the following requirement for these modes:

1o

1 þ ðt=2Þli o1. 1  ðt=2Þli

(52)

As it was said before, for lio0 condition (52) is always satisfied, and condition (51) leads to the following time step restriction: tp

2 1 ¼ . j10lmin j j5lmin j

(53)

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Variant 2: There exist modes with liX0, i.e. nondamping modes. In this case damping modes are not so important, and we may only require that their numerical analogues should also damp. However, for numerical analogues of non-damping modes we shall require both a positive sign and a non-damping character: li o0 ) 1o

1 þ ðt=2Þli o1; 1  ðt=2Þli

li X0 )

1 þ ðt=2Þli X1. 1  ðt=2Þli

These conditions are satisfied at to

2 , max lþ i i

where lþ i are positive eigenvalues of matrix R. To increase the quality of description for the most quickly growing mode, we shall require tp

1 . 5max lþ i

(54)

i

Depending on the situation (variant 1 or variant 2), it is necessary to use either condition (53), or condition (54). Both conditions may be merged in a single condition: tp

1 , j 5max li j

(55)

i

where the maximum is taken not in the modulus, but in the immediate value. The disadvantage of implicit schemes is the necessity to solve the system of linear equations to find ~ pnþ1 . But in our case this is a system of the dimension 2  2, and it may be solved analytically. As a result, the growth of operations per time step in comparison with explicit scheme becomes insignificant. But in stiff situations, when there are damping modes, for which |li| is two or more orders of magnitude higher than jmaxi li j, the implicit scheme (50) will be essentially more effective than the explicit scheme (47). That is why in our numerical method the implicit approximation of source terms as given by Eq. (50) has been selected. 3. Approximation of Favre averaged Navier–Stokes equations

Rodionov scheme (GKR) is used. This scheme was described in our previous paper [4]. Here we shall give only a brief description of this scheme and of its new features. For the predictor step calculations are performed along Eq. (4), where the convective fluxes are calculated by formulas ~n Þ. The function F~conv ðPÞ ~ is given by ~i1=2 ¼ F~conv ðP like F i1=2 n ~ (12), and Pi1=2 are calculated by the following formulas: !n ~ q P n n ~ ~ P DLiþ1=2 xi , iþ1=2;j;k ¼ Pi;j;k þ qxi ijk !n ~ qP ~n ~n P DR ð56Þ i1=2;j;k ¼ Pi;j;k  i1=2 xi . qxi ijk

Derivatives in (56) are calculated using a nonlinear limiter function, which analyses the distribution of parameters along the index axis and limits the values of gradients to obtain the monotonic property of the scheme:     Di1=2 ðPm Þn Diþ1=2 ðPm Þn qPm n  limiter ; , (57) qxi i;j;k Di1=2 xi Diþ1=2 xi where Pm is the m-component of the vector of the flow ~ As a limiter function, we use a classical parameter P. ‘‘minmod’’ function [62,63] or a function as proposed by van Leer [37]. In our previous paper, a limiter function was separately ~ In particular, it applied to each component of the vector P. was applied to each component of the velocity vector ~ ¼ ðu; v; wÞ, as to independent values. V Now we have introduced an important modification. It is obvious that the velocity vector is an indivisible physical object. That is why the limiter function is now applied to the modulus of a velocity vector. It is not applied to the ~. Indeed, a variation of the direction of the direction of V ~ at constant jV ~j does not require the variation of vector V gasdynamic parameters at a given point and therefore does not produce any non-physical oscillations of these parameters during the calculation. The new algorithm consists of three stages: At the first stage we determine the directions of the velocity vectors at opposite sides of the (i; j; k) cell: ~ ei1=2;j;k ¼

3.1. Special features of an approximation of convective fluxes In Section 2, model equations were used to illustrate our decisions. Now let us describe the scheme for the solution of the full Favre averaged Navier–Stokes system (4)+(12). In accordance with the principle of inheritance (see the beginning of Section 2), the terms, which pass to the Navier–Stokes equations from the Euler equations (conconv vective fluxes F~ in (12)) are approximated in the same way as in the inviscid case—the Godunov–Kolgan–

93

~0 V i1=2;j;k 0

~ jV i1=2;j;k j

; ~ eiþ1=2;j;k ¼

~0 V iþ1=2;j;k 0

~ jV iþ1=2;j;k j

,

where ~0 ~ ~0 ~0 ~ ~0 V i1=2 ¼ Vi;j;k  DVi;j;k ; Viþ1=2 ¼ Vi;j;k þ DVi;j;k , ~ 0 þ DR V ~0 Þ, ~0 ¼ 1ðDL V DV i;j;k 2 ~ ~ ~0 ¼ Di1=2 V DR xi ; DR V ~0 ¼ Diþ1=2 V DL xi . DL V i1=2 Di1=2 xi Diþ1=2 xi iþ1=2 0

~ j ¼ 0, we take ~ e ¼ 0. If jV

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

94

At the second stage the limiter function is applied to the lengths of vectors:    1 0 ! ~ ~ D D V   ~ i1=2 iþ1=2 V qjVj @ A. ¼ limiter ; qxi Di1=2 xi Diþ1=2 xi i;j;k

At the third stage, knowing the velocity directions, we calculate the vectors of velocity at the opposite sides: 2 3 ! ~ ~i1=2 ¼ ~ ~i;j;k j  qjVj 5 V ei1=2;j;k 4jV DR i1=2 xi , qxi i;j;k 2 3 ! ~ qj Vj ~iþ1=2 ¼ ~ ~i;j;k j þ V eiþ1=2;j;k 4jV DLiþ1=2 xi 5. qxi i;j;k

Finally we calculate the velocity gradient as follows: ! ~i1=2 ~iþ1=2  V ~ V qV ¼ R . qxi Di1=2 xi þ DLiþ1=2 xi i;j;k

For scalar values (p; T; q; o) we apply the standard formula (57). Our experience shows that the use of this algorithm for velocity vectors results in an essential improvement of the scheme quality. Any dependency of the numerical solution on the turn of a reference frame disappears. Non-physical near-wall entropy layers are getting thinner. Strong discontinuities become less sensitive to the sharp turns of grid lines. Also the smoothing of shock waves diminishes while the convergence to steady state accelerates. The scheme becomes more stable in critical situations. In this paper the values of gradients, which are used in the GKR scheme and are calculated using limiter functions, will be called ‘‘convective’’ gradients. It is worth noting that ‘‘convective’’ gradients are calculated only with first-order approximations—this is enough to obtain second-order approximations in calculation of flow parameters. For the corrector step calculations are again performed using Eq. (4), but now convective fluxes are calculated by

~conv ðP ~nþ1=2 Þ. To obtain values of the formula F~i1=2 ¼ F i1=2 flow parameters at the cell border (i+1/2, j, k) at time layer (n+1/2), we initially find flow parameters in the cell centres ~nþ1=2 ¼ ð1=2ÞðP ~n þ P ~ Þ. Then at time layer (n+1/2): P i;j;k

i;j;k

i;j;k

we assume that at the time layer (n+1/2) the parameters are linearly distributed within each cell, and in addition we take !nþ1=2 !n ~ ~ qP qP  . qxi qxi ijk

ijk

Then at the cell borders we have discontinuities in parameter distributions, and parameters at the cell borders ~nþ1=2 are calculated from the solution of the Riemann P i1=2

problem about the decay of an arbitrary discontinuity [34,64]: ~L ; P ~R Þ, ~nþ1=2 ¼ decayðP P iþ1=2;j;k where !n ~ qP DLiþ1=2 xi , qxi ijk !n ~ ~R ¼ P ~nþ1=2  qP P DR iþ1;j;k iþ1=2 xi . qxi ~L ¼ P ~nþ1=2 þ P i;j;k

ð58Þ

iþ1;j;k

In our method we use the exact iterative self-similar solution of the 1D Riemann problem about the decay of discontinuity between two half-infinite uniform flows of an inviscid gas (Fig. 17) [34]. The use of this procedure provides high reliability and quality of results. Switching to time-averaged Navier–Stokes equations, the question arises: is it possible to use for turbulent flow of a viscous gas the solution of the 1D self-similar solution of the Riemann problem, obtained in the framework of Euler equations? It may be shown that in the GKR scheme the solution of the Riemann problem is considered at the very initial

Fig. 17. Typical structure of the exact solution of the Riemann problem.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

moments after the decay of the discontinuity (t-tn+1/2+0). This is due to this fact that we may use the solution of the self-similar Riemann problem in a situation, when flows from both sides of the cell border are not uniform. Let us note that in the inviscid limit (m-0) the structure of the discontinuity decay arises instantly. At the initial time t ¼ tn+1/2 we have a discontinuity of gas parameters at the cell border, but at any indefinitely small t-tn+1/2+0 we already have two waves and the contact discontinuity (Fig. 17). One may assume that at small (ttn+1/2) we may neglect the actions of gradual processes—molecular and turbulent diffusion and production/dissipation of turbulence. Therefore we may expect that initially the flows from both sides of the discontinuity interact in accordance with the laws of continuous medium, and the ‘‘inviscid’’ structure of discontinuity decay arises; and only after that this structure begins to be smoothed by diffusion, and turbulence begins to develop. Consequently, in the framework of the basic assumptions of the GKR scheme, we may neglect the action of viscosity and turbulence and use solution of the 1D selfsimilar Riemann problem for an inviscid gas. In fact, the use of the 1D Riemann problem implies that the flow parameters vary only along the normal to the cell border. For such a flow, the tangential-to-the-border ~t ¼ V ~  ðV~ ~nÞ~ component of gas velocity ðV nÞ and also turbulence parameters (q and o) are passively transported values. Passive values have no influence on the structure of the solution of the Riemann problem. They are transmitted without variation along the trajectories of gas volumes. Therefore, if in the solution of the Riemann problem the contact discontinuity passes to the right from the cell ~t , q and o at the border (as in Fig. 17), then the values of V cell border after decay of the discontinuity should be equal to their values from the left of the cell border before the decay; if the contact discontinuity passes to the left of the cell border, we should take their values after the decay from the state right of the cell border before the decay. 3.2. Approximation of diffusive fluxes Now we shall consider the approximation of terms in the system (4)+(12), which are connected with molecular and turbulent transport. These terms appear in scheme (4) as diffusive fluxes through the borders of grid cell; each diffusive flux is a linear combination of the following expressions:   m mT qf l þ , Pr PrT qxm where fl are parameters from the set f~ ¼ ðT; u; v; w; q; oÞ. The model of decay of discontinuity between two semiinfinite flows cannot be used in the approximation of these terms: in our model, the solution of the Riemann problem has piece-wise-constant character, and almost everywhere qfl/qxm ¼ 0. In real flow, there are no decays of discontinuities at the cell borders. A decay of discontinuity is only a

95

model that provides the correct distribution of the information, which comes to the cell border along characteristic lines from the flows from both sides of the border. This model cannot serve for diffusive fluxes. Let us consider for example the border (i+1/2, j, k) (see Fig. 11). We shall use the curvilinear coordinate system, connected with grid lines (xi, xj, xk) (see Section 1.2). The derivatives at the border will be calculated by formula: qf~ qf~ qxi qf~ qxj qf~ qxk ¼ þ þ . qxm qxi qxm qxj qxm qxk qxm

(59)

Approximations of the metric coefficients are described in Section 1.2—formulas (6)–(9). Therefore, now we should describe the approximation of the derivative along a grid line, which intersects the border, ðqf~=qxi Þ, and also approximations of derivatives along tangential-to-border directions qf~=qxj and qf~=qxk . On the basis of ideas from Section 2.1, we shall take for ðqf~=qxi Þ the approximation that would convert to the scheme (23) and (24) for 1D diffusion equation (19) on a uniform grid. In the case of a uniform grid along the grid direction xi, this is the central-difference approximation: ! Diþ1=2 f~ qf~  . Diþ1=2 xi qxi iþ1=2

But we should obtain the second-order approximation. This formula provides a second-order approximation only at a uniform grid, when the cell border is placed at the centre of the grid line xi element that connects the centres of cells (i, j, k) and (i+1, j, k) in the point ðxi ÞC ¼ ð1=2Þ½ðxi Þi;j;k þ ðxi Þiþ1;j;k . On a non-uniform grid, this formula should be corrected. This may be done as follows: ! ! !   qf~ qf~ q2 f~ ¼ þ ðx Þ  ðx Þ i iþ1=2 i C qxi qxi qx2i ðx Þ iþ1=2 ðxi ÞC i C ! !  qf~ q2 f~ 2 þ O ðDxi Þ ¼ þ qxi qx2i ðx Þ ðxi ÞC i C   1 L 2 R  Diþ1=2 xi  Diþ1=2 x þ O ðDxi Þ . 2 For the calculation of ðq2 f~=qx2i Þðxi ÞC we may use the values qf~=qxi in the centres of cells, adjacent to this border. The following approximation arises:  .  ! ~ Diþ1=2 f~ Diþ1=2 qf qxi qf~  þ Diþ1=2 xi qxi Diþ1=2 xi iþ1=2   1  DLiþ1=2 xi  DR x ð60Þ iþ1=2 i . 2 For a calculation of derivatives along tangential-toborder directions, we also may use the values of derivatives

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

96

in the cell centres: 0 ! ! qf~ qf~ @  interpolation qxj qxj iþ1=2

! qf~ qxj !

0

 interpolation@ iþ1=2

qf~ qxk

; i;j;k

1

A ; DLiþ1=2 xi ; DR iþ1=2 xi ,

iþ1;j;k

qf~ qxk

But we have full sets of derivatives along all grid directions in the cell centres—they are calculated with second-order approximation by formula (62). So, the natural idea appears—to replace formula (60) for the derivative in the intersecting-the-border direction by the same formula as for tangential derivatives (61):

!

@f~ @xk

! qf~ qxi

! , i;j;k

iþ1=2

1

A ; DLiþ1=2 xi ; DR iþ1=2 xi .

0

! ! qf~ qf~ @ ,  interpolation ; qxi qxi i;j;k iþ1;j;k ! DLiþ1=2 xi ; DR iþ1=2 xi .

ð63Þ

ð61Þ

iþ1;j;k

For closing approximations (59)–(61) we should describe the means of calculating the derivatives in the cell centres. Formulas (60) and (61) provide second-order approximations only under the condition that derivatives in the cell centres are also approximated with second order. In our numerical method the derivatives in the cell centres are calculated using the following formula:    Di1=2 f l Diþ1=2 f l qf l  interpolation ; ; qxm i;j;k Di1=2 xm Diþ1=2 xm  Di1=2 xm Diþ1=2 xm ; . ð62Þ 2 2

This formula also provides a second-order approximation, but together with (61) it gives an absolutely uniform algorithm for the calculation of all derivatives. Unfortunately, an attempt to use formula (63) leads to unexpectedly bad results. As an example, in Fig. 18(a) one may see the profile of boundary layer at the wall, which was calculated using formula (63), and in Fig. 18(b), the profile obtained using formula (60). It is easy to see that formula (63) leads to a step-like form of the profile; however, formula (60) achives a smooth profile. For an explanation let us apply formulas (60) and (63) for an approximation of the 1D scalar model Eq. (19) on a uniform grid.

On a uniform grid, this formula converts to a central difference. The derivatives, calculated in this way, are used only in the ‘‘diffusive’’ part of the scheme. Below we shall call them ‘‘diffusive’’ gradients. During calculation, we also have to calculate and store in the memory one more set of gradients, ‘‘convective’’ gradients of parameters, which are used in the reconstruction of flow distribution over the cells before the solution of the Riemann problem (see Section 3.1). ‘‘Convective’’ gradients are calculated using different formulas, which include limiter functions, (57). These formulas cannot be used in the approximation of diffusive fluxes, because they are asymmetric and give values of gradients only with first order of approximation. For the corrector step, the values of the parameters f~ in (60)–(62) are calculated at the time layer (n+1/2). For the predictor step it would be possible to take parameters from the time layer n, as well as in scheme (23). But the predictor is an intermediate step, which should give a solution only with the 1st order approximation in time; this fact allows us to exclude additional calculations of gradients and take for the predictor step the values of gradients from the time layer (n1/2), which were calculated for the corrector of the previous time step. (Difficulties with the first time step do not occur, because the value of the first time step is always taken to be equal to zero.) The approximation of diffusive fluxes through the cell border (60) and (61) is not uniform: for the derivative in the intersecting-the-border direction and for tangential derivatives we use different formulas.

Fig. 18. Influence of the method for derivative approximation on the form of the boundary layer profile: (a) formula (63) and (b) formula (60).

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Formula (60) generates the following approximation of equation (19): uinþ1  uni ui1  2ui þ uiþ1 ¼D (64) t h2 and formula (63) generates the following approximation: uinþ1  uni ui2  2ui þ uiþ2 ¼D . (65) t 4h2 We may see the first disadvantage of using formula (63): it doubles the stencil of the scheme and, consequently, diminishes its resolution. Scheme (64) resolves details of the solution with wavelengths lX2h, scheme (65), with lX4h. As a consequence, scheme (65) is more dissipative. But scheme (65) has a more serious disadvantage: it does not include the adjacent cells (i1) and (i+1). As a result, in the calculation on a uniform grid the cells with even numbers and the cells with odd numbers do not interact at all, and the solution within these two sets of cells develops independently. That is why we have obtained a step-like profile in Fig. 18(a). Thus, in scheme (65) there is no even–odd interaction. However, scheme (64) provides such interaction. That is the reason, why we have realized in our numerical method the approximation of diffusive terms, based on formula (60), and have rejected formula (63). Finally, we also have to calculate the molecular and turbulent viscosity at the border (i1/2, j, k), and also the term (2/3)rq2dij, that is included in the turbulent fluxes of momentum. According to formulas (17) and (10), for this purpose one should know the full set of parameters at the ~iþ1=2 ¼ ðT; u; v; w; p; q; oÞiþ1=2 . It would be cell border, P possible to take these parameters from the solution of the calculation of the decay of an arbitrary discontinuity. But the Riemann problem solution is principally an asymmetric procedure that is determined by the convective interaction of fluxes from both sides of the cell border. This fact contradicts the isotropic character of diffusion, which gives a symmetrical propagation of perturbations in all directions. Sometimes the use of values from the Riemann problem leads to the following situation: at one time step ~iþ1=2 is taken from the left side of the border, but the value P at the next time step, from the right side. As a result, nonphysical self-excited oscillations arise. So, we use in diffusive fluxes through the cell border the following values of flow parameters: ~diff ¼ ðT; u; v; w; p; q; oÞdiff P iþ1=2 iþ1=2 ~i;j;k ; P ~iþ1;j;k ; DL xi ; DR xi Þ: ¼ interpolationðP iþ1=2 iþ1=2 ð66Þ ~iþ1;j;k in (66) are ~i;j;k ; P For the predictor step, the values P taken from the time layer n; for the corrector step, from the time layer (n+1/2). Parameters at the cell centres at the time layer (n+1/2) are calculated as in the GKR scheme, as half the sum of the value at the time layer n and of the value at the time layer (n+1), which was obtained for the predictor.

97

It is necessary to note that the parameter q appears both in convective and in diffusive parts of the vector of fluxes ~ (12). In these parts we use through the cell border F different values for the parameter q. In the convective flux of the parameter q (term Qq in (12)) and in the expression for the total energy of gas E (that appears in the convective flux of energy) we use the value of q from the solution of the Riemann problem at this boundary. However, in the formula for turbulent viscosity (10), in the formula for Reynolds stresses (13) and in the formula for fluxes of turbulence kinetic energy (16) we use the value of q from ~diff (66). the interpolated vector P iþ1=2

The distance to the wall (yW)i+1/2 is also calculated by an interpolation from the centres of adjacent cells. Values of yW in the cell centres should be calculated anyway, because they appear in the function Fwall in the source terms (18). These values are calculated by formula ðyW Þi;j;k ¼ jð~ ri;j;k  ~ rW Þ~ nW j, where ~ rW is the radius vector of the centre of the cell border on the wall surface, which is closest to the cell (i, j, k), and ~ nW is normal to this cell border. 3.3. Approximation of source terms For the source terms we use the local-implicit approximation that was selected in Section 2.4. Here the whole realization of this scheme will be described. For source terms we use expressions (18), where the coefficients Aq, Bq, Ao, Bo depend on the values of ‘‘diffusive’’ gradients, and the coefficients Cq, Co are constant. For the predictor step we use the values of ‘‘diffusive’’ gradients from the previous step, from time layer (n1/2), and an explicit approximation of source terms with first order in accuracy in time is used:

 n n1=2 1 n n1=2 n SðqÞ ¼ r Aq þ Bq þ C q o qn , on h i n1=2 n1=2 n þ Bo o þ C o ðon Þ2 . SðoÞn ¼ rn Ao For the corrector step, new values of ‘‘diffusive’’ gradients are calculated, using parameters from the time layer (n+1/2). Then with the use of new ‘‘diffusive’’ gradients, new values of coefficients are calculated Aqnþ1=2 ; Bqnþ1=2 ; Aqnþ1=2 ; Bqnþ1=2 . Parameters of turbulence at the time layer (n+1) are obtained from solution of the following system of nonlinear equations, which are based on (3): 8  < rnþ1 qnþ1 ¼ rn qn  t DF ðqÞnþ1=2 =V þ t SðqÞn þ SðqÞnþ1 ; 2  : rnþ1 onþ1 ¼ rn on  t DF ðoÞnþ1=2 =V þ 2t SðoÞn þ SðoÞnþ1 ;

(67) nþ1=2

nþ1=2

; DF ðoÞ are differences of fluxes where DF ðqÞ through the cell borders (fluxes are multiplied by the vectors of the cell border area, see (12)); V is the cell volume; in the formulas for S(q)n, S(q)n+1, S(o)n, S(o)n+1 we used the new values of coefficients Aq ; Bq ; Ao ; Bo .

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

98

Scheme (67) is an implicit scheme, based on (50); contrary to (44), " # " # r q SðqÞ ~ ~ p¼ and W ¼ . r o SðoÞ The essential difference of scheme (67) from implicit schemes, which were considered in Section 2.1, is that (67) is a point-implicit scheme: equation system (67) is not connected with unknown parameters at the time layer (n+1) in adjacent cells. As a result, there is no need to solve a giant system of nonlinear equations for all cells of the computational grid simultaneously; in each cell the system of two nonlinear equations is solved. The nonlinear system (67) is solved approximately, using one iteration of the Newton method (that is equivalent to a linearization (44)). This approximation is justified by the fact that the (qo)-model itself is only an approximate semi-empirical mathematical model for the description of ~ =q~ turbulence. An inversion of the matrix R ¼ ðqW pÞni;j;k is performed analytically. Initially with the use of new values of the coefficients Aq ; Bq ; Ao ; Bo the following values are recalculated:

 n nþ1=2 1 n nþ1=2 n þ Bq þ C q o qn , SðqÞ ¼ r Aq on h i nþ1=2 n n 2 þ B o þ C ðo Þ SðoÞn ¼ rn Anþ1=2 . o o o After that, the values RHS1 ¼ SðqÞn  DF ðqÞnþ1=2 = V RHS2 ¼ SðoÞn  DF ðoÞnþ1=2 =V are determined. Finally, the values rnþ1 qnþ1 and rnþ1 onþ1 are calculated by the following formulas: (

t ½a22 RHS1  a12 RHS2 ; rnþ1 qnþ1 ¼ rn qn þ det t a11 RHS2 ; rnþ1 onþ1 ¼ rn on þ det

where a11 ¼ 1(t/2)l1, a22 ¼ 1(t/2)l2, a12 ¼ ðt=2Þ½C q  ðAqnþ1=2 =ðon Þ2 Þqn , det ¼ a11a22 and the eigenvalues of the matrix R are equal to l1 ¼

Aqnþ1=2 þ Bnþ1=2 þ C q on ; q on

nþ1=2 n l2 ¼ Bnþ1=2 þ 2C o o . o

(68) (The matrix R is calculated at the time layer n instead of the layer (n+1/2), because the values of parameters at the time layer (n+1/2), obtained by a simplified procedure for the predictor step, are less reliable (in conditions of exponentially varying modes) than the values of parameters at the time layer n. But values of the coefficients Aq ; Bq ; Ao ; Bo algorithmically may be calculated only at the time layer (n+1/2); so, we assume that these coefficients are constant during the whole time step and equal to Aqnþ1=2 ; Bnþ1=2 ; Aqnþ1=2 ; Bqnþ1=2 . (This allows us to calculate q the values of S(q)n and S(o)n also with the use of values Aqnþ1=2 ; Bnþ1=2 ; Aqnþ1=2 ; Bqnþ1=2 ). q

3.4. Stability condition of the scheme Let us consider a quasi-1D flow, in which all parameters vary only along the coordinate axis n, but along coordinate axes t1 and t2, perpendicular to n axis, the flow is constant: q/qt1 ¼ q/qt2 ¼ 0. For this case the equation system (4)+(12) may be transformed to the following form: 8 ! > > qz1 qz1 q2 z1 mT m 1 > > qt þ V n qn ¼ r þ A1 ; > Cp T þ Cp T qn2 > Pr Pr > T > > > > > qz2 qz2 mþmT q2 z2 > > ; > qt þ V n qn ¼ r qn2 > > > 2 > qz3 qz mþm q z > > þ V n qn3 ¼ r T qn23 ; > > qt <   qz4 qz4 q 2 z4 mT (69) 1 þ V ¼ m þ þ SðqÞ q n > r r ; qt qn qn2 Pr > T > > >   2 > > qz qz q z m > > qt5 þ V n qn5 ¼ r1 m þ PrTo qn25 þ SðoÞ > r ; > T > > > > qz6 qz mþm q2 z > þ ðV n þ cÞ qn6 ¼ 43 r T qn26 þ A2 ; > > qt > > > > > qz7 þ ðV n  cÞ qz7 ¼ 4 mþmT q2 z7  A2 ; : qt

qn

3

r

qn2

where z1 ¼ Cp lnR TR ln P, z2 ¼ RVt1, z3 ¼ Vt2, z4 ¼ q, z5 ¼ o, z6 ¼ u+ (dp/rc), z7 ¼ u (dp/rc) are the values that become Riemann invariants in the case of an inviscid flow. It is easy to see that in the case, when m ¼ mT ¼ 0, S(q) ¼ 0 and S(o) ¼ 0, the values z1, z2, z3, z4, z5 are constant along the characteristic lines dn/dt ¼ Vn (trajectories of flow volumes), whereas the values z6 and z7 are constant along the characteristics dn/dt ¼ Vn+c and dn/dt ¼ Vnc (trajectories of acoustic perturbations). One may also see that in the case of turbulent flow of a viscous flow the information about zi propagates over the space not only due to convection (along the characteristics), but also due to molecular and turbulent diffusion. In addition, even in an inviscid case the values z4 and z5 are not constant along characteristic lines due to the presence of source terms. The additional terms A1, A2 have the following form: !"    # 1 m mT 1 qT 2 q 1 qp þ þR A1 ¼ Cp r PrC p T PrC p T T qn qn p qn T "  2  2   # 1 m þ mT 4 qV n qV t1 qV t2 2 þ þ þ , T r 3 qn qn qn !   4 m þ mT q 1 qp 1 m mT c q2 T þ þ A2 ¼  3 r qn rc qn r PrC p T PrC p T T qn2 T "  2  2   # c m þ mT 4 qV n qV t1 qV t2 2 þ þ þ . CpT r 3 qn qn qn (the simplifying assumptions mEconst and mTEconst were made in the derivation of Eq. (69)). Now it is important for us that each equation of system (69) is similar to Eq. (37), if we exclude from our consideration (approximately) the cross-diffusive terms

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

A1, A2. So, we may use the stability condition (42) that was derived in Section 2.3 for an approximation of Eq. (37). Let us describe the realization of this stability condition for an arbitrary cell (i, j, k). When the Riemann problems about the decay of discontinuity are solved for the cell borders (i1/2, j, k) and (i+1/2, j, k), we know the speeds of the extreme waves, which propagate from these two borders inside the cell— Si1/2 and Si+1/2. Assuming that these are plane waves, propagating along the normal vectors to borders, let us require that in one time step these plane waves should not reach the centres of opposite borders (this is classical CFLcondition). For this purpose, the time step should not exceed the value: tconv

  

 ð~ riþ1=2  ~ ri1=2 Þ~ ni1=2  ð~ riþ1=2  ~ ri1=2 Þ~ niþ1=2   ¼ min  ;   . S i1=2 Siþ1=2

After that we shall calculate the diffusive restriction for the time step, taking the diffusion coefficient to be equal to the maximal of diffusion coefficients in system (69): tdiff ¼

1 4 

j~ riþ1=2  ~ ri1=2 j2

( 1 r max

4 3 ðm

þ mT Þ;

m PrC p T

þ

mT Pr

Cp T T

;mþ

). mT q PrT

;m þ

mT Pro T

Then using formula (42) we may calculate the restriction tNS i . In an analogous way, we calculate the time step restrictions along others grid lines: tNS and tNS j k . For 3D tasks, the time step restriction is calculated from these 1D restrictions using the following formula (see e.g. [64]): 1 . tNS ¼   NS . NS 1 ti þ 1 tNS j þ 1 tk In accordance with the conclusions from Section 2.4, our approximation of the source terms is stable under condition (55). In this condition, one should substitute the eigenvalues from (68). As a result, we obtain the time step restriction, connected with the source terms: tsource ¼

1 . j5maxðl1 ; l2 Þj

Then the final time step restriction is calculated by formula: tmax ¼ minðtNS ; tsource Þ. 4. Numerical boundary conditions 4.1. General principles of the numerical boundary condition formulation Calculations are performed on a multi-block regular computational grid. At the boundaries between the blocks the grid lines may be either continuous or discontinuous. Each block of computational domain is a curvilinear

99

hexahedron (possibly, with degenerated boundaries). The numerical method described above nominally has secondorder approximation in all variables. Its stencil includes 25 cells: (i, j, k), (i71, j, k), (i, j71, k), (i, j, k71), (i72, j, k), (i, j72, k), (i, j, k72), (i71, j71, k), (i71, j, k71), (i, j71, k71). When the stencil, corresponding to the cell (i, j, k), intersects the boundary of the block of the computational domain, the problem of numerical boundary condition arises: we should in some way fill up the missing values of flow parameters, which correspond to the stencil points that fall outside the boundary of the block of the computational domain. One way to solve this problem is to store two layers of fictitious out-of-boundary cells in the computer memory, which are filled in accordance with physical properties of each boundary of the block. An alternative way is to store in memory not only values of flow parameters, but also values of their gradients along grid lines. In this case it is enough to use only one layer of out-of-boundary cells. The advantage of these two approaches is that after the filling of the out-of-boundary cells the computation in all inner cells of the block is performed with the same algorithm. In our numerical method, another approach is realized. In the cell centres we also store the values of flow parameters and their gradients along grid lines. But out-of-boundary cells are introduced only for some types of boundary condition, in the cases, when it is natural and useful. In other cases they are not used. It means that in the first near-boundary layer of cells the calculations are performed with a special algorithm. For the sake of brevity, let us consider only the block boundary that is intersected by the axis xi; at that, this axis is directed from this boundary inside the block. Then the near-boundary cell has the number i ¼ 1, and its border, which lies at the block boundary, has the number i ¼ 1/2. The algorithm for the calculation of the flux through this boundary F~1=2 has the following form: Before the predictor step, the following actions are performed: 1. Initially, in the centre of the near-boundary cell (i ¼ 1) the ‘‘convective’’ derivatives of flow parameters are calculated along the grid lines, which intersect the boundary. Along the grid directions, tangential to the boundary, derivatives are calculated as usual, by formula (57). The algorithm for the calculation of the derivative in the direction of the axis i depends on the type of boundary condition and will be described below. 2. Subsequently, for all types of boundaries at the cell border i ¼ 1/2 the ‘‘diffusive’’ derivatives of the flow parameters are calculated, in the following order: At the beginning, all gradients along the grid directions are calculated. This procedure depends on the type of boundary condition. Then they are converted into gradients along Cartesian axes (x, y, z) by using formula (59). Finally, these gradients along

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

100

axes (x, y, z) may be corrected in the way, specific for each boundary condition (see below). The exception is the condition ‘‘Connect’’, for which the ‘‘diffusive’’ gradients at the cell border i ¼ 1/2 are calculated directly in Cartesian coordinates (x, y, z). 3. Finally, at the border i ¼ 1/2 two sets of flow ~conv Þ is used parameters are calculated. The first set ðP 1=2;j;k for the calculation of convective fluxes through this border. It is obtained by an extrapolation from the centre of the near-boundary cell (i ¼ 1) using ‘‘convective’’ derivatives, which were calculated by step 1:

~conv ¼ P ~n  P 1=2;j;k i;j;k

!n ~ qP qxi

as in step 3, but uses information from the time layer (n+1/2) instead of layer n. ~1=2 for the After that it is possible to calculate the flux F corrector step. It remains to describe the specific details of the realization of each boundary condition, the methods for the calculation of the following values:

 

DR 1=2 xi .



1;j;k

~diff Þ is used for calculation The second set ðP 1=2;j;k of diffusive fluxes through this border. The way of calculating this set depends on the individual type of boundary condition and will be described below. This information is sufficient for the calculation of ~1=2 in the predictor step. the flux F Between the predictor step and the corrector step, the following activities are performed: 4. Initially we calculate the flow parameters adjoining the border i ¼ 1/2 from the inner side, using the derivative ~ i Þn , which was calculated before the predictor ðqP=qx 1;j;k step:

~R Þ1=2 ¼ ðP

~nþ1=2 P 1;j;k



!n ~ qP qxi

DR 1=2 xi .

1;j;k

~R Þ1=2 we calculate the 5. Then using the obtained value ðP flow parameters adjoining to the border i ¼ 1/2 from ~L Þ1=2 . The algorithm for the the outer side ðP ~L Þ1=2 depends on the type of boundary calculation of ðP condition and will be described below. 6. Then in the centre of a near-boundary cell (i ¼ 1) the ‘‘diffusive’’ derivatives of the flow parameters along the grid lines are calculated. Along the tangential-toboundary grid directions the derivatives are calculated as usual—by formula (62). Derivatives along the axis i are calculated in different ways, depending on the type of boundary condition. 7. After that, the values of ‘‘diffusive’’ derivatives of flow parameters at the border i ¼ 1/2 are calculated. Algorithmically, this step coincides with the step 2, but uses information from the time layer (n+1/2) instead of layer n. 8. Finally, at the border i ¼ 1/2 two sets of flow ~conv is now parameters are calculated. The set P 1=2;j;k calculated by a solution of the Riemann problem: ~L Þ1=2 ; ðP ~R Þ1=2 Þ. The second set ~conv ¼ decayððP P 1=2;j;k diff

~ ðP 1=2;j;k Þ is calculated algorithmically in the same way

 

‘‘convective’’ derivatives in near-boundary cell at time ~ i Þn ; layer n; ðqP=qx 1;j;k ‘‘diffusive’’ derivatives along the i-axis in the nearnþ1=2 boundary cell at the time layer ðn þ 1=2Þ; ðqf~=qxi Þ1;j;k ; the full set of ‘‘diffusive’’ derivatives at the border i ¼ 1/2 at the current time layer, ðqf~=@xi Þ1=2;j;k ; ðqf~=@xj Þ1=2;j;k and ðqf~=@xk Þ1=2;j;k ; flow parameters adjoining the border i ¼ 1/2 from the ~L Þ1=2 ; outer side at the time layer ðn þ 1=2Þ; ðP flow parameters for the calculation of diffusive fluxes ~diff . at the border i ¼ 1/2 at the current time layer P 1=2;j;k

Below we shall describe for the numerical scheme only those boundary conditions that has been used in the simulation of the flow in ETW. 4.2. Boundary condition on the boundary between adjacent blocks In the calculation of the ETW flow, two variants of boundaries between adjacent blocks are used. In the first variant, grid lines are continuous at the boundary between adjacent blocks (see Fig. 19(a)). The corresponding numerical boundary condition is called ‘‘Joint’’. In the formulation of this boundary condition, one layer of out-of-boundary cells is introduced. These outof-boundary cells coincide with near-boundary cells of the adjacent block. Out-of-boundary cells are filled with flow parameters and their gradients from the near-boundary cells of the adjacent block. As a result, the boundary between the blocks becomes ‘‘invisible’’ for the numerical algorithm. In the second variant, the grid lines are not continuous at the boundary between blocks (see Fig. 19(a)). In this case, we use the boundary condition ‘‘Connect’’. Parameters at the outer side of the boundary are determined by an interpolation of parameters from near-boundary cells of the adjacent block. Let us describe the condition ‘‘Connect’’ in more detail. Among the cells of the adjacent block, we find the cell with the centre closest to the centre of the border i ¼ 1/2. The centre of that cell will be denoted as X. We assume that the grid line xi, which passes through the centre of the nearboundary cell i ¼ 1 and through the centre of the border i ¼ 1/2, also passes through the point X, as if the point X ~ i Þn has the index i ¼ 0. Then ðqP=qx 1;j;k is calculated by

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

101

Fig. 19. Boundary conditions named ‘‘Joint’’ and ‘‘Connect’’: (a) a grid near such boundaries; (b) the Mach number field in a calculation of supersonic flow around a ‘‘wedge–flat body’’ combination with the boundary condition ‘‘Connect’’.

the formula: !

  ðPm Þn1;j;k  ðPm ÞnX D3=2 ðPm Þn qPm n ¼ limiter ; . qxi 1;j;k D3=2 xi j~ r1=2  ~ rX j þ D R 1=2 xi Analogous actions are performed in the adjacent block. As a result, after step 1 (see Section 4.1) both in the centre of cell i ¼ 1 and in the one of cell X (in the adjacent block) the ‘‘convective’’ derivatives in all grid directions are ~L Þ1=2 , in the cell X ‘‘index determined. To obtain ðP coordinates’’ are introduced. These coordinates in the centre of the cell X take on the values (iX, jX; kX), corresponding to the number of the cell X in the adjacent block, in the centres of the cell borders, values (iX71/2, jX; kX), (iX, jX71/2; kX), (iX, jX; kX71/2), and in the cell vertices, values (iX71/2, jX71/2; kX71/2) (8 variants). Using the position of the border centre i ¼ 1/2 inside the cell X, we determine the values of ‘‘index coordinates’’ that correspond to this point (iX+Di, jX+Dj; kX+Dk), where the increments Di, Dj, Dk belong to the range ~L Þ1=2 is determined by the [1/2; 1/2]. Then the value ðP formula !n !n ~ ~ q P q P nþ1=2 ~L Þ1=2 ¼ P ~ ðP þ ðhi ÞX Di þ ðhj ÞX Dj X qxi qxj X X !n ~ qP þ ðhk ÞX Dk, qxk X

where (hi)X, (hj)X, (hk)X are the distances between opposite borders of the cell X along the grid lines xi, xj, xk; e.g. L R L ðhi ÞX ¼ ðDR i1=2 xi þ Diþ1=2 xi ÞX ¼ 2ðDi1=2 xi ÞX ¼ 2ðDiþ1=2 xi ÞX (here it is to be taken into account that the cell centre is

placed at the central point of the segment that connects the centres of opposite borders). For the calculation of the ‘‘diffusive’’ derivatives nþ1=2 ðqf~=qxi Þ1;j;k , the values of the flow parameters at the time layer (n+1/2) in the point X are used: 0 nþ1=2 !nþ1=2 nþ1=2 ~ f~1;j;k  f~X qf ¼ interpolation@ ; qxi j~ r1=2  ~ rX j þ DR 1=2 xi 1;j;k 1 nþ1=2 j~ r1=2  ~ rX j þ DR x D3=2 f~ D x i 1=2 3=2 i A ; . ; D3=2 xi 2 2 The boundary condition ‘‘Connect’’ is the only case for which the ‘‘diffusive’’ derivatives at the border i ¼ 1/2 are directly calculated along Cartesian axes (x, y, z). For this purpose, the derivatives along the grid directions are initially converted into derivatives along the Cartesian axes (x, y, z) both in the cell 1 and in the cell X (in the adjacent block) using formula (59), and subsequently interpolated at the border i ¼ 1/2: ! qf~ ¼ interpolation qxm 1=2;j;k

0

qf~ @ qxm

! X

qf~ ; qxm

1   R ; ~ r1=2  ~ rX ; D1=2 xi A.

! 1;j;k

An interpolation is applied to the derivatives along the Cartesian axes (x, y, z) instead of the derivatives along grid lines, because the grid directions in the cell (l, j, k) and in the cell X (in the adjacent block) may be very different. Flow parameters that are used in the expressions for diffusive fluxes (in the formula for turbulent viscosity and

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

102

so on), are calculated at the border by formula  ~diff ¼ interpolation P ~X ; P ~1;j;k ; P 1=2  j~ r1=2  ~ rX j þ DR 1=2 xi ; D3=2 xi (for the predictor step the parameters from time layer n are used, for the corrector step from layer (n+1/2)). This boundary condition ‘‘Connect’’ allows achieving good solutions even in the case, when a computational grid contains ‘‘overlaps’’ and ‘‘holes’’ in the region of joining the adjacent blocks (Fig. 19(b)). 4.3. The boundary condition ‘‘symmetry’’ We consider the flow in the ETW test section as symmetrical with respect to the longitudinal vertical plane. This simplification allows us to calculate only half of the flow field. At the symmetry plane, a special boundary condition is used. It would be equivalent to the condition ‘‘joint’’, if we would calculate the flow from both sides of the symmetry plane. This boundary condition is also used on the elements of the ETW solid walls, where the cells of the computational grid are too large to generate a classical boundary layer with a no-slip condition on the wall. In such situations this boundary condition simulates a solid wall with a slip flow. Let us consider an imaginary cell i ¼ 0, placed from outside the boundary symmetrically to a near-boundary cell i ¼ 1. The flow in this imaginary cell is symmetric to the flow in the cell i ¼ 1 at the current time layer: ~0 ¼ V ~1;j;k  2ðV ~1;j;k ; ~ V nÞ~ n, p0 ¼ p1,j,k, T0 ¼ T1,j,k, q0 ¼ q1,j,k, o0 ¼ o1,j,k. Derivatives in the near-boundary cell are calculated using the same formulas as for the inner cells: !   ðPm Þn1;j;k  ðPm Þ0 D3=2 ðPm Þn qPm n ¼ limiter ; , qxi 1;j;k D3=2 xi 2DR 1=2 xi   qf m nþ1=2 ¼ interpolation qxi 1;j;k 

! nþ1=2 ðf m Þ1;j;k  ðf m Þ0 D3=2 ðf m Þnþ1=2 R D3=2 xi ; ; D x ; . 1=2 i D3=2 xi 2 2DR 1=2 xi

~L Þ1=2 are obtained by symmetric The parameters ðP ~R Þ1=2 . reflection of the flow with parameters ðP For a calculation of ‘‘diffusive’’ derivatives at the border i ¼ 1/2 it is necessary to use the same algorithm as for the cell border inside the block, see Section 3.2. As the cells i ¼ 0 and 1 are symmetric with respect to the boundary, for the calculation of the derivative ðqf~=qxi Þ1=2 formula (60) simplifies to ! f~  f~ qf~ ¼ 1 R 0. qxi 2D1=2 xi 1=2

For derivatives along tangential directions one may use formulas (61). Parameters for the ‘‘diffusive’’ fluxes are ~diff ¼ ð1=2ÞðP ~0 þ P ~1;j;k Þ. calculated by formula P 1=2

4.4. The boundary condition ‘‘TW_const’’ On the walls of the ETW test section, a no-slip boundary condition with given surface temperature Twall is imposed. This boundary condition supposes the near-boundary cell to be placed in the region of an essential effect of viscosity, i.e. the Reynolds number ~t Þ1;j;k jDR xi r1;j;k jðV 1=2 m1;j;k

o 1.

~t Þ1;j;k ¼ V ~1;j;k  ðV ~1;j;k ; ~ (Here ðV nÞ ~ n is the tangential-towall component of the velocity in the near-boundary cell.) Additionally, it is assumed that the flow in near-boundary cell is slow and practically parallel to the wall, the velocity profile is linear, and the density is practically a constant. On the wall itself the mean velocity of the flow and its turbulent fluctuations should be equal to zero because of the absence of slipping. The turbulence parameter o is assumed to be constant in the near-boundary cell.4 Then we may define the following values for the flow ~W ¼ ~ parameters on the solid wall: V 0, TW ¼ Twall, T wall pW ¼ T p1;j;k , qW ¼ 0, oW ¼ o1,j,k. Arranging these 1;j;k ~W and f~W , we may calculate values into vectors P derivatives in the direction of axis i in the centre of the near-boundary cell as follows: !n ~ qP qxi

1;j;k

¼

~n  P ~n P 1;j;k W DR 1=2 xi

,

!nþ1=2 qf~ ¼ interpolation qxi 1;j;k

0

1 nþ1=2 nþ1=2 nþ1=2 DR xi D3=2 xi f~1;j;k  f~W D3=2 f~ 1=2 A. ; @ ; ; D3=2 xi 2 2 DR 1=2 xi

‘‘Diffusive’’ derivatives on the wall in direction of axis xi are calculated by the following formulas: ! ~1;j;k  0 qT  ~ V T 1;j;k  T wall qV ¼ ; ¼ , R qxi qxi 1=2 DR D1=2 xi 1=2 xi 1=2     q1;j;k  0 qo qq ¼ ; ¼ 0. qxi 1=2 qxi 1=2 DR 1=2 xi

4 Approaching the wall, the characteristic linear size of turbulent eddies diminishes (it cannot be larger than the distance to the wall), and the characteristic frequency of turbulence is growing. It tends to a limit value, which is connected with the minimal possible size of turbulent eddies (Kolmogorov microscale [8]). That is why the eventual variant of a boundary condition for the parameter o-(qo/qn)W ¼ 0 [16].

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

For the ‘‘diffusive’’ derivatives along tangential directions we use the formulas: ! ~ qV qxj



qq qxj



~ qV qxk

¼ 1=2



¼ 1=2

qq qxk

!



¼ 0; 

1=2

¼ 0; 1=2

   qT qT ¼ ¼ 0, qxj 1=2 qxk 1=2         qo qo qo qo ¼ ; ¼ . qxj 1=2 qxj 1;j;k qxk 1=2 qxk 1;j;k

After conversion of ‘‘diffusive’’ derivatives into derivatives along Cartesian axes (x, y, z) the gradient of parameter o is corrected to nullify the derivative along the normal to the boundary: ro ¼ ro  qo=qn ~ n, where qo=qn ¼ ðro ~ nÞ. In the calculation of convective fluxes at the border the ~conv and P ~diff are Riemann problem is not solved. Both P 1=2;j;k 1=2;j;k ~W . taken to be equal to P 4.5. The boundary condition ‘‘Riemann’’ To investigate the influence of the WT walls, calculations of the flow around the model in the ETW test section are compared to calculations of the model, placed into an ~1 (i.e. an infinite uniform airflow with given parameters P isolated model). In this case it is supposed that a free gas flow across the outer boundary of the computational domain is allowed. The formulation of a boundary condition at such boundaries depends on the flow direction (inflow/outflow) and on the Mach number along the normal to the boundary. Initially, we estimate at the boundary the velocity of gas along the normal to the boundary and the speed of sound. For this purpose we approximately solve (with the use of the non-iterative Roe algorithm [65]) the Riemann problem for the decay of discontinuity between parameters of undisturbed flow at infinity and parameters at the centre of near-boundary cell at the current time layer, and zero approximation to parameters at the boundary is obtained: ~ð0Þ ¼ decay ðP ~1 ; P ~1;j;k Þ. From this approximation, we P

invariant is transmitted from inside the computational domain to the boundary, then for the predictor step it is calculated using parameters at the centre of the near~1;j;k , and for the corrector step, using boundary cell P ~ parameters ðPR Þ1=2 . When the values of all Riemann invariants are determined, we find the first approximation ~ð1Þ . With the use of P ~ð1Þ to parameters at the boundary P 1=2

1=2

we calculate the derivatives along the axis i in the centre of near-boundary cell: !  ð1Þ ~m Þn ðPm Þn1;j;k  P1=2 D3=2 ðP qPm n ¼ limiter ; , qxi 1;j;k D3=2 xi DR 1=2 xi 0 1 !nþ1=2 nþ1=2 ð1Þ nþ1=2 f~1;j;k  f~1=2 D3=2 f~ DR xi D3=2 xi qf~ 1=2 A. ; ¼ interpolation@ ; ; qxi D3=2 xi 2 2 DR 1=2 xi



1;j;k

~L Þ1=2 are taken to be equal to P ~ð1Þ . The Parameters ðP 1=2 calculation of the ‘‘diffusive’’ derivatives of flow parameters at the border i ¼ 1/2 begins with an analysis of flow velocity along the normal to boundary Vn. If the gas flows inside the computational domain, then all ‘‘diffusive’’ derivatives along all directions are taken to be equal to zero (condition ‘‘Riemann’’ is based on the assumption that the outer flow is uniform). But if gas flows out from the computational domain, then the ‘‘diffusive’’ derivatives are determined by extrapolation from the inner cells of computational domain in the following way. For the ð1Þ calculation of ðqf~=qxi Þ we use again the parameters f~ : 1=2

! qf~ qxi

1=2

1=2

0 1 ð1Þ R f~1;j;k  f~1=2 D3=2 f~ DR 1=2 xi D1=2 xi þ D3=2 xi A @ ; ¼ extrapolation ; ; D3=2 xi 2 2 DR 1=2 xi

and in the calculation of derivatives along tangential grid directions we use the derivatives in the cell centres i ¼ 1 and i ¼ 2, which are already known ! qf~ qxj

1=2

exclude the values of Vn and c. After that, depending on values of Vn and c, the analysis of Riemann invariants is performed by the approximation of quasi-1D flow along the normal to the boundary. Let us introduce designations: ~ n-unit vector of outer normal to boundary, and n-variable in direction of this vector. We further assume that along characteristics dn=dt ¼ V n (gas trajectories) the values of the Riemann invariants z1 ¼ Cp ln T–R ln P, z2 ¼ Vt1, z3 ¼ Vt2, z4 ¼ q and z5 ¼ o are transmitted, and along the acoustic characteristics dn/dt ¼ Vn+c and dn/ dt ¼ Vnc (trajectories of small, acoustic perturbations) accordingly the values of the Riemann invariants z6 ¼ R R V n þ dp=rc  V n þ 2c=g  1 and z7 ¼ V n  dp=rc  V n  2c=g  1 are transmitted (these approximate expressions are strictly exact in isentropic flow only). If the Riemann invariant is transmitted from outside to the boundary of the computational domain, it is calculated using parameters of the undisturbed flow at infinity. If the

103

qf~ qxk

0

1=2

! 1=2

1 ! ! ~ ~ q f q f A ¼ extrapolation@ ; ; DR 1=2 xi ; D3=2 xi , qxj qxj 1;j;k 2;j;k 0 1 ! ! ~ ~ q f q f A ¼ extrapolation@ ; DR 1=2 xi ; D3=2 xi . qxk qxk 1;j;k

2;j;k

For the calculation of parameters for the diffusive fluxes, ð1Þ ~ P1=2 are used. 4.6. The boundary conditions at the inlet and outlet sections of ETW An important problem in the numerical flow simulation in the ETW test section is a correct formulation of boundary conditions at the inlet and outlet sections of the test section of the WT. We assume that inlet and outlet sections are perpendicular to the axis x (tunnel centreline). For the considered flow regimes of ETW, the flow is subsonic both at the inlet and at the outlet. From a purely

ARTICLE IN PRESS 104

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

mathematical viewpoint this means that from 5 flow parameters—Mach number M, total pressure Pt, total temperature Tt and turbulence parameters q, o—one should impose the values of 4 parameters at the inlet and of 1 parameter at the outlet. But from the practical viewpoint, a user should have the chance of obtaining the given flow regime, i.e. to impose all 5 parameters of the flow at the control point, which is placed at the centreline of the ETW test section side wall, not far from the test section inlet (x ¼ 646 mm). In fact, the user knows one ‘‘superfluous’’ parameter at the inlet and does not know anything about the outlet section. In practical operation, the control of the flow regime is accomplished by a variation of total pressure, total temperature and mass-flow-rate in the closed circuit of the ETW WT (that includes compressor and adaptive second throat). But in a calculation it is practically impossible to simulate this regulation directly (too complex geometry and too long calculations). As a whole, the boundary conditions at the inlet and outlet should provide the realization of the user-defined flow regime near the inlet section in minimal time. As a result of prolonged work, devoted to the choice and optimization of boundary conditions, the authors have selected the following combined approach: At the inlet section, all flow parameters are imposed— (M, Pt, Tt, q, o, v ¼ w ¼ 0). Formally (from the viewpoint of formulas and algorithm of calculations) this condition is equivalent to the ‘‘Riemann’’ boundary condition for the case, when gas flows into computational domain with supersonic speed. But in practice, as a rule, Mo1. It seems that such boundary condition contradicts subsonic character of flow (for subsonic flows, any information about a parameter is propagating upstream along the characteristic dx/dt ¼ uc). But, as one may see from the description of the ‘‘Riemann’’ boundary condition, the given parameters of outer flow are used only for a determination of ~ð1Þ . parameters from the outer side of the boundary P 1=2 After that, all parameters at the boundary are determined from the solution of the Riemann problem about the decay of discontinuity between flows from both sides of the boundary. The solution of the Riemann problem automatically uses part of the information from the inner nearboundary cell. As a result, the boundary condition remains mathematically correct. From the other point, one may hope that an imposition of the whole set of parameters from the outer side of the boundary should stimulate more effectively the convergence of solution to the required regime. At the outlet section, the pressure pexit is imposed. From a formal algorithmical viewpoint, this boundary condition almost coincides with the ‘‘Riemann’’ boundary condition for the case when gas flows out from a computational domain with subsonic speed; the only difference is that in the outer flow, adjoining the outlet boundary from the outer side, we impose pexit instead of the Riemann invariant z7. But the value of pexit is determined by a

special procedure. It is adjusted in a way that instant mass-flow rate through the outlet section of ETW coincides with the imposed mass-flow rate through the ETW inlet section. The outlet section is covered by a grid, consisting of tetragonal borders of the inner near-boundary cells. We shall numerate these borders by indices (j, k). Let us consider for some border (j, k) the system of equations for the parameters of flow, adjoining to the boundary from the ~ð1Þ : outer side P 1=2;j;k

8 pð1Þ > > 1=2;j;k ¼ pexit; > >     > ð1Þ > > 2c 2c > u þ ¼ u þ > g1 1=2;j;k g1 in ; > > > <  ð1Þ   p p ¼ g r rg in ; > 1=2;j;k > > > ð1Þ > > > v1=2;j;k ¼ vin ; wð1Þ > 1=2;j;k ¼ win ; > > > > : qð1Þ oð1Þ 1=2;j;k ¼ qin ; 1=2;j;k ¼ oin ;

(70)

where the index ‘‘in’’ corresponds to the flow parameters, adjoining the boundary from the inner side. The solution of this system will give us the values of the parameters ~ð1Þ P 1=2;j;k ¼ f ðpexit Þ. Except for the pressure, these parameters will be different in various points of the outlet section, because parameters in the inner near-boundary cells, generally speaking, are not constant across the section (especially during the process of convergence to stationary state). To determine pexit, we require the full mass-flow rate through the outlet section to coincide with the given massflow rate through the inlet section: o Xn ð1Þ r1=2;j;k ðpexit Þ uð1Þ 1=2;j;k ðpexit Þ DAj;k ¼ re ue Ae , j;k

where the index ‘‘e’’ corresponds to the given parameters at the inlet section, DAj, k is an area element of the outlet section, Ae is the area of the inlet section. It is obvious that (70) is an algebraic equation for the unknown value of pexit. This equation may be converted to the form a x  b xgþ1=2  1 ¼ 0, 1=g

where x ¼ pexit and the coefficients a, b are defined by a¼

1 X ðrin Þj;k ðuin þ 2cin =g  1Þj;k DAj;k , 1=g re ue Ae j;k ðp Þj;k in

X ðrin Þj;k ðcin Þj;k DAj;k 2 b¼ . gþ1=2g ðg  1Þre ue Ae j;k ðp Þj;k in

This transcendental equation has two roots: x1 and x2. The higher root corresponds to subsonic outflow through the outlet section. 5. Numerical modelling of the ETW test section The ETW sited in Cologne (Germany) is a pressurized cryogenic facility using nitrogen as test gas. Its ability to

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

independently vary tunnel pressure (120–450 kPa), temperature (310–110 K) and Mach number (0.15–1.35) allows performing pure Reynolds number, pure aeroelastic and pure Mach number investigations up to real flight conditions. The test section features dimensions 2.4  2  9 m3 (width  height  length). For full model testing floor and ceiling exhibit 6 slots each generating a total open area ratio of about 6%. Mass flow bypassing the model through the slots will return further downstream in the re-entry area. The test section itself is surrounded by a plenum of 10 m diameter. The classical support for transport aircraft models is the straight sting as shown on Fig. 21. The challenge for any numerical simulation of the real flow field is the existing geometrical complexity of slots, slats and the re-entry area (see Fig. 20). Obviously, the slot flow plays a major role at the development of the main flow despite the small dimensions of the slots themselves. Hence, a correct representation of the geometry is essential for the quality of he calculations. This will

Fig. 20. The re-entry area (looking upstream).

Fig. 21. Aircraft model on straight sting support in ETW.

105

strongly affect grid generation as the global test section dimensions as well as tiny details have to be covered simultaneously (Fig. 21). On the way to a numerical treatment of an aerodynamic problem the first step is always the geometrical representation of the boundaries. To limit the number of grid nodes in the final mesh it has to be carefully analysed which details have to be correctly represented and were simplifications are tolerable. Unfortunately, this decision can often not be made a priori and, hence has to be subject for specific studies. In the present case only printed design drawings were available to describe the tunnel hardware in detail. From that basis cross-sectional figures have been developed as shown in detail in Fig. 22 for a station x ¼ 3677 mm downstream of the inlet of the test section (i.e. the model area). The existing right–left and top–bottom symmetry revealed extremely helpful as one can then reduce the area to be considered to a quarter of the tunnel only. In the frame of the performed approach the lower left quarter has been selected. Referring to the remark above it should be emphasized that the slot width is 25 mm compared to the tunnel width of 2400 mm. Further on the slot channels have been aerodynamically optimized during the detailed design phase of the facility featuring a more complex geometry when expanding from their minimum width to about 90 mm inside the wall structure. In this area they contribute to the development of low-energy flow as we will see later. By nature, not all tiny details could be included in the final model. So some smoothing of the slot channel contour had to be applied (i.e. neglecting the slot radii; see close-up in Fig. 22). Generally, the test section was approximated as a rectangular box of 2000 mm height and 2400 mm width. As for standard tunnel operations regarding full models the top and bottom walls are set to a divergence of 0.551 each, the used coordinate system is referenced to the floor while the x-axis coincidences with the tunnel centreline. The resulting tunnel floor with its supporting structure is presented in Fig. 23. The flow mixing zone in the re-entry area is characterized by the saw-tooth or wedge-type structure (see also Fig. 20), profiled slots, ribs and other structural elements. Further specific elements are the reentry flaps located downstream of the test section. These flaps are going to be used to control the amount of flow bypassing the test section by propagating through the slot channels. The arrangement can nicely be identified in Fig. 20 including their hinged attachment. They can only be operated block-wise (3 flaps form a block) and their setting angles (nose down) have been optimized during the calibration phase of the facility. A sketched side view is given by Fig. 24. Attention is to be drawn to the front radius allowing a smoother return of the slot flow. The flaps exposed to the shaped re-entry side-walls feature an adapted contour but leave a small gap, which could not be modelled in the numerical representation. The re-entry area also houses the circular sector with the integrated sting-boss for an attachment of model supports.

ARTICLE IN PRESS 106

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 22. Sketch of the test section cross section at x ¼ 3677 mm (only the lower left quarter is shown).

Fig. 23. Numerical representation of the slotted floor of ETW.

Its downstream extension called ‘‘centre-body’’ divides the second throat duct in two symmetrical parts before they merge at the inlet of the high-speed diffuser. Using the original design drawings similar representations have been developed for the other tunnel areas like side-walls and plenum chamber. The latter one has a diameter of 10 m exhibiting a series of structural elements for the sake of accessibility and operation of tunnel components. It was concluded that there would be no flow in this area anymore and, hence, a neglect of these elements in the numerical model would be appropriate.

After meshing, the simplest numerical grid consists of 27 blocks. While one block represents the plenum chamber, the others are distributed over the test section, re-entry and slots. For example, each slot contains its appropriate block. One block is devoted to the mixing zone and another one covers the sting. Additional blocks are arranged upstream of the test section inlet and downstream of its exit. They are required to formulate the boundary conditions for this task. Fig. 25 presents an example of a numerical grid generated for a cross-section x ¼ const. To accelerate the calculation and to economize computer RAM, the numerical grid in different adjacent blocks is not joined (like a chimera condition). In principle, such an arrangement will result in accuracy losses but assuming the local gradients of flow parameters in the conjunction area to be small resulting accuracy losses will be small, too. With respect to the existing complexity of the problem to calculate the flow field around a test article in a slotted wall WT including walls the selected approach consisted of 3 major steps: (1) The viscous calculation in and around a single isolated slot channel. The relevance of this approach seemed acceptable as the existing hardware does not allow for any lateral flow exchange between the channels. (2) The calculation of the flow field around a centreline probe (short axial probe, SAP) attached to the sing

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

107

Fig. 24. Sketch of a re-entry flap in the mixing zone.

Fig. 26. Simplified model of a single slot.

Fig. 25. Conjunction of grid nodes at a slot–main flow interface.

boss. This probe features a body of revolution, hence a major simplification in its geometry. (3) The calculation of the flow field around a single straight sting supported wing-body configuration at transonic speeds. Available WT and model data could be used for comparative checks between CFD and experimental results for all selected configurations. 6. The flow behaviour within a single slot channel Preparatory numerical investigation revealed the flow near the ETW slots to be unsteady. This might be due to different effects such as a re-entry influence acting upstream, the presence of the plenum chamber and an unsteadiness of mixing layers. The selected approach simplifies the reality as given in Fig. 20 by considering the flow behaviour in a single isolated slot-channel as presented in Fig. 26. The geometrical model accurately represents the local test section floor with its 0.551 divergence, the slot channel beneath, the local re-entry area including the re-entry flap, the plate located beneath the slot channel (see Fig. 23) and part of the surrounding plenum. Preparatory calculations had investigated the depth of the plenum to be taken into account for a realistic representation.

The meshing is developed to describe a selected part of the geometry with maximal accuracy. It exhibits compressing the nodes near mixing layers and solid walls. Unfortunately, it is impossible to perform a perfect flow simulation near the re-entry, so the ‘‘non-slip’’ condition is formulated on re-entry flap. The calculation is performed in the frame of the RANS method described above. Regarding turbulence the performed calculations relied on the application of the (qo)-model. To be able to compare CFD and experiments the boundary conditions were taken from experiments: M ¼ 0.78, total pressure P0 ¼ 197 Kpa, total temperature T0 ¼ 300 K. Turbulence parameters at the inlet of the computational domain are equal to q ¼ 2.7 m/sec, o ¼ 200 Hz. With respect to the comments above it sounds fair to assume a lateral symmetry. So, this condition is formulated for all lateral vertical planes restricting the computational domain. The computational domain has been extended by a 5 m long buffer zone upstream of the physical domain to allow for a correct simulation of the boundary layer formed on the corresponding wall of ETW. Nevertheless, the boundary layer thickness at the entrance of the test section differs from its real value. Boundary conditions based on an analysis of the ‘‘Riemann invariants’’ are formulated at the inlet and exit of the computational domain. It should also be noted that the flow parameters at the inlet of the computational domain are not equal to the ones at WT Control Point (x ¼ 646 mm). For example, the Mach number is defined for this calculations as M ¼ 0.78, but in the experiment M ¼ 0.753 had been measured at the Control Point. This difference is produced

ARTICLE IN PRESS 108

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 27. Measured and calculated boundary layer profiles at different axial positions.

by a complicated process of stabilizing the unsteady process described as ‘‘starting the WT’’. This problem may be solved in the future by using the new algorithms that are not available now. The present status only allows a qualitative comparison of boundary layer profiles obtained by CFD and experiment. Considering Fig. 27 showing boundary layer profiles in sections x ¼ 2810 and x ¼ 4970 mm, the velocity is presented as non-dimensional values obtained by dividing locally measured values by the corresponding ones of the Control Point Vinf. It has to be noted that Vinf in the experiment was different from the one in CFD study (see above). Globally calculated and experimental profiles match fairly, nevertheless, essential differences exist. For x ¼ 2810 the ‘‘experimental’’ boundary layer is thinner than the ‘‘CFD’’ prediction (see Fig. 27). Further downstream at x ¼ 4970, the experimental profile reveals an s-shape characteristic for a separated boundary layer but CFD is in contradiction still indicating attached flow. The effect of slots in the walls of WTs is well known and has been subject of numerous experimental investigations all suffering the complex 3d flow structure combined with the restricted accessibility. As a sample, let us consider Mach number flow fields in different cross-sections corresponded to the beginning and end of a slot as well as in the re-entry zone (see Fig. 28). The analysis reveals the flow in the test section above the slots to be non-uniform. Different layers can be identified: The mixing layer shows its maximal intensity directly above the slot. Here it strongly penetrates the main flow in the WT. The slot flux intensity decreases when laterally moving away from the slot and the mixing layer starts merging the boundary layer (between slots) on top of the slat. It is interesting to note that the thickness of the mixing layer is comparable with the lateral distance between two slots (i.e. the width of

the computational domain indicated in Fig. 28) and that the boundary layer thickness is about 1/3 of this distance. This conclusion should be taken into account when trying to assess the ‘‘integral thickness of the boundary layer’’ on the surface of a slotted wall. The flow inside the slots channels is also characterized by high its nonuniformity. Fig. 29 presents velocity vectors documenting the existence of numerous vortices resulting in a strong generation of turbulence and, hence, providing an essential upstream influence inside the slot channel. Fig. 30 presents a summary of calculated flow details in a 3d reconstruction. It clearly illustrates the effect of the reentry zone and the flap on the slot flow development. This is due to the presence of strong pressure gradients caused by the channel expansion combined with a damping effect of re-entry flap. This process is unsteady by its nature and, hence, the given ‘‘flow reconstruction’’ represents a snapshot, while the real process is strongly time dependant. This fact may be better represented by a comparison of longitudinal and transversal velocity components along the slot at two distinct moments (Fig. 31). Note the two charts do not coincide, which may be explained by the existence of reverse flow inside the slots. Such reverse flow typically starts disturbing the mixing layer (at the boundary between test section and slot), hence, forming the main sources of unsteadiness, the primary subject under investigation. The performed investigations document the regular character of static pressure oscillations at different locations on the slotted surface. Fig. 32 presents the time-wise behaviour of the static pressure at different axial positions (x ¼ const.). For all considered cases, the graphs demonstrate a quasi-periodic process. It means that large-scale turbulence (including coherent structures) plays an essential part in forming flow near slot.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

109

Fig. 28. Mach number contours in cross-sections 1–3 (test section/slot-channel/plenum) (section 1–x ¼ 1400 mm, section 2–x ¼ 5680 mm, section 3–x ¼ 6950 mm).

7. The test section flow with the installed short axial probe (SAP) The SAP is a body of revolution holding to axial rows of pressure taps for a fast calibration of the Mach number and an assessment of the axial pressure gradient in the test section of ETW. A view of the SAP installed in the tunnel is given by Fig. 33. Its axial composition representing a straight model support can be taken from Fig. 34 (in millimetres). Being installed on the ETW model cart dedicated to full model testing its apex is located 1270 mm downstream of the test section inlet, while the point of model rotation is at x ¼ 3677 mm. The numerical representation in the frame of the considered CFD work has been done in accordance with Fig. 34 leading to the meshed arrangement presented in Fig. 35. It is shown attached to the model supporting sting-boss integrated in the circular sector. All slots and reentry area are taken into account. In this numerical model we were facing the same problem as already mentioned

above: the WT Control Point does not coincide with a grid node resulting in a non-perfect prescription of the correct WT operating condition. To overcome this problem the approach described in [3,4] has been applied because of the good results achieved during earlier calibration activities. A special boundary condition ‘‘Inlet-Exit’’ had been developed and permitted to calculate the flow around the SAP with acceptable accuracy. To understand this approach, we consider a simplified sample for flat channel with dimensions corresponding to the test section of ETW. The floor is also declined by 0.551 corresponding to the selected wall divergence in ETW to minimize the axial pressure gradient in the test section. Selected boundary conditions prescribe non-shock parameters in each section of the wind tunnel. Then, iterations are performed to adjust local parameters of the flow and to fulfil the ‘‘Inlet-Exit’’ mass flow law. It is assumed that all slots in the walls are closed and the Mach number M=0.95 at the inlet of the computational domain. Using the boundary condition ‘‘Inlet-Exit’’, one may obtain M=0.951 at Control Point.

ARTICLE IN PRESS 110

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 29. Flow structure in the cross-section x ¼ 6950 mm (re-entry).

This is considered a very promising result, as alternative boundary conditions may give errors up to 50% and more. In reality, this problem of ‘‘regime’’ is solved using more complex procedures including global iterations and specific calibrations. Along the SAP the flow field was analysed for different regimes with close and open wall slots. For each selected variant two additional cases have been calculated corresponding to the closed (no outflow from the slot channel into the re-entry) and open re-entry. The model was set to zero angle of attack. The used numerical grid was of medium quality and consisted of 100 cells along the model plus 10 for the representation of the sting boss. Some calculations have been performed using the approach based on Euler equations [3,4]. To take viscosity effects into account we used well-known procedure of geometry correction. A three-step procedure has been implemented: (1) inviscid flow calculation; (2) geometry correction; (3) inviscid calculation for corrected geometry. As a result, the effective ‘‘throat’’ of slot becomes smaller and, in some cases, even choking of a slot may occur, see Fig. 36. The solid line with symbols represents the original slot contour before any correction (x ¼ const.). The isolated symbols correspond to the corrected contour. When performing preliminary calculations, the boundary layer, obviously, has no pre-history and starts growing from the beginning of slots (x ¼ 250 mm). In reality the slot blocking may be more severe due to an essential influence of the boundary layers generated on the walls upstream of the test section. Fig. 37 presents two important cases:

A down view on a slot shape before and after correction. It demonstrates again that the slot boundary layer essentially changes the geometry of especially the slot width. At their downstream ends slots are practically ‘‘closed’’ and hence, the influence of the plenum chamber on the main flow diminishes. On the other hand the described correction procedure does not take into account many essential factors. For example, this approach will not work in the case of separated boundary layer, etc. The ‘‘only’’ benefit of this method is probably to lead to ‘‘fast’’ results and in many cases to achieve estimates useful for practical purposes. It is easy to calculate the effective open area of slots to correct the mass flow across. A more complicated sample is presented in Fig. 38. Here the axial distribution of the static pressure along the SAP surface is compared between CFD and experiment. All data are given as relative values by taring to the reference pressure, in our case the pressure in WT Control Point Pref. The calculation has been performed for closed and open slots in the floor and ceiling. Their geometry is corrected for the boundary layer thickness displacement according to the assessment described above. With respect to the complexity of the problem the numerical and experimental data are in a good agreement (do1% at xo6000). For closed slots, the discrepancy between CFD and experimental data (experimental data represent open slots) increases up to 10%. This clearly demonstrates the need to include the slots in a numerical simulation when striving for realistic and reliable data. An essential difference between CFD and experimental data (more than 10%) is to be noticed in re-entry area (x46000). Obviously, this is a result of the selected simplified numerical approach (Euler+boundary layer), which is not capable of covering the complex flow structures in that area. Later on, results obtained in this zone will have been excluded from an analysis. The presence of WT walls affects the flow around a test article and vice versa, known as ‘‘wall interference’’. To investigate this phenomenon, the reading from static pressure taps on the ETW walls has been made available (about 400 taps). They were placed along definite ‘‘lines’’ on the centre of the slots, i.e. between the slots. Each ‘‘line’’ is marked as shown in Figs. 39 and 40. One should pay attention to additional ‘‘lines’’ downstream of the test section named ROL 0, SBL 1, RBL 1, RBL 3 indicated in Fig. 40 and used for measurements in the re-entry zone and further downstream. The calculations were performed for different tunnel configurations. It was decided to start with closed slots but an open and closed re-entry (see Fig. 41). If the boundary layer is not taken into account, one may observe discrepancies between CFD and experimental data of about 5% (for Xo6000). Although it sounds strange it has been concluded that the flow development in the re-entry had only little effect on the quality of results obtained in the test section (compare solid and dash-dot lines in Fig. 41). The assessed difference was in the order of less than 1%.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 30. Reconstruction of the flow within and at the exit of a slot channel.

Fig. 31. Distributions of longitudinal and vertical velocities along slot.

111

ARTICLE IN PRESS 112

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 32. Time dependencies of the static pressure at some points on the lateral boundary of the computational domain (y ¼ 0.439 m, z ¼ 1.0 m).

Fig. 33. The short axial probe (SAP) installed in the test section of ETW.

Fig. 34. Sketch of the SAP.

In a second approach the flow in test section with open slots was investigated (see Fig. 42). It is evident that a consideration of the slots provides smoothing the pressure distribution along WT walls in comparison with the previous case (close slots). For example, open slots (X ¼ 6000) generate an increasing static pressure along ‘‘line’’ TBL 4 (compare solid and dash-dotted lines). The relative pressure comes closer to the measured value of P/Pref ¼ 1 (0oXo6000). For the case of open slots the

boundary layer seems not to influence the result above. So, it has been concluded that that the dominant part on flow development near the slots is taken by the shear layers and not by the boundary layers. The results given in Fig. 42 confirm this conclusion (compare solid and dotted lines). The differences received from a comparison of the pure EULER calculation and a consideration of the boundary layer is assessed negligible compared to the global difference gained for open and closed slots.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

113

Fig. 35. Meshed test section including SAP in ETW (one-quarter of test section only).

section of ETW, then extended to about 4 m. Three model supports have been studied. They had different geometries and inclination angles (01, 2.51, 51). Fig. 44 presents CFD and experimental data for the configurations with a nonzero angle. The numerical results for the pressure coefficient Cp achieved by a RANS calculation are in a good agreement with the experimental ones. The only region with visible differences is the most downstream part of the model in the area of 4100oxo4624. There, the experiments reveal a higher level of compression, which exceeds the ones predicted by CFD up to DCp ¼ 0.004. It has to be noted that the geometry in that area is extremely complex featuring a cavity inside the body, which has been accurately modelled in the simulation. Nevertheless, the numerical results demonstrate a level of accuracy achieved with the applied code which may be sufficient for practical purposes. Fig. 36. Applied correction methodology for the slot geometry to account for viscous effects.

Another sample documenting the benefits of an application of CFD in combination with WT investigations is presented by the configuration presented in Fig. 43. The test article features a body of revolution supported by different stings targeting for an experimental assessment of sting interference. This model holds 4 ‘‘lines’’ of pressure taps (top/bottom and right/left) to measure the static pressure. Data from the ‘‘central line’’ along the horizontal plane Y ¼ 0 are used in this article. The model axis coincides with the centreline of the tunnel. While the body is always kept in its position, different support systems have been attached; hence, the primary area of interest is the near field around the model rear. The nose of the body is located about 2 m downstream of the inlet of the test

8. Flow calculations on a wing-body configuration in the test section of ETW In the final paragraph of this article, flow calculations are presented for a wing-body configuration with a straight sting support inside the ETW test section. The sting diameter is 910 mm, while its tapered part starts at X ¼ 4861. The centre of model rotation in ETW is at X ¼ 3677 due to the radius of the sector, where the model support is attached to. Top/bottom wall divergence is equal to d ¼ 0.551, combined with a re-entry flap divergence of g ¼ 701. The investigated model as well as the main elements of the WT regarding the CFD simulation are given in Fig. 45. Two series of calculations have been performed: The first series reflects the ETW test conditions, the second one corresponds to free flight conditions. Hence, taring the results to the later is allowing the assessment of the WT effects on the aerodynamic

ARTICLE IN PRESS 114

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 37. A correction of the slot the geometry to consider viscous boundary layer effects in EULER calculations.

Fig. 39. Arrangement of pressure taps in a typical cross-section of the ETW test section.

Fig. 38. Comparison of CFD and experimental data for the static pressure distribution along the SAP surface. M ¼ 0.78, p0 ¼ 197 KPa, T0 ¼ 300 K.

performances of the considered model. The calculations have been performed both using methodologies [3,4] and the method described in this paper. The WT control system is linked to the reading of a total pressure and a total temperature probe located in the settling chamber, while the static pressure is picked up from a side wall tap named the Control Point at (X, Y, Z)ref=(646; 1000; 1200). Ideally, the customer wants to test his article at specifically defined conditions at the model reference point x=3677 on the centreline. Therefore, during the main calibration phase an experimentally assessed link between the conditions at the Control Point and the Model Reference Point (PMR) will be established (e.g. by using the SAP). The achieved correlation is going to be used for model testing, as the presence of a model does not allow for any reliable measurement at PMR anymore, by nature. While the conditions at the Control

Point are considered to be free of any upstream effect generated by the model the flow field around it is more or less disturbed by the presence of the surrounding walls , which is known as wall interference. To obtain the required values of flow parameters at the model location, several iterations with recalculations of the flow near the model and with a subsequent correction of the inlet parameters at the computational domain are necessary. Usually, two major parameters are given as initial ones: (1) Mach number Mref at PMR and (2) the model angle of attack a. The considered configuration underwent a series of tests in ETW over a wide range of parameters. The comparison of CFD and experimental data is based on the surface pressures gathered from the pressure-plotted model. The layout presented in Fig. 46 shows the arrangement of pressure taps organized in lines P1–P7 parallel to the centreline of the tunnel (and the model axis). Let us consider some comparisons of static pressures along those lines given in Fig. 47. The solid line corresponds to results obtained by the current RANS method, the dotted line shows formerly obtained results from the method of [3,4] and the symbols refer to experimental data. The main difference to be noticed is the disappearance of the shock

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

115

Fig. 40. Arrangement of ‘‘pressure lines’’ holding series of static pressure taps along the walls.

Fig. 41. Comparison of calculated and measured static pressure along the WT walls (closed slots). M ¼ 0.78, p0 ¼ 197 KPa, T0 ¼ 300 K.

on the upper side when going for the RANS solution, hence, achieving a much better match with the experiment. Obviously, this is mainly due to the consideration of viscous effects on the model surface by introducing the boundary layers. To get a feeling for the sensitivity of the solution regarding flow parameters and in the light of the relevant experimental results the flow field has been calculated for only slightly different free stream Mach numbers (DM=0.002) using the method of [3,4] (see Fig. 48). The resulting influence on shock position and intensity is clearly seen. A more severe change of the Mach number DM ¼ 0.01 (i.e. the accuracy of the regime definition in the

Control Point in this calculation) even shifts the shock by 10% of chord and reduces its intensity up to 20%. This may be an indicator for a peak design of the wing. The calculated effect of WT walls on the wing pressure distribution is impressively documented in Fig. 49. For the considered case the presence of WT walls leads to a shock displacement and a reduction of its intensity. A similar conclusion may already be drawn from the results applying method [3,4], which are presented in Fig. 50. It underlines the finding that often a qualitative result may be achieved, even regarding the flow in such a complex environment as a WT, by an application of simpler and faster methods.

ARTICLE IN PRESS 116

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 42. Comparison of calculated and measured static pressures along the WT walls (opened slots). M ¼ 0.78, p0 ¼ 197 KPa, T0 ¼ 300 K.

Fig. 43. Body of revolution attached to an ETW support system. Fig. 45. The wing-body configuration on its support as tested in ETW.

Fig. 44. Comparison of CFD and experimental data of static pressure distribution on the model surface. M ¼ 0.85, a ¼ 51, p0 ¼ 125 KPa, T0 ¼ 300 K.

Finally, an incremental comparison of calculated (RANS) aerodynamic coefficients is provided in Table 3. The condition M ¼ 0.85, a ¼ 1.56 (free flow) is used as reference. Note all given data are established dividing the actual by the corresponding reference values. As seen in Table 3 the aerodynamic coefficients differ strongly. But one can see some regularity. For example, influence of WT walls results in increasing CD coefficients, in the case of positive incidence angles. As already mentioned above one crucial problem regarding the match between CFD and experimental results is the correct representation of the Mach number at the model position, which is presently handled by a comparison of the conditions at the Control Point in the tunnel. Unfortunately, even a perfect agreement between CFD and test data at the Control Point will not automatically guarantee correct results around the model. The applied iterative method for the adjustment at this location leads to long calculation times. In the

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

117

Fig. 46. ‘‘Experimental pressure lines’’ on the model wings.

Fig. 47. Comparison of calculated and experimental data. M ¼ 0.85, a ¼ 1.56. (solid ¼ RANS, dashed ¼ method [3,4], symbols=experimental data).

9. Conclusions

Rodionov (GKR) for the description of convective fluxes, the explicit modified central-difference approximation of diffusive fluxes and the point-implicit approximation of source terms in equations for turbulence parameters. To overcome the problem of small time steps, the local time stepping is used for the solution of stationary tasks by a stabilization method and for the solution of unsteady tasks the fractional time stepping. Developing the numerical method resulted in some important findings:

A numerical method for the solution of time-averaged Navier-Stokes equations is proposed. It features a 2nd order approximation in all variables and includes an explicit monotone scheme of the Godunov–Kolgan–

(1) The estimation of the characteristic speed of the propagation of information due to diffusion is proposed. Applying this estimation, the generalized variant of the Courant–Friedrichs–Levi stability condition is

described cases, two to three loops to correct the parameters in the Control Point (by changing parameters at WT entrance) have been performed. As a result, the accuracy of the adjustment has been estimated as difference in Mach number DM ¼ 0.015. In some cases this might be sufficient for a practical implementation but generally it is still a subject for future modifications and improvements of the code (Table 4).

ARTICLE IN PRESS 118

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

Fig. 48. Effect of Mach number variations on the pressure distribution along line P1 for a=1.56 (calculations by method [3,4]).

Fig. 49. Calculated pressure distribution (RANS) for free flight and model inside the ETW M=0.85, a=1.56 (‘‘Line’’ P5).

proposed, which takes into account an interaction of convection and diffusion. (2) The criterion for a comparison of amplitude errors of numerical schemes is proposed. This criterion shows that the quality of an implicit scheme for the approximation of diffusive fluxes in the range of its applicability is much less than the quality of an explicit scheme in its working range.

(3) A new method for determination of gradients in GKR schemes is proposed, in which the limiter function is applied to the modulus of the velocity vector, but not to its direction. This method essentially improves the quality of the numerical solution. (4) A modification of the central-difference approximation of derivatives in diffusive fluxes is proposed.

ARTICLE IN PRESS S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

119

Acknowledgements The authors are grateful to the member of the Russian Academy of Sciences Prof. V. Ya. Neyland for taking the scientific leadership in the work described in this paper.

References

Fig. 50. Method [3,4] data for a configuration calculated in ETW and in free flow M=0.85, a=1.56 (‘‘Line’’ P5).

Table 4 Table regime

ETW

Free

M

a

CD

CL

CD

CL

0.85 0.7 0.89 0.85 0.85 0.85

1.56 1.56 1.56 3 0 2

0.9935 0.8608 1.4466 1.4466 0.7379 0.6893

1.0212 0.8167 1.1276 1.1276 0.5983 0.0747

1 0.8900 1.3398 1.3398 0.7573 0.7314

1 0.8249 1.1455 1.1455 0.5919 0.0899

It maintains the 2nd order approximation on strongly non-uniform grids and provides the interaction of cells with even and odd numbers. (5) It could be shown that in the approximation of source terms in equations for turbulence parameters the pointimplicit approximation is more reliable and sometimes more effective than the explicit approximation. A time step restriction for the point-implicit approximation is proposed on the basis of the analysis of exponentially varying modes of solution. (6) From various combinations of boundary conditions for the inlet and outlet sections of the WT test section, the best combination has been selected. It allows to obtain the given flow regime in the region of the inlet section. Flow parameters, which adjoin the entrance section from outside, are imposed inflexibly. At the outlet section, the pressure is imposed. The value of pressure is chosen to obtain the given mass-flow rate.

The proposed method is now available for a practical application but additional modifications may be required.

[1] Bosniakov SM. Approach to real aircraft 3D calculation based on time-dependent Euler equation system with boundary layer taking into account. Research Bulletin no. 6, Warsaw University of Technology, 1997. [2] Bosniakov S, Glazkov S, Matyash S, Neyland V, Remeev N, Yatckevich N. Numerical modelling flow-field around model in the big transonic wind tunnel using monotone 2nd order Godunov-type scheme. In: Wind tunnels and wind tunnel test techniques, Proceedings, Cambridge, UK, April 1997. p. 14–6. [3] Bosniakov S. Experience in integrating CFD to the technology of testing models in wind tunnels. Progr Aerospace Sci 1998;34:391–422. [4] Neyland V, Bosniakov S, Glazkov S, Ivanov A, Matyash S, Mikhailov S, et al. Conception of electronic wind tunnel and first results of its implementation. Progr Aerospace Sci 2001;37(2):121–45. [5] Jameson A. A perspective on computational algorithms for aerodynamic analysis and design. Progr Aerospace Sci 2001;37(2):197–243. [6] Loytzansky LG. Mechanics of fluids and gases. Moscow: Nauka; 1970 [in Russian]. [7] Piomelli U. Large-eddy simulation: achievements and challenges. Progr Aerospace Sci 1999;35:335–62. [8] Wilcox DC. Turbulence modeling for CFD. 2nd ed. DCW Industries; 1998. [9] Wilcox DC. Reassessment of the scale determining function for advanced turbulence models. AIAA J 1988;19(2):1299–310. [10] Jones WP, Launder BE. The prediction of laminarization with a twoequation model of turbulence. Int J Heat Mass Transfer 1972;15:301–14. [11] Launder BE, Sharma BI. Application of the energy dissipation model of the turbulence to the calculation of flow near a spinning disc. Lett Heat Mass Transfer 1974;1:131–8. [12] Lam CKG, Bremhorst KA. Modified form of the model predicting wall turbulence. J Fluids Eng 1981;103:456–60. [13] Menter FR. Improved two-equation turbulence models for aerodynamic flows. NASA TM-103975, 1992. [14] Menter FR, Langtry RB, Likki SR, Suzen YB, Huang PG, Vo¨lker S. A correlation based transition model using local variables Part 1, Part 2—model formulation ASME-GT2004-53452, ASME-GT2004-53454 ASME TURBO EXPO 2004, Vienna, Austria. [15] Rodi W, Scheuerer G. Scrutinizing the k–e turbulence model under adverse pressure gradient conditions. ASME J Fluids Eng 1986; 108:174–9. [16] Marvin JG. Turbulence modeling for hypersonic flows. The Third Joint Europe/US Short Course in Hypersonics, RWTH AachenUniversity of Technology, October 4, 1990. [17] Gerolymos GA, Sauret E, Vallet I. Effect of boundary-layer conditions in shock-wave/turbulent-boundary layer interaction computation. In: 33rd fluid dynamic conference, AIAA Paper 2003-3465, Orlando, FL, USA, 23–26 June 2003. [18] Gerolymos GA, Vallet I. Wall-normal-free near-wall Reynolds-stress closure for 3-D compressible separated flows. AIAA J 2001;39:1833–42. [19] Pope SB. An explanation of the turbulent round-jet/plane-jet anomaly. AIAA J 1978;16(3):279–81. [20] Krasheninnikov SJ, Mironov AK, Pavlyukov EV, Shenkin AV, Zhitenev VK. Mixer-ejector nozzles: acoustics and thrust characteristics. Int J Aeroacoust 2005;4(3–4):267–88. [21] Cambier L, Gazaix M. An efficient object-oriented solution to CFD complexity. AIAA Paper 2002-0108, 2002.

ARTICLE IN PRESS 120

S. Bosnyakov et al. / Progress in Aerospace Sciences 44 (2008) 67–120

[22] Boerstoel JW, Kassies A, Kok JC, Spekreijse SP. ENFLOW, a fullfunctionality system of CFD codes for industrial Euler/Navier– Stokes flow computations. In: Second international symposium on Aeron. Science and Techn., Jakarta, 1996 (also NLR-TP96286, 1996). [23] Eliasson P, Nordstr.om J, Torngren L, Tysell L, Karlsson A, Winzell B. Computations and measurements of unsteady pressure on a delta wing with an oscillating flap. In: Proceedings of the ECCOMAS, 1996 conference. [24] Hoffren J, Siikonen T, Laine S. Conservative multiblock Navier– Stokes solver for arbitrarily deforming geometries. J Aircraft 1995;32:1342–50. [25] Heinrich R, Kalitzin N. Numerical simulation of three dimensional flows using the chimera-technique. In: Proceedings of the 11th AG STAB/DGLR symposium, 1998. [26] Boniface J-Ch, Guillen Ph, Le Pape M-C, Darracq D, Beaumier Ph. Development of a chimera unsteady method for the numerical simulation of rotorcraft flowfields. AIAA Paper 98-0421, 1998. [27] Sillen M. Turbulent flow modeling using EARSM on parallel computers. TASK Q 2001;5(2):247–60. [28] Gacherieu C, Collercandy R, Larrieu P, Soumillon S, Tourette L, Viala S. Navier–Stokes calculations at aerospatiale matra airbus for aircraft design. In: Proceedings ICAS 2000 congress. [29] Kalitzin G, Gould ARB, Benton JJ. Application of two-equation turbulence models in aircraft design. AIAA Paper 96-0327, 1996. [30] Amato M, Paparone L. On the application of the multizone modeling and the local grid refinement technique for high-lift airfoil analysis. In: Proceedings of the ECCOMAS 1994 conference. [31] Bell HJ, Schairer ET, Hand LA, Mehta RD. Surface pressure measurements using luminescent coating. Annu Rev Fluid Mech 2001;33:155–206. [32] Vos J, Rizzi A, Darracq D, Hirschel E. Navier–Stokes solvers in European aircraft design. Progr Aerospace Sci 2002;38(38):601–97. [33] Marchuk GI. Methods of computational mathematics. Moscow: Nauka; 1980 [in Russian]. [34] Godunov SK, Zabrodin AV, Ivanov MYa, Kraiko AN, Prokopov GP. Numerical solution of multidimensional problems of gasdynamics. Moscow: Nauka; 1976 [in Russian]. [35] Lifshitz Yu.B, Sorokin AM, Shagaev AA, Plozky AI. Methods of calculation of flows around elements of aircrafts at transonic speeds. Part II: methods for grid construction. TsAGI Review no. 688, Russia, 1989 [in Russian]. [36] Jameson A. Transonic potential flow calculations in conservation form. In: Proceedings of the AIAA second computational fluid dynamic conference, Hartford, 1975. p. 148–61. [37] van Leer B. Towards the ultimate conservative difference scheme. Part V; a second order sequel to Godunov’s method. J Comput Phys 1979;32(1). [38] Eliasson P. Numerical validation of a half model high lift configuration in a wind tunnel. AIAA paper 2007-262, 2007. [39] Tysell L, Berling T, Eneroth P. Adaptive grid generation for 3D unstructured grids. In: Num. grid generation in computational fluid simulation, Proceeding of the sixth conference, London, England, July 6–9, 1998. [40] Melber-Wilkending S, Heidebrecht A, Wichmann G. A new approach in CFD supported wind tunnel testing. In: 25th international congress of the aeronautical sciences, ICAS 2006-3.4.2. [41] Kroll N, Rossow C, Schwamborn D, Becker K, Heller G. MEGAFLOW—a numerical flow simulation tool for transport aircraft design. ICAS 2002, Toronto, 1.5–10.5, 2002. [42] Heidebrecht A. A numeric far field model for support interference studies in a slotted wall wind tunnel (European Transonic Windtunnel). JAERO92, Special Issue Paper 581, 2006.

[43] Collercandy R, Marques B, Lory J, Dbjay S, Espiau L. Application of CFD for wall and sting effects. HiReTT report HIRETT WP2.231102003, Airbus France, October 2003. [44] Karlson K, Sedin J. Numerical design and analysis of optimal slot shapes for transonic test sections—axisymmetric flows. J Aircraft 1981;18(3):168–75. [45] Sedin J, Agrell N. Computation of wind tunnel flow in transonic slotted-wall test sections. ICAS-98-3, 11, 1–2. [46] Favre A. Equations des gaz turbulents compressibles I, II. J Me`canique 1965;3:361–90; 1965;4:391–421. [47] Coakley T. Turbulence modeling methods for the compressible Navier–Stokes equations. AIAA-83-1693, 1983. [48] Coakley T, Hsieh T. Comparison between implicit and hybrid methods for the calculation of steady and unsteady inlet flows. AIAA-85-1125, 1985. [49] Boussinesq J. Essai sur la theorie des eaux courantes. Paris Memories presentees par diverses savants a l’Acad d Sci, 1877. [50] Dash S, Weilersteen G, Vaglio-Laurin R. Compressibility effects in free turbulent shear flows, TR-75-1436, AFOSR, 1975. [51] Donaldson C, Rosenbaum H. Calculation of turbulent shear flows through closure of the Reynolds equations by invariant modelling. Aero Res. Assoc. of Princeton Report, 1968, p. 127. [52] Launder B, Reece G, Rodi W. Progress in the development of a Reynolds-stress turbulence closure. J Fluid Mech 1975;68:537–66. [53] Abramovich GN. Applied gasdynamics, vol. 1. Moscow: Nauka; 1991 [in Russian]. [54] Tikhonov AN, Samarsky AA. Equations of mathematical physics. Moscow: Nauka; 1972 [in Russian]. [55] Kogan MN. Dynamics of rarefied gas. . Moscow: Nauka; 1967 [in Russian]. [56] Khonkin AD. About paradox of infinite speed of perturbations propagation in gydrodynamics of viscid heat–conducting medium and in equations of hydrodynamics of fast processes. Article in the book Aeromechanics, devoted to 60th anniversary of Academician V.V. Struminsky. Moscow: Nauka, 1976 [in Russian]. [57] Fletcher K. Computational methods in fluid dynamics, vol. 1. Moscow: Mir; 1991 [in Russian]. [58] Ladyzhenskaya OA. Mathematical questions of dynamics of viscid incompressible fluid. Moscow: Nauka; 1970 [in Russian]. [59] Pervaiz MM, Baron JR. Spatiotemporal adaptation algorithm for two-dimensional reacting flows. AIAA J 1989;27(10). [60] Shokin Yu.I, Yanenko PN. Method of differential approximation: application to gasdynamics. Syberian Division, Nauka: Novosibirsk; 1985 [in Russian)]. [61] Ortega J, Poul W. Introduction in numerical methods for solution of differential equations. Moscow: Nauka; 1986 [in Russian]. [62] Kolgan VP. Application of the principle of minimal values of derivative to construction of finite-difference schemes for calculation of discontinuous solutions in gasdynamics. TsAGI Sci Notes 1972;3(6):68–77 [in Russian]. [63] Harten A. High resolution schemes for hyperbolic conservation laws. J Comput Phys 1983;49(3). [64] Kulikovskii AG, Pogorelov NV, Semenov AYU. Mathematical aspects of numerical solution of hyperbolic systems. Monographs and surveys on pure and applied mathematics, vol. 188. Boca Raton, FL: Chapman & Hall/CRC; 2001. [65] Roe PL. Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys 1981;43:357–72.