Numerical and physical aspects of large eddy simulation of turbulent flows

(1)

NUMERICAL AND PHYSICAL ASPECTS

OF LARGE EDDY SIMULATION OF

TURBULENT FLOWS

MARCO LAFEBER

Delft University Press

TR diss

1558

(2)

/

x

rl NUMERICAL AND PHYSICAL ASPECTS

V i ' OF LARGE EDDY SIMULATION OF

v>

TURBULENT FLOWS

<J(SL ^

(3)

NUMERICAL AND PHYSICAL ASPECTS

OF LARGE EDDY SIMULATION OF

TURBULENT FLOWS

Proefschrift

Ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus, Prof.dr. J . M . Dirken, in het openbaar te verdedigen

ten overstaan van een commissie aangewezen door

het College van Dekanen

op dinsdag 8 september 1987 te 16.00 uur door

MARCO LAFEBER geboren te Leiden werktuigkundig ingenieur

Delft University Press / 1987

(4)

Dit proefschrift is goedgekeurd door de promotoren:

PROF.DR.IR. G. OOMS

(5)

STELLINGEN

1 De keuze van het differentie schema is van

wezenlijk belang voor succesvol voorspellen van het

gedrag van turbulente stromingen met behulp van de

"large eddy simulation" methode.

2 De immer dalende prijs/prestatie verhouding van

digitale rekenmachines verhoogt in toenemende mate

de aantrekkelijkheid tot het gebruik van

rekenintensieve programma's voor het voorspellen

van het gedrag van turbulente stromingen.

3 Voor engineering doeleinden is het k-G model zo gek

nog niet.

4 "Large eddy simulation" kan in belangrijke mate

bijdragen aan de bestudering van coherente

structuren.

5 Tweefasen stroming door een horizontale pijp kan op

kleine schaal realistisch worden nagebootst.

6 Het bij fabricage onnauwkeurig positioneren van

snijplaatjes op "Polycrystalline Diamond Compact"

boorbeitels kan ongewenste trillingen veroorzaken

tijdens het boren.

7 De ontwikkeling van het bewustzijn van het belang

van veiligheid dient een essentieel onderdeel te

vormen van iedere technische opleiding.

8 Prinsjesdag moet drie weken worden vervroegd.

9 Zij die ambiëren te promoveren naast een volledige

hier geen verband mee houdende dienstbetrekking,

dienen in overweging te nemen dat dit een eentonige

wijze van vrijetijdsbesteding tot gevolg zal

hebben.

10 Met stellingen verzinnen,

Kan men niet vroeg genoeg beginnen;

Want is het proefschrift zelve klaar,

Dan staat je hoofd er niet meer naar.

(6)

(7)

SUMMARY

Basic physical and numerical aspects of large eddy simulation are discussed in this work. The applicability of various numerical methods are examined and results are presented from theoretical considerations as well as from numerical experiments. Problems near solid walls are investigated and a more generally applicable proposal for wall shear stress calculation is presented and tested.

A review of large eddy simulation research, as it has been carried out over the past two decades, revealed that the emphasis lay mainly on improvement of subgrid-scale modelling. Applications were limited to turbulent flow problems of relatively simple nature and geometry. In order to apply large eddy simulation to problems of engineering interest, numerical methods and models that can readily be used for this purpose have to be selected. A choice is made between two different approaches to separate the large from the small scale motions : the filter- and the volume balance procedure. Though it is shown that for second order accurate finite difference schemes both methods result in basically identical discretized equations, preference is given to the volume balance procedure, since it allows more flexibility with respect to the geometry of the field of flow. Also, no fundamental problems arise when applying rectangular or non-equidistant grids. Derivation of formal relations between surface- and volume-averaged quantities allows calculation of truncation errors arising from central and upwind differencing methods as well as from various time advancement methods.

Near solid walls problems arise when velocity gradients have to be calculated that originate from diffusive terms. An empirically based method to calculate these gradients has been generalized, whereby greater flexibility can be obtained with respect to geometry and local conditions of surface roughness and (mean) flow direction. Subgrid-scale modelling near solid walls poses a problem because of the absence of locally isotropic behaviour for the unresolved motions. A method proposed by Schumann to split the subgrid-scale stresses into locally isotropic and inhomogeneous parts is applied. This allows grid elements near solid walls to be relatively coarse, hence limiting both the total number of elements as well as the number of time steps required.

To solve the pressure field the commonly applied fast Poisson solvers cannot be used for non-rectangular flow geometries. The slower but more

(12)

general iterative CG and preconditioned CG procedures do not suffer from this problem. Application of both procedures to laminar and turbulent Poiseuille flow (channel flow problem) demonstrated that they can be used as relatively fast iterative methods.

Furthermore, turbulent channel flow has been used to test the applicability of various numerical methods to large eddy simulation. With a relatively coarse grid randomly chosen initial conditions evolved into fully developed channel flow. Good agreement with measured data was found for various statistical quantities, demonstrating the potential of large eddy simulation as a (future) engineering tool. Application of finer grids leads to distinct increases in computing time requirements but improves the results, in particular near solid walls.

(13)

1 - INTRODUCTION

1 . 1 - Methods of investigating turbulent flow behaviour

Nearly all flows encountered in nature and technology are turbulent. From a practical point of view the study of turbulence is therefore extremely important. A mathematical description of both laminar and turbulent flows is given by well known basic equations of motion. Unfortunately it is highly unlikely that a general solution for these non-linear equations will ever be found. To study the behaviour of turbulent flows, experiments can be carried out by which the quantities of interest are measured or particular flow patterns are observed. However, not all quantities can be measured, e.g. those involving pressure fluctuations. Furthermore, isolation of the effect of a single parameter is usually not possible without affecting other parameters simultaneously. Measurements cannot always be carried out, for instance if the domain of the flow is either too large or too small, hidden inside machinery, or perhaps even radioactive. On the other hand, the fluid itself can be dangerous to deal with, e.g. poisonous, or have extreme temperatures. For such cases the use of a scale model and/or a different fluid can be a solution to allow experiments to be carried out under well defined and controlled conditions. Unfortunately this is not always attractive for economic reasons, since building a scale model and taking accurate measurements by well qualified people is usually very expensive. A more serious argument against the use of scale models for certain flow problems is that they cannot always be scaled down properly (e.g. when temperature gradients are involved). Similarity between the real and the scaled flow then becomes questionable.

The introduction of digital computer systems made attempts possible to simulate flow behaviour by solving the basic equations with numerical methods. In the past decades the popularity of this approach to fluid dynamics grew along with the increase of memory capacity and computational speed. In contrast to laminar flow problems, numerical simulation of turbulent flows cannot be carried out by simply discretizing the governing equations and solving them on a certain grid. This is caused by the fact that turbulence is essentially three dimensional and contains many length scales simultaneously. With increasing Reynolds number the length scales of the smallest eddies in the flow become smaller and smaller. Consequently, the amount of data that has to be handled in an attempt to describe all the

(14)

length scales that occur, increases with the Reynolds number. Even the largest of the modern computers do not (yet?) have the required speed and memory capacity available to handle this amount of data, except for turbulent flows with relatively low Reynolds numbers (sec. 1.2). For approximate computations at high Reynolds numbers to be carried out, various methods can be applied to model the influence of fluctuations on the calculated mean flow field. In the next section several of these methods will be discussed and placed in perspective with respect to the subject of this work : the method of large eddy simulation.

1.2 - Numerical simulation of turbulent flows

A number of books and articles have been published that give excellent reviews and discussions on the impressive number of numerical simulation methods available for turbulent flows (e.g. Reynolds and Cebeci [1976], Launder and Spalding [1972], Bradshaw et al. [1981]). For obvious reasons these methods will not all be discussed here. Instead a comparison will be made between various "levels" at which turbulent flow behaviour is currently simulated.

The basic level comprises the numerical simulation methods that are at present amongst the most popular to solve flow problems for engineering purposes. Their predictions are based on time-averaged properties of turbulence. Since these properties vary only gradually in space, the use of coarse grids can generally provide enough accuracy. Furthermore only a limited number of time steps is usually required to converge from initial conditions to the final-solution. However, a price has to be paid for neglecting the fluctuations : while averaging the basic equations of motion correlations involving fluctuating velocities (and temperatures) appear. They are additional unknowns for which no equations can be derived without again introducing new unknowns. This is called the closure problem of turbulence. Closure models are required to give additional relations for the correlations, so that there are as many equations as unknowns. Numerous attempts to construct such models have been undertaken. Most of them give good results for (the classes of) flow problems for which they have been designed. However, the reliability of the results obtained with such models decreases with increasing complexity of- e.g. boundary conditions.' Experimental data is required to support this basic level of simulation.

(15)

This limits the potential of these methods to replace or assist experiments meant to study the (fundamental) behaviour of turbulence.

The highest level of numerical simulation uses the basic equations of motion without any averaging applied. These so-called direct simulations (DS) (or full turbulent simulations (FTS)) are in principle preferable to any other since no modelling assumptions are required. From a theoretical point of view this level is an excellent replacement for experiments. Even the smallest motions are simulated and detailed information can be obtained for any combination of quantities at any position and time. Unfortunately, the presence of small motions causes practical problems. The characteristic length scale for the smallest motions is given by the Kolmogorov microscale

n (A circumflex ( ) denotes that the variable is dimensional). The relation

between r) and the length scale t of the largest eddies is given by

