12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

WALL-MODELLED LES USING HIGH- AND LOW-ORDER CFD CODES:APPLICATION TO A FLAT-PLATE BOUNDARY LAYERTimofey Mukha
SimEx/FLOW, Engineering Mechanics
KTH Royal Institute of Technology
SE-100 44 Stockholm, Sweden
tmu@kth.se

Philipp Schlatter
SimEx/FLOW, Engineering Mechanics
KTH Royal Institute of Technology
SE-100 44 Stockholm, Sweden
pschlatt@mech.kth.se

ABSTRACT

These software vary greatly in the numerical approach to solving the LES equations, OpenFOAM being based on finite
volume discretization and Nek5000 on the spectral element
method. The latter has the advantage of being high order,
whereas using finite volumes offers the benefit of generality
with respect to the shape of the computational grid cells.
In OpenFOAM, wall modelling for LES has been previously developed in Mukha et al. (2019). But for Nek5000 thedevelopment is ongoing, and we are the first to demonstrate theapplication of WMLES to a developing boundary layer using
this solver. The presented comparison allows to highlight boththe advantages and the difficulties of conducting WMLES in ahigh-order code.

Results from wall-modelled large-eddy simulations ofa zero-pressure-gradient flat-plate turbulent boundary layer
in the range Reθ ≈ [4000, 12000] are reported. The simulations are performed using a low- and a high-order code:
OpenFOAM® and Nek5000, respectively. For the latter, such
simulations have previously not been reported in the literature. Structured hexahedral meshes are used, with two levels of refinement. As an important aspect in the wall modelling methodology, we use a temporally varying wall viscosity in order to enforce the wall shear stress. An equivalent
inflow generation procedure is used for both codes, allowing
for a more fair comparison. Results from Nek5000 simulations are generally more accurate. Both the skin friction andthe profiles of velocity statistics are in good agreement withreference data. For Nek5000, this is an important milestone in
the development of wall modelling capabilities for this solver.
The results from OpenFOAM simulations exhibit a significant
over-prediction of the skin friction, which has not been previously reported in the literature. Further investigation of thesimulation methodology is necessary to find the cause of the
problematic behaviour.

SIMULATION SETUP
A summary of the simulation parameters is provided ini ×Table 1. The computational domain is a box of size 150δ99
i × 20δ i , where δ i is the thickness of the inflow bound20δ99
99
99
ary layer. The applied boundary conditions are as follows.
The Synthetic Eddy Method (SEM) (Jarrin et al., 2006; Poletto et al., 2013) is used to generate the inflow turbulence,corresponding to a flat-plate boundary layer at Reθ ≈ 4000.
The profiles of the mean velocity, turbulent kinetic energy andturbulent length scale necessary for the SEM are obtained fromthe DNS of Schlatter & Örlü (2010). To avoid introducing velocity fluctuations that are too large, we found it necessary to
avoid generating the synthetic eddies at y+ / 20, where y is
the wall-normal coordinate. This location is below the firstoff-wall grid node, which is located at y+ ≈ 37 on the finest
grid. The inflow generator is adopted from Hufnagel et al.
(2018), with the implementation available in Nek5000. For
the OpenFOAM simulations, a dummy Nek5000 simulation
with a just few elements in the streamwise direction was used.
In this simulation, velocity signals were probed at all the locai . This data
tions of the OpenFOAM inflow meshes for y < 3δ99
was then converted to an appropriate format for OpenFOAM
to read from disk. As a result, both solvers use exactly thesame inflow generation procedure.
At the top boundary a mixed boundary condition is used,

