A multiscale approach to hybrid RANS/LES wall modeling within a high-order discontinuous Galerkin scheme using function enrichment

TL;DR

This paper introduces a multiscale hybrid RANS/LES wall modeling method using function enrichment within a high-order discontinuous Galerkin scheme, enabling coarse near-wall meshes and improving simulation accuracy and efficiency.

Contribution

It provides a rigorous derivation of a new multiscale turbulence model and a DG discretization that seamlessly combines RANS and LES in a single equation with independent shape functions.

Findings

Achieves grid independence and superior accuracy in separated flows.

Provides a two-order magnitude speed-up over wall-resolved LES.

Effectively prevents log-layer mismatch in boundary layers.

Abstract

We present a novel approach to hybrid RANS/LES wall modeling based on function enrichment, which overcomes the common problem of the RANS-LES transition and enables coarse meshes near the boundary. While the concept of function enrichment as an efficient discretization technique for turbulent boundary layers has been proposed in an earlier article by Krank & Wall (J. Comput. Phys. 316 (2016) 94-116), the contribution of this work is a rigorous derivation of a new multiscale turbulence modeling approach and a corresponding discontinuous Galerkin discretization scheme. In the near-wall area, the Navier-Stokes equations are explicitly solved for an LES and a RANS component in one single equation. This is done by providing the Galerkin method with an independent set of shape functions for each of these two methods; the standard high-order polynomial basis resolves turbulent eddies where the mesh is sufficiently fine and the enrichment automatically computes the ensemble-averaged flow if the LES mesh is too coarse. As a result of the derivation, the RANS model is consistently applied solely to the RANS degrees of freedom, which effectively prevents the typical issue of a log-layer mismatch in attached boundary layers. As the full Navier-Stokes equations are solved in the boundary layer, spatial refinement gradually yields wall-resolved LES with exact boundary conditions. Numerical tests show the outstanding characteristics of the wall model regarding grid independence, superiority compared to equilibrium wall models in separated flows, and achieve a speed-up by two orders of magnitude compared to wall-resolved LES.

Tables

Table 1. Table 1: Overview of simulation cases for the turbulent channel flow: Investigation of several Reynolds numbers (Figure  5 ), comparison of the new approach to delayed-detached-eddy simulation (DDES) including function enrichment 16 (Figure  6 ), grid stretching factors and mesh aspect ratios (Figure  7 ), transition from wall-resolved to wall-modeled LES (Figure  8 ), weighting of the enrichment function with a constant or linear basis (Figure  9 ), and the performance comparing the wall model with wall-resolved LES (Figure  10 ). The number of cells per direction i 𝑖 i is denoted N i ​ e subscript 𝑁 𝑖 𝑒 N_{ie} and the normalized dimensions of the first off-wall cell is defined as Δ ​ y 1 ​ e + = Δ ​ y 1 ​ e ​ u τ / ν Δ superscript subscript 𝑦 1 𝑒 Δ subscript 𝑦 1 𝑒 subscript 𝑢 𝜏 𝜈 \Delta y_{1e}^{+}=\Delta y_{1e}u_{\tau}/\nu in wall-normal direction and Δ ​ { x , z } e + = Δ ​ { x , z } e ​ u τ / ν Δ superscript subscript 𝑥 𝑧 𝑒 Δ subscript 𝑥 𝑧 𝑒 subscript 𝑢 𝜏 𝜈 \Delta\{x,z\}_{e}^{+}=\Delta\{x,z\}_{e}u_{\tau}/\nu in wall-parallel direction. The number of LES grid points per direction i 𝑖 i is N i = ( k + 1 ) ​ N i ​ e subscript 𝑁 𝑖 𝑘 1 subscript 𝑁 𝑖 𝑒 N_{i}=(k+1)N_{ie} with k = 4 𝑘 4 k=4 in all cases.

Fig.	Case	$N_{1e}\times N_{2e}\times N_{3e}$	$R{e}_{\tau }$	$l$	$\gamma$	$\Delta y_{1e}^{+}$	$\max (\Delta x_e^{+},\Delta z_e^{+})$	$\frac{\Delta y_{1e}^{+}\left(k+1\right)}{\max (\Delta x_e^{+},\Delta z_e^{+})}$
5	$ch395\_N12\times 8\times 12\_k4l0\_\gamma 0.8$	$12\times 8\times 12$	$395$	$0$	$0.8$	$76$	$207$	$1.84$
	$ch950\_N12\times 8\times 12\_k4l0\_\gamma 1.6$	$12\times 8\times 12$	$950$	$0$	$1.6$	$91$	$497$	$0.92$
	$ch2000\_N24\times 8\times 24\_k4l0\_\gamma 2.2$	$24\times 8\times 24$	$2,000$	$0$	$2.2$	$96$	$524$	$0.92$
	$ch5200\_N48\times 16\times 48\_k4l0\_\gamma 2.05$	$48\times 16\times 48$	$5,200$	$0$	$2.05$	$114$	$681$	$0.84$
6	$ch950\_N{16}^3\_k4l0\_\gamma 0.001$	$16\times 16\times 16$	$950$	$0$	$0.001$	$119$	$373$	$1.60$
	$ch950\_N{16}^3\_k4l0\_\gamma 0.001\_DDES$16	$16\times 16\times 16$	$950$	$0$	$0.001$	$119$	$373$	-
7	$ch950\_N{16}^3\_k4l0\_\gamma 0.001$	$16\times 16\times 16$	$950$	$0$	$0.001$	$119$	$373$	$1.60$
	$ch950\_N{16}^3\_k4l0\_\gamma 0.8$	$16\times 16\times 16$	$950$	$0$	$0.8$	$85$	$373$	$1.14$
	$ch950\_N{16}^3\_k4l0\_\gamma 1.2$	$16\times 16\times 16$	$950$	$0$	$1.2$	$59$	$373$	$0.79$
	$ch950\_N{16}^3\_k4l0\_\gamma 1.5$	$16\times 16\times 16$	$950$	$0$	$1.5$	$42$	$373$	$0.56$
	$ch950\_N{16}^3\_k4l0\_\gamma 1.0$	$16\times 16\times 16$	$950$	$0$	$1.0$	$72$	$373$	$0.97$
	$ch950\_N16\times 16\times 32\_k4l0\_\gamma 1.0$	$16\times 16\times 32$	$950$	$0$	$1.0$	$72$	$373$	$0.97$
8	$ch950\_N{24}^3\_k4\_\gamma 2.25$	$24\times 24\times 24$	$950$	-	$2.25$	$10$	$249$	$0.20$
	$ch950\_N{24}^3\_k4l0\_\gamma 1.75$	$24\times 24\times 24$	$950$	$0$	$1.75$	$19$	$249$	$0.38$
	$ch950\_N{24}^3\_k4l0\_\gamma 1.5$	$24\times 24\times 24$	$950$	$0$	$1.5$	$27$	$249$	$0.54$
	$ch950\_N{24}^3\_k4l0\_\gamma 1.25$	$24\times 24\times 24$	$950$	$0$	$1.25$	$36$	$249$	$0.72$
	$ch950\_N{24}^3\_k4l0\_\gamma 1.0$	$24\times 24\times 24$	$950$	$0$	$1.0$	$46$	$249$	$0.92$
9	$ch950\_N{16}^3\_k4l0\_\gamma 1.2$	$16\times 16\times 16$	$950$	$0$	$1.2$	$59$	$373$	$0.79$
	$ch950\_N{16}^3\_k4l1\_\gamma 1.2$	$16\times 16\times 16$	$950$	$1$	$1.2$	$59$	$373$	$0.79$
10	$ch950\_N12\times 8\times 12\_k4l0\_\gamma 1.6$	$12\times 8\times 12$	$950$	$0$	$1.6$	$91$	$497$	$0.92$
	$ch950\_N{24}^3\_k4\_\gamma 2.25$	$24\times 24\times 24$	$950$	-	$2.25$	$10$	$249$	$0.20$
	$ch950\_N{32}^3\_k4\_\gamma 2.0$	$32\times 32\times 32$	$950$	-	$2.0$	$10$	$187$	$0.27$

Table 2. Table 2: Simulation cases and resolutions of the periodic hill flow. The cases at R ​ e H = 10 , 595 𝑅 subscript 𝑒 𝐻 10 595 Re_{H}=10{,}595 and 19 , 000 19 000 19{,}000 use 32 × 16 × 16 32 16 16 32{\times}16{\times}16 cells with varying grid stretching in vertical direction and, at R ​ e H = 37 , 000 𝑅 subscript 𝑒 𝐻 37 000 Re_{H}=37{,}000 , a finer mesh with 64 × 32 × 32 64 32 32 64{\times}32{\times}32 grid cells is additionally considered. The polynomial degrees are k = 4 𝑘 4 k=4 and l = 0 𝑙 0 l=0 for all simulation cases, and the number of grid points per direction is k + 1 𝑘 1 k+1 in each cell. The addition N ​ W ​ M 𝑁 𝑊 𝑀 NWM stands for ‘no wall model’. The separation and reattachment lengths x 1 , sep subscript 𝑥 1 sep x_{1,\mathrm{sep}} and x 1 , reatt subscript 𝑥 1 reatt x_{1,\mathrm{reatt}} correspond to the zero-crossings of the skin friction.

Case	$N_{e1}\times N_{e2}\times N_{e3}$	$N_1\times N_2\times N_3$	$R{e}_H$	$\max \left(\Delta y_{1e}^{+}\right)$	$x_{1,\mathrm{sep}}/H$	$x_{1,\mathrm{reatt}}/H$
$ph10595\_stretch1$	$32\times 16\times 16$	$160\times 80\times 80$	$10,595$	$86$	$0.32$	$4.23$
$ph10595\_stretch2$	$32\times 16\times 16$	$160\times 80\times 80$	$10,595$	$73$	$0.23$	$4.24$
$ph10595\_stretch3$	$32\times 16\times 16$	$160\times 80\times 80$	$10,595$	$60$	$0.19$	$3.98$
$ph10595\_stretch2\_NWM$	$32\times 16\times 16$	$160\times 80\times 80$	$10,595$	-	$0.29$	$3.99$
$KKW\_DNS$31	-	$896\times 448\times 448$	$10,595$	-	$0.19$	$4.51$
$ph19000\_stretch1$	$32\times 16\times 16$	$160\times 80\times 80$	$19,000$	$138$	$0.42$	$2.85$
$ph19000\_stretch2$	$32\times 16\times 16$	$160\times 80\times 80$	$19,000$	$113$	$0.35$	$3.60$
$ph19000\_stretch3$	$32\times 16\times 16$	$160\times 80\times 80$	$19,000$	$95$	$0.25$	$3.64$
$ph19000\_stretch2\_NWM$	$32\times 16\times 16$	$160\times 80\times 80$	$19,000$	-	$0.32$	$1.99$
$RM\_EXP$55	-	-	$19,000$	-	-	$3.94$
$ph37000\_stretch3$	$32\times 16\times 16$	$160\times 80\times 80$	$37,000$	$156$	$0.37$	$2.78$
$ph37000\_stretch3\_fine$	$64\times 32\times 32$	$320\times 160\times 160$	$37,000$	$78$	$0.28$	$3.38$
$ph37000\_stretch3\_NWM$	$32\times 16\times 16$	$160\times 80\times 80$	$37,000$	-	-	-
$RM\_EXP$55	-	-	$37,000$	-	-	$3.76$
$CM\_WMLES$8	-	$128\times 64\times 64$	$37,000$	-	-	$2.3$
$CM\_WMLES\_fine$8	-	$256\times 128\times 128$	$37,000$	-	-	$2.8$

Table 3. Table 3: Location of support points y i + superscript subscript 𝑦 𝑖 y_{i}^{+} and corresponding values ψ i subscript 𝜓 𝑖 \psi_{i} for the efficient evaluation of the enrichment function.

$i$	$y_i^{+}$	${\psi }_i$
$0$	$0.0$	$0.0$
$1$	$5.0$	$4.88298776233176$
$2$	$11.0$	$8.91824406645381$
$3$	$24.0$	$12.3978516813118$
$4$	$59.0$	$15.1875389926298$
$5$	$144.0$	$17.4177125619900$
$6$	$361.0$	$19.6484300823042$
$7$	$946.0$	$21.9930107788854$
$8$	$2517.0$	$24.3778307011372$