(Tennekes, Lumley [1972]) :

$/£

- R

-3/4 _{; i . D}

where R is the Reynolds number with respect to £ . If the dimensions of the -, o

mean flow field are of the order t and the sizes ft of the grid elements are

about equal to r\ it can easily be verified that the number of grid elements

LHN, needed to discretize the field of flow, results from eqn (1.1) as :

LMN - R, 9/4 (1.2)

In figure 1.1 the above relation is shown.

\ —

figure 1.1 relation between number of grid

(16)

The dotted part of the line is of no meaning with respect to eqn (1.2) since it lies in the range of low R numbers where the flow will generally be laminar. It may.be clear from this figure that DS methods are only possible at relatively low Reynolds numbers, whereby the required number of grid elements is already extremely large. Even larger is the number of data involved, since at least three velocities and one pressure variable have to be stored for each grid element. The largest of the present generation of super-computers are still by far too small for practical application of DS. This can already be concluded even without discussing the calculation time that would be required to simulate a representative time interval for such an instationary flow problem. Furthermore, the question has been raised whether all of the detailed information produced by a direct simulation is relevant for other than

At an intermediate

theoretical purposes.

level a method is under development that needs a closure model, but a much more universal applicable one than those used at the basic level. This method is called large eddy simulation (LES). As the name already indicates, only the large eddies (or resolvable scale motions) are calculated, wherea's the small eddies (or subgrid-scale motions) are modelled by a closure assumption. The distinction between large and small eddies is typically made at length scales of the order of magnitude of the grid elements. Ferziger and Leslie [1979] summarize the properties of the large and the small eddies as shown in table 1.1.

large eddies

produced from mean flow

transporting (mass, momentum, energy) flow dependent

anisotropic long lived hard to model

small eddies

produced by large eddies dissipating

~ universal ~ isotropic short lived easier to model

table 1.1 - properties of large and small eddies (from Ferziger and Leslie [1979])

The flow dependent large eddies are directly affected by the boundary conditions. They are therefore the most difficult ones to model with methods of the basic level mentioned above. LES avoids this problem by explicit computation of these motions. Since tne small eddies dissipate energy from

(17)

large ones, a so-called subgrid-scale (SGS) model is needed that has to account for this physical energy cascade process. The sizes of the grid elements, and thereby the distinction between large and small eddies, has to be chosen to be small enough for the unresolved subgrid-scale motions to behave statistically in a nearly isotropic manner (sec. 1.3). If this is the case the subgrid-scale motions can be modelled independently of the flow geometry. The parameters involved will then become real constants.

Large eddy simulation is not only an intermediate level method as far as modelling is concerned, but also with respect to required data storage capacity and computing time. Compared to the basic level methods more grid elements are needed to separate the small from the large eddies. Furthermore, the calculation has to be time-dependent, with time steps small enough to calculate the motions of even the smallest resolvable scale eddies. Consequently, large numbers of grid elements and time steps will be required to reach statistically steady state solutions. Fortunately these numbers are by far not as large as required for direct simulations. Mainframe computers allow extensive LES experiments, whereby the present capacities of the largest allow simulations using grids with over half a million of elements (sec. 1.4).

Since LES lies intermediate between the basic and the highest level of turbulent flow simulation, it has been used to provide valuable information to adjust model parameters for basic level methods (Schumann [1973]). On the other hand DS calculations have been used to test subgrid-scale models (Clark et al. [1977], McMillan and Ferziger [1979]).

1.3 - Basic aspects of large eddy simulations

The aspects that basically need attention in order to perform a large eddy simulation stem from various origins (table 1.2). As far as theory is concerned, fluid mechanics provides the equations that describe the large scale motions as well as the model for the energy drain caused by the (unresolved) small scale motions. With respect to this model, fluid mechanics can give an estimate for the grid element sizes in relation to the Reynolds number of the flow. Since nearly isotropic properties are required for the subgrid-scale eddies, the sizes of the grid elements have to be

(18)

theory software hardware comparative data ,, fluid mechanics \ numerical methods s LES program ^ plot package * "small" computer ^ mainframe computer

y DS or fine(r)-grid LES calculations

\

^ experiments

table 1.2 - basic aspects of large eddy simulation

small enough for the unresolved eddies to behave as so-called locally isotropic turbulence (see e.g. Monin and Yaglom [1979] for a detailed discussion of this concept). The adjective "local" refers to the small scales whose strain-rate field itself may be expected to be isotropic in the mean (Tennekes and Lumley [1972]). For a turbulent flow to have local isotropy at the smallest scales the Reynolds number has to be large enough. Schumann [1973] recommends R >2*10 and grid side length H of the order L /30 for the unresolved motions to behave locally isotropic (except for wall adjacent grid cells). However, results reported by e.g. Schumann [1973,1975a], Kwak et al. [1975] and Horiuti and Kuwahara [1982] as well as calculations presented in chapter 6 of the present work show that coarser grids can also give good results.

The software required for large eddy simulation must of course contain discretized equivalents of the basic equations and models as provided by fluid mechanics. The numerical methods used to discretize and solve these equations are of essential importance. Successful large eddy simulations can only be carried out if the numerical methods allow the solutions to be obtained without being disturbed by numerical instabilities or large truncation errors. For certain discretization methods the errors can become too large to be acceptable, e.g. when so-called artificial viscosity effects overshadow the subgrid-scale model terms (sec. 3.5).

