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

NEAR-WALL MODEL FOR LARGE EDDY SIMULATION THATINCORPORATES RESOLVED NON-EQUILIBRIUM ANDREYNOLDS-NUMBER EFFECTS

Kevin Griffin
Center for Turbulence Research
Stanford University
481 Panama Mall, Stanford, California 94305 USA
kevinpg@stanford.edu

Lin FuDepartment of Mechanical and Aerospace Engineering
Department of Mathematics
Center for Ocean Research in Hong Kong and Macau (CORE)
The Hong Kong University of Science and Technology
Clear Water Bay, Kowloon, Hong Kong
linfu@ust.hk

ABSTRACT

this ODE becomes

In wall-modeled large-eddy simulations (WMLES), thenear-wall model plays a significant role in predicting the skinfriction, although the majority of the boundary layer is resolved by the outer large-eddy simulation (LES) solver. In
this work, we aim at developing a new ordinary differential equation (ODE)-based wall model, which is as simple as
the classical equilibrium model yet capable of capturing nonequilibrium effects and low Reynolds number effects. The
proposed model reformulates the classical equilibrium modelby parameterizing the mixing-length function in terms of theboundary-layer shape factor. The performance of the newmodel is validated by predicting a wide range of canonicalflows with the friction Reynolds number between 200 and
5200, and the Clauser pressure-gradient parameter between 0.3 and 4. Compared to the classical equilibrium wall model,
remarkable error reduction in terms of the skin friction prediction is obtained by the new model. Moreover, since the
new model is ODE-based, it is straightforward to be deployedfor predicting flows with complex geometries and therefore
promising for a wide range of applications.

τ+
dU +
=
,
dy+
1 + ℓ+(1)

The mixing length can be redefined as ℓ+n to account for thewall-normal variation of the total shear stress, which leads to
1
dU +
=
,
dy+
1 + ℓ+n(2)

The mixing length is parameterized as



y+ 2
+
ℓ+
=
κy
1
−
exp
−(
)
n
A+

(3)

Cabot & Moin (2000) chose the damping coefficient A+ = 17.
For near-equilibrium flows, Cabot’s model is quite successful since the log law coefficients κ and B are relatively
robust. However, for flows where non-equilibrium effects
are significant, such as the adverse-pressure-gradient boundary layers shown in Figure 1(a), the log intercept constant B
exhibits non-universality. The log slope also exhibits nonuniversality but to a lesser extent. Therefore, it is assumed
to be a constant in this work.
In this work, we permit the damping function to dependon the interior LES solution in order to capture the variationin the log intercept in non-equilibirium flows. The motivation
for this decision is that the interior partial differential equation(PDE) solver can capture non-equilibrium effects that are neglected by the ODE wall model. The exact dependence that is
advanced is informed by data analysis and is discussed next.

MODEL
We begin by considering the definition of the total shear
stress τ/ρ = (ν + νt )(∂U/∂ y), where ν is the kinematic viscosity, ∂U/∂ y is the mean shear, and νt is the eddy viscosity,
which Cabot (1995) defines as νt = ℓuτ , where ℓ is thepmixing
length (specified later), the wall friction velocity uτ = τw /ρ,the wall shear stress τw = τ|y=0 , and ρ denotes the density.
This definition of the total shear stress can be interpreted as an
ordinary differential equation (ODE) for the mean streamwise
velocity profile U[y]. After non-dimensionalizing by uτ and ν,1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
H = δ ∗ /θ , where the displacement thickness δ ∗ is defined as

Dependence of wall model closure on the LES
A database is assembled from various existing highfidelity simulations (either Direct Numerical Simulations,
DNSs, or Wall-Resolved Large-Eddy Simulations, WRLESs).
Included are zero-pressure-gradient boundary layers (ZPGBLs) (Sillero et al., 2013; Spalart, 1988; Eitel-Amor et al.,2014), fully developed channel and pipe flows (Lee & Moser,
2015; Wu & Moin, 2008), adverse-pressure-gradient boundary layers (APGBLs) with five different streamwise-varying
pressure-gradient conditions (Bobke et al., 2017), and two airfoil flows with specified angle of attack AoA= 0◦ , 5◦ respectively (Tanarro et al., 2020; Vinuesa et al., 2018).
For each profile in this database, we solve the ODE given
in Eq. 2 for various values of the damping function A+ to determine the best fit value of A+ , which is defined as the valueof A+ that minimizes the least square error between the reference data and the modeled velocity profile (the solution of theODE) over the domain y ∈ [0, 0.1δ ], where δ is a measure of
the boundary layer thickness.
The best fit values of the damping function A+ are plotted+versus the pressure gradient parameter α = ∂∂ Px+ in Figure 1(a)and versus the boundary layer shape factor H in Figure 2(b).
Note that H = δ ∗ /θ , where δ ∗ denotes the boundary layer
displacement thickness and θ denotes the boundary layer momentum thickness. Clearly, there is a much higher correlationof A+ with H than with α. This motivates the following modeling choice
A+ [H] = C0 +C1 H

δ∗ =

Z δ

1−

0

