Implementation and Validation of an Algebraic Wall Model for LES in Nek5000

Turbulent flows are most often wall-bounded, rendering the treatment of the wall essential. In this work, the near-wall layer is modelled in large eddy simulations which enables simulating high Reynolds number flows. An algebraic wall model has been implemented in the spectral element code Nek5000. It consists of an approximate boundary condition that relates the wall shear stress to the velocity measured close to the wall, on the upper edge of the first spectral element. The wall shear stress model approximates the law of the wall for hydraulically smooth cases. The model is applied to channel flow cases at Reτ=1000\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_{\tau }=1000$$\end{document} and at Reτ=5200\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_{\tau }=5200$$\end{document}. The wmLES results obtained with the present implementation are seen to compare very well with those of reference direct numerical simulations in the resolved region. They also remain remarkably close to the reference results for a large part of the under-resolved region; which is not necessarily the case when using low order implementations and even other types of high order discretizations, as found in the literature. Various parameters are studied such as the time averaging, the height of the near wall under-resolved element, and the mesh requirements. The obtained results indicate that accurate results can be obtained with Nek5000 at a reduced cost thanks to this newly implemented wmLES model. This work provides the necessary guidelines for simple flows, and it will serve as a first basis for simulating more complex flows at high Reynolds numbers.


Introduction
Performing a turbulent flow simulation is computationally costly due to the wide range of scales that it encompasses. With the vast majority of turbulent flows being wall-bounded, the present work aims to improve today's tools, unlocking the capabilities for simulating more complex flows. Even with the increasing power of today's large scale high performance computers (HPC), the computational cost remains a major criterion for high-fidelity, scale-resolving simulations.
Direct numerical simulation (DNS) of the Navier-Stokes equations is extremely costly and is only computationally affordable for moderate Reynolds number flows. Results from DNS will be used in this work as a reference, to assess the accuracy of the modelling presented. Most often, the quantities of interest are integral quantities (such as the aerodynamic forces exerted on an object), or the profiles of mean velocity and turbulent fluctuations, the mean turbulence intensity, the dynamic wake behaviour, etc. In contrast, the bulk of the computational cost in DNS lies in resolving the smallest scales that are not the main contributors.
As a remedy, large eddy simulation (LES) filters out the smallest scales, and models their effects on the large scale using a sub-grid scale (SGS) model. While LES is much less costly than DNS and hence allows simulating higher Reynolds number flows, it is still limited mainly by the cost of resolving the near-wall region. Modelling this region further reduces the cost significantly, and it makes it possible to simulate high Reynolds number flows. Previous research Choi and Moin (2012); Yang and Griffin (2021) suggest that the number of grid points for wall-resolved LES (wrLES) scales as Re 13∕7 , while it scales as Re 1 for wall-modelled LES (wmLES). The use of wall modelling is therefore mandatory to enable the simulation of high Reynolds number flows using LES.
WmLES reduces the computational cost of LES by modelling the inner layer ( y ≲ 0.15 where y is the wall-normal distance, and is the boundary layer thickness), and connecting it to the LES used for the rest of the flow domain. A simple way to picture this is to consider the region in the boundary layer (here assumed not near separation) and above the inner layer ( 0.15 ≲ y ≤ 1 ): this region is constrained by the pressure gradient driving the flow, and by the shear stress propagated from the wall. As such, the inner layer itself is not of much interest, and one can replace it with the wall shear stress that would have been produced if it was fully resolved. What is called the "model" is the computation of this shear stress replacement, and the "connection" to the resolved LES, is this shear stress imposed as a boundary condition.
Wall modelling for LES can be classified in two main types: the hybrid RANS-LES approach which computes this region in Reynolds Averaged Navier-Stokes (RANS) formulation, and the wall-stress models which uses an equation as model; typically an algebraic equation. Both types sample the velocity information from the LES at a distance from the wall as an input (in the resolved region), and they output the wall shear stress w , used as boundary condition for the LES. The distinction between the resolved scales from the LES and the modelled region is not made by a formal interface. It is here taken as the edge of the first element of the high order method used here (that uses a 7th order polynomial). The present work thus focuses on the implementation of a simple algebraic wall-stress model in a high order spectral element method. It will be seen that this framework also provides better results in the upper part of the under-resolved region when compared to lower order implementations, and even to other types of high-order discretization found in the literature; hence providing more accurate results at a lower cost.