The result of a large eddy simulation is not a single "steady state" flow field. Instead, series of successive velocity and pressure vectors are calculated. They are perhaps best presented by a computer generated movie showing the eddies of the resolved flow field (or any desired quantity or

(19)

correlation) changing as time progresses. Since this presentation technique cannot be used for written reports a plot package is required by which results of a single time step as well as the evolution of calculated statistical quantities can be shown.

Both the large eddy simulation program as well as the plot package used in this work have been developed on a relatively small computer (PDP 11/34). Once certain parts of the software had been completed they were transferred to an IBM 3083-JX1 mainframe computer. This was found to be efficient because of the easy and immediate access to the small computer. Interactive programming and immediate availability of printed and plotted results lightened the often laborious testing and debugging work. At present only super-computers like e.g. CYBER 205 or CRAY I allow relatively fast LES calculations with grid side lengths within the previously mentioned order L /30. Fortunately slower mainframe computers, like the IBM 3083-JX1 used for this work, can already produce promising results and allow numerous valuable tests on models and (numerical) methods.

For quantitative testing of the numerical results laboratory experiments can be simulated. Measured quantities can be used for comparison with those retrieved from calculated velocity (and pressure) fields. Qualitative tests can be carried out by using (proven) reliable results from either direct simulations or fine(r)-grid large eddy simulations. The influence of different parameter values or model changes can then be evaluated against data with smaller numerical errors and, in case of DS, without errors caused by closure models.

1.4 - Review of large eddy simulation research

In this section a review is given of large eddy simulation research from its earliest beginning till the present state of the art. The aim is primarily to present various flow types, geometries, models and basic approaches that have been used so far. It is neither intended to be complete in every detail nor to mention every article or paper that has been published on the subject.

The history of LES goes back to the early sixties when the meteorologists Lilly [1962] and Smagorinsky [1963] introduced the idea of the eddy viscosity to be a function of local large scale dynamics. Lilly [1967] is mentioned to be the first to present explicit use of grid averaging, and

(20)

thereby the separation of eddies into resolvable- and subgrid-scale ones (Herring [1979]). Deardorff [1970] tested this (so far only meteorological) approach to what he called "an interesting case of laboratory turbulence : plane Poiseuille flow (channel flow) driven by a uniform pressure gradient". With 6,720 grid elements (24*14*20) he achieved results that agreed well with experimental data of Laufer [1950] and even better with those more recently obtained by Comte-Bellot [1965]. Subsequently he extended the turbulence model to the study of neutral and unstable planetary boundary layers using up to 32,000 (40*40*20) grid elements (Deardorff [1972]). Temperature was introduced as an additional variable. When treating a stable layer overlying a buoyantly active convective region he found that calculation of the eddy viscosity in the usual way, by employing the so-called Smagorinsky model, produced unrealistic smoothing of the vertical temperature field in the transition layer (Deardorff [1972a]). Instead of the Smagorinsky model, subgrid transport equations were used that permitted much improved treatment of the motions in the aforementioned transition layer (Deardorff [1973]). The price to be paid was an increase in calculation time and required computer memory, as well as more complex models that needed more closure constants. Deardorff concluded that adverse effects of truncation errors appear to be more serious than uncertainties in assumptions made by the subgrid transport equations modelling. He therefore suggested that the second order energy conserving scheme be replaced by a fourth order differencing scheme for the convective terms.

The procedures used by Deardorff [1970] were improved by Schumann [1973,1975a]. His volume balance procedure uses conservative grid volume averages of the equations of motion, averaging over grid volumes fixed in space (sec. 2.5). Apart from computation of Deardorff's [1970] channel flow problem, annular flows could be simulated because grid volume averaging allowed the use of anisotropic grids. Introduction of an inhomogeneous part next to a locally isotropic part for the subgrid-scale stresses made treatment of flows with non-zero time-mean strain components possible. This proved to be of importance near solid walls where it allowed the use of coarse grid resolutions. The maximum number of elements used was 65,536 for

it

a 64*32*32 grid. Grotzbach [1977] further extended Schumann's work at Karlsruhe by including temperature as an additional variable. He also improved modelling to allow locally different surface roughness, which enabled him to study the effect of secondary flows in a channel.

At Stanford University a program of development and exploration of LES was started at 1973 (Reynolds et al. [1976]). The objective was to put a

(21)

sound foundation under the method by computing simple flows first. The first (theoretical) contribution was made by Leonard [1974]. He presented a conceptual framework to mathematically separate the large from the small eddies by means of a convolution filter (sec. 2.3). Unlike Schumann's grid volume averaging, such filtering does not satisfy some of the Reynolds

conditions (Monin and Yaglom [1979]), like uu=uu and uu'=0 (here ( ) denotes a filtered variable and ( ) ' a subgrid-scale one). The remaining additional non-zero terms (usually referred to as the Leonard stresses and cross term stresses) were in theory found to be responsible for a non-negligible part of the total dissipation (Leonard [1974]). The first flows considered by the Stanford group were simple decaying isotropic turbulence and homogeneous

3 3

turbulence (Kwak et al. [1975]). Both with 16 and 32 grid elements these flows could be simulated using models that proved to be independent of size and number of grids elements. Confidence in the validity of the concept of LES was thereby provided. Results from direct simulations with 64 grid elements were used by Clark et al. [1979] and McMillan and Ferziger [1979]

3 3

to test subgrid-scale models on 8 respectively 16 grids for homogeneous turbulence. From these tests the non-negligible influence of the Leonard-and cross term stresses was confirmed when applying fourth order differencing schemes (Clark et al. [1977]). Next the channel flow problem as mentioned in the foregoing was given major attention by the Stanford group

(Moin et al. [1978], Kim and Moin [1980]). Up to 516,096 (63*64*128) grid elements have been used by Moin and Kim [1982] for a numerical investigation of this flow problem. Particular attention was given to flow structures in the vicinity of the wall. The results obtained were in very good agreement with measurements. They could even be used as a data-basis for detailed investigations of bursting processes and hairpin vortices (e.g. Moin and Kim [1985] and Kim and Moin [1986]).

At Queen Mary College in the UK large eddy simulation has been subject of investigation since 1975. The QMC group first studied subgrid-scale models. Isotropic flow was used to investigate the dependence of the model on filter shape, grid-scale spectrum and grid anisotropy (Leslie and Quarini [1979]). Love and Leslie [1979] and Love [1980] report results on subgrid modelling studies using Burgers' equation. They found the standard Smagorinsky model with or without the Leonard variant to be superior to more complex models derived from it. Antonopoulos-Domis [1981] examined the concept of filtering

3 3

for homogeneous isotropic turbulence on 16 and 32 grid elements. It was found that Schumann's volume balance procedure is capable of producing good

(22)

agreement between simulation and laboratory experiments. Inclusion of the Leonard term was not found to be necessary. Furthermore the SGS effects on the resolvable scale motions can successfully be modelled by the Smagorinsky model and the model parameter was found to be independent of the mesh size. The QMC group therefore prefered to use averaging methods such as the volume balance procedure, which makes the Leonard (and cross term) stresses automatically zero (Leslie [1982]). Large eddy simulation of a passive scalar in isotropic turbulence was investigated by Antonopoulos-Domis [1981a]. The results were in good agreement with experiments, thereby confirming the aforementioned conclusions on averaging and modelling.