U
Ue


dy.

(5)

and the momentum thickness θ is defined as
Z δ

θ=
0

U
Ue



U
1−
dy.
Ue

(6)

The rationale for the hypothesis that A = A[H] is based on thefollowing two observations. First, non-equilibrium effects directly modify the the boundary-layer shape factor H Tamaki
et al. (2020). Second, the dominant contribution to the shapefactor in WMLES comes from the outer PDE solver, which
captures non-equilibrium (history) effects by construction, as
demonstrated in Griffin & Fu (2020). These observations imply that by correlating the inner wall model with the shapefactor, the history effects captured by the outer solver can beleveraged by the wall model.
As show in Figure 2, a much better collapse is observed
when plotted versus the shape factor H, than versus the innerand outer pressure-gradient parameter (with a coefficient ofdetermination R2 = 0.90 compared to 0.73).
Note that these regressions only apply to the present mixing length model and damping function given in Eq. (3), forfully turbulent, incompressible flows with zero wall penetration. Although the optimal choices of the regression coefficients may be different for other mixing-length models, thesuitability of correlating A+ with H and Reτ may still hold ingeneral.

(4)

where C0 = 27 and C1 = 7.4 (according to the least squres
regression of the data in Figure 2(b)). Fitting of B or A+ versus
α was considered by Huffman & Bradshaw (1972); Granville(1989); Johnstone et al. (2010); Nickels (2004); Duprat et al.
(2011), but fitting versus H has not been considered before andis better motivated by the present database.
In equilibrium flows, such as the channel, pipe, and zeropressure gradient boundary layer, the law of the wall has
needed only small modifications to optimally fit these velocity
profiles. However, in pressure-gradient boundary layers, relatively large changes in the logarithmic intercept are required(see Figure 1(a)). Most of the classical mixing-length-basedmodels for non-equilibrium flows rely on a non-dimensionalpressure-gradient parameter as an input.
As discussed above, the classical pressure-gradient parameters α or β in combination with Reynolds number Reτ do
not uniquely and completely define the boundary layer velocity profile. The reason is that these pressure-gradient parameters are unaware of the spatial (or temporal in a Lagrangian
sense) history of the flow, see e.g., Johnstone et al. (2010)(concerning α) and Bobke et al. (2017) (concerning Reτ andβ ). Only the local pressure gradient effect is taken into accountinstead of the integrated effect on the flow in the streamwise(or temporal) dimension.
However, these integrated history effects are significant in
non-equilibrium flows. The most straightforward solution forincorporating the boundary layer history effects into the wallmodel is to employ a PDE-based wall model, but this leads
to a significantly increased computational cost. Instead, we
propose that the same objective can be achieved by correlating A+ with the boundary-layer shape factor H for the ODEbased wall model and hypothesizing that A+ (H) is a universal
function. For boundary layers, the shape factor is defined as

Prediction the wall shear stress
In this a priori study, the data at the matching location(y = 0.1δ ) is provided from DNS or WRLES, such that any
resulting errors can be attributed to the wall model instead ofto the matching data. The relative error ετw is defined as thedifference between the wall shear stress τw computed from thewall model and that from DNS or WRLES.
As shown in Figure 3(a), the relative error of the wall
shear stress from the classical Cabot’s model is as large as
17% for cases with strong pressure gradients. Meanwhile,the error from the new model, as shown in Figure 3(b), is
typically less than 2%, with a maximum of 5%. Cases withthe strongest pressure gradients have the largest errors for theclassical model and the most remarkable error reductions byswitching to the new model. For very few cases there is a tiny
error increase of about 1%, which can be attributed to the fitting errors evident in Figure 2(b).
The resulting velocity profile from the classical andpresent models are shown in Figure 1(b). The classical model
incorrectly predicts the equilibrium solution in this flow. On
the other hand, the new model correctly predicts a shift in thelog law that is in agreement with the reference WRLES data.

CONCLUSION
Whereas most classical stress-based wall models assume
a universal value for the mixing-length damping coefficient
A+ , the new method correlates A+ with the boundary-layer
shape factor H. The proposed correlation of A+ [H] makes asubstantial improvement to the prediction of the velocity profile and the wall shear stress for a large range of Reynolds numbers and pressure gradient conditions. By depending on H, the212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Tamaki, Yoshiharu, Fukushima, Yuma, Kuya, Yuichi &Kawai, Soshi 2020 Physics and modeling of trailing-edge
stall phenomena for wall-modeled large-eddy simulation.
Phys. Rev. Fluids 5 (074602).
Tanarro, Á., Vinuesa, Ricardo & Schlatter, P. 2020 Effect ofadverse pressure gradients on turbulent wing boundary layers. J. Fluid Mech. 883 (A8).
Vinuesa, R., Negi, P. S., Atzori, M., Hanifi, A., Henningson,
D. S. & Schlatter, P. 2018 Turbulent boundary layers aroundwing sections up to Rec =1,000,000. Int. J. Heat Fluid Flow
72, 86–99.
Wu, Xiaohua & Moin, Parviz 2008 A direct numerical simulation study on the mean velocity characteristics in turbulent
pipe flow. J. Fluid Mech. 608, 81–112.