3
The presented framework, with its investigation of channel flows, falls within the validity of the wall model used; and even more so for the high Reynolds number case of Re = 5200 that has a better separation between the inner scales and the outer scales. Wall modelled LES of more complex flows, such as in the presence of an adverse pressure gradient will have to be considered as next step. In such flow, the present simple algebraic model for wmLES can only be used if the influence of the pressure gradient does not persist close to the wall; and here more precisely close to the edge of the first element. This is indeed the case in high Reynolds number boundary layers as long as the complement function (i.e., the so-called "wake function" of Coles that adds to the logarithmic profile) has not reached the edge of the first element; see also the investigation by Frère et al. (2018). It is hence not expected that the present simple model will properly allow to simulate a boundary layer very near separation or at separation (recalling that the velocity profile at separation is solely the wake function, without law of the wall part). A more general analysis on wmLES for complex turbulent flows can also be found from Bose and Park (2018).
For reasons that will become clear to the reader, the code of choice used here is Nek5000, which is a high order spectral element code. Wall modelling for LES is well documented for widely used discretization methods such as finite differences, with available guidelines in Larsson et al. (2016), or as recently proposed for a discontinuous Galerkin method in Frère et al. (2017). These guidelines cannot be used as such in the present code, due to the large differences in the discretization methods, which has much influence on the wall model behaviour. A dedicated analysis is carried out for the present spectral element code, and specific guidelines can be established from these results. As the reader will notice from the results of this paper, the modelled region is highly sensitive to the discretization method, emphasizing the need for the present analysis. In parallel to the present paper, a similar implementation to wmLES has been performed by Argone National Laboratory with successful results in Pal et al. (2022). The authors include a thermal law-of-the wall, complementing well the present paper. The main difference in the present work compared to the above-mentioned publication is the use of a lower resolution (analysis on a less expensive computation) ad at much higher Reynolds numbers.
The code will be presented in some details in Sect. 2, also describing its advantages which will also lead to good results in part of the near-wall under-resolved element. Then the formulation of the algebraic model will be exposed, as well as the related constraints for its use. The used algorithm is iterative, and its performance will be assessed. An alternate method will also be proposed based on other research from Maheu et al. (2012). As mentioned previously, the ability of thge present wmLES to reproduce the results of a DNS will serve as comparison. For this, the test case of a channel flow at Re = 1000 will be described, as well as the convergence criteria, and the criteria for error quantification. In addition to those results, and in an attempt for generalization, a channel flow case at Re = 5200 will serve as a validation of the method at higher Reynolds numbers. The two main regions will be considered: the under-resolved inner part and the resolved outer part. A sensitivity analysis on the location of the boundary separating those parts will also be performed on the case at Re = 1000 , which constitutes the main guideline for creating a mesh for wmLES in this framework. Then a sensitivity analysis on the resolution of the outer layer of the same case will outline one main drawback of the spectral element method, while still providing a close comparison to the DNS. The comparison with the higher Reynolds number case will support the generalizability of the method. Finally, we will conclude with the guidelines on the placement of the boundary between the resolved region and the under-resolved first element, the numbers of degrees of freedom (DOFs) required, and a summary of the obtained results.