Schumann's volume balance procedure offers the possibility to use anisotropic grids. This was utilized by Friedrich and Su [1982] and Schmitt and Friedrich [1984] who applied curvilinear coordinates to study

boundary-layer flow over a curved wall with up to 65,536 (32*32*64) grid elements. One of the French contributions to large eddy simulation was made by Baron [1982]. He simulated the channel flow problem with 65,536 grid elements (64*32*32) as well as homogeneous turbulence maintained by a

3

uniform mean shear with a 32 grid. Baron also studied turbulent flow induced by a jet in a cavity using 83,200 elements (52*40*40). At present this is geometrically the most complex problem to which large eddy simulation techniques have been applied.

The channel flow problem also received attention in Japan. Horiuti and Kuwahara [1982] used a rather coarse grid with 5,376 elements (16*16*21). Nevertheless their results were found to be in very good agreement with experimental data. Theoretical work on LES was carried out by Yoshizawa [1982]. He derived the Smagorinsky model from a statistical viewpoint under the assumption of small convection compared to production and dissipation in the SGS turbulent energy transport equation.

Mason and Callen [1986] used turbulent channel flow simulations for studies on the magnitude of the subgrid-scale eddy coefficient (Smagorinsky constant).

Application of large eddy simulation techniques to the type of flow problems that actually gave birth to the method were reported by the Dutch meteorologists Nieuwstadt and Brost [1986] and Nieuwstadt, Brost and van Stijn [1986]. They studied the decay of turbulence in a convective boundary layer in which the surface heat flux is suddenly turned off. This type of flow, which occurs for example around sunset, is difficult to study directly in the atmosphere as well as in the laboratory. A large eddy model proved to be able to simulate this decay realistically.

(23)

Other non-isothermal research was carried out by Eidsen [1985] in the US. He successfully simulated turbulent natural convection (the Rayleigh-Bénard problem).

1.5 - M m of the present work

As can be concluded from the historical review presented in the previous section, experience in large eddy simulation techniques has gradually been gained over the past 20 years. Major contributions originate from research carried out in the US, West-Germany and the UK. In the early eighties, when the present large eddy research project was started at the Delft University of Technology, no known effort had been made in the Netherlands to catch up with this promising method for turbulent flow simulation. The initial intention was to use LES as a method for prediction of flow behaviour in climatised rooms. However during the course of investigation it became apparent that it is essential to study relatively "simple" geometries first in order to recognize and understand the specific physical and numerical problems of large eddy simulation. The problems of physical nature are mainly related to models that need to be generalized. This to make them applicable to flows of a more complicated geometrical structure than that of a channel with a distinct mean flow direction. The numerical problems are also strongly related to this increased complexity. The boundary conditions for instance may be much more difficult to define in such a way that they are physically realistic without causing numerical instabilities. The pressure solution method is of iterative nature and bound to become a very time consuming part of the calculation. A method that solves the pressure field accurately and quickly in a general geometry is required for most engineering flow problems.

The aim of the present work can be specified as : "To investigate numerical and physical aspects of large eddy simulations of turbulent flows, thereby focussing on :

- averaging methods used to derive the basic equations of LES, their practical applicabilities and mutual relations

- the applicability of various numerical methods for discretization of the basic equations of LES and their influence on SGS modelling

- physical and numerical problems near solid walls and generalization of wall shear stress computation

(24)

- fast and accurate calculation of velocity and pressure fields, using numerical methods that can be applied in a general geometry

- validation of theory by numerical experiments."

A computer program has been developed that allows testing of models and numerical methods. Though the conceptual framework for this program allows simulations of flows enclosed by more complex boundaries than for a channel, the present investigation will be focussed on the latter problem. As may have been concluded from the previous section, this problem already carries a considerable amount of difficulties with it that are typical for large eddy simulations. Furthermore, it is an excellent problem to study when starting a large eddy simulation investigation project, for the results obtained with new models and (numerical) methods can be compared to results already available from both empirical and numerical origin.

(25)

2 - BASIC EQUATIONS OF LARGE EDDY SIMULATION

2.1 - Introduction

As was already pointed out in the previous chapter, the small scale motions have to be separated from the large scale ones. For this purpose two approaches are in use : the so-called "filter procedure" and the "volume balance procedure". Both approaches will be discussed in this chapter and the basic equations for large eddy simulation that result after application will be compared. The small scale motions that cannot be computed explicitly have to be modelled. Several models are discussed and a choice is made for one of them that is used to perform the calculations reported in chapter 6.

The following symbols are applied for the variables present in the basic equations. Time is denoted by t, space coordinates by either x,y,z or x ,x ,x and velocities by u,v,w or u ,u ,u . Pressure is represented by p. A Cartesian coordinate system is used throughout this work unless stated otherwise. Einstein's summation convention will be applied to repeated lower indices i,j,k=l,2,3 only (so not to upper indices ! ) .

(26)

2.2 - Basic equations of motion

For isothermal flows of isotropic, incompressible Newtonian fluids with constant physical properties, the basic equations of motion (Monin, Yaglom [1979]) are the equation of continuity (a circumflex ( ) denotes that the variable is dimensional) :

ail

—A = 0 (2.1)

3x.

i

and the Navier-Stokes (NS) equations :

,A , 2 A 3u. , 3 u. ,A A 1 A 0 , A A , A X O P / . -, r, ^ \ / -, -, \ p = - p (u.U.) + p. - —*- (i = l,2,3) (2.2) O ,A * 0 ,A 1 D O -A ,A ,A 9t 3x . 3x .3x . 3x. 3 D 3 1

Together with appropriate boundary conditions these four equations are sufficient to describe both laminar and turbulent flows in three space dimensions and in time.

Dimensionless equivalents of the above equations can be obtained with the aid of the following scales :

length scale : L _o velocity scale : u (2.3) time scale : t = L /u o o o A A A 2 pressure scale : p = p u o *o o

The choice for the length- and velocity scales L and u will be made such o o

that they are characteristic for the flow problem considered (see sec. 6.2 and 6.3.1). The dimensionless basic equations of motion read :

3u. ^ = 0 (2.4) i 3u. a2u. T Ti = - ^ - ( u . u . ) + l - ^ - f^- (i=l,2,3) (2.5) 3t 3x . l ] R 3x . 3x . 3x. 3 3 3

1-where R is the Reynolds number, expressed in terms of the characteristic scales as :

(27)

R = -2-2 (2.6)

A

V O

The kinematic viscosity v is a property of the fluid.

2.3 - Averaging procedures applied in large eddy simulation

Two approaches are used by LES researchers to separate the small scale from the large scale motions. In the "filter procedure" proposed by Lilly [1967] and unified by Leonard [1974], the large scale (or resolvable scale) field a(x,t) is defined by means of a convolution filter :

a(x,t) = ƒ G(x-x") a(x",t) dx" (2.7) flow

volume

Several filter functions G have been used. Most of the earlier workers (e.g. Lilly [1967] and Deardorff [1970-1973]) have used a box (or top-hat) filter:

G(x)

1/&1 at |x-x | < §A

A 1 ^ A ( 2_8 )

0 otherwise

The filter width A. is of the size h of the grid elements. The Stanford A

group (e.g. Kwak et al. [1975], Ferziger [1977,1981], Mansour et al. [1979]) prefer to use a Gaussian filter :