INTRODUCTION
In the context of large-eddy simulation (LES) wall modelling typically refers to introducing additional equations that
couple the flow quantities in the boundary layer’s outer region
with the value of the wall-shear stress. In principle, this allows
to leave the inner region of the turbulent boundary layer (TBL)
unresolved by the grid, leading to more favourable scaling of
the grid size with the Reynolds number of the flow. While,in spite of these savings, wall-modelled LES (WMLES) remains computationally expensive, a surge of interest in WMLES could be observed in the recent years, with many newcontributions to the literature. For a review of the topic thereader is referred to (Bose & Park, 2018; Larsson et al., 2016).
In this work, we present results from WMLES of a flatplate turbulent boundary layer. The focus is on comparing thepredictive accuracy of two solvers: OpenFOAM and Nek5000.

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
assigning the streamwise velocity to the freestream value U0 ,the lateral velocity to 0, and using a homogeneous Neumanncondition for the vertical velocity.
At the outflow, a homogeneous Neumann condition for
the velocity is used and a homogeneous Dirichlet condition
for the pressure.
At the wall, the boundary condition is determined by thewall model. Since wall modelling is the primary concern of
this work, it is discussed separately in more detail in the following section.
Unfortunately, the subgrid-scale (SGS) models available
for selection in both codes are not the same. For this reason,in OpenFOAM, the WALE model (Nicoud & Ducros, 1999) isused, whereas in Nek5000 the Vreman (2004) model is employed. However, both models are algebraic and therefore
similar in terms of computational complexity. In the case of
OpenFOAM, WALE has been shown to be a good choice for
WMLES in previous works Rezaeiravesh et al. (2019). Conversely, for Nek5000 the model selection is still actively researched. Preliminary results from channel flow simulations
using the Vreman model have shown good agreement with reference data.
For OpenFOAM, the selection of the scheme for interpolating the convective fluxes is important. Here, a scheme linearly blending linear interpolation with second-order upwind
interpolation is used. The weight of the linear scheme is 0.75.
This scheme is second-order accurate but quite diffusive dueto the upwinding. It would be possible to use linear interpolation only, however the blended scheme is known to keep simulations stable even for highly unstructured grids. The OpenFOAM setup is thus representative of what one would use in
an industrial flow simulation. For time-integration, a secondorder backward-differencing scheme is used.
In Nek5000, the most important choice is the order of theinterpolating polynomials. Here, we use order 7, which is awell-tested selection. This refers to the basis used for the velocity field. For pressure, the polynomials are lower by twoorders as part of the pressure-velocity coupling approach referred to as PN -PN−2 , which allows to avoid spurious modesin the pressure field (see Deville et al. (2012) for details). Different time-integration schemes are applied to different termsof the momentum equation. For the viscous term, a thirdorder backward-differencing scheme is used, whereas thirdorder extrapolation is applied to the convective and source
terms, see Karniadakis et al. (1991).
Finally, in both codes adaptive time-stepping was used tomaintain a CFL number of ≈ 0.5. Time-averaging was peri /U , which is equivalent to about 27 eddy
formed over 700δ99
0
turnover times.
Two meshes were employed for each code, referred to as
M1 and M2. In the region occupied by the boundary layer,
the meshes are isotropic, and the resolution is defined by thei : 1 and 2 for M1
number of Nek5000 elements covering δ99
and M2, respectively. The same resolution is maintained up to
i after which the grid is coarsened in the wall-normal
y = 4δ99
direction. In OpenFOAM, the corresponding grids are constructed considering that each Nek5000 element contains 7 intervals along a single edge. Therefore, 7 and 14 cells coveri in the two OpenFOAM grids. Additionally, we make use
δ99of the fact that OpenFOAM supports unstructured grids byicoarsening the grid in the wall-parallel direction at y = 4δ99
for the denser grid. An illustration of the M2 grid for bothcodes is given in Figure 1. It should be noted that even thei by
M2 grid is coarse as one would ideally want to resolve δ99
at least 4 elements (or 28 cells with OpenFOAM) across the

Figure 1. A portion of the computational mesh M2 at the inlet of the domain OpenFOAM (left) and Nek5000 (right).