Numerical Methodology
The code used in the present work is the spectral element code Nek5000 developed by Argonne National Laboratory. This is an open-source code, providing users the license to modify the source code as needed. As any spectral code, it presents the advantage of computing very efficiently partial derivatives. Using a meshing method, the domain can be discretized in elements with different sizes according to the local accuracy needed. Nek5000 combines both of these characteristics by discretizing the domain in finite elements, and solving each element using a spectral method.
This spatial discretization is based on the spectral element method (SEM) formulated in Patera (1984). It is a high-order weighted residual technique using Gauss-Lobatto-Legendre polynomials (described in Deville et al. (2002)) on each element. The order of accuracy of the code is dictated by the order of the polynomials used as a basis function in the discretized elements. The main advantage of a higher order method is to produce smaller dispersion errors as depicted in Wang et al. (2013) for a similar cost compared to lower order methods.
More specifically related to Nek5000, previous work by Argonne Laboratory (2013) has demonstrated that the code exhibits minimal numerical dispersion and dissipation errors, taking into perspective the current performance of lower-order CFD codes. The solution and data are represented in terms of Nth-order tensor-product polynomials within each hexahedral elements, as described in Argonne Laboratory (2021).
The efficiency of the code comes from the above-mentioned factors, resulting in reduced cost per degree of freedom, and exponentially fast convergence. It is built as a highly parallelized code, and its impressive scaling capabilities have been demonstrated in Offermans et al. (2016).
The equations are, in the LES formulation (all field thus being LES fields, here written omitting to use any complex notation, for simplicity): with the fluid density, u the velocity, t the time, = ∇u + ∇u T the molecular shear stress tensor with the fluid dynamic viscosity, sgs the modelled subgrid-scale (SGS) stress (see later), p the pressure (including any no-modelled trace of the SGS tensor), and f a forcing (such as an added mean pressure gradient used to produce the flow in a channel, see later). This formulation is known as the stress formulation, and is required to prescribe the wall shear stress as boundary condition.
The boundary condition used on the wall to impose the shear stress is referred in Nek5000 as 'sh '. As it is used in the x direction, it is specified with the following code TRX = -(u_tau**2)*rho, with rho being defined from param(1) (referring to the .par file).
The SGS model is an implementation of an "approximate deconvolution model" (ADM), and for which the most influential parameter is a "relaxation parameter" that provides the It is user-defined using the variable filterWeight; a default value of 10 provided an excellent match with a DNS, as shown in Schlatter et al. (2004) at Re ≃ 210 , and is used here. The model also uses some high-pass filtering to build the SGS model; it is thus not a simple eddy-viscosity SGS model. The proposed study will use the 7th order polynomial, as previous research in Fischer (2013) indicates that this choice provides the best performance for a wide range of applications. Elements are loosely coupled, meaning that the velocity on the boundary nodes of an element is shared with the neighbouring element, making the velocity continuous between two elements (unlike the Discontinuous Galerkin method). The velocity gradient is however discontinuous at the interface between the two elements. Implications from this characteristic will become clear when examining the velocity profile or the turbulent fluctuation profiles as the grid becomes coarser: the discontinuities in the gradient at the element boundaries will create a Gibbs phenomenon, Fang et al. (2011), as manifested by spikes in the gradient or in the velocity fluctuations. As this phenomenon reduces to zero as the mesh converges to a DNS mesh, it can therefore be used as a qualitative measurement of the local mesh quality without the need of a reference solution (which is not available for most flows). For a quantitative assessment of this behaviour, a detailed analysis is given by Fang et al. (2011).
A clear distinction must of course be made between an element and a node, see Fig. 1.

Wall-Model Formulation
A wall-stress model samples the streamwise velocity u || i at each node i close to the wall, and it uses the model to deduce the corresponding wall shear stress w i , which is then used as a boundary condition for the LES, as shown in Fig. 2. The sampled velocity should be far enough from the wall to be in the resolved region, while being sufficiently close for the law of the wall to be valid ( y ≲ 0.15 ). More specifically, the sampled velocity is here taken at the last node of the first element. Moreover, the code parallelization is done element-wise; therefore the most efficient implementation is to impose the shear stress based of the velocity within a single element, so as to avoid global transfer of data between elements.
The model chosen to approximate the law of the wall is that of Reichardt (3), as also used in Frère et al. (2017): where u + || = u || u and y + = y u . We used the values = 0.38 and C = 4.1 as suggested in Osterlund et al. (2000). Note that, following Lee and Moser (2015), the most recent values are ≃ 0.383 and C ≃ 4.25 ; see also the analysis in Winckelmans and Duponcheel (2021). The law of the wall characterizes the flow near the wall, using the time-averaged streamwise velocity; whereas the LES inherently solves the fluctuating velocity field. This causes a mismatch at the interface between the wall layer and LES domain, which can be accommodated by averaging either in time or in space. The performance of the wall model using either of these two approaches has been compared. Spatial averaging in the spanwise direction element-wise has be tested, but it does not appear to make a significant improvement. Averaging over the full spanwise extent could also have been considered in view of reducing the time constant C t needed to obtain a similar averaging. However, the guidelines sought for are meant to be also applicable to inhomogeneous flows; hence it was not investigated. Time averaging was observed to be sufficient, and with a minimal effect on the computational cost.
The following moving averages were explored : • The Simple Moving Average (SMA) which is the average of the last few sampled points. It requires to define the averaging time and is not very reactive. • The Weighted Moving Average (WMA) which provides the possibility to apply more weight to the last values so as to increase the reactivity of the system. • The Exponential Moving Average (EMA) which is the averaging of all points and using an exponential weight, thus also favouring the last value. It is simple to implement and light on memory (no need to store other previous data) and is mathematically equivalent to a first order low-pass filter. The time constant is chosen in relation to a physical quantity. This choice also results in a reactive system.
One should understand the reactivity of the system here as the inverse of the time needed for the moving average to output a changed input. A steady flow will eventually converge to its steady state, but the transition period to reach steady state will be shorter with a higher reactivity, hence less computationally costly. For a wmLES, the EMA (4) appears to be an appropriate method for averaging the streamwise velocity, whereby: where i is the node index, n is the time-step index, = Δt∕(C t + Δt) with Δt the simulation time-step, and C t is the time constant of the filter. With this formulation, the time constant is chosen in relation to a physical timescale, such as C t ∝ u .