G(x) = ( - ^ )3 / 2 exp (-7x2/A2) (2.9)

"AA

with 7=6 and A»-2h. According to McMillan et al. [1979] the choice of the filter has only a minor influence on the results of the LES calculations and is "a matter of the user's preference".

In contrast to the filter procedure the "volume balance procedure" of Schumann [1973] obtains the resolvable variables by averaging over volumes that are fixed in space. This results in both surface and volume averaged variables, defined at fixed locations in the grid (as will be shown in section 2 . 5 ) . Sometimes (e.g. Antonopoulos-Domis [1981]) the volume balance procedure is called the "filter procedure" whereas what is here (and more

(28)

commonly) named "filter procedure" is referred to as "prefilter procedure". The volume balance procedure has mainly been used by the Karlsruhe group (e.g. Schumann [1973,1975a,1980], Grotzbach [1977,1979]). In the present work this procedure will also be applied.

In section 2.8 a relation will be shown between the final sets of discretized basic equations that result from both averaging procedures. For this purpose the averaged basic equations that result after application of either of the two procedures will be needed; they are derived in sections 2.5 and 2.6. Since the derivation of the volume averaged basic equations requires more detailed information about the finite difference grid system that is used, this will first be discussed in the next section.

2.4 - Finite difference grid system

The finite difference grid system employed in this work is the so-called "staggered grid" (Harlow, Welch [1965]). In this grid the velocities are calculated for points located on the faces of the control volumes, whereas pressure is attached to points at the centre. Figure 2.1 shows the locations for u, v, w and p in an element of a (Cartesian) staggered grid.

figure 2.1 - staggered grid element

An advantage of the use of this grid is that no approximations are necessary for discretization of the continuity equation, as will be demonstrated in section 2.5. Furthermore, the staggered grid leads to a particularly simple discretization of the Poisson equation that needs to be solved to obtain the pressure field (sec. 5.3). No spatial pressure oscillations occur as found with non-staggered grids (see Patankar [1980] for an extensive discussion of

(29)

the staggered grid and its particular advantages). This is very important for LES calculations, where sources for "non-physical" oscillations need to be avoided, since they can at least disturb or even gradually destabilize a complete calculation. Consequently, it is not surprising that the staggered grid is very popular in large eddy simulations and is, with a few exceptions, always employed.

2.5 - Volume balance procedure applied to the basic equations

In order to examine Schumann's volume balance procedure a (rectangular) grid element V is considered. This element is fixed in the staggered grid as shown in figure 2.2 and has a content Av=AxAyAz.

figure 2.2 - volume V within staggered grid

In general the volume average of a continuous function a(x,t) for volume v, with coordinates of its corner points as shown in figure 2.2, is defined by:

a(t) = 1 z yt x n t e AxAyAz ! I I a(x,t) dx dy dz (2.10) z yw x s 'b w

Applying this definition to a partial derivative (with respect to one of the space coordinates), a first order centered difference quotient results. By way of example this is demonstrated for 3a(x,t)/3x :

(30)

z y x

^

a (

**5 < t ) = ^ z " '" ' '" i ï * ^ »**

d x d y dz ZS yD XW Zn yt ^ ^ ƒ ƒ [ a ( xe, y , z , t ) - a ( xw, y , z , t ) ] dy d z (2.11) S "'t)

The surface average of function a ( x , t ) for a surface A = A y A z w i t h normal in x-direction and contained in the aforementioned volume V is defined by :

z y^ _n

Jt

xa(x,t) = T^TZ I I a(x,t) dy dz _AyAz (2.12) s b

With this definition eqn (2.11) can be written as :

iJ

a ( X

' ° = to

[

" e

( t )

" %

( t ) ]

=

6

x "

( t )

(2.13)

Applying the volume balance procedure, with respect to a volume V that coincides with a staggered grid element, to the partial derivatives of the continuity equation (2.4), results in :

6. u. (2.14)

No approximations are necessary for this equation if the (surface average) velocities 1u are taken to be equal to the components in x. direction as defined in the staggered grid (sec. 2.4). From here on the velocities u stored in the staggered grid will therefore be defined as surface average velocities. In figure 2.3 an example is shown for velocity u (t) representing the average value of u(x ,y,z,t) for surface A.

^r^

" A - 2 . ^ X

0 > >

X*»Xp

(31)

Note that the surface averaged velocities u (x ,t) (i,j=l,2,3) are D i

continuous functions of x. and t only.

Let us now turn to the NS equations (2.5), and apply the volume balance procedure as defined by eqn (2.10). The volumes V. used for the averaging differ from the volume V that was used for the continuity equation. For the i'th NS equation they are staggered in x.-direction with respect to V, so as to surround the grid point for the velocity component in this direction. The pressure points now lie on the faces of V. with normal in x.-direction.

l i

Figure 2.4 shows the staggered grid volumes V , V and V . x y z VY V2

-O

—

w / 1

y^~

'

o

V

y

s

yf

y

e.

hi

v

x

y

o

"y

y

,y

-»-x

figure 2.4 - staggered grid volumes V , V and V x y z

By applying eqns (2.10) and (2.13) to eqn (2.5) the volume averaged NS equations can easily be derived as :

3u. . : 3u.

—A = - ó.^u.u.) + % o. (T-2-) - S / p

at D I D R 3 OX . l i=l,2,3) (2.15)

x. (i=l,2,3) and t are independent variables and the surface averaged . velocity gradients with respect to x . in the diffusive terms involve

integrations with respect to coordinates x. and x whereby i.k^j. Consequently eqn (2.15) can be written as :

V.

a

x

u .

1 3t n

3

j

ü.

6.3(u.u.) + r 6. - —i - fi/p D I D R D 3 X . i (i=l,2,3) (2.16)

(32)

No further approximations were necessary for the continuity equation (2.14) that resulted after applying the volume balance procedure. This is unfortunately not the case with equations (2.16). Their terms need further approximations, in order to express them in terms of the averaged quantities stored in the staggered grid. Some of these approximations are straightforward and will be discussed in chapter 3. Others however need models based on physical reasoning and empirical information. One of the latter type of approximations concerns the non-linear convective terms. They can only partly be expressed in terms of surface averaged velocities. This is done by separating the velocities u.(x,t) into a part 3u.(x.,t) that is

the average for surface A and a part u!(x,t) that represents the fluctuations with respect to u. for that surface A :

u.(x.t) = 3ü.(x.,t) + -V(x,t) x e 3A (2.17)

l - l 3 l

-In this way the fluctuations u!(x,t) constitute the subgrid-scale motions at very specific and well defined surfaces A in the discretized flow field. The convective terms in eqn (2.16) may now be written as :

