K. Hanjali´ , Y. Nagano and S. Jakirli´ (Editors)
K. Hanjali´ , Y. Nagano and S. Jakirli´ (Editors)
Wall-Functions and Boundary Layer Response to Pul-
sating and Oscillating Turbulent Channel Flows
D. Panara1 , M. Porta2 ,R. Dannecker1, and B. Noll1
DLR Deutsches Zentrum f¨ r Luft- und Raumfahrt e. V.
Institut f¨ r Verbrennungstechnik, Stuttgart, Germany, [email protected],
[email protected], [email protected]
CERFACS Centre Europeen de Recherche et de Formation Avancee en Calcul Scientiﬁque,
Toulouse, France, [email protected]
Abstract — Some limitations of the classical Wall-Function approach for the near-wall boundary layer
treatment in LES and URANS are presented for isothermal oscillating and pulsating channel ﬂows.
Despite their simple geometry, pulsating and oscillating ﬂows are interesting unsteady ﬂow test cases
representative of many industrial components. A special attention will be focused on the unsteady wall-
shear stress prediction since it is also an indirect measure of the unsteady wall-heat transfer. A correct
evaluation of the unsteady wall heat transfer is in fact critical, for example, in combustion chamber
applications when ﬂow unsteadiness due to ﬂame instability occurs and in general in each unsteady ﬂow
situation where thermo-acoustic phenomena play an important role.
With the continuous growing of the computational resources, the use of unsteady numerical
simulations is becoming a strategic tool in designing complex industrial components. Many
examples of that can be found in turbomachinery and combustion chamber applications. Com-
mercial CFD softwares are widely used in industry and a critical factor for the design develop-
ment resides in an optimal balance between numerical costs and prediction accuracy. Unsteady
Reynolds averaged Navier-Stokes (URANS) and, in the last years, Large Eddy Simulations
(LES) have been used increasingly by industry for the study of critical components. For in-
ternal turbulent ﬂow simulations, in order to save computational time, the wall-boundary layer
resolution is often neglected. In both, the LES and URANS approach, the tendency is to model
the inner-part of the boundary layer by means of so-called wall functions. The use of wall func-
tions signiﬁcantly reduces the number of computational points required in the boundary layer
region and for this reason it is widely used. However, it is questionable whether in pulsating and
unsteady ﬂow conditions the hypothesis under which the wall functions have been developed
are still valid. A deeper investigation on the accuracy of the use of wall functions for unsteady
turbulent ﬂow applications is therefore needed.
Despite their simple geometry, pulsating channel ﬂows are representative of many interesting
industrial conﬁgurations. In the present work, the accuracy of the use of wall functions in
LES and URANS is investigated by means of such test cases. Our attention is focused on the
unsteady wall-shear stress prediction since it is important also as an indirect measure of the
unsteady wall heat transfer.
In the present work two very different CFD research solvers are used. OpenFoam  has been
chosen for the URANS calculations and AVBP  for the LES calculations. These two codes
differ not only for the discretisation schemes used but also for the mathematical formulation
2 Turbulence, Heat and Mass Transfer 5
of the problem. A brief description of the codes and governing equations will be given in the
next sections. For the near wall treatment, both codes employ a similar wall function approach.
Details about the LES code and its wall functions implementation are given in section (2.) and
(2.1.). Details of the URANS solver are given in section (3.) and the implementation of wall
functions is brieﬂy described in section (3.1.). Detail on code validations and their results for
oscillating and pulsating ﬂow are given in . In section (4.) the most interesting validation
results are reported. Finally in the last part of the paper, the URANS results are discussed to-
gether with comparable LES pulsating channel ﬂow simulations and experimental results using
characteristic non-dimensional parameters. Main attention is paid on wall-shear stress and its
phase shift with respect to the free stream velocity.
2. LES Solver
DNS and LES computations object of this work, were performed using the AVBP code de-
veloped by CERFACS and IFP (Institut Francais du P´ trole). AVBP is a parallel CFD code
that solves the laminar and turbulent compressible Navier-Stokes equations in two and three
space dimensions on unstructured and hybrid grids. The data structure of AVBP employs a cell-
vertex ﬁnite-volume approximation. The basic numerical methods are based on a Lax-Wendroff
or a Finite-Element type low-dissipation Taylor-Galerkin discretisation in combination with a
linear-preserving artiﬁcial viscosity model. In this paper only the former was used, the study of
the inﬂuence of the numerical scheme on the solution remains open for future works. The time
discretisation is explicit making use of a Runge-Kutta multi-stage time stepping. For turbu-
lent compressible ﬂows, AVBP solves the LES formulation of the Navier-Stokes equations. In
the LES approach, the governing equation are ﬁltered in space before discretising and solving.
Additional unresolved terms appear in the convective ﬂuxes. For the Reynolds stresses we have:
τij t = −¯(uiuj − ui uj )
ρ ˜˜ (1)
The over-bar represents a ﬁltered numerically resolved quantity and the tilde represents a mass-
weighted Favre average.
ρ and ui represent respectively the density and the ith component of the velocity vector u. The
unresolved sub-grid scale (SGS) terms are generally closed using the following formulation:
τij t = 2¯νt Sij − τll t δij
1 ∂ ui ∂ uj
˜ ˜ 1 ∂ uk
Sij = ( + )− δij (3)
2 ∂xj ∂xi 3 ∂xk
In our computations two different closure are used: the classical Smagorinsky model 
νt = (CS ∆)2 2Sij Sij (4)
and the WALE model 
(sd sd )3/2
νt = (Cw ∆)2 (5)
(Sij Sij )5/2 + (sd sd )5/4
where Cw and CS are model constants ( Cw = 0.4929 and CS = 0.18 ), ∆ is the characteristic
ﬁlter length and
D. Panara et al. 3
Figure 1: Typical velocity proﬁle near the wall and notation used for near-wall quantities.
1 2 1 2
sd = (˜ij + gji) − gkk δij
ij g ˜2 ˜ (6)
where gij denotes the resolved velocity gradient. The WALE model was developed for wall
bounded ﬂows in an attempt to recover the scaling laws of the wall without using the wall
2.1. Wall Functions Implementation
The wall law implemented in AVBP is presented in detail by Schmitt  and is here shortly
described. As mentioned above AVBP uses a cell-vertex scheme. All quantities are stored at
the cell-corners. For the calculation of the viscous ﬂuxes AVBP needs the shear stresses at the
cell boundary and heat ﬂuxes. Imposing the appropriate values of velocity and temperature
at the boundary plus the correct ﬂuxes, using the wall-law formulations, constrains the ﬂow
too much and leads to oscillatory solutions. The strategy used by Schmitt is then to impose the
wall-shear stress τw (and heat ﬂux qw ) at the boundary using the wall-function approach without
ﬁxing the value of velocity and temperature at the cell corners (u1 and T1 in Figure (1)). Only
the normal component of the velocity at the wall is imposed to vanish for continuity reasons.
This is equivalent, as shown in Figure (1), to imagine the real wall boundary shifted by a small
distance δw away from the computational domain. Assuming that the shift is small compared
to the distance between the wall and the ﬁrst point in which the wall-function is evaluated
(δw ≪ yw ), it can be neglected when computing the wall distance. The wall shear stress is then
imposed at the boundary following the steps in Table (1).
3. URANS Solver
The URANS calculation where performed using the OpenFoam solver. OpenFoam (Open Field
Operation and Manipulation) is a CFD toolbox that uses ﬁnite volume numerics to solve systems
of partial differential equations ascribed on any 3D unstructured mesh of polyhedral cells. The
top-level code used for our computations is the standard OpenFoam solver turbFoam, a transient
solver for incompressible, turbulent ﬂow. turbFoam solves the URANS equation for a turbulent
ﬂuid ﬂow using a robust, implicit, pressure-velocity, iterative algorithm based on the PISO
scheme  (Pressure-Implicit with Splitting of Operators).
4 Turbulence, Heat and Mass Transfer 5
Step 1 Compute uτ iteratively from Eq. (7) or (8) with y + = uν y and u+ =
Input values: u = u2 , ν = ν(T1 ), y = ∆y.
Step 2 Compute τw from uτ =
with ρ = ρ1 .
Step 3 Apply τw and advance ﬂow equations.
Step 4 Set normal component of u1 to zero and go to Step 1.
Table 1: Working principle of the wall-function boundary condition.
For the URANS simulations, two series of near wall treatments have been employed: the Wall
Function approach as explained in section (3.1. together with a k-ǫ High Reynolds Number
model and a resolved boundary layer approach using the k-ǫ Low Reynolds Number model
from Launder and Sharma .
3.1. Wall Functions Implementation
The use of wall functions is based on two important assumptions:
1. The validity of the universal law of the wall
u+ = y + 0 < y+ < 5 (7)
u+ = ln y + + C y + > 40 (8)
2. The assumption of turbulent local equilibrium in the near wall region:
P = ǫ; P = νt ( ∂u )2
where P and ǫ represent the production and dissipation term of the turbulent kinetic en-
ergy transport equation. The turbulent kinetic viscosity νt in the k-ǫ turbulent model
formulation depends on the solution of the two additional transport equations for k and ǫ.
In the above hypothesis, assuming a known distribution of k and ǫ near the wall by solving their
transport equation, it is possible to obtain an expression for τw that can be used as boundary
condition for the solution of the momentum equation.
4. Code Validations
4.1. LES Code Validation
In order to validate the LES code in a turbulent oscillating channel ﬂow application the numer-
ical results are compared with a DNS from Spalart and Baldwin . Two series of near wall
treatments have been employed: the Wall Function approach as explained in section (2.1. and a
so-called Wall-Normal Resolved approach .
For the Wall-Function computations an equally-spaced grid in the three axis direction is used
and a wall function treatment is employed in the near wall region as explained in (2.1. For the
closure of the SGS terms, the classical Smagorinsky model has been used.
D. Panara et al. 5
In the Wall-Normal Resolved computations, the used grid is equally spaced in two of the three
axis dimensions, but stretched in the wall-normal direction in order to resolve the boundary
layer up to a wall unit value (y + ) of the order of 2. The wall unit values are computed using the
maximum velocity value in the cycle. For the closure of the SGS terms, the WALE model has
been used in order to reproduce the asymptotic behavior of the turbulence at the wall .
In the following we report just the main ﬁndings on the behavior of the wall shear stress (τw ).
Details on velocity proﬁles, velocity ﬂuctuations and pressure ﬂuctuations are reported in Pa-
nara et al. . In Figure (2) τw is displayed vs. phase and compared with the DNS results.
The Wall-Function computations seem to better reproduce the amplitude of the wall shear stress
oscillations. We notice a little phase shift between the peak value of τw that is not present in
the Wall-Normal Resolved results. The Wall-Normal Resolved τw predictions seem to underes-
timate the peak value of the wall-shear stress but seem to be more in phase with the DNS data.
The overall behavior of both approaches is quite good for this oscillation frequency.
The most striking discrepancy has been found on the pressure ﬂuctuation prediction. Unfortu-
nately, no direct DNS data are available. The only pressure ﬂuctuations related quantity reported
by Spalart and Baldwin  is the pressure term of the Reynolds-stress transport equation:
1 ′ ∂P ′ ∂P ′
Pterm = Π12 = − u − v′ (10)
ρ ∂y ∂x
In Figure (3) Pterm is reported compared with the DNS data. Only the phase φ = π has been
showed for shortness but similar trends are observed in all the other phases. The ﬁgure shows
that the Wall-Normal Resolved results agree with the DNS data. The Wall-Function results
are not reported in the graph since out of scale. The discontinuous results are probably due to
numerical noise related to the discretisation scheme employed.
In conclusion for the frequency investigated, the differences between the two near wall ap-
proaches are small concerning the wall-shear stress prediction. The pressure ﬂuctuations in-
stead are underestimate by the Wall-Function computation and this could be a critical limitation
for thermo-acoustic applications.
6 Turbulence, Heat and Mass Transfer 5
Figure 2: Wall-shear stress vs. phase, Wall-Normal Resolved and Wall Function Computation.
Figure 3: Pressure term, Reynolds stress budget.
D. Panara et al. 7
4.2. URANS Code Validation
The URANS code is validated against the experimental data of Tardu et al.. The pulsation
frequency is given in relation to the so called dimensionless viscous Stokes layer thickness (ls )
ls = (11)
ω+ = (12)
Where u2 is the skin friction velocity of the relative steady channel case. Two series of compu-
tations were performed using two different near wall approaches. In the High Reynolds Number
(HR) computations, the High Reynolds Number k-ǫ turbulence model is used together with the
wall function wall treatment as explained in (3.1.). For the Low Reynolds Number (LR) com-
putations, the Low Reynolds number k-ǫ turbulence model from Launder and Sharma  has
been used with a wall reﬁned mesh. The points in the direction normal to the channel wall are
stretched using a simple grading algorithm. The instantaneous values of y + are varying during
the unsteady computations but are always between 10 and one. The number of grid points due
to the boundary layer resolution is consequently higher respect to the non stretched HR grid. At
the inlet a turbulent inlet proﬁle is pulsated and a ﬁxed static value of pressure is prescribed at
In Figure (4) the near-wall velocity proﬁles in different phases compared with the experimental
results for ls = 8.1,auc = 0.64 and Uc = 0.169 m/s is reported (See Tardu et al. for
parameter deﬁnitions). The high Reynolds number and low Reynolds number boundary layer
treatments are indicated with HR and LR. The ﬁgure shows a ﬁrst limitation of the wall law
approach. In this ﬂow regime, the magnitude of the pulsations determines the ﬂow reversal
close to the solid walls. This is not captured by the HR model even though there is a quite good
agreement between the HR and LR model far from the wall.
Computing the wall shear stress according to the wall-function formulation and using the fol-
lowing for the LR model:
τw = ρν
¯ ¯ | (13)
a phase shift with respect to the velocity centreline of −8◦ and 33◦ as shown in Figure (5) has
been obtained. Besides the calculations done for ls = 8.1, a variety of further simulations was
done for other values of ls . In Figure (5) the value of wall-shear stress phase shift is shown
as function of the dimensionless viscous Stokes layer thickness ls . The numerical results are
compared with the experimental results reported by Tardu  (symbols in black). The graph
shows the incapability of the wall law approach to predict the experimental wall-shear stress
phase shift. The red and blue vertical lines separate the regimes of quasi-laminar and quasi-
steady boundary layer behavior. Depending on the pulsation frequency, different boundary
layer regimes are experienced. In the quasi-steady regime, the turbulence has time to relax to
the local (in time) equilibrium. The ﬂow can be studied as a succession of steady states and
the wall function assumption seems to be still valid in this ﬂow conditions as shown in Figure
(5). With increasing frequency the turbulence production and dissipation start to show a phase
lag. In this situation a change in amplitude and phase of the wall-shear stress in respect of the
8 Turbulence, Heat and Mass Transfer 5
Figure 4: Instantaneous velocity proﬁles in the presence of reverse ﬂow. Uc = 16.9 cm/s,
auc = 0.64, ls = 8.1.
outer velocity is measured. A Stokes layer, where the effects of the outer ﬂow oscillations are
conﬁned, occurs. The thickness of the Stokes Layer decreases with increasing forcing frequency
and in the inertia dominated or ( quasi-laminar ) regime the Stokes layer resides completely
within the viscous sub layer. In this case a ﬂow solution can be obtained combining the laminar
Stokes solution in the laminar sub-layer with a turbulent plug ﬂow far from the wall.
5. LES and URANS Near Wall Numerical Predictions in Turbulent Pul-
For the comparison between URANS and LES near-wall numerical prediction in turbulent pul-
sating ﬂow we designed a LES pulsating channel case in order to meet the value of ls consid-
ered by Tardu et al.. The computational domain and boundary conditions are analogous to
the validation of the oscillating case. The source term has been modiﬁed as follows:
− = Kosc sin ωt + K (14)
Kosc in ﬁrst approximation has been evaluated using the laminar analytical solution for oscillat-
ing ﬂows (Kosc = Uc ω). K has been chosen, as a ﬁrst approximation, to balance the wall mean
where h is the channel height and in our case the size of the cubic numerical domain chosen.
For the case of ls = 14.14 we considered a value of Uc = 70 m/s and an height of the channel
h = 0.006 m. We obtained then using steady channel correlation a value of τw = 13.88 N/m2
D. Panara et al. 9
and uτ = 3.48 m/s with a channel Reynolds number of 98800. The consequent pulsation
frequency can be sought from the ls deﬁnition and it is around 1140 Hz.
The amplitude of the velocity pulsation (Auc in the Tardu et al. notation) has been chosen equal
to 20 m/s. For different values of ls only the oscillating frequency has been changed in the
evaluation of the source terms. All the computations are then made with the same channel size
(h = 0.006 m) and the computed values of Uc and Auc remain always around 70 and 20 m/s.
For comparison, a DNS channel of h = 0.0015 m with same values of Uc and Auc has also ˜
been computed. In order to match the nominal condition of ls = 14.14 the pulsation frequency
has been set to 1500 Hz. The DNS grid consists of a wall-normal stretched grid with 73x73x9
The results are shown in Figure (5). In comparison to the DNS calculation the wall shear
stress phase shift is underestimate by the Wall-Normal Resolved approach. The Wall-Function
approach instead underestimate the phase shift and is also quite off from the experimental data.
The ﬁgure shows clearly the limitation of the wall function approach not only in URANS but
also in LES. The Wall-Normal Resolved approach is much more accurate and in agreement with
the experimental data.
Figure 5: Wall shear stress phase shift dependence on pulsation frequency
The use of wall-functions for LES and URANS has been investigated paying special attention
to the wall shear stress phase shift. The LES computations were validated using the DNS
results from Spalart and Baldwin  on oscillating ﬂows showing a good agreement between
the Wall-Function and Wall-Normal Resolved approaches. When the two near wall modeling
were tested in pulsating conditions discrepancies have been found in the wall-shear stress phase
shift predictions. Similar results have been obtained also for URANS calculations using as a test
case the experimental results from Tardu et al.. It is interesting at this stage to point out that
10 Turbulence, Heat and Mass Transfer 5
the oscillations in the case of Spalart and Baldwin can be considered to be in the quasi-steady
regime. The ls parameter in this case has been computed using the maximum skin friction
velocity during the period of oscillation. We do not expect strong phase shift effects in cases
with such a large values of ls and for this reason the validation results are in accordance with
our ﬁndings. The Wall-Normal Resolved computations and the Low-Reynolds turbulent model
seem to capture the unsteady effects of pulsation on the wall-shear stress phase shift. The use of
wall functions is accurate only in cases in which the oscillations are well above the quasi-steady
regime. In all the other cases the use of wall-functions in URANS and LES computations is
questionable especially in applications for which phase lags can play an important role such as
the prediction of thermo-acoustic instabilities.
This work was carried out within the Marie Curie Research Training Network FLUISTCOM (
Fluid-Structure Interaction for Combustion Systems ) under FP6, European Commission, DG
 H. G. Weller, G. Tabor, H. Jasak and C. Fureby. A Tensorial Approach to Computational
Continuum Mechanics using Object Orientated Techniques. In Computers in Physics.,
 V. Moureau, G. Lartigue, Y. Sommerer, C. Angelberger, O. Colin and T. Poinsot. Numer-
ical Methods for Unsteady Compressible Multi-Component Reacting Flows on Fixed and
Moving Grids. In Journal of Computational Physics, 202-2:710-736, 2005.
 P. R. Spalart and B. S. Baldwin. Direct Simulation of a Turbulent Oscillating Boundary
Layer. In Turbulent Shear Flows 6, pp. 417-440, Springer-Verlag Berlin Heidelberg, 1989.
 S. F. Tardu, G. Binder and R. F. Blackwelder. Turbulent Channel Flow with Large-
Amplitude Velocity Oscillations. In J. Fluid Mech., 267:109-151, 1994.
 J. Smagorinsky. General Circulation Experiments with the Primitive Equations, I. The
Basic Experiment. In Monthly Weather Review, 91-3:99-164, 1963.
 F. Ducros, F. Nicoud and T. Poinsot. Wall-Adapting Local Eddy-Viscosity Models for
Simulations in Complex Geometries. In Numerical Methods for Fluid Dynamics VI,
(Edited by M. J. Baines , pp. 293-299, ICFD, 1998.
 R. I. Issa. Solution of the Implicitly Discretised Fluid ﬂow Equations by Operator-
Splitting. In Journal of Computational physics, 62:40-65, 1985.
 P. Scmitt. Simulation aux grandes echelles de la combustion etag´ e dans les turbines a
´ ´ e `
gaz et son interaction stabilit´ - polluants-thermique - TH/CFD/05/45 Institut National
Polytechnique de Toulouse PhD Thesis, 2005.
 B. E. Launder and B. I. Sharma. Application of the Energy Dissipation Model of Tur-
bulence to the Calculation of Flow Near a Spinning Disc. In Letters in Heat and Mass
Transfer, 1:131-138, 1974.
 D. Panara, M. Porta, T. Schoenfeld. LES and URANS Unsteady Boundary Layer Strate-
gies for Pulsating and Oscillating Turbulent Channel Flow Applications. In Proceedings
ECCOMAS CFD , TU Delft, The Netherlands, 2006.
- Related pdf books
- DRAFT: LARGE-EDDY SIMULATION AND CONJUGATE HEAT TRANSFER ...
- Format des données SAFRAN et scénarios climatiques désagrégés ...
- Approximating Hessians in multilevel unconstrained optimization
- Institute of Aeronautics and Astronautics Inc All rights reserved
- Towards Large Eddy Simulation in InternalCombustion Engines
- The Augmented Block Cimmino Distributed method
- Analyse th¶eoriqueet simulationnum¶erique
- As shown by Sausen, the level of scientific understanding of the ...
- A two-step chemical scheme for kerosene-air
- Mixed-precision preconditioners in parallel domain ...
- EDF Electricité de France Fonctionnement et modélisation
- Numerical Analysis of the Tip Leakage Flow Field in a Transonic
- CARACTERISATION ACOUSTIQUE D’UNE PLAQUE´ MULTIPERFOREE´
- Incomplete factorization based on a HID Hybrid Solver
- Using Numerical Values in Sparse Matrix Reorderings for ARMS ...
- Thermoacoustic stability of a helicopter gas turbine combustor