Description of the Wall Model Solver
As expressed previously, the wall model solver takes as input the sampled streamwise averaged velocity u || at a distance y from the wall, and it outputs the wall shear stress u . To do so, it attempts to solve the law of the wall, formulated by the Reichardt's equation (3). As there is no direct analytical solution to find u , it is common practice to use a numerical method as in Liefvendahl et al. (2017). To solve it numerically, the equation is rearranged in the form of f (u ) in equation (5), where the root of the function is the solution to u : Most published implementations use the Newton-Raphson method, which start from one initial guess and uses the slope of the curve to reach the root iteratively. Alternatively, as suggested in Maheu et al. (2012), the process to find the value of wall can be performed using a lookup table rather than an iterative process. This method was used further in Frère et al. (2018) by using the quasi-linear relationship between u + y + and y+ to use the lookup table efficiently. Nevertheless, the presented method uses an iterative process based on the convergence of two initial guesses, with indistinguishable performance hit on the total time/iteration. For this channel flow test case, each node of the wall is calling this algorithm independently, converging to machine precision in seven steps on average. The cumulative time for this process on all near wall nodes sums to only 0.4% of the total time required for a timestep, which is imperceptible. With such performance, there is no need for a lookup table method. Additionally, a lookup table is only accurate up to its pre-computed resolution, whereas the proposed method reaches machine precision.

Time Averaging Requirements
Averaging in time such that the velocity profile appears to be sufficiently converged statistically can take a long computational time, and would slow down the ability of the model to self-adapt to the steady state convergence. It is of utmost importance to reduce the time constant C t used for averaging. Previous studies such as Yang et al. (2017);Fowler et al. (2022) suggest that this time averaging helps significantly the computation. This present paper investigates it in the framework described in previous sections.
In case of a steady flow, the effect of this time constant can only be seen in the settling time; which is the time to reach steady state from the initial conditions. On the other hand, when simulating an unsteady flow, such a delay should be limited as the shear (4) u i n = u i n−1 + u i n − u i n−1 , stress boundary condition will need time for the averaged sampled velocity to adapt to the unsteadiness of the flow. To investigate this, several cases have been explored with C t ranging from 0 to a large value of 2 u in the case of Re = 1000 . The case at C t = 2 u provides good agreement with the DNS, but at the cost of a long settling time. A more detailed comparison between the wmLES and the DNS will be explored in the following sections. Without averaging ( C t = 0 ), the law of the wall cannot handle the high fluctuations of the non-averaged streamwise LES velocity to determine the wall shear stress, leading to a divergence. This is most likely due to the Gibbs phenomenon at the element boundary creating high fluctuations in the sampled velocity, hence affecting the convergence of the algorithm.
It appears that the smallest time constant C t = 0.0025 u provides the best compromise. Such small time constant keeps most of the turbulent quantities unaffected, and only filters out the highest frequency content, which is most likely due to the Gibbs phenomenon on the derivatives. The scales of motion with higher frequency than 1 C t are therefore filtered out. The final solution is statistically equivalent to the case with C t = 2 u , validating the hypothesis of this independence for a steady flow. The time constant is then equivalent to 5 time-steps ( 5 Δt ). As we will see in Sect. 9, this is persistent at higher Reynolds Re = 5200 , at which the same value of C t is found to be near optimal.
A qualitative comparison is given between the instantaneous and averaged flows in Fig. 3, as well as for the sampled plane at y + = 153 in Fig. 4 where the averaged velocity is used as input to the wall model. The time constant used here is C t = 0.0025 u . We see that the difference between the instantaneous and averaged velocity field is small. This implies that the law of the wall is quite robust to the input quantity, yet it still requires a minimal time averaging. Notice that since the averaging time required can remain very short, this also opens the possibility of applying the present methodology to unsteady flows. In the present case, the time averaging has no impact on the final solution.