3 & J(\x.u.) = 6 . (^ü. + V ) ( V + V . ) 3 1 ] 3 1 1 ] 3 3 ._ ■_ ] ~ : ] . ._ 3 . = 6. ( V - V + 3u Du' + V - V + W . ) (2.18) ] i ] 1 3 1 3 1 3

From the definition of surface averaging over a surface A fixed in space (eqn (2.12)), it is easy to prove that :

■V = ju. and ju! = 0 (2.19)

1 1 1

It can also easily be shown that :

V ju . = V ju. = 0 (i,j=l,2,3) (2.20)

With eqns (2.19) and (2.20) eqn (2.18) can be simplified to :

6 .j( u . u . ) = 6 Vju . + 6. Vju ' . (2.21)

(33)

The last term on the right-hand side of the above equation contains subgrid-scale velocities. It can only be approximated with the aid of models based on physical reasoning and empirical information. This will be discussed in section 2.7. Another approximation of this type is required for the velocity gradients 3 u./3x. that occur in the diffusive terms in case they have to be calculated for surfaces that coincide with solid walls. Section 4.2 is devoted to this problem.

By substitution of eqn (2.21) into eqn (2.16) the volume averaged NS equations are obtained :

V

3 u. ._ ._ 3- : 3 u.

—rr

3

- = - S.fV-V) - 6. W

. +

z 6. - — - - b.

X

p (i=l,2,3) (2.22)

ot 3 1 3 3 1 3 R 3 3 x . 1

2.6 - Filter procedure applied to the basic equations The filter procedure defines the resolvable field by :

u = ü + u' (2.23) The overbar denotes the filter defined by eqn (2.7), separating the large

scale motions from the subgrid-scale ones that are represented by u'. Substitution of eqn (2.23) into the basic equations of motion (2.4) and (2.5) and applying the filter (2.7) to the resulting terms in order to eliminate as much of the subgrid-scale variables as possible, results in :

du. ^ = 0 (2.24) 1 2-3u- a ~ - 3 u- a" TT^ = - TT- K_{ot 3x . 1 1} + "')<"■_{1 i R ox . 3x . ox.} + "'■) + Ï a, at " l^T (1=1,2,3) (2.25) J D 3 1

In the same way as for eqn (2.18) the convective terms C. . can be decomposed

3

C. . = - — (u. + u!)(u. + u'. ) = - — ( u.u. + u.u'. + u!u. + u!u'. ) (2.26)

(34)

Leonard [1974] first drew attention to the fact that none of the terms on the right-hand side of eqn (2.26) can generally be omitted. Only if the u remain constant over a certain averaging volume or surface (as is the case with the volume balance procedure) one can prove that :

u.u = u.u. = 0 and u.u. = u.u. (2.27)

1 3 X 3 l 3 l 3

Usually the filtered convective terms C. . are written as :

ID

C. . = - ê - (Ü.Ü. + L. . + T. .) (2.28) i] 3x 1 3 l] ID

where u u is the true resolvable scale inertial term, while :

'1 D

L. . = ü.u . - u.u . (2.29)

13 I D I D

is known as the "Leonard stress". Finally :

T . = u u ' + u ! u . + u ' u ' . (2.30) ID i j I D I D

is the subgrid-scale stress. Leonard [1974] reports that the L. .-terms should not be neglected. This has been confirmed by research carried out by a.o. Clark et al. [1979] and Moin et al. [1982]. The same applies to the first and second term on the right-hand side of (2.30) (Clark et al. [1979]). Kwak et al. [1975] derived approximations for the latter terms, usually called "cross term stresses" CR , as well as for the Leonard stress

ID L . . : ID CR. . = u . u'. + u : u . = - r A? (u V2 u . + ü . V2 ü . ) ( 2 . 3 1 ) ID i ] i 3 2 4 A l D D 1 L . . = u . u . - ü . u . ~ ^ A ^ V2 ( u . u . ) ' ( 2 . 3 2 ) ID I D 1 3 2 4 A I D 2

(where V is the Laplace operator). Both approximations have frequently been used in conjunction with the filter procedure.

(35)

2.7 - Subgrid-scale modelling

As was already mentioned in section 2.5 approximations are required for some of the terms of the averaged NS equations (2.22) or (2.25). To these terms belong the subgrid-scale stresses; they can only be calculated with models based on physical reasoning and empirical information. For the volume balance procedure the SGS stresses read :

jr. . = ju:ju'. (2.33a)

I D i :

and for the filter procedure :

T. . = u!u'. (2.33b)

ID I D

Smagorinsky [1963] suggested to relate the SGS stress r. . to the resolvable velocity field by an eddy viscosity model. This suggestion has ever since been followed with considerable success in LES calculations. For the filter procedure the proposed relation reads :

T. . = - 2 v S. . + ~ 6. . u'u' (2.34a)

13 T i] 3 13 k k

When applied to the volume balance procedure we have :

K. . = - 2 ^ S. . + \ 6. . V V (2.34b)

ID T I D 3 i : k k

The s t r a i n tensor S. . i s calculated from :

iD

S. . = \ ( 6 Xü. + 6. jü . ) ( 2 . 3 5 a )

I D 2 j l 1 3

which is the finite difference analogue for S. ., defined by : 3u. 3u .

S

i3 = i

(

i ^

+

^ >

(2

"

35b)

D i

The second term on the right-hand side of eqn (2.34a,b) is generally not explicitly calculated, but added to the pressure term to result in a virtual

— * pressure (indicated by p ) :

(36)

= p + t u!u' (2.36a)

3 ï ï and

p = p + - u: u! (2.36b) 3 x i

Smagorinsky [1963] suggested the following model for v :

»_ = ( C A )2 (2 S. . S. , )1 / 2 (2.37)

T S A ij i]

where C i s the s o - c a l l e d Smagorinsky c o n s t a n t (and A as for eqn ( 2 . 8 ) ) . This "Smagorinsky model" can be d e r i v e d by e q u a t i n g the s u b g r i d - s c a l e p r o d u c t i o n and d i s s i p a t i o n in homogeneous t u r b u l e n c e . Kwak e t a l . [1975] used t h e v o r t i c i t y u. i n s t e a d of t h e s t r a i n t e n s o r S . . and c a l l e d the

i ID

r e s u l t i n g r e l a t i o n the " v o r t i c i t y model" :