Equations

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial \mathbf{u}}{\partial t}+$$

$${\mathcal{F}}^c(\mathbf{u})=\mathbf{u}\otimes \mathbf{u},{\mathcal{F}}^{\nu }(\mathbf{u})=2\nu \mathbf{ϵ}(\mathbf{u}),$$

$$\mathbf{u}={\mathbf{u}}_0\text{ in }\Omega \times \{0\}.$$

$$\mathbf{u}={\mathbf{g}}_{\mathbf{u}}=0\text{ on }\partial {\Omega }^D\times [0,\mathcal{T}].$$

$${\mathbf{u}}_H$$

$$p_H$$

$$a_1+a_2$$

$$⟨{\mathbf{u}}_H⟩$$

$$\nabla \cdot {\mathbf{u}}_H$$

$$\frac{\partial {\mathbf{u}}_H}{\partial t}+$$

$${\mathbf{u}}_H$$

$$p_H$$

$$a_1+a_2+a_3$$

$$⟨{\mathbf{u}}_H⟩$$

$${\mathbf{\tau }}_H=$$

$$+a_1{a}_2\left(\bar{\mathbf{u}}-⟨\mathbf{u}⟩\right)\otimes \left(\bar{\mathbf{u}}-⟨\mathbf{u}⟩\right)+a_1{a}_3\left(\bar{\mathbf{u}}-⟨\bar{\mathbf{u}}⟩\right)\otimes \left(\bar{\mathbf{u}}-⟨\bar{\mathbf{u}}⟩\right)+a_2{a}_3\left(⟨\mathbf{u}⟩-⟨\bar{\mathbf{u}}⟩\right)\otimes \left(⟨\mathbf{u}⟩-⟨\bar{\mathbf{u}}⟩\right),$$

$${\mathbf{\tau }}_{\mathrm{LES}}$$

$${\mathbf{\tau }}_{\mathrm{RANS}}$$

$${\mathbf{\tau }}_3$$

$${\mathbf{\tau }}_H=$$

$$=$$

$${\mathbf{\tau }}_H\approx ⟨{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }⟩-⟨{\bar{\mathbf{u}}}^{\prime }\otimes {\bar{\mathbf{u}}}^{\prime }⟩.$$

$${\mathbf{\tau }}_H$$

$$\approx -2{\nu }_t\left(\mathbf{ϵ}\left(⟨\mathbf{u}⟩\right)-\mathbf{ϵ}\left(⟨\bar{\mathbf{u}}⟩\right)\right)$$

$$=-2{\nu }_t\mathbf{ϵ}\left(\tilde{\mathbf{u}}\right)$$

$${\nu }_t=l_{\mathrm{mix}}^2|\mathbf{ϵ}\left(\tilde{\mathbf{u}}\right)|$$

$$l_{\mathrm{mix}}=\kappa y_{\mathrm{DES}}\left(1-\exp \left(-y^{+}/A^{+}\right)\right),$$

$$y_{\mathrm{DES}}=\min \left(y,C_{\mathrm{DES}}h\right)\approx {\left(\frac{1}{y^{\beta }}+\frac{1}{{\left(C_{\mathrm{DES}}h\right)}^{\beta }}\right)}^{-1/\beta },$$

$$\nabla \cdot \left(\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)$$

$$\frac{\partial \left(\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)}{\partial t}+$$

$${\mathcal{F}}^{{\nu }_t}\left(\tilde{\mathbf{u}}\right)\cdot \mathbf{n}=0\text{ on }\partial {\tilde{\Omega }}^I\times [0,\mathcal{T}].$$

$$\mathcal{C}\left(q,\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)=$$

$$\mathcal{M}\left(\mathbf{v},\frac{\partial \left(\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)}{\partial t}\right)+{\mathcal{F}}^c\left(\mathbf{v},\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)+\mathcal{P}\left(\mathbf{v},p_H\right)-{\mathcal{F}}^{\nu }\left(\mathbf{v},\bar{\mathbf{u}}\right)-{\mathcal{F}}^{\nu +{\nu }_t}\left(\mathbf{v},\tilde{\mathbf{u}}\right)=$$

$$\mathcal{C}\left(q,\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)=$$

$$\mathcal{M}\left(\bar{\mathbf{v}},\frac{\partial \left(\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)}{\partial t}\right)+{\mathcal{F}}^c\left(\bar{\mathbf{v}},\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)+\mathcal{P}\left(\bar{\mathbf{v}},p_H\right)-{\mathcal{F}}^{\nu }\left(\bar{\mathbf{v}},\bar{\mathbf{u}}\right)-{\mathcal{F}}^{\nu +{\nu }_t}\left(\bar{\mathbf{v}},\tilde{\mathbf{u}}\right)=$$

$$\mathcal{M}\left(\tilde{\mathbf{v}},\frac{\partial \left(\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)}{\partial t}\right)+{\mathcal{F}}^c\left(\tilde{\mathbf{v}},\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)+\mathcal{P}\left(\tilde{\mathbf{v}},p_H\right)-{\mathcal{F}}^{\nu }\left(\tilde{\mathbf{v}},\bar{\mathbf{u}}\right)-{\mathcal{F}}^{\nu +{\nu }_t}\left(\tilde{\mathbf{v}},\tilde{\mathbf{u}}\right)=$$

$$\mathcal{M}\left(\tilde{\mathbf{v}},\frac{\partial \left(\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)}{\partial t}\right)+{\mathcal{F}}^c\left(\tilde{\mathbf{v}},\bar{\mathbf{u}}+\tilde{\mathbf{u}}\right)+\mathcal{P}\left(\tilde{\mathbf{v}},p_H\right)-{\mathcal{F}}^{\nu +{\nu }_t}\left(\tilde{\mathbf{v}},\tilde{\mathbf{u}}\right)=l\left(\tilde{\mathbf{v}}\right)$$

Full text

\jnlcitation\cname

, , and , \ctitleA multiscale approach to hybrid RANS/LES wall modeling within a high-order discontinuous Galerkin scheme using function enrichment, \cjournalInt. J. Numer. Meth. Fluids, \cvol2017;00:1–32.

††articletype: Research Article

1 Introduction

The simulation of turbulent boundary layers poses two fundamental challenges: (i) the mean velocity gradient in the viscous sublayer can become very high in engineering applications; this gradient has to be resolved or modeled for the prediction of the skin friction, and (ii) eddies of various spatial and temporal scales have to be resolved or modeled for the computation of the mean velocity profile (and thus the skin friction). Classical hybrid RANS/LES wall modeling is popular in research and industry as it promises a vast reduction in computational cost compared to wall-resolved LES through modeling of the near-wall turbulence (cf. (ii)), while being potentially much more accurate than RANS in separated flows. Yet the velocity gradient is usually resolved by the scheme (cf. (i)). A complete review of wall modeling approaches and RANS/LES methods would be beyond the scope of this article, so we refer to a series of review articles in this field1, 2, 3, 4, 5. Within the high-order DG method, detached-eddy simulation models16, 6, a hybrid RANS/ILES model7, and wall-stress models based on equilibrium wall functions8, 9, 10 have been developed.

The major challenge in hybrid RANS/LES methods is the transition region from RANS to LES11. The problem is that “LES content has to be generated in the outer part of the boundary layer”4 since the momentum transfer is fully modeled in RANS but relies on the turbulent structures in LES. If the model parameters are carefully tuned within the transition region, the results may be quite convincing12; an error in the sum of the resolved and modeled Reynolds stress would otherwise produce a so-called log-layer mismatch in equilibrium boundary layers11. However, existing zonal RANS/LES models often rely heavily on empiricism, and grid-dependence can be an issue. Enhancements, including artificial forcing methods or synthetic turbulence, may improve the results obtained for attached boundary layers4.

1.1 A multiscale approach to wall modeling via function enrichment

In this work, we propose a novel approach to hybrid RANS/LES wall modeling, which resolves the velocity gradient at the wall with coarser meshes (cf. (i)), models the near-wall turbulence (cf. (ii)), and additionally solves the problem of the RANS–LES transition. The basic idea is that, in the vicinity of the wall, the solution is decomposed into a RANS and an LES component according to the ancient strategy divide et impera. The two components overlap in this near-wall layer, and together they represent the velocity solution. The decoupling is achieved by applying the RANS eddy viscosity on the RANS degrees of freedom only and allows the energy-carrying turbulent structures to fully develop inside the near-wall layer, which entirely circumvents the aforementioned problem of transition at the RANS/LES interface. This structure of the eddy viscosity term is a direct result of applying a new three-level hybrid RANS/LES filter, which is defined based on Germano’s framework of additive filters13. Despite the fact that the two components are expressed and modeled separately, the full Navier–Stokes equations are computed in the whole boundary layer in a single equation, with exact boundary conditions. This ensures robustness in challenging flow conditions such as separated flows featuring a high adverse pressure gradient. The velocity gradient is resolved, although the first off-wall node for the LES may be placed in a range $y_{1}^{+}\sim 1\ldots 100$, which is possible due to the interaction of the RANS and LES scales and due to the use of a tailored spatial discretization via function enrichment.

Function enrichment in the context of wall modeling was recently proposed by Krank and Wall for the standard (continuous) FEM14 and for the discontinuous Galerkin method (DG)16, 15; see also the monograph 17 for further details. The concept of function enrichment utilizes the basic characteristic of the Galerkin method: The user chooses a set of shape functions, and the method automatically finds the best match among those available shape functions. For the LES solution, we take classical polynomial shape functions. However, the RANS solution is based on a discrete solution space tailored for this specific application. The key idea is to construct the shape functions by using a wall-law as an enrichment function, which enables the resolution of the high velocity gradient at the no-slip boundary with very few degrees of freedom. The general principle of function enrichment is universal and has been successfully employed in a broad spectrum of other computational disciplines18, 19, 20.

Once the solution is decomposed into a RANS and an LES component, each with different shape functions, we employ a second beneficial feature of the Galerkin method: The weighting functions have the same structure as the solution functions, allowing an explicit decomposition of the governing equations into separate equations for the RANS and the LES components – although the problem is still solved in one single equation – according to the variational multiscale methodology21, 22. The scale separation allows appropriate turbulence modeling for each of the scales, a RANS turbulence model for the RANS equation and an LES turbulence model for the LES scale, exactly as they are derived through the three-level hybrid RANS/LES filter. We note that this is a fundamental difference to the application of function enrichment in detached-eddy simulation (DES) 16, where the domain is clearly separated into the inner RANS region and the outer LES region, with a ‘grey area’ in between, since the RANS eddy viscosity acts on the polynomial velocity component as well in that study; these differences are detailed in the subsequent section.

1.2 Comparison to other hybrid RANS/LES methodologies

This composition of the solution relates the present method to the nonlinear disturbance equations (NLDE)23, 24, 25, which have been proposed in order to reconstruct the energy-carrying turbulent structures (LES) around an averaged mean flow (RANS), with the primary application in the field of aeroacoustics. In the NLDE, the RANS and LES components overlap in a part of the domain, as they do in the present method. Several features of our approach distinguish between the methods:

• •

In NLDE, the RANS equations are first solved for the mean flow and the turbulent fluctuations are reconstructed in a subsequent step, whereas the RANS and LES solutions are computed simultaneously in a single equation in the present model.

• •

In NLDE, the LES represents fluctuations around a mean flow, while the LES and RANS are each a variable fraction of one solution in the present methodology.

• •

The NLDE usually require the LES to be wall-resolved, whereas in our model the RANS provides the LES with the necessary dissipation if the LES is underresolved, ie, not all energetic scales have to be computed in the LES.

The multiscale wall model is further compared to two of the classical hybrid RANS/LES wall models, the original DES approach by Spalart et al26, which was recently applied with wall modeling via function enrichment in Krank et al16 (in the delayed DES27 (DDES) variant), and the method of blending subgrid models through a weighted sum by Baggett28. Both models follow the idea of two spatially separated zones, a RANS layer and an LES bulk flow, with a transition region in between. They compare to the proposed method as follows:

• •

DES achieves the blending by limiting the wall distance function in the Spalart–Allmaras model by a characteristic grid length scale such that the RANS model acts as a one-equation LES model in the bulk flow, see, eg6 for an application with DG. The present method applies the idea of limiting the wall distance of the RANS model by a characteristic cell length scale in order to account for the resolved Reynolds stresses.