Test Case at Re = 1000
A periodic channel flow is simulated in order to evaluate the new wall model implemented in Nek5000, with a domain size of 2 streamwise, spanwise, and 2 in the wall-normal direction, with being the half channel height. Periodicity is applied in the two homogeneous directions (streamwise and spanwise). The flow is driven by a constant flow rate imposed in the streamwise direction. The Reynolds number based on and u is Re = 1000 . The wmLES results are compared to those of the DNS database from Lee and Moser (2015).
As introduced in Sect. 2, the mesh quality influences the results in a very distinguishable way. With a coarse mesh, the velocity fluctuations exhibit a spike at the elements boundaries (Gibbs phenomenon Fang et al. (2011)) due to the spectral element method. This provides a natural framework for the refinement of the mesh. Consequently, the mesh in the streamwise and spanwise directions has been refined such that the Gibbs phenomenon becomes small compared to the local variations within each element. This results in 22 and 16 elements in the streamwise and spanwise direction respectively. With the dimensions of the domain, this leads to 0.286 per element streamwise, and 0.196 per element spanwise. With the 7th order polynomial order used, it leads to an average of 0.040 and 0.028 between the nodes, respectively. So, the resolution is about 24.5 and 35.7 nodes per , respectively. The mesh in the wall-normal direction is not given here, as it will be explored in more details in the following sections with a parametric study. This framework to refine the mesh provides similar results to the guidelines from Larsson et al. (2016), validating the approach.
The initial conditions are composed of a mean streamwise velocity profile equal to the Reichardt equation and perturbed using sine waves in the three directions to initiate turbulence generation.
The results shown hereafter are first averaged in space (in the two homogeneous directions). Within each element, a weighted average based on the weight of each node is used, while each element is averaged without weights. Additionally, the averaging procedure uses the channel's geometrical symmetry to accelerate statistical convergence by averaging between the top and the bottom half channels. Averaging is then also done in time after the transition period has passed, until a statistical convergence of at least 99.5% is achieved for all quantities of interest. The criteria used to quantify the statistical convergence is the average of the values contained in the profiles of u ′ u ′ , v ′ v ′ and w ′ w ′ . The values are shown in Fig. 5 for one of the following test case, with an exponential fit exposing the convergence; the same methodology being applied to all other cases. In order to compare the proposed results with the reference DNS, an error is defined. The normalized root-mean-square error (NRMSE) formulated in Eq. (6) will be used to compare the results obtained by wmLES to the reference DNS: with the DNS data being linearly interpolated to the location of the wmLES nodes. As the DNS data are of high resolution, a linear interpolation is sufficiently accurate.
Quantifying the match of Re between the wmLES and DNS is given in an attempt for transparency to the reader. For this, the error on the shear stress will be reported as a percentage. This is computed using the dimensionless ratio u ∕u mean , where u mean is the mean flow velocity across the channel. u is averaged across both walls, both homogeneous directions and in time for ensure proper convergence. The percentage is computed as the difference of this parameter between the wmLES and the DNS, divided by the DNS.

