ERCOFTAC Special Interest Group on “Quality and Trust in Industrial CFD”
Best Practice Guidelines
Editors: Michael Casey and Torsten Wintergerste Fluid Dynamics Laboratory Sulzer Innotec
[email protected] [email protected]
Version 1: January 2000
© ERCOFTAC 2000
page 1 of 94
© ERCOFTAC 2000
page 2 of 94
CONTENTS 1. 1.1. 1.2. 1.3. 1.4. 1.5. 2. 2.1. 2.2. 2.3. 2.4. 2.5. 3. 3.1. 3.2. 3.3. 3.4. 3.5. 3.6. 3.7. 3.8. 4. 4.1. 4.2. 4.3. 4.4. 4.5. 4.6. 4.7. 5. 5.1. 5.2. 5.3. 5.4. 5.5. 6. 6.1. 6.2.
INTRODUCTION PURPOSES OF THIS DOCUMENT SCOPE BACKGROUND TO THE PROBLEMS OF CFD BACKGROUND TO THIS DOCUMENT STRUCTURE OF THIS DOCUMENT SOURCE OF ERRORS AND UNCERTAINTIES IN CFD SIMULATIONS OVERVIEW OF CFD SIMULATIONS IMPORTANCE OF ERRORS AND UNCERTAINTIES IN ENGINEERING APPLICATIONS ERRORS AND UNCERTAINTIES CATEGORISATION OF ERRORS AND UNCERTAINTIES DEFINITIONS OF VERIFICATION, VALIDATION AND CALIBRATION NUMERICAL ERRORS, CONVERGENCE AND ROUND-OFF ERRORS GLOBAL SOLUTION ALGORITHM CONVERGENCE ERRORS ROUND-OFF ERRORS GRIDS AND GRID DESIGN SPATIAL DISCRETISATION ERRORS TEMPORAL DISCRETISATION ERRORS FOR TIME-DEPENDENT SIMULATIONS TYPES OF IMPLEMENTATION OF BOUNDARY CONDITIONS SPECIAL REMARKS ON THE FINITE-ELEMENT METHOD TURBULENCE MODELLING RANS EQUATIONS AND TURBULENCE MODELS CLASSES OF TURBULENCE MODELS NEAR-WALL MODELLING WEAKNESSES OF THE STANDARD K-ε MODEL WITH WALL FUNCTIONS INFLOW BOUNDARY CONDITIONS UNSTEADY FLOWS LAMINAR AND TRANSITIONAL FLOWS APPLICATION UNCERTAINTIES GEOMETRICAL UNCERTAINTIES BOUNDARY CONDITIONS INITIAL CONDITION AND INITIAL GUESS UNCERTAINTIES WITH REGARD TO STEADY FLOW, SYMMETRY, AND PERIODICITY PHYSICAL PROPERTIES USER ERRORS GENERAL COMMENTS CONTROL OF THE WORKING PROCESS
© ERCOFTAC 2000
5 5 5 5 6 6 8 8 8 9 9 10 11 11 12 13 14 16 17 18 19 20 20 21 24 27 28 28 28 29 29 29 31 32 32 33 33 33 page 3 of 94
6.3. 7.
TRAINING REQUIREMENTS FOR CFD USERS CODE ERRORS
7.1. 7.2. 8.
37 37
VALIDATION AND SENSITIVITY TESTS OF CFD MODELS
39
SELECTION OF TEST CASES FOR VALIDATION DESIGN AND USE OF SENSITIVITY TESTS
39 40
EXAMPLES OF APPLICATION OF BEST PRACTICE GUIDELINES
9.1. 9.2. 9.3. 9.4. 9.5. 9.6. 9.7. 9.8. 9.9. 10.
41
INTRODUCTION TO THE TEST CASES TEST CASE A : 2-D TRANSIENT SCALAR BUBBLE CONVECTING AT 45 DEGREES TEST CASE B : T-JUNCTION BETWEEN MAIN AND AUXILIARY PIPE TEST CASE C : NATURAL CONVECTION FLOW IN A SQUARE CAVITY TEST CASE D : SUDDEN PIPE EXPANSION TEST CASE E : TRANSONIC AIRFOIL RAE2822 TEST CASE F : ENGINE VALVE TEST CASE G : LOW SPEED CENTRIFUGAL COMPRESSOR (LSCC) TEST CASE H : TURBULENT FLOW IN A MODEL OUTLET PLENUM
41 41 43 47 51 56 59 65 71
SUGGESTIONS AND NEED FOR FURTHER WORK
10.1. 10.2. 10.3. 11.
37
GUIDELINES FOR THE CODE DEVELOPER AND VENDOR GUIDELINES FOR THE CODE USER
8.1. 8.2. 9.
35
77
REVISION OF BEST PRACTICE GUIDELINES EXTENSION OF THE BEST PRACTICE GUIDELINES APPLICATION PROCEDURES
77 77 77
CHECKLIST OF BEST PRACTICE ADVICE FOR INDUSTRIAL CFD
11.1. 11.2. 11.3. 11.4. 11.5. 11.6. 11.7. 11.8. 11.9. 11.10. 11.11. 11.12. 11.13. 11.14.
78
GUIDELINES ON THE TRAINING OF CFD USERS GUIDELINES ON PROBLEM DEFINITION GUIDELINES ON SOLUTION STRATEGY GUIDELINES ON GLOBAL SOLUTION ALGORITHM GUIDELINES ON VALIDATION OF MODELS GUIDELINES ON TURBULENCE MODELLING GUIDELINES ON DEFINITION OF GEOMETRY GUIDELINES ON GRIDS AND GRID DESIGN GUIDELINES ON BOUNDARY CONDITIONS GUIDELINES ON THE SOLUTION OF DISCRETISED EQUATIONS USING A CFD CODE GUIDELINES ON ASSESSMENT OF ERRORS GUIDELINES ON INTERPRETATION GUIDELINES ON DOCUMENTATION GUIDELINES ON COMMUNICATION WITH CODE DEVELOPER
78 78 78 79 79 80 82 82 83 86 88 89 89 89
12.
ACKNOWLEDGEMENTS
91
13.
REFERENCES
92
© ERCOFTAC 2000
page 4 of 94
1. Introduction 1.1. Purposes of this document This document is intended as a practical guide giving best practice advice for achieving high-quality industrial Computational Fluid Dynamics (CFD) simulations using the Reynolds-averaged NavierStokes (RANS) equations. The advice given is primarily aimed at less experienced CFD users and is summarised as a checklist in chapter 11. The summary of best practice advice is also considered to be relevant for managers of scientific and engineering projects involving CFD simulations and for inspectors in regulatory bodies who have to examine submissions involving CFD simulations. It will also assist experienced CFD users, especially if they are moving to a new application area. The document also provides a useful compendium of relevant information on the most important issues relevant to the credibility of CFD simulations, especially with regard to the most common sources of errors and uncertainties in CFD. For each aspect considered, simple statements of advice are given which provide clear and generally accepted guidance for the user of CFD in industrial applications. The statements are understandable without elaborate mathematics, and are generally preceded by a short section providing a discussion of the issues involved, with references to texts where more detailed discussion of the issues may be found. A basic level of scientific or engineering knowledge is expected of the reader, and in particular some elementary knowledge of fluid dynamics and numerical methods is assumed (such as may be attained by reading one of the general texts given in the list of references in chapter 13).
1.2. Scope These guidelines cannot hope to be exhaustive. They have been written and edited making use of extensive consultation with CFD code vendors, code developers, academic experts and code users and this gives the guidelines wide support. It is intended that they offer roughly those 20% of the most important general rules of advice that cover roughly 80% of the problems likely to be encountered. The document is not part of a formal quality assurance (QA) management system but it nevertheless addresses most of the issues that a formal QA system for CFD simulations would need to include. Users who follow the advice given can be expected to avoid the most common pitfalls in CFD simulations. The technical content of the guidelines is limited to single-phase, compressible and incompressible, steady and unsteady, turbulent and laminar flow with and without heat transfer. The guidelines do not cover combustion, two-phase flows, flows with radiation, non-Newtonian flows, supersonic and hypersonic flows with strong shocks and many other more complex flow situations. The guidance given is relevant to many mechanical, aeronautical, automotive, power, environmental, medical and process engineering applications. More detailed guidelines and application procedures are needed for the remaining specific problems and to extend these guidelines to more complex flows (see chapter 10).
1.3. Background to the problems of CFD Over the last ten years CFD using the RANS equations has become a standard industrial simulation tool for the design, analysis, performance determination and investigation of engineering systems involving fluid flows. This development has been driven by the ready availability of robust in-house and commercial CFD software and by the massive increase in affordable computer speed and memory capacity, leading to a steady reduction in the costs of simulations compared to prototype and model experiments. Relatively few new technical developments in fluid flow modelling for industrial applications, in particular for turbulence modelling, have been made in this period. The fundamental problem of CFD simulations lies in the prediction of the effects of turbulence, which, using the words of Lamb [1895] of more than a hundred years ago, still remains “the main outstanding difficulty of our subject”. At a scientific level, turbulence is one of the great unsolved problems of nonlinear computational physics. To simulate turbulence directly by Direct Numerical Simulation (DNS) one needs to be able to capture the time and length scales of all the characteristic structures of the flow, from the energy-carrying large scales to the small dissipative scales. As these vary by several orders of magnitude, a typical spatial scale of 10-5 to 10-6 of the size of the computational domain in each co-ordinate direction has to be resolved. For engineering problems this is beyond the capacity of present or foreseeable computers, needing an increase of at least several orders of magnitude in computing power before it can become a general tool. Even on the largest possible computers, DNS © ERCOFTAC 2000
page 5 of 94
will remain a research tool for simple geometries at low Reynolds number for at least the next decade and probably longer. Large Eddy Simulation (LES), in which the large turbulent vortices are captured by the computational grid and the fine-scale turbulent motions are modelled by a so-called sub-grid model, is less intensive in computational resources and closer to engineering application. This is still unlikely to become a common engineering tool during the next ten years except, perhaps, for bluff body flows with low-frequency and large-scale oscillations (for example, external flows around buildings). In both academic and industrial circles there is little doubt that CFD on the basis of the RANS equations with a suitable turbulence model will form the basis of most engineering calculations for many years to come. Even when LES or DNS become practical, RANS solutions will still be used for the initial design explorations, just as simpler tools are used in design today, so that a way to deal with uncertainty arising from the turbulence modelling will still be needed. In addition to the physical difficulties of modelling the effects of turbulence, there are many other sources of error in CFD simulations (see chapter 2). A major issue is the accuracy of the numerical discretisation in CFD simulations on grids that are not good enough to produce grid independent solutions and this has led to extensive discussion of the credibility of many CFD simulations. The difficulties are amply demonstrated by the many CFD validation exercises involving blind test cases, where only sufficient information is made available to the participants to allow a CFD model to be set up and run, but the full test results are not available. The results of such exercises can be highly userdependent even when the same CFD software with the same models is being used. Recent examples are the ASME blind test case of the transonic compressor rotor NASA Rotor 37 (Strazisar and Denton [1995]), the EU-funded EMU Project on modelling of atmospheric dispersion near buildings (Hall and Cowan [1998]) and the ERCOFTAC workshop on draft tube flows (ERCOFTAC/IAHR [1999]). The greatest benefit arising from such validation exercises has often been to identify the many causes for differences in simulation results for what is nominally the same calculation, and to give insight into the range of results obtainable. One of the major causes of this problem is that for complex geometries in many industrial applications, the practical constraints on memory capacity and computer power lead to results that are not grid independent. The question naturally arises as to what are the best procedures to ensure that the results of a CFD simulation are accurate and credible. There are currently no standard best practice guidelines that can be used to achieve and to confirm that the best possible numerical accuracy has been achieved. There is an urgent need for such a code of best practice from the knowledge currently available to experienced industrial CFD users, code developers and academic experts to overcome this shortcoming.
1.4. Background to this document The newly-formed ERCOFTAC Special Interest Group (SIG) on “Quality and Trust in Industrial CFD” has identified that the production of best practice guidelines would be of great value to many of its members. With funding obtained from its members, it has commissioned the Fluid Dynamics Laboratory of Sulzer Innotec to organise the writing and editing of the ERCOFTAC Best Practice Guidelines. The work has been fully supported by CFD users in many industrial companies, by academic experts and by many CFD code developers and code vendors (see acknowledgements in chapter 12). Following a survey of users to identify their needs, an initial version of the Best Practice Guidelines was produced by the Fluid Dynamics Laboratory of Sulzer Innotec during May 1999. The initial version was then reviewed and criticised by an invited panel of CFD experts at a workshop held in Zurich, Switzerland on 21 and 22 June 1999. Following the workshop an interim version was produced during July 1999. This was circulated in August 1999 to participants of the workshop and to members of the ERCOFTAC SIG on “Quality and Trust in Industrial CFD” for further editing. Following the feedback on this, a draft of the final version was prepared for the ERCOFTAC SIG meeting on October 27 1999 in Florence, Italy. Discussion of this with academic reviewers led to changes in chapter 4 of the final version of January 2000.
1.5. Structure of this document The following structure has been adopted for the ERCOFTAC Best Practice Guidelines. Chapter 2 gives a general view of CFD, and provides a categorisation and description of the types of errors and
© ERCOFTAC 2000
page 6 of 94
uncertainties that can arise. This includes some definitions of terms, such as validation, verification, calibration, error and uncertainty. In the chapters 3 to 7 that follow, each of the different sources of error and uncertainty are examined and discussed and some guidelines on their effect and elimination are given. In view of the fact that the validation of models and the use of sensitivity tests are commonly accepted ways of dealing with some of the problems, an additional chapter on these issues has been added in chapter 8. The sections in the individual chapters overlap to a certain extent as similar issues often arise under different headings. Examples of the issues arising in the best practice guidelines as applied to some specific test cases are given in chapter 9. Suggestions for the revision and extension of the guidelines are given in chapter 10. A summary of the most important guidelines in the form of a checklist is provided as chapter 11, whereby this is an edited version of the guidelines in the earlier chapters organised roughly in the chronological order of the activities needed during a typical CFD simulation. Acknowledgements are given in chapter 12 and references in chapter 13.
© ERCOFTAC 2000
page 7 of 94
2. Source of errors and uncertainties in CFD simulations 2.1. Overview of CFD simulations CFD makes use of computer simulation to obtain an approximate solution to the governing equations of fluid flow. The solution is always approximate because only discretised versions of the continuum transport equations for fluid flow and energy transfer can be solved numerically. Moreover, in turbulent flows, the effects of turbulence cannot be represented in a mathematically accurate sense, but are modelled by approximate theories. Because it is not exact, CFD is used together with more traditional techniques, such as flow measurement or analytical methods, for the investigation of fluid flows. A typical CFD calculation can be broken down into a number of important steps and activities, such as: •
Training of CFD users.
•
Problem definition.
•
Selection of solution strategy.
•
Choice of numerical procedure.
•
Validation of models.
•
Selection of turbulence model.
•
Definition of geometry.
•
Generation of computational grid.
•
Definition of boundary conditions and physical properties.
•
Definition of initial guess or initial condition in unsteady flows.
•
Solution of numerical equations.
•
Assessment of errors and solution accuracy.
•
Post-processing and visualisation of simulation results.
•
Analysis and interpretation of simulation results.
•
Documentation and archiving of results.
•
Communications with the code developer.
Errors and uncertainty can arise in almost all of these steps and these are discussed in the sections below.
2.2. Importance of errors and uncertainties in engineering applications The requirement on the accuracy of a simulation depends to a large extent on the purpose of the analysis, and can vary from the need to know in which direction the flow is most likely to move, to accurate performance estimations for complex machines. Engineering flows are generally analysed with a view to identifying weak features of a design so that a component may be improved to provide economic gains from better product competitiveness and functionality. Even an inaccurate simulation can be of use in an engineering design, provided that the error bounds on the predicted parameter can be defined. For example, if the calculations show that the upper limit of a numerical quantity is far below an appropriate reference level of concern, then engineering decisions may be safely made. In the absence of quantitative accuracy the engineer often tends to examine the predicted flow-field in qualitative terms to assess global flow structures and trends (for example, incidence of flow onto blading, presence and extent of flow separations, existence of strong secondary flows, etc.). The details of the complex flow-fields provided by CFD, even when not perfectly accurate, are in this way an important source of insight into design improvements and also into the design of experiments and the most appropriate location of any instrumentation.
© ERCOFTAC 2000
page 8 of 94
Unfortunately, there are no standard ways of defining the error bounds on a CFD simulation. Nevertheless the CFD analyst should in the documentation of the results try to comment on the expected accuracy of the simulations.
2.3. Errors and uncertainties The deficiencies or inaccuracies of CFD simulations can be related to a wide variety of errors and uncertainties. The recent publication of the AIAA Guide for the Verification and Validation of Computational Fluid Dynamics Simulations, (AIAA [1998]), provides useful definitions of error and uncertainty in CFD along the following lines: Error:
A recognisable deficiency that is not due to lack of knowledge.
Uncertainty:
A potential deficiency that is due to lack of knowledge.
These rather philosophical definitions can be made clearer by examples. Typical known errors are the round-off errors in a digital computer and the convergence error in an iterative numerical scheme. In these cases, the CFD analyst has a reasonable chance of estimating the likely magnitude of the error. Unacknowledged errors include mistakes and blunders, either in the input data or in the implementation of the code itself, and there are no methods to estimate their magnitude. Uncertainties arise because of incomplete knowledge of a physical characteristic, such as the turbulence structure at the inlet to a flow domain or because there is uncertainty in the validity of a particular flow model being used. An error is something that can be removed with appropriate care, effort and resources, whereas an uncertainty cannot be removed as it is rooted in lack of knowledge.
2.4. Categorisation of errors and uncertainties The structure of this guideline makes use of a simple categorisation of errors and uncertainties given in a lecture by Ferreira and Scheuerer [1997] at the first European meeting of the ERCOFTAC SIG on “Quality and Trust in Industrial CFD” in Paris. Ferreira and Scheuerer suggested some clear and logical distinctions between types of errors occurring in a CFD computation, which with slight modification for the present purposes, are explained in the sections below and form the basis for the main chapters of this document.
2.4.1. Model uncertainty These are the uncertainties due to the difference between the real flow and the exact solution of the model equations. This includes errors due to the fact that the exact governing flow equations are not solved but are replaced with a simplified model of reality. The most well-publicised errors in this category are the errors from turbulence modelling, but other model errors may occur. Examples would be the simplification of an equation of state of a real gas to that of an ideal gas, the assumption of incompressible flow when compressibility effects and strong heat transfer occur, the neglect of nonNewtonian viscous effects, or the simplification of a complex combustion process to a few simple equations. In short, the model errors and uncertainties can be described as the uncertainties which arise because we are in fact solving the wrong equations.
2.4.2. Discretisation or numerical error These errors arise due to the difference between the exact solution of the modelled equations and a numerical solution with a limited resolution in time and space. For consistent discretisation schemes, the greater the number of grid cells, the closer the results will be to the exact solution of the modelled equations, but both the fineness and the distribution of the grid points affect the result. This type of error arises in all numerical methods and is related to the approximation of a continually varying parameter in space by some polynomial function for the variation across a grid cell. In first order schemes, for example, the parameter is taken as constant across a certain region. In short, discretisation errors arise because we do not find an exact solution to the equations we are trying to solve, but numerical approximations to them.
2.4.3. Iteration or convergence error These errors occur due to the difference between a fully converged solution on a finite number of grid points and a solution that is not fully converged. The equations solved by CFD methods are usually iterative, and starting from an initial approximation to the flow solution, iterate to a final result. This
© ERCOFTAC 2000
page 9 of 94
should ideally satisfy the imposed boundary conditions and the equations in each grid cell and globally over the whole domain, but if the iterative process is incomplete then errors arise. In short, convergence errors arise because we are impatient or short of time or the numerical methods are inadequate and do not allow the solution algorithm to complete its progress to the final converged solution.
2.4.4. Round-off errors Round-off errors are due to the fact that the difference between two values of a parameter is below the machine accuracy of the computer. This is caused by the limited number of computer digits available for storage of a given physical value.
2.4.5. Application uncertainties Inaccuracy is also introduced because the application is complex and precise data needed for the simulation is not always available. Examples of this are uncertainties in the precise geometry, uncertain data or models that need to be specified as boundary conditions (such as turbulence properties at an inlet or a tabulated equation of state) and uncertainties as to whether the flow is likely to be steady or unsteady.
2.4.6. User errors Many errors also arise from mistakes and carelessness of the user. Such errors generally decrease with increasing experience of the user, but in the nature of things cannot be completely eliminated as “to err is human”. This error is often described by the popular jibe “garbage in, garbage out”. Sometimes, however, garbage in does not lead to (terrible) garbage out because a genuine user error does not always have a significant effect on the results.
2.4.7. Code errors Errors also occur due to bugs in the software, unintended programming errors in the implementation of models or compiler errors on the computer hardware being used. Such errors are often difficult to find, as CFD software is highly complex, typically involving hundreds of thousands of lines of code for a commercial product. Computers are very unforgiving. Even a relatively simple typing error that might easily be overlooked on this page, such as an “i” for a “j” in a single word, can have disastrous consequences when incorporated into a line of code.
2.5. Definitions of verification, validation and calibration In discussions of CFD errors and uncertainties it is useful to make some clear distinctions between the meaning of the terms validation, verification and calibration. The definitions used in these guidelines follow closely the similar definitions given in the AIAA guide [1998], Roache [1998], Rizzi and Vos [1998] and Fisher and Rhodes [1996]: Verification:
Procedure to ensure that the program solves the equations correctly.
Validation:
Procedure to test the extent to which the model accurately represents reality.
Calibration:
Procedure to assess the ability of a CFD code to predict global quantities of interest for specific geometries of engineering design interest.
Calibration is sometimes also used to describe the process of adjusting the values of the coefficients of a turbulence model to provide better agreement with experimental data, but this is not its meaning in this document. In many industrial engineering cases, the distinction between validation and calibration becomes blurred.
© ERCOFTAC 2000
page 10 of 94
3. Numerical errors, convergence and round-off errors The physical problem being solved by an analysis of the fluid flow can be mathematically described by the equations for conservation of mass, momentum and energy. These partial differential equations (PDEs) are based on the assumption that the fluid can be described as a continuous medium. The techniques used to solve the problem replace the PDEs by a set of algebraic equations by breaking down the physical domain into a large number of discrete control volumes, called elements or cells. Within these cells algebraic relationships describe how the flow variables, such as velocity, temperature or pressure, vary locally with the space co-ordinates. For instance, a quadratic variation across a cell results in a formally third order scheme, linear variation in a second order scheme and constant behaviour in a first order scheme, at least on regular Cartesian grids. The order indicates up to which term of a Taylor expansion series the discretised system of equations is consistent with the original set of PDEs. How the basic volume elements are defined mathematically is particular to the numerical method chosen. Examples of such methods are the finite-difference, finite-volume or finite-element methods, each of which has its own specific advantages and disadvantages. The discretisation in time follows a similar scheme to that in space but respects the fact that time exhibits a clearly distinct forward direction. Discretisation is probably the most crucial source of error for the accuracy of numerical fluid flow simulations and needs to be carefully analysed by the user of CFD-codes. However, other sources of error may be related to the method used to calculate the fluxes at the cell boundaries and to solve the resulting system of linear equations. As the numerical procedures are iterative the question of convergence needs additional attention.
3.1. Global solution algorithm The discretised set of equations can be obtained with various solution procedures which include both pressure-based and density-based methods (for a review, see Ferziger and Peric [1972], Fletcher [1991], Patankar [1980] or Hirsch [1991]). The solution algorithms make use of numerous tuning parameters, such as artificial time-steps, under-relaxation, etc., to improve the convergence behaviour and robustness of the code. The field of application of a code and the modelling technique included influence the choice of the numerical method and the solution procedure. For example, computation of highly compressible supersonic and hypersonic flow fields often requires the use of a different solution technique than that used for an incompressible flow simulation. In principle, the solution of a well converged simulation is independent of the numerical method and the solution algorithm chosen, provided the grid is fine enough.
3.1.1. Guidelines on global solution algorithm Check the adequacy of the solution procedure with respect to the physical properties of the flow by reference to the code user manual, for instance compressible or incompressible flow. ¾ As a first step in the simulation process, use the parameters that control convergence of the solution algorithm (e.g. relaxation parameters, damping factors or time steps) as suggested by the CFD-code vendor or developer. ¾ If it is necessary to change these parameters to aid convergence, do not change too many parameters in one step, as it then becomes difficult to analyse which of the changes have influenced the convergence. In case of persistent divergence check the advice given in the sections on boundary conditions, grid discretisation and convergence errors, before changing the solution algorithm. ¾ If a steady solution has been computed and there is reason to be unsure that the flow really is steady, then carry out an unsteady simulation with the available steady flow field as the initial condition, preferably with a small perturbation imposed. Examination of the time-development of the physical quantities in the locations of interest will identify whether the flow is steady or not.
¾
© ERCOFTAC 2000
page 11 of 94
3.2. Convergence errors Iterative algorithms are used both for steady state solution methods and for procedures to obtain an accurate intermediate solution at a given time step in transient methods. Progressively better estimates of the solution are generated as the iteration count proceeds. There is no convergence theory for the solution of the discrete RANS equations. Hence convergence cannot be enforced by theoretical means, but is based on empirical criteria. Driving all the residuals in all equations plus the residuals of the integral balances down to machine accuracy is the most generally accepted convergence criterion. This level of convergence is an ideal that, for practical reasons related to time constraints, is not used for most engineering computations. Instead the level of convergence is most commonly evaluated based on average values of the residuals as described below. However, some methods measure the normalised change in variable values between successive iterations, or on values of globally integrated parameters, such as lift coefficient or heat transfer coefficient, or on time/iteration signals of a physical quantity at a monitor point, which is located at a suitably chosen location in the flow domain. In some situations the iterative procedure will not converge, but either diverges or remains at a fixed and unacceptable level of error, or oscillates between alternative solutions. Careful selection and optimisation of control parameters (such as damping factors, relaxation factors, or time-steps) may be needed in these cases to ensure that a converged solution can be found wherever possible. A procedure which proved successful for one class of problems may not be successful for a different problem.
3.2.1. Residuals Residuals are 3D fields associated with a conservation law, such as conservation of mass or momentum. They indicate how far the present approximate solution is away from exact cancellation of flux balances in each cell. Usually, the residuals are normalised by dividing by a reference value times a reference flux. The reference value may be one of the following: •
Maximum value of the related conserved quantity.
•
Average value of the related conserved quantity.
•
Range (i.e. maximum minus minimum) value of related conserved quantity.
•
Inlet flow of the related conserved quantity.
The reference flux may be: •
The local central finite difference of finite volume coefficient
•
Inlet mass flow rate of a related quantity
Convergence is usually monitored on the basis of one representative number, a so-called norm (usually a p-norm) characterising the residual level in the three-dimensional flow field. This single value may be: •
The maximum of the absolute values (p=∞).
•
The sum of absolute values (p=1)
•
The square root of the sum of squared values (P=2).
The large number of variants makes it difficult to give precise statements how to judge convergence and at which residual level a solution may be considered converged. In principle, a solution is converged if the level of round-off error is reached. Special care is needed in defining equivalent levels of convergence if different codes are used for comparison purpose.
3.2.2. Guidelines on convergence Carefully define solution-sensitive target quantities for the integrated global parameters of interest and select an acceptable level of convergence based on the rate of change of these target parameters (such as heat load, mass flow, lift, drag, and moment forces on a body). Convergence of a simulation should not be assessed purely in terms of the achievement of a particular level of residual error. ¾ For each class of problem you are interested in, carry out a test of the effect of converging to different levels of residual on the integrated parameter of interest (this can be a single calculation
¾
© ERCOFTAC 2000
page 12 of 94
¾ ¾
¾
¾ ¾
¾
that is stopped and restarted at different residual levels). This test demonstrates at what level of residual the parameter of interest can be considered to have converged and identifies the level of residual that should be aimed at in similar simulations of this class of problem. Be aware that different codes have different definitions of residuals. Always check the convergence on global balances (conservation of mass, momentum and energy) where possible, such as the mass flow balance at inlets and outlets and at intermediate planes within the flow domain. Check the maxima and minima of all variables which may be limited by the code (for example k and ε which should always be positive), and the number of cells where any limitation on values is active. Check not only the residual itself but also the rate of change of the residual with increasing iteration count. Monitor the solution in at least one selected point in a sensitive area to see if the region has reached convergence, but be aware that this alone is insufficient as a criterion for convergence as the solution may be "stuck" far from convergence but appear converged if solution variables are unchanging. For calculations that are proving difficult to converge, the following advice may be helpful: • • • • • • • • • •
¾
Use more robust numerical schemes during the first (transient) period of convergence and switch to more accurate numerical schemes as the convergence improves. Reduce parameters controlling convergence, such as the under-relaxation parameters or the CFL1 number. If the solution is heavily under-relaxed, increase the relaxation factors at the end to see if the solution holds. Check whether the grid quality in areas with large residual has any effect on the convergence rate. Look at the residual distribution and associated flow field for possible hints, e.g. regions with large residuals or unrealistic velocity levels. When the residuals stagnate at too high a value, try to further converge the problem with a double precision version of the flow solver. Check whether switching from a steady-state to a time-accurate calculation has any effect. Consider using a different initial condition for the calculation. Check the numerical and physical suitability of the boundary conditions. Examine the choice of variable on which the residual is based, for example density may be a perfectly reasonable choice for transonic flow calculations, but is certainly not for low speed subsonic flows.
Guidelines and recommendations to the code developers: • • • •
CFD codes should make available the maximum possible information to judge convergence. This includes residuals for every conserved quantity. The residuals should be dimensionless and the definitions of the residuals should be clearly defined in the handbook. The codes should supply information on the spatial distribution of residuals. To avoid the confusion of CFD users, CFD code developers should try to adopt one commonly accepted standard definition of the residual.
3.3. Round-off errors In situations where the small arithmetical differences between two large numbers become relevant, cancellation due to round-off may lead to severe errors. To avoid a negative influence of round-off errors, it is common practice to calculate pressure relative to a reference value. Some codes also use a residual form of the governing equations to combat round-off errors. Examples where round-off errors are known to be of significance are: •
Low Reynolds number turbulence models with large exponential terms.
1
The CFL number is generally defined as CFL=∆t (v+a) / ∆x, where ∆t is the time step, ∆x the local cell size, v the local velocity and a the sonic velocity. For incompressible flows computed with an incompressible solver, the CFL number is defined by CFL=∆t v / ∆x. © ERCOFTAC 2000
page 13 of 94
• • • • • •
Flows with density driven buoyant forces with small density and temperature differences. High aspect ratio grids with large area ratios on different sides of the grid. Conjugate heat transfer. Calculations of scalar diffusion with low concentrations of one species. Low Mach number flows with a density based solver. Flows with large hydrostatic pressure gradients.
3.3.1. Guidelines on round off errors Check the influence of round-off errors by performing the same calculation with a single precision and a double precision version of the flow solver, and compare target parameters. ¾ Always use the 64-bit representation of real numbers (double precision on common UNIX workstations), where possible. ¾ Developers are recommended to use the 64-bit representation of real numbers (REAL*8 in FORTRAN) as the default settings for their CFD code.
¾
3.4. Grids and grid design The computational grid represents the geometry of the region of interest. The computational domain needs to be discretised using grid cells that should provide an adequate resolution of the geometrical and expected flow features. Two approaches exist to capture the geometrical details of the domain: •
A regular (and structured) grid is created that is large enough to incorporate the complete geometry. Regions extending over the boundary are then cut away. This may result in two different ways of representing curved domain surfaces: • By a sequence of steps, whereby the elements close to the surface keep their regular shape. • By a sequence of flat facets, of which at least the two ends lie on the geometrical surface. The boundary elements then have an irregular shape.
•
Body-fitted grids where the cell surfaces follow a curved domain surface as a sequence of flat facets. The cells are designed to keep a regular shape. In structured meshes the grid lines follow the surface. This type is most frequently used.
In the case of body-fitted grids the interior of the domain must be built up to satisfy the geometrical constraints imposed by the domain boundary. Several kinds of mesh topology are available: •
Structured grid: The points of a block are addressed by a triple of indices (ijk). The connectivity is straight-forward because cells adjacent to a given face are identified by the indices. Cell edges form continuous mesh lines which start and end on opposite block faces. Cells have a hexahedral shape.
•
Unstructured grid: Meshes are allowed to be assembled cell by cell freely without considering continuity of mesh lines. Hence, the connectivity information for each cell face needs to be stored in a table. The most typical cell shape is the tetrahedron, but any other form including hexahedral cells is possible. These grids may or may not have matching cell-faces. Special cases of unstructured grids are: • Block structured grid: For the sake of flexibility the mesh is assembled from a number of structured blocks attached to each other. Attachments may be regular, i.e. cell faces of adjacent blocks match, or arbitrary (general attachment without matching cell faces). • Chimera grid: Structured mesh blocks are placed freely in the domain to fit the geometrical boundaries and to satisfy resolution requirements. Blocks may overlap, and instead of attachments at block boundaries information between different blocks is transferred in the overlapping region. • Hybrid grid: These grids combine different element types, i.e tetrahedra, hexahedra, prisms and pyramids.
The grids must be fine enough to provide an adequate resolution of the important flow features, as well as geometrical features. This may be achieved by local grid refinement. Unstructured meshes are especially well suited for this purpose. If block structured grids are used local refinement results in block attachments with dissimilar number of grid lines. Some CFD codes provide algorithms to adapt the grid resolution locally according to numerical criteria from the flow solution, such as gradient information or other error estimators.
© ERCOFTAC 2000
page 14 of 94
The accuracy of the simulation usually increases with increasing number of cells, i.e. with decreasing cell size. However, due to limitations imposed by the increased computer storage and run-time some compromise in mesh size is nearly always inevitable. In addition to grid density, the quality of a mesh depends on various criteria such as the shape of the cells (aspect ratio, skewness, warp angle or included angle of adjacent faces), distances of cell faces from boundaries or spatial distribution of cell sizes. The introduction of special topological features such as O-grids or C-grids and care taken to locate block-interfaces in a sensible manner can help to improve the overall quality of a block-structured mesh. Unstructured meshing techniques may take advantage of prism layers with structured sub-meshes close to domain boundaries.
3.4.1. Guidelines on grids and grid design ¾
¾
¾
¾
¾
¾
¾
¾
¾
¾
¾ ¾
Choose a suitable global topology of the mesh (for example, a scalable grid topology where the grid can be refined by parametric changes without degrading the grid quality) to help satisfy the specific code's requirements with regard to skewness, aspect ratio and expansion ratios, as outlined below. Ensure that the extent of the computational domain has been chosen to capture relevant flow and geometrical features, and if necessary examine the sensitivity of the calculation to the choice of computational domain. Assess which geometrical features can be omitted, for example, those that have dimensions below that of the local grid size. In areas where local detail is needed, then consider the use of local grid refinement to capture fine geometrical details. If grid refinement is used the additional grid points should lie on the original geometry and not simply be a linear interpolation of more grid points on the original coarse grid. The use of distributed losses or porosity might also be considered to take into account highly obstructed zones or fine detail of some obstructions (for example, porous filters and packings), whereby the code manual should be consulted for the interpretation of the turbulent fields obtained. Avoid highly skewed cells. For hexahedral cells or prisms the included angles between the grid lines should be optimised in such a way that the angles are approximately 90 degrees. Angles with less than 40 or more than 140 degrees often show a deterioration in the results or lead to numerical instabilities. Tetrahedra should tend to have their four angles equal. Avoid highly warped cells, that is cells with large deviations from co-planar faces. Warp angles (measured between the surface normals of triangular parts of the faces) greater than 75 degrees can lead to serious convergence problems and deterioration in the results. Avoid non-orthogonal cells near boundaries. The angle between the grid lines and the boundary of the computational domain (the wall or the inlet and outlet boundaries) should be close to 90 degrees. This requirement is stronger than the requirement given for the grid angles in the flow field far away from the domain boundaries. Avoid the use of tetrahedral elements in boundary layers. Prismatic or hexahedral cells are to be preferred because of their regular shape and their ability to adjust in accordance with the near-wall turbulence model requirements. Avoid aspect ratios that are too high. Away from boundaries, ensure that the aspect ratio (the ratio of the sides of the elements) is not too large. This aspect ratio should be typically not larger than 20 to 100 (depending on the flow solver) in important regions of the flow domain but may be larger in non-critical regions, whereby precise values should be supplied in the code manual. The accuracy and possible convergence difficulties in such cases depends greatly on the flow direction. Near walls this restriction may be relaxed and indeed it can be beneficial to have high aspect ratio in the boundary layers. Observe any specific requirements on mesh expansion ratios. The specific code's requirements for cell mesh stretching or expansion ratios (rates of change of cell size for adjacent cells) should be observed. The change in mesh spacing should be continuous and mesh size discontinuities be avoided, particularly in regions of large changes. Note that some codes offer automatic grid adaptation techniques, but that depending on the software being used these might not improve the grid quality (skewness, aspect ratio). Use a finer and more regular mesh in critical regions with high flow gradients or with large changes, such as shocks, regions with high shear, regions with significant changes in geometry or
© ERCOFTAC 2000
page 15 of 94
¾ ¾ ¾
¾
where suggested by error estimators. Make use of local refinement of the mesh in these regions, in accordance with the selected turbulence wall modelling. Ensure high geometric precision of the periodic grid interface when using periodic boundary conditions. Avoid the presence of arbitrary mesh coupling, non-matching cell faces, grid refinement interfaces or extended changes of element types in the critical regions of high flow gradients. Check that the assumptions made when setting up the grid with regard to critical regions of high flow gradients and large changes agree with the result of the computation, and rearrange grid points if found to be necessary. For each class of problem, make use of a grid dependency study to analyse the suitability of the mesh and to give an estimate of the numerical error in the simulation. Ideally, at least three significantly different grid resolutions should be used. Strictly, one should double the grid twice in each direction and then apply Richardson extrapolation (Roache [1997]) to be really sure. If this is not feasible, apply selective local refinement of the grid in critical flow regions of the domain to allow greater factors, or try to compare different order of spatial discretisations on the same mesh.
3.5. Spatial discretisation errors Spatial discretisation errors are primarily concerned with the approximation of convective terms in the governing transport equations. Different numerical methods evaluate the fluxes at either the same grid locations as the transported quantities or somewhere in between (co-located or staggered grids). In both cases, an algebraic approximation of the spatial functions is required to calculate the gradients at these locations. This approximation is called the differencing scheme in finite volume or difference methods or the basis function in finite element methods. The accuracy of the scheme depends on the form of the algebraic relationship and on the number and location of participating grid points (sometimes known as the stencil). The spatial discretisation or truncation error equals the difference between the scheme and the exact formulation based on a Taylor expansion series. A formally second order scheme is consistent with the exact formulation up to the term with a power of two, a third order scheme also takes into account the next higher term. For some discretisation schemes, the formal order of accuracy is not preserved on irregular meshes, where it reduces by one, or even more on cells with discontinuous sizes. Reducing the cell size by introducing a finer grid has the biggest impact on the accuracy of the solution if higher order schemes are applied. Halving the elements in all directions using a third order scheme will reduce the numerical error by a factor of 8, while this factor is only 2 with a first order scheme. It should, however, be borne in mind that the order of a method has no direct connection to the accuracy of a solution on a given grid. For instance lower order upwind methods can yield better solutions on coarse grids than higher order schemes. However, the higher order central schemes will approach the exact solution faster when the grid is refined. If the solution of the physical problem considered is smooth, and exhibits only small gradients, even a first order scheme can do a good job, especially if the flow is oriented along the grid lines. A first order scheme is usually not at all suitable for general engineering applications involving complex flows with large gradients and thin boundary layers. The large truncation error introduced by the first order scheme is known as numerical viscosity or diffusivity as it gives rise to artificial diffusion fluxes, which may be much stronger than the real molecular or turbulent contributions. On the other hand, higher order schemes suffer from a different more obvious problem, namely the appearance of a characteristic wavy pattern with a wavelength of two cell sizes in the neighbourhood of steep gradients. These so-called wiggles are caused by dispersion, i.e. waves with different wave lengths are not transported with the same speed. Dispersion is most prominent in central differencing schemes for finite volume methods and quadratic basis function schemes for finite element methods. Higher order upwind schemes are less prone to it. In some cases these wiggles can be considered a benefit, as their size is an indicator of the mesh error. If necessary, this problem may be remedied using special (non-linear) monotone schemes, such as total variation diminishing (TVD) or shockcapturing schemes (Hirsch [1991]). Due to their capability to resolve steep gradients while avoiding dispersion effects they are frequently applied in supersonic flows with shock waves, and for the transport of scalar quantities with weak molecular diffusion.
© ERCOFTAC 2000
page 16 of 94
3.5.1. Guidelines on spatial discretisation Avoid the use of first order schemes. The use of methods of higher order (at least second) is recommended for all transported quantities. ¾ First order schemes may be acceptable under the following circumstances. It may be necessary to use a first order scheme at the start of a calculation as it is likely to be more robust, but as convergence is approached a second order or higher scheme should be used. From a practical point of view it will sometimes be acceptable to use a first order scheme for the transported turbulence variables (which are often dominated by local production and dissipation). For large engineering calculations a second order scheme may be impractical because of the available computing resources. ¾ Try to give an estimate of the discretisation error in the simulation by applying a mesh refinement study (or if this is not possible by mesh coarsening) and by examining the effect of changes caused by the use of different order of spatial discretisations on the same mesh. ¾ If available in the code, make use of the calculation of an error estimator (which may be based on residuals, on the difference between two solutions of different order of accuracy or on wiggle size for a high order numerical scheme). ¾
3.6. Temporal discretisation errors for time-dependent simulations Purely steady flow fields with the time-derivative equal to zero are only a special case of the timedependent equations; in general, fluid flows are transient. The sources for this time-dependent behaviour are: •
External transient or non-transient forces.
•
Transient boundary conditions, moving walls (e.g. the fluttering of an airfoil).
•
Vortex stretching, a three-dimensional phenomenon due to the non-linear term of the governing equations, which also gives rise to the fluctuating nature of turbulence.
•
Vortex shedding.
The computation of steady turbulent flow is the most common kind of simulation for the industrial use of CFD. In these cases the Reynolds-averaged flow is steady while the average turbulent quantities account for the time-dependence of the turbulent fluctuations. In this way, the RANS equations represent the temporal average of the flow. However, the RANS-equations also allow the time-dependent Reynolds-averaged flow fields to be computed, based on the assumption that the temporal average of the turbulent quantities is not affected by the global unsteadiness. This is physically realistic when the time scale of the turbulence is much smaller than the time scale of the mean unsteadiness. A limitation is, however, that the typical turbulence frequencies are much higher than typical externally imposed or internally generated mean-flow frequencies. If the turbulence and mean flow frequencies are in the same band, ensemble averaging loses its sense, and LES or DNS methods should be employed. A time-dependent simulation is always needed if the scale of eddies or vortices is large and is comparable to the dimensions of the geometry (e.g. the computation of a vortex-street behind bluff bodies like a car or a building). In some turbomachinery design applications special models involving averaging-techniques with mixing planes between the blade rows are used to remove the time dependence of the flow due to the relative motion of the blades. Such techniques allow the flow to be computed using common steady state solution methods. If an accurate spatial discretisation is applied, flows which are physically time-dependent will often fail to converge using a steady-state method. However, if integrated with a large time step, even physically unsteady flows can be converged to a steady state (for example, a cylinder in a cross flow). Very often convergence problems with a steady simulation can be interpreted as a hint that the flow is unsteady and a time-stepping scheme would be appropriate. On the other hand, symmetry boundary conditions may impose a steady flow, although it would be transient in reality. If the complete geometry including both sides of the symmetry plane were used the velocity field would oscillate perpetually. Averaging the solution over a long time interval would lead to a symmetrical field, which, however, may differ from the steady state solution obtained with flow symmetry actually imposed.
© ERCOFTAC 2000
page 17 of 94
The temporal discretisation scheme provides an approximation of the time derivative. Most CFD codes offer implicit first and second order schemes, which are most effective in terms of computer memory and stability requirements. Low-storage higher-order Runge-Kutta methods are also available. The order of the scheme and the choice of the time-step influence the size of the amplitude and the phase error, the two components of the temporal discretisation error. To improve time-accuracy self-adaptive time-stepping procedures (such as predictor-corrector methods) can be used. The choice of the time step depends on the time scales of the flow being analysed. If time steps are too large the simulation might fail to capture important flow and mimic unphysical steady behaviour. It is therefore advisable to start with relatively small time steps, and corresponding CFL number, even though this is not required from the point of view of numerical stability. Some CFD codes use a timestepping scheme for steady state simulations.
3.6.1. Guidelines on temporal discretisation Second order accuracy is recommended in both space and time, as the overall solution accuracy is determined by the lower order component of the discretisation. For time dependent flows the time and space discretisation errors are strongly coupled. Hence finer grids, smaller time steps or higher order schemes are required (in both space and time). ¾ Check the influence of the order of the temporal discretisation by analysis of the frequency and time-development of a quantity of interest (e.g. the velocity in the main flow direction). ¾ Check the influence of the time-step on the results. ¾ Ensure that the time-step is adapted to the choice of the grid and the requested temporal size by resolving the frequency of the realistic flow and ensure that it complies with specific stability requirements. ¾
3.7. Types of implementation of boundary conditions Two types of boundary conditions and combinations of them are most commonly encountered. The Dirichlet condition specifies the distribution of a physical quantity over the boundary at a given time step and the Neumann condition defines the distribution of its first derivative. Fluid flow boundary conditions are more complex than in simple pure heat conduction simulations due to the coupling of the velocity fields with the pressure distribution. Users have normally no control on the spatial discretisation in the neighbourhood of boundaries. The CFD code developer should ensure that the implementation in the boundary region retains the overall accuracy of the numerical scheme. There is common consent that good practice for outflow boundaries is to set the convective derivative normal to the boundary face equal to zero and to combine this with a stream-wise extrapolation of transported quantities. Outflow and inflow boundaries bring about the following difficulties: •
Non-physical reflection of outgoing information back into the domain (Giles [1990]).
•
Difficulties in providing information about the properties of the fluid which may inadvertently enter the domain from the outside, and some CFD codes prevent fluid from entering the domain through outflow boundaries.
•
Difficulties may also arise if open boundary conditions are placed in regions of high swirl, large curvatures or pressure gradients.
3.7.1. Guidelines on the implementation of boundary conditions Ensure that appropriate boundary conditions are available for the case being considered. For swirling flows at an outlet consult the code manual to check that appropriate boundary conditions have been implemented (for example, radial equilibrium of pressure field instead of constant static pressure). Special non-reflecting boundary conditions are sometimes required for outflow and inflow boundaries where there are strong pressure gradients (Giles [1990]). ¾ Check whether the CFD code allows inflow at open boundary conditions. If inflow cannot be avoided at an open boundary then ensure that the transported properties of the incoming fluid including turbulence properties are properly modelled. ¾
© ERCOFTAC 2000
page 18 of 94
3.8. Special remarks on the finite-element method The best practice guidelines are based mainly on the finite-volume method which is the most common numerical method in commercial CFD-codes of industrial interest. Nevertheless a non-negligible number of CFD-codes based on the finite-element method are used by the CFD-community and this section gives additional advice that is specific to these. In the finite-element method (FEM) the domain is broken into a set of discrete finite elements that are usually unstructured. Their special feature is that the Navier-Stokes equations are multiplied by a weight (or basis) function before they are integrated over the entire domain. The weight function is defined locally in each element and it can be of higher order. The approximation is substituted into a weighted integral of the conservation laws. The result is a set of non-linear algebraic equations. An advantage of the early finite-element methods was the ability to deal with complex geometries where the grids can be easily refined in the area of interest. It can also be shown that the finiteelement methods have optimal properties for certain types of equations. An additional advantage is the mathematical background, which allows the development of interesting numerical possibilities, such as error estimators based on the residuals, and the use of higher order schemes. The principal drawback is that one needs to calculate and store the basis functions for each type of element that might appear in the grid (even if it only appears once). Moreover, there might be some difficulties related to conservation when using weak formulations (global conservation would be satisfied but local conservation would not). In addition the matrices of the linearised equations are not as well structured as in finite volume methods and this makes it more difficult to find efficient solution algorithms, although these are evolving rapidly through the use of iterative methods with pre-conditioners. For further information on the finite-element methods in fluid mechanics see Girault and Raviart [1986], Oden et. al [1975] and Zienkiewicz and Taylor [1991]. Most of the guidelines on the numerical methods mentioned above are also valid for the finite-element method. Some aspects may be either more or less important with the finite-element method than with the finite-volume method. The main differences between the guidelines for finite element and finite volume simulations are described in the next section.
3.8.1. Specific guidelines on finite-element methods in CFD ¾
¾ ¾ ¾ ¾ ¾
The requirements for the quality of the mesh of a finite-element model are less stringent than those for finite volume methods (e.g. expansion rate, aspect ratio) but there might be problems of local conservation of conserved variables if the quality of the mesh is too low. The code manual should define any special requirements of the mesh. The distribution of elements needs to be able to resolve the flow variation. In particular, the grid lines should be roughly aligned with the flow direction near to walls. Be aware that the grid might affect the flow at the degrees of freedom located on sharp corners (more so than in finite volume methods). The shape of the elements should be acceptable. Minimise the mesh bandwidth. Finite-element computations need more computer memory than Finite-Volume methods on structured grids with same number of nodes. Accurate FEM solutions can be obtained with a minimal mesh by using a quadratic basis function.
© ERCOFTAC 2000
page 19 of 94
4. Turbulence modelling Most flows of practical engineering interest are turbulent, and turbulent mixing then usually dominates the behaviour of the flow. Turbulence plays a crucial part in the determination of many relevant engineering parameters, such as frictional drag, heat transfer, flow separation, transition from laminar to turbulent flow, thickness of boundary layers, extent of secondary flows, and spreading of jets and wakes. The turbulent states which can be encountered across the whole range of industrially relevant flows are rich, complex and varied. After a century of intensive theoretical and experimental research, it is now accepted that no single turbulence model can span these states and that there is no generally valid universal model of turbulence. A bewildering number and variety of models have appeared over the years, as different developers have tried to introduce improvements to the models that are available. The extremely difficult nature of this endeavour caused Bradshaw [1994] to refer to turbulence as “the invention of the Devil on the seventh day of creation, when the Good Lord wasn't looking”. The available turbulence models can be roughly divided into three main categories, with the following subdivisions: •
Linear eddy viscosity models • Algebraic (or zero-equation) models • One-equation models (with ordinary and/or partial differential equations) • Two-equation models
•
Reynolds stress or second moment closure models • Full Reynolds stress transport models • Algebraic Reynolds stress models
•
Non-linear eddy viscosity models (NLEVM) • These can be sub-divided in the same way as linear eddy viscosity models (above) but are usually implemented as two-equation models • A certain sub-class of NLEVM are equivalent to algebraic Reynolds stress models
Within each of these categories there are a wide variety of different models and options available (see below). The choice of which turbulence model to use and the interpretation of its performance (i.e. establishing bounds on key predicted parameters) is a far-from-trivial matter. A set of application procedures is evidently required which documents the performance of various turbulence models across a broad class of flow regimes and for different applications. A full categorisation of this type is beyond the scope of the present guidelines. Instead the general features and broad limitations of different classes of model will be discussed and guidance will be given on the practical deployment of the turbulence model most commonly used in industrial practice, the standard k-ε model. A fuller introduction to the subject of turbulence and turbulence modelling can be obtained by consulting standard reference texts on the subject, such as Launder and Spalding [1972], Cebeci and Smith [1974], Rodi [1993], Tennekes and Lumley [1972], Menter et al. [1997], Wilcox [1998], Hanjalic [1994], Leschziner [1990] and Launder [1989].
4.1. RANS equations and turbulence models Turbulent flows contain many unsteady eddies covering a range of sizes and time scales. For flows in industrial applications, the effects of turbulence are nearly always examined using the RANS equations. These are developed from the time-dependent three-dimensional Navier-Stokes equations which are averaged in such a manner that unsteady structures of small sizes in space and time are eliminated and become expressed by their mean effects on the flow through the so-called Reynolds or turbulent stresses. These stresses must be interpreted in terms of averaged variables in order to close the system of equations. This requires the construction of a mathematical model known as a turbulence model. If heat transfer is important, the time-dependent energy equation is averaged in a similar manner giving rise to turbulent heat fluxes which account for the transport of energy due to the action of the unre-
© ERCOFTAC 2000
page 20 of 94
solved turbulent motions. The turbulence model must be extended to include a representation of these energy fluxes in order to fully close the system of equations2. Because models are based on different assumptions, all available turbulence models have limitations which depend on the modelling strategy.
4.2. Classes of turbulence models 4.2.1. Linear Eddy viscosity models The simplest turbulence modelling approach rests on the concept of a dynamic turbulent viscosity, µT. This relates the turbulent stresses appearing in the RANS equations to the averaged velocity gradients (i.e. the rate of strain) in direct analogy to the classical interpretation of viscous stresses in laminar flow by means of the fluid viscosity, µ. Thus for example, in a shear layer where the dominant velocity gradient is ∂u/∂y (u is the averaged velocity in the principal direction of flow and y is the crossstream co-ordinate) the turbulent shear stress is given as µT·∂u/∂y. The viscous analogy can be extended to the interpretation of the turbulent energy fluxes by the use of the so-called Reynolds extended analogy. A turbulent Prandtl number is introduced, PrT, so that the energy fluxes are related to time-averaged temperature gradients via a turbulent diffusivity of energy given as µT/(ρ·PrT). Unfortunately, unlike their viscous counterparts, the turbulent viscosity and turbulent Prandtl number are not fluid properties but functions of the state of turbulence. From dimensional considerations, µT/ρ is proportional to V·L, where V is a velocity scale and L is a length scale of the larger turbulent motions (called the mixing length in so-called mixing length models). Both the velocity scale V and the length scale L are determined by the state of turbulence, and, over the years, various prescriptions for V and L have been proposed.
4.2.2. Algebraic (or zero-equation) models The simplest prescription of V and L is with the so called algebraic (or zero-equation) class of models. These assume that V and L can be related by algebraic equations to the local properties of the flow. This is fairly straightforward for simple flows but can often be difficult in geometrically complex configurations (see Cebeci and Smith [1974], Baldwin and Lomax [1978] and the discussion in Wilcox [1998]). For example, in a wake or free shear layer V is often taken as proportional to the velocity difference across the flow and L is taken as constant and proportional to the width of the layer. In a boundary layer close to the wall V is given as L·∂u/∂y (or L·Ω where Ω is the magnitude of the vorticity) and L is related to the wall-normal distance from the wall (y-direction). The outer part of the boundary layer is treated in a similar manner to a wake. The turbulent Prandtl number is given a constant value except very close to a wall where molecular effects become important. Algebraic models of turbulence have the virtue of simplicity and are widely used with considerable success for simple shear flows such as attached boundary layers, jets and wakes. For more complex flows where the state of turbulence is not locally determined but related to the upstream history of the flow, a more sophisticated prescription is required. Algebraic turbulence models are based on boundary layer concepts, like shear layer thickness, distance from a wall and velocity differences across the layer. These quantities cannot easily be computed in a Navier-Stokes method and introduce significant additional uncertainty into the computation, as the solutions are dependent on implementation details. This is the main reason why these models are not normally recommended for general applications of RANS methods.
4.2.3. One-equation models The one-equation models attempt to improve on the zero-equation models by using an eddy viscosity that no longer depends purely on the local flow conditions but takes into account where the flow has come from, i.e. upon the flow history. The majority of approaches seek to determine V and L separately and then construct µT/ρ as the product of V and L. Almost without exception, V is identified with 2
If the transport of one or more scalar contaminants or species is important, the appropriate timedependant concentration equations are treated in a similar way to the energy equation. The turbulent mass fluxes which then arise are usually modelled following the same approach adopted for the energy fluxes. Thus, for the sake of simplicity and clarity, the text will focus on the treatment of the energy equation. © ERCOFTAC 2000
page 21 of 94
k1/2, where k is the kinetic energy per unit mass of fluid arising from the turbulent fluctuations in velocity around the averaged velocity. A transport equation for k can be derived from the Navier-Stokes equations and this is the single transport equation in the one-equation model. This is closed (i.e. reduced to a form involving only calculated variables) by introducing simple modelling assumptions thereby furnishing a prescription for the velocity scale V which accounts at least partially for non-local effects. In shear-layer type flows, and especially in regions close to a wall, it is often possible to algebraically prescribe L with reasonable confidence. In geometrically complex configurations the prescription of L is difficult because it depends on non-local quantities, such as the boundary layer thickness, displacement thickness, etc. and introduces similar uncertainties as in an algebraic turbulence model. Spalart and Allmaras [1992] have devised an alternative formulation of a one-equation model specifically for aerodynamic flows, which determines the turbulent viscosity directly from a single transport equation for µT and this model is proving quite successful for practical turbulent flows in external airfoil applications. The model is not well-suited for more general flows, as it leads to serious errors even for simple shear flows (round jet).
4.2.4. Two-equation models For general applications, it is usual to solve two separate transport equations to determine V and L, giving rise to the name two-equation model. In combination with the transport equation for k, an additional transport equation is solved for a quantity which determines the length scale L. This class of models is the best known and the most widely used in industrial applications since it is the simplest level of closure which does not require geometry or flow regime dependent input. The most popular version of two equation models is the k-ε model, where ε is the rate at which turbulent energy is dissipated to smaller eddies (Launder and Spalding [1974]). This is included in almost every commercial and in-house code and is a de-facto standard in industrial applications. A modelled transport equation for ε is solved and then L is determined as Cµ k3/2/ε where Cµ is a constant. Several commercial codes also include variants of the standard k-ε model, which improve predictions under special circumstances as explained in 4.4.1 below. The second most widely used type of two equation model is the k-ω model, where ω is a frequency of the large eddies (Wilcox [1998]). A modelled transport equation for ω is solved and L is then determined from k1/2/ω. The k-ω model performs very well close to walls in boundary layer flows, particularly under strong adverse pressure gradients (hence its popularity in aerospace applications). However it is very sensitive to the free stream value of ω and unless great care is taken in setting this value, spurious results can be obtained in both boundary layer flows and free shear flows. The k-ε model is less sensitive to free stream values but is often inadequate in adverse pressure gradients. Many different variants are to be found in the literature to circumvent this problem. An interesting option has been proposed by Menter [1994a, 1994b] through a model which retains the properties of k-ω close to the wall and gradually blends into the k-ε model away from the wall. This model has been shown to eliminate the free stream sensitivity problem without sacrificing the k-ω near wall performance. The performance of two-equation turbulence models deteriorates when the turbulence structure is no longer close to local equilibrium. This occurs when the production of turbulence energy departs significantly from the rate at which it is dissipated at the small scales (i.e. ε), or equivalently when dimensionless strain rates (i.e. absolute value of the rate of strain times k/ε) become large. Various attempts have been made to modify two equation turbulence models to account for strong non-equilibrium effects. For example, the so-called SST (shear stress transport) variation of Menter's model [1993, 1996] leads to marked improvements in performance for non-equilibrium boundary layer regions such as those found behind shocks or close to separation. However, such modifications should not be viewed as a universal cure. For example SST is less able to deal with flow recovery following reattachment.
4.2.5. Reynolds stress or second moment closure models The eddy viscosity models described above are based on the eddy viscosity concept and assume that the turbulent stresses are linearly related to the rate of strain by a scalar turbulent viscosity, and that the principal strain directions are aligned to the principal stress directions. This is reasonable for fairly © ERCOFTAC 2000
page 22 of 94
simple states of strain, especially when the model constants have been carefully calibrated from similar classes of flows, but may prove inadequate for modelling complex strain fields arising from the action of swirl, body forces such as buoyancy or extreme geometrical complexity. Under such circumstances a more subtle relationship between stress and strain must be invoked. The so-called Reynolds stress transport models (RSTM) dispense with notion of turbulent viscosity, and determine the turbulent stresses directly by solving a transport equation for each stress component. This requires the solution of six additional coupled equations, together with an equation for ε to provide a length scale (Launder and Spalding [1972], Rodi [1993], Launder et al. [1975], Leschziner [1998], and Speziale [1987a]). These methods are also sometimes referred to as second moment closure models (SMC) or second order closure models as the equations can be considered as second moments of the momentum equations. In a similar way, the turbulent heat fluxes can be determined directly by solving three extra equations, one for each flux component thereby removing the notion of a turbulent Prandtl number. This type of model can handle complex strain and, in principle, can cope with non-equilibrium flows. However, it is complex and expensive to compute. A much simplified form of second moment closure can be derived from the RSTM which reduces the number of coupled differential equations to two (as in two-equation turbulence models) whilst retaining the subtle interaction between stress and strain. This form, known as an algebraic Reynolds stress model (ARSM), is generated by assuming the transport (convection plus diffusion) of Reynolds stress (RS) in the RSTM equations is equal to the transport of turbulence energy multiplied by the factor (RS/k). The result is an implicit algebraic relationship for the Reynolds stresses which is closed by solving transport equations for k and a lengthscale determining quantity such as ε (Rodi [1993], Launder [1984] and [1989]). RSTM require boundary conditions for each of the stress components which can be difficult to specify. Also, both RSTM and ARSM can give rise to problems of convergence. Such models are offered in most commercial codes nowadays, but have not yet been widely adopted as an industrial tool, except for very particular types of flow for which extensive calibration has been carried out and where the standard k-ε model is known to do a poor job (see section 4.4).
4.2.6. Non-linear eddy viscosity models An alternative, somewhat simpler approach for dealing with complex strain is provided by the nonlinear eddy viscosity class of models, often referred to as NLEVM (see for example, Apsley et al. [1997]). These models retain the idea that the turbulent stresses can be algebraically related to the rate of strain (i.e. time averaged velocity gradients), but higher order quadratic and cubic terms are included which feature additional coefficients which can then be calibrated to experimental data. Such models are gaining in popularity since they involve the same number of equations as two equation models and thus are computationally efficient. They often provide improved results in flow impingement/reattachment areas where linear methods mostly fail. A sub-class of these models are can be shown to be formally equivalent to algebraic Reynolds stress models (Apsley et al. [1997], Gatski and Speziale [1993]). However, unlike ARSM, the stress components appear only on the left-hand side of the equation, obviating the need for an implicit solution process which, as mentioned above, often gives rise to convergence problems. For this reason models of this type are known as explicit algebraic Reynolds stress models (EARSM), and are gaining popularity as an effective way of implementing ARSM.
4.2.7. Low Reynolds number models Some turbulence models are valid for the turbulent flow region, but fail in the viscous sub-layer close to the wall. Various so-called low-Reynolds number versions of the k-ε and RSTM models have been proposed incorporating modifications which remove this limitation (see for example, Patel et al. [1985] and Wilcox [1998]). Alternatively the standard k-ε and RSTM models can be used in the interior of the flow and coupled to a one-equation model which is used to resolve just the wall region, a so-called two-layer model (see the review by Rodi [1991]).
4.2.8. Other models Turbulence modelling is an active area of research and many other models have been developed which are not referred to here. Some of these are quite promising, such as the k-ε-υ2 model of Durbin [1995], which is intermediate between the RSTM and the k-ε models. For further details the interested © ERCOFTAC 2000
page 23 of 94
reader is referred to the very extensive literature on this subject (starting with the references given in chapter 13).
4.2.9. Guidelines on turbulence modelling ¾
¾
¾ ¾
¾
Ensure that low numerical and convergence errors have been achieved in turbulent flow simulations. The relevance of turbulence modelling only becomes significant in CFD simulations when other sources of error, in particular the numerical and convergence errors, have been removed or properly controlled. Clearly no proper evaluation of the merits of different turbulence models can be made unless the discretisation error of the numerical algorithm is known, and grid sensitivity studies become crucial for all turbulent flow computations. Be aware that there is no universally valid general model of turbulence that is accurate for all classes of flows. Validation and calibration of the turbulence model with test data are necessary for all applications3. If possible, examine the effect and sensitivity of results to the turbulence model by changing the turbulence model being used. When using a particular turbulence model, check the published literature with regard to the known weaknesses of the model. The weaknesses of the standard k-ε model (Launder and Spalding [1974]), which is the most commonly used model in industrial applications, are listed below together with some indications of possible palliative actions which might be fruitfully considered. Select an appropriate near-wall model. Decide whether to use a wall function method, in which the near-wall region is bridged with wall functions (see the discussion below), or a low Reynolds number model, in which the flow structure in the viscous sub-layer is resolved. This decision will be based on the available resources and the requirements for resolution of the boundary layer. The validity of the wall function approach or the use of a low Reynolds number model should be examined for the flow configuration under study. Wall function methods are not valid in the presence of separated regions and/or strong three dimensional flows.
4.3. Near-wall modelling In wall attached boundary layers, the normal gradients in some flow variables become large as the wall distance reduces to zero. A large number of mesh points packed close to the wall is required to resolve these gradients. Furthermore, as the wall is approached, turbulent fluctuations are suppressed and eventually viscous effects become important in the region known as the viscous sub-layer. This modified turbulence structure means that many standard turbulence models (see summary given above) are not valid all the way through to the wall. Thus special wall modelling procedures are required.
4.3.1. Wall functions This is the procedure most commonly used in industrial practice. The difficult near-wall region is not explicitly resolved with the numerical model but is bridged using so-called wall functions (Launder and Spalding [1974] and Wilcox [1998]). Attention will focus on smooth walls, the alternative treatment of rough walls being mentioned as a special case. In order to construct these functions the region close to the wall is characterised in terms of variables rendered dimensionless with respect to conditions at the wall. The wall friction velocity uτ is defined as (τw/ρ)1/2 where τw is the wall shear stress. Let y be normal distance from the wall and let U be time-averaged velocity parallel to the wall. Then the dimensionless velocity, U+ and dimensionless wall distance, y+ are defined as U/uτ and y·ρ·uτ /µ respectively. If the flow close to the wall is determined by conditions at the wall then U+ can be expected to be a universal (wall) function of y+ up to some limiting value of y+. This is indeed observed in practice, with a linear relationship between U+ and y+ in the viscous sub-layer, and a logarithmic relationship, known as the law of the wall, in the layers adjacent to this (so-called log-layer). For rough walls, this law of the wall must be modified by scaling y on the equivalent roughness height, z0 (i.e. y+ is replaced by y/z0 ), and also by adjusting the coefficients (Schlichting [1979]). The y+-limit of validity depends on external factors 3
Calibration in this sentence refers to testing the ability of a CFD code to predict global quantities of interest for specific geometries of engineering design interest (see definition in chapter 2) and not to tuning the coefficients of a turbulence model, which cannot be recommended for a novice user. © ERCOFTAC 2000
page 24 of 94
such as pressure gradient and the penetration of far field influences. In some circumstances the range of validity may also be affected by local influences such as buoyancy forces if there is strong heat transfer at the wall. Turbulent quantities, such as the turbulence velocity (k1/2), dissipation and length scales, when treated in the same way, also exhibit a universal behaviour. For example, in the loglayer, under equilibrium conditions k/uτ2 is a constant, C. It is of interest to note this allows an alternative definition of dimensionless wall distance, denoted y* where y* is given by y ρ(C-1 k)1/2/µ . The two forms are equivalent (i.e. y+ = y*) when the turbulence structure is in equilibrium (see for details Chieng and Launder [1980] and Grotjans and Menter [1998]). The role of y* in the application of wall functions is discussed further below. The standard wall functions are valid for smooth walls, but can be modified to allow for rough surfaces by adjustment of the constants in the law of the wall (see Tennekes and Lumley [1972]). If a rough wall is being modelled the wall distance in the law of the wall is non-dimensionalised with an equivalent roughness height, and appropriate values must be selected from the data or literature. These universal functions can be used to relate flow variables at the first computational mesh point, displaced some distance y from the wall, directly to the wall shear stress without resolving the structure in between. A similar universal, non-dimensional function can be constructed which relates the temperature difference between the wall and the mesh point to heat flux at the wall (Launder and Spalding [1974]). This can be used to bridge the near-wall region when solving the energy equation. The constraints to be observed when using wall functions are: a) The calculated flow must be consistent or nearly consistent with the assumptions made in arriving at the wall function equations, i.e. an attached two-dimensional Couette flow with small pressure gradients, local equilibrium of turbulence (production rate of k equals its dissipation rate), and a constant near-wall stress layer. Applying wall functions outside this application range, e.g. for general three-dimensional flows, separating flows, swirling flows and flows along spinning walls, may lead to inaccuracies. b) The value of y+ at the first mesh point remains within the limit of validity of the wall functions. Unfortunately, it is not always possible to satisfy these constraints everywhere within a flow and then certain practical advantages (rarely improvements in accuracy) can be gained from basing wall distance appearing in the wall-functions on y* rather than y+ . For example, at flow re-attachment or at point of impingement, τw vanishes and so y+ is zero irrespective of the value of y (an example of far field influences penetrating through to the wall reducing the y+ limit of validity to zero). However, y* remains large due to the high values of k convected towards the wall, facilitating the (apparent) satisfaction of constraint b). Also, in this particular case the use of y* rather than y+ results in a qualitative improvement in the prediction of wall heat transfer. Wherever the conditions of constraint a) are met, y+ and y* are equivalent. Many commercial codes employ y* in the application of wall functions as the default4. Standard wall functions are one of the biggest sources of misconceptions in turbulent flow computations, even for experienced users. Their purpose is to bridge the extremely thin viscous layer near the surface. They do not free the user from the need to adequately resolve the turbulent portion of the boundary layer. Adequate resolution means at least eight to ten points in the turbulent regime, if the accurate prediction of boundary layer effects, like wall friction, heat transfer and separation is intended. In addition to the lower limit on y+, which ensures that the first (integration) point does not fall into the viscous sublayer, there are also upper limits due to the requirement of adequate boundary layer resolution. These limits are the more restrictive the lower the running-length of the boundary layer is. For boundary layers with a Reynolds number lower than Re=105 it is difficult to satisfy the requirement that the first (integration) point has to lay outside the viscous sublayer and at the same time to achieve a sufficient boundary layer resolution. The most sensible approach to satisfy the logarithmic profile assumption and to resolve the boundary layer is to place the first (integration) point slightly above the lower limit recommended by the code developers (typically around y+ = 20).
4
Throughout Chapter 4, y+ will be used to represent either y·ρ·uτ /µ or y* without distinction, depending upon the wall function implementation adopted. © ERCOFTAC 2000
page 25 of 94
To generate a grid with a pre-specified y+ distribution is difficult, as y+ depends on the solution, which is not known during the grid generation phase. As many engineers compute similar flow in similar geometries at similar Reynolds numbers, previous computations can serve as a guideline. Standard wall functions are clearly a weakness of present day CFD methods, as they require a very careful placement of the near wall grid points. In addition, they do not allow for a consistent grid refinement as the solutions deteriorate if the near wall grid is refined below the log-law limit. Also they are rarely accurate in all regions of the flow.
4.3.2. Specific guidelines on wall functions ¾
¾
¾
¾
¾
Check the lower limit of y+. In the commonly used applications of wall functions, the meshing should be arranged so that the values of y+ at all the wall-adjacent integration points is only slightly above the recommended lower limit given by the code developers, typically between 20 and 30 (the form usually assumed for the wall functions is not valid much below these values). This procedure offers the best chances to resolve the turbulent portion of the boundary layer. It should be noted that this criterion is impossible to satisfy close to separation or reattachment zones unless y+ is based upon y*. Check the upper limit on y+. In the case of moderate Reynolds number, where the boundary layer only extends to y+ of 300 to 500, there is no chance of accurately resolving the boundary layer if the first integration point is placed at a location with the value of y+ of 100. Check the resolution of the boundary layer. If boundary layer effects are important, it is recommended that the resolution of the boundary layer is checked after the computation. This can be achieved by a plot of the ratio between the turbulent to the molecular viscosity, which is high inside the boundary layer. Adequate boundary layer resolution requires at least 8-10 points in the layer. Exercise care when calculating the flow using different schemes or different codes with wall functions on the same mesh. Cell centred schemes have their integration points at different locations in a mesh cell than cell vertex schemes. Thus the y+ value associated with a wall-adjacent cell differs according to which scheme is being used on the mesh. Check that the correct form of the wall function is being used to take into account the wall roughness. An equivalent roughness height and a modified multiplier in the law of the wall must be used.
4.3.3. Near-wall resolution with low Reynolds number models As already mentioned, a universal near-wall behaviour over a practical range of y+ may not be realisable everywhere in a flow. Under such circumstances the wall-function concept breaks down and its use will lead to significant errors, particularly if wall friction and heat transfer rates are important. The alternative is to fully resolve the flow through to the wall by the selection of a suitable low Reynolds number wall resolving model, but the cost of the solution is around an order of magnitude greater than when a wall function is used because of the extra mesh involved Many turbulence models can be validly used for this purpose, others cannot. For example, the k-ω two-equation model can be deployed through to the wall as can the one-equation k-L model (e.g. Wolfshtein [1969]). The standard k-ε and RSTM models cannot. Various so-called low-Reynolds number versions of the k-ε and RSTM models have been proposed incorporating modifications which remove this limitation (see, for example, the reviews in Patel et al. [1985] and Wilcox [1998]). Alternatively the standard k-ε and RSTM models can be used in the interior of the flow and coupled to the k-L model which is used to resolve just the wall region (see the review by Rodi [1991]). Whatever modelling approach is adopted, a sufficient number of mesh points must be packed into a very narrow region adjacent to the wall in order to capture the variation in the flow variables.
4.3.4. Specific guidelines on low Reynolds number models Choose carefully the value of y+ at the first node adjacent to the wall. This should to be less than 4, and preferably close to unity. ¾ Depending on the Reynolds number, ensure that there are between five and ten mesh points between the wall and y+ equal to 20, which likely results in thirty to sixty points inside the boundary layer for adequate boundary layer resolution. ¾
© ERCOFTAC 2000
page 26 of 94
4.4. Weaknesses of the standard k-ε model with wall functions Despite the great variety of turbulence modelling options available to the user, the standard k-ε model with wall functions, as set out by Launder and Spalding [1974], remains the work-horse of industrial computation. It is therefore of value to catalogue the major weaknesses associated with this model in practical application and, where possible, indicate palliative actions which might be fruitfully considered. These are listed below. The advisory actions are drawn from an extensive literature on the subject and should not to be viewed as definitive cures. There is no a priori guarantee that the advice given will yield significantly improved results, and careful validation is always appropriate. The manuals of commercial and inhouse codes may proffer alternative and equally effective advice, and many commercial codes will include alternatives to the standard k-ε model. Where the action given below involves a modification or adjustment to the standard k-ε model, this should be regarded as a specific palliative for the weakness under consideration and will not usually prove of general benefit (and may even be worse).
4.4.1. Guidelines on weaknesses of the standard k-ε model ¾
¾
¾ ¾
¾
The turbulent kinetic energy is over-predicted in regions of flow impingement and re-attachment leading to poor prediction of heat transfer and the development of boundary layer flow around leading edges and bluff bodies. The high turbulence levels predicted upstream of a stagnation point are transported around the body and the real boundary layer development is swamped by this effect. Kato and Launder [1993] have proposed an ad-hoc modification (called the Yap correction) to the equation for the production term of k appearing in the k and ε equation, which is designed to tackle this problem. The RNG k-ε model (Yakhot et al. [1992]) includes a modification to the transport equation for ε and may also improve predictions in this area. The problems depend on the free-stream values of k and ε and may not always occur. Highly swirling flows are often poorly predicted due to the complex strain fields and regions of recirculation in a swirling flow are often under-estimated. These deficiencies require turbulence models with improved physics, such as non-linear k-ε models, algebraic Reynolds stress models, or full Reynolds stress models with low Reynolds number capability limitation (for example, the reviews in Patel et al. [1985] and Wilcox [1998]). Mixing is often poorly predicted in flows with strong buoyancy effects or high streamline curvature. Again this deficiency requires improved physics, see comment above. Flow separation from surfaces under the action of adverse pressure gradients is often poorly predicted. The real flow is likely to be much closer to separation (or more separated) than the calculations suggest. The SST version of Menter's k-ω based, near wall resolved model (Menter [1993, 1996]) can offer a considerable improvement. Flow recovery following re-attachment is often poorly predicted. If possible, avoid the use of wall functions in these regions, but low Reynolds number models may also have some problems. A promising possibility is the use of a length-scale limiting device (Ince and Launder [1995]).
The far-field spreading rates of round jets are predicted incorrectly. The use of non-linear k-ε models should be investigated for these problems (Apsley et al. [1995]). ¾ Turbulence driven secondary flows in straight ducts of non-circular cross section are not predicted at all. Linear eddy viscosity models cannot capture this feature. Use RSTM or non-linear eddy viscosity modelling. ¾ Dilation effects on turbulence production, which may become important for compressible flows, are not accounted for. These may be important in mixing layers at moderate Mach numbers. Sarkar et al. [1989] and Zeman [1990] have proposed modifications to the k-ε model that correct this deficiency but the database is small. ¾
Laminar and transitional regions of flow cannot be calculated with the standard k-ε model but nonlinear modelling may help in this respect. Transition is an active area of research in turbulence modelling. Some simple practical advice is given in section 4.7 below. ¾ The vertical heat flux in a buoyancy affected horizontal boundary layer, or in a thermal plume, is poorly represented by an eddy viscosity model which assumes a turbulent Prandtl number to calculate the turbulent transport of heat. RSTM models including the solution of equations for the three heat fluxes is currently the only option if this is important. ¾
© ERCOFTAC 2000
page 27 of 94
4.5. Inflow boundary conditions The use of a turbulence model (other than an algebraic model) requires that turbulence properties at a domain inlet region need to be specified. Verified quantities should be used as inlet boundary conditions for k and ε, because their magnitude can significantly influence the results. If there are no data available, then the influence of the choice should be examined by sensitivity tests with different simulations. The specification of turbulence properties is an important aspect of application uncertainties, and is discussed in more detail in the section 5.2.3.
4.6. Unsteady flows The use of the RANS equations in an unsteady flow is valid, provided that the time-scale of the unsteadiness is sufficiently far removed from the turbulence scales. If the unsteadiness is provoked from an external source (such as an internal combustion engine, wake and blade interaction, pitching airfoils) the turbulence model does not necessarily interact with the frequency and amplitude of the unsteadiness. However, in unsteady flows all properties of a turbulence model are put under severe scrutiny, with the result that the predictive capability of a model is generally significantly reduced in comparison with the equivalent steady state conditions. This is clearly observed in all publications on oscillating airfoils, for instance, where large discrepancies with experimental data are observed even at moderate angular amplitudes (Ekaterinaris and Menter [1994]). If the unsteadiness is self-induced, such as vortex shedding from a bluff body or highly separated flow from an airfoil at high incidence, then the difficulties with turbulence modelling are much worse. There is presently no known turbulence model able to handle correctly highly stalled airfoils. See also section 3.6 on temporal discretisation errors.
4.7. Laminar and transitional flows The distinction between laminar, transitional and turbulent flow is difficult for a number of applications. Sometimes the flow appears in different states depending on the location of the area of interest, for example, the flow in an inlet of a machine can be laminar whereas the flow inside the machine is turbulent. Another example is the flow over an airfoil, which is normally laminar at the leading edge undergoes transition to a turbulent state downstream and is turbulent in the wake behind. The general problem of the transition from laminar to turbulent flow, that is the computation of the onset of transition as well as the length of the transition zone is a subject of fundamental research. It cannot presently be included in general industrial CFD-computations. The simplest way around this problem is to calculate the flow as a turbulent one or to prescribe the location of transition if it is available from experiment. If a k-ε model with a low Reynolds number behaviour is applied, then the turbulent kinetic energy is approximately zero in the nominal laminar flow regimes, and in this case the models are able to predict a transitional behaviour. However, much care has to be exercised, as this does not ensure that the physics of the transition is respected by the turbulence model. The physics can be very complex and highly different according to the type of transition, such as by-pass transition, separation induced transition, etc. (see, for example Mayle [1991] and Savill [1996]).
4.7.1. Guidelines on transitional flows Use experimental data to check whether the flow contains extensive regions of laminar or transition flow which would be incorrectly estimated by the k-ε model with wall functions, or poorly predicted by low Reynolds number k-ε models. ¾ Model transition by intervention in the code to switch the turbulence model on or off at predetermined locations determined by experiment. ¾
© ERCOFTAC 2000
page 28 of 94
5. Application uncertainties 5.1. Geometrical uncertainties In many industrial and engineering problems, the geometry of the object to be simulated is extremely complex and requires much effort to specify it exactly for a computer simulation. There are many sources of error which can arise in this process, such as: •
Changes in geometry that have occurred during the design process have been neglected.
•
CAD geometry definition is insufficiently complete for flow simulation. Some surfaces and curves may not meet at the intended end point locations due to different levels of accuracy in different parts of the CAD model. Other curves may be duplicated.
•
The geometry of a tested component may get modified during the testing procedure, and these modifications may not have been added to the original drawings.
•
The geometry may not be manufactured within the tolerances as shown on the drawing, particularly with regard to fine flow features, such as the rounded shape of turbomachinery leading edges, or symmetrical features.
•
The change of shape due to flow loading or other mechanical loading may not have been taken into account, for example, the untwist of turbomachinery blade rows due to centrifugal forces.
•
The effective geometry of the surface may have changed during use due to wear, erosion or dirt fouling.
•
Small details of the geometry may have been omitted, such as tip clearance gaps in turbomachinery, roughness on the walls, welding fillet radii, details of furniture in building flow simulations, etc.
•
The co-ordinate system and units used in the CAD system may be different from that used in the CFD code (rotational direction).
5.1.1. Guidelines on geometrical uncertainties Check and document that the geometry of the object being calculated is the geometry as intended taking into account items such as design modifications, wear, tolerances, deflection due to loading, fouling, etc. If the geometry is substantially altered or deformed by the aerodynamic, mechanical or thermal loading, then use a some structural or mechanical calculation to determine the exact geometry. ¾ Check that the geometry is defined in the correct co-ordinate system and with the correct units requested by the CFD-code. ¾ Check that the transfer of geometrical data from a CAD system to a CFD-system, or grid generation software, does not involve loss of surface representation accuracy and that the final surface grid conforms to the intended geometry. CAD-systems often define the geometry in millimetres and this must be converted to SI-units if the code assumes that the geometry information is in these units (which is commonly done by most codes). It may be necessary to clean up the CAD geometry so that surfaces and curves meet at the intended end point locations. Visual display of the geometry helps here. ¾ Consider whether local details of the geometry can be omitted or require grid refinement. In general, it is not necessary to explicitly include geometrical features that have dimensions below that of the local grid size provided that they are taken into account in the modelling (e.g. roughness in wall layer, porous elements or detail of pipe-work in process plant explosion calculations). In areas where local detail is needed, such as in the neighbourhood of fine edges or small clearance gaps, then grid refinement in the areas with fine details should be used,. ¾
5.2. Boundary conditions In many real applications, there is a frequent difficulty in defining some of the boundary conditions at the inlet and outlet of a calculation domain in the detail that is needed for an accurate simulation. A typical example is the specification of the turbulence properties (turbulence intensity and length scale) at the inlet flow boundary, as these are seldom available for a new design configuration. Other examples are the specification of the boundary layer velocity profile on the walls at an inlet, or the precise distribution of a species concentration at an inlet boundary. The user needs to be aware of these problems and needs to develop a good feel for the certainty or uncertainty of the boundary conditions © ERCOFTAC 2000
page 29 of 94
which are imposed. This can best be achieved if the user knows and understands the application he is calculating. Additional uncertainties can arise because boundary condition data that needs to be specified is inconsistent with the model being used.
5.2.1. General guidelines on boundary conditions Examine the possibilities of moving the domain boundaries to a position where the boundary conditions are more readily identified, are well-posed and can be precisely specified. ¾ Check whether upstream or downstream obstacles (such as bends, contractions, diffusers, blade rows etc.) outside of the flow domain are present which might have significant effects on the flow distribution. Often information on the components upstream or downstream of the domain is lacking or not available at the beginning of a project. ¾ For each class of problem that is of interest, carry out a sensitivity analysis in which the boundary conditions are systematically changed within certain limits to see the variation in results. Should any of these variations prove to have a sensitive effect on the simulated results and lead to large changes in the simulation, then it is clearly necessary to obtain more accurate data on the boundary conditions that are specified. ¾ Be aware that the code may have default boundary conditions for regions of the domain boundaries where the user has specified nothing.
¾
5.2.2. Guidelines on inlet conditions ¾
Carry out a sensitivity analysis in which the key inlet boundary conditions are systematically changed within certain limits. Depending on the problem, the key parameters that might be examined are: • • • •
¾
Inlet flow direction and magnitude. Uniform distribution of a parameter or a profile specification, for example a uniform inlet velocity (slug flow) or an inlet velocity profile. Physical parameters. Turbulence properties at inlet.
Particular care is needed with the specification of inlet mass flow as a boundary condition for compressible flow-fields with supersonic regions, as the flow may be choked and unable to pass the mass flow specified.
5.2.3. Guidelines on specification of turbulence quantities at an inlet Use verified quantities from experiments as inlet boundary conditions for turbulent kinetic energy k and dissipation ε, if these are available, as the magnitude can significantly influence the results. If there are no experimental data available, then the values need to be specified using sensible engineering assumptions (see below), and the influence of the choice should be examined by sensitivity tests with different simulations. ¾ Specify values of the turbulent kinetic energy k which are appropriate to the application. These values are often specified through a turbulence intensity level Tu, which is defined by the ratio of the fluctuating component of the velocity to the mean velocity. In external aerodynamic flows over airfoils the turbulence level is typically Tu = 0.003 (0.3%). In atmospheric boundary layer flows the level can be two orders of magnitude higher (Tu = 0.30 (30%)) and details of the actual boundary layer profiles are needed. In internal flows the turbulence level of Tu = 0.05 to 0.10 (5 to 10%) is usually appropriate. ¾
¾
Specify values of the turbulent length scale, as an equivalent parameter for the dissipation ε, that are appropriate to the application. For external flows remote from boundary layers a value determined from the assumption that the ratio of turbulent and molecular viscosity µT/µ is between 1 and 10 is a reasonable guess. For internal flows a constant value of length scale derived from a characteristic geometrical feature can be used (e.g. 1 to 10% of the hydraulic diameter for internal flows). For simulations in which the near-wall region is modelled with a low Reynolds number model the length scale should be based on the distance to the wall and be consistent with the internal modelling in the code.
¾
If more sophisticated distributions of k and ε are used, check that these are consistent with the velocity profile. An inconsistent formulation may lead to an immediate unrealistic reduction of the turbulence quantities after the inlet.
© ERCOFTAC 2000
page 30 of 94
For RSTM models the stresses themselves need to be specified. If these are not available, as is often the case, make an assumption of isotropic flow conditions (that is with the normal stresses given by 2/3 k and with zero shear stresses). ¾ For inlets that are representative of fully developed pipe or channel flows, algebraic profiles might be used.
¾
Check the consistency of the definitions of k and ε by making a plot of the ratio of turbulent to molecular viscosity µT /µ. Note that the ratio of turbulent to laminar viscosity depends on Tu, L, k and ε as µT/µ ~ k2·ε or µT/µ ~ Tu·L. ¾ In cases where problems arise, move the inflow boundary sufficiently far from the region of interest so that a natural inlet boundary layer can develop.
¾
5.2.4. Guidelines on outlet conditions ¾
¾ ¾
¾
¾
¾
Place open boundary conditions (outflow, pressure) as far away from the region of interest as possible and avoid open boundaries in regions of strong geometrical changes or in regions of recirculation. If the latter point cannot be achieved, consult the code manual to see what boundary condition is considered best. Take special care with the orientation of outlet planes with regard to the mean flow, especially when the boundary condition consists of a constant static pressure profile. Select the boundary conditions imposed at the outlet to have only a weak influence on the upstream flow. Extreme care is needed when specifying flow velocities and directions on the outlet plane. The most suitable outflow conditions are weak formulations involving specification of static pressure at the outlet plane. Particular care should be taken in strongly swirling flows where the pressure distribution on the outlet boundary is strongly influenced by the swirl. It is therefore not acceptable to specify constant static pressure across the outlet. Be aware of the possibility of inlet flow inadvertently occurring at the outflow boundary, which may lead to difficulties in obtaining a stable solution or even to an incorrect solution. If it is not possible to avoid this by relocating the position of the outlet boundary in the domain, then one possibility to avoid this problem is to restrict the flow area at the outlet, provided that the outflow boundary is not near the region of interest. If the outflow boundary condition allows flow to re-enter the domain, appropriate Dirichlet conditions should be imposed for all transported variables. If there are multiple outlets, impose either pressure boundary conditions or mass flow specifications depending on the known quantities.
5.2.5. Guidelines on solid walls Care should be taken that the boundary conditions imposed on solid walls are consistent with both the physical and numerical models used (for example, adiabatic walls, local heat fluxes, free-slip boundaries, roughness, moving or rotating walls). In cases with moving or rotating walls the boundary conditions need to be specified consistent with the motion of the walls. ¾ Great care is needed if roughness on the wall is not negligible. If no detailed information is available, significant levels of uncertainty can arise through incorrect specification of equivalent roughness on solid walls. ¾
5.2.6. Guidelines on symmetry and periodicity planes ¾
If symmetry or periodicity planes cross the inlet or outlet boundaries then care should be taken to specify inlet or outlet variables that are consistent with these. Symmetry and periodicity planes assume that the gradients perpendicular to the plane are either zero (for symmetry) or determined from the flow field (periodicity).
5.3. Initial condition and initial guess The iterative solution methods calculate the flow from an initial estimate of the flow field. The initial guess can influence the convergence and, in some special cases, on the converged solution. In the case of transient calculations, the initial condition might be a valid solution of the flow equations. This can either be the trivial solution, i.e. all quantities set to zero, or a converged steady-state solution.
© ERCOFTAC 2000
page 31 of 94
5.3.1. Guidelines on initial guess or initial condition Be aware that in special cases the initial guess can influence the converged flow field e.g. buoyant flows, flows with bifurcations or flows with a hysteresis depending on the direction of the changes. ¾ Care should be taken to impose realistic conditions that are consistent with the flow equations as the initial condition of an unsteady-state flow calculation (for example, use a steady-state flow simulation as the initial condition for a gas release simulation).
¾
¾
Special attention should be directed to the initial guess of k and ε, especially with low Reynolds number models.
5.4. Uncertainties with regard to steady flow, symmetry, and periodicity A symmetric steady computation, or a computation with periodic boundary conditions, is often carried out in order to reduce the computing time and memory required. There are many applications where the nominal geometry is symmetric but the flow is asymmetric, and the flow field can be asymmetric even in the case of symmetry of the geometry (for example, an oblate spheroid at very high incidence). This can be an important factor in predicting the detail of the dynamical behaviour of fluid flows. Asymmetry can also be forced if the real inflow conditions are not geometrically perfectly symmetrical or some minor distortions are present within the inlet flow. Due to the physical temporal instability of the flow (for example, in a Karman vortex street or the periodic attachment and detachment of a hot jet below a cold ceiling) or due to time-dependent boundary conditions the flow field can be unsteady. This effect should be carefully examined because the flow solvers can often compute a spurious steady solution of the flow field which is in contradiction to the physics. In cases with very strong unsteady effects within the flow field the solution algorithm does not always converge to a steady solution.
5.4.1. Guidelines on uncertainties with steady flow, symmetry and periodicity ¾ ¾
¾ ¾ ¾
Check carefully whether the geometry is symmetric or whether a geometrical distortion or disturbance in the inlet conditions is present which can trigger asymmetric solutions. Estimate the Reynolds-number of the inflow and check whether the flow could be asymmetric, turbulent and/or unsteady (e.g. by sources or literature). Coupling between competitive physical phenomena might also be considered as a potential source of unsteadiness. After obtaining a steady solution, switch to the transient mode and check whether the solution remains stable. Switch to the transient mode if it is difficult to get a converged steady solution, especially if there is an oscillation of the residuals. In case of doubt, the simulation should be unsteady and without symmetry assumptions as boundary conditions.
5.5. Physical properties The commonly used commercial CFD-codes require the specification of the physical properties. The codes typically solve the equations in dimensional form.
5.5.1. Guidelines on physical properties Ensure that the correct physical properties are specified (e.g. water, air, ideal gas). If constant property assumptions are being used make sure that they are valid. ¾ Ensure that the specified value is the correct one (e.g. dynamic and not kinematic viscosity). ¾ Ensure that the values are specified in the correct units. ¾ Check whether user coding is being used to specify non-constant properties via tabulated values, as this can be an additional source of uncertainty (e.g. interpolation, extrapolation, user coding errors, inconsistency with base code). ¾
© ERCOFTAC 2000
page 32 of 94
6. User errors 6.1. General comments In CFD the human factor plays an important role, as the results depend to a large extent on the competence and expertise of the user. It is worthwhile spending a few words on this rather embarrassing aspect of CFD, as it is one of the prime causes of uncertainty in the results of CFD simulations. This may help to avoid some, if not all, of the most easily avoidable mistakes in the future. Several factors may give rise to user errors: •
Lack of attention to detail, sloppiness, carelessness, mistakes and blunders.
•
Too optimistic and uncritical use of CFD, thanks to the high accessibility through simple interactive graphical user interfaces in commercial software, and the convincing and seductive power of the colourful visualisations.
•
Lack of experience so that the user is unaware of a technical difficulty or unaware that critical information is missing.
•
Unfamiliarity with a particular CFD code or code version, and the tacit assumption that certain parameter settings are equivalent to those in a code or code version with which the user is more familiar.
While the first two points are associated with the user’s attitude and personal disposition the remaining points refer to the question of experience and training.
6.2. Control of the working process Many mistakes are made by mere lack of attention to detail, or because the user is not aware of factors that can give rise to them. The best way to deal with these issues is for the user to have a clear checklist of issues that can arise which helps to ensure that all relevant problem areas have been dealt with. This becomes most important if the user has limited experience. A formal management Quality Assurance (QA) system with checklists can help to support the inexperienced user to produce quality CFD simulations. It has been noted by Roache [1998], however, that a CFD project can meet all formal QA requirements and still be of low quality (or flatly erroneous). On the other hand high quality work can be done without a formal QA system. The guidelines given below provide examples of the sorts of issues that should be dealt with in a formal QA management system. The issues covered are based on the process of carrying out a CFD simulation as outlined in chapter 2.
6.2.1. Guidelines on problem definition ¾
Give careful thought at the outset of a CFD project to the requirements and objectives of the simulation. The following points should be considered: • •
• • •
• • •
Are the objectives of the simulation clearly defined? Is a CFD simulation really appropriate or is a simple analysis or an experiment needed? Circumstances where CFD is generally considered appropriate are when the flow physics or geometry are sufficiently complex that simpler modelling approaches are invalid, or where experimental data is lacking and cannot be obtained because it is too costly, dangerous or impractical. What local/global quantities are needed from the simulation? What are the requirements on accuracy? What level of validation is necessary? Is this a routine application, where validation and calibration has already been carried out on similar flow fields, and where only relatively small changes can be expected from earlier similar simulations? Or is it a non-routine application, where little earlier validation work has been done? What are the important flow physics involved (heat transfer, radiation, steady, unsteady, compressible, single phase, laminar, turbulent, transitional, internal, external, etc.)? What is the area of primary interest (domain) for the flow calculation? Is the geometry well defined, and can symmetry or periodic boundary conditions be applied?
© ERCOFTAC 2000
page 33 of 94
• • • •
Are the inlet and outlet boundary conditions clearly specified? What sensitivity tests should be carried out to assess the accuracy of the simulation? What level of computational resources is needed for the simulation (memory, disk space, cpu time, serial or parallel run because of time requirement) and are these available? What are the documentation/reporting requirements?
6.2.2. Guidelines on solution strategy ¾
Translate the problem definition into a solution strategy involving issues and questions that have been addressed in this document, such as: • • • • • • • •
Mathematical and physical models, numerical schemes. Solution algorithm, pressure or density based solution method. Turbulence and other physical models. Available code/solver. Grid design and computational mesh. Boundary conditions. Convergence criterion. Validation and sensitivity test requirements.
6.2.3. Guidelines on user errors in the code-handling ¾
Try to minimise errors in the handling of the code by the use of a formal check list or by letting another CFD analyst check through the code input data. This will minimise the potential source of user errors when implementing the solution strategy with a particular code. The types of questions which should be considered are: • • • • • • • • • • • • • • •
Is the correct version of the code being used? The default parameter settings or alternatives should be reviewed. Have default parameters been changed which may affect the solution? Are the user defined routines or other software coupled to the code compatible with the current version of the code? Have the boundary conditions not only been properly defined, but also properly applied? Has the appropriate system of units been used? Is the geometry correct? Are all parts of the mesh correctly matched and not overlapping, especially after importing grid parts from different sources? Are the correct physical properties specified? Have the intended physical and mathematical models been used (e.g. gravity forces, rotation, user defined functions)? Have the specified models been changed from those in an earlier code version? Have any user-defined subroutines been used, and have these been checked with regard to programming errors, units, and references to pre-defined variables? Has the appropriate convergence criterion been defined and used? Is there a reasonable maximum iteration number defined, especially in connection with a cputime limit for batch jobs? Have any warnings been generated by the software, and have they been properly understood and acted upon? Have the necessary files from other computer systems been transferred correctly?
6.2.4. Guidelines on interpretation ¾
Make full use of the wide variety of visual tools available in modern CFD packages to interpret the flow-field that has been calculated. A good physical understanding of the application and experience with similar simulations are the best guide as to which plots to make. Options include: • • • • • •
Contours of parameters on arbitrary cross-sections through the flow. Contours of parameters on wetted surfaces in flow, such as wall pressures. Special plots involving projections for specific applications, such as blade-to-blade or meridional views in turbomachinery blade rows. Line diagrams of parameters. Three-dimensional plots of streamlines. Vector plots.
© ERCOFTAC 2000
page 34 of 94
• • • •
Three-dimensional views of iso-surfaces, such as regions of low pressure or high mach number. Surface streak lines. Use of identification methods to visualise vortical structures, such as line integral convolution (LIC) Video sequences of unsteady simulations.
Produce plots of parameters that are relevant to the physical understanding of the flow (such as pressure fields, velocities, entropy, losses, wall shear stresses, regions of reverse flow, etc.), and also of parameters that might help in the interpretation of the quality of the simulation (such as the ratio of turbulent to laminar viscosity, grid skewness, grid aspect ratio, residuals, etc.). ¾ Don’t be seduced into believing that the solution is correct just because it has converged and produced high-quality colour plots (or even seductive video presentations) of the CFD simulations. Make sure that an elementary interpretation of the flow-field explains the fluid behaviour and that the trends of the flow analysis can be reconciled with a simple view of the flow. ¾ Be aware of the possible errors from bugs in the pre- and post-processors (interpolation, extrapolation, averaging, etc.). Make sure that the mean values of engineering parameters derived from the simulation are computed consistently (e.g. mass-average values, area-average values, timeaverage values). Calculation of local and mean engineering parameters with external postprocessing software may be inconsistent with the solution method of the code used (e.g. calculating shear stresses from the velocities, calculating shear stresses using nodal values instead of wall functions). Check that any test data used for comparison with the simulations is also computed in the same way as the data from the simulation. ¾ Consider whether the interpretation of the results, and any decisions made on the basis of them, are within the accuracy of the computation.
¾
6.2.5. Guidelines on documentation Keep good records of the simulation with clear documentation of assumptions, approximations, simplifications, geometry and data sources. ¾ Organise the documentation of the calculations so that another CFD expert can follow what has been done. ¾ Be aware that the level of documentation required depends strongly on the customer's requirements as defined in the problem definition.
¾
6.3. Training requirements for CFD users The growth in the use of CFD codes and the trend for them to become rich packages with lots of alternative modelling options, steadily increases the risk of user errors. This trend is reinforced by the ease of use of modern computer codes with simple graphical user interfaces making them available for inexperienced users. Although efforts taken to simplify the usage of CFD codes are most welcome, careful training with realistic exercises should still be considered as the starting point of any user’s CFD career. The theoretical part of the training should focus on fundamental modelling features, their underlying assumptions and their limitations. The same information is also a central part of a good user documentation. Unlike linear finite element stress analysis, CFD still requires expertly trained users for good results. In situations where non-experienced users have to be used, some restriction on their freedom to adjust critical parameters might be advisable, and they should be limited to simulations of routine types. Depending on the CFD software, additional training on grid generation may be advisable.
6.3.1. Guidelines on the training of CFD users ¾
A CFD user for non-routine applications should have good training and knowledge in classical fluid mechanics (including experimental, analytical and computational aspects), a broad understanding of numerical methods as applied to fluid flows, and detailed knowledge of the application being examined. This means that the user should be able to understand the limitations of the models used (e.g. turbulence, radiation, buoyancy driven flows, numerical models, etc.) and be aware of their relevance when being applied to a particular problem.
© ERCOFTAC 2000
page 35 of 94
The training and education requirement for more routine applications can be less stringent, provided that clear guidelines or procedures have been established for the use of the code being used. An example of a routine application would be the simulation of a standard component in a design environment where many previous designs have been calculated and only relatively small changes in geometries and boundaries conditions occur. ¾ In both routine and non-routine applications, training on the use of the specific CFD code with the solution of realistic exercises is needed. ¾ Training should be repeated regularly, and the user should frequently read the release notes and manuals issued by the code supplier, as these are a good source of educational information.
¾
© ERCOFTAC 2000
page 36 of 94
7. Code errors The success of a code leads to it becoming used more widely. As the user-base expands, there are increasing demands for more options and the code becomes more and more complex. As it can deal with more difficult problems, there is again an expansion of its use. In the end it is inevitable that code errors will be discovered by the many users who outnumber the developers by an order of magnitude and have a much wider range of applications and test cases than the code developers themselves. The size and complexity of large CFD software packages inevitably mean that code errors (bugs) may still be present in the software even if it has been in use and development for many years. For an illuminating description of some errors that were unobserved for many years, see Roache [1998]. He quotes an example where a Fortran code used the symbol MACH for Mach number, which is automatically interpreted as an integer. The mistake went undiscovered as all of the validation and verification cases were for integer values of the Mach number! The painstaking but straightforward process of verification provides a means of checking that the code faithfully reproduces the model approximations incorporated in the algorithms being programmed. The main problem associated with code verification is that the accuracy of a code can never be formally demonstrated for all possible conditions and applications, and for all possible combinations of valid code input options. In fact it can never be proven that a code is correct in that sense, as at any time a new bug may be found.
7.1. Guidelines for the code developer and vendor ¾ ¾ ¾
¾
¾ ¾
¾
The code developer or code vendor needs to demonstrate that he has applied stringent methods of quality control of the software development and maintenance. Verification of the code is to be carried out by the code vendor or developer, and he should provide the necessary information on the verification process for the user. The code developer or code vendor should maintain and publish a databank of verification test cases that are used for testing. The cases should include simple code verifications tests (e.g. that solutions are independent of co-ordinate systems). The code developer or code vendor should provide documentary evidence of the verification tests that the software has undergone, which should include clearly details of the code options which are used during testing. For all new versions of the code, new compilers and all new machines a standard set of verification test cases should be repeated. Code vendors and code developers should supply a list of known bugs and errors in each version of the code (hot-line, password secured web-page). This list should demonstrate that the number of bugs reduces as the code matures. The code developers should try to include warning notices and guidance for the user in the output. For example, when basic rules on grid generation (expansion ratios, skew, etc.) are being broken, when important specific default options are being overruled by the code input data or when the near wall grid is inconsistent with the turbulence modelling.
7.2. Guidelines for the code user Read the release notes issued by the code supplier, as these are a good source of information on changes to the code to remove existing bugs. ¾ Recognise that codes can only be validated and verified for a class of problems involving specific variables. If the application involves an area where the code is not fully verified there is more risk of code errors. ¾ Set up a suite of test cases that can be run on new code releases to provide an independent check on the code and to highlight any changes between releases (for example in default parameters). ¾ Check that the release version of the code is compatibile with any user defined routines, if applicable. ¾
© ERCOFTAC 2000
page 37 of 94
When a code error is suspected, communicate this to the code vendor or developer as soon as possible, especially if no list of known bugs has been published. Other users may then profit from this experience. If the bug is well-known a solution or work-around may already be available. ¾ In any communications with the code developer or code vendor about a suspected program error, provide a short concise description of the problem and all the necessary input data files so that the error can be reproduced. In cases where commercial sensitivity precludes this, special arrangements will need to be made. ¾ Contact the code developer if the meaning and definition of variables in the manuals are unclear and if support is needed in the implementation of special features. ¾
© ERCOFTAC 2000
page 38 of 94
8. Validation and sensitivity tests of CFD models 8.1. Selection of test cases for validation Validation is the process of ensuring that the CFD model being solved is a good representation of the real world. The need for careful validation should be clear from the earlier discussions, as the simplifying assumptions made in the CFD models mean that there is no guarantee that the results will agree with experiments. CFD validation is undertaken by systematic comparison of CFD results to reliable and appropriate experimental data or algebraic solutions, when available. This is a highly time-consuming process, as it requires careful selection of suitable experimental data, assessment of its accuracy, careful and systematic simulation of the test case and analysis of the results. Verification is limited to simple test cases likely to be performed by code developers. In industrial flows verification, validation and calibration, as defined in chapter 2, cannot be entirely separated.
8.1.1. Guidelines on validation ¾
¾
¾
¾ ¾
¾
¾
Before using a CFD code incorporating a specific model for a particular application, validate it against test data for a similar application with similar flow structures and flow physics. A design application, where accurate performance and flow-field data is needed, will require more detailed validation than an application involving a preliminary investigation of a flow. Validation of engineering cases should be aimed to check the most prominent and relevant flow features, rather than to capture every detail in the calculation. For the improved validation of complex flows, consider the taxonomy of the flow and validate the models being used on relevant simpler flows. For example, the validation of a code and turbulence model for a radial pump application should include not only a similar radial pump test case, but also test cases demonstrating the typical flow structures found in the pump, such as decelerating flow in a cascade of blades, flow in a rotating channel, secondary flow in bends, and leading and trailing edge flows. Simple test cases used in validation will allow finer grids than those used in real applications so that results of calculations at the typical grid resolution expected in the application should also be investigated. Carefully check the quality and accuracy of any measurements used as a validation test case, to provide confidence in the experiment and the quality of the experimental data. The experiments should attempt to capture the relevant boundary conditions and initial conditions where appropriate. If the experimental results are found to be sensitive to the boundary conditions and experimental set-up then this should be reported. The experimental data should include measurements of the quantity which is of interest and needed for CFD validation. This is best done by close collaboration between experimentalists and CFD analysts. If collaboration between experimentalists and CFD analysis is not possible care should be taken that the physically relevant quantities have been measured and that the quality and accuracy of the data is known. The validation data set should be independent of any data set used to construct the CFD model. Carry out the validation simulations in accordance with the guidelines outlined in this document. In particular, one of the validation calculations should use about the same mesh size as will be used in the real application. A simple test of validity should be carried out by examining elementary flow features of physical relevance. The global mass balance should be checked, the total temperature should not vary if no heat or work is being added to the flow, and in areas away from solid boundaries the total pressure may often be assumed to be constant. Predictions of the entropy in a flow field can also be a helpful validation check. When using CFD for a design application, examine the flow of a well understood and well tested baseline design comprising the same important flow features as a means of calibrating the expected accuracy of the code to this application. The reliability of predictions of a new design can then be assessed by taking into account the calibration of the code with the baseline design. Such a calibration allows a more accurate estimate of the predicted performance of the new design to be made by consideration only of the predicted differences and trends between the new design and the baseline design. In this way CFD can be used to rank different designs with similar flow structures.
© ERCOFTAC 2000
page 39 of 94
¾
Document the validation efforts. In any CFD calculations which are carried out without undergoing a thorough validation process (which is sometimes necessary as there may be no suitable measurements, or the time available allows no alternative), then attention should be drawn in the documentation of the results to the speculative nature of the simulation, with comments on why the user believes that the results are correct.
8.2. Design and use of sensitivity tests The objective of a sensitivity test is to allow the user to examine the effect on the predictions of the variation of any particular parameter in the simulation. In the first instance this technique can be used to examine the predicted flow over a range of certain flow parameters and boundary conditions whose exact value may be ill-defined, such as the turbulence level at the inlet to the flow domain or the inlet flow angle to a turbomachinery blade row. In addition such sensitivity analyses can be carried out to examine the effect of various model parameters, empirical constants, physical properties, code default options, discretisation schemes, grid types and mesh sizes.
8.2.1. Guidelines on use of sensitivity tests Before using a CFD code and model for any specific application, carry out and document tests to determine the sensitivity of the results to the key parameter settings. The key issues to be examined are sensitivity to the grid, discretisation scheme, turbulence model and boundary conditions. Sensitivity analyses can be very time-consuming and this means that in complex cases they may be limited in scope. Because of this, it may be satisfactory to carry out sensitivity tests not on each particular application, but for similar classes of applications requiring the use of the same models and the same code settings. ¾ Base the sensitivity studies on experience of the fluid dynamics of the case under consideration. Rather than changing parameters arbitrarily, physical understanding of the flow should be used whenever possible to identify the parameters most likely to exhibit sensitivity in the results. ¾ Document clearly the results of a sensitivity analysis, to show to what extent the simulation results depend strongly or not at all on the parameter settings. Unsubstantiated assertions as to the insensitivity of the results to certain factors should not be made, but the argument should be demonstrated by reference to quantitative examples.
¾
© ERCOFTAC 2000
page 40 of 94
9. Examples of application of best practice guidelines 9.1. Introduction to the test cases At the workshop held in Zurich (21 and 22 June 1999), a small working group decided on the following conditions with regard to the test cases that should be included in the ERCOFTAC Best Practice Guidelines: •
The cases should be selected to show each of following important physical effects: strong swirl, heat transfer, flow in rotating components (turbomachinery), compressible flow, and unsteady flow.
•
The cases should individually show two or three of the following effects of relevance to the best practice guidelines: numerical errors, grid dependency, spatial discretisation, temporal discretisation convergence, turbulence model, wall function, application uncertainty (boundary conditions and geometry).
•
The test cases should include some comparison with experimental data or to an analytic solution.
Various code vendors and code developers offered a selection of test cases for inclusion in these guidelines. From discussions at the workshop, and in accordance with the conditions outlined above, a selection from the various test cases was made. It is important to note that the presence of an example from a particular code vendor in a particular application area is not intended to provide any endorsement of the acceptability of the code for this particular purpose. Likewise, the absence of an example from a particular code vendor in a certain application area provides no statement on the unsuitability of the code for this particular application. The literature references which are of direct relevance to the specific test case are given at the end of the test case, but references of a more general nature are given in chapter 13.
9.2. Test case A : 2-D transient scalar bubble convecting at 45 degrees Test case calculated with the finite volume Code STAR supplied by Computational Dynamics, London.
9.2.1. Introduction The purpose of this example is to demonstrate the effects of the order of temporal discretisation. A first-order scheme (Euler Implicit) is compared with a second-order scheme (Crank-Nicholson). This example substantiates the best practices described in section 3.6. Additionally, the effects of the order of spatial discretisation are illustrated on the same case. A firstorder scheme is compared with the second-order scheme above using Crank-Nicholson temporal discretisation. This example illustrates the points in section 3.5. A perfectly circular scalar "bubble" is initially positioned at the bottom left of a square domain, in a fixed constant velocity field directed at 45 degrees towards the top right of the domain. The bubble circular shape should be preserved as it convects across the mesh. The ability of the numerical scheme to preserve this perfectly circular shape is a function on the temporal and spatial order of discretisation. High temporal accuracy is important for any applications which rely on time-dependent propagation phenomena, such as pressure pulses, and especially when the propagation occurs across mesh lines, as with vortical flow structures. Therefore, applications such as acoustic analyses or large-eddy simulation (LES) rely on higher-order capabilities or else the discretisation practice will introduce unphysical effects (phase change, or distorted eddies) which may be indistinguishable from the true physics. High spatial accuracy is necessary to preserve the sharpness of bubble interface. To summarise, three cases are presented; Case 1 - First order temporal, second order spatial discretisation Case 2 - Second order temporal, first order spatial discretisation Case 3 - Second order temporal, second order spatial discretisation
© ERCOFTAC 2000
page 41 of 94
9.2.2. Geometry and boundary conditions The bubble is defined by setting a scalar value to unity inside the bubble and zero outside at time = 0.0. The bubble diameter is 0.5m. The 2-dimensional domain extends 3m x 3m, the left and bottom faces of which are defined as inlet boundaries at u = v = 0.5m/s. The top and right faces are defined as outlets.
9.2.3. Grid and time step The mesh contains 300 x 300 cells, and is advanced at a time step size of 0.005 seconds. This corresponds to a calculation CFL number (based on convection time) of 0.25.
9.2.4. Features of the simulation For simplicity, only the scalar transport is solved for. The scalar is passive and the flow domain is unchanging in time (u = v = 0.5m/s).
9.2.5. Results The results of the simulations with different temporal and spatial resolutions are shown in Figure 1. In Case 1, (row 1) it is noticeable that the shape of the bubble distorts; it becomes elongated on the propagation direction, and the interface becomes smeared in this direction. This is a feature of lowerorder temporal discretisation. In contrast, Case 2 (row 2) preserves the shape of the bubble (higher-order temporal discretisation) by smearing the bubble interface. This is a feature of lower-order spatial discretisation. Finally, Case 3 (row 3) retains the shape and the sharpness of the bubble.
9.2.6. Conclusions The importance of higher order temporal and spatial discretisation for preservation of the phase and sharpness of a propagating flow features has been demonstrated. The first-order method distorts and smears, whereas the second-order method preserves phase and sharpness. This example also illustrates the strong coupling between spatial and temporal discretisation. For accurate solutions in time dependent calculations it is necessary to use higher-order discretisation in both time and space.
Figure 1: Evolution of the scalar bubble for different temporal and spatial discretisations. Upper row: 1st order temporal, 2nd order spatial discretisation. Middle row: 2nd order temporal, 1st order spatial discretisation. Lower row: 2nd order temporal, 2nd order spatial discretisation. © ERCOFTAC 2000
page 42 of 94
9.3. Test case B : T-Junction between main and auxiliary pipe Test case calculated with the finite volume Code ESTET Version 3.2 (Mattéi et al. [1993]) supplied by EDF.
9.3.1. Introduction This case is related to the isothermal swirling flow developing in a T-junction between a high Reynolds number main pipe and an auxiliary pipe at very small and zero mass flow rates. The motivation is to understand and model the hydraulic behaviour of various auxiliary lines connected to the primary circuit of Pressurised Water Reactors. Experimental data have been obtained on a hydraulic mock-up using flow visualisation and particle image velocimetry (PIV). They show an helicoidal (corkscrew) flow extending along the auxiliary pipe. This case illustrates the importance of the choice of the turbulence model k-ε eddy viscosity model and a Reynolds stress turbulence model) for computing this type of swirling flow and a grid dependency study shows the determining influence of the mesh refinement on the quality of the results.
9.3.2. Geometry and boundary conditions A sketch of the experimental set-up is shown in figure 1. The cross-section surface is Sm = 1.303 10-2 m2. The dead leg diameter D is 100 mm. Its dimensionless length H/D is 20. The origin is located at the main pipe centre. Auxiliary pipe (inlet) R = 50 mm
VA 2000 mm VM
Main pipe (inlet)
outlet
500 mm
1500 mm
100 mm
R = 130 mm 100 mm
Figure 1: Sketch of the mock-up.
Main pipe characteristics Flow-rate (m3/s)
0.120
Bulk velocity (m/s)
VM
9.2
Auxiliary pipe characteristics Bulk velocity (m/s) VA / VM
VA
0.092
0.046
0.023
0.
1%
0.5 %
0.25 %
0.%
Table 1: Flow rates At the inlets (main pipe and auxiliary pipe), Dirichlet conditions are used for all variables. In the main pipe, the incoming flow is fully developed. Hence, the mean velocity and turbulent quantities were obtained from preliminary periodic pipe computations with the same bulk velocity and the same crosssection. For the auxiliary pipe, due to the very small flow-rates, it is not necessary to describe so accu© ERCOFTAC 2000
page 43 of 94
rately the incoming flow. The mean velocity was taken constant and turbulent quantities were estimated from usual experimental laws for the friction velocity. In fully developed pipe flows with zeroroughness, and for Reynolds numbers Re (based on the bulk velocity VA and the hydraulic diameter D) ranging between 5,000 and 30,000, the friction velocity uτ can be determined from uτ = VA ((0.3164/Re0.25)/8)0.5. The turbulent kinetic energy k can then be estimated from uτ as k = uτ2/0.3 and the dissipation as ε = uτ3/(0.42 D 0.1). As for Reynolds stresses, the diagonal components are taken equal to 2/3 k whereas the extra-diagonal components are set to zero (assuming isotropic turbulence). Flow rates and corresponding bulk velocities are given in Table 1. Note that the bulk flow rate remains constant and the auxiliary flow rate varies from zero to 1% of the bulk flow rate. At the outlet (main pipe), zero gradient conditions were applied for all transported variables. Wall functions were used at solid walls.
9.3.3. Grid The equations are discretized on a 3D semi-staggered grid. Velocity is located at the vertices and pressure is cell-centred. When the k-ε model is used, the turbulent variables are located at the vertices, whereas they are cell-centred with the Reynolds stress model. The mesh is single block structured in Cartesian orthogonal coordinates. Curved boundaries are represented by special features of the code which enables the treatment of slanted boundary cells (prisms, tetrahedrons, etc.). A grid dependence study has been carried out with meshes of 100 000, 400 000, and 1 500 000 nodes for each turbulence model and each auxiliary mass flow rate considered (except for zero mass flow rate, for which only the finest grid has been used). Aspect ratios (ratio of the longest to the smallest edge in each cell) have been kept approximately between 1 and 2.
9.3.4. Features of the simulation The numerical techniques are based on finite difference and finite volume discretizations. The algorithm for the solution of Navier-Stokes equations relies on a segregated velocity-pressure formulation. The advective terms are treated by a method of characteristics. The trajectory is approximated by a second order Runge Kutta scheme with a third order 3D Hermitian polynomial for interpolation. Diffusion with explicit and implicit source terms for dynamic variables is solved implicitly. For the computation of the velocity components, a third step is required in order to prescribe the mass conservation, leading to a Poisson equation for the pressure increment. In order to avoid non physical oscillations of the pressure field and the associated difficulties in obtaining a converged solution, a variant of the Rhie and Chow interpolation is used (Deutsch et al., [1997]). The turbulence models that have been used are: - the standard high Reynolds number k-ε model proposed by Launder and Spalding [1974] - a high Reynolds number Reynolds stress model (second moment closure) with isotropic dissipation, turbulent self-transport modelled by the usual gradient diffusion following Daly and Harlow [1970] and the pressure-strain correlation term recommended by Launder [1989], i.e. consisting of the sum of Rotta's return to isotropy and isotropisation of production terms. The fluid is incompressible and its viscosity is ν = 1.03 10-6 m2/s. The corresponding Reynolds number based on the main pipe bulk velocity and the diameter of the auxiliary pipe is 895 000. Default initial conditions implemented in the code have been used for each simulation. The flow investigated being unsteady, the convergence had to be determined from time averages. Hence, once the vortical motion was fully established, time average values of the variables were estimated in the dead leg. The vortex penetration was then determined by visualising secondary velocities. The calculations were stopped once the value for the mean vortex penetration was stabilised.
9.3.5. Results Figure 2 compares the length of vortex penetration obtained by numerical simulation with experimental measurements. The length of the vortex penetration has been determined numerically by visualising the time-averaged velocity field obtained in different planes of the dead leg. This method is consistent © ERCOFTAC 2000
page 44 of 94
with the way the experiments were carried out (visualisation of averaged results from one thousand instantaneous velocity fields obtained using PIV system). It is also important to note that experiments and calculations were conducted by the same team, thus making comparisons easier and less subject to interpretation errors. The results for three different flow rates are presented. Significant improvement is obtained when refining the mesh. The prediction of the vortex penetration in the dead leg using Reynolds stress model on the finest mesh is in very good agreement with the experiments. The differences between medium and fine meshes illustrate the necessity to use very fine grids with the Reynolds stress model. It might be necessary to perform yet another simulation on an even finer mesh to achieve full grid independence. With the k-ε model, on the contrary, the computed vortex penetration is not satisfactory, even on the finest mesh. Moreover, the very small difference between the results obtained on medium and fine meshes indicates that the use of even more refined meshes would be unlikely to bring significant improvement.
H/D 10
Experimental SMC fine mesh SMC medium mesh SMC coarse mesh
8 6
k-ε fine mesh k-ε medium mesh k-ε coarse mesh
4 2 0,2
0,4
0,6
0,8
1
VA/VM (%)
Figure 2: Length of vortex penetration for different flow rates in the auxiliary pipe. Comparison between Reynolds stress model (SMC) and k-ε model. Coarse, medium and fine meshes have been used. Figure 3 presents the rotation velocity at different sections in the auxiliary pipe obtained by both numerical and experimental simulations. According to tests by Robert [1992], the strength of the rotation velocity in the auxiliary line is very sensitive to geometric details of the T-junction. Small differences between the mock-up geometry and the numerical mesh could generate differences on the rotation velocity. To validate independently the vortex decay in the auxiliary line, numerical and experimental results are compared using dimensionless quantities. The rotation velocities obtained from numerical and experimental simulations are scaled by the corresponding rotation velocities near to the T-junction taken at H/D = 3. Using Reynolds stress model, the vortex decay is very close to the experimental result, whereas it is overestimated with the k-ε model.
© ERCOFTAC 2000
page 45 of 94
dimensionless rotation velocity 1.00
SMC k-ε
0.75
Experimental 0.50
0.25
0.00 2
4
6
8
H/D
10
Figure 3: Maximum of rotation velocity scaled by the corresponding value at the section H/D=3. Comparison between Reynolds stress model (SMC), k-ε results and experimental data. Figures 4 and 5 help to visualise the flow features. They represent the vortex penetration in the dead leg for VA/VM = 0.25 % (k-ε model and RSTM respectively). The positive vertical velocity component is displayed in black and the negative in grey: black trajectories represent fluid particles coming from the main pipe and grey trajectories represent the fluid particles coming from the top of the auxiliary pipe. Velocity vectors are plotted at three different sections of the auxiliary pipe H/D=2, 5 and 7 from the results of fine, medium and coarse meshes.
fine mesh
Va
medium mesh
Coarse mesh
H/D= 9
H/D= 7
H/D= 5
H/D=3
Vm
Figure 4: k-ε model prediction of vortex penetration in the dead leg: VA/VM = 0.25 %. © ERCOFTAC 2000
page 46 of 94
fine mesh
Va
medium mesh
Coarse mesh
H/D= 9
H/D= 7
H/D= 5
H/D=3
Vm
Figure 5 : Reynolds stress model prediction of vortex penetration in the dead leg: VA/VM = 0.25 %.
9.3.6. Conclusions Numerical calculations performed on the T-junction with coarse (100 000), medium (400 000) and fine (1 500 000) meshes were used to test the influence of the grid. The results illustrate that a thorough mesh refinement is necessary before trying to draw conclusions from any numerical simulation. Indeed, on too coarse a mesh, results from the Reynolds stress model might turn out to be significantly worse than those obtained with a k-ε model. The numerical results obtained with the Reynolds stress turbulence model on the finest mesh confirm the presence of the swirling flow structure observed experimentally in the dead leg. On the other hand, with the k-ε eddy viscosity model, the vortical penetration in the dead leg is not satisfactorily reproduced. This demonstrates that Reynolds stress modelling might be a very interesting alternative to less sophisticated turbulence models, such as the k-ε model, especially to capture swirling flows.
9.3.7. References Mattéi J.D., Simonin O., "Logiciel (1993), “Reference manual : ESTET - Manuel théorique de la version 3.1 - Tome 1 - Modélisation physiques" , EDF report HE-44/92/038/A. Deutsch E., Montanari P., Mallez C. (1997), "Isothermal study of the flow at the junction between an auxiliary line and primary circuit of pressurised water reactors", Journal of hydraulic Research Vol. 35, No. 6, 799-811. Robert M. (1992), "Corkscrew flow pattern in piping system dead legs", NURETH-5, September, Salt Lake City, USA.
9.4. Test case C : Natural convection flow in a square cavity Test case calculated with CASTEM2000 v.99 developed by CEA.
9.4.1. Introduction The test case chosen is the computation of the steady-state, natural convection flow of a compressible fluid (air) in a square cavity, with differentially heated vertical walls. It is customary to assume incompressible flow (the local Mach numbers are extremely low), constant flow properties and to model buoyancy forces with the so-called Boussinesq approximation (Gray and Giorgini, [1976]). Under this
© ERCOFTAC 2000
page 47 of 94
hypothesis, density is assumed constant in the flow equations except in the gravity force term, where small variations due to temperature differences are modelled in the form g(ρ-ρ0 )/ρ0 ≈ g∆T/T0 . Most natural convection solvers used in the industry are indeed based on this model, but these assumptions are only valid for small temperature differences, and may be the source of substantial errors in heat transfer coefficients. For large temperature differences, compressibility effects come into play even if the local Mach numbers remain extremely small so that it becomes necessary to use a compressible model to compute the flow. This test case will illustrate the following guidelines: • •
the effect of the flow model on the heat transfer results (modelling errors); the grid-sensitivity of results (numerical errors).
9.4.2. Geometry and boundary conditions The geometry is a square cavity, of length L, with thermally insulated horizontal walls, and differentially heated vertical walls. The right wall is maintained at a “cold” temperature Tc, the left wall at a “hot” temperature Th. No slip boundary conditions are assumed on all solid walls.
9.4.3. Grids The grids used consist of a series of structured quadrilateral meshes, slightly clustered near the solid boundaries where stronger gradients occur. In order to conduct a grid-convergence study, the following meshes were used: 20x20, 40x40, 80x80, and 160x160.
9.4.4. Features of the simulation The cavity is filled with air, assumed to behave as a perfect gas, with Prandtl number Pr=0.71 and the ratio of specific heats γ equal to 1.4. The viscosity µ and the thermal conductivity λ are assumed to depend only on temperature, according to Sutherland’s law,
µ(T ) T = ∗ µ∗ T
3/ 2
T∗ + S T +S
with T* =273 K, S=110.5 K and µ* =1.68 10-5 kg/m/s. We shall consider here a Rayleigh number Ra=106, for which the flow is known to be laminar. We recall that the Rayleigh number is defined as: 2
gρ0 (Th − Tc )L3 Ra = Pr T0µ 2 Uniform initial conditions are assumed, with air at atmospheric pressure Po=101325 Pa, and a temperature To= (Th + Tc)/2 = 600 K. Two temperature ratios ε = (Th - Tc) / (Th + Tc) will be considered, ε = 0.01 for which the Boussinesq flow model is valid, and ε = 0.2 for which a compressible flow model is required. With these definitions, expressions for the hot and cold temperatures Th and Tc are respectively Th = To (1 + ε) and Tc = To (1 - ε). The numerical scheme used for the computations is an elliptic finite element solver based on an asymptotic approximation of the Navier-Stokes equations at low Mach numbers in which the pressure is split into an average time-dependent thermo-dynamic pressure which obeys the perfect gas law, and a hydrodynamic fluctuation pressure which obeys a Poisson pressure equation just as for incompressible flow. The discretisation is based on bilinear quadrilateral elements (Galerkin method), and the solutions are obtained by time-marching to steady state using a time-accurate scheme.
9.4.5. Numerical results Small temperature difference case The main objective of the simulation is to compute the heat transfer to the walls. Hence, the figures below show distribution plots of the local Nusselt number along the vertical walls. The Nusselt number is a non-dimensional number defined as:
Nu ( y ) =
© ERCOFTAC 2000
L ∂T λ λO (Th − Tc ) ∂ x
wall
page 48 of 94
In this first computation, we perform a grid convergence study for the Ra = 106 and ε = 0.01 case. Table 1 shows the evolution of the average, maximum and minimum Nusselt number computed on the right wall as the mesh is progressively refined. In all cases, steady-state was obtained with a drop of the momentum residual of more than 10 orders of magnitude. A reference solution correspond to gridconverged Boussinesq solutions, which assume incompressible flow and constant flow properties. In that case, a symmetric solution is obtained so that identical heat transfer distributions are obtained on the walls.
ε=0.01
20x20
40x40
80x80
160x160
8.379 (5%) 1.018 (4%) 15.046 (14%)
8.760 (0.7%) 0.982 (0.3%) 16.999 (3%)
8.822 (<0.1%) 0.981 (0.2%) 17.433 (0.6%)
8.832 (<0.1%) 0.981 (<0.2%) 17.510 (<0.2%)
Nu_min Nu_max
Ref. Incomp. Sol. 8.826 0.979 17.536
Table 1: Heat transfer results for ε = 0.01 case; grid convergence study and comparison with reference solutions; the figures in brackets represent the relative error with respect to Le Quéré et al. [1985]. The grid convergence study shows that the solution varies sensibly up to the 80x80 mesh, and very little when the mesh is refined further (see Fig. 1). Thus, one can assume that the solution on the 160x160 mesh is nearly grid converged. For that solution, the results are in close agreement (relative error in heat transfer of the order of 0.2%) with reference solutions found in the literature, which correspond to incompressible flow models with the Boussinesq approximation. Thus, discrepancies between the solutions computed on the coarse meshes and the reference solutions can be attributed to numerical errors, stemming from the use of too coarse meshes or excessively dissipative schemes. On the other hand, on fine enough meshes, the incompressible flow models and the more general low Mach number compressible flow models yield extremely close results (and within the band of reference solutions). Thus, for small temperature differences, the modelling errors (errors due to the choice of flow model) are negligible compared to numerical errors.
“Large” temperature difference case
In this case, we consider a larger temperature difference, corresponding to ε = 0.2. Note that this parameter is still quite small compared to 1, and that it would be tempting to use an incompressible flow model with the Boussinesq approximation. However, we shall see that even for this moderately large temperature difference, and in spite of the extremely low Mach numbers, a compressible flow formulation is required to solve the problem. Table 2 shows the grid converged heat transfer results obtained on the fine mesh (160x160) (validated against other compressible solvers by Paillère et al., [1999]), alongside the reference incompressible flow results. Note that the average Nusselt number has remained almost constant (as was observed in Chenoweth et al., [1986], Le Quéré et al., [1992] and Paillère et al., [1999]). However, the minimum and maximum Nusselt numbers are substantially different, as may be observed from the distribution plots in Fig. 2, which display a noticeable non-symmetric character attributed to compressibility effects. The maximum local Mach number in the flow was found to be of the order of 3.4 x10-4 . Relative errors in heat transfer prediction between the compressible and the incompressible flow models are of the order of 5%. Thus, for large temperature differences, the incompressible flow model with Boussinesq approximation and the low Mach number compressible flow model yield substantially different results. The difference can be attributed to modelling errors, as it is known (but unfortunately often forgotten) that the former model is only valid for small temperature differences in Gray et al. [1976].
© ERCOFTAC 2000
page 49 of 94
ε=0.2
160x160
Ref. Incomp. Sol.
8.817
8.826 (0.1%) Nu_min 0.930 0.979 (5%) Nu_max 16.869 17.536 (4%) Table 2: Heat transfer results for ε = 0.2 case ; the figures in brackets represent the relative error between the reference incompressible solution (Le Quéré et al., [1985]) and the ‘grid converged’ compressible solution
9.4.6. Conclusions This example illustrates two sources of error that can affect a CFD simulation, numerical error due to insufficient grid refinement and modelling error due to inappropriate choice of flow model. When analysing flow with heat transfer and very low Mach numbers, care must be used in the choice of flow model (i.e extremely low Mach does not imply incompressible flow), and grid-sensitivity assessed by comparing results on grids of different size.
9.4.7. References Chenoweth, D. and Paolucci, S. (1986), “Natural Convection in an Enclosed Vertical Layer with Large Temperature Differences”, J. Fluid Mech., Vol. 169, pp. 173-210. Gray, D.D. and Giorgini, A., (1976) “The Validity of the Boussinesq Approximation for Liquids and Gases”, Int. J. Heat Mass Transfer, Vol. 15, pp. 545-551, 1976 Le Quéré, P., Masson, R. and Perrot, P. (1992), “A Chebyshev Collocation Algorithm for 2D NonBoussinesq Convection”, J. Comput. Phys., 57, pp. 320-335. Le Quéré, P. and Alziary de Roquefort, T. (1985), “Computation of Natural Convection in Twodimensional Cavities with Chebyshev Polynomials”, J. Comput. Phys., Vol. 57, pp.210-228. Paillère, H. and Magnaud, J.P. (1998), “A Finite Element Elliptic Flow Solver for Low Mach Number Compressible Flows”, pp. 419-424, Proc. 10th Int. Conf. on Finite Elements in Fluids, Tucson, Arizona, January 5-8, 1998. Paillère, H., Viozat, C., Kumbaro, A. and Toumi, I. (1999), “Comparison of Low Mach Number Models for Natural Convection Problems”, EUROTHERM Seminar No. 63 ‘Single and Two-Phase Natural Circulation`, Genoa, Italy, 6-8 September 1999 20
15
Nu(y)
10
right wall 20x20 right wall 40x40 right wall 80x80 right wall 160x160
5
0
0
0.25
0.5
0.75
1
y/L
Figure 1: Grid convergence study for ε = 0.01 case.
© ERCOFTAC 2000
page 50 of 94
20
15
Nu(y)
10 epsi=0.2, right wall epsi=0.2, left wall epsi=0.01, right wall epsi=0.01, left wall
5
0
0
0.25
0.5
0.75
1
y/L
Figure 2: Effect of temperature difference on heat transfer results
9.5. Test case D : Sudden Pipe Expansion Test case calculated with finite element code FEAT, Version 3.1.4 (FEAT User Guide, 1997), supplied by M. Rabbitt of British Energy.
9.5.1. Introduction The test case has been included to provide information about heat transfer at flow re-attachment, and to compare predictions with mean flow and turbulence experimental data in a spreading shear layer. The focus of the presented work is the sensitivity of the results to mesh errors and the accuracy of the turbulence model applied for global predictions of the flow-field and heat transfer. Results are presented for the k-ε turbulence model with wall functions. In addition, a modification to the turbulence model is suggested to improve prediction of heat transfer in a region of flow re-attachment. The geometry of the test problem is a sudden pipe expansion, as illustrated in Figure 1. Turbulent flow approaching the expansion from the inlet pipe separates at the step corner forming a shear layer which finally attaches to the larger pipe wall at a distance downstream of the expansion. A recirculating region of flow is contained between the shear layer, step and larger pipe wall and a secondary recirculation is set up in the corner. In the heat transfer experiments, a constant heat flux was applied to the larger pipe side wall, the other walls being insulated, and the inlet flow was at uniform temperature. The experimental data of Szczepura [1985], and Baughn et. al [1984] were used as control data for the velocity and thermal fields respectively.
9.5.2. Geometry and boundary conditions The problem consists of calculating the velocity and turbulence fields at a Reynolds’ number of 892000, using inlet boundary conditions specified at x=-0.25 in Figure 1. Having obtained the flow field, the thermal field is calculated, assuming a Prandtl number of 0.7. The simulation has been carried out assuming a two-dimensional rotational symmetric flow field. Figures 3 show the variation at inlet for the velocity, turbulence energy and turbulence viscosity (µeff=k1/2l, l is the length scale Cµk3/2/ε). This corresponds to a Reynolds number, based on the centreline speed and downstream pipe diameter, of 892000.
© ERCOFTAC 2000
page 51 of 94
Once a k-ε solution is obtained, the thermal energy equation is solved, applying a constant heat flux on the downstream pipe wall to provide a bulk temperature at the outlet of 1.0. The local bulk temperature therefore increases linearly along the pipe.
9.5.3. Grid The mesh design as a layout of blocks or patches is shown in Figure 2. The grid design is parametric, in that the parameter N controls the number of elements in both directions for each patch, as described in Table 1, which also indicates the squeeze factors chosen. The squeeze factor shown provides the ratio of the dimension of an element in one direction in relation to its neighbour. Additionally, to preserve the thickness of the wall element, a constant thickness layer of 0.005m is positioned next to the walls. Table 2 gives the number of elements in each mesh. Each element consists of eight nodes, with four corner and four mid-side nodes.
9.5.4. Features of the simulation The finite element code uses continuous nodal quadratic basis functions for the velocity, temperature, turbulence speed (q=k1/2) and turbulence frequency (f=ε/k). For pressure, continuous nodal linear and piecewise element constant basis functions are used together. The numerical scheme is formally fourth order accurate. The standard high-Reynolds number k-ε turbulence model is solved using the q-f formulation (Smith [1984]). Constant values are supplied for Cµ (0.09) c1ε (1.44) and c2ε (1.92). Wall functions are assumed. The effective viscosity required to provide the log-law of the wall for the tangential velocity is supplied to obtain all velocity components in the near wall element. The turbulence variables are obtained by assuming a constant shear stress normal to the wall. The length scale correlation l=Cµ1/4κy, where y is the distance normal to the wall, is reproduced in equilibrium boundary layers, but when conditions deviate from this transport effects are also incorporated. For the thermal energy equation, turbulent thermal fluxes are represented by an effective conductivity, kT=Cpµt/σT, where σT is the turbulent Prandtl number. In the wall element, similarity between boundary layers for velocity and temperature is assumed via the Jayatilleke function (Jayatilleke [1969]). Separate calculations are also undertaken with the assumption that the turbulent length scale can grow at a maximum rate of 0.25×y, where y is the wall distance. This improves the prediction of heat transfer in flow re-attachment regions, as shown in the results section. To solve, initial values are first provided for the turbulence speed, q, and length scale, l. These determine the eddy viscosity, to enable a single Newton-Raphson iteration of the Reynolds averaged Navier-Stokes (RANS) equations. The turbulence variables are then updated, by a series of NewtonRaphson iterations of the q and f equations. These values are relaxed with previous values to obtain an eddy viscosity to plug into the RANS equations. This process is cycled until the maximum change divided by the range for all variables, between successive cycles, is less than 0.003
9.5.5. Results Results showing the sensitivity of the simulation to mesh errors are analysed first, followed by results addressing turbulence modelling issues. Two reference points are chosen to demonstrate sensitivity of axial velocity and effective viscosity to mesh errors, point R1 located at x=0.05, r=0.25, and point R2, at x=1.0, R=0.1. The maximum Nusselt number was also found on each mesh. Table 3 shows values of the axial velocity at point R2, effective viscosity at point R1 and maximum Nusselt number on each of the meshes. It can be shown from these values that the error relative to the solution on the N=16 grid decays with the fourth order of the grid size (N4), which is consistent with the fourth order accuracy of the solution method being used. The wiggles produced from the high order numerical scheme were analysed to indicate mesh errors; the wiggle size is proportional to mesh Reynolds number. The maximum local error indicated on the N=4 mesh is 38%, it is 12% on the N=8 mesh and 3% on the N=16 mesh. Average errors are much smaller than this. The location of maximum error is near the corner.
© ERCOFTAC 2000
page 52 of 94
The experimental value for maximum Nusselt number on the wall is 932. Consequently, the standard modelling over-predicts the heat transfer by more than a factor of 2. Restricting the rate at which the length scale can grow away from a wall to 0.25×y provides a maximum Nusselt number of 1200. This is a significant improvement. The inclusion of such a length scale limiting device is advisable for all reattaching flows, and other simulations have shown that it can also be recommended for all buoyancy driven flows (e.g. Ince and Launder [1995]). Figures 4 and 5 show centre-line profiles of the turbulence energy and mean axial velocity in comparison with experimental values (calculated without the length scale limiter). The predictions give reasonable agreement with experiment for the gross features in the flow, apart from in the re-circulation region, where turbulence levels are significantly under-predicted, and near the centre-line edge of the shear layer, due to the poor prediction of the spreading rate of the shear layer. The global features of the flow prediction are not very sensitive to the use of the length scale limiter.
9.5.6. Conclusions The finite element scheme used in this example makes use of quadratic basis functions and this provides a fourth order accurate numerical scheme. The results demonstrate that the mesh errors decrease with the fourth order of the mesh spacing or grid size (N4). The standard k-ε turbulence model used seriously overestimates the heat transfer in a reattachment region by a factor of 2. A length scale limiting device is essential to provide a reasonable prediction of heat transfer in flow re-attachment regions. The predictions give reasonable agreement with experiment for the gross features in the flow, apart from in the re-circulation region, where turbulence levels are significantly under-predicted, and near the centre-line edge of the shear layer, due to the poor prediction of the spreading rate of the shear layer.
9.5.7. References FEAT User Guide Version 3.1.0 (1997), Nuclear Technology Branch, British Energy, Barnett Way, Barnwood, Gloucester, GL4 3RS Szczepura, R.T. (1985), “Flow characteristics of an axisymmetric sudden pipe expansion: results obtained from the turbulence studies rig. Part 1, Mean and turbulence velocity results”, CEGB Report TPRD/B/0702/N85. Baughn, J.W., Hoffman, M.A., Takahashi, R.K., and Launder, B.E. (1984), “Local heat transfer downstream of an abrupt expansion in a circular channel with constant heat flux”, J. Heat Transfer 106. Smith, R.M. (1984), “A practical method of two-equation turbulence modelling using finiteelements”, Int.J.Num.Meth.Fluids, 4, 1984, pp. 321-336. Jayatilleke, C.L.V. (1969), “The influence of Prandtl number and surface roughness on the resistance of the laminar sublayer to momentum and heat transfer”, Prog in Heat and Mass Transfer, 1, Pergamon Press. Ince, N.Z and Launder, B.E. (1995), “Three-dimensional and heat-loss effects on turbulent flow in a nominally two-dimensional cavity”, Int. J. Heat and Fluid Flow, 16, pp. 171-177.
Axial elements N
3N/2
4N
1,8,9, 10
2,3,4
5,6,7
Radial elements 3N/2 all
Axial squeeze factor 60 1/2 10 1,2,3, 4
5,6,7
8,9,1 0
Radial squeeze factor 10 1/20 1/40 3,6,9
1,2,5, 8
4,7,10
Table 1: Mesh design, showing patch details
N=4 534
N=6 1170
N=8 2052
N=10 3180
N=12 4554
N=14 6174
N=16 8040
Table 2: Number of elements in the various meshes
© ERCOFTAC 2000
page 53 of 94
Variable U1(R2) VEFF(R1) Nu(max)
N=4 0.921 0.138 2328.2
N=6 0.934 0.135 2341.0
N=8 0.939 0.124 2342.0
N=10 0.941 0.111 2341.3
N=12 0.941 0.108 2342.8
N=14 0.942 0.104 2343.8
N=16 0.942 0.104 2343.2
Table 3: Values of axial velocity, effective viscosity and maximum Nusselt number on the various meshes
X=6.0 Centre line
X=3.0
X=0.05
X=-0.25 Step corner, x=0.0, r=0.257 Figure 1: Pipe Expansion Test Problem Geometry
X=0.4
X=2.8
X=6
4
7
10
3
6
9
2
5
8
R=0.4 R=.257
1
R=0 Figure 2: Mesh Design, Layout of Patches
© ERCOFTAC 2000
page 54 of 94
axial velocity 1.20E+00 1.00E+00
radius
8.00E-01 6.00E-01 4.00E-01 2.00E-01 0.00E+00 0.00E+00
6.00E-02
1.20E-01
1.80E-01
2.40E-01
radial velocity (solid), k (dashed) 0.01 0.009 0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0 0
50
100
150
200
250
300
radius
turbulent viscosity 7.00E-04 6.00E-04 5.00E-04 4.00E-04 3.00E-04 2.00E-04 1.00E-04 0.00E+00 0.00E+0 0
5.00E-02
1.00E-01
1.50E-01
2.00E-01
2.50E-01
3.00E-01
radius
Figure 3: Inlet Profiles for Velocity,Turbulence Energy and Turbulence Viscosity
© ERCOFTAC 2000
page 55 of 94
Centre-Line k 0.04 0.035 0.03
k
0.025 Exp k
0.02
Pred k
0.015 0.01 0.005 0 0
1
2
3
Axial Location
Figure 4: Centre-Line Turbulence Energy, k, Prediction vs Experiment Centre-Line Axial Velocity 1.2
Axial Velocity
1 0.8 Exp U
0.6
Pred U
0.4 0.2 0 0
1
2
3
Axial Location
Figure 5: Centre-Line Axial Velocity, Prediction vs Experiment
9.6. Test case E : Transonic airfoil RAE2822 Test case calculated with the finite volume Code STAR supplied by Computational Dynamics, London.
9.6.1. Introduction The purpose of this example is further to illustrate the effects of spacial discretisation. Two formally second-order accurate schemes are demonstrated. The first uses a gradient-based advection scheme which can be termed 2nd-order Upwind Differencing, and the second is a similar scheme based on a more compressive monotone advection algorithm. The first scheme essentially uses a compact stencil as a result of the upwinding practice that it employs while the more compressive scheme uses a wider stencil and therefore it resembles a gradient-based central differencing scheme. This example is used further to illustrate spatial discretisation practices described in section 3.5 The RAE2822 airfoil, Mach 0.731 at 2.51 degrees incidence and a Reynolds Number of 6.5 x106, is a well known benchmark exhibiting supersonic expansion and shock compression. Predictions of the
© ERCOFTAC 2000
page 56 of 94
pressure coefficient at the surface of the airfoil, lift and drag are compared with experimental data( Cooke et al. [1979]). In such cases which contain rapid expansions (leading edge suction) and compression (suction surface shock), non-compressive schemes have the effect of smearing the high-gradient features. Therefore, the prediction of critical parameters, such as lift, are adversely affected. Compressive discretisation practices are necessary.
9.6.2. Geometry and boundary conditions The computational domain for 2-dimensional airfoil calculations needs to extend between 10 to 20 chords around the airfoil. This requirement is to negate the effect of 'upwash' (which is significant for 2D lifting bodies) at the far-field boundaries. Upstream conditions are defined as Mach 0.731 at 2.51 degrees to the airfoil chord.
9.6.3. Grid In Figure 1 the computational mesh is shown. It contains 12288 cells, using a C-grid construction with 192 cells wrapped around the airfoil (128 of which are on the surface) and 64 in normal to surface direction spaced appropriately to enable a low Reynolds number turbulence model to be used (see section 4.3.3). For this grid, y+ for the near-wall cells range between 1.5 and 5.0, the higher part of the range occurring around the suction peak.
9.6.4. Features of the simulation The simulations are performed as fully turbulent, compressible flow calculations. The boundary layers are captured using a low Reynolds number two-equation eddy-viscosity k-ε model. Convergence is judged according to the stabilisation of the lift and drag.
9.6.5. Results The figures 2 and 3 below show the Cp (Pressure coefficient) distribution on the airfoil surface, Mach number and near-surface mesh distributions. The compressive scheme captures well the high leading edge suction gradient and upper surface shock. In contrast the low-compression scheme positions the shock correctly, but fails to capture the sharpness of both the leading edge suction and the shock. The effect is to underpredict the upper surface suction, and reduce the predicted level of lift. In both cases, drag is overpredicted. The two-equation eddy-viscosity k-ε model is known to overpredict turbulence viscosity (section 4.4.1).
Non-compressive Compressive Experiment
CL(lift) CD(drag) 0.693 0.0279 0.733 0.0184 0.743 0.0127
9.6.6. Conclusions In this example, non-compressive and compressive disretisation practices have been employed to predict the flow around an RAE2822 airfoil, and results compared with experimental data. The compressive scheme has been shown to be necessary for capturing sharp flow gradients such as leading edge suction and shock compression.
9.6.7. References Cooke P., McDonald, M. & Firmin, M. (1979), "Airfoil RAE2822 - Pressure Distribution and Boundary Layer Wake Measurements", AGARD AR-138.
© ERCOFTAC 2000
page 57 of 94
1.26 1.17 1.08 0.99 0.90 0.81 0.72 0.63 0.54 0.45 0.36 0.27 0.18 0.09 0.00 Mach No.
Figure 1: Computational mesh of the airfoil simulation
-1.5
-1.0
Cp -0.5
0.0
0.5
1.0
1.5 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
chord
Figure 2: Pressure distribution on the surface of the airfoil for the compressive scheme.
© ERCOFTAC 2000
page 58 of 94
-1.5
-1.0
Cp -0.5
0.0
0.5
1.0
1.5 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
chord
Figure 3: Pressure distribution on the surface of the airfoil for the non-compressive scheme.
9.7. Test Case F : Engine Valve Test case calculated with FLUENT version 5.1 supplied by Fluent.
9.7.1. Introduction The fluid flow inside an automobile engine cylinder is important for combustion efficiency and engine performance. Understanding this flow is crucial in the design of the cylinder, the valve, and the inlet port. This example has been included to show the effects of boundary condition uncertainty, discretisation and convergence level in a realistic three-dimensional industrial geometry where important recirculation phenomena occur. The geometry and problem set-up for this case refers to the experimental study of a generic model engine valve (Chen et al, [1995]) in which a mixture of liquid phase turpentine and tetraline was used to represent the working fluid.
9.7.2. Geometry and boundary conditions The geometry consists of the cylinder connected to the inlet port by means of the valve passage, as shown in Figure 1. The port axis is offset from the cylinder axis by 4.00 mm in the x direction and 21.87 mm in the y direction, and it is elevated from the horizontal plane at an angle of 40 degrees. The cylinder diameter is 93.65 mm, the inlet port diameter is 46 mm, and the valve diameter is 43 mm. For more details about the geometry of the problem refer to Figure 1 of Chen et al. [1995]. For the inlet, the turbulence length scale was taken to be the diameter of the port (46 mm). In the absence of experimental measurements, the turbulence intensity was initially assumed to be 10%. This value was later reduced to 1% so that the effect of boundary condition uncertainty could be illustrated. The inlet velocity was calculated from the known mass flow rate (of 1.379 kg/s) and is 0.928 m/s. A zero gauge pressure boundary condition was used at the outlet. The wall boundaries have been given a no-slip condition.
9.7.3. Grid A tetrahedral grid with 176461 cells, was used. The concentration of cells is higher in the vicinity of the valve, due to the complexity of the flow in this region. The maximum skewness of the cells is 0.875. Skewness is defined as (Optimal-Cell-Size – Cell Size)/(Optimal-Cell-Size) where optimal-cell-size is the size of an equilateral cell with the same circumradius. According to this definition of skewness, a
© ERCOFTAC 2000
page 59 of 94
value of 0 indicates an equilateral cell (best) and a value of 1 indicates a completely degenerate cell (worst).
9.7.4. Features of the simulation The fluid used in the experimental work was a mixture of turpentine and tetraline. Constant values were used for the fluid properties: • •
ρ = 894 kg/m3 µ = 0.00152875 kg/(m·s)
Density Viscosity
The Reynolds number based on the diameter of the inlet port is Re = 24 970. The RNG k−ε turbulence model with standard wall functions was employed in the calculation. The first order discretisation scheme was used for momentum and turbulence equations. The calculation was then repeated with second order discretisation to illustrate sensitivity to the discretisation method. All the variables were initialised with a zero value prior to calculation. Results have been produced at two levels of convergence to illustrate its impact on accuracy.
9.7.5. Results For each of the best practice issues, figures have been produced in which FLUENT predictions can be compared with experimental data. The plots show the velocity component in z direction at two different values of the z co-ordinate. Plots are also shown for turbulent viscosity (although no experimental data is available for this quantity). Where ‘best practice’ is applied the predicted velocities are seen to be in good agreement with the measurements from Chen et al. [1995].
Z X Y Flow Through an Engine Inlet Valve Grid
Jul 26, 1999 FLUENT 5.1 (3d, segregated, rngke)
Figure 1: Grid of the computational domain. •
Discretisation Method
Calculations were performed using both first and second order discretisation. The predicted velocity profiles obtained in both cases are compared with experimental measurement in figures 2 and 3. It can be seen that the second order scheme yields results that are in slightly closer agreement with the experimental data points. However, for this application it would appear that the accuracy of the first order results might be sufficient for most engineering purposes.
© ERCOFTAC 2000
page 60 of 94
1st order 2nd order measurements
0.00e+00 -2.00e-01 -4.00e-01 -6.00e-01
Z Velocity [m/s]
-8.00e-01 -1.00e+00 -1.20e+00 -1.40e+00 -1.60e+00 -0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
Y Position [m]
Z-Velocity at z=+10 mm
Oct 06, 1999 FLUENT 5.2 (2d, segregated, lam)
Figure 2: Discretisation – Z Component of Velocity at z= +10 mm 1st order 2nd order measurements
1.50e+00
1.00e+00
5.00e-01
Z Velocity [m/s]
0.00e+00
-5.00e-01
-1.00e+00
-1.50e+00 -0.03 -0.02 -0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Y Position [m]
Z-Velocity at z=-5 mm
Oct 06, 1999 FLUENT 5.2 (2d, segregated, lam)
Figure 3: Discretisation – Z Component of Velocity at z = -5mm
© ERCOFTAC 2000
page 61 of 94
•
Convergence
The iterative calculation was monitored and stopped when the ‘maximum residual’ value had fallen to 1.0E-2 and overall mass imbalance was 1.4%. [‘residual’ is defined in this case as the error associated with the discretised equation, summed over all computational cells and normalised by the sum of the centre coefficients of the equation. The term ‘maximum’ refers to the ‘residual’ for the transport equation showing the largest ‘residual’ value.] The situation described above represents the ‘not converged’ result shown in figures 4 and 5. Further iterations were then performed to obtain the ‘fully converged’ result shown in figures 4 and 5. Full convergence was deemed to have occurred when there was no further noticeable change in the monitored variables, ‘maximum residual’ value was at 1.0E-5 and overall mass imbalance had fallen to 0.002%. Figures 4 and 5 clearly show the error that can result if the user relies on an insufficiently converged solution. not converged fully converged measurements
1.50
1.00
0.50
Z Velocity (m/s)
0.00
-0.50
-1.00
Y Z
-1.50 -0.040
X
Flow through an Engine Inlet Valve Z Velocity at z=-5mm, x=0mm
-0.020
0.000
0.020
0.040
0.060
0.080
Position (m)
Jul 22, 1999 FLUENT 5.1 (3d, segregated, rngke)
Figure 4: Convergence – Z Component of Velocity at z = -5 mm •
Boundary Condition Uncertainty
As already noted, no experimental data on inlet turbulence intensity was available. To examine the effect of this uncertainty two cases with different values of turbulence intensity were compared. The comparison has been carried out by using the code’s default value of 10% in the first case and 1% in the second case. The velocity plots in Figures 6 and 7 show the results. For this application the assumption about inlet turbulence intensity is seen to have negligible effect on the predicted velocity field down stream of the narrow valve opening. However, other transported quantities are affected. This can be seen in figure 8 where predicted turbulent viscosity is compared for the two cases.
© ERCOFTAC 2000
page 62 of 94
not converged fully converged measurements
0.00 -0.20 -0.40 -0.60
Z Velocity (m/s)
-0.80 -1.00 -1.20 -1.40
Y Z
-1.60 -0.025
-0.020
X
Flow through an Engine Inlet Valve Z Velocity at z=+10mm, x=0mm
-0.015
-0.010
-0.005
0.000
0.005
0.010
0.015
0.020
Position (m)
Jul 22, 1999 FLUENT 5.1 (3d, segregated, rngke)
Figure 5: Convergence - Z Component of Velocity at z = +10mm
9.7.6. Conclusions This example has illustrated the effects of discretisation, level of convergence and boundary condition uncertainty on the quality of CFD prediction. It has been shown that the use of a second order scheme can yield more accurate results than those obtained using first order methods. However, in this particular application, the difference is very small and might not warrant the additional computational expense. The study of convergence error has clearly illustrated the need to ensure a properly converged solution and shows that unconverged solutions can be very misleading. The effect of boundary condition uncertainty has also been demonstrated. For this application it has been shown that uncertainty about inlet turbulence levels does not have any significant impact on the quality of velocity predictions within the system. However, other transported quantities have been shown to be more sensitive.
9.7.7. References Chen, A., Lee, K.C., Yianneskis, M., and Ganti, G. (1995), “Velocity Characteristics of Steady Flow Through a Straight Generic Inlet Port”, International Journal for Numerical Methods in Fluids, 21:571590.
© ERCOFTAC 2000
page 63 of 94
10% Intensity 1% Intensity measurements
1.50
1.00
0.50
Z Velocity (m/s)
0.00
-0.50
-1.00
Y Z
-1.50 -0.030 -0.020 -0.010 0.000 0.010 0.020 0.030 0.040 0.050 0.060 0.070 0.080
Position (m)
X
Flow though an Engine Inlet Valve Z Velocity at z=-5mm, x=0mm
Jul 22, 1999 FLUENT 5.1 (3d, segregated, rngke)
Figure 6: Boundary Conditions - Z Component of Velocity at z = -5mm 10% Intensity 1% Intensity measurements
0.00 -0.20 -0.40 -0.60
Z Velocity (m/s)
-0.80 -1.00 -1.20 -1.40
Y Z
-1.60 -0.025
X
Flow though an Engine Inlet Valve Z Velocity at z=+10, x=0
-0.020
-0.015
-0.010
-0.005
0.000
0.005
0.010
0.015
0.020
Position (m)
Jul 22, 1999 FLUENT 5.1 (3d, segregated, rngke)
Figure 7: Boundary Conditions - Z Component of Velocity at z = +10mm
© ERCOFTAC 2000
page 64 of 94
Intensity = 1% Intensity = 10% 7.00e-02
6.00e-02
5.00e-02
4.00e-02
Turbulent Viscosity [kg/m-s]
3.00e-02
2.00e-02
1.00e-02
0.00e+00 -0.03 -0.02 -0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
Y Position [m]
Turbulent Viscosity at Z = -5mm
Oct 06, 1999 FLUENT 5.2 (2d, segregated, lam)
Figure 8: Boundary Conditions – Turbulent Viscosity at z = -5mm
9.8. Test case G : Low speed centrifugal compressor (LSCC) Test case calculated with Code Fine/Turbo Version 3.1, supplied by NUMECA International.
9.8.1. Introduction The prediction of the flow and performance in turbomachinery blading is one of the major uses of CFD in industrial applications. The simulation given below demonstrates the use of CFD for a radial compressor calculation, and highlights the important effects of mesh size, turbulence modelling, and the geometric modelling of the clearance gap above the impeller vanes on both the flow and the performance predictions. The test case is the NASA Low Speed Centrifugal Compressor (LSCC), which is an experimental facility designed to duplicate the essential flow physics of high speed subsonic centrifugal compressor flow fields in a large size, low-speed machine. A complete description of the facility can be found in Hathaway et al. [1993].
9.8.2. Geometry and boundary conditions The compressor has an open back-swept impeller, followed by a vaneless diffuser. The impeller has 20 full blades with a backsweep of 55 deg. The inlet diameter is 0.870m and the inlet blade height is 0.218 m. The exit diameter is 1.524 m and exit blade height is 0.141m. The tip clearance between the blade tip and the shroud is 2.54mm, and is constant from the impeller inlet to exit. The design tip speed is 153 m/s, with a rotation speed of 1862 rpm and a design mass flow of 30 kg/s. Some aspects of the geometry can be seen in figure 1. The computational domain for these simulations has three blocks. The first block covers the flow passage, extending from 40% meridional shroud length upstream of the impeller to 15% impeller tip radius downstream of the impeller in the radial direction. The second block represents the region behind the blunt trailing edge of the blades. The third block occupies the space in the clearance gap above the tip of the blades. For convenience, a blunt blade tip is meshed, even though the real tested blade tips are rounded as reported by Chriss et al. [1996]. At the inlet plane, measured profiles of absolute flow angles and total pressure, uniform total temperature, turbulent kinetic energy k and dissipation ε are imposed. The inlet values of k and ε are estimated in the following way,
© ERCOFTAC 2000
page 65 of 94
k inlet = 15 . (TuV z ) 2 , and 2 2 / µ t = C µ k inlet / 10µ l ε inlet = C µ k inlet
where Cµ=0.09. The ratio of the turbulent viscosity to the laminar viscosity is set to 10 whereas the turbulent intensity Tu is set to 2%, these values being standard for turbulent internal flow computations. As in the experiments, the hub wall of the impeller is moving with the rotor blade, while the diffuser hub is stationary. The whole shroud is stationary and all solid wall boundaries are set adiabatic. Periodic conditions are enforced along the boundaries upstream and downstream of the impeller as well as in the gap.
9.8.3. Grid Three different H-grids, obtained with the auto-grid generation software IGG/AutoGrid of NUMECA, are used. The first mesh, Mesh-1 consists of 49x61x129, 13x61x33 and 13x13x65 points, in the tangential, radial and streamwise directions, for blocks 1, 2 and 3 respectively, resulting in 422,735 points in total. Hence there are 13 lines over the gap height and 13 lines across the blade profile. The second mesh, named Mesh-2 hereafter, has 61 points, instead of 49 in Mesh-1, over a pitch in Block 1 and 73 points, instead of 61 in Mesh-1, from hub to shroud in Blocks 1 and 2. The total grid number of Mesh-2 is 616,739. Particular care has been taken in order to ensure sufficient resolution in the end walls and blade boundary layers by controlling the position of the mesh points close to the solid surfaces. The grid location near the blade surfaces for Mesh-1 and for Mesh-2 produces a y+ value for the first cell center less than 1. A third mesh named Mesh-3, coarser than the two first counts 31x37x129, 13x37x33 and 13x17x65 points with in total 178201 mesh points. The y+ value at the first inner cell is about y+ =2.
9.8.4. Feature of the simulation The time-dependent Reynolds averaged Navier-Stokes equations, are solved with the finite volume flow solver EURANUS of the FINE/TURBO environment, with either the algebraic turbulence model of Baldwin-Lomax or two-equation k-ε models for the turbulence closure. The flow solver is based on a structured multiblock, multigrid approach, including non-matching block boundaries and incorporates various numerical schemes, based on either a central or an upwind space discretization. For the present calculations, the Baldwin-Lomax algebraic turbulence model and the linear k-ε model of Yang and Shih [1993], modified by Khodak and Hirsch [1996], are selected. In Yang and Shih 's model, the wall damping function used in the eddy viscosity is based on Ry = (k1/2y/ν). The calculations were performed with a second-order centered scheme, with second and fourth order artificial dissipation terms and a W-cycle multigrid technique. The numerical procedure applied a four-stage Runge-Kutta scheme, coupled to a local time stepping and implicit residual smoothing for convergence acceleration. The Reynolds number is based on the tip radius (1.524 m) and tip speed (153 m/s) and is equal to Re=1.610+7 whereas the Prandtl number is set to Pr=0.7. Thermodynamic properties are defined by Cp=1004.5 J/kg/oK and γ=1.4. For all computations, the convergence criteria is based on a mass flow error less than 0.1% and 4 to 5 order of magnitude reduction in the RMS residual.
9.8.5. Results Table 1 lists the cases of CFD results presented in this report. There are three different meshes, Mesh-1, Mesh-2 and Mesh-3, two different tip gap sizes, Design Gap and 50% Design Gap, and two turbulence models, Baldwin-Lomax (B-L) and linear k-ε models. Design flow rate Design Gap 50% Design Gap Off-design flow rate Design Gap
Mesh-1 B-L *
k-ε * * *
Mesh-2 k-ε *
Mesh-3 k-ε *
*
Table 1: Test cases definition
© ERCOFTAC 2000
page 66 of 94
At design flow rate of 30 kg/s, there are four calculated operating points for the design gap case, one computed with Mesh-2, two points computed with Mesh-1 and one point computed with Mesh-3, while for the 50% gap case, there is only one point computed with Mesh-1 and the Baldwin-Lomax turbulence model. At off-design flow rate, there are two calculated operating points for the design gap case with the k-ε turbulence model respectively on Mesh-1 and Mesh-2. On the coarser mesh Mesh-3, a single k-ε computation is carried on at design flow rate and tip gap. In many radial compressor simulations there is some considerable uncertainty in the running tip clearance of the blades, due to the effects of centrifugal and bending stresses in the blades, temperature distortions of the casing and manufacturing tolerances. The sensitivity to these effects have been examined by running a single simulation with a smaller tip gap. For the design gap case, the gap height is constant from the leading edge to the trailing edge while for the 50% design gap case the gap height is linearly reduced from the design height at the leading edge to half of it at the trailing edge. Figure 2 shows the effect of grid density, geometry and turbulence models on overall absolute total pressure ratio. We notice that the CFD data obtained with the k-ε model on Mesh-2 is in excellent agreement with the experimental data. For the coarser mesh, Mesh-1, the k-ε solution underpredicts slightly the pressure ratio. Reducing the tip gap size by 50% at the trailing edge, the pressure ratio increases about 0.4%. The computation on the coarse mesh Mesh-3 gives the lowest pressure ratio prediction with an underestimation of about 1%. The effect of turbulence model is shown at the design flow rate of 30 kg/s where the Baldwin-Lomax model gives a lower pressure ratio than the k-ε model on Mesh-1 and Mesh-2, putting the calculated operating point further away from the experiments. Figure 3 presents grid density (Mesh-1 and Mesh-2) and geometry influence for the k-ε model at the design flow rate. A comparison of reduced static pressure measured at six sections from 5% to 98% span is shown. As seen, in all the computations, with different meshes and different gaps, apparent difference in the predicted pressure at the spanwise sections can only be found downstream of the mid-chord. The best agreement with the experimental data is obtained on the fine mesh, Mesh-2. The 98% spanwise distribution clearly shows an effect of gap size. Figure 4 shows the effect of grid density, and turbulence model on the throughflow velocity at design conditions. A comparisons of computed throughflow isolines and the experimental results at the experimental station 165 are presented. This station is located at 6% of meridional shroud length upstream of the impeller exit. The throughflow velocity is normalized with the exit tip speed. We observe on Mesh-1 that the k-ε model reproduces better the isolines trend and in particular the location of the maximum throughflow velocity near the suction side. As expected, the Mesh-2 computation gives better results than those obtained on the coarser mesh, Mesh-1, showing the necessity to have sufficient grid density in the spanwise and circumferential directions. This grid density effect is also seen from the computation on Mesh-3 where the maximum velocity location is shifted towards the hub, showing a flow pattern similar to the prediction of the Baldwin-Lomax model on Mesh-1.
9.8.6. Conclusions The 3D viscous flow in the NASA LSCC impeller has been computed with different meshes, different tip gap sizes, and different turbulence models to demonstrate the sensitivity of the simulations to these effects. It is shown that compared to a low Reynolds k-ε turbulence model, the algebraic BaldwinLomax turbulence model underpredicts the global pressure ratio and is not able to predict the flow pattern with the same accuracy, showing the necessity to use precise turbulence models for quantitative simulations in turbomachines. It is also shown that grid density in the spanwise and circumferential directions play an important role in the accuracy of the global performance data such as the pressure ratio as well as on local flow features. The sensitivity to the uncertainties in the geometry is demonstrated by a 50% tip gap reduction where the gap height is linearly reduced from the design height at the leading edge to half of it at the trailing edge. The tip gap reduction slightly increases the global pressure ratio and affects the pressure distribution mainly on the pressure side near the shroud (in the region where a tip leakage vortex could be expected to lie).
© ERCOFTAC 2000
page 67 of 94
9.8.7. References Chriss, R. M., Hathaway, M. D., and Wood, J. R., (1996), "Experimental and Computational Results from the NASA Lewis Low-Speed Centrifugal Impeller at Design and Part-Flow Conditions", J. of Turbomachinery, Vol. 118, pp.55-65. Hathaway, M. D., Chriss R. M., Wood, J. R., and Strazisar, A. J., (1993), "Experimental and Computational Investigation of NASA Low-Speed Centrifugal Compressor Flow Field", J. of Turbomachinery, Vol. 115, pp.527. Khodak, A., and Hirsch, Ch., (1996), "Second Order Non-Linear k-e Models with Explicit Effect of Curvature and Rotation", Computational Fluid Dynamics'96, Proceeding of the Third ECCOMAS Computational Fluid Dynamics Conference, 690-696. Yang, Z., and Shih, T. H., (1993), "New Time Scale Based on k-e Model for Near-Wall Turbulence", AIAA J., Vol. 31, No.7, pp.1191-1198. S. Kang and Hirsch Ch., (1999) “Effect of flow rate on the development of three-dimensional flow in NASA LSCC impeller, based on numerical solutions, ISABE Paper No. 99-7225.
Figure 1: Geometry and grid for the LSCC impeller
© ERCOFTAC 2000
page 68 of 94
Figure 2: Effect of grid density, geometry and turbulence model on total pressure ratio
Figure 3: Effect of grid density and geometry on static pressure distribution at six sections from 5% to 98% span (circle: experiments; Full line: Mesh2, Design Gap; Dashed line: Mesh1, Design gap; Dashed dotted line: Mesh1, 50% gap)
© ERCOFTAC 2000
page 69 of 94
Baldwin-Lomax Mesh-1
K-E Mesh-2
K-E Mesh-1
K-E Mesh-3
Experiments Figure 4: Comparisons of computed throughflow isolines with the experimental results at the station located at 6% of meridional shroud length upstream of the impeller exit
© ERCOFTAC 2000
page 70 of 94
9.9. Test Case H : Turbulent flow in a model outlet plenum Test case calculated with CFX-TASCflow, Version 2.9, from AEA Technology Engineering Software.
9.9.1. Introduction The selected test case involves the modelling of two-dimensional, turbulent flow in a model nuclear reactor outlet plenum. Experimental results for this case were reported by Boyle and Golay (1983). It features the following phenomena which are challenging for turbulence models: • •
Curved shear layers Impingement
The assumption of two-dimensionality allows the concepts of iteration, truncation, and solution (discretization) error estimation, the assessment of boundary condition uncertainty and, ultimately, turbulence model performance to be demonstrated within reasonable computing time and memory requirements. The methods and algorithms used in this report can be directly transferred to three-dimensional applications.
9.9.2. Geometry and boundary conditions The model geometry is shown in Fig. 1. The model plenum is 228.6 mm high and 304.8 mm wide, with a chimney width of 57.2 mm. The depth of the model is 60.8 mm, yielding an aspect ratio based on width of only 0.2. The outlet duct has a width of 19.1 mm. Its length was set to four outlet duct widths to avoid recirculating flow at the outlet plane.
Figure 1: Test case geometry Because of the symmetry of the geometry, inlet conditions, and data, the flow in one half of the geometry was simulated. A symmetry boundary condition was used along the centre line. The steady-state inlet mass flow rate in the experiments of Boyle and Golay (1983) correspond to a Reynolds number of 70000, based on chimney width, H. The inlet boundary conditions are compiled in Table 1. At the outlet, the average static pressure, P, was set to zero. The no-slip condition in conjunction with the scalable wall function method developed by Grotjans and Menter (1998) was applied along all solid walls.
© ERCOFTAC 2000
page 71 of 94
Description x-Velocity y-Velocity z-Velocity Turbulence intensity Length scale
Symbol U V W Tu L
Value 0.0 0.823 0.0 0.0252 0.00295
Unit m/s m/s m/s m
Table 1: Inlet Boundary Conditions
9.9.3. Grid A sequence of hexahedral numerical grids was used for the calculations. The finer grids were derived from the coarse grids by successive division of the flux element sizes in the x - and y -direction. Figure 2 shows the basic grid structure using the coarsest grid. Details on the grid size and quality of all the grids are summarised in Table 2.
Flux elements Min. grid angle Max. grid angle Max. aspect ratio Max./min. volume ratio
Grid 0 25 x 35 77.8 102.0 12.6 142
Grid 1 50 x 70 77.8 102.0 13.6 184
Grid 2 100 x 140 77.8 102.0 6.6 216
Grid 3 150 x 210 77.8 102.0 5.4 228
Grid 4 200 x 280 77.8 102.0 5.4 234
Table 2: Grid quality parameters
Figure 2: Grid structure
9.9.4. Features of the simulation The code used allows the choice of various discretisation schemes for the convective terms in the model transport equations. In the present calculations the Modified Linear Profile Skew-Upwind method in combination with the Physical Advection Correction scheme was employed to discretise the momentum equations. This scheme has a spatial truncation error of second order.
© ERCOFTAC 2000
page 72 of 94
The experimental fluid was water. Constant properties are assumed in the calculations, with the density, ρ, set to 1000 kg/m3, and the molecular viscosity, µ, set to 6.725x10-4 Pa-s. Initial values used to start the coarse grid calculation were U = 0.0 m/s, V = 0.1 m/s, W = 0.0 m/s, P = 0.0 Pa, Tu = 0.025 , and L = 2.5 mm. Calculations on finer grids were started from interpolated coarse grid solutions. The standard k − ε turbulence model is used in this initial study. Future work will build on this first component to investigate the performance of different turbulence models. However, as turbulence model uncertainties stand at the end of the uncertainty chain, it is necessary to assess other sources of uncertainties before being able to make any scientific statement on the turbulence model performance. The first step in uncertainty or error quantification is to establish target quantities. The target quantities are problem and case dependent. They can be whole fields, norms, or integrated quantities like forces, or efficiencies. In this test case the target quantities are: • The efficiency, η . Efficiencies are often very sensitive to changes in computational parameters, making them good candidates for error and uncertainty analysis. • Profiles of the normal velocity, V , at a line corresponding to y = 113.4 mm. •
Profiles of the turbulent kinetic energy, k , at a line corresponding to y = 113.4 mm.
In this test case, attention is focussed on the target quantity efficiency, η . The efficiency is defined as the ratio of mass-averaged total pressure at the inlet to that at the outlet:
η=
Pt , 2 Pt ,1
.
9.9.5. Results Iteration Error Iteration errors arise from insufficient convergence criteria. The convergence criterion used is that the maximum normalized residual for the discrete momentum and mass conservation equations fall below a specified value. It is up to the code user to ensure that this convergence criterion is sufficient. A procedure for determining if a convergence criteria is sufficient is demonstrated below. This procedure is independent of the software used and also independent of the particular definition of the residuals used for measuring convergence. The procedure consists of monitoring target quantities as a function of the value of the given convergence criteria. For this case, the efficiency, η , as the most sensitive target quantity, is used to determine the necessary convergence level. When η ceases to change with the value of the convergence criterion, one can deem the solution to be converged. Conv. Criteria -2
1x10 5x10-3 1x10-3 5x10-4 1x10-4 5x10-5 1x10-5 5x10-6 1x10-6
η
∆η , %
0.682 0.636 0.615 0.610 0.513 0.514 0.515 0.515 0.515
32.4 23.5 19.4 18.4 0.0 0.0 0.0 0.0 0.0
Table 3: Iteration Error A detailed study has been done for Grid 1. Table 3 shows the efficiencies for a range of convergence levels. The data show that a convergence criterion corresponding to a maximum normalized residual of no more than 1x10-4 is necessary to reduce iteration errors to insignificant values. It is interesting
© ERCOFTAC 2000
page 73 of 94
that the results first level off at values of approximately 0.610, and then decrease by almost 20 % to the final value of 0.515. Table 4 shows the change in efficiencies for the other grids when going from a convergence criterion of 1x10-3 to a convergence criterion of 1x10-6. Conv. Criteria 1x10-3 5x10-4 1x10-4 5x10-5 1x10-5 5x10-6 1x10-6
Grid 0 0.711 0.701 0.575 0.574 0.574 0.575 0.575
Grid 1 0.615 0.610 0.513 0.514 0.515 0.515 0.515
Grid 2 0.459 0.460 0.489 0.493 0.498 0.498 0.499
Grid 3 0.443 0.475 0.492 0.494 0.496
Grid 4 0.399 0.435 0.474 0.487 0.493 0.495 0.495
Table 4: Iteration errors for various grids It must be emphasized, that different target quantities will generally have different convergence criteria for arriving at the same iteration error. Efficiencies, which are associated with viscous losses, usually require tougher convergence criteria than, for instance, pressure forces. Truncation Error Truncation, or round-off, errors should be checked to ensure that the accuracy of the chosen computer is sufficient for the problem. Round-off errors can be quantified by performing the same calculation with single and double precision executables and by comparing the results for the target quantities. Table 5 shows the evolution of efficiencies predicted for a single and double precision version of the code on Grid 3. Conv. Criteria 1x10-3 1x10-4 5x10-5 1x10-5 1x10-6
Single Prec. 0.443 0.475 0.492 0.494 0.494
Double Prec. 0.514 0.485 0.485 0.495 0.495
Table 5: Comparison of efficiency for single and double precision calculations Both versions converge to the same result for this case. It is, however, evident that conclusions on the influence of machine accuracy and round-off problems should only be made once the problem is sufficiently well converged. Solution Error Solution error is defined as the difference between a solution on a grid with finite grid width and the `exact' solution of the very same problem on a grid with infinitely small grid width. In most practical cases the exact solution is not known. Hence solution error can only be estimated, see Roache (1997). A method based on Richardson extrapolation is used in this report. The first target quantity to which the solution error estimation method is applied, is the efficiency, η . Table 6 shows the predicted efficiencies for the various grids as a function of the number of grid nodes in the y -direction. The estimated efficiency error, Eη , is defined by (Roache, 1997):
Eη = where
e r −1
(
p
)
e = (η fine − η coarse ) / η fine , r is the grid refinement ratio, and p=2 is the order of the discretiza-
tion.
© ERCOFTAC 2000
page 74 of 94
Grid 0 Grid 1 Grid 2 Grid 3 Grid 4
Ny
η
35 70 140 210 280
0.575 0.515 0.493 0.494 0.495
Eη , % 3.88 1.49 -0.16 -0.26
Table 6: Efficiency solution errors This solution error analysis indicates that a grid number of more than 140 is required for an assessment of efficiency which is independent of the computational grid. Figure 3 shows the velocity profiles at y = 113.4 mm. The calculated velocity profiles converge to a single `grid-independent’ solution. The results obtained on Grids 3 and 4 are almost indistinguishable. The agreement with data is satisfactory. Large discrepancies appear for the turbulent kinetic energy, k, shown in Fig. 4. The standard k-ε turbulence model over predicts the level of turbulent kinetic energy because it is insensitive to streamline curvature. This problem will be addressed by application of a second-moment closure turbulence model.
Figure 3: Vertical velocity distribution
Figure 4: Turbulent kinetic energy profile
9.9.6. Conclusions and outlook The results and analysis presented for this test case demonstrate a means for determining required levels of convergence from a suitable convergence criteria. The application of an approach for estimating solution error has been shown as a means for assessing the necessary grid resolution for grid independence of a solution. A basis for subsequent analysis of the performance of alternative turbulence models in predicting the flow phenomena particular to this case has been established. Application of a second-moment closure turbulence model will be reported in Scheuerer and Stokes (1999).
© ERCOFTAC 2000
page 75 of 94
9.9.7. References Boyle, D.R., and Golay, M.W. (1983), ‘Measurement of a Recirculating, Two-dimensional, Turbulent Flow and Comparison to Turbulence Model Predictions. I: Steady-State Case’, J. Fluids Eng., Vol. 105, pp. 439-445. Grotjans, H., Menter F. R. (1998), `Wall Functions for General Application CFD Codes’, Proc., ECCOMAS ’98, 4th European Computational Fluid Dynamics Conference, Athens, eds. Papailiou, K. D. et al., J. Wiley & Sons, pp. 1112 – 1117. Scheuerer, G., Stokes, J. A. (1999), `Application of the ERCOFTAC Best Practice Guidelines to Turbulence Model Performance Assessment’, AEA Technology GmbH, Internal Report AEAT/TR-9926.
© ERCOFTAC 2000
page 76 of 94
10. Suggestions and need for further work 10.1. Revision of best practice guidelines The guidelines have been written and edited making use of extensive consultation with CFD code vendors, code developers, academic experts and code users and this gives the guidelines wide support. The editors have tried to ensure that the guidelines include those 20% of the most important general rules of advice that cover roughly 80% of the problems likely to be encountered. In the nature of things, certain useful guidelines will have been forgotten and others will have been included which are superfluous. In addition CFD is a rapidly developing subject and technical developments will lead to inevitable changes in the most useful guidelines. Because of this, it is intended that the guidelines will be revised after they have been in use for a certain period of time. Readers with suggestions for corrections and improvements should forward them directly to the ERCOFTAC co-ordination centre at [email protected] from where they will be passed on the ERCOFTAC SIG on Quality and Trust in Industrial CFD, or should directly contact the editors by email ([email protected] and [email protected]) who are interested in any feedback.
10.2. Extension of the best practice guidelines This first edition of the best practice guidelines is limited to single-phase, compressible and incompressible, steady and unsteady, turbulent and laminar flow with and without heat transfer. The guidelines need extension to cover specific types of flow not dealt with in the current edition. In particular, separate guidelines on the computation of the following aspects would be of interest: •
Combustion
•
Multi-phase flows (particles, drops, sprays, films, free surfaces, bubbles, cavitation)
•
Flows with radiation
•
Non-Newtonian flows
•
Hypersonic and supersonic flows with strong shocks and compressibility
•
Fluid-structure interaction
•
Multi-physics problems
10.3. Application procedures In addition to extending the guidelines to cover different types of flow, the guidelines could be extended to provide application procedures with detail information on how to achieve the best engineering prediction in specific application areas. The objective here would be to examine the performance of various physical models on specific test cases, to offer advice on how to achieve the best results in specific cases and to suggest what level of accuracy is to be expected in key predicted parameters. The list of possibilities for this work is very large, but specific examples with major industrial impact might be taken from the following areas: •
External flows on airframes and wings.
•
Turbomachinery component performance prediction.
•
Heat transfer and cooling of gas turbine blades.
•
Internal and external building flows.
•
Pollutant dispersion problems.
•
Automotive engine flows.
•
Internal automotive cooling flows.
•
External automotive flows.
•
Internal flows for the power industry.
© ERCOFTAC 2000
page 77 of 94
11. Checklist of best practice advice for industrial CFD This section contains a checklist of the best practice guidelines outlined in the chapters given above. The list is organised differently to the body of the report, in that the structure roughly follows the important chronological activities normally required to complete a CFD-project (as given in section 2.1). This structure should help the reader to rapidly check the quality of a CFD simulation.
11.1. Guidelines on the training of CFD users A CFD user for non-routine applications should have good training and knowledge in classical fluid mechanics (including experimental, analytical and computational aspects), a broad understanding of numerical methods as applied to fluid flows, and detailed knowledge of the application being examined. This means that the user should be able to understand the limitations of the models used (e.g. turbulence, radiation, buoyancy driven flows, numerical models, etc.) and be aware of their relevance when being applied to a particular problem. ¾ The training and education requirement for more routine applications can be less stringent, provided that clear guidelines or procedures have been established for the code being used. An example of a routine application would be the simulation of a standard component in a design environment where many previous designs have been calculated and only relatively small changes in geometries and boundaries conditions occur. ¾ In both routine and non-routine applications, training on the use of the specific CFD code with the solution of realistic exercises is needed. ¾ Training should be repeated regularly, and the user should frequently read the release notes and manuals issued by the code supplier, as these are a good source of educational information.
¾
11.2. Guidelines on problem definition ¾
Give careful thought at the outset of a CFD project to the requirements and objectives of the simulation. The following points should be considered: • •
• • •
• • • • • • •
Are the objectives of the simulation clearly defined? Is a CFD simulation really appropriate or is a simple analysis or an experiment needed? Circumstances where CFD is generally considered appropriate are when the flow physics or geometry are sufficiently complex that simpler modelling approaches are invalid, or where experimental data is lacking and cannot be obtained because it is too costly, dangerous or impractical. What local/global quantities are needed from the simulation? What are the requirements on accuracy? What level of validation is necessary? Is this a routine application, where validation and calibration has already been carried out on similar flow fields, and where only relatively small changes can be expected from earlier similar simulations? Or is it a non-routine application, where little earlier validation work has been done? What are the important flow physics involved (heat transfer, radiation, steady, unsteady, compressible, single phase, laminar, turbulent, transitional, internal, external, etc.)? What is the area of primary interest (domain) for the flow calculation? Is the geometry well defined, and can symmetry or periodic boundary conditions be applied? Are the inlet and outlet boundary conditions clearly specified? What sensitivity tests should be carried out to assess the accuracy of the simulation? What level of computational resources is needed for the simulation (memory, disk space, cpu time, serial or parallel run because of time requirement) and are these available? What are the documentation/reporting requirements?
11.3. Guidelines on solution strategy ¾
Translate the problem definition into a solution strategy involving issues and questions that have been addressed in this document, such as: • • •
Mathematical and physical models, numerical schemes Solution algorithm, pressure or density based solution method. Turbulence and other physical models.
© ERCOFTAC 2000
page 78 of 94
• • • • •
Available code/solver. Grid design and computational mesh. Boundary conditions. Convergence criterion. Validation and sensitivity test requirements
11.4. Guidelines on global solution algorithm Check the adequacy of the solution procedure with respect to the physical properties of the flow by reference to the code user manual, for instance compressible or incompressible flow. ¾ As a first step in the simulation process, use the parameters that control convergence of the solution algorithm (e.g. relaxation parameters, damping factors or time steps) as suggested by the CFD-code vendor or developer. ¾ If it is necessary to change these parameters to aid convergence, do not change too many parameters in one step, as it then becomes difficult to analyse which of the changes have influenced the convergence. In case of persistent divergence see sections on boundary conditions grid discretisation and convergence errors. ¾ If a steady solution has been computed and there is reason to be unsure that the flow really is steady, then carry out an unsteady simulation with the available steady flow field as the initial condition, preferably with a small perturbation imposed. Examination of the time-development of the physical quantities in the locations of interest will identify whether the flow is steady or not.
¾
11.5. Guidelines on validation of models ¾
¾
¾
¾ ¾
¾
¾
Before using a CFD code incorporating a specific model for a particular application, validate it against test data for a similar application with similar flow structures and flow physics. A design application, where accurate performance and flow-field data is needed, will require more detailed validation than an application involving a preliminary investigation of a flow. Validation of engineering cases should be aimed to check the most prominent and relevant flow features, rather than to capture every detail in the calculation. For the improved validation of complex flows, consider the taxonomy of the flow and validate the models being used on relevant simpler flows. For example, the validation of a code and turbulence model for a radial pump application should include not only a similar radial pump test case, but also test cases demonstrating the typical flow structures found in the pump, such as decelerating flow in a cascade of blades, flow in a rotating channel, secondary flow in bends, and leading and trailing edge flows. Simple test cases used in validation will allow finer grids than those used in real applications so that results of calculations at the typical grid resolution expected in the application should also be investigated. Carefully check the quality and accuracy of any experimental data used as a validation test case, to provide confidence in the experiment and the quality experimental data. The experiments should attempt to capture the relevant boundary conditions and initial conditions where appropriate. If the experimental results are found to be sensitive to the boundary conditions and experimental set-up then this should be reported. The experimental data should include measurements of the quantity which is of interest and needed for CFD validation. This is best done by close collaboration between experimentalists and CFD analysts. If collaboration between experimentalists and CFD analysis is not possible care should be taken that the physically relevant quantities have been measured and that the quality and accuracy of the data is known. The validation data set should be independent of any data set used to construct the CFD model. Carry out the validation simulations in accordance with the guidelines outlined in this document. In particular, one of the validation calculations should use about the same mesh size as will be used in the real application. A simple test of validity should be carried out by examining particular flow features. The global mass balance should be checked, the total temperature should not vary if no heat or work is being added to the flow, and in areas away from solid boundaries the total pressure may often be assumed to be constant. Predictions of the entropy in a flow field can also be a helpful validation check. When using CFD for a design application, examine the flow of a well understood and well tested baseline design comprising the same important flow features as a means of calibrating the ex-
© ERCOFTAC 2000
page 79 of 94
pected accuracy of the code to this application. The reliability of predictions of a new design can then be assessed by taking into account the calibration of the code with the baseline design. Such a calibration allows a more accurate estimate of the predicted performance of the new design to be made by consideration only of the differences and trends between the new design and the baseline design. ¾ Document the validation efforts. In any CFD calculations which are carried out without undergoing a thorough validation process (which is sometimes necessary as there may be no suitable measurements, or the time available allows no alternative), then attention should be drawn in the documentation of the results to the speculative nature of the simulation, with comments on why the user believes that the results are correct.
11.6. Guidelines on turbulence modelling 11.6.1. General guidelines on turbulence modelling ¾
¾
¾ ¾
¾
Ensure that low numerical and convergence errors have been achieved in turbulent flow simulations. The relevance of turbulence modelling only becomes significant in CFD simulations when other sources of error, in particular the numerical and convergence errors, have been removed or properly controlled. Clearly no proper evaluation of the merits of different turbulence models can be made unless the discretisation error of the numerical algorithm is known, and grid sensitivity studies become crucial for all turbulent flow computations. Be aware that there is no universally valid general model of turbulence that is accurate for all classes of flows. Validation and calibration of the turbulence model with test data are necessary for all applications5. If possible, examine the effect and sensitivity of results to the turbulence model by changing the turbulence model being used. When using a particular turbulence model, check the published literature with regard to the known weaknesses of the model. The weaknesses of the standard k-ε model (Launder and Spalding [1974]) ,which is the most commonly used model in industrial applications, are listed below together with some indications of possible palliative actions which might be fruitfully considered. Decide whether a wall function method, in which the near-wall region is bridged with wall functions, or a low Reynolds number model, in which the flow structure in the viscous sub-layer is resolved, is to be used. This decision will be based on the available resources and the requirements for resolution of the boundary layer. Wall function methods are not valid in the presence of separated regions and/or strong three dimensional flows. The validity of the wall function approach or the use of a low Reynolds number model should be examined for the flow configuration under study.
11.6.2. Guidelines on wall functions Check the lower limit of y+. In the commonly used applications of wall functions, the meshing should be arranged so that the values of y+ at all the wall-adjacent integration points is only slightly above the recommended lower limit given by the code developers, typically between 20 and 30 (the form usually assumed for the wall functions is not valid much below these values). This procedure offers the best chances to resolve the turbulent portion of the boundary layer. It should be noted that this criterion is impossible to satisfy close to separation or reattachment zones because y+ will approach zero there, independent of the distance to the wall. ¾ Check the upper limit on y+. In the case of moderate Reynolds number, where the boundary layer only extends to y+ of 300 to 500, there is no chance of accurately resolving the boundary layer if the first integration point is placed at a location with the value of y+ of 100. ¾ Check the resolution of the boundary layer. If boundary layer effects are important, it is recommended that the resolution of the boundary layer is checked after the computation. This can be achieved by a plot of the ratio between the turbulent to the molecular viscosity, which is high inside the boundary layer. Adequate boundary layer resolution requires at least 8-10 points in the layer. ¾
5
Calibration in this sentence refers to testing the ability of a CFD code to predict global quantities of interest for specific geometries of engineering design interest (see definition in chapter 2) and not to tuning the coefficients of a turbulence model, which cannot be recommended for a novice user. © ERCOFTAC 2000
page 80 of 94
Exercise care when calculating the flow using different schemes or different codes with wall functions on the same mesh. Cell centred schemes have their integration points at different locations in a mesh cell than cell vertex schemes. Thus the y+ value associated with a wall-adjacent cell differs according to which scheme is being used on the mesh. ¾ Check that the correct form of the wall function is being used to take into account the wall roughness. An equivalent roughness height and a modified multiplier in the law of the wall must be used.
¾
11.6.3. Guidelines on low Reynolds number models The value of y+ at the first node adjacent to the wall should be less than 4, and preferably close to unity. ¾ Depending on the Reynolds number, there should be between five and ten mesh points between the wall and y+ equal to 20, which likely results in thirty to sixty points inside the boundary layer for adequate boundary layer resolution.
¾
11.6.4. Guidelines on weaknesses of the standard k-ε model ¾
¾
¾ ¾
¾
The turbulent kinetic energy is over-predicted in regions of flow impingement and re-attachment leading to poor prediction of heat transfer and the development of boundary layer flow around leading edges and bluff bodies. The high turbulence levels predicted upstream of a stagnation point are transported around the body and the real boundary layer development is swamped by this effect. Kato and Launder [1993] have proposed an ad-hoc modification (called the Yap correction) to the equation for the production term of k appearing in the k and ε equation, which is designed to tackle this problem. The RNG k-ε model (Yakhot et al. [1992]) includes a modification to the transport equation for ε and may also improve predictions in this area. The problems depend on the free-stream values of k and ε and may not always occur. Highly swirling flows are often poorly predicted due to the complex strain fields and regions of recirculation in a swirling flow are often under-estimated. These deficiencies require turbulence models with improved physics, such as non-linear k-ε models, algebraic Reynolds stress models, or full Reynolds stress models with low Reynolds number capability limitation (for example, the reviews in Patel et al. [1985] and Wilcox [1998]). Mixing is often poorly predicted in flows with strong buoyancy effects or high streamline curvature. Again this deficiency requires improved physics, see comment above. Flow separation from surfaces under the action of adverse pressure gradients is often poorly predicted. The real flow is likely to be much closer to separation (or more separated) than the calculations suggest. The SST version of Menter's k-ω based, near wall resolved model (Menter [1993, 1996]) can offer a considerable improvement. Flow recovery following re-attachment is often poorly predicted. If possible, avoid the use of wall functions in these regions, but low Reynolds number models may also have some problems. A promising possibility is the use of a length-scale limiting device (Ince and Launder [1995]).
The far-field spreading rates of round jets are predicted incorrectly. The use of non-linear k-ε models should be investigated for these problems (Apsley et al. [1995]). ¾ Turbulence driven secondary flows in straight ducts of non-circular cross section are not predicted at all. Linear eddy viscosity models cannot capture this feature. Use RSTM or non-linear eddy viscosity modelling. ¾ Dilation effects on turbulence production, which may become important for compressible flows, are not accounted for. These may be important in mixing layers at moderate Mach numbers. Sarkar et al. [1989] and Zeman [1990] have proposed modifications to the k-ε model that correct this deficiency but the database is small. ¾
Laminar and transitional regions of flow cannot be calculated with the standard k-ε model but nonlinear modelling may help in this respect. Transition is an active area of research in turbulence modelling. Some simple practical advice is given in section 4.7 below. ¾ The vertical heat flux in a buoyancy affected horizontal boundary layer, or in a thermal plume, is poorly represented by an eddy viscosity model which assumes a turbulent Prandtl number to calculate the turbulent transport of heat. RSTM models including the solution of equations for the three heat fluxes is currently the only option if this is important. ¾
© ERCOFTAC 2000
page 81 of 94
11.6.5. Guidelines on transitional flows Use experimental data to check whether the flow contains extensive regions of laminar or transition flow which would be incorrectly estimated by the k-ε model with wall functions, or poorly predicted by low Reynolds number k-ε models. ¾ Model transiton by intervention in the code to switch the turbulence model on or off at predetermined locations determined by experiment ¾
11.7. Guidelines on definition of geometry Check and document that the geometry of the object being calculated is the geometry as intended taking into account items such as design modifications, wear, tolerances, deflection due to loading, fouling, etc. If the geometry is substantially altered or deformed by the aerodynamic, mechanical or thermal loading, then use a some structural or mechanical calculation to determine the exact geometry. ¾ Check that the geometry is defined in the correct co-ordinate system and with the correct units requested by the CFD-code. ¾ Check that the transfer of geometrical data from a CAD system to a CFD-system, or grid generation software, does not involve loss of surface representation accuracy and that the final surface grid conforms to the intended geometry. CAD-systems often define the geometry in millimetres and this must be converted to SI-units if the code assumes that the geometry information is in these units (which is commonly done by most codes). It may be necessary to clean up the CAD geometry so that surfaces and curves meet at the intended end point locations. Visual display of the geometry helps here. ¾ Consider whether local details of the geometry can be omitted or require grid refinement. In general, it is not necessary to explicitly include geometrical features that have dimensions below that of the local grid size provided that they are taken into account in the modelling (e.g. roughness in wall layer, porous elements or detail of pipe-work in process plant explosion calculations). In areas where local detail is needed, such as in the neighbourhood of fine edges or small clearance gaps, then grid refinement in the areas with fine details should be used,. ¾
11.8. Guidelines on grids and grid design ¾
¾
¾
¾
¾
¾
¾
Choose a suitable global topology of the mesh (for example, a scalable grid topology where the grid can be refined by parametric changes without degrading the grid quality) to help satisfy the specific code's requirements with regard to skewness, aspect ratio and expansion ratios, as outlined below. Ensure that the extent of the computational domain has been chosen to capture relevant flow and geometrical features, and if necessary examine the sensitivity of the calculation to the choice of computational domain. Assess which geometrical features can be omitted, for example, those that have dimensions below that of the local grid size. In areas where local detail is needed, then consider the use of local grid refinement to capture fine geometrical details. If grid refinement is used the additional grid points should lie on the original geometry and not simply be a linear interpolation of more grid points on the original coarse grid. The use of distributed losses or porosity might also be considered to take into account highly obstructed zones or fine detail of some obstructions (for example, porous filters and packings), whereby the code manual should be consulted for the interpretation of the turbulent fields obtained. Avoid highly skewed cells. For hexahedral cells or prisms the included angles between the grid lines should be optimised in such a way that the angles are approximately 90 degrees. Angles with less than 40 or more than 140 degrees often show a deterioration in the results or lead to numerical instabilities. Tetrahedra should tend to have their four angles equal. Avoid highly warped cells, that is cells with large deviations from co-planar faces. Warp angles (measured between the surface normals of triangular parts of the faces) greater than 75 degrees can lead to serious convergence problems and deterioration in the results. Avoid non-orthogonal cells near boundaries. The angle between the grid lines and the boundary of the computational domain (the wall or the inlet- and outlet-boundaries) should be close to 90 de-
© ERCOFTAC 2000
page 82 of 94
¾ ¾
¾
¾ ¾
¾ ¾ ¾
¾
grees. This requirement is stronger than the requirement given for the grid angles in the flow field far away from the domain boundaries. Avoid the use of tetrahedral elements in boundary layers. Avoid aspect ratios that are too high. Away from boundaries, ensure that the aspect ratio (the ratio of the sides of the elements) is not too large. This aspect ratio should be typically not larger than 20 to 100 (depending on the flow solver) in important regions of the flow domain but may be larger in non-critical regions, whereby precise values should be supplied in the code manual. The accuracy and possible convergence difficulties in such cases depends greatly on the flow direction. Near walls this restriction may be relaxed and indeed it can be beneficial to have high aspect ratio in the boundary layers. Observe any specific requirements on mesh expansion ratios. The specific code's requirements for cell mesh stretching or expansion ratios (rates of change of cell size for adjacent cells) should be observed. The change in mesh spacing should be continuous and mesh size discontinuities be avoided, particularly in regions of large changes. Note that some codes offer automatic grid adaptation techniques, but that depending on the software being used these might not improve the grid quality (skewness, aspect ratio). Use a finer and more regular mesh in critical regions with high flow gradients or with large changes, such as shocks, regions with high shear, regions with significant changes in geometry or where suggested by error estimators. Make use of local refinement of the mesh in these regions, in accordance with the selected turbulence wall modelling. Ensure high geometric precision of the periodic grid interface when using periodic boundary conditions. Avoid the presence of arbitrary mesh coupling, non-matching cell faces, grid refinement interfaces or extended changes of element types in the critical regions of high flow gradients. Check that the assumptions made with regard to critical regions of high flow gradients and large changes when setting up the grid agree with the result of the computation, and rearrange grid points if found to be necessary. For each class of problem, make use of a grid dependency study to analyse the suitability of the mesh and to give an estimate of the numerical error in the simulation. Ideally, at least three significantly different grid resolutions should be used. Strictly, one should double the grid twice in each direction and then apply Richardson extrapolation (Roache [1997]) to be really sure. If this is not feasible, apply selective local refinement of the grid in critical flow regions of the domain to allow greater factors, or try to compare different order of spatial discretisations on the same mesh.
11.9. Guidelines on boundary conditions 11.9.1. General guidelines on boundary conditions ¾
¾
¾ ¾
¾
Ensure that appropriate boundary conditions are available for the case being considered. For swirling flows at an outlet consult the manual to check that appropriate boundary conditions have been implemented (for example, radial equilibrium of pressure field instead of constant static pressure). Special non-reflecting boundary conditions are sometimes required for outflow and inflow boundaries where there are strong pressure gradients (Giles [1990]). Check whether the CFD code allows inflow at open boundary conditions. If inflow cannot be avoided at an open boundary then ensure that the transported properties of the incoming fluid including turbulence properties are properly modelled. Examine the possibilities of moving the domain boundaries to a position where the boundary conditions are more readily identified, are well-posed and can be precisely specified. Check whether upstream or downstream obstacles (such as bends, contractions, diffusers, blade rows etc.) outside of the flow domain may be present which might have significant effects on the flow distribution. Often information on the components upstream of the inlet boundary or downstream of the outlet boundary of the domain are lacking or not available at the beginning of a project. For each class of problem that is of interest, carry out a sensitivity analysis in which the boundary conditions are systematically changed within certain limits to see the variation in results. Should any of these variations prove to have a sensitive effect on the simulated results and lead to large
© ERCOFTAC 2000
page 83 of 94
changes in the simulation, then it is clearly necessary to obtain more accurate data on the boundary conditions that are specified. ¾ Be aware that the code may have default boundary conditions for regions of the domain boundaries where the user has specified nothing.
11.9.2. Guidelines on inlet conditions ¾
Carry out a sensitivity analysis in which the key inlet boundary conditions are systematically changed within certain limits. Depending on the problem, the key parameters that might be examined are: • • • •
¾
Inlet flow direction and magnitude. Uniform distribution of a parameter or a profile specification, for example a uniform inlet velocity (slug flow) or an inlet velocity profile. Physical parameters. Turbulence properties at inlet.
Particular care is needed with the specification of inlet mass flow as a boundary condition for compressible flow-fields with supersonic regions, as the flow may be choked and unable to pass the mass flow specified.
11.9.3. Guidelines on specification of turbulence quantities at an inlet Use verified quantities from experiments as inlet boundary conditions for turbulent kinetic energy k and dissipation ε, if these are available, as the magnitude can significantly influence the results. If there are no experimental data available, then the values need to be specified using sensible engineering assumptions (see below), and the influence of the choice should be examined by sensitivity tests with different simulations. ¾ Specify values of the turbulent kinetic energy k which are appropriate to the application. These values are often specified through a turbulence intensity level Tu, which is defined by the ratio of the fluctuating component of the velocity to the mean velocity. In external aerodynamic flows over airfoils the turbulence level is typically Tu = 0.003 (0.3%). In atmospheric boundary layer flows the level can be two orders of magnitude higher (Tu = 0.30 (30%)) and details of the actual boundary layer profiles are needed. In internal flows the turbulence level of Tu = 0.05 to 0.10 (5 to 10%) is usually appropriate. ¾
¾
Specify values of the turbulent length scale, as an equivalent parameter for the dissipation ε, that are appropriate to the application. For external flows remote from boundary layers a value determined from the assumption that the ratio of turbulent and molecular viscosity µT/µ is between 1 and 10 is a reasonable guess. For internal flows a constant value of length scale derived from a characteristic geometrical feature can be used (e.g. 1 to 10% of the hydraulic diameter for internal flows). For simulations in which the near-wall region is modelled with a low Reynolds number model the length scale should be based on the distance to the wall and be consistent with the internal modelling in the code.
If more sophisticated distributions of k and ε are used, check that these are consistent with the velocity profile. An inconsistent formulation may lead to an immediate unrealistic reduction of the turbulence quantities after the inlet. ¾ For RSTM models the stresses themselves need to be specified. If these are not available, as is often the case, make an assumption of isotropic flow conditions (that is with the normal stresses given by 2/3 k and with zero shear stresses). ¾ For inlets that are representative of fully developed pipe or channel flows, algebraic profiles might be used. ¾
Check the consistency of the definitions of k and ε by making a plot of the ratio of turbulent to molecular viscosity µT /µ. Note that the ratio of turbulent to laminar viscosity depends on Tu, L, k and ε as µT/µ ~ k2·ε or µT/µ ~ Tu·L. ¾ In cases where problems arise, move the inflow boundary sufficiently far from the region of interest so that a natural inlet boundary layer can develop.
¾
© ERCOFTAC 2000
page 84 of 94
11.9.4. Guidelines on outlet conditions ¾
¾ ¾
¾
¾
¾
Place open boundary conditions (outflow, pressure) as far away from the region of interest as possible and avoid open boundaries in regions of strong geometrical changes or in regions of recirculation. If the latter point cannot be achieved, consult the code manual to see what boundary condition is considered best. Take special care with the orientation of outlet planes with regard to the mean flow, especially when the boundary condition consists of a constant static pressure profile Select the boundary conditions imposed at the outlet to have only a weak influence on the upstream flow. Extreme care is needed when specifying flow velocities and directions on the outlet plane. The most suitable outflow conditions are weak formulations involving specification of static pressure at the outlet plane. Particular care should be taken in strongly swirling flows where the pressure distribution on the outlet boundary is strongly influenced by the swirl. It is therefore not acceptable to specify constant static pressure across the outlet. Be aware of the possibility of inlet flow inadvertently occurring at the outflow boundary, which may lead to difficulties in obtaining a stable solution or even to an incorrect solution. If it is not possible to avoid this by relocating the position of the outlet boundary in the domain, then one possibility to avoid this problem is to restrict the flow area at the outlet, provided that the outflow boundary is not near the region of interest. If the outflow boundary condition allows flow to re-enter the domain, appropriate Dirichlet conditions should be imposed for all transported variables. If there are multiple outlets, impose either pressure boundary conditions or mass flow specifications depending on the known quantities.
11.9.5. Guidelines on solid walls Care should be taken that the boundary conditions imposed on solid walls are consistent with both the physical and numerical models used (for example, adiabatic walls, local heat fluxes, free-slip boundaries, roughness, moving or rotating walls). In cases with moving or rotating walls the boundary conditions need to be specified consistent with the motion of the walls. ¾ Great care is needed if roughness on the wall is not negligible. If no detailed information is available, significant levels of uncertainty can arise through incorrect specification of equivalent roughness on solid walls. ¾
11.9.6. Guidelines on symmetry and periodicity planes ¾
If symmetry or periodicity planes cross the inlet or outlet boundaries then care should be taken to specify inlet or outlet variables that are consistent with these. Symmetry and periodicity planes assume that the gradients perpendicular to the plane are either zero (for symmetry) or determined from the flow field (periodicity).
11.9.7. Guidelines on initial guess or initial condition Be aware that in special cases the initial guess can influence the converged flow field e.g. buoyant flows, flows with bifurcations or flows with a hysteresis depending on the direction of the changes. ¾ Care should be taken to impose realistic conditions that consistent with the flow equations as the initial condition of an unsteady-state flow calculation (for example, use a steady-state flow simulation as the initial condition for a gas release simulation).
¾
¾
Special attention should be directed to the initial guess of k and ε, especially with low Reynolds number models.
11.9.8. Guidelines on uncertainties with steady flow, symmetry and periodicity Check carefully whether the geometry is symmetric or whether a geometrical distortion or disturbance in the inlet conditions is present which can trigger asymmetric solutions. ¾ Estimate the Reynolds-number of the inflow and check whether the flow could be asymmetric, turbulent and/or unsteady (e.g. by sources or literature). Coupling between competitive physical phenomena might also be considered as a potential source of unsteadiness. ¾ After obtaining a steady solution, switch to the transient mode and check whether the solution remains stable. ¾
© ERCOFTAC 2000
page 85 of 94
Switch to the transient mode if it is difficult to get a converged steady solution, especially if there is an oscillation of the residuals. ¾ In case of doubt, the simulation should be unsteady and without symmetry assumptions as boundary conditions. ¾
11.9.9. Guidelines on physical properties Ensure that the correct physical properties are specified (e.g. water, air, ideal gas). If constant property assumptions are being used make sure that they are valid. ¾ Ensure that the specified value is the correct one (e.g. dynamic and not kinematic viscosity). ¾ Ensure that the values are specified in the correct units. ¾ Check whether user coding is being used to specify non-constant properties via tabulated values, as this can be an additional source of uncertainty (e.g. interpolation, extrapolation, user coding errors, inconsistency with base code). ¾
11.10. Guidelines on the solution of discretised equations using a CFD code 11.10.1. Guidelines on convergence ¾
¾
¾ ¾
¾
¾ ¾
¾
Carefully define solution-sensitive target quantities for the integrated global parameters of interest and select an acceptable level of convergence based on the rate of change of these target parameters (such as heat load, mass flow, lift, drag, and moment forces on a body). Convergence of a simulation should not be assessed purely in terms of the achievement of a particular level of residual error. For each class of problem you are interested in, carry out a test of the effect of converging to different levels of residual on the integrated parameter of interest (this can be a single calculation that is stopped and restarted at different residual levels). This test demonstrates at what level of residual the parameter of interest can be considered to have converged and identifies the level of residual that should be aimed at in similar simulations of this class of problem. Be aware that different codes have different definitions of residuals. Always check the convergence on global balances (conservation of mass, momentum and energy) where possible, such as the mass flow balance at inlets and outlets and at intermediate planes within the flow domain. Check the maxima and minima of all variables which may be limited by the code (for example k and ε which should always be positive), and the number of cells where any limitation on values is active. Check not only the residual itself but also the rate of change of the residual with increasing iteration count Monitor the solution in at least one selected point in a sensitive area to see if the region has reached convergence, but be aware that this alone is insufficient as a criterion for convergence as the solution may be "stuck" far from convergence but appear converged if solution variables are unchanging. For calculations that are proving difficult to converge, the following advice may be helpful: • • • • •
Use more robust numerical schemes during the first (transient) period of convergence and switch to more accurate numerical schemes as the convergence improves. Reduce parameters controlling convergence, such as the under-relaxation parameters or the CFL6 number. If the solution is heavily under-relaxed, increase the relaxation factors at the end to see if the solution holds. Check whether the grid quality in areas with large residual has any effect on the convergence rate. Look at the residual distribution and associated flow field for possible hints, e.g. regions with large residuals or unrealistic velocity levels.
6
The CFL number is generally defined as CFL=∆t (v+a) / ∆x, where ∆t is the time step, ∆x the local cell size, v the local velocity and a the sonic velocity. For incompressible flows computed with an incompressible solver, the CFL number is defined by CFL=∆t v / ∆x. © ERCOFTAC 2000
page 86 of 94
•
¾
When the residuals stagnate at too high a value, try to further converge the problem with a double precision version of the flow solver. • Check whether switching from a steady-state to a time-accurate calculation has any effect. • Consider using a different initial condition for the calculation. • Check the numerical and physical suitability of the boundary conditions • Examine the choice of variable on which the residual is based, for example density may be a perfectly reasonable choice for transonic flow calculations, but is certainly not for low speed subsonic flows. Guidelines and recommendations to the code developers: • • • • • •
• • • • • • •
CFD codes should make available the maximum possible information to judge convergence. This includes residuals for every conserved quantity. The residuals should be dimensionless and the definitions of the residuals should be clearly defined in the handbook. The codes should supply information on the spatial distribution of residuals. To avoid the confusion of CFD users, CFD code developers should try to adopt one commonly accepted standard definition of the residual. Round-off errors In situations where the small arithmetical differences between two large numbers become relevant, cancellation due to round-off may lead to severe errors. To avoid a negative influence of round-off errors, it is common practice to calculate pressure relative to a reference value. Some codes also use a residual form of the governing equations to combat round-off errors. Examples where round-off errors are known to be of significance are: Low Reynolds number turbulence models with large exponential terms. Flows with density driven buoyant forces with small density and temperature differences. High aspect ratio grids with large area ratios on different sides of the grid. Conjugate heat transfer. Calculations of scalar diffusion with low concentrations of one species. Low Mach number flows with a density based solver. Flows with large hydrostatic pressure gradients.
11.10.2. Guidelines on round off errors Check the influence of round-off errors by performing the same calculation with a single precision and a double precision version of the flow solver, and compare target parameters. ¾ Always use the 64-bit representation of real numbers (double precision on common UNIX workstations), where possible. ¾ Developers are recommended to use the 64-bit representation of real numbers (REAL*8 in FORTRAN) as the default settings for their CFD code. ¾
11.10.3. Guidelines on spatial discretisation Avoid the use of first order schemes. The use of methods of higher order (at least second) is recommended for all transported quantities. ¾ First order schemes may be acceptable under the following circumstances. It may be necessary to use a first order scheme at the start of a calculation as it is likely to be more robust, but as convergence is approached a second order or higher scheme should be used. From a practical point of view it will sometimes be acceptable to use a first order scheme for the transported turbulence variables (which are often dominated by local production and dissipation). For large engineering calculations a second order scheme may be impractical because of the available computing resources. ¾ Try to give an estimate of the discretisation error in the simulation by applying a mesh refinement study (or if this is not possible by mesh coarsening) and by examining the effect of changes caused by the use of different order of spatial discretisations on the same mesh. ¾ If available in the code, make use of the calculation of an error estimator (which may be based on residuals or on the difference between two solutions of different order of accuracy or on wiggle size for a high order numerical scheme).
¾
© ERCOFTAC 2000
page 87 of 94
11.10.4. Guidelines on temporal discretisation Second order accuracy is recommended in both space and time, as the overall solution accuracy is determined by the lower order component of the discretisation. For time dependent flows the time and space discretisation errors are strongly coupled. Hence finer grids, smaller time steps or higher order schemes are required (in both space and time). ¾ Check the influence of the order of the temporal discretisation by analysis of the frequency and time-development of a quantity of interest (e.g. the velocity in the main flow direction). ¾ Check the influence of the time-step on the results. ¾ Ensure that the time-step is adapted to the choice of the grid and the requested temporal size by resolving the frequency of the realistic flow and ensure that it complies with specific stability requirements. ¾
11.11. Guidelines on assessment of errors 11.11.1. Guidelines on user errors in code handling ¾
Try to minimise errors in the handling of the code by the use of a formal check list or by letting another CFD analyst check through the code input data. This will minimise the potential source of user errors when implementing the solution strategy with a particular code. The types of questions which should be considered are: • • • • • • • • • • • • • • •
Is the correct version of the code being used? The default parameter settings or alternatives should be reviewed. Have default parameters been changed which may affect the solution? Are the user defined routines or other software coupled to the code compatible with the current version of the code? Have the boundary conditions not only been properly defined, but also properly applied? Has the appropriate system of units been used? Is the geometry correct? Are all parts of the mesh correctly matched and not overlapping, especially after importing grid parts from different sources? Are the correct physical properties specified? Have the intended physical and mathematical models been used (e.g. gravity forces, rotation, user defined functions)? Have the specified models been changed from those in an earlier code version? Have any user-defined subroutines been used, and have these been checked with regard to programming errors, units, and references to pre-defined variables? Has the appropriate convergence criterion been defined and used? Is there a reasonable maximum iteration number defined, especially in connection with a cputime limit for batch jobs? Have any warnings been generated by the software, and have they been properly understood and acted upon? Have the necessary files from other computer systems been transferred correctly?
11.11.2. Guidelines on use of sensitivity tests Before using a CFD code and model for any specific application, carry out and document tests to determine the sensitivity of the results to the key parameter settings. The key issues to be examined are sensitivity to the grid, discretisation scheme, turbulence model and boundary conditions. Sensitivity analyses can be very time-consuming and this means that in complex cases they may be limited in scope. Because of this, it may be satisfactory to carry out sensitivity tests not on each particular application, but for similar classes of applications requiring the use of the same models and the same code settings. ¾ Base the sensitivity studies on experience of the fluid dynamics of the case under consideration. Rather than changing parameters arbitrarily, physical understanding of the flow should be used whenever possible to identify the parameters most likely to exhibit sensitivity in the results. ¾ Document clearly the results of a sensitivity analysis, to show to what extent the simulation results depend strongly or not at all on the parameter settings. Unsubstantiated assertions as to the in-
¾
© ERCOFTAC 2000
page 88 of 94
sensitivity of the results to certain factors should not be made, but the argument should be demonstrated by reference to quantitative examples.
11.12. Guidelines on interpretation ¾
Make full use of the wide variety of visual tools available in modern CFD packages to interpret the flow-field that has been calculated. A good physical understanding of the application and experience with similar simulations are the best guide as to which plots to make. Options include: • • •
¾
¾
¾
¾
Contours of parameters on arbitrary cross-sections through the flow. Contours of parameters on wetted surfaces in flow, such as wall pressures. Special plots involving projections for specific applications, such as blade-to-blade or meridional views in turbomachinery blade rows. • Line diagrams of parameters. • Three dimensional plots of streamlines. • Three-dimensional views of iso-surfaces, such as regions of low pressure or high mach number. • Surface streak lines. • Use of identification methods to visualise vortical structures, such as line integral convolution (LIC) • Video sequences of unsteady simulations. Produce plots of parameters that are relevant to the physical understanding of the flow (such as pressure fields, velocities, entropy, losses, wall shear stresses, regions of reverse flow, etc.), and also of parameters that might help in the interpretation of the quality of the simulation (such as the ratio of turbulent to laminar viscosity, grid skewness, grid aspect ratio, residuals, etc.). Don’t be seduced into believing that the solution is correct just because it has converged and produced high-quality colour plots (or even seductive video presentations) of the CFD simulations. Make sure that an elementary interpretation of the flow-field explains the fluid behaviour and that the trends of the flow analysis can be reconciled with a simple view of the flow. Be aware of the possible errors from bugs in the pre- and post-processors (interpolation, extrapolation, averaging, etc.). Make sure that the mean values of engineering parameters derived from the simulation are computed consistently (e.g. mass-average values, area-average values, timeaverage values). Calculation of local and mean engineering parameters with external postprocessing software may be inconsistent with the solution method of the code used (e.g. calculating shear stresses from the velocities, calculating shear stresses using nodal values instead of wall functions). Check that any test data used for comparison with the simulations is also computed in the same way as the data from the simulation. Consider whether the interpretation of the results, and any decisions made on the basis of them, are within the accuracy of the computation.
11.13. Guidelines on documentation Keep good records of the simulation with clear documentation of assumptions, approximations, simplifications, geometry and data sources. ¾ Organise the documentation of the calculations so that another CFD expert can follow what has been done. ¾ Be aware that the level of documentation required depends strongly on the customer's requirements as defined in the problem definition.
¾
11.14. Guidelines on communication with code developer 11.14.1. Guidelines for the code developer and vendor The code developer or code vendor needs to demonstrate that he has applied stringent methods of quality control of the software development and maintenance. ¾ Verification of the code is to be carried out by the code vendor or developer, and he should provide the necessary information on the verification process for the user. ¾
© ERCOFTAC 2000
page 89 of 94
¾
¾
¾ ¾
¾
The code developer or code vendor should maintain and publish a databank of verification test cases that are used for testing. The cases should include simple code verifications tests (e.g. that solutions are independent of co-ordinate systems). The code developer or code vendor should provide documentary evidence of the verification tests that the software has undergone, which should include clearly details of the code options which are used during testing. For all new versions of the code and all new machines a standard set of verification test cases should be repeated. Code vendors and code developers should supply a list of known bugs and errors in each version of the code (hot-line, password secured web-page). This list should demonstrate that the number of bugs reduces as the code matures. The code developers should try to include warning notices and guidance for the user in the output. For example, when basic rules on grid generation (expansion ratios, skew, etc.) are being broken, when important specific default options are being overruled by the code input data or when the near wall grid is inconsistent with the turbulence modelling.
11.14.2. Guidelines for the code user ¾ ¾
¾
¾ ¾
¾
¾
Read the release notes issued by the code supplier, as these are a good source of information on changes to the code to remove existing bugs. Recognise that codes can only be validated and verified for a class of problems involving specific variables. If the application involves an area where the code is not fully verified there is more risk of code errors. Set up a suite of test cases that can be run on new code releases to provide an independent check on the code and to highlight any changes between releases (for example in default parameters). Check that the release version of the code is compatibility with any user defined routines, if applicable. When a code error is suspected, communicate this to the code vendor or developer as soon as possible, especially if no list of known bugs has been published. Other users may then profit from this experience. If the bug is well-known a solution or work-around may already be available. In any communications with the code developer or code vendor about a suspected program error, provide a short concise description of the problem and all the necessary input data files so that the error can be reproduced. In cases where commercial sensitivity precludes this, special arrangements will need to be made. Contact the code developer if the meaning and definition of variables in the manuals are unclear and if support is needed in the implementation of special features.
© ERCOFTAC 2000
page 90 of 94
12. Acknowledgements The material presented in this document is a real team effort and thanks are due to the many people and organisations that have co-operated in its preparation, but in particular to the management and staff of the code developers and code vendors who have encouraged and supported the work, to the academic experts who have corrected some elementary mistakes of the editors, and the many CFD users from industry who have freely shared their experience in the use of CFD. The production of the guidelines was organised by the Fluid Dynamics Laboratory of Sulzer Innotec, Ltd. in a contract issued by the ERCOFTAC SIG (Special Interest Group) on Quality and Trust in Industrial CFD. The major organisational and editorial work in association with the document and its preparation was carried out by the project leader at Sulzer Innotec, Dr. Michael Casey, assisted by Torsten Wintergerste. The first version of the best practice guidelines was produced in May 1999 with major input from the following persons: Dr. G. Scheuerer, Dr. F. Menter, Dr. J. Ferreira and Dr. S. Phillipson of AEA Technology, Dr. F. Mendonca of Computational Dynamics, Dr. M. Rabbitt of British Energy, Prof. M. A. Leschziner of University of London, Prof. I. Castro of the University of Surrey, Dr. C. Lea of HSE, Dr. S. Gilham of WSA, Dr. N. Hakimi of Numeca, Dr. P. Wild of Fluent and Dr. M. Casey, Dr. E Lang, Dr. P. Holbein, Dr. M. Wehrli and T. Wintergerste of Sulzer Innotec. The first version was reviewed and improved at a technical Workshop held in Zurich on June 21 and 22, 1999. The following invited persons were present, giving coverage of both academic and industrial points of view from a range of applications, levels of experience, and European countries: Dr. M. Casey, Dr. E. Lang, Dr. P. Holbein, Dr. M. Wehrli, T. Wintergerste, Miss E. Wiederkehr of Sulzer Innotec, Dr. E. Benz (ABB), Dr. J. Vos (EPFL), Prof. A. Hutton (DERA), Dr. C. Carey (Alstom), Dr. A. Kelsey (HSE), Dr. N. Rhodes (Mott Macdonald), Dr. S. Gilham (WSA), Dr. P. Tucker (NAFEMS), Dr. P. Andreasson (Vattenfall Utveckling), Prof. Ch. Hirsch (VUB), Prof. W. Rodi (Karlsruhe), Dr. F. Klimetzek (Daimler Chrysler Cars), Dr. W. Haase (Daimler Chrysler Aerospace), Dr. F. Archambeau (EDF), Dr. H. Paillere (CEA), Dr. F. Mendonca (CD), Dr. S. Phillipson (AEA), Dr N. Hakimi (Numeca) and Dr. P. Wild (Fluent). Many other persons offered their services for the workshop but could not take part either because numbers were limited or because the date was unsuitable. An interim version based on the input from the workshop was produced by Sulzer Innotec for final editing, with particularly high quality and helpful input to chapter 3 from Dr. Marc Wehrli (Sulzer Innotec) and to chapter 4 from Prof. Tony Hutton (DERA). Helpful revisions to the interim version were received from Dr. P. Andreasson (Vattenfall Utveckling), Dr. C. Lea and Dr. A. Kelsey (HSE), H.-D. Grein (NEC), Dr. P. Holbein (Sulzer Innotec), Dr. H. Pailliere (CEA), Dr. M. Rabbitt (British Energy) Dr P. Drtina and Dr. H. Keck (Sulzer Hydro), Prof. Ch. Hirsch (VUB), Prof. I. Castro (Surrey University), Prof. W. Rodi (Karlsruhe University), Dr. W. Haase (Daimler Chrysler Aerospace), Dr. F. Archambeau and Dr. N. Gay (EDF), Dr. S. Phillipson, Dr. G. Scheuerer and Dr. F. Menter (AEA). A draft of the final version was prepared for the ERCOFTAC SIG meeting in Florence at the end of October 1999. After this meeting, some additional revisions of the chapter on turbulence modelling were received from Prof. Ian Castro, Dr. Georg Scheurer, Prof. Mike Leschziner and Prof. Tony Hutton, and these were incorporated into the final version. The test cases were calculated by the code vendors and code developers whose names appear in the sections of chapter 9. Sulzer Innotec and the ERCOFTAC SIG extend many thanks to the management and staff of the code developers and code vendors who have encouraged and supported the work in this way. Individual thanks go to Dr. F. Mendonca, Dr. F. Archambeau, Dr. H. Pailliere, Dr. M. Rabbitt, Dr. N. Hakimi, Dr. P. Wild, Dr. G. Scheuerer and Dr. J. Stokes
© ERCOFTAC 2000
page 91 of 94
13. References AIAA (1998), “AIAA guide for the verification and validation of computational fluid dynamics simulations”, AIAA G-077-1998. Apsley, D., Chen, W-L., Leschziner, M. & Lien, F-S. (1997), “Non-linear eddy-viscosity modelling of separated flows”, IAHR J. Hydraulic Research, vol. 35, pp. 723-748. Baldwin, W.S. & Lomax, H., (1978), “Thin-layer approximation and algebraic model for separated turbulent flows”, AIAA Paper 87-257. Bradshaw, P. (1994), “Turbulence: the chief outstanding difficulty of our subject”, Experiments in Fluids, vol. 16, pp. 203 – 216. Cebeci, T. & Smith, A.M.O., (1974), “Analysis of turbulent boundary layers”, Series in Appl. Math. & Mech., Vol. XV, Academic Press. Chieng, C.C. & Launder, B.E. (1980), “On the calculation of turbulent heat transfer downstream from an abrupt pipe expansion“, Journal of Heat Transfer, vol. 3, pp. 189-207 Daly B. J., Harlow F. H., (1970), “Transport Equations in Turbulence”, Physics of fluids, Vol 13, n° 11, 2634-2649. Durbin, P. A., (1995), “Separated flow computations with the k-ε-υ2 model” AIAA Journal, Vol 33, No. 4, pp. 659-664. Ekaterinaris, J. A., and Menter, F. R., (1994), “Computation of oscillating airfoil flows with one- and two-equation turbulence models”, AIAA Journal, Vol. 32, No.12, pp. 2359-2365. ERCOFTAC/IAHR, (1999), “Workshop on draft tube flows”, June 20 - 23 Porjus Hydro Power Center, Sweden. Ferreira, J.C. & Scheuerer, G. (1997), “Quality assurance in industrial CFD calculations”, Presentation at the First European Meeting of the ERCOFTAC Special Interest Group on Quality and Trust in Industrial CFD, La Defense, Paris. Ferziger, J.H. & Peric, M. (1997), “Computational methods for fluid dynamics”, Springer, Berlin. Fisher, E.M. & Rhodes, N. (1996), “Uncertainty in computational fluid dynamics”, Proc. Mech. Eng., vol. 210, Part C: Journal of Mech. Eng. Sci., pp. 91 – 94. Fletcher, C.A.J. (1991), “Computational techniques for fluid dynamics, vol. I & II”, Springer, Berlin. Gatski, T.B. & Speziale, C.G. (1993), “On explicit algebraic stress models for complex turbulent flows“, J. Fluid Mech. Vol. 254, pp. 59 – 78. Giles, M. B., (1990), “Non-reflecting boundary conditions for euler equation calculations”, AIAA Journal, Vol. 28, No. 12, pp. 2050-2058. Girault, V. & Raviart, P. A., (1986), “Finite element methods for Navier-Stokes equations. Theory and algortihms”, Springer Series in Computational Mathematics, Vol. 5. Grotjans, H. & Menter, F.R., (1998), “Wall Functions for General Application CFD Codes”, ECCOMAS 98, Fourth European Computational Fluid Dynamics Conference, Athens. Hall, R.C. & Cowan, I.R. (1998), “Modelling of atmospheric dispersion near buildings”, NAFEMS Int. Journal of CFD Case Studies, vol. 1, April 1998, pp. 7 – 18. Hanjalic, K. (1994), “Advanced turbulence closure models: A view of current status and future prospects“, Int. J. Heat and Fluid Flow, pp. 178-203. Hirsch, C. (1991), “Numerical computation of internal and external flows, vol. I & II”, Wiley, New York. Ince, N.Z and Launder, B.E. (1995), “Three-dimensional and heat-loss effects on turbulent flow in a nominally two-dimensional cavity”, Int. J. Heat and Fluid Flow, 16, pp. 171-177. Kato, M. & Launder , B.E. (1993), “The modelling of turbulent flow around stationary and vibrating square cylinders”, Proc. 9th Symp. on Turbulent Shear Flows, Kyoto, Japan. Lamb, M. (1895), “Hydrodynamics”, University Press, Cambridge. Launder B. E., (1989) “Second-moment closure : present...and future”, Intl J. Heat Fluid Flow 10, 282-300. Launder, B.E. & Spalding, D.B. (1972), “Mathematical models of turbulence”, Academic Press. Launder, B.E. & Spalding, D.B. (1974), “The numerical computation of turbulent flow”, Comp. Meth. In Appl. Mech. And Engng., 3, pp 269-289. Launder, B.E. (1984), “Second-moment closure: methodology and practice”, Chapter in Turbulence Models and their Applications, Vol. 2, Collection de la Direction des Etudes et Recherches d’Electricite de France, Eyrolles, Paris. Launder, B.E. (1989), “Second moment closure: Present and future? “, Int. J. Heat and Fluid Flow, 10, p. 282.
© ERCOFTAC 2000
page 92 of 94
Leschziner, M. (1990), “Modelling engineering flows with Reynolds stress turbulence closure“, J. Wind Engineering and Industrial Aerodynamics, vol. 35, no. 21. Leschziner, M. A., (1998), “The computation of turbulent engineering flows with turbulent transport closures”, UMIST Report, TFD/98/03, to be published by Springer. Mayle, R. E., (1991), “The role of laminar turbulent transition in gas turbine engines” ASME Journal of Turbomachinery, Vol.113, pp. 506-537. Menter, F.R. (1993), “Zonal two equation k-ω turbulence models for aerodynamic flows”, AIAA Paper 93-2906. Menter, F.R. (1994a), “Two-equation eddy-viscosity turbulence models for engineering applications”, AIAA Journal, vol. 32, No 8. pp. 1598-1605. Menter, F.R. (1994b), “Eddy viscosity transport equations and their relation to the k-ε model”, NASA-TM-108854. Menter, F.R. (1996), “A comparison of some recent eddy-viscosity turbulence models”, Transactions of the SAME, vol. 118, pp. 514 – 519. Menter, F.R., Grotjans, H., Unger, F., (1997), “Numerical Aspects of Turbulence Modelling for the Reynolds Averaged Navier-Stokes Equations”, VKI Lecture Series 1997-02 on Computations Fluid Dynamics, Rhode Saint Genese. Oden, J.T., Zienkiewicz, O.C., Gallagher, R.H. & Taylor, C. (1975), “Finite elements in Fluids, vol. I & II”, John Wiley and Sons, New York. Patankar, S.V. (1980), “Numerical heat transfer and fluid flow”, McGraw-Hill, New York. Patel, V. C., Rodi, W. & Scheuerer, G. (1985), “Turbulence models for near-wall and low Reynolds number flows: a review”, AIAA Journal, vol. 23, No. 9, pp1308-1319. Rizzi, A. & Vos, J. (1998), “Toward establishing credibility in computational fluid dynamics simulations”, AIAA Journal, vol. 36, no. 5, pp. 668 – 675. Roache, P.J. (1972), “Computational fluid dynamics”, Hermosa Publishers, Albuquerque. Roache, P.J., (1997), “Quantification of Uncertainty in Computational Fluid Dynamics”, Annual Rev. Fluid Mech., Vol. 29, pp.123-160 Roache, P.J. (1998), “Verification and validation in computational science and engineering.”, Hermosa Publishers, Albuquerque. Rodi, W. & Spalding, D.B. (1970), “A two-parameter model of turbulence and its application to free jets”, Wärme- und Stoffübertragung, vol. 3, p. 85. Rodi, W. (1981), “Progress in turbulence modelling for incompressible flows”, AIAA Paper 81-45, St. Louis. Rodi, W., (1991), “Experience with two-layer models combining the k-ε model with a one-equation model near the wall”, AIAA paper 91-0216. Rodi, W. (1993), “Turbulence models and their application in hydraulics”, 3rd edition, Balkema, Rotterdam. Rodi, W. (1997), “Comparison of LES and RANS calculations of the flow around bluff bodies”, Journal of Wind Engineering and Industrial Aerodynamics, vol.69-71, pp. 55-75. Sarkar, S., Erlebacher, G., Hussaini, M.Y. & Kreiss, H. O. (1989), “The analysis and modelling of dilatational terms in compressible turbulence”, ICASE Report 89-79, Univ. Space Research Assoc. Hampton, VA, USA. Sarkar, S. & Speziale, C.G. (1990), “A simple non-linear model for the return to isotropy in turbulence”, Physics of Fluids A, vol. 2, pp. 84-93. Savill, A. M (1996), “One point closures for transition”, chapter 6 in “Turbulence and Transition”, eds P.H. Alfredsson and D. Henningston, Kluwer Academic. Schlichting, H. (1979), “Boundary Layer Theory”, McGraw-Hill. Spalart, P.R. & Allmaras, S.R. (1992), “A one-equation turbulence model for aerodynamic flows”, AIAA Paper 92-0439. Speziale, C.G. (1987a), “Second-order closure models for rotating turbulent flows”, Q. Appl Math., vol. 45, pp. 69-71. Speziale, C.G. (1987b), “On non-linear k-L and k-ε models of turbulence”, Journal of Fluid Mechanics, vol. 178, pp. 459-475. Strazisar, A.J. & Denton, J.D. (1995), “CFD code assessment in turbomachinery: a progress report”, Global Gas Turbine News, May/June 1995, pp. 12 – 14. Tennekes, H. & Lumley, J.L. (1972), “A first course in turbulence”, MIT Press, Cambridge. Latest edition 1994. Wilcox, D.C. (1998), “Turbulence modelling for CFD”, DCW Industries, Inc. Wolfshtein, M.W. (1969), “The velocity and temperature distribution in a one-dimensional flow with turbulence augmentation and pressure gradient”, Int. J. Heat and Mass Transfer, 12, pp 301-312.
© ERCOFTAC 2000
page 93 of 94
Yakhot, V., Orszag, S. A., Tangham, S., Gatski, T.B., and Speciale, C. G., (1992), “Development of turbulence models for shear flows by a double expansion technique”, Physics of Fluids A4 (7). Zeman, O. (1990), “Dilatational dissipation: the concept and application in modelling compressible mixing layers”, Physics of Fluids A, Vol. 2, No. 2, pp 178-188. Zienkiewicz, O.C. & Taylor, R.L. (1991), “The finite element method, Volume 2: Solid and fluid mechanics dynamics and non-linearity”, McGraw-Hill, New York.
© ERCOFTAC 2000
page 94 of 94