• •

The method of weighted sum blends the subgrid stress tensor from a RANS model at the wall into an LES model in the bulk flow via an explicitly defined blending function. This approach has been refined by Germano13 and Rajamani and Kim29, who suggested to blend LES and RANS filtering operators instead of subgrid models, and it is this filter that will be extended to a three-level filter in the present work.

• •

The major difference of the proposed approach to both the DES and the weighted sum method is, however, that a marginally/underresolved LES, which is not directly affected by the RANS eddy viscosity, extends to the no-slip boundary. It is this difference that avoids nonphysical velocity fluctuations, such as described by Baggett28, within the present method. A direct qualitative and quantitative comparison of DES using function enrichment16 and the present multiscale wall model will be presented in Section 5.1.

1.3 Outline

The wall model has been implemented in a high-performance discontinuous Galerkin incompressible flow solver, which was previously applied for computation of DNS and LES30, 31, DES16, and RANS15. The wall model may be implemented in a continuous Galerkin (standard FEM) framework as well, but the discontinuous Galerkin method provides the simplest formulation for enrichments including adaptation. The flow solver is implemented using fast state-of-the-art matrix-free algorithms.

The remainder of this paper is organized as follows. In Section 2, we derive consistent governing equations for the RANS and LES components in one equation and show how these equations may be solved in the framework of the variational multiscale method. In Section 3, the concept of function enrichment for turbulent boundary layers is revisited. Section 4 gives details on the solution procedure and variational forms of the present discontinuous Galerkin methodology. In Section 5, numerical examples demonstrate the excellent mesh independence of the present wall model and show how the full consistency of the model enables accurate results under separated flow conditions. Conclusions close the article in Section 6.

2 A multiscale approach to wall modeling

The primary innovation of this work is the derivation of a consistent multiscale RANS/LES wall modeling framework for application to wall modeling via function enrichment. Wall modeling via function enrichment has been introduced by Krank and Wall14 and represents a highly efficient discretization method problem-tailored for turbulent boundary layers. The topic of near-wall subgrid modeling in the hybrid RANS/LES context has recently been addressed by extending the wall-function-enriched RANS solver based on the Spalart–Allmaras model15 to DES26. While that work allows the use of comparably coarse grids in wall-normal direction in the inner layer as compared to classical DES, the issue of the RANS–LES transition persisted. However, the Galerkin method in conjunction with function enrichment enables a much more general turbulence model, which avoids the common problem of the RANS–LES transition entirely; the derivation of this model is the topic of this section.

We derive a consistent multiscale framework for application to wall modeling via function enrichment and define a suitable subgrid model. Further essential ingredients to the wall modeling approach are the discrete function spaces including function enrichment, which are recited in Section 3, as well as crucial modifications to the discontinuous Galerkin formulation of the spatial discretization, which are presented in Section 4.

In the following, the solution of the incompressible Navier–Stokes equations in a layer near the wall is split into two parts, an eddy-resolving component (LES) and an ensemble-averaged component (RANS). The equations are derived by applying a new three-level filter to the Navier–Stokes equations, which is defined based on Germano’s framework for additive filtering. The implications of this filter are that both the RANS eddy viscosity term and the definition of the eddy viscosity variable depend only on the RANS component, thus the LES component can evolve independently and is limited solely by its own resolution. Further, the RANS and the LES components are not spatially separated, but overlap in the near-wall layer.

In this section, we begin with a presentation of the governing equations in Section 2.1. Subsequently, the scale separation is introduced by a three-level filter and the differences to the classical hybrid RANS/LES filter are discussed in Section 2.2. Finally, Section 2.3 shows how the variational multiscale method is employed to solve for the RANS and LES solutions using weak forms and appropriate solution as well as weighting functions.

2.1 Incompressible Navier–Stokes Equations

We consider the incompressible Navier–Stokes equations in conservation form:

[TABLE]

with the velocity vector $\bm{u}=\left(u_{1},u_{2},u_{3}\right)^{T}$, the pressure $p$, the body force $\bm{f}=\left(f_{1},f_{2},f_{3}\right)^{T}$, the three-dimensional spatial domain $\Omega$, and the simulation time $\mathcal{T}$. The convective and viscous fluxes are:

[TABLE]

where the kinematic viscosity is denoted $\nu$ and the rate-of-deformation tensor $\bm{\epsilon}=1/2\left(\nabla\bm{u}+(\nabla\bm{u})^{T}\right)$.

At $t=0$, an initial velocity field $\bm{u}_{0}$ is specified as

[TABLE]

Periodic as well as no-slip Dirichlet boundary conditions on solid walls $\partial\Omega^{D}=\partial\Omega$ close the problem with

[TABLE]

2.2 Three-level hybrid RANS/LES filter

We propose a novel hybrid RANS/LES filter based on Germano’s framework of additive filtering13. The filter is constructed such that an underresolved LES solution overlaps with the RANS solution in the near-wall area. In this subsection, we first recite the classical two-level additive hybrid RANS/LES filter13 and then use that framework to define a new three-level filter. In applying the new filter to the incompressible Navier–Stokes equations, we obtain a set of equations for an eddy-resolving velocity variable and a statistical velocity variable, and their sum represents the final velocity solution. These equations have formally seven unknowns, three LES velocity components, three RANS velocity components, and one pressure variable. In practice, the LES and RANS components will be solved by taking a different set of shape functions tailored for each of the two scales (see Section 2.3 and 3).

The derivations address the near-wall area of a computational domain. Figure 1 gives an overview of the general concept and the variables defining the wall model. The near-wall region is denoted $\widetilde{\Omega}\subset\Omega$; the baseline LES solver is used outside this area. Coupling conditions defining the transition from the near-wall region to the outer region at the interface $\partial\widetilde{\Omega}^{I}$ are also discussed in this section.

2.2.1 The classical two-level additive hybrid RANS/LES filter

Germano’s hybrid RANS/LES filter13 defines the velocity and pressure as a weighted sum of the LES-filtered and RANS-averaged quantities, yielding

[TABLE]

Herein, $\bar{(\cdotp)}$ denotes an LES-filtering operator and $\langle\cdotp\rangle$ a statistical (ensemble-averaging) operator. As an LES filter, implicit grid filtering is commonly considered, which may also be interpreted as a variational projection in the context of the Galerkin method32. The continuous blending functions $0\leq a_{1}\leq 1$ and $a_{2}=1-a_{1}$ may be spatially dependent and are predefined. The filter fulfills two basic requirements, namely

[TABLE]

and equivalently for the pressure. Applying this filter to the governing equations (1)–(2) and using the rule of permutation of the derivative and the filtering operator13, we obtain

[TABLE]

where $\bm{\tau}_{H}$ is the corresponding subgrid stress tensor, which is given by Germano13. It is noted that the chain rule applies in all spatial derivatives if $a_{1}$ is spatially dependent.

We argue that the construction of this filter promotes some of the issues associated to hybrid RANS/LES methods, such as in the transition region. For any blending factor $a_{1}<1$, the resolved LES scales in $\bm{u}_{H}$ are not allowed to evolve with their full magnitude as they are multiplied by a factor smaller than one, which induces a damping of the resolved turbulence. Therefore, the smallest resolved scales would be larger than the given grid filter size and the method does not exhibit optimal efficiency in resolving turbulent scales. It is one of the primary incentives of this work to construct a numerical method which allows the LES scales to evolve independently, only limited by the coarseness of the LES filter size. This is achieved by proposing a new formulation of the hybrid RANS/LES filter in the following, which allows the use of $a_{1}=1$ throughout.

2.2.2 A new three-level additive hybrid RANS/LES filter

We modify the classical additive hybrid RANS/LES filter by adding a third term, yielding

[TABLE]

for the filtered velocity and pressure. Herein, the third filtering operator $\langle\bar{(\cdotp)}\rangle$ defines a hierarchical filter by successively applying an LES filter and a statistical filter. This particular construction enables the choice of the constant factors $a_{1}=1$, $a_{2}=1$, and $a_{3}=-1$ throughout the near-wall region. The third component subtracts the statistical average of the resolved variables from the balance. Due to $a_{1}=1$, resolved turbulent structures are not damped artificially by the method. The eddy resolving component $\bar{\bm{u}}$ may be significantly underresolved as we place the first off-wall node in the region of $y_{1}^{+}\sim 1...100$, so the energy-carrying scales do not have to be resolved in the near-wall area. The resulting filtering operator fulfills the requirements:

[TABLE]

and equivalently for the pressure.

The present hybrid filter may be related to the NLDE, where full LES and RANS solutions overlap in a part of the domain23, 24, 25. The NLDE is formally derived using a hierarchical scale separation operator, applying both an LES filter and a RANS operator, such that $\langle\bar{\bm{u}}\rangle=\langle{\bm{u}}\rangle$ holds by definition. This characteristic stands in contrast to the present approach, since the LES filter size is usually so coarse that the mean flow is not captured sufficiently by the LES scale and we have $\langle\bar{\bm{u}}\rangle\neq\langle{\bm{u}}\rangle$.

When this approach is applied to the incompressible Navier–Stokes equations, the derivatives of $a_{i}$ with respect to the spatial coordinates vanish and the subgrid tensor becomes13:

[TABLE]

with the LES and RANS subgrid tensors

[TABLE]

where $(\cdotp)^{\prime}$ denotes the fluctuating component of the respective quantity, ie $\bm{u}=\langle\bm{u}\rangle+\bm{u}^{\prime}$. In order to simplify the subgrid terms in the following, we define a new velocity variable summarizing the statistical velocity contributions, $\widetilde{\bm{u}}=\langle\bm{u}\rangle-\langle\bar{\bm{u}}\rangle$.

Remark: The quantity $\widetilde{\bm{u}}$ plays a major role in the remainder of this article, as we will explicitly solve for this component in the final numerical method and the velocity solution is considered in the form $\bm{u}_{H}=\bar{\bm{u}}+\widetilde{\bm{u}}$, as a sum of an LES and a RANS component.

Inserting the constants $a_{i}$ in Equation (15) and application of the relation (14), $\langle\bm{u}\rangle=\langle\bar{\bm{u}}\rangle+\widetilde{\bm{u}}$, yields

[TABLE]

In the near-wall area, we assume that Reynolds stresses dominate, allowing the simplification

[TABLE]

Herein, the LES subgrid stresses play a minor role in the vicinity of the wall, as barely any turbulence is resolved due to an overly coarse LES mesh and usually almost the entire solution is represented by $\widetilde{\bm{u}}$. In the present work, we consider the implicit LES capabilities of the numerical method30, 31, where the numerical upwind fluxes introduce an appropriate amount of dissipation in the high-frequency range, so any LES subgrid tensor is omitted in the following. Further remarks on how to include an explicit LES model are given in Section 2.3. As a model for the remaining subgrid terms, we employ an eddy viscosity closure

[TABLE]

with the RANS eddy viscosity $\nu_{t}$. The relevant velocity component, on which the eddy viscosity model acts, is the RANS component $\widetilde{\bm{u}}$, which is a direct result of the formal definition of the three-level filter, as indicated in Equation (22).

Remark: The earlier elaborations on using a constant blending factor of $a_{1}=1$ manifest themselves in the numerical method: the LES degrees of freedom are not directly impacted by the RANS eddy viscosity, so they are not damped.

For simplicity and due to the limitation to the inner layer in the present application, we employ Prandtl’s algebraic mixing length eddy viscosity model including van Driest’s damping function33

[TABLE]

with the mixing length

[TABLE]

where $A^{+}=26$, the norm $|\bm{\epsilon}\left(\widetilde{\bm{u}}\right)|=\sqrt{2\bm{\epsilon}\left(\widetilde{\bm{u}}\right):\bm{\epsilon}\left(\widetilde{\bm{u}}\right)}$, and a modified wall distance $y_{\mathrm{DES}}$. We further note that the Reynolds stresses in Equation (21) include the difference between the full and resolved stress tensor $\langle{\bm{u}}^{\prime}\otimes{\bm{u}}^{\prime}\rangle-\langle\bar{\bm{u}}^{\prime}\otimes\bar{\bm{u}}^{\prime}\rangle$. It would be possible to explicitly compute all terms involving the resolved scales in Equation (21) $\left(-\langle\bar{\bm{u}}^{\prime}\otimes\bar{\bm{u}}^{\prime}\rangle\right)$ by averaging over homogeneous directions and subtracting the result from the eddy viscosity tensor in the philosophy of Medic et al34. In this work we account for the resolved Reynolds stress tensor by considering the classical DES approach of Spalart et al26, which limits the wall distance $y$ with the characteristic grid length scale $h$ times a model constant $C_{\mathrm{DES}}$:

[TABLE]

and we use the latter modification in order to avoid kinks in the residual. For $C_{\mathrm{DES}}$, Spalart et al26 specify a value of the order of unity and a careful calibration using turbulent channel flow yields $C_{\mathrm{DES}}=0.85$ with the solver presented in Section 4. We further employ the constant $\beta=6$ and take $h$ as the wall-normal width of the element $\Delta y_{e}$ divided by $k+1$, $h=\Delta y_{e}/(k+1)$, with $k$ the polynomial degree of the discrete function space of $\bar{\bm{u}}$ (see Section 3.1) and $k+1$ nodes in each spatial direction. We note that a factor of $k$ instead of $k+1$ is sometimes used in the DG context, see, eg, Wurst et al6, but the factor $k+1$ seems to be more consistent with the analysis carried out Moura et al35.

Remark: Taking the wall-normal cell width in Equation (25) stands in contrast to the recommendation in DES of using $\max(\Delta x,\Delta y,\Delta z)$36 as a characteristic cell length. The wall-normal grid size is chosen in this work as we require the flow to be fully turbulent at the interface $\partial\widetilde{\Omega}^{I}$, where the RANS component ends, and the resolved LES stresses are assumed to behave approximately universal in wall-normal direction. The requirement that turbulent eddies have to be resolved by the first off-wall cell limit the cell aspect ratio $\max(\Delta x,\Delta z)/\Delta y_{1}$ in the same way as for wall-stress models. The specific cell aspect ratio requirements are analyzed in detail in Section 5.1.

As the final result we get for the hybrid-filtered Navier–Stokes equations

[TABLE]

where the components of the velocity solution have been expanded in $\bm{u}_{H}=\bar{\bm{u}}+\widetilde{\bm{u}}$. We solve for these two velocity components, an eddy-resolving velocity component $\bar{\bm{u}}$ and a statistical velocity component $\widetilde{\bm{u}}$, explicitly in our numerical method. The idea is that the method computes turbulent eddies where the LES mesh is sufficiently fine. If the mesh is too coarse to resolve the near-wall flow, the method automatically promotes the RANS modes and computes the flow, or parts of it, in a statistical sense. Since the resolution of the pressure is of much less relevance in turbulent boundary layers, the pressure is kept as a single variable filtered with the hybrid RANS/LES filter.

Remark: A straightforward discretization of Equations (26) and (27) with the same function space for $\bar{\bm{u}}$ and $\widetilde{\bm{u}}$ would not be possible, as there are formally seven unknowns (three components each in $\bar{\bm{u}}$ and $\widetilde{\bm{u}}$ as well as $p_{H}$) whereas solely four equations are available. This underdetermination requires the user to make further assumptions on the composition of the solution, in particular that the two components of the velocity field are described with finite numerical resolution and using linearly independent shape functions. An efficient description of this basis should take further knowledge about turbulent boundary layers into account and our choice is presented in Section 3. The general framework of weak forms in a multiscale context, where different shape functions are used for particular scales, is discussed in the subsequent section.

Outside of the near-wall layer $\widetilde{\Omega}$, only the LES scale $\bar{\bm{u}}$ is considered, and the coupling condition on the interface $\partial\widetilde{\Omega}^{I}$ in Figure 1 is $\bm{u}_{H}=\bar{\bm{u}}$. Additionally, we require that the viscous flux by the eddy viscosity term $\bm{\mathcal{F}}^{\nu_{t}}\left(\widetilde{\bm{u}}\right)$ becomes zero in normal direction to avoid a kink in the solution, giving the interface condition

[TABLE]

Since $\nu_{t}$ is in general not zero, this condition is equivalent to $\bm{\epsilon}\left(\widetilde{\bm{u}}\right)\cdotp\bm{n}=\bm{0}$.

2.3 Variational multiscale formulation

The goal of the present subsection is twofold: We explain, on an abstract level, how the solvability of Equations (26) and (27) is enabled despite the seemingly underdetermination of the equation system. This is done by choosing weighting functions for the weak form which are of the same structure as the velocity components and by taking a linearly independent functions for the two scales $\bar{\bm{u}}$ and $\widetilde{\bm{u}}$ with finite dimension each. An efficient set of discrete linearly independent function spaces via function enrichment is presented in Section 3. Furthermore, we explain how the viscous operator in Equation (27) has to be modified in order to enable a physically suitable RANS solution to develop.

For this purpose, we derive an abstract weak form of Equations (26) and (27) in the standard procedure by multiplication of these equations with appropriate weighting functions $\bm{v}\in\mathcal{V}^{\bm{u}}$ as well as $q\in\mathcal{V}^{{p}}$ and integration over the whole spatial domain $\Omega$. As a result we obtain the variational form of the incompressible Navier–Stokes equations as

[TABLE]

Herein, the terms correspond to the transient (mass) term $\mathscr{M}$, convective term $\mathscr{F}^{c}$, pressure term $\mathscr{P}$, the rearranged viscous terms $\mathscr{F}^{\nu}$ and $\mathscr{F}^{\nu+\nu_{t}}$, and the right-hand-side term $l$, as well as the velocity-divergence term of the continuity equation $\mathscr{C}$. All terms are bilinear regarding $\bm{v}$ as well as $\bar{\bm{u}}$ and $\widetilde{\bm{u}}$ except for the convective term, which is nonlinear in the second slot. The detailed definitions of these operators will be elaborated in the context of the solver description in Section 4.

The classical Bubnov–Galerkin method suggests to take weighting functions of the same space as the solution functions. We choose appropriate solution functions $\bar{\bm{u}}\in\mathcal{V}^{\bar{\bm{u}}}$ and $\widetilde{\bm{u}}\in\mathcal{V}^{\widetilde{\bm{u}}}$, assume direct sum decomposition of the underlying spaces $\mathcal{V}^{{\bm{u}}}=\mathcal{V}^{\bar{\bm{u}}}\oplus\mathcal{V}^{\widetilde{\bm{u}}}$, ie, $\left(\bar{\bm{u}}+\widetilde{\bm{u}}\right)\in\mathcal{V}^{{\bm{u}}}$, and we require $\mathcal{V}^{\bar{\bm{u}}}\cap\mathcal{V}^{\widetilde{\bm{u}}}=\{\bm{0}\}$ for linear independence. Employing the same basis for the weighting functions, we obtain two components for $\bm{v}$, $\bar{\bm{v}}=\mathcal{V}^{\bar{\bm{u}}}$ and $\widetilde{\bm{v}}=\mathcal{V}^{\widetilde{\bm{u}}}$, and we have $\bm{v}=\bar{\bm{v}}+\widetilde{\bm{v}}$. By inserting this ansatz into the weak forms (29) and (30), the momentum equation may be decomposed into two equations resembling the two scales of the solution according to the classical variational multiscale paradigm21, 22: an LES scale, weighted with $\bar{\bm{v}}$, and a RANS scale, weighted with $\widetilde{\bm{v}}$, yielding

[TABLE]

The equation for the continuity equation remains unchanged. Through the weighting of Equations (32) and (33) with linearly independent function spaces, the equation system is well posed and resolves the underdetermination discussed in the previous section.

The equation for the LES scale, Equation (32), is essentially a standard LES approach on a background convective velocity $\widetilde{\bm{u}}$, similar to NLDE24, except for an additional eddy viscosity term based on $\widetilde{\bm{u}}$, providing the LES scale with physical dissipation in regions where the energy-carrying scales are not sufficiently resolved. Potential LES subgrid models contained in $\bm{\tau}_{\mathrm{LES}}$ would also be added to this equation.

Equation (33) represents the RANS scale, which resolves a variable fraction of the averaged velocity. It may be observed that the LES solution couples into the RANS equation in (33) and this coupling should be investigated in more detail. Considering again the NLDE methodology as a reference, the RANS equations are well defined without the additional LES source terms. However, the definition of the eddy viscosity in this work (Equation (23)) takes the resolved Reynolds stresses of the LES arising in the convective term into account, such that the full convective flux should be included, and thus the transient term. In contrast, there is no physical justification for the viscous LES part. Indeed, our numerical tests revealed that the method tends to promote an unphysical RANS mode, such as a RANS solution directed adverse to the primary flow direction even in attached boundary layers, if the viscous LES term is included in the RANS equation. Further arguments to be considered are that we expect the viscous RANS flux to be dominant in the underresolved LES region $\left(|\bm{\mathcal{F}}^{\nu}\left(\bar{\bm{u}}\right)|\ll|\bm{\mathcal{F}}^{\nu+\nu_{t}}\left(\widetilde{\bm{u}}\right)|\right)$. Additionally, numerical stability issues according to a coercivity analysis (in Section 4.3) suggest that this term would significantly restrict the range of application if taken into account. We conclude that the viscous LES term should be canceled from the RANS equation according to

[TABLE]

in replacement of Equation (33). It is noted that this modification solely has an impact on one component of the equations and therefore primarily represents a measure for helping the method to find a physically suitable composition of $\bar{\bm{u}}$ and $\widetilde{\bm{u}}$ within the solution space. Similar modifications in the framework of the variational multiscale method are frequently considered regarding scale separation into large and small resolved/unresolved scales in LES turbulence modeling, see eg the review article37.

In this section, a consistent set of governing equations has been derived, which decomposes the incompressible Navier–Stokes equations into a RANS scale and an LES scale. We have argued that a linearly independent function space is required for these two components in order to guarantee solvability of the equation system. In the next section, we show how the concept of function enrichment may be used to construct a highly efficient function space for $\bar{\bm{u}}$ and $\widetilde{\bm{u}}$.

3 RANS and LES velocity components via function enrichment

The concept of function enrichment was introduced by Krank and Wall14 within the standard FEM as a way to economize the computation of turbulent boundary layers. This approach has since been applied to the discontinuous Galerkin method for DES16 and the RANS equations by Krank et al15; the present section is closely related to the latter articles. The basic idea is to construct additional shape functions using a wall-law as enrichment function, which enable the use of very coarse meshes in the near-wall region along with exact boundary conditions since the function space is capable of resolving the sharp boundary layer gradient present in the mean velocity. In the current work, this capability is used to resolve the RANS scale according to Section 2 in a layer near the wall with very few degrees of freedom. The high-order polynomial component additionally represents eddies where the mesh is sufficiently fine. Since these basis functions are linearly independent, they may be used to approximate the RANS and LES scales in the variational multiscale method described in Section 2.3. The major idea is summarized in Section 3.1 and the enrichment function employed in this article is presented in Section 3.2.

3.1 Function enrichment

We consider a tessellation of the three-dimensional spatial computational domain $\Omega_{h}\subset\mathbb{R}^{3}$ into $N_{e}$ nonoverlapping hexahedral finite elements $\Omega_{e}$, given as $\Omega_{h}=\bigcup_{e=1}^{N_{e}}\Omega_{e}$. The subscript $(\cdotp)_{h}$ relates the respective variable to a characteristic element size $h$. The discrete velocity solution is denoted as $\bm{u}_{h}$ in the following and, within the first element layer near the no-slip wall, it consists of two parts, a high-order polynomial component $\bar{\bm{u}}_{h}\left(\bm{x},t\right)\in\mathcal{V}_{h}^{\bar{\bm{u}}}=\left(\mathcal{V}_{h}^{p}\right)^{3}$ and an enrichment component $\widetilde{\bm{u}}_{h}\left(\bm{x},t\right)\in\mathcal{V}_{h}^{\widetilde{\bm{u}}}$, yielding

[TABLE]

These two components are of course related to the respective RANS and LES contributions of the solution as introduced in the previous section. In each element, we use a polynomial of degree $k$ to represent the LES component of both $\bar{\bm{u}}_{h}$ and $p_{h}$, for example for the pressure:

[TABLE]

where $P_{k}$ are all polynomials up to the tensor degree $k$ and the solution allows for discontinuities at the element boundaries. Written in terms of a finite element expansion, the standard polynomial component in each element becomes