Sensitivity to the First Element Height
The boundary between the resolved and under-resolved regions is non-trivial in most cases. In the current framework, the mesh resolution above the first element is sufficiently resolved to capture all the scales of motion, and are therefore considered as a resolved region. The first element adjacent to the wall is unresolved near the wall (it does not satisfy the no slip condition at the wall), and it must be considered unresolved for the entire element (as the solution is computed with one single polynomial for the entire element). This creates a distinctive boundary between the resolved and unresolved regions on the domain. One could further investigate what would happen if the boundary between the resolved and unresolved regions was placed somewhere within the upper part of the first element rather than at its upper edge; this could be investigated in future work.
The height of this first element provides an adjustable size of under-resolved region, as defined by the mesh. The input velocity is sampled at the end of the first element, which is an efficient method due to the way the solution is stored numerically. This section will focus on the effect that this size parameter has on the final solution.
Three simulations were performed at Re = 1000 , varying the size of this first element, while maintaining six elements within the half channel height . The mesh resolution of the outer (resolved) layer has a negligible impact on the solution as it is fully resolved. Starting from a uniform distribution y u ∕ , the following stretching law (also used in Duponcheel et al. (2014)) is applied in the wall-normal direction: Figure 6 illustrates the resulting wall-normal distribution obtained for = 1.5 , 1.0, and 0.4, corresponding to y = 0.063 , 0.10 and 0.15 for the height of the first element respectively; that is y + = 63 , 104 and 153 when expressed in wall units. Note that they are all located within the law of the wall region.
The results from this setup are compared to the reference DNS in Fig. 7 for the mean velocity profile and in Fig. 8 for the mean fluctuation profiles. For transparency to the reader, all the under-resolved nodes within the first element are shown, even though those in the lower part of the element are expected to behave poorly (recall that the no slip condition is not enforced at the wall; what is enforced is the wall shear stress). In general, the wmLES results are found to be in good agreement with the DNS both for the mean velocity and mean fluctuation profiles.
In Fig. 7, the individual nodes within each element are shown; even those in the underresolved region (except the node at y + = 0 due to the logarithmic scale). It is expected for a wmLES to provide accurate results in the resolved region, while accepting that it performs poorly in the under-resolved zone. The closeness between the wmLES and DNS confirms the validity of the implemented method. The reader is invited to focus on the boundary between the resolved and under-resolved regions, at y + = 63 , 104, and 153. It appears that even inside the under-resolved zones (below those three thresholds), the mean profile is still very well captured down to y + ≃ 10 ; only the first two nodes differ significantly.  This is quite remarkable; given that, in the literature, the wmLES results are not shown in the under-resolved region; or shown but dismissed. Due to the high-order character of the spectral method, a good part of the first element is still able to behave in a way that is similar to the resolved case, at a fraction of the computational cost.
Qualitatively, the turbulent fluctuations ( u ′ u ′ , v ′ v ′ , w ′ w ′ ) show negligible spikes at element boundaries, which means that the mesh is here refined enough to accurately capture the mean and turbulent flow fields.
In general, the agreement between the wmLES and the DNS is quite good, and for all three cases. Similarly as for the mean velocity profile, the fluctuating quantities are expected to provide inaccurate results in the under-resolved region, and perform well in the resolved zone. Here again, the individual nodes within each element are shown. A good part of the unresolved region appears to also capture the fluctuating quantities. The part that differs significantly from the DNS is that of the first three nodes of the first element, where both the streamwise ( u ′ u ′ ) and spanwise ( w ′ w ′ ) fluctuations do not settle to 0 at the wall. The fluctuations in the wall-normal direction appear to match entirely the DNS in the under-resolved region. This is due to the boundary condition restricting the normal velocity to 0 at the boundary.
The difference between each wmLES cases is imperceptible, which means that the first element height does not appear to have a significant impact on the solution (as long as it remains within the region where the law of the wall is valid, which was the case here). As a quantitative analysis, we use the NRMSE as a measure of the total error between the wmLES and the DNS. The under-resolved region should be omitted from this calculation, as is it not a region of interest since this region is commonly accepted to be inaccurate. Table 1 provides the NMRSE for these three cases, and only excluding the first 3 nodes of the first element. These results show small errors both in the mean velocity ( ≲ 1% ) and in the mean fluctuations ( < 8% ) profiles. A trend is visible with errors increasing slightly with the first element height, and mostly on the mean velocity and u ′ u ′ profiles; suggesting that using a first element height comparable to the end of the law of the wall region (here y + ≃ 150 ) is not necessarily the best choice. The errors are nevertheless quite acceptable in all three cases. In addition to the error on the profiles ( u , u ′ u ′ , v ′ v ′ , w ′ w ′ ), the relative error on the wall shear stress between the wmLES and the DNS is quantified with the dimensionless ratio u ∕u mean .
Given the relative insensitivity of the results with respect to this first element height (as long as it remains within the law of the wall region), the wall model used here constitutes a quite flexible approach.