^T = (CvAfi)2 (ü± Ü J1 7 2 ( 2 . 3 8 )

Compared to the Smagorinsky model of eqn (2.37) it has the advantage of being zero in irrotational regions. However, no significant differences have been found between applications of the two models (Kwak et al. [1975], Ferziger [1977], McMillan et al. [1979]).

Schumann [1973] relates v to the subgrid-scale kinetic energy Ë'

rather than to S. .. In essence his model reads :

\ = C

2

V

/ 2 j

Ê'

1 / 2

(2.39)

A is the size of the surface over which the averaging takes place and its square root is consequently a characteristic length scale. The characteristic velocity is the square root of the SGS kinetic energy

j

E' = ^

(u

±

-

j

u\)

2

=

\ [ u

2

- ( V )

2

] (2.40)

Ë' is evaluated from a separate and quite complicated transport equation for Ê' (Schumann [1973,1975a]). Fortunately Schumann [1975a] found that application of the Smagorinsky model of eqn (2.37) did not cause

(37)

considerable changes to his results (choosing A =(AxAyAz) ). The model used in this work will therefore be based on the relatively simple Smagorinsky model :

%= ( CSV/ 2)2 (2 S..S..)1/2 (2.41)

The value used for C is discussed in section 6.5.

The subgrid-scale models described in this section were derived for turbulence with local isotropic behaviour. In regions where the mean (resolvable) flow field is inhomogeneous, such as near solid walls, the behaviour is anisotropic. Special care is required for these regions. Chapter 4 is devoted to problems that arise near solid walls with sec. 4.3 in particular to SGS modelling.

2.8 - Relation between volume balance and filter procedure

The surface averaged variables of the volume averaged NS equations (2.22) can be approximated by volume averaged ones with the aid of eqns (A.3-5) derived in appendix A : V._ V ._ 3 xu. V ._ V. V. 3 "V V._ —T T * = - 6. [ V Dü. + V u ' . ] + :: 6. ~Z—i - 6. 1p + HOT (2.42) 3t ] i ] l 3 R 3 OX . l 2

where HOT represents truncation errors of order Ax and smaller. The partial derivatives in the filtered NS equations (2.25) can be discretized with central difference quotients, resulting in :

3u. _ _ 3u.

— - = - 6 . [Ü.Ü. + U!u'. + L. . + CR. . ] + ~ 6 . " T ^ - S.p + HOT ( 2 . 4 3 )

3 t D I D I D ID ID R D °x • i

There is a close resemblance between eqns (2.42) and (2.43). Both equations contain variables that are identically connected to similar grid elements (provided that a staggered grid is applied for the filtered equations). Apart from the subgrid-scale terms, the only difference is the definition of the variables (volume averaged versus filtered). Consequently, the resulting computer program codes for the resolvable parts of both equations will be identical. In fact the filter inherent in the volume

(38)

balance procedure is the box or top-hat one given by eqn (2.8) with A =h . Differences occur for the subgrid-scale terms, since L. . and CR.. in eqn

ID 13

(2.43) are not present in eqn (2.42). Substitution of the approximation formulas for L. . and CR. ., given by eqns (2.31) and (2.32), into the centered difference terms of eqn (2.43) results in :

6. [L. . + CR. .] = 6. [ -7 A^ V2(ü.ü.) - rr A2 (ü.V2ü. + Ü.V2ü. ) ] (2.44)

3 13 13 3 24 A i n 24 A 1 3 3 x

Working out the central difference for a staggered grid element V . results

2 ^ in a term of order Ax.. This means that it is not necessary to calculate the

L. . and CR. . terms explicitly when applying the volume balance procedure, since they turn out to be of the order of magnitude of the truncation errors gathered in HOT. It can therefore be concluded that when second order accurate schemes are applied in filter procedure calculations, the volume balance procedure and the filter procedure are in fact identical. This means that the Smagorinsky model, which proved to be a good model in conjunction with the filter procedure, can equally well be used with the volume balance procedure. Schumann [1975a] came to a similar conclusion after performing numerical experiments with the volume balance procedure where the Smagorinsky model had been adopted.

(39)

3 - THE BASIC EQUATIONS IN DISCRETIZED FORM

3.1 - Introduction

The terms that occur in the volume averaged NS equations (2.22) need approximations both in space and time. Obviously, spatial approximations are required for the partial derivatives that remained in the diffusive terms. However, spatial approximations are also necessary to express the volume averaged velocities in surface averaged ones, and vice versa for pressure. Moreover, surface averaged velocities that are not located at positions defined in the staggered grid, have to be approximated in terms of those that do. Two different schemes will be studied to approximate the resolvable part of the convective terms: the central and the upwind difference scheme. With respect to the approximations in time three different three-time-level methods will be applied : the Leapfrog method, the DuFort Frankel Leapfrog method and the Adams Bashforth method. An Euler method is used for the diffusive and the upwind difference approximated resolvable scale convective terms. The maximum allowed time steps required for numerical stability will be estimated. The orders of magnitude of the truncation errors are mutually compared and related to the terms that model the subgrid-scale motions.

For brevity we confine ourselves to the first of the volume averaged NS equations (2.22) and consider only two-dimensional examples. It will not be difficult to generalize the findings to all components of the equations in a three-dimensional space domain.

The instationary-, resolvable scale convection-, diffusion- and (virtual) * pressure terms of eqn (2.22) are replaced by symbols I., R. ., D. . and P..

l i] 13 1

Leaving out the subgrid-scale stress terms this means that eqn (2.22) can be written as :

3

I. = E [ - R. . + D. . ] - P* (i=l,2,3) (3.1)

i j = 1 1 3 1 3

Figure 3.1 shows a cross-section through a volume V . The surface-averaged velocities and volume surface-averaged pressures in and around this volume are used in this chapter to discretize the terms of the above equation. The indices used for the velocities and pressures in this figure indicate their positions in the grid.

(40)

-*u_ o

"

U

w

J>\

W

- ► n o -►

-—x

figure 3.1 - velocities and pressures in the vicinity of volume V (dotted)

3.2 - Spatial approximations

3.2.1 - Instationarv term

The instationary term I of the volume averaged first NS equation reads

V

i = * —u

x at (3.2)

It contains a velocity component u that has been averaged over volume V .

V _ x

This velocity u has to be expressed in terms of a surface averaged velocity u, since velocities u are only in this way defined in the staggered grid. Using eqn (A.3), derived in appendix A, the desired relation can easily be found :

x-T 3 " P Ax-T

*x = I T

+ A l

x

(3.3a)

where the truncation error Al equals

AI

X = ii

A

*

2 ' r r ï ' p

+ H0T

3t 3x

(41)

3.2.2 - Resolvable scale convective terms central difference approximation

The convective terms R and R that are considered in the 2d example of xx xy

the first volume averaged NS equation read :

R = 6 (Xü Xü ) = T1 (Xü Xü - Xü Xü ) (3.4a)

xx x Ax e e w w

R = 6 (yv yü ) = T1 (yv_yü\ - yv yü\ ) (3.4b)

xy y Ay t t b b

The velocities Xu and yv in eqn (3.4) are not located in the staggered grid

at planes that actually contain these components (see fig. 3.1). Furthermore the velocities u in eqn (3.4b) are averaged over surfaces with normal in y-instead of x-direction. The latter problem can be overcome with eqn (A.6) (derived in appendix h) that relates averaging over perpendicular planes

contained in volume V. Using this equation and Taylor series the following central difference approximations for R and R can easily be derived :

rr xx xy Rxx - ± < \ < W - % < W * + A Rx x ( 3-5 a ) R

xy = ife < ^t

(

V % > -

"h < W 5 + AR

xy

(3.5b)

with : X- 1 .X- X- . X- 1 , X- X- . U

i

=

2

( U

E

+

V

U

w

=

2

( U

W

+

V

yv~ = i (yv + yv ) Yv , = i (yv + Yv ) t 2 k E P; b 2 v EB B'

The truncation errors AR and AR read : xx xy

AR

*xx = ' J *** i J < " ~T>*

P + H0T (3

-

5c)

3x 2 2 X -»„ « 2 j _ i 3 , y - 3 u , 1 3 , y - 3. v , ^ AR

xy =

A x {

Ü By"

( V

^ i >

P

' i a?

( U 3 x

2

}

P

}

-

Ay2

6 i ^

a

rl

ü

)p

+ H0T ( 3

-

5 d

>

3y

(42)

upwind difference approximation

The upwind difference scheme (Roache [1976], Patankar [1980]) uses one sided rather than space-centered differencing. Backward differences are used when the velocities are positive and vice versa; for R the approximation reads:

1 ( , iX- i. , i i. x-Rxx = ito { " < V l Uwl' UW ' ( V l "wl' UP

+ (

X

=

g

-|

X

=

5

D %

+

(

X

V|

X

u

?

D % } + AR

xx

(3.6a)

with approximated values for the velocity components at w and e (figure 3.1) as :

x- 1 ,x- x- . x- 1 ,x- , x- .

Uw = 2 ( UW + V Ui = 2 ( UP + UE>

The upwind approximation for R is :

B = _L_ { _ (yv~+|yv~|) xü - (yv~-|yv~|) Xü xy 2Ay x ( E ' B'' UB ( b ' b>' UP + (yvr|yv£|) % + (yvE+|yvE|) Xüp } + A Rx y (3.6b) with : yv ~ = ~ (yv + Yv ) yv ~ = - (yv + Yv ) b 2 * B EB' t 2 V P E;

The truncation errors AR and AR can be calculated in the same way as for xx xy

the central difference approximation, using eqn (A.6) and Taylor series expansions. Considering the case wnere all velocities are positive, we find:

1 x - A - d2 Xü , 1 . , 3Xu \ 2 1 , 2 3 , x - 32 Xü , A Rxx = 2 UP ^ < — T > P + 2 A x 'ix^P " Ï A x to ( U TT*? + H 0 T 3x 3x ( 3 . 6 c ) AO I y~ A , 92 yü> 1 . ,3yv, ,3yu> 1 * 2 3 ,y- 32 yü , A R

x y = 2

V

P

Ay ( 3y2

>

P +

i

Ay (

iF

j

P

(

^'

,

P " 6

Ay

^ <

v

- ^ T ' p

" 8

^

f, <» ^*

P

+ H0T

<

3

^>

3x

(43)

3.2.3 - Diffusive terms

The diffusive terms D and D that are considered in the 2d example of

xx xy r

the first volume averaged NS equation read :

x— x—

D

xx = i Z x t <*£>» -<**!>. J

(3

-

7a)

The partial derivatives in the above equations can be approximated using

Taylor series expansions and eqn (A.6) derived in appendix A, resulting in :

D = Z - ^ [Xu_ - 2XÜD + Xü „ ] + AD (3.8a)

XX R . 2 E P WJ XX

Ax

D

xy = i ^2" t % "

2

%

+

" B

] + A D

xy

(3

"

8b

'

with for the errors :

AD

XX

= - i i ¥ <

a

7i

u

>p

+ HOT (3

-

8c)

3x

*

D

xy = - Ii ¥ ^ ' P * it ^ 4

y

^ > P ♦ HOT ,3.Bd,

3y 3y 8x

If the partial derivatives of eqn (3.7a,b) have to be calculated for planes coinciding with solid walls, the approximations derived in the foregoing can no longer be used. For these situations a special treatment is required that will be discussed in section 4.2. Other types of boundary conditions, required to approximate velocities on or near edges of the flow domain, are discussed in section .5.4.

(44)

3.2.4 - Pressure term

The pressure term P :

* 1 ,x-« x-*. ,, _v

Px = to ( PE " PW> ( 3-9 )

contains surface averaged pressures. They are replaced by volume averaged ones in order to obtain a scalar type variable for pressure that is the same for all three volume averaged NS equations. Eqn (A.5) derived in appendix A is used for the aforementioned replacements, resulting in :

K - Kx <% -

v

p ; >

+

< <

3

— >

with for the error term AP x

a3 x-*

AP* = - zr Ax^ (d P )p + HOT (3.10b)

X 24 , , 3 P 3x

3.3 - Time advancement

For a start the discussion will be confined to the spatially approximated 2d instationary convection diffusion equation as a part of the first volume averaged NS equation. The resolvable part of this equation reads :

-,x-3 u

3t - R XX - R + D + D _{xy xx xy} (3.11)

where the abbreviations used on the right-hand side are defined by eqns (3.5.) or (3.6) and (3.8). In order to approximate the time derivative, integration with respect to time t has to be carried out over a certain interval At. In the following sub-sections methods to perform this integration will be discussed. These methods can only result in numerically stable calculation procedures if the interval At is smaller than a certain critical time step. In section 3.4 formulas are presented by which these critical time steps can be estimated for the various methods described in the present section. Their validity is checked by numerical experiments discussed in chapter

(45)

6-3.3.1 - Euler/Leapfrog method

A method applied in LES calculations by e.g. Schumann [1973] and Grotzbach [1977] is the Leapfrog method with respect to the resolvable scale convective terms and the Euler method for diffusion. This three-time-level Euler/Leapfrog (E/L) method calculates the unknown velocity u (upper right indices denote the time step in the calculation) from :

x-n+1 _ x-n-1

"P ~ "P = - Rn - Rn + D""1 + D""1 + AT (3.12a)

2 At xx xy xx xy E/L

The spatial discretization for the R and the D terms are obtained by central difference approximations (eqns (3.5) and (3.8)). The error AT , made by this approach can be calculated using Taylor series expansions :

a3 xu 3D 3D

A TE/L = " 6 A t < —- f) U - 2 A t t(^ f »n + '"if ^ + H 0 T ( 3-1 2 t ) at

3.3.2 - DuFort Frankel Leapfrog method

The DuFort Frankel Leapfrog (DFFL) method (Roache [1976]) uses Leapfrog with respect to the resolvable scale convective terms. The diffusive terms D and D are now calculated according to DuFort and Frankel's [1953] _{xx xy} method : _n+l 1 ,x-n , x-n x-n+1 x-n-1. ,, ,, . °XX = . 2 ( UE + UW " UP " UP > ( 3-1 3 a ) R Ax _n+l 1 ,x-n , x-n x-n+1 x-n-1. , , .,,,, Dxy = D . 2 ( UT + UB " UP - UP ) ( 3-1 3 b ) R Ay

The error AT made by the approximated convection diffusion equation : Dr FL

x-n+1 _ x-n-1

"P ~ "P = - R" - R" + Dn + 1 + D "+ 1 + AT pT (3.14a)

2 At xx xy xx xy DFFL can be estimated in a similar way as for the Leapfrog scheme :

(46)

AT

DpFL

= - J (£)

2

( ^ V - J At

2

< ^ , »

+

HOT (3.14b)

ot dt

3.3.3 - Euler/Adams Bashforth method

The Euler/Adams Bashforth (E/AB) method (Roache [1976]) is also a three-time-level method. It calculates the unknown velocity u from :

x-n+1 _ x-n

"P " "P = -f(Rn + Rn ) + ^ R " "1 + R " -1) + Dn + D" + A TW a R (3.15a)

At 2 xx xy 2 xx xy xx xy E/AB using Adams Bashforth for convection and Euler for diffusion. The R and D terms are again given by eqns (3.5) and (3.8). The error AT . made when applying this scheme is found to be :

A TE / A B = ^ [ ( ^ f ) " + < ^ f ) " ] - j A t [ ( ^ , « + < ^ f , » ] + HOT (3.15b)

ot dt

3.3.4 - Euler/upwind method

The Euler/upwind difference approximation (E/u) for the resolvable convective terms (sec. 3.2.2) allows the use of a two-time-level approximation for the integration in time. This is possible without loss of numerical stability, as will be shown in the next section. The advantage over the methods discussed in the previous sub-sections is that only two velocity vector-fields instead of three have to be stored at the same time. With R and R defined by eqn (3.6) the unknown velocity u is calculated from :

x-n+1 x-n

U P " "P = - R " - Rn + Dn + Dn + AT , (3.16)

At xx xy xx xy E/u

The error AT . can be calculated using Taylor series expansions as : ,2

Numerical and physical aspects of large eddy simulation of turbulent flows