[TABLE]

where $N_{B}^{k}$ is the $B^{\mathrm{th}}$ shape function of degree $k$ and $\bar{\bm{u}}_{B}$ is the corresponding nodal value. The enrichment may analogously be written in terms of an FE expansion of tensor degree $l$, multiplied by an enrichment function $\psi$ according to

[TABLE]

this expansion corresponds to the general definition of $\widetilde{\bm{u}}_{h}$ in $\widetilde{\Omega}_{h}\subset{\Omega}_{h}$:

[TABLE]

and we set $\widetilde{\bm{u}}_{h}=\bm{0}$ outside $\widetilde{\Omega}_{h}$. As indicated, $\widetilde{\Omega}_{h}$ corresponds to the RANS/LES zone in Figure 1, and, in the present higher order context (eg $k=4$), this layer is chosen equivalent to the width of the first off-wall cell. If $k\leq 2$, two cell layers may be considered. This construction allows the numerical method to resolve solutions similar to the enrichment function with very coarse meshes. At the same time, full consistency is maintained and the method has sufficient flexibility in flow situations where the solution differs from the wall-law through the high-order polynomial component. In the context of boundary-layer enrichments, we therefore employ a wall function as enrichment function, which will be introduced in the subsequent section. It is noted that the blending of function spaces is not an issue within the discontinuous Galerkin method since neighboring elements may have nonconforming shape functions by default, due to the weak element coupling. As a consequence, ramp functions used within the standard FEM14 are not necessary. The resulting decomposition of $\bm{u}_{h}$ into the two components is illustrated in Figure 2. The discrete pressure variable is adequately represented by a high-order polynomial as high gradients are not expected in the pressure distribution.

3.2 Enrichment function

In our previous work14, 15, 16, we employed Spalding’s law as the enrichment function, which has proven to be suitable both for RANS, DES, and LES. Preliminary investigations of the present wall model have indicated the significant potential of using the particular wall function which is consistent with the eddy viscosity model presented in Section 2.2 in the viscous sublayer and the buffer layer ($y^{+}<30$). For example, the predictions of the wall shear stress could be enhanced in accuracy and robustness. This wall function has been proposed by van Driest33 and fulfills the boundary conditions $u^{+}(y=0)=0$ and $\frac{\partial u^{+}}{\partial y^{+}}|_{y=0}=1$. The enrichment function is defined as

[TABLE]

with $u^{+}=\psi$, and we take $\kappa=0.41$ and $A^{+}=26$. Since the enrichment function has to be evaluated on every quadrature point of the numerical method in a layer of elements near the wall, and an efficient evaluation may not be obvious, we describe a possible approach in Appendix A.

As the enrichment function is universal with respect to the wall coordinate $y^{+}=y\sqrt{\tau_{w}/\rho}/\nu$, given the wall shear stress $\tau_{w}$ and the density $\rho$, the wall coordinate scales the enrichment function in wall-normal direction to match the local wall shear stress. We have developed an adaptation algorithm14, 15 which allows for spatial and temporal variations of the flow in the wall model by discretizing the wall shear stress $\tau_{w}$ based on linear continuous shape functions $N_{B}^{c,m}$ ($m=1$): The wall shear stress is computed on each node $B$ in $N^{c,m}$ similar to15 with

[TABLE]

where $\bm{\tau}_{w}$ is the wall shear stress vector, which may be computed as $\bm{\tau}_{w}=\rho\nu\left(\nabla\bm{u}_{h}\right)\cdotp\bm{n}$. Note that we compute the integral for each component of the numerator and subsequently take the vector norm. The nodal values are interpolated using

[TABLE]

The quantity $\tau_{w,h}$ includes a spatial coarsening of the wall shear stress distribution similar to14. This coarsening is physically motivated since the wall function represents a relation for the mean quantities, ie, the mean wall shear stress gives a mean velocity profile. Without this coarsening, one would observe a statistical overprediction of $\tau_{w}$. The wall distance is likewise defined using a finite element expansion where the weights $y_{B}$ are given as the shortest distance between the current and the closest wall node,

[TABLE]

which facilitates the evaluation of the enrichment shape functions and their spatial gradients. The $y^{+}$-variable in the mixing length eddy viscosity model in Equation (24) is also computed based on these definitions of $\tau_{w,h}$ and $y_{h}$.

The wall shear stress is updated at the beginning of every time step using the velocity solution of the previous time step for $\tau_{w,h}$, resulting in a new FE space in every time step. A cell-wise $L^{2}$-projection is employed to project the solution vector onto the new function space with

[TABLE]

where $n$ is the current time step number. Supplementary solution vectors required by the time integration scheme are projected as well or recomputed from the velocity solution $\bm{u}_{h}^{n,\mathrm{new}}$, if possible. Instead of using a wall shear stress based on the solution of the previous time step, it would also be possible to use a higher order estimate such as an extrapolated wall shear stress based on a series of previous time steps for improved accuracy. Since the ‘lag’ in $\tau_{w,h}$ is small in this work due to the semi-explicit formulation of the scheme, such an extrapolation is not considered.

Finally, the present adaptation algorithm also allows to switch off the wall model when not needed, ie, when the standard polynomial component is sufficient to resolve the necessary scales according to15. In particular, one observes problems in conditioning when the enriched element spans a $y^{+}$-range smaller than approximately 15 wall units, depending on the respective choice of $k$ and $l$, as the standard and enrichment shape functions become close to linearly dependent when $\tau_{w}\to 0$, yielding $\psi\to 0$. This issue is circumvented by switching off the enrichment and the complete multiscale wall model if all quadrature points of the wall-nearest cell lie in $y^{+}<30$. A parameter study regarding the latter quantity, using turbulent channel flow with the setup applied in Figure 8, yielded almost indistinguishable results for values in the range of 20 to 30. In case $y^{+}$ becomes larger than 30 at any quadrature point of this cell in a subsequent time step, the wall model is switched on again by taking the enrichment in the solution of the new time step into account and an appropriate RANS solution develops automatically within a few time steps.

4 Application to a high-order discontinuous Galerkin solver

The present multiscale wall model may in principle be implemented in any discontinuous Galerkin flow solver. In order to exemplify the approach, we employ our high-performance matrix-free incompressible code INDEXA30. An extension to compressible flow, where the energy equation has to be discretized in addition, is straightforward since, in adiabatic boundary layers, high gradients are not expected in the energy variable, allowing it to be treated in the same way as the pressure variable by the high-order polynomial space. The present section gives an overview of the temporal discretization (Section 4.1) and the spatial operators of the variational form (Section 4.2), and presents a numerical stability analysis of the viscous multiscale term (Section 4.3). Finally, we discuss possibilities for the efficient implementation of the enriched elements (Section 4.4).

4.1 Temporal dual splitting scheme

The incompressible Navier–Stokes equations are integrated in time using the high-order semi-explicit dual-splitting scheme by Karniadakis et al38, which decomposes the governing equations into three sub-steps. In the first step, the convective term is treated explicitly as

[TABLE]

yielding the first intermediate velocity $\hat{\bm{u}}$. Herein, $\Delta t$ denotes the time increment, $n$ the time step, $J$ the order of accuracy, and $\alpha_{i}$, $\beta_{i}$, as well as $\gamma_{0}$ time integration constants of the extrapolation and backward-differentiation (BDF) schemes. The explicit character of this step restricts the time step size according to the Courant–Friedrichs–Lewy (CFL) condition. The time step size is chosen adaptively in order to maximize the time increments in each step with the accurate CFL condition15, and the time integration parameters are calculated to match variable time increments. We take a value for the Courant number of $\mathrm{Cr}=0.14$ and a second order accurate scheme ($J=2$) for all examples considered in this article. Unlike $k$-adaptivity, we have not observed time step restrictions through the enrichment.

In the second step of the time integration scheme, a Poisson equation is solved for the pressure $p^{n+1}$ at the time $t^{n+1}$,

[TABLE]

We adopt the high-order boundary conditions for the Poisson problem38,

[TABLE]

on $\partial\Omega^{D}$, where we have made use of the boundary conditions $\bm{g}_{\bm{u}}=\bm{0}$ and $\nu_{t}=0$. Only the solenoidal viscous part is taken into consideration38, and the vorticity $\bm{\omega}$ is computed as in Krank et al15, 30. We further note that the present dual-splitting scheme allows equal-order interpolation of velocity and pressure39. With the pressure available, the velocity field is made divergence-free in a local projection step, given as

[TABLE]

resulting in the second intermediate velocity $\hat{\hat{\bm{u}}}$. Finally, the viscous effects are taken into account in a Helmholtz-like equation,

[TABLE]

yielding the velocity solution $\bm{u}^{n+1}$ at time $t^{n+1}$.

4.2 Galerkin formulation

Variational formulations are derived for all steps of the time integration scheme (Equations (44)–(48)) and the weak forms are related to the symbolic variational formulations discussed in Section 2.3. These weak forms are adapted from15, 30 with major modifications of the viscous term summarized in the following. They consist of element-wise volume integrals and integrals over interior faces $\partial\Omega_{e}^{\Gamma}=\partial\Omega_{e}^{-}\cap\partial\Omega_{e}^{+}$ between two adjacent elements $\Omega_{e}^{-}$ and $\Omega_{e}^{+}$. Unit normal vectors point outwards of the respective cells and are denoted $\bm{n}_{\Gamma}^{-}$ and $\bm{n}_{\Gamma}^{+}$ such that $\bm{n}_{\Gamma}^{-}=-\bm{n}_{\Gamma}^{+}$ holds. Jumps in the solution at the interfaces between two neighboring elements are abbreviated with $[\phi]=\phi^{-}-\phi^{+}$ and $\llbracket\phi\rrbracket=\phi^{-}\otimes\bm{n}_{\Gamma}^{-}+\phi^{+}\otimes\bm{n}_{\Gamma}^{+}$. Likewise, an averaging operator defines the mean value $\{\!\{\phi\}\!\}=w^{-}\phi^{-}+w^{+}\phi^{+}$ with the weights $w^{-}=w^{+}=1/2$ unless specified otherwise. Boundary conditions are imposed by prescribing external values for $\phi^{+}$ and $\nabla\phi^{+}$ on $\partial\Omega^{D}$ and interface conditions on $\partial\widetilde{\Omega}^{I}$. The weak form of the multiscale model further makes use of the definitions $\widetilde{\Omega}_{e}=\Omega_{e}\cap\widetilde{\Omega}$ (enriched cell interior), $\partial\widetilde{\Omega}_{e}^{I}=\partial\Omega_{e}\cap\partial\widetilde{\Omega}^{I}$ (face on $\partial\widetilde{\Omega}^{I}$, exactly one face enriched), and $\partial\widetilde{\Omega}_{e}=\partial\Omega_{e}\cap\widetilde{\Omega}\backslash\partial\widetilde{\Omega}_{e}^{I}$ (both adjacent faces enriched). The superscripts $(\cdotp)^{-}$ and $(\cdotp)^{+}$ are omitted in the remainder of this section if possible without ambiguity. The subsequent weak forms employ the common notation for $L^{2}$ inner products abbreviating volume integrals with $(a,b)_{\Omega_{e}}=\int_{\Omega_{e}}ab\ \mathrm{d}\Omega$ for scalars, $(\bm{a},\bm{b})_{\Omega_{e}}=\int_{\Omega_{e}}\bm{a}\cdotp\bm{b}\ \mathrm{d}\Omega$ for vectors, and $(\bm{a},\bm{b})_{\Omega_{e}}=\int_{\Omega_{e}}\bm{a}:\bm{b}\ \mathrm{d}\Omega$ for tensors. The notation for face integrals is analogous.

We derive element-wise weak forms for each term of the Navier–Stokes equations by multiplying the respective term with a weighting function and integrating over one element volume:

• •

Mass term: The mass term is obtained without further modification as

[TABLE]

Neighboring cells are not connected due to the absence of face terms.

• •

Convective term: The weak form is integrated by parts, yielding

[TABLE]

for the convective term. The convective flux is defined via the local Lax–Friedrichs flux

[TABLE]

and we take the maximum of $\Lambda=\max(\lambda^{-},\lambda^{+})$ across the interface of the largest eigenvalue of the flux Jacobian (see Reference30)

[TABLE]

• •

Pressure gradient: The pressure gradient is integrated by parts as suggested in Krank et al30, yielding

[TABLE]

• •

Velocity divergence: The velocity divergence is integrated by parts according to30 as well and results in

[TABLE]

• •

Viscous multiscale term: We discuss a suitable viscous term in view of the modifications introduced in Section 2.3. The baseline viscous implementation of our solver is an interior penalty method in symmetric form in the standard incompressible flow solver30 and in nonsymmetric form in the enriched solver15. The nonsymmetric variant of the interior penalty method has the advantage of being stable with very low requirements on the penalty parameter40, which is beneficial in the case of nonpolynomial shape functions, such as the present boundary layer enrichment, where the derivation of inverse estimates on the fly is not economical.

The viscous term considered in this work follows Reference15 and the impact of the LES solution $\bar{\bm{u}}$ on the equations for the RANS scale weighted with $\widetilde{\bm{v}}$ are canceled as suggested in Section 2.3. We propose the following formulation of the viscous multiscale term:

[TABLE]

where all terms are expanded into the individual contributions from the RANS and LES scale. Herein, the row ($\widetilde{\bm{v}},\widetilde{\bm{u}}$) represents the RANS contribution to the RANS equation, the row ($\widetilde{\bm{v}},\bar{\bm{u}}$) the LES contribution to the RANS equation, the row ($\bar{\bm{v}},\widetilde{\bm{u}}$) the RANS contribution to the LES equations and finally the last row the LES contributions to the LES scale. The columns correspond to the respective volume terms, the adjoint face terms, the standard consistency face terms, and the penalty terms. We note the skew-symmetry of the standard consistency and adjoint face terms as well as the symmetry of the penalty terms.

In Equation (55) we do not entirely cancel the LES impact on the RANS equations but shift the adjoint term at the inner enriched faces $\partial\widetilde{\Omega}_{e}$ from the row ($\bar{\bm{v}},\widetilde{\bm{u}}$) to the row ($\widetilde{\bm{v}},\bar{\bm{u}}$). This modification is consistent, since it relies on the discontinuity of the velocity solution at element boundaries, and is required in order to guarantee a stable scheme through skew-symmetric face terms. The skew-symmetry is beneficial in the context of the coercivity analysis in the subsequent section. In contrast, on the interface between enriched and nonenriched cells $\partial\widetilde{\Omega}^{I}_{e}$, the opposite face terms are considered, namely the adjoint term in row ($\bar{\bm{v}},\widetilde{\bm{u}}$) and the skew-symmetric standard consistency term in row ($\widetilde{\bm{v}},\bar{\bm{u}}$) again in order to allow the fulfillment of the coercivity argument. The resulting Galerkin formulation is in agreement with the interface conditions $\bar{\bm{u}}^{+}=\bar{\bm{u}}^{-}+\widetilde{\bm{u}}^{-}$ and $\widetilde{\bm{u}}^{+}=\bm{0}$ on $\partial\widetilde{\Omega}_{e}^{I}$ (ie, $\llbracket{\bm{u}}\rrbracket=\left(\bar{\bm{u}}^{-}+\widetilde{\bm{u}}^{-}-\bar{\bm{u}}^{+}\right)\otimes\bm{n}_{\Gamma}$ and $\llbracket\widetilde{\bm{u}}\rrbracket=\widetilde{\bm{u}}^{-}\otimes\bm{n}_{\Gamma}$), where the left face $(\cdotp)^{-}$ is enriched and the right face $(\cdotp)^{+}$ is not enriched, without loss of generality. However, the condition $\bm{\epsilon}\left(\widetilde{\bm{u}}^{-}\right)\cdotp\bm{n}=\bm{0}$ is not exactly fulfilled as we have $\bm{\mathcal{F}}^{\nu+\nu_{t}}\left(\widetilde{\bm{u}}^{-}\right)\cdotp\bm{n}_{\Gamma}^{-}=\{\!\{\bm{\mathcal{F}}^{\nu}\left(\bar{\bm{u}}\right)\}\!\}\cdotp\bm{n}_{\Gamma}$ on $\partial\widetilde{\Omega}_{e}^{I}$, but the error is considered to be small since kinks have not been observed in the solution at that boundary. As an alternative, our numerical tests have shown that the standard consistency term in row ($\widetilde{\bm{v}},\bar{\bm{u}}$) may be neglected, which yields a better fulfillment of the latter interface condition while problems in stability have not been observed with this modified formulation. The numerical examples shown in Section 5 yet employ the variant with proven stability.

We further note that the material parameter $\nu$ is applied to the interior penalty term in the whole domain. As interior penalty stabilization parameter $\tau_{\mathrm{IP}}$ we use the definition by Hillewaert41 (see also References15, 30), since this choice has yielded favorable results in underresolved turbulent flows (ILES) in Krank et al30. On the Dirichlet boundary, we have taken measures to enhance the accuracy of the weakly enforced no-slip condition15. In the present work, we increase the interior penalty parameter $\tau_{\mathrm{IP}}$ on all Dirichlet boundary faces by a factor of $10$, which has a similar effect while being more solver-friendly in 3D. The weights of the averaging operators with spatially varying material parameter $\nu_{t}$ included in the $\bm{\mathcal{F}}^{\nu+\nu_{t}}$-terms are given through harmonic weighting15, 42:

[TABLE]

On $\partial\widetilde{\Omega}^{I}$, all terms including the eddy viscosity vanish, making a consideration of the varying material law unnecessary.

• •

Velocity div-div penalty: The present scheme requires a div-div penalty operator for stabilization of mass conservation, see Reference30 for a detailed discussion. The corresponding operator reads

[TABLE]

We note that the contribution of this term vanishes if the velocity field is exactly divergence-free, as it relies on the continuity residual $\nabla\cdotp{{\bm{u}}}$. The stabilization parameter $\tau_{\mathrm{D}}$ is given as30

[TABLE]

where $\|{\bm{u}}^{n}\|$ is the volume-averaged velocity magnitude of the previous time step, $h$ the element length based on the cube root of the element volume and Cr the global Courant number.

• •

Pressure Laplace term: The pressure Poisson equation further requires the discretization of a Laplace term. We use the symmetric interior penalty method43, yielding

[TABLE]

For the numerical flux $\bm{\mathcal{P}}^{*}$, we include an interior penalty term on internal faces and use the pressure boundary conditions on the Dirichlet boundary (46) as

[TABLE]

The same interior penalty parameter definition as for the viscous term is employed41 and the vorticity vector $\bm{\omega}_{h}\in\mathcal{V}^{{\bm{u}}}_{h}$ is precomputed via $L^{2}$-projection30.

Finally, these weak operators allow a compact definition of the dual-splitting scheme in element-wise Galerkin form. The explicit convective step reads

[TABLE]

The pressure Poisson equation is

[TABLE]

and as boundary condition for $\hat{\bm{u}}_{h}$ on $\partial\Omega^{D}$ we employ the interior value $\{\!\{\hat{\bm{u}}_{h}\}\!\}=\hat{\bm{u}}_{h}^{-}$. The element-wise projection step is obtained by adding the div-div penalty term (57) to the left-hand side, resulting in

[TABLE]

where boundary conditions on the pressure variable are again applied using the internal value $\{\!\{p_{h}^{n+1}\}\!\}=p_{h}^{n+1,-}$ on $\partial\Omega^{D}$. The viscous step may be written as

[TABLE]

Boundary conditions for all remaining quantities on no-slip walls are specified using the so-called mirror principle44, defining the exterior velocity contribution with $\bm{u}_{h}^{+}=-\bm{u}_{h}^{-}+2\bm{g}_{\bm{u}}=-\bm{u}_{h}^{-}$ on $\partial\Omega^{D}$.

4.3 Coercivity analysis

We sketch a coercivity analysis of the viscous term in order to investigate the stability of our numerical scheme. The Lax–Milgram theorem ensures solvability of the variational Helmholtz problem (64) if there is a constant $C\geq 0$ such that

[TABLE]

holds. Herein, we disregard the mass terms, as they always yield a positive contribution on the left-hand side of the inequality and the Helmholtz equation degenerates to an $L^{2}$ projection if $-\mathscr{F}^{\nu}_{e}\to 0$.

Due to their skew-symmetric construction, the face terms on $\partial\widetilde{\Omega}_{e}$ in Equation (55) cancel each other when the solution and weighting functions are equal. Solely the interior penalty face terms remain, which have an entirely positive contribution. Therefore, the following inequality holds for enriched cells,

[TABLE]

leaving only the volume terms for further consideration. The first and the last volume term are each symmetric and thus positive while the nonsymmetric second term is problematic as it can yield a negative contribution. In the following, we derive an estimate for this term in order to guarantee that the sum of all terms in the bilinear form is bounded from below by $C\left(\|\bar{\bm{v}}\|^{2}+\|\widetilde{\bm{v}}\|^{2}\right)$.

Young’s inequality

[TABLE]

for any $\varepsilon>0$ applied to the second volume term in (66) results in:

[TABLE]

Herein, we have introduced the modified viscosity term $\widehat{(\nu+\nu_{t})}$. Inserting this relation in inequality (66) yields an estimate for $-\mathscr{F}^{\nu}_{e}$ based on the two symmetric terms,

[TABLE]

allowing the conclusion that the Lax–Milgram theorem is fulfilled under the following conditions:

[TABLE]

Upon rearrangement, these inequalities directly result in the condition for coercivity of the viscous multiscale term:

[TABLE]

This means that the amount of viscous dissipation introduced in the LES scale by the RANS solution has to be limited if ${\nu+\nu_{t}}>4\nu$ for reasons of stability. We apply the modified eddy viscosity

[TABLE]

in the volume term of the row ($\bar{\bm{v}},\widetilde{\bm{u}}$) in (55), which comes along with a minor limitation of the application range of the wall model. This aspect will be investigated in detail in Section 5.1. We anticipate at this point that the width of the first off-wall cell should not exceed a $y^{+}$-range of approximately 120 wall units in the statistical data. The relation (72) guarantees the stability of the scheme in outliers due to turbulent fluctuations. We emphasize that this behavior is related to the particular discrete method used for the viscous multiscale term in the present work and we encourage mathematicians to develop alternative formulations in order to weaken or remove this limitation.

For completeness, we remark on the relations for the case with the viscous LES term according to Equation (33), ie, without the modification introduced in Equation (34) and therefore including all terms in the rows ($\widetilde{\bm{v}},\bar{\bm{u}}$) and ($\bar{\bm{v}},\widetilde{\bm{u}}$) in $\widetilde{\Omega}_{e}$ as well as on $\partial\widetilde{\Omega}_{e}$ (in Equation (55)). The nonsymmetric term on the right hand side of inequality (66) would get a material factor of $2{(2\nu+\nu_{t})}$ instead of $2{(\nu+\nu_{t})}$, resulting in the condition $\widehat{(\nu+\nu_{t})}\leq 2\sqrt{\nu(\nu+\nu_{t})}-\nu$ for the material parameter of the volume term in the row ($\bar{\bm{v}},\widetilde{\bm{u}}$). This condition limits the growth of the eddy viscosity variable if ${\nu+\nu_{t}}>\nu$, which would represent a too strong restriction in practical applications.

4.4 Implementation

The Galerkin formulations of the dual splitting scheme, Equations (61)–(64) are integrated by numerical quadrature, yielding a matrix formulation for each sub-step. The matrix forms are similar to the ones described in earlier studies15, 30. The numerical evaluation and integration is performed using the high-performance sum factorization kernels by Kronbichler and Kormann45 available in the deal.II finite element library46. Herein, the nonenriched cells are evaluated with Gaussian quadrature rules of appropriate order such that all integrals are evaluated exactly on affine cells ($k+1$ quadrature points for linear terms and $\left\lfloor\frac{3k}{2}\right\rfloor+1$ points for nonlinear terms). Enriched cells pose higher requirements on the accuracy of the quadrature rules due to the presence of the nonpolynomial shape functions. Therefore, $8^{3}$ to $12^{3}$ quadrature points have been used for the examples in the present paper, where more points have been employed if the enriched cells span a wider $y^{+}$-range. Face integrals are evaluated with the same number of quadrature points per direction if at least one of the two adjacent cells is enriched. In an earlier work15, we have shown how an object-oriented programming language with template capabilities (eg C++) may be used to implement the present boundary layer enrichment with minimal computational overhead in nonenriched cells and zero extra cost in nonenriched simulations in one single implementation of the Navier–Stokes solver.

All sub-steps are implemented in a matrix-free manner, including the iterative solvers for the global Poisson and Helmholtz problems as well as the local projection solver according to Reference30; further details on the implementation of the multigrid solver of the Poisson problem are given in Reference47. In particular, we have developed an algorithmic toolbox for implementing enrichments in such a matrix-free context at very little additional cost. For example, these algorithms allow a much faster application of the inverse mass matrix than the approach via LU factorization chosen in Krank et al15. Furthermore, most of the additional cost of the higher quadrature rules used in enriched cells can be circumvented. Since a detailed description and performance analysis of these algorithms is beyond the scope of this article, it is the subject of a separate publication48; see also the monograph 17.

5 Numerical examples

The wall modeling approach is validated with turbulent channel flow (Section 5.1) as well as flow past periodic hills (Section 5.2). These two benchmark examples provide insight in the performance of the wall model both in attached and separated boundary layers. We follow our earlier recommendations of using a polynomial degree of $k=4$ for the discretization of the LES scale, since this choice yields a good compromise between accuracy and time-to-solution within the present solver15, 30, 16, 31. In addition, a polynomial of $4^{\mathrm{th}}$ degree inside the wall layer is capable of resolving sufficient turbulence within a single layer of cells. The influence of the degree of the polynomial weighting of the enrichment function, $l$, is investigated and guidelines for the use of the wall model are formulated.

5.1 Turbulent channel flow

As a first validation example, we consider turbulent flow in a rectangular channel of the dimensions $2\pi\delta{\times}2\delta{\times}\pi\delta$ in streamwise, vertical, and spanwise direction, respectively, with the channel half-width $\delta$. Periodic boundary conditions are applied in the streamwise and spanwise direction and no-slip boundary conditions at the top and bottom wall. The grid is graded towards the walls by the hyperbolic mapping given as $f$: $[0,1]\to[-\delta,\delta]$:

[TABLE]

with the mesh stretching parameter $\gamma$. The wall model is active in the first off-wall element layer and a typical mesh is shown in Figure 3. The initial conditions are solely applied on the polynomial degrees of freedom of the solution and the RANS component develops in the course of the simulation once $y^{+}>30$ at any quadrature point of a wall-layer cell according to Section 3.2. The flow is driven by a constant pressure gradient derived from the nominal flow quantities and the results are normalized using the numerical wall shear stress in the friction velocity $u_{\tau}=\sqrt{\tau_{w}/\rho}$. We consider the friction Reynolds numbers $Re_{\tau}=u_{\tau}\delta/\nu$ in accordance to DNS reference data at $Re_{\tau}=395$49, $Re_{\tau}=950$50, $Re_{\tau}=2{,}000$51, and $Re_{\tau}=5{,}200$52. One snapshot of the turbulent flow is visualized in Figure 4 via contours of the velocity magnitude and turbulent vortex structures are made visible using iso-surfaces of the Q-criterion. All simulation parameters and discretization cases are summarized in Table 1. Statistics were averaged in a time interval corresponding to approximately 60–80 flow-through times based on a fixed time interval.

The results of a first application of the wall model including all four Reynolds numbers are depicted in Figure 5. Overall, excellent agreement is observed both for the mean velocity, the RMS velocities and the Reynolds shear stress. Furthermore, it is apparent that the results do not depend on the Reynolds number. Regarding the mean velocity, the enrichment represents the whole mean velocity and the time-averaged LES solution is almost zero in the wall layer. Therefore, the LES result is only visualized for the lowest Reynolds number and not considered in the remainder of this section. The Reynolds shear stress is computed based on the resolved fluctuations, which explains the small gap to the reference data in the near-wall zone. In particular the $u^{\prime+}$ and $w^{\prime+}$ curves exhibit small ticks at the element interfaces, which is a typical result for coarse meshes in DG and has been observed in an earlier study as well30. These ticks vanish with increasing resolution30.

Since the primary incentive for the development of the present multiscale wall model was the deficiencies of DES in the hybrid RANS/LES transition region, we directly compare the two simulation methodologies. To this end, we consider DDES including wall modeling via function enrichment16 for the spatial discretization, such that the exact same meshes can be used for both models. Figure 4 shows a qualitative comparison of two simulations at $Re_{\tau}=950$; they deviate drastically. DDES operates in RANS mode in the inner layer up to $y/\delta\approx 0.05$ and the eddy viscosity acts on the polynomial velocity component as well, such that no vortices are resolved in the inner layer. This stands in contrast to the present multiscale wall model, which computes turbulent motions in the inner boundary layer as well. In DDES, turbulent eddies evolve outside of the inner layer, but the flow generally behaves differently due to the use of the LES eddy viscosity subgrid model given through the Spalart–Allmaras model in LES mode. These two simulations are further compared quantitatively through velocity statistics in Figure 6. As it is expected, DDES overpredicts the mean velocity in the outer layer as a result of a log-layer mismatch caused by the RANS–LES transition. Such a log-layer mismatch is not observed with the present model. The turbulent stresses show that the transition in the DDES model extends until approximately $y/\delta=0.4$, whereas the present multiscale wall model makes much better use of the eddy-resolving capability of the mesh and resolves the important turbulent motions beyond $y/\delta=0.05$.

The wall model including its application range and robustness with regard to grid dependence is analyzed in a systematic manner in the following. In Figure 7, five grid stretching factors are investigated using the same number of cells in each spatial direction. In addition, the number of cells is doubled in the $x_{3}$-direction in one case. The results show that the wall model generally exhibits an excellent robustness regarding the choice of the mesh. The dependence of the cell aspect ratio may be analyzed by considering the measure ${\Delta y_{1e}^{+}(k+1)}/{\max(\Delta x_{e}^{+},\Delta z_{e}^{+})}$, which quantifies the available number of grid points in the wall-parallel direction to resolve a turbulent motion of the size of $\Delta y_{1e}^{+}$; this quantity is included in Table 1 for all simulation cases. Remember that we require sufficient turbulence to be resolved at a distance from the wall $\Delta y_{1e}^{+}$ such that the RANS model can be ‘switched off’, so the cell aspect ratio should not exceed a certain limit. It is noted that other wall modeling approaches, such as wall stress models, in principle have the same requirement if not even more stringent. The factor ${\Delta y_{1e}^{+}(k+1)}/{\max(\Delta x_{e}^{+},\Delta z_{e}^{+})}$ varies in the range $0.56$ to $1.60$. The simulation case with the highest aspect ratio, $ch950\_N16^{3}\_k4l0\_\gamma 1.5$, shows a minor overprediction of the mean velocity where the wall model ends, allowing the conclusion that ${\Delta y_{1e}^{+}(k+1)}/{\max(\Delta x_{e}^{+},\Delta z_{e}^{+})}$ should not go below the limit $\sim 0.5$. According to Jiménez 53, the streamwise size of the energetic eddies in the log-layer is approximately $5y$ so they are well-resolved with $2.5$ grid points at the height $\Delta y_{1e}^{+}$, where the RANS layer ends.

The coercivity analysis in Section 4.3 requires the clipping of the eddy viscosity beyond $\widehat{\nu+\nu_{t}}>4\nu$ in the nonsymmetric volume term of the viscous operator. A normalization of this relation may be obtained by division with the viscosity $\widehat{\nu}_{t}^{+}=\widehat{\nu+\nu_{t}}/\nu>4$. This relation is universal in $y^{+}$-units. It is thus possible to investigate the clipping of $\widehat{\nu}_{t}^{+}$ for one flow configuration and to transfer the conclusions to other simulations. Given that the maximum values of $\widehat{\nu}_{t}^{+}$ grow with the width of the wall layer, it is sufficient to find the upper limit of the application range (cf. Equation (72)). In Figure 7, $\widehat{\nu}_{t}^{+}$ is clipped on approximately $11\%$ of the quadrature points of all wall-layer cells for the case $ch950\_N16^{3}\_k4l0\_\gamma 0.001$. Although the results are excellent for this case, we do not recommend the application of the wall model significantly beyond this $y^{+}$-range in order to guarantee the reliability of the results. In conclusion, the wall model should not extend beyond approximately $120$ wall units in the statistical data. We stress at this point that stability of the numerical method is guaranteed irrespective of the application of the wall model and this $y^{+}$-limit is only due to the result quality, which may degrade for thicker wall layers.

The transition from wall-modeled to wall-resolved LES is investigated in Figure 8. At least $24^{3}$ cells are required for sufficiently resolving the energetic scales in the case without wall model, based on resolution guidelines by Chapman54 and our experience with the present scheme30 (with $k+1$ grid points per cell, we have $\Delta y_{1}^{+}\approx\Delta y_{1e}^{+}/(k+1)\approx 2$, $\Delta z_{e}^{+}\approx\Delta z^{+}/(k+1)\approx 20$). In Figure 8, the grid stretching of the wall resolved LES is successively reduced until the wall model is fully active. The wall-resolved simulation slightly overpredicts the mean velocity and the turbulent fluctuations agree well with the reference. Regarding the case $ch950\_N24^{3}\_k4l0\_\gamma 1.75$, the wall model is taken into account in the algorithm, but the model is switched off during the entire simulation due to the condition described in Section 3.2. The results in Figure 8 are equivalent to the wall-resolved case. For the case $ch950\_N24^{3}\_k4l0\_\gamma 1.5$ the wall model is activated dynamically in a temporally varying fraction of the cells and the enrichment shape functions constitute a small part of the mean velocity. The mean velocity profile agrees well with the DNS, even better than the wall-resolved case. For smaller grid stretching factors, the wall model is fully switched on after the initial laminar-turbulent transition and the enrichment solution is equivalent to the mean velocity profile. The mean velocity in these simulation cases is slightly underpredicted compared to the reference data and the turbulent fluctuations agree well with the reference data. From this numerical test it may be concluded that the wall model is well capable of handling the full range from wall-resolved to wall-modeled LES.

Until this point, all simulation cases have used a constant weighting of the enrichment function ($l=0$) and we put forward reasons for this choice. Figure 9 compares two simulation cases, one with $l=0$ and one with $l=1$, where the latter represents a linear weighting of van Driest’s law. The lines are on top of each other, including the enrichment solution and the fluctuations. This suggests that, for attached boundary layers, there is no need for a weighting of the enrichment with more than a constant factor. Furthermore, the enrichment using linear functions requires additional degrees of freedom ($24$ instead of $3$ per cell) and, most importantly, results in approximately three times the number of linear iterations in the GMRES solver of the Helmholtz equation due to worse condition numbers. Regarding nonequilibrium boundary layers, a linear weighting may yield slightly better results due to the additional flexibility within the RANS solution, but the higher computational cost is not justified.

We conclude the present section of the turbulent channel flow with a performance comparison of wall-resolved and wall-modeled LES. On the one hand, the wall model allows much coarser meshes near the wall, which also comes along with significantly larger time steps due to the explicit time integration, and the lower grid anisotropy yields fewer iterations in the expensive pressure Poisson solver. On the other hand, the wall model requires additional computational effort due to the larger number of quadrature points in the enriched cells and supplementary algorithmic steps. Figure 10 compares the results of a wall-modeled case with the above coarse wall-resolved simulation as well as an additional finer wall-resolved case. The wall-modeled simulation exhibits better agreement with the DNS than the coarse wall-resolved simulation, but worse than the fine wall-resolved case. The computational cost is measured in number of processor cores times wall clock time of the entire simulation and the result is normalized by the wall-modeled simulation cost. The wall-modeled case ($ch950\_N12{\times}8{\times}12\_k4l0\_\gamma 1.6$) gives a cost of 1, the coarse LES simulation ($ch950\_N24^{3}\_k4\_\gamma 2.25$) yields a factor of $77$, and the fine wall-resolved calculation ($ch950\_N32^{3}\_k4\_\gamma 2.0$) a cost of $158$. Comparing the coarse wall-resolved calculation and the wall-modeled calculation in more detail, the use of the wall model reduces the number of time steps by a factor of $5.7$ (for the same simulation time) and reduces the number of cells by a factor of 12. In addition, the wall-resolved simulation requires approximately 22 Poisson solver iterations instead of 9 in the wall-modeled case to yield the same relative accuracy in the iterative solver due to the higher mesh stretching, which overcompensates the extra cost of the wall model. These results allow the conclusion that the wall model promises a speed-up compared to wall-resolved LES by a factor of approximately two orders of magnitude and we expect an even larger benefit for higher Reynolds numbers.