Sensitivity to the Outer Region Mesh
The computational cost and accuracy of the results due to the mesh size in the inner layer y ≲ 0.15 was studied in the previous section at Re = 1000 . The focus is now on determining how sensitive the results are relatively to the mesh resolution used in the outer layer y ≳ 0.15 . Based on the previous section's results, it appears that the results are insensitive to the first element height (within the validity of the law of the wall). To reduce the   . 9 Mesh for the sensitivity analysis to the number of DOFs used computational cost, the first element height is therefore set to the largest value of y = 0.15 limited by the validity of the law of the wall, and is kept constant for all cases. Three meshes are explored, as illustrated in Fig. 9, where the half channel contains 6, 5, or 4 elements in the wall-normal direction. A more consistent way to evaluate the computational cost of the cases is to use the degrees of freedom (DOFs, number of nodes on which the solution is computed). Using the 7th order polynomial, the three cases result in 42, 35 and 28 DOFs in the wall-normal direction for the half channel.
The mean velocity profile given in Fig. 10 shows the individual nodes, as in the previous section. The overall shape of the profile is properly captured by all three mesh sizes. Regarding the nodes within the first element, they again appear to provide results close to the DNS, even though they are under-resolved. This demonstrates that the mean profile is well captured, even with the coarser mesh in the outer layer.
Regarding the mean fluctuation profiles presented in Fig. 11, all three cases overall appear to match closely the DNS. Some differences are seen between the three mesh sizes, as the coarser one is not able to capture the physical behaviour as well as the finer meshes. The larger elements size in the bulk of the domain result inevitably in the introduction of a Gibbs phenomenon at the element boundaries, especially visible for the 4 elements case. As introduced previously, this phenomenon gives an indication of the mesh quality, and can be used to easily analyse qualitatively if the mesh needs to be refined. We notice a spike in the fluctuations (Fig. 11) in the centre of the channel, for all three cases, irrespective of the size of the element. This is likely due to the fact that the channel centre corresponds to the edge of two elements. The mesh in the bulk of the domain should be further refined there, where the solver must capture a change in the gradient (between the bottom and top halves of the channel). Indeed, as the location of this change of gradient (from positive to negative) is varying over time yet remaining near the centre, it might be difficult for the solver to capture it properly with the present mesh configuration (recall that the present method only enforces the continuity of the velocity fields at the interface between two elements; it does not enforce the continuity on the derivatives). It would also be interesting to run a case where the centre of the channel is fully covered by an element; this is left to future investigations. While the under-resolved region spans from y + = 0 to y + ≃ 150 , this region is still shown in the graphs and the error in the under-resolved region is again fairly low for such mwLES.
For a more quantitative comparison, the NRMSE error between the DNS and wmLES is reported in Table 2. As in the previous section, those errors exclude the first three nodes of the first element closest to the wall. The NRMSE are sometimes lower for coarser resolution, but this is due to the fact that, in all cases, the mean fluctuations are under-estimated by the wmLES, while the averaged value is artificially improved thanks to the amplified spikes. The case with 4 elements in the half channel height must therefore be considered as not sufficiently resolved to provide an accurate solution to the problem. The case with 5 elements can be considered as acceptable. Here again the relative error on the wall shear stress between the wmLES and the DNS is quantified using the dimensionless ratio u ∕u mean .

