values at a location sufficiently far downstream of the inlet (Re
θ≈ 1300).
The sampled values were then rescaled with eddylicious using the rescaling procedure defined by the LWS method [26] to approximately match the Re
θat the inflow of the first TBL simulation. The obtained values were then used as inflow for the second TBL simulation, see figure 6.6. This set-up mimics the weak recycling method defined in [26], however, instead of obtaining the inflow values dynamically, they are computed using a database from a precursor TBL simulation. Apart from the inflow generation method the set-up of the second TBL simulation is identical to that of the first one.
Figure 6.7 shows the evolution of the friction coefficient, c
f, and the shape factor, H = δ
∗/θ, as a function of Re
θ. For c
fit is clear that the adaption length for the proposed method is very short, while the rescaling-based inflow generation leads to a dip in c
fbefore recovery happens. For the shape factor, the overall agreement with DNS is worse for both methods.
The under-prediction of H can, however, be anticipated based on the results for channel flow shown in the previous section. Indeed, for the chosen grid resolution the produced mean velocity profiles were flatter than the DNS data, which corresponds to an under-predicted shape factor. Regarding the adaption length, it is clear that for the proposed method the behaviour of H stabilises faster. The curve corresponding to the LWS rescaling exhibits a steep decline over a significant part of the domain, which is not consistent with the behaviour of the DNS data.
Figure 6.7: Evolution of the friction coefficient, c
f, and the shape factor, H, as a function of Re
θ.
Figure 6.8 shows the mean streamwise velocity at three Re
θfor which
DNS data are provided in [36]. The overall accuracy of the profiles is in line
with what could be expected based on the grid resolution study for channel
flow. An important result is that at Re
θ≈ 1420 the profiles corresponding
to the two inflow generation methods are already in excellent agreement
with each other. Noticeably, in the inner layer the agreement is fairly good
6.2. Turbulent Boundary Layer 47 even at Re
θ≈ 1000, which reflects the fact that channel flow and the TBL have a similar flow structure in that region.
Figure 6.8: Profiles of the streamwise component of velocity at Re
θ≈ 1000, Re
θ≈ 1420, Re
θ≈ 2000. The profiles are shifted by hui/u
+τ= 5 upwards for increasing Re
θ.
The overall conclusion is that the correct flow behaviour is recovered at Re
θsomewhere in the range 1200-1300. The corresponding range in terms of the momentum thickness of the inflow TBL is 270θ
in-370θ
in. For the rescaling method, the adaption length is of comparable size. Based on these results it is possible to conclude that the performance of the proposed inflow generation method is on par with other precursor-based approaches.
Combined with its simplicity and robustness, this makes it an attractive
alternative to methods based on weak recycling.
Chapter 7
Summary and Future Work
The work presented in this thesis can be seen as consisting of three parts.
The first part is the development of the methodology for computing LES of wall-bounded flows using a general-purpose CFD solver. To that end, turbulent channel flow was chosen as a model problem. Various aspects of the simulation set-up and modelling choices where analysed such as the con-struction of the computational grid and the SGS modelling. Several further developments are possible here. One is to consider the possibility of using unstructured grids. While this is not needed for most academic applications where the geometry of the flow is commonly simple, unstructured grids are the standard in industrial applications. A different possibility is to consider nested (embedded) structured grids. The fact that OpenFOAM can use arbit-rary polyhedral cells makes this idea especially attractive. The benefit of using a nested grid would be the ability to make the stream- and spanwise resolutions of the grid a function of y. This would lead to large savings in computational resources, especially for higher Re-number flows.
The second major contribution of this work is the software package eddylicious. Multiple ways of expanding the package are possible. This includes adding support for other CFD solvers, implementing more inflow generation methods, and supporting arbitrary inlet geometry. Perhaps most importantly, further promotion of the package among the CFD community has to be done, since the implementation of the above mentioned function-ality is only possible as a collaborative effort.
The third part of this thesis is the work on the proposed inflow genera-tion method. A framework for using channel flow as a precursor simulagenera-tion for generating inflow for a TBL simulation was presented. While the meth-odology is established, it is necessary to see how the method performs for a wider range of flows. Of particular interest for our research group are flows with a separating boundary layer. As mentioned in chapter 4, numerous
49
issues with strong recycling methods such as the one proposed here remain
unsolved. In particular, no concrete guidelines exist for how the grids of the
precursor and the main simulations should be related besides for the general
conclusion that it is best that they coincide. Also of interest is to test the
range of applicability of the LWS method and its variations. For instance,
the difference in the Re-number between the rescaled flow and the inflow
that can be allowed is still not known [19].
Bibliography
[1] Y. Bentaleb, S. Lardeau, and M. Leschziner, Large-eddy simulation of turbulent boundary layer separation from a rounded step, Journal of Turbulence 13 (2012), no. 4.
[2] D. R. Chapman, Computational aerodynamics development and outlook, AIAA Journal 17 (1979), no. 12, 1293–1313.
[3] H. Choi and P. Moin, Grid-point requirements for large eddy simulation:
Chapman’s estimates revisited, Physics of Fluids 24 (2012), no. 1, 30–
35.
[4] Y. M. Chung and H. J. Sung, Comparative study of inflow conditions for spatially evolving simulation, AIAA Journal 35 (1997), no. 2, 269–274.
[5] E. De Villiers, The potential of large eddy simulation for the model-ing of wall bounded flows, Ph.D. thesis, Imperial College of Science, Technology and Medicine, 2006.
[6] D. Dietzel, D. Messig, F. Piscaglia, A. Montorfano, G. Olenik, O. T.
Stein, A. Kronenburg, A. Onorati, and C. Hasse, Evaluation of scale resolving turbulence generation methods for large eddy simulation of turbulent flows, Computers & Fluids 93 (2014), 116–128.
[7] E. R. Van Driest, On turbulent flow near a wall, Journal of the Aero-nautical Sciences 23 (1956), no. 11, 1007–1011.
[8] J. H. Ferziger and M. Peric, Computational methods for fluid dynamics, 2002.
[9] C. Fureby and F. F. Grinstein, Large eddy simulation of high-Reynolds-number free and wall-bounded flows, Journal of Computational Physics 181 (2002), no. 1, 68–97.
[10] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, A dynamic subgrid-scale eddy viscosity model, Physics of Fluids A: Fluid Dynamics 3 (1991), no. 7, 1760.
51
[11] F. F. Grinstein, L. G. Margolin, and W. J. Rider, Implicit large eddy simulation: computing turbulent fluid dynamics, Cambridge University Press, 2007.
[12] A. H. Herbst, P. Schlatter, and D. S. Henningson, Simulations of tur-bulent flow in a plane asymmetric diffuser, Flow, Turbulence and Com-bustion 79 (2007), no. 3, 275–306.
[13] K. Horiuti, Large eddy simulation of turbulent channel flow by one-equation modeling, Journal of the Physical Society of Japan 54 (1985), no. 8, 2855–2865.
[14] R. I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, Journal of Computational Physics 62 (1986), no. 1, 40–65.
[15] N. Jarrin, Synthetic inflow boundary conditions for the numerical simu-lation of turbulence, Ph.D. thesis, The University of Manchester, 2008.
[16] H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows, Ph.D. thesis, Imperial College of Sci-ence, Technology and Medicine, 1996, p. 394.
[17] J. W. Jewkes, Y. M. Chung, and P. W. Carpenter, Modifications to a turbulent inflow generation method for boundary-layer flows, AIAA Journal 49 (2011), no. 1, 247–250.
[18] H. Kanchi, K. Sengupta, and F. Mashayek, Effect of turbulent inflow boundary condition in LES of flow over a backward-facing step using spectral element method, International Journal of Heat and Mass Trans-fer 62 (2013), no. 1, 782–793.
[19] A. Keating, U. Piomelli, E. Balaras, and H.-J. Kaltenbach, A priori and a posteriori tests of inflow conditions for large-eddy simulation, Physics of Fluids 16 (2004), no. 12, 4696.
[20] W.-W. Kim and S. Menon, An unsteady incompressible Navier–Stokes solver for large eddy simulation of turbulent flows, International Journal for Numerical Methods in Fluids 31 (1999), no. 6, 983–1017.
[21] J. Larsson, Blending technique for compressible inflow turbulence: Al-gorithm localization and accuracy assessment, Journal of Computa-tional Physics 228 (2009), no. 4, 933–937.
[22] H. Le, P. Moin, and J. Kim, Direct numerical simulation of
turbu-lent flow over a backward-facing step, Journal of Fluid Mechanics 330
(1997), 349–374.
Bibliography 53 [23] M. Lee and R. D. Moser, Direct numerical simulation of turbulent
chan-nel flow up to Re
τ≈ 5200, Journal of Fluid Mechanics 774 (2015), 395–415.
[24] D. K. Lilly, A proposed modification of the germano-subgrid-scale clos-ure method, Physics of Fluids A 4 (1992), no. 3, 633–635.
[25] K. Liu and R. H. Pletcher, Inflow conditions for the large eddy sim-ulation of turbulent boundary layers: A dynamic recycling procedure, Journal of Computational Physics 219 (2006), no. 1, 1–6.
[26] T. S. Lund, X. Wu, and K. D. Squires, On the generation of turbulent inflow conditions for boundary layer simulations, Journal of Computa-tional Physics 140 (1998), 233–258.
[27] C. Meneveau, T. S. Lund, and W. H. Cabot, A Lagrangian dynamic subgrid-scale model of turbulence, Journal of Fluid Mechanics 319 (1996), 353–385.
[28] N. Nikitin, Spatial periodicity of spatially evolving turbulent flow caused by inflow boundary condition, Physics of Fluids 19 (2007), no. 9, 1–4.
[29] T. Nishikawa, Y. Yamade, M. Sakuma, and C. Kato, Application of fully-resolved large eddy simulation to KVLCC2, Journal of the Japan Society of Naval Architects and Ocean Engineers 16 (2012), 1–9.
[30] J. M. ¨ Osterlund, Experimental studies of zero pressure-gradient turbu-lent boundary layer flow, Ph.D. thesis, Royal Institute of Technology, 1999.
[31] S. B. Pope, Turbulent flows, Cambridge University Press, 2000.
[32] F. T. Pronk, A comparison of inflow generation methods for large-eddy simulation, Master thesis, Delft University of Technology, 2012.
[33] S. Rezaeiravesh, M. Liefvendahl, and C. Fureby, On grid resolution requirements for LES of wall-bounded flows, Proceedings of ECCOMAS Congress 2016, 2016.
[34] H. Rusche, Computational fluid dynamics of dispersed two-phase flows at high phase fractions, Ph.D. thesis, Imperial College of Science, Tech-nology and Medicine, 2002.
[35] P. Sagaut, Large eddy simulation for incompressible flows,
Springer-Verlag, 2006.
[36] P. Schlatter and R. ¨ Orl¨ u, Assessment of direct numerical simulation data of turbulent boundary layers, Journal of Fluid Mechanics 659 (2010), 116–126.
[37] U. Schumann, Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli, Journal of Computational Physics 18 (1975), no. 4, 376–404.
[38] J. Sillero, J. Jim´enez, R. D. Moser, and N. P. Malaya, Direct simulation of a zero-pressure-gradient turbulent boundary layer up to Re
θ= 6650, Journal of Physics: Conference Series 318 (2011), no. 2.
[39] M. P. Simens, J. Jim´enez, S. Hoyas, and Y. Mizuno, A high-resolution code for turbulent boundary layers, Journal of Computational Physics 228 (2009), no. 11, 4218–4231.
[40] J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, and D. Mavriplis, CFD vision 2030 study: A path to revolu-tionary computational aerosciences, Tech. report, NASA, 2014.
[41] J. Smagorinsky, General circulation experiments with the primitive equations, Monthly Weather Review 91 (1963), no. 3, 99–164.
[42] P. R. Spalart, Direct simulation of a turbulent boundary layer up to R
θ= 1410, Journal of Fluid Mechanics 187 (1988), 61–98.
[43] P. R. Spalart, W. H. Jou, M. Strelets, and S. R. Allmaras, Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach, Advances in DNS/LES, vol. 1, 1997, pp. 4–8.
[44] P. R. Spalart, M. Strelets, and A. Travin, Direct numerical simulation of large-eddy-break-up devices in a boundary layer, International Journal of Heat and Fluid Flow 27 (2006), no. 5, 902–910.
[45] A. Spille-Kohoff and H.-J. Kaltenbach, Generation of turbulent inflow data with a prescribed shear-stress profile, DNS/LES - Progress and Challenges, Proceedings of the third AFOSR International Conference on DNS/LES, 2001, pp. 319–326.
[46] B. Stevens, C.-H. Moeng, and P. P Sullivan, Large-eddy simulations of radiatively driven convection: Sensitivities to the representation of small scales, Journal of the Atmospheric Sciences 56 (1999), no. 23, 3963–3984.
[47] G.R. Tabor and M. H. Baba-Ahmadi, Inlet conditions for large eddy
simulation: A review, Computers & Fluids 39 (2010), no. 4, 553–567.
Bibliography 55 [48] X. Wu, Inflow turbulence generation methods, Annual Review of Fluid
Mechanics 49 (2017), no. 1.
[49] X. Wu, K. D. Squires, and T. S. Lund, Large eddy simulation of a spatially-developing boundary layer, Supercomputing, 1995. Proceed-ings of the IEEE/ACM SC95 Conference, 1995.
[50] A. Yoshizawa, A statistically-derived subgrid model for the large-eddy simulation of turbulence, Physics of Fluids (1958-1988) 25 (1982), no. 9, 1532–1538.
[51] A. Yoshizawa and K. Horiuti, A statistically-derived subgrid-scale
kin-etic energy model for the large-eddy simulation of turbulent flows,
Journal of the Physical Society of Japan 54 (1985), no. 8, 2834–2839.
Paper I
Large-Eddy Simulation of Turbulent Channel Flow
Timofey Mukha, Mattias Liefvendahl
Department of Information Technology Uppsala University
Box 337, SE-751 05 Uppsala, Sweden
Technical report 2015-014 May 2015 ISSN 1404-3203
1 Introduction 3
2 Turbulence modelling and numerical methods 5
2.1 Governing equations . . . 5 2.2 Large-eddy simulation . . . 5 2.3 Subgrid stress modelling . . . 7 2.4 Discretization and interpolation methods . . . 8 2.5 Solver algorithm . . . 9
3 Simulation case set-up 11
3.1 The continuous problem and the physical parameters . . . 11 3.2 Mesh generation . . . 13 3.3 Boundary conditions . . . 14 3.4 Modelling the pressure gradient . . . 14 3.5 The simulation campaign . . . 15
4 Results 16
4.1 Global flow quantities . . . 16 4.2 Mean velocity profiles . . . 17 4.3 Velocity fluctuations . . . 18 4.4 Turbulent shear stress . . . 20 4.5 Skewness and flatness . . . 22 4.6 Vorticity fluctuations . . . 25 4.7 Two-point velocity correlations . . . 26
5 Conclusions 30
A Definitions and methods of computation for statistical quantities 34 A.1 Averaging . . . 34 A.1.1 Temporal averaging . . . 34 A.1.2 Spatial averaging . . . 35 A.2 Higher order statistical moments . . . 35 A.3 Two-point correlations . . . 36
2
Chapter 1
Introduction
Results and analysis are reported from a simulation campaign, using large-eddy simulation (LES), of fully developed turbulent channel flow. Channel flow is a classical model problem for the investigation of wall-bounded turbulence, and it has been extensively studied using both ex-perimental and computational methods, [11, 7]. The flow takes place between two infinite parallel planes, and it is driven by a constant pressure gradient. Due to the particularly simple geometry, specially designed spectral methods can be applied to the problem, [9, 7, 10, 3, 5], which provides extremely good accuracy and allows for direct numerical simulation (DNS), i.e. without turbulence modelling, up to moderate Reynolds numbers (at least Reτ> 2000).
The purpose of the present study is to evaluate the predictive capability of the numerical methods and turbulence modelling (a particular LES-subgrid model is tested) implemented in the general-purpose CFD-software OpenFOAM1. The subgrid model is based on the ideas first proposed in [14], and relies on a subgrid viscosity which is computed using an additional transport equation for the subgrid kinetic energy, see section 2.3 for a detailed description. One Reynolds number is investi-gated, Reb= 13 350, which roughly corresponds to Reτ= 395, see the discussion in section 3.1, for the parameters of the problem and the definition of these numbers. Corresponding simulations are carried out on three grid refinement levels in order to investigate grid convergence, and how the subgrid model behaves with varying grid size. The coarsest grid contains 135 000 cells and the finest grid 8.65·106cells. The results obtained are naturally not as accurate as those obtained by spectral methods, at the same computational cost. But, since the methods employed here can be applied to general configurations/geometries, it is essential to have a good understanding of their predictive capability for wall-bounded turbulence, and channel flow is an ideal case for this investigation.
A wide range of turbulence statistics are computed, and are included in the mesh refinement study. Due to the geometry/symmetry of the problem, all statistics only depend on the wall-normal coordinate. The resulting profiles are presented for the following statistical quantities. The mean (streamwise) velocity and all components of the Reynolds stress tensor, hence including the level of turbulent fluctuations of the velocity components. Higher order statistical moments, skewness and flatness, are also computed for all three velocity components. Profiles of the vorticity fluctuations are included as well. For all these quantities, the statistics are computed using time series at a given spatial location. In addition to this, the two-point correlations of velocity fluctuations have been
1Disclaimer: this study is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM Rand OpenCFD Rtrade marks.
3
statistics and how this depends on the mesh resolution.
This study should be useful for CFD-users of OpenFOAM, and similar software, and provides detailed information concerning grid requirements, and other settings, for problems where turbulent boundary layers play an important role. All of the data presented in the report is made available online2along with complete OpenFOAM simulation case set-up files.
2https://bitbucket.org/lesituu/channel flow data
4
Chapter 2
Turbulence modelling and numerical methods
The mathematical model and numerical methods employed in the study are described in this chapter. This includes a description of the underlying governing equations, the approach used towards turbulence modelling, the numerical methods employed for discretizing the equations, and the algorithm employed by the solver program to take the pressure-velocity coupling into account.
2.1 Governing equations
The momentum equation for an incompressible viscous fluid is the (incompressible) Navier-Stokes equation which, in the absence of external forces, is given by,
∂u
∂t+∇ · (u ⊗ u) = −1
ρ∇˜p + ν∆u. (2.1)
Here u is the velocity of the fluid, ρ is the density, ˜p is the pressure, and ν is the kinematic viscosity.
The incompressibility of the fluid is expressed by,
∇ · u = 0. (2.2)
Equations (2.1) and (2.2) constitute the governing equations which fully describe the fluid flow.
Note that ρ and ν are assumed constant. Below, p = ˜p/ρ will be referred to as the “pressure” for brevity. Both the notation, (u, v, w), and the notation, (u1, u2, u3) will be used for the Cartesian components of the velocity vector.
2.2 Large-eddy simulation
The need for a different approach in turbulent fluid flow simulation rather then attempting to solve (2.1)-(2.2) directly stems from the fact that resolving all the relevant scales of turbulent motion requires both spatial and temporal resolutions that are incompatible with the computational resources currently available.
5
computational meshes possible as well as increasing the time-step. Formally, the separation of scales is achieved via a filtering operation which is mathematically expressed as a convolution of the relevant flow field with a chosen filter kernel,
φ(x, t) = ZZZ∞
−∞
φ(x, t)G (x− ξ, ∆) d3ξ, (2.3)
where G is the kernel and ∆ is the filter’s cut-off width, which is a parameter that defines the size of the scales that are filtered out. Throughout the report we will use the over-bar to denote filtered variables. The unresolved part of the flow field, which is left out after the filtering is defined as,
φ00(x, t) = φ(x, t)− φ(x, t). (2.4)
The length scales associated with the unresolved part are commonly referred to as subgrid scales (SGS).
A large variety of filters with different properties is described in the literature. For an in-depth discussion of the most important such filters, we refer to, [11, 13]. In spite of the large number of filters available, most of them are difficult to apply in a general-purpose CFD-code. The filter most commonly used in conjunction with finite volume discretization is the top-hat filter, see page 99 of [16],
G (x− ξ, ∆) =
(1/∆3, |x − ξ| ≤ ∆/2
0, otherwise. (2.5)
It is evident from (2.5) that the filtering gives a value which is an average over a rectangular volume
∆3. A common choice for ∆ is the cubic root of the volume of the (local) computational cell,
∆ =p3
∆x∆y∆z, (2.6)
where ∆x, ∆y and ∆z are the cell-sizes along the corresponding coordinate axes. This choice of ∆ conveniently makes φ equal to the average value of φ in the computational cell. This implies that no explicit filtering needs to be performed during the computational procedure, instead the filtering is built into the discretization method itself. Due to this attractive property the combination of the top-hat filter (2.5) and the ∆ defined as (2.6) is used in the channel flow simulations described in this report.
By formally applying the filtering operation to the continuity equation (2.2) and Navier-Stokes equations (2.1) it is possible to derive conservation laws for the filtered flow variables. When doing this, it is convenient to use indicial tensor notation, which leads to more compact expressions.
Due to the linearity of the continuity equation, applying the filtering is straightforward. The form of the equation remains unchanged,
∂ui
∂xi
= 0. (2.7)
This result implies that the SGS velocity field u00is also solenoidal. Filtering the Navier-Stokes equations results in the following.
∂ui
∂t + ∂
∂xj
(uiuj) =−∂p
∂xi
+ ν ∂2ui
∂xj∂xj
(2.8)
6
The main complication here is that the advection term in (2.8) cannot be expressed in terms of ui. The common way to address this issue is to introduce the so called SGS stress tensor B, the components of which are defined by,
Bij= uiuj− uiuj. (2.9)
Inserting this into equation (2.8) leads to,
∂ui
∂t + ∂
∂xj
(uiuj) =−∂p
∂xi−∂Bij
∂xj
+ ν ∂2ui
∂xj∂xj
. (2.10)
In order to close this system of equations, B has to be modelled.
2.3 Subgrid stress modelling
A rich variety of approaches to modelling B has been developed. A thorough description of many of them can be found in [13]. Only a small amount of the models proposed in the literature has been implemented in general-purpose CFD packages. The reason is that the implementation of some of the models is often difficult, or impossible, to combine with the general framework of the code or the discretization practices the code is based on.
A common approach to SGS-modelling it to employ the Boussinesq assumption, which is the hy-pothesis that the SGS-stress can be modelled in a structurally similar way to the viscous stress, [13].
The analogous idea is also used in most RANS turbulence models. Mathematically it can be ex-pressed as,
B =1
3Tr(B)I + νsgs ∇u + ∇Tu
, (2.11)
where Tr(B) stands for the trace of the tensor B, I is the identity matrix, and νsgsis the so called SGS viscosity which, in turn, is to be computed from the filtered velocity field.
Assuming (2.11), the task is then to obtain a way of calculating νsgs. In order to be able to do that, one has to adopt the hypothesis, that a characteristic length scale lsgsand time scale tsgs
are sufficient to describe the subgrid scales, [13]. Then, based on dimensional grounds, the SGS viscosity can be calculated as
νsgs∼lsgs2 tsgs
= usgslsgs, (2.12)
where usgsis the corresponding velocity scale. A natural choice for lsgsis the filter cut-off width ∆.
The choice of usgsis less obvious and different models use different approaches.
Here we use a model based on solving a transport equation for the subgrid turbulent kinetic energy ksgs, which was proposed independently by several researchers, [2, 6, 14, 18, 19, 15]. The natural choice for the characteristic velocity scale is then usgs=p
ksgs. The transport equation for ksgsis,
∂ksgs
∂t +∂uiksgs
∂xi
= 2νsgs|Dij|2− Ce
k3/2sgs
∆ + ∂
∂xi
νsgs
∂ksgs
∂xi
+ ν∂2ksgs
∂xi∂xi
, (2.13)
7
νsgs= Ck∆ ksgs, (2.14) where, Ck= 0.094, is another model constant.
Physically the four terms on the right-hand side of (2.13) represent, respectively, the production of turbulence by the resolved scales, turbulent dissipation, turbulent diffusion, and viscous dissi-pation. More details on the derivation of (2.13) and the employed modelling assumptions can be found on page 128 in [13]. Note that, ksgs= Tr(B)/2, and therefore can be used for calculating the isotropic subgrid stresses.
When defined as (2.4), νsgsdoes not exhibit correct behaviour in the limit y→ 0, where y is the distance from a wall. To rectify this problem a damping function can be employed. In the simulations presented here the van Driest damping function is used. It has the following form.
f = κ C∆
1− e−y+/A+y
(2.15)
Here, κ = 0.41, is the von Karman constant, C∆= 0.158, and, A+= 26.
2.4 Discretization and interpolation methods
As many other contemporary CFD codes, OpenFOAM is based on the finite volume method for disretizing partial differential equations. The method itself, and its application to the Navier-Stokes equations, are thoroughly described in a large number of monographs dedicated to numerical methods and fluid flow modelling. In particular the books, [16, 1], provide a good introduction to the method in the context of CFD, and a detailed OpenFOAM-oriented discussion can be found in, [4, 12, 17]. In this section only a short overview of the method will be presented.
The finite volume method is based upon dividing the computational domain into many small non-intersecting polyhedra called control volumes (CVs). Different approaches exist, but in OpenFOAM all variables are stored at the centroids of the control volumes. The value at the centroid represents the whole CV. This can be easily shown to be a second order accurate approximation, see [1].
Derivation of the descretized form of the equations begins with integrating the original equations over a control volume and a time interval ∆t. Then, where it is applicable, the Gauss theorem is used to convert volume integrals into surface integrals.
The procedure can be illustrated using a convection-diffusion equation for some quantity φ,
∂φ
∂t+∇ · (uφ) = ∇ · (Γ∇φ), (2.16)
which is integrated to, Zt+∆t
t
∂
∂t Z
CV
φ dV + I
SCV
φu· n dS − I
SCV
Γ∇φ · n dS
dt = 0. (2.17)
The equation above is an integral form expressing the conservation of the quantity φ. If this is summed over all the control volumes a conservation law for the whole domain is obtained. As a consequence, the finite volume method is intrinsically conservative, which is one of its most attractive properties.
8
In order to get algebraic equations, the volume integral in (2.17) is approximated as, φPVCV, while surface integration is approximated as a sum over the faces of the control volume:
I
SCV
φu· n dS ≈X
f
(uf· n) Sfφf, (2.18)
I
SCV
Γ∇φ · n dS ≈X
f
(∇fφ· n) SfΓf. (2.19)
The index f in the above equations stands for the value in the centroid of the face with the corresponding index. These values are unknown and have to be interpolated using the values in the centroids of the cells.
Choosing spatial interpolation and time-marching schemes has to be done with care in order to have a good balance between accuracy and stability. In the simulations presented here a second order backward differencing scheme was used for time marching. To calculate the time derivative in (2.17), the scheme uses the unknown value of φ from the current time-step, φn, and the values from the two preceding time-steps, φn−1and φn−2,
∂φi
∂t ≈
3
2φn− 2φn−1+12φn−2
∆t , (2.20)
where ∆t is the time-step. The time integrals of the convective and diffusive terms in (2.17) are approximated by the (to be determined) values at the current step, which makes the time-marching scheme fully implicit. Additional information on the second order backward differencing scheme can be found in [4]. The OpenFOAM-keyword for this scheme is backward.
The time-step ∆t was chosen to be small enough to keep the Courant number below one. For meshes M1 and M2 the value of, ∆t = 0.2s, was sufficient, but for the finest mesh M3 it had to be lowered down to 0.08s.
The face-fluxes of momentum were calculated using a linear interpolation scheme (referred to as linear in OpenFOAM). The same scheme was also used to evaluate the values of the gradients in the centroids of the faces.
The face-fluxes of the subgrid scale turbulent kinetic energy were calculated using a TVD inter-polation scheme based on upwind and central differencing, (linearLimited 1 in OpenFOAM). The scheme is based on a flux limiter of the form, max(min(2r, 1), 0), where r is the ratio of successive gradients.
2.5 Solver algorithm
The solver pimpleFoam, provided as part of OpenFOAM, was used to solve the equations derived in sections 2.2 and 2.3. The solver is capable of treating a variety of different flow problems and employ different types of turbulence modelling approaches. Here we give a short overview of the principal components of the underlying algorithm. More information can also be found in [4].
9