inner ODE-based wall model can leverage additional information from the outer PDE-based solver. As a result, the newmodel incorporates an integral measure of the streamwise andtemporal history of the flow and is accurate in non-equilibrium
scenarios, while retaining similar computational efficiency asclassical equilibrium models.

ACKNOWLEDGEMENTS
Kevin Griffin acknowledges support from the National Defense Science and Engineering Graduate Fellowship, the Stanford Graduate Fellowship, and NASA grant No.
NNX15AU93A. Lin Fu acknowledges the funding from the
Guangdong Basic and Applied Basic Research Foundation(No. 2022A1515011779) and the funding from CORE as a
joint research center for ocean research between QNLM andHKUST. We wish to acknowledge helpful feedback from P.
Moin.

REFERENCES
Bobke, A., Vinuesa, Ricardo, Örlü, Ramis & Schlatter, P. 2017
History effects and near equilibrium in adverse-pressuregradient turbulent boundary layers. J. Fluid Mech. 820,
667–692.
Cabot, William H. 1995 Large-eddy simulations with wallmodels. Cent. Turbul. Res. Annu. Res. Briefs pp. 41–50.
Cabot, William H. & Moin, Parviz 2000 Approximate wallboundary conditions in the large-eddy simulation of highReynolds number flow. Flow, Turbul. Combust. 63, 269–
291.
Duprat, C., Balarac, G., Métais, O., Congedo, P. M. &Brugière, O. 2011 A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient. Phys.
Fluids 23 (1).
Eitel-Amor, Georg, Örlü, Ramis & Schlatter, Philipp 2014
Simulation and validation of a spatially evolving turbulent
boundary layer up to Reθ = 8300. Int. J. Heat Fluid Flow
47, 57–69.
Granville, P. S. 1989 A modified van driest formula for the
mixing length of turbulent boundary layers in pressure gradients. J. Fluids Eng. Trans. ASME 111 (1), 94–97.
Griffin, Kevin Patrick & Fu, Lin 2020 A new ODE-based turbulence wall model accounting for pressure gradient andReynolds number effects .
Huffman, G. David & Bradshaw, Peter 1972 A note on vonKármán’s constant in low Reynolds number turbulent flows.
J. Fluid Mech. 53 (1), 45–60.
Johnstone, Roderick, Coleman, Gary N. & Spalart, Philippe R.
2010 The resilience of the logarithmic law to pressure gradients: evidence from direct numerical simulation. J. Fluid
Mech. 643, 163–175.
Lee, Myoungkyu & Moser, Robert D. 2015 Direct numericalsimulation of turbulent channel flow up to Reτ = 5200. J.
Fluid Mech. 774, 395–415.
Nickels, T. B. 2004 Inner scaling for wall-bounded flows subject to large pressure gradients. J. Fluid Mech. 521, 217–
239.
Sillero, Juan A., Jiménez, Javier & Moser, Robert D. 2013
One-point statistics for turbulent wall-bounded flows at
Reynolds numbers up to δ + =2000. Phys. Fluids 25,
105102.
Spalart, Philippe R. 1988 Direct simulation of a turbulent
boundary layer up to Rθ = 1410. J. Fluid Mech. 187, 61–98.

3

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

Figure 1: (a) The mean streamwise velocity U + plotted versus the wall-normal coordinate y+ for wall-bounded flows withvarious values of the streamwise pressure gradient parameter β . β = 4.2, 2.3, 1.1 correspond to WRLESs of APGBLs by
Bobke et al. (2017), β = 0 is the DNS of a ZPGBL at Reτ = 2000 by Sillero et al. (2013), and β = −0.1 is the DNS ofa channel flow at Reτ = 5200 by (Lee & Moser, 2015). (b) The data from the APGBL with β = 4.2 is reduplicated in(b) (dashed line). Also shown are wall-modeled velocity profiles, which are constructed by solving Cabot’s model (redcircles) and the present model (green diamonds) using matching data from the WRLES sampled at y = 0.1δ . In bothpanels, the log law is shown for reference (κ = 0.41 and B = 5.2).

Figure 2: Distributions of the best-fit value of the wall-model damping coefficient A+ plotted versus (a) the inner pressure-gradient parameter α and (b) the boundary layer shape factor H. The green lines denote the least-squares regressions of
the best-fit coefficients with the corresponding abscissas. The symbol color indicates Reτ of the reference data and thesymbol type indicates the flow type, i.e., channel flows (diamonds), ZPGBLs (circles), pipe flows (squares), APGBLs(triangles), and airfoil flows (pentagrams).

Figure 3: Distributions of the relative error ετw between the wall stress predicted by the well-resolved simulations and thatby the wall-modeled simulations (from Cabot’s model (a) and the present model (b)) versus the underlying shape factor
H. The symbol color indicates Reτ of the reference data and the symbol type indicates the flow type, i.e., channel flows(diamonds), ZPGBLs (circles), pipe flows (squares), APGBLs (triangles), and airfoil flows (pentagrams).

4