Case at a Higher Reynolds Number
This section showcases the results of the above-mentioned wall model at a significantly higher Reynolds number, as a further assessment of the newly implemented wall model. The highest Reynolds number DNS database currently available and of high quality is that Re = 5200 , by Lee and Moser (2015), which is used as reference. The goal here is to use the guidelines extracted from the previous sections at Re = 1000 , and push the limit of the model to a much more challenging case to highlight that good accuracy can still be obtained, even on a very coarse mesh.
To this end, a mesh with acceptable Gibbs phenomenon is chosen. The case from Re = 1000 with 4 elements in the half-channel height has noticeable Gibbs spikes; yet it remains a fair compromise between accuracy and computational cost. For the higher Re = 5200 , we first investigate a case with uniform resolution and using 20 elements in the half channel height: the first element thus ends at y∕ = 0.05 , which corresponds to y + = 260 thus already within the log layer. A case with 10 elements will be considered in second (see later). Using the same outer domain dimensions as the previous case ( 2 streamwise, spanwise, and 2 in the wall-normal direction), the resolution is doubled relatively to the case at Re = 1000 : thus 44 and 32 elements streamwise and spanwise, respectively. This is also done so as to mitigate the Gibbs phenomenon on the mean fluctuation diagnostics.
The time averaging needed was once again explored at this higher Reynolds number, and the same value of C t = 0.0025 u is found to be near optimal; which supports the generalizability of that choice.
The mean velocity profile is provided in Fig. 12. As previously, all nodes are shown (also the under-resolved nodes of the first element, the node closest to the wall still being far away from the laminar sublayer).
The mean profile is well captured, and even in the under-resolved region of the first element. This is due to the fairly high resolution of this wmLES; even the value at the first node, that lies in the transition region, is still close to the reference. The velocity fluctuations are shown in Fig. 13. The wmLES results match with good accuracy the DNS reference, However, there is still some Gibbs phenomenon (mostly visible close to the wall, and especially for w ′ w ′ ). This illustrates again a main drawback of such high order spectral element method: since the fields are represented by a high order polynomial in each element, and since only the continuity on the field is enforced between elements (and not on its derivatives), the statistics on the fluctuations of any LES field will display a Gibbs phenomenon. The only way to reduce that is to increase the resolution.
The fluctuations in the first under-resolved element are also provided, for completeness. They are of course not expected to be accurate.
The differences between the wmLES and the DNS are reported in Table 3, using the NMRSE as for the previous cases. The first element was omitted from the calculation as it would bias the comparison. The relative error on the shear stress is also provided in the table.   As a second test, we investigate a case using 10 elements in the half channel height: the first element thus ends at y∕ = 0.10 , which corresponds to y + = 520 ; well within the log layer. The resolution is also taken back to 22 and 16 elements streamwise and spanwise, as for the wmLES at Re = 1000.
The mean velocity profile is also provided in Fig. 12. The mean profile appears well captured in the resolved region; and not so well in the under-resolved region of the first element. The velocity fluctuations are provided in Fig. 13. As expected, the Gibbs phenomenon is more important. Apart from that, the u ′ u ′ and v ′ v ′ fluctuations are similar to those obtained at the twice finer resolution; larger differences are however seen on the w ′ w ′ fluctuations.
The present results at the higher Reynolds number of Re = 5200 indicate that the present wmLES framework can indeed be used for a wide range of applications. The results show an error of < 1% on the mean velocity profile and < 9% on the velocity fluctuations; which are similar to what was obtained for the cases at Re = 1000.

Conclusion
In this paper, the capabilities of combining wall-modelled LES (wmLES) with the high order spectral element code Nek5000 are presented, reducing significantly the computational cost. An algebraic wall-stress model is implemented in Nek5000, and assessed with using the results of reference direct numerical simulation databases on channel flows: first at Re = 1000 , and then at Re = 5200 . The velocity in the near-wall resolved region is used to solve a "law of the wall" model for the wall shear stress, and this stress is then imposed at the lower edge of the under-resolved near wall element. It is solved numerically using an efficient iterative process which has an imperceptible impact on the computational time of each time-step. Averaging in time has also been investigated. It is seen that it can be significantly reduced; this opens the way to using the framework to perform wmLES of more complex, highly unsteady, flows.
The boundary between the resolved and under-resolved regions is defined by the height of the first element. It was shown that this parameter does not significantly affect the solution, which makes this a quite flexible approach (as long as it is placed in the region where the law of the wall is valid: y ≲ 0.15). Regarding the outer region mesh requirements, the solution behaves similarly to a wrLES. No further limitations on the wall model are introduced. The resolved region compares well to the reference DNS with errors on the mean velocity profile below 1%, and errors on the mean fluctuation profiles of 7-9% for the worst case, for the two Reynolds numbers tested. The results however highlight a drawback of such high order spectral element method: since the fields are represented using a high order polynomial in each element, and since only the continuity on the field is enforced between elements (and not on its derivatives), the statistics on the fluctuations of the LES fields display a Gibbs phenomenon.
Another interesting observation is that the wmLES results for the mean profile in the under-resolved first element do not deviate much from the reference DNS results; and even down to the first node when the LES resolution is such that the first node lies in the transition layer. This good behaviour is likely due to the spectral element method (which has already proved to be very accurate and efficient in DNS). Most of the deviations, in comparison to the DNS, are observed on the fluctuating components, and of course even more so in the under-resolved first element.
The reduced cost of LES with the use of this simple wall model will enable simulations of high Reynolds flows, and it provides promising perspectives for the use of wmLES to simulate more complex flows, and for a potential industrialization of the approach.