whole range of Reθ . However, the corresponding grid size becomes prohibitively large. In the region occupied by the TBLi ) the M1 and M2 meshes have, respectively, ≈ 4 · 106(y < 4δ99
and ≈ 32 · 106 degrees of freedom (DOFs). The latter refers to
the number of cells for OpenFOAM and number of nodes for
Nek5000.
While the resolution of the outer scales is defined by thegrid construction strategy, the resolution in wall units is an outcome of the simulation. This result is shown in Figure 2. The
top plot shows the distribution of y+
1 . Clearly, all meshes arevery coarse by this measure. Since ∆x+ = ∆z+ it suffices toconsider one of these quantities. The bottom plot in Figure 2
shows ∆x+ , and to make the plot clear only a small range inx is depicted. The modulation of the values by the change in
the inner length scale can be inferred from the behaviour ofy+
1 . One clearly sees the large effect of the predefined distribution of nodes within the element in Nek5000. On the M2
grid, ∆x+ is as low as 40 at the inter-element boundaries. It
is, however, not entirely clear how the effective resolution of
the Nek5000 grid should be defined. One could argue that the
maximum spacing is the relevant measure. Regardless, if one
were to decrease the grid size by a factor of 2 (which is desirable for better resolution of the outer scales), ∆x+ would startgetting close to that of a wall-resolved simulation. This highlights the necessity of considering high Reynolds number flow
for WMLES.
Table 1.

2

Simulation parameters.

Parameter

ValueiInflow TBL thickness, δ99

1m

Freestream velocity, U0

1 m/s

Kinematic viscosity, ν

2.86 · 10−5 m2 /sInflow Re-number, Reθ≈ 4060

Domain lengthi
150 m = 150δ99

Domain heighti
20 m = 20δ99

Domain widthi
20 m = 20δ99

Time averaging length

i /U
700 s = 700δ990iGrid size, y < 4δ99

M1: ≈ 4 · 106 , M2: ≈ 32 · 106

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

175
150
y1+

the rate of strain, the necessary τw value is recovered.
In this work, we use the latter approach, which is now
presented in more detail. Let νt be the value of SGS viscosity
at the wall, and Si j the components of the rate-of-strain tensor.
Then, the magnitude of the wall shear stress is computed as

M1, OF
M2, OF
M1, NEK, 2nd node
M1, NEK, 3rd node
M2, NEK

200

125

q
2 + S2|τw | = (ν + νt ) Sxy
zy

100

(1)

75

It follows, that given |τw | from the wall model and the valuesof Si j as computed at the wall, we can set νt according to

50
0

50

100

150

q
2 + S2 − ννt = |τw |/ Sxy
zy

x

(2)

300

to prescribe |τw | at the boundary.
In OpenFOAM, this approach has been well tested in several works, e.g. Mukha et al. (2019). Conversely, in Nek5000
it has not been previously applied, and the comparative benefits with respect to the Neumann condition is a topic we are
actively investigating, but will not focus on here. In the context of finite differences, the two approach have been analysedby Bae & Lozano-Durán (2021).
The import role of h in the accuracy of the predicted τw
has been initially highlighted in Kawai & Larsson (2012) andlater confirmed in other works, e.g. Rezaeiravesh et al. (2019);
Frère et al. (2017). The consensus is that, independent of theunderlying numerics, sampling from the first off-wall solution
point leads to suboptimal results. More generally, sampling
has to take place from a node where the velocity signal is accurate since any errors in u propagate directly into τw . In thiscontext, it is observed that in the first off-wall node the error
in u is always significant, see Kawai & Larsson (2012).

250

∆x+

200
150
100
50
148.00 148.25 148.50 148.75 149.00
x+
Figure 2. The size of the grid in wall units: y+
1 (top) and ∆x(bottom). Note that ∆z+ = ∆x+ .

WALL MODELLING
The role of the wall model is to provide the correct valueof the wall shear stress τw for each boundary node. This task
can be split in two: obtaining the sought value of τw and ingesting it into the simulation by means of a suitable boundary
condition.
The general procedure for obtaining τw is to use a relationship coupling the stress with the LES solution at somepoint away from the wall. The distance to this point is an important parameter, which we will refer to as h. The simplest
wall model can be constructed using the law of the wall, suchas as Spalding’s law (Spalding, 1961), which is used here. The
law provides an algebraic expression connecting the streamwise velocity with the wall shear. Generally, this can be expressed as F(u+ , y+ ) = 0. At each wall node, the value of u isthen sampled at distance y = h, these values are inserted into
F and the equation is solved for the friction velocity uτ using
Newton-Raphson iteration. The procedure is repeated at eachtime-step.
The value of τw is formed as a product of viscosity andthe shear components of the rate of strain tensor evaluated at
the boundary. Therefore, one can generally distinguish two
main approaches to setting the obtained stress as a boundary
condition. One is to use a Neumann boundary condition for thevelocity field, effectively resulting in a slip velocity at the wallof such magnitude that the resulting rate of strain produces
the desired value of τw . Alternatively, one can leave the noslip condition for velocity intact and instead manipulate theviscosity at the wall in such a way that the when it multiplies

RESULTSThe figures in this section present results from 5 simulations. OpenFOAM simulations are denoted as OF in the legend, and Nek5000 simulations with NEK. For each code, asimulation on each of the meshes M1 and M2 is performed
with h set to the 2nd and 3rd off-wall node, respectively. This
corresponds to approximately the same distance from the wall(see Figure 5). Additionally, a Nek5000 simulation on the M1
grid with h set to the 3rd node is performed to highlight theimportance of the selection of h.
To demonstrate the level of resolution of turbulent structures, Figure 5 shows a snapshot of the streamwise velocity
taken at the outlet of the M2 grids. Qualitatively, the turbulence in the outer layer appears well-represented in both
Nek5000 and OpenFOAM simulations. However, one should
keep in mind that at the outlet the grid size is smallest relativeto δ99 . Even so, the Nek5000 solution clearly exhibits numerical artefacts at element boundaries. Particularly the near-wallelement is easily distinguished in the solution. This indicatesthat we can expect an improvement of the results using finer
grids.
We begin the quantitative analysis of the results by looking at the performance of the wall modelling, which manifests
itself in the values of the skin friction coefficient, c f . The obtained values are shown in Figure 4. For reference, we use
the wall-resolved LES (WRLES) data by Eitel-Amor et al.
(2014), which, however, covers only a part of the Reθ range ofthe WMLES. Therefore, the power law by Rezaeiravesh et al.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 3.

Streamwise velocity distribution at the outlet. Simulations on mesh M2, Nek5000 (top) and OpenFOAM (bottom)(2016) is also depicted. The agreement between the WRLES
and the power law is quite good.
The results from OpenFOAM simulations exhibit a significant over-prediction of c f , although on the M2 grid the results
are improved at high Reθ . The accuracy is worse to what waspreviously reported at lower Re-number simulations in Mukha
et al. (2021). There are methodological differences, which
could explain that, particularly the fact that in Mukha et al.
(2021) an unstructured grid adapted to the growth of the TBL
was employed. The latter implies that the resolution of the
TBL remains constant with x and that h is adapted to remain
at the same fraction of the TBL thickness. It may therefore be
possible to improve the OpenFOAM results, and further investigation is required.
The curves from Nek5000 simulations are more faithful
to the reference power law. A very significant effect of h canbe observed, however. The choice between 2nd or 3rd off-wall
grid node on the M1 grid results in a shift of the c f values, withsignificant under-prediction when the 2nd node is used. By
contrast, the effect of grid resolution is more or less negligible
per se—as long as an appropriate h is selected.
The profiles of the mean streamwise velocity are analysed
next, see Figure 5. The location corresponds to Reθ ≈ 8183. It
is instructive to see whether the expected connection between
the accuracy of c f and the velocity prediction at the sampling
point can be found. Indeed, we see that the OpenFOAM curves
over-predict the reference profile, including the points corresponding to h. Looking at the Nek5000 profiles from grid M1
it is clear that the 3rd node is indeed in better agreement withthe reference than the second, the latter lying below the curve.
This is consistent with the c f values. However, it is likely thatthe error in hui is only part of the explanation for the c f errors. For example, the Nek5000 simulation on the M2 grid also
over-predicts the reference and the accuracy is only marginallybetter than that of the OpenFOAM value on the same grid. Yet
the difference in c f much larger. Finally, in inner scaling, thevelocity values at the sampling points coincide perfectly withthe reference log-law, as expected.
It is worth noting that the value of h has more or less no

M1, OF
M2, OF
M1, NEK, 2nd node
M1, NEK, 3rd node
M2, NEKPower law
Ref.

3.6
3.4
cf · 103

3.2
3.0
2.8
2.6
2.4
4000

6000

8000
Reθ

10000

12000

Figure 4. Development of the skin friction coefficient, c f .
Reference data from Eitel-Amor et al. (2014).

influence on hui, the Nek5000 profiles on the M1 grid overlap. So, the coupling between the sampled velocity and c f isunidirectional. This has previously been observed in OpenFOAM, see Rezaeiravesh et al. (2019), however, one could
expect a stronger coupling in Nek5000 since the velocity values are more strongly coupled via the underlying polynomial
basis.
The accuracy of the Nek5000 solutions are clearly superior. However, the solution on the M2 mesh is worse than on
M1 indicating that some fortuitous error cancellation may beat play. Interestingly, the plot in inner scaling shows no evidence of log-layer mismatch for the OpenFOAM solution, thediscrepancy with the reference is instead further away from thewall.
Finally, we consider the second-order statistics, particularly the mean turbulent kinetic energy, hki and the mean turbulent shear stress hu0 v0 i. The profiles are shown in Figure 6.
The accuracy of the OpenFOAM profiles is very similar to thatreported in previous works, e.g. Mukha et al. (2019). The

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
shear stress is well-predicted excluding small non-physicalpeak near the wall. For hki one observes the typical pattern ofover-prediction near the wall and then under-prediction closer
to the free stream. In the Nek5000 solutions one clearly see
the boundaries of the elements: they coincide with the spikes
in hki. Clearly, on the M1 mesh the TBL is covered by 2 elements at the considered streamwise position. This resolution
is too coarse to properly resolve the turbulent fluctuations andthey are over-predicted. Note that even here the effect of h is
negligible. On the M2 mesh the results are better and close inaccuracy to that of OpenFOAM.

1.0

hui/U00.8
0.6
M1, OF
M2, OF
M1, NEK, 2nd node
M1, NEK, 3rd node
M2, NEK
Ref.

0.4
0.2

CONCLUSIONS

0.0
0.0

0.5

The most import outcome of this study is a successful
WMLES of a spatially-evolving boundary layer in Nek5000.
It is clear that, with a careful selection of simulation parameters, good agreement with reference data can be obtained. The
simulations showed that the sampling height h should especially be chosen with care. The fact that Nek5000 outperforms
OpenFOAM in terms of accuracy on equivalent grids is very
encouraging. Using viscosity to enforce the wall shear stressappears to be a robust way to transmit the values predicted by
the wall model into the simulation.
Several methodological question remain open: i) The necessary grid resolution to obtain convergence of statistics in theouter layer; ii) Length-scale computation for the subgrid viscosity, and selection of the model; iii) Selection of h, particularly its non-uniform distribution in space. All of these are
currently the subject of our research.
For OpenFOAM, an important step is to try to apply thesame unstructured meshing strategy as in (Mukha et al., 2021)
to see if that alleviates the issues with c f predictions. Further,simulations on denser grids have to be performed to see if better agreement with the mean velocity profile can be obtained.