The wall model has been thoroughly investigated regarding its application limits, mesh dependence and performance. The major results are that the model can compute the whole range from wall-resolved to wall-modeled LES and is extremely robust regarding different meshes as well as aspect ratios. However, the wall model should not be used beyond $y^{+}=120$ to guarantee accurate results. Given this limitation, the wall model has accelerated the simulations of turbulent channel flow by a factor of approximately two orders of magnitude.

5.2 Flow over periodic hills

As a second benchmark example, we consider flow over periodic hills at the Reynolds numbers $Re_{H}=10{,}595$, $Re_{H}=19{,}000$, and $Re_{H}=37{,}000$ based on the hill height $H$ and the bulk velocity $u_{b}$. This test case is a popular benchmark for wall models due to its simplicity regarding simulation setup and boundary conditions on the one hand and complexity with respect to turbulence phenomena and modeling on the other hand. For example, common RANS models fail to predict the separation from the curved hill crest15, 56. In addition, the nonequilibrium boundary layer within the recirculation bubble necessitates the consideration of all terms in the Navier–Stokes equations in order to yield accurate results. As an example, the simulations using an equilibrium wall model within a high-order DG method do not converge to the reference data8. In contrast, the results obtained with wall modeling via function enrichment for LES in a standard FEM code14 are very promising, since all terms of the Navier–Stokes equations are taken into consideration.

For this validation case, several reference databases exist. The lower Reynolds number has been computed by LES using multiple codes57, 58, experiments have been carried out55, 58, and the database has recently been extended by a fully resolved DNS using the standard version of the present solver30 at a higher polynomial degree ($7^{\mathrm{th}}$ order of accuracy) and 180 million grid points31. This reference data may be downloaded from the public repository59. Therein, it was found that the DNS is in closer agreement with the experiments than the LES data and predicts the reattachment length at $4.51H$, slightly shorter than the previous LES references. Except for the spatial discretization and temporal order of accuracy, we use the same simulation setup as for the DNS here. Reference data for the higher Reynolds numbers is available via experiments55.

The computational domain for the present simulations is of the size $9H{\times}3.036H{\times}4.5H$ in streamwise, vertical, and spanwise direction, respectively, and the lower wall is smoothly curved in the surrounding of the hill. Since the baseline solver used in this study has already been investigated in detail in the context of wall-resolved LES31 using four polynomial degrees and three mesh refinement levels, we focus on the performance of the wall model in the present work. For most of the simulation cases, a mesh consisting of $32{\times}16{\times}16$ cells of degree $k=4$ is selected, resulting in $160{\times}80{\times}80$ nodes, and the enrichment with $l=0$ is included in the first off-wall element layer, see Figure 11. A single simulation case is presented, in which the number of cells is doubled in each spatial direction, considering the highest Reynolds number. The grid is graded towards the walls and we investigate three grid stretching factors in order to show the influence of the width of the wall layer. According to Figure 12, the wall layer and thus the first off-wall cells span a range in $y^{+}$-units up to $86$ for the lowest, $138$ for the medium, and $156$ for the highest Reynolds number, with peaks near the hill top and minima near the separation and reattachment points. The mesh is mapped onto the exact hill geometry using the same polynomial ansatz as for the LES scale via manifold techniques available in the deal.II library46. The enrichment degrees of freedom constitute equal or less than $0.1\%$ of the overall number of degrees of freedom. Further, the enrichment and thus the multiscale wall model are switched off dynamically if the boundary layer is sufficiently resolved. The wall model is therefore effectively not taken into account behind the hill crest and in the region of the flow reattachment. More cells are switched off in the lower Reynolds number cases; the active cells in one snapshot are depicted in Figure 11. The statistical quantities are averaged during a simulation time corresponding to 61 flow-through times. The velocity field of the fully turbulent flow and turbulent vortex structures are visualized in Figure 13. All discretization cases are summarized in Table 2, including the labels employed in the subsequent figures. The investigations of the wall model include a comparison to one simulation for each Reynolds number, which uses the exact same mesh but without the wall model and is labeled with the addition $NWM$ (no wall model) in Table 2.

We begin the discussion of the results with the friction and pressure coefficients for the lowest Reynolds number. They are defined as

[TABLE]

and for the reference pressure $p_{\mathrm{ref}}$ we take the value at $x_{1}=0$ at the upper wall. The pressure coefficient is multiplied by the density since the kinematic pressure is used in this work. The curves are compared to the DNS reference data in Figure 14. Overall, good agreement of all curves with the DNS is observed and the largest error occurs near the hill top. Even the case without wall modeling is in acceptable agreement with the reference data, although the peak in the skin friction is predicted less distinct in that case. Also, the separation and reattachment lengths in Table 2, given as the zero-crossings of the skin friction, are in good agreement with the reference. The reattachment lengths of the cases $ph10595\_stretch1$ and $ph10595\_stretch2$ at $\sim 4.25H$ are slightly underpredicted in comparison to the DNS reference of $4.51H$. Regarding the case $ph10595\_stretch3$, the reattachment length is predicted even shorter as $3.98H$. With respect to the latter, a possible source of error could be the coarser mesh in the shear layer within the LES region, resulting in a more dissipative behavior of the numerical scheme and thus affecting the length of the recirulation zone.

The time-averaged streamwise velocity, vertical velocity, and Reynolds shear stress are depicted in Figure 15 at ten stations. The curves essentially lie on top of each other, including the wall-modeled simulations, the underresolved case, and the DNS. From this fact we conclude that the wall model shows an outstanding performance for nonequilibrium boundary layers, which confirms the results obtained for wall modeling via function enrichment within the continuous FEM14 and in DES16 as well as the robustness of the method with respect to mesh aspect ratios observed in the previous section.

The fact that even the underresolved case without the wall model yields a good agreement with the reference data motivates an application of the same meshes to a higher Reynolds number. The mean streamwise velocity, vertical velocity, and Reynolds shear stress for $Re_{H}=19{,}000$ are shown in Figure 16. The wall model exhibits results of similar quality as for the lower Reynolds number and the sensitivity regarding the grid stretching is slightly higher. However, the results for the case without wall model do not agree with the reference data at all. In particular, the reattachment length with $1.99H$ instead of $3.94H$ is significantly underpredicted.

Finally, we assess whether the wall model gives any added value in separated flow conditions in comparison to the simpler wall stress model by Carton de Wiart and Murman within a high-order DG method8. Two meshes are employed, the grid with the highest mesh stretching considered previously and a refined variant with twice as many cells in each spatial direction. These meshes are quite similar to the baseline and fine grid used by Wiart and Murman8, labeled as $CM\_WMLES$ and $CM\_WMLES\_fine$, respectively, in the following. Slightly more points are used herein, the approach in the reference8 exhibits an order of accuracy of 8 in comparison to the $4^{\mathrm{th}}$ order accuracy of the present method, however. A full data set including reference data is only available for the mean streamwise velocity, which is presented in Figure 17. The coarser case using wall modeling via function enrichment exhibits solutions of similar quality as the coarse case with the equilibrium wall model. A possible explanation for this agreement may be that some of the relevant scales in the separation region may not be sufficiently resolved, which would also explain the delayed separation length of the case $ph37000\_stretch3$ of $0.37H$. Comparing the refined simulation cases, the equilibrium model does not show a significant improvement in comparison to the coarse cases while the results of the present wall model are almost converged to the experimental reference data. An investigation of the reattachment lengths confirms these observations, as the simulation $ph37000\_stretch3\_fine$ shows the best agreement considering this quantity. This superiority of the function-enrichment-based wall model in separated flow conditions is due to the full consistency of the method and the consideration of all terms of the Navier–Stokes equations. The case without wall modeling does not show a separation region at all as the near-wall region is significantly underresolved.

In summary, the present multiscale wall model enables an accurate computation of separated boundary layers, is robust with regard to mesh sensitivities and exhibits much more accurate results compared to an equilibrium wall stress model.

6 Concluding remarks

In this work, we have developed a new turbulence modeling approach for the wall modeling concept of function enrichment within the high-order discontinuous Galerkin method. Based on a rigorous derivation of the modeling terms using Germano’s framework of additive filtering, a RANS model is applied in a thin layer near the wall. Unlike existing wall modeling approaches, the RANS and an underresolved LES solution overlap inside the near-wall layer. This composition allows the LES solution to develop within the wall layer such that the typical issue of the RANS–LES transition is avoided. The method is also capable of resolving the velocity gradient and thus the wall shear stress in the viscous sublayer with coarse meshes, where the first off-wall cell extends up to $120$ wall units. This is achieved by using the recently proposed technique of function enrichment14, 15, 16, 17 inside the inner boundary layer. The wall model has exhibited outstanding characteristics in attached and separated boundary layers, does not show a log-layer mismatch and has reduced the cost of a benchmark simulation by approximately two orders of magnitude.

The turbulence model used in the original publication on wall modeling via function enrichment within a standard FEM flow solver14 consists of a residual-based approach, which is supported by a structural LES model (multifractal subgrid scales) in the bulk flow. Therein, the idea was that the residual-based method provides numerical dissipation, which is both appropriate to stabilize the scheme and to model the unresolved energetic scales at the wall. Although the data has not been directly compared in this article, we argue based on personal experience with both models that the present methodology with a clean separation of the RANS and LES scales including turbulence modeling (the concept of divide et impera, as noted in the introduction) yields a more reliable and robust approach. Moreover, the present data has been compared with results of a detached-eddy simulation hybrid RANS/LES approach that uses function enrichment in the near-wall area16. That method showed problems in the hybrid RANS/LES transition region, whereas the present method does not.

Acknowledgments

Financial disclosure

The research presented in this paper was partly funded by the German Research Foundation (DFG) under the project “High-order discontinuous Galerkin for the EXA-scale” (ExaDG) within the priority program “Software for Exascale Computing” (SPPEXA), grant agreement no. KR4661/2-1 and WA1521/18-1.

Conflict of interest

The authors declare no potential conflict of interests.

Appendix A Numerical evaluation of van Driest’s law

The enrichment function $\psi$ given as van Driest’s wall-law in Equation (39) has to be computed on each quadrature point during evaluation of the weak forms. The wall function necessitates the evaluation of an integral in the interval $[0,y^{+}]$. In order to avoid computing the integral over the entire interval each time the wall function is evaluated, we tabulate the wall function $\psi_{i}$ at several discrete values of the wall coordinate $y^{+}_{i}$. When the integral is evaluated at run time, solely a small slice of the interval has to be computed, namely $[y^{+}_{i},y^{+}]$ with $i$ such that $y^{+}_{i}\leq y^{+}<y_{i+1}^{+}$:

[TABLE]

This remaining integral is evaluated employing a nine-point Gaussian quadrature rule (with a polynomial order of accuracy of 17) irrespective of the width of the interval in our code. The fixed selection allows for efficient SIMD instructions that include quadrature points in several elements each with different parameters. The tabulated quantities have been computed using 32 digit precision in Matlab and are listed in Table 3 for $\kappa=0.41$ and $A^{+}=26$ and yield a guaranteed relative accuracy for $\psi$ of $1e-12$ up to $y^{+}=5{,}000$. The rapid change of the $y_{i}^{+}$ values shows that the near-wall area necessitates a denser distribution of the support points, while the integrand varies more slowly beyond $y^{+}>1{,}000$, such that only a few quadrature nodes are sufficient. As the wall function is only computed once per time step on each quadrature point and stored, the evaluation cost is negligible. Aiming at fine-tuning the approach, one may consider more support points $y_{i}^{+}$ and a Gaussian rule with fewer support points.

Besides the enrichment function itself, the spatial derivative $\mathrm{d}\psi/\mathrm{d}y^{+}$ is also required for evaluating the weak forms. This derivative is given by the integrand in Equation (74). The derivatives with respect to the spatial location vector $\bm{x}$ are computed starting from Equation (35) by successive application of the chain rule as it is detailed Krank and Wall14.