1.0
y/δ99

30.0
Log-law

27.5
25.0

hu+ i

22.5
20.0
17.5
15.0
12.5
10.0
101

102

103
y+

Figure 5. Mean velocity profiles at Reθ ≈ 8183 in outer scaling (top) and inner scaling (bottom). Stars show the locationof the sampling point. Reference data from Eitel-Amor et al.
(2014). The log-law parameters are κ = 0.387 and B = 4.21
where selected according to Yamamoto & Tsuji (2018).

REFERENCES
Bae, H. J. & Lozano-Durán, A. 2021 Effect of wall boundary conditions on wall-modeled large-eddy simulation in a
finite-difference framework. Fluids 6 (112).
Bose, S. T. & Park, G. I. 2018 Wall-modeled large-eddy simulation for complex turbulent flows. Annual Review of Fluid
Mechanics 50, 535–561.
Deville, E. H., Fischer, P. F. & Mund, E.H 2012 High-Order
Methods for Incompressible Flow. Cambridge University
Press.
Eitel-Amor, G., Örlü, R. & Schlatter, P. 2014 Simulation andvalidation of a spatially evolving turbulent boundary layer
up to Reθ = 8300. International Journal of Heat and Fluid
Flow 47, 57–69.
Frère, A., de Wiart, C. C., Hillewaert, K., Chatelain, P. &Winckelmans, G. 2017 Application of wall-models to discontinuous Galerkin LES. Physics of Fluids 29, 085111.
Hufnagel, L., Canton, J., Örlü, R., Marin, O., Merzari, E. &Schlatter, P. 2018 The three-dimensional structure of swirlswitching in bent pipe flow. Journal of Fluid Mechanics
835, 86–101.
Jarrin, N., Benhamadouche, S., Laurence, D. & Prosser, R.
2006 A synthetic-eddy-method for generating inflow conditions for large-eddy simulations. International Journal of
Heat and Fluid Flow 27 (4), 585–593.
Karniadakis, G. E., Israeli, M. & Orszag, S. A. 1991 Highorder splitting methods for the incompressible Navier-

1.4
M1, OF
M2, OF
M1, NEK, 2nd node
M1, NEK, 3rd node
M2, NEK
Ref.

1.2

y/δ99

1.0
0.8
0.6
0.4
0.2
0.0
−2

0

2

4

6

8

103 · hu0 v 0 i/U02 , 103 · hki/U02

Figure 6. Profiles of the mean turbulent kinetic energy, hki,and the turbulent shear stress, hu0 v0 i at Reθ ≈ 8183. Reference
data from Eitel-Amor et al. (2014).

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
conditions for les. Flow, Turbulence and Combustion 91 (3),519–539.
Rezaeiravesh, S., Liefvendahl, M. & Fureby, C. 2016 On grid
resolution requirements for LES of wall-bounded flows. In
ECCOMAS Congress 2016. Crete, Greece.
Rezaeiravesh, S., Mukha, T. & Liefvendahl, M. 2019 Systematic study of accuracy of wall-modeled large eddy simulation using uncertainty quantification techniques. Computers
and Fluids 185, 34–58.
Schlatter, P. & Örlü, R. 2010 Assessment of direct numerical simulation data of turbulent boundary layers. Journal of
Fluid Mechanics 659, 116–126.
Spalding, D. B. 1961 A single formula for the ‘law of the wall’.
Journal of Applied Mechanics 28 (3), 455–458.
Vreman, A. W. 2004 An eddy-viscosity subgrid-scale modelfor turbulent shear flow: Algebraic theory and applications.
Physics of Fluids 16 (10), 3670–3681.
Yamamoto, Y. & Tsuji, Y. 2018 Numerical evidence of logarithmic regions in channel flow at Reτ = 8000. PhysicalReview Fluids 3 (1), 012602.

Stokes equations. Journal of Computational Physics 97 (2),414–443.
Kawai, S. & Larsson, J. 2012 Wall-modeling in large eddy
simulation: Length scales, grid resolution, and accuracy.
Physics of Fluids 24 (1), 015105.
Larsson, J., Kawai, S., Bodart, J. & Bermejo-Moreno, I. 2016
Large eddy simulation with modeled wall-stress: Recent
progress and future directions. Mechanical Engineering Reviews 3 (1), 1–23.
Mukha, T., Bensow, R.E. & Liefvendahl, M. 2021 Predictive accuracy of wall-modelled large-eddy simulation on
unstructured grids. Computers & Fluids 221, 104885.
Mukha, T., Rezaeiravesh, S. & Liefvendahl, M. 2019 A library for wall-modelled large-eddy simulation based on
OpenFOAM technology. Computer Physics Communications 239, 204–224.
Nicoud, F. & Ducros, F. 1999 Subgrid-scale stress modellingbased on the square of the velocity gradient tensor. Flow,
Turbulence and Combustion 62 (3), 183–200.
Poletto, R., Craft, T. & Revell, A. 2013 A new divergence freesynthetic eddy method for the reproduction of inlet flow

6