The cost of LES in the near-wall region of attached turbulent boundary layers scales as the square of the friction Reynolds number thus limiting LES to moderate Reynolds numbers. Wall stress boundary conditions are frequently used to alleviate this resolution requirement, but commonly used models are shown to perform poorly at high Reynolds numbers even in turbulent channel flow. Techniques from optimal control theory are used to find wall stresses that yield much better results in turbulent channel flow at high Reynolds numbers than existing models even on extremely coarse grids. In this approach, an optimal control strategy is used in which the objective is to force the outer LES towards a desired solution by using the wall stress boundary conditions as control. The optimal wall stresses are not necessarily physical, rather they are whatever is necessary to overcome the numerical and modeling errors present in the near-wall region to yield the correct mean velocity profile. Furthermore, the optimal control strategy generates reference data for comparing and deriving new wall models. Using linear stochastic estimation it is shown that the dynamically relevant part of the optimal wall stresses can be predicted from the local velocity field. A wall model derived from linear stochastic estimation yields good mean flow predictions in LES of turbulent channel flow on a 32 3 uniform grid for friction velocity Reynolds numbers from 640 to 20000.

Introduction

The large number of grid points required in the near-wall region of attached turbulent boundary layers is the chief obstacle to applying large eddy simulation (LES) to many flows of engineering interest. Current subgrid scale models do not accurately model the subgrid scale (SGS) Reynolds stresses. [START_REF] Baggett | Resolution requirements in large-eddy simulations of shear flows[END_REF][START_REF] Jiménez | LES: Where we are and what we can expect[END_REF] Thus, if the LES is to include the near-wall region, the filter width has to be sufficiently small so that most of the Reynolds stresses are carried by resolved motions. This requires the filter width to scale as a fixed fraction of the local turbulent integral scales which are proportional to the distance from the wall. Baggett, et al. [START_REF] Baggett | Resolution requirements in large-eddy simulations of shear flows[END_REF] calculated that the number of grid points required for accurate LES of a turbulent boundary layer with a resolved near-wall region scales as N ∼ Re 2 τ , where Re τ is the friction Reynolds number.

To alleviate these prohibitive near-wall resolution requirements, the no-slip boundary condition is usually replaced by approximate boundary conditions. Approximate boundary conditions attempt to account for the effects of the unresolved near-wall region on the outer flow. Deardorff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF] was the first, in 1970, to employ this approach in simulations of turbulent channel flow at infinite Reynolds number. In the most common approach, the wall stresses are modeled and the transpiration velocity is set to zero. The latter approach was first introduced by Schumann, [START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF] in 1975, who assumed the streamwise wall stress was in phase with the streamwise velocity at the first off-wall grid point. Many variations and improvements on this basic idea exist and will be reviewed in the next section.

Nearly all existing wall stress models are designed to reproduce first order statistics, in the spirit of wall functions for Reynolds averaged Navier-Stokes simulations (RANS). Typically, wall stress models attempt to reproduce the logarithmic law of the wall in the mean. Fluctuations about the mean are usually modeled in an ad hoc manner. Even models for the wall stresses that attempt to incorporate more physics, such as the two-layer model to be discussed in the next section, while sometimes improving the prediction of the mean wall stress do not necessarily yield wall stress fluctuations that will improve the prediction of the outer LES (see Figure 2 in Section 2). In fact, it is difficult to know a priori what is required of a wall stress model to ensure the accuracy of the outer LES.

The difficulty in specifying good approximate boundary conditions stems from the fact that in the near-wall region the filter width (i.e. the resolution) is larger than a significant fraction of the energy containing eddies that would be present in an actual flow. The usual LES approximation that the subgrid scale model has to account only for the more isotropic small scales is invalid at the filter widths that are typically maintained in the near-wall region. In fact, in LES of a turbulent attached boundary layer in which the first off-wall grid point is in the logarithmic region, the filter width will always be larger than the local turbulent integral scales which are proportional to the distance from the wall. Thus, subgrid scale modeling errors will dominate the near-wall region. The problem of wall stress modeling is inherently coupled to the subgrid scale model used and the numerical resolution. In this work, we propose a strategy for finding wall stress models that account for this crucial linkage.

Since there is no known subgrid scale model for the near-wall region when the filter width is larger than the turbulent integral scales, the information that must be supplied by a wall stress model to the outer LES is essentially unknown. To overcome this difficulty, we use an optimal control strategy in which the objective is to force the outer LES towards a desired solution by using the wall stress boundary conditions as control. This strategy does not result in a practical wall model for LES because the solution must already be known, but it does produce reference data that can be used to compare or derive new wall models. In fact, we use the data from an optimally controlled simulation of turbulent channel flow to derive a practical wall stress model as will be shown below.

In the next section, we present a brief review of existing wall stress models and also mention some other strategies that have been proposed for alleviating the near-wall resolution requirements. In Section 3 the details of the numerical scheme for the channel flow LES are presented. The optimal control strategy is described in Section 4 and results are presented in Section 5. The reference data from the optimally controlled simulation is used to derive a practical wall model via Linear Stochastic Estimation in Section 6 where the results from the new model are also presented. We conclude with some remarks in Section 7.

Wall Stress Models

Nearly all attempts to account for the effects of the near-wall region have focused on modifying the no-slip boundary conditions. This is usually done by replacing the no-slip boundary with conditions on the wall stresses, that is, on the gradients of the tangential velocities at the wall. Typically the wall-normal velocity is set to zero. There have also been attempts to supply the velocities on some plane parallel to the wall, so that only the portion of the flow in which the LES can be trusted is computed, but these approaches have met with limited success. [START_REF] Baggett | Some modeling requirements for wall models in large eddy simulation[END_REF][START_REF] Jiménez | Approximate lateral boundary conditions for turbulent simulations[END_REF][START_REF] Nicoud | Boundary condition for les away from the wall[END_REF] The first wall model or approximate boundary condition was introduced by Deardorff, 3 who simulated a channel flow at infinite Reynolds number. He provided constraints on the wall-normal second derivatives of the wall-parallel velocities that ensured the logarithmic law was satisfied on average. The normal velocity at the wall was equal to zero. Unlike most subsequent wall models, Deardorff's model did not specify the wall stresses directly.

The first wall stress model is due to Schumann, 4 who used a staggered grid system to simulate finite Reynolds number channel flow. In this formulation, the wall stresses and the normal velocity at the wall are required (see Figure 1). Schumann assumed that the streamwise (spanwise) wall stress is in phase with the streamwise (spanwise) velocity at the first off-wall grid point. The constant of proportionality was obtained from the logarithmic law by requiring the mean wall stress to balance the mean pressure gradient which therefore has to be known a priori. The normal velocity at the wall was again set to zero as it is in all existing wall stress models for staggered grid systems.

A number of improvements to the basic idea of Schumann, that the wall stress is a simple deterministic function of the velocity at the first off-wall point, have been proposed by Grötzbach, [START_REF] Grötzbach | Direct numerical and large eddy simulation of turbulent channel flows[END_REF] Piomelli, et al, [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] Werner and Wengle 10 and Hoffmann and Bennoci. [START_REF] Hoffmann | Approximate wall boundary conditions for large eddy simulations[END_REF] Wall stress boundary conditions are also frequently used for LES of atmospheric boundary layers. For instance, Mason and Callen 12 enforce an instantaneous logarithmic law in each grid volume adjacent to the wall.

The above cited models have in common that a simple algebraic relationship is assumed between the wall stress and the velocity at the first off-wall points. In this paper, we take the shifted model of Piomelli, et al 9 as representative of these algebraic models. It provides the approximate boundary condition as:

τ w 12 (x, z) = u(x + ∆ s , y 1 , z) u τ w V (x, 0, z) = 0 (1)
τ w 32 (x, z) = w(x + ∆ s , y 1 , z) u τ w
where ∆ s is an empirical displacement, related to the mean orientation of the vortical structures near the wall. Moreover • denotes the plane average and the mean shear stress in the streamwise direction < τ w > is computed assuming a log law. We have found little difference in the results when using this model or some of the others mentioned in the previous paragraph at the high Reynolds numbers considered in this work.

More recently, a two-layer approach has been employed in which the three dimensional unsteady boundary layer equations are integrated on an embedded near-wall grid to estimate the wall stresses. [START_REF] Balaras | Two-layer approximate boundary conditions for large-eddy simulations[END_REF][START_REF] Cabot | Large-eddy simulations with wall models[END_REF] The two tangential components of the velocity are computed first in the inner region with no-slip boundary condition at the wall, the outer Dirichlet condition being provided by the LES computation in the outer region. This approach is still cost effective as compared with a fully resolved LES because the normal component of the velocity in the inner region is computed as a consequence of the divergence-free constraint. Thus the pressure field in the inner region is not computed through a Poisson solver but imposed from the outer region. Consequently the normal velocities in the inner and outer regions do not match at the interface between the two layers. This approach incorporates more physics in the wall model and it has met some success in the computation of the flow behind a backward-facing step. [START_REF] Cabot | Near-wall models in large eddy simulations of flow behind a backward-facing step[END_REF] However, as illustrated in the next paragraph, it does not produce better results than simpler algebraic wall models for coarse LES.

Figure 2 shows typical mean velocity profiles from LES of a periodic channel flow with the shifted model of Piomelli, et al. [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] and the two-layer model provided by Cabot. [START_REF] Cabot | Large-eddy simulations with wall models[END_REF] The simulation shown in Figure 2a is for a friction Reynolds number 640. The channel dimension is (2π h, 2h, 2π h/3) where h is the channel half height. The grid is uniform with 32x33x32 cells -see sections 3 and 5 for more details on the numerical method and the simulation. The mean wall shear stress is fixed so that the first velocity point matches the logarithmic law-of-the-wall shown in the figure . Clearly the algebraic model succeeds rather well in producing the shear stresses needed to recover the expected log law (the results from the two-layer model are very similar). Figure 2b shows however that the with the two-layer model. The grid is uniform with 32x33x32 points.

: < u > + = 2.41 ln y + + 5.2 same model is not as accurate when the Reynolds number is higher, say Re τ = 4000 instead of 640 (same domain, same grid). Also shown in Figure 2b is the result of a simulation with the two-layer model. [START_REF] Cabot | Large-eddy simulations with wall models[END_REF] The mean gradient is correct after the first few grid points but an artificial boundary layer develops near the wall and the log law intercept is under-predicted. This simulation illustrates the need for improved wall models that can cope with the high Reynolds numbers and necessarily coarse grids typical of many engineering simulations. In the following, we will compare our results to the shifted model (1) since it gives results comparable to the two-layer model at high Reynolds numbers.

The wall stress models mentioned above all attempt to produce the correct mean wall stress and often include ad hoc terms to account for the wall stress fluctuations. It is also possible to use data from a simulation in which the near-wall region is resolved to find a model which reproduces the wall stresses corresponding to a given velocity field. In an attempt to derive such a model, Bagwell [START_REF] Bagwell | Improved approximation of wall shear stress boundary conditions for large eddy simulation[END_REF][START_REF] Bagwell | Stochastic estimation of near wall closure in turbulence models[END_REF] approximated the conditional average of the wall stresses given the velocities on a plane above the wall by a Linear Stochastic Estimate (LSE) of the form:

τ w i2 = L ij * u j (•, y, •)
, where i = 1, 3, and j = 1, 2, 3 -see section 6 for further details. The convolution coefficients L ij were estimated using data from a direct numerical simulation (DNS) of turbulent channel flow at low Reynolds number, Re τ = 180. [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] The resulting model worked well as a wall stress boundary condition in a Re τ = 180 channel simulation where both the mean flow and turbulence intensities were well predicted. However, extrapolating the convolution coefficients to higher Reynolds numbers where wall models are required is difficult. Bagwell 17 only had some success in re-scaling the correlation data to Reynolds number 640. In this case the linear stochastic model performed marginally better than the above shifted model. A possible reason for the limited success met by Bagwell is that the linkage between numerical and subgrid scale model errors and wall model was not recognized. In the following we use an approach similar to Bagwell's one [START_REF] Bagwell | Stochastic estimation of near wall closure in turbulence models[END_REF] except for the way the reference data is generated (see section 4).

Providing approximate boundary conditions to the outer LES is not the only way to account for effects of the near-wall region. An alternative is to refine the grid in the wall-normal direction so that the steep wall-normal gradients typical of boundary layers can be captured on the grid, at least for the mean flow. Of course, if the grid is not refined in the directions parallel to the wall, then the correct near-wall physics cannot be represented on the grid and the subgrid scale model must be modified in the near-wall region to account for the unresolved motions. The simplest approach is to augment the SGS model with a RANS eddy viscosity model in the vicinity of the wall. This idea was originally proposed by Schumann 4 who used a mixing length eddy viscosity to supplement the Smagorinsky eddy viscosity in the near-wall region. Moin and Kim [START_REF] Moin | Numerical investigations of turbulent channel flow[END_REF] and Sullivan et al. [START_REF] Sullivan | A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows[END_REF] have explored similar approaches. More recently, Spalart et al. [START_REF] Spalart | Comments on the feasibility of les for wings, and on a hybrid rans/les approach[END_REF] have advocated the use of the one-equation Spalart-Allmaras eddy viscosity model as an SGS model. The length scale in the destruction term is modified so that the eddy viscosity crosses over from the usual Spalart-Allmaras RANS eddy viscosity near the wall to a proposed LES eddy viscosity, similar to that of Smagorinsky, away from the wall. Spalart, et al. call this approach "Detached-Eddy Simulation" since it is intended to be used in regions, such as separated regions, in which only eddies that are detached from the surface must be resolved for accurate LES of the flow away from the wall.

In what follows we will use only wall stress boundary conditions to account for the unresolved motions in the near-wall region. In the next section we outline the numerical method used for this study before giving the details of our approach.

Numerical Method

In this study a second-order accurate finite difference scheme is used to discretize the LES equations on a staggered grid system. [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF] Given the simple geometry considered (periodic channel flow, see Section 5), more accurate (spectral) methods could have been used. However, these highly accurate methods are not flexible enough to handle industrial applications with complex geometries (e.g. flow around an airfoil), where both low-order numerics for simplicity and wall modeling for high-Reynolds number boundary layers are needed. A staggered grid system is used to avoid the decoupled pressure-velocity mode as well as the prescription of a boundary condition for the pressure. The time integration is a third-order Runge-Kutta scheme for all the convection and diffusive terms. The diffusive terms in the normal direction to the wall are not treated implicitly since the grid is uniform in that direction. Periodic conditions are imposed in the two directions x 1 and x 3 (or x and z) parallel to the walls so that the Poisson equation can be solved efficiently using a FFT-based Poisson solver.

The subgrid scale model is the Smagorinsky model with the coefficient determined by the planeaveraged dynamic procedure. [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF] Unless otherwise stated, all quantities are nondimensionalized by the friction velocity, u τ , and channel half-height, h. The channel walls are at y = ±1. The skin friction Reynolds number is then defined as Re τ = u τ h/ν. When the mean flow is converged in the statistical sense, the mean streamwise pressure gradient is equal to the mean wall stress, that is,

-∂ p /∂x = τ w = 1.
Since the first grid point is outside the viscous layer in this study, the classical no-slip boundary condition for the velocity components is replaced by a set of approximate boundary conditions. More precisely, the transpiration velocity is set to zero and the two shear stresses τ w 12 and τ w 32 are provided to the momentum equation in the x 1 -and x 3 -direction respectively. The sketch in Figure 1 shows the location of the variables and boundary conditions in the staggered grid system. The wall normal direction is x 2 (or y) while u i (or u,v,w) and p are the velocity components and the pressure. Prescribing zero transpiration velocity is not fully justified theoretically but this choice has been made in view of the difficulties met in previous studies when this assumption is not made. [START_REF] Baggett | Some modeling requirements for wall models in large eddy simulation[END_REF][START_REF] Jiménez | Approximate lateral boundary conditions for turbulent simulations[END_REF][START_REF] Nicoud | Boundary condition for les away from the wall[END_REF] Whether or not the optimal framework used in this study can be extended to the more general case with non zero transpiration velocity will be considered elsewhere. The shear stresses τ w 12 and τ w 32 are prescribed either from one of the wall stress models discussed in Section 2 or the optimal control strategy.

Optimal formulation

In the process of deriving a new wall model for LES, it is crucial to keep in mind that:

1. The objective is to produce approximate boundary conditions able to handle very large Reynolds numbers. In this respect, using DNS data as a guide may not be the most judicious choice since these data are only available for low to moderate Reynolds numbers.

2. In any coarse grid LES where the first grid point is within the logarithmic region, the turbulent integral length scale (L ≡ κy) is less than half the grid spacing (∆y = y). As a consequence both subgrid modeling and numerical errors are important. Therefore the approximate boundary condition should compensate for these errors if the correct mean profile is to be obtained in a coarse grid computation. In this case the best shear stresses to supply are not the physical ones, as from a DNS.

It follows that the reference data used to compare or derive new wall models should be obtained from a high Reynolds number LES simulation on a coarse grid. Of course, such a simulation requires a good model for the near-wall region in the first place or a prohibitive number of grid points. The optimal control framework is applied here to conduct such a simulation without a priori knowledge of the necessary wall stress boundary conditions. The case of a channel flow with constant pressure gradient is considered (∂ p /∂x = 1). The objective is to optimize the shear stresses τ w 12 and τ w 32 in order to minimize a given cost function. As stated earlier, we only consider the case with zero transpiration velocity. The mathematical formulation is detailed in the following subsections.

State equation

The problem considered is governed by the unsteady, incompressible, filtered Navier-Stokes equations as well as the divergence-free constraint which arises from continuity. The governing equations read:

∂u i ∂t + ∂u i u j ∂x j = - ∂p ∂x i + δ i1 + ∂ ∂x j (ν + ν t ) ∂u i ∂x j + ∂u j ∂x i (2) ∂u j ∂x j = 0,
where the δ i1 represents the constant streamwise pressure gradient. Note that no specific notation is used to describe the spatial filtering associated with the LES formulation. Each variable in the previous and subsequent equations should be understood as a low-pass filtered version of the actual variable (e.g. u i ≡ u i , where the overbar stands for the filtering operator). Equation ( 2) is valid for any subgrid model based on the Boussinesq assumption. The boundary conditions are:

∂u ∂y = 1 ν w τ w 12 v = 0 (3) ∂w ∂y = 1 ν w τ w 32 or ∂u ∂y n = 1 ν w φ u v n = 0 (4) ∂w ∂y n = 1 ν w φ w
where the subscript n stands for the outward normal to the walls (y when y = +1 and -y when y = -1). Also ν w is the wall value of the total dynamic viscosity ν + ν t . Since the wall-normal velocity, v, is zero in the present study, the eddy viscosity ν t is zero at the walls and ν w = ν. From Eqs. ( 3) and ( 4), (φ u , φ w ) = (τ w 12 , τ w 32 ) at y = +1 and (φ u , φ w ) = -(τ w 12 , τ w 32 ) at y = -1. In the classical optimal control procedure the objective is to reduce the given cost function for some period of time. This method has been proven to be efficient. [START_REF] Abergel | On some control problems in fluid mechanics[END_REF] However, this is a very expensive procedure in terms of storage and manipulation of many 3D fields over the entire period under consideration. [START_REF] Choi | Feedback control for unsteady flow and its application to the stochastic burgers equation[END_REF] We therefore make use of a more affordable sub-optimal procedure in which the state equation is first discretized in time, then a control procedure is used to minimize the cost function over a short period of time (the time step) at each time step. [START_REF] Bewley | Optimal feedback control of turbulent channel flow[END_REF] This method does not necessarily provide the 'best' answer but it is much more cost effective than the optimal strategy. The equation of state (2) is therefore discretized in time by assuming an implicit discretization:

u n+1 i + 2β∆t ∂P ∂x i + ∂u i u j ∂x j - ∂ ∂x j (ν + ν t ) ∂u i ∂x j + ∂u j ∂x i n+1 = RHS n (5) -2β∆t ∂u n+1 j ∂x j = 0
The boundary conditions, Eq. ( 4), apply to Eq. ( 5) and the second order Crank-Nicolson method is recovered if the parameter β is chosen to be β = 1 4 . Note that in Eq. ( 5) only the terms involving the solution at the current time step n + 1 are written explicitly. This is because in the sub-optimal framework developed in Bewley et al, [START_REF] Bewley | Optimal feedback control of turbulent channel flow[END_REF] only the terms at time n + 1 in the semi-discretized state equation are used in the optimization process. The terms which depend only on the variables at the previous time step n are gathered in the generic notation RHS n and disappear in the analytical development -see Eq. (12). Note that there is a slight inconsistency between the time integration used to solve the state equation (explicit third-order Runge-Kutta, Section 3) and time integration supposed in the process of defining the control problem (implicit Crank-Nicolson). The optimal control problem and the resulting optimization algorithm proposed in Section 4.4 may therefore be considered as an approximation of the 'exact' control problem. However, in view of the results shown in Sections 5 and 6, and the small time steps used for the explicit scheme, this approximation is sufficient to provide valuable insight for identifying a practical wall model for coarse grid LES.

Cost function

The goal of the sub-optimal approach is to provide numerical shear stresses to the flow solver so that the overall solution is consistent with what is expected in a channel flow. In this particular case, the mean velocity profile should be well approximated by a logarithmic profile. Therefore the cost function is chosen to be a measure, at each time step, of the difference between the actual mean velocity profile and the target profile u + ref = 1 κ ln y + + C. This difference can be expressed as:

δ u (y) = 1 A (u -u ref ) dx dz (6)
where A stands for the channel area in the homogeneous plane. While a logarithmic mean velocity profile has been chosen as the target profile, in fact, any target profile could be used. Notably, a more realistic shape could be used in the centerline region. The logarithmic profile is expected in the 'near wall' region if the grid spacing is large enough so that the first grid point belongs to the log region. In the same way one can define an error in the w mean velocity profile as:

δ w (y) = 1 A (w -w ref ) dx dz (7)
where the reference velocity in the z direction is obviously w ref = 0 in the present case. A functional of the form :

J (φ) = +1 -1 δ u (y) 2 + δ w (y) 2 dy (8)
provides a measure of the difference between the actual and the target mean velocity profiles. The control parameter φ is defined as φ = (φ u , φ w ). A second term related to the 'cost' of the control is added to build the final cost function to be minimized:

J (φ) = +1 -1 δ u (y) 2 + δ w (y) 2 dy
Error in mean profiles

+ α A w φ 2 u + φ 2 w dx dz Cost (9)
The second term is quadratic in φ and was introduced in previous studies on drag reduction to monitor the energy input by the control system and reduce the cost of the control. [START_REF] Bewley | Optimal feedback control of turbulent channel flow[END_REF] In the present study the quadratic term is introduced to avoid the numerical instabilities which would arise if the imposed shear stresses become too large during the optimization process. The parameter α is tuned to balance the effects of the two terms in J . Note that in the particular case of the channel flow considered in the present study, the mean velocity < w > in the spanwise direction tends naturally towards zero over a sufficiently long time of integration. The term δ w is kept in the definition of J for symmetry and also to show that other parallel flows could be investigated with the present strategy. More generally, the cost function could include all velocity components.

Adjoint problem

The gradient of the cost function J with respect to the control parameter φ is estimated by using the Fréchet differential 27 defined for any functional F as:

DF Dφ φ = lim ǫ→0 F(φ + ǫ φ) -F(φ) ǫ ( 10
)
where φ is an arbitrary direction. From Eq. ( 9) the derivative of the cost function J is:

DJ Dφ φ = Ω 2δ u A U + 2δ w A W dx dy dz + 2α A w φ u φu + φ w φw dx dz (11)
where U and W denote the Fréchet derivatives of u and w respectively. An adjoint problem must be formulated to estimate the gradient of J since the derivatives U and W are unknown. The first step is to take the derivative of the semi-discrete state equation (5) to obtain:

U i + 2β∆t ∂P ∂x i + U j ∂u i ∂x j + u j ∂U i ∂x j - ∂ ∂x j (ν + ν t ) ∂U i ∂x j + ∂U j ∂x i = 0 (12)
-2β∆t ∂U j ∂x j = 0 with the boundary condition:

∂U 1 ∂y n = ∂U ∂y n = 1 ν w φu U 2n = V n = 0 (13) ∂U 3 ∂y n = ∂W ∂y n = 1 ν w φw
The right hand side term in Eq. ( 12) is now zero since the flow field at time step n does not depend on the shear stresses to be imposed for the current time step. Therefore the superscript 'n + 1' has been dropped for clarity. Note also that the Fréchet derivative of the eddy viscosity was supposed to be zero, viz. Dν t /Dφ = 0, to obtain Eq. ( 12). Moreover, this system of equations is linear in the variables U i and P, where P is the Fréchet derivative of pressure. Therefore it can be written in the form:

A Θ = 0 ( 14
)
where A is a linear operator acting on the vector Θ = (U i , P) T .

To bypass the resolution of the differential problem Eq. ( 14) with unknown boundary conditions (13), an adjoint operator A * is formulated by considering the equation

< A Θ, Ψ >=< Θ, A * Ψ > +BT, ( 15
)
where < a, b > denotes the inner product defined as the integral over the volume of the dot product of the two terms a and b and Ψ is the adjoint state vector Ψ = (η i , π) T . Practically, the adjoint operator is formed by using successive integrations by parts to turn all the spatial derivatives acting on Θ to derivatives acting on Ψ. Some boundary terms arise during this process, which are contain in the term BT of Eq. ( 15). It is straightforward to show that the adjoint operator A * acting on the adjoint state vector Ψ is such that

A * Ψ =    η i + 2β∆t ∂π ∂x i + η j ∂u j ∂x i -u j ∂η i ∂x j -∂ ∂x j (ν + ν t ) ∂η i ∂x j + ∂η j ∂x i -2β∆t ∂η j ∂x j (16)
and that the boundary terms are:

BT = 2β∆t w Pη 2n -ν w η i ∂U i ∂y n -U i ∂η i ∂y n + η i ∂U 2n ∂x i -U i ∂η 2n ∂x i dx dz (17)
Note that several terms arising from the integration by parts have not been included in BT because of the zero transpiration assumption, viz. u 2 = U 2 = 0. From Eq. ( 14), the relation ( 15) defining the adjoint operator reduces to

< A * Ψ, Θ >= -BT (18)

Gradient estimate

We now have the liberty to choose boundary condition and right hand side terms for the adjoint problem such that the relation (18) can be utilized to compute the gradient of J . By comparing Eqs. ( 18) and ( 11) it appears that a judicious choice for the definition of the adjoint problem is:

A * Ψ =      2δu A 0 2δw A 0      (19)
with boundary conditions at the wall:

∂η 1 ∂y n = ∂η 3 ∂y n = η 2n = 0 (20)
In doing so, Eq. ( 18) can be re-written as

DJ Dφ φ = 2β∆t w η 1 φu + η 3 φw dx dz + 2α A w φ u φu + φ w φw dx dz (21)
Since Eq. ( 21) is valid for all directions φ = ( φu , φw ), the gradient of J may be extracted:

27        DJ Dφ = 2β∆t η 1,w + 2α A φ u if φ = (φ u , 0) DJ Dφ = 2β∆t η 3,w + 2α A φ w if φ = (0, φ w ) (22)
where η i,w stands for the value of η i,w at the wall. A control procedure using a simple gradient algorithm at each time step may now be proposed such that:

φ n+1 k+1 = φ n+1 k -µ DJ (φ n+1 k ) Dφ (23)
Note that the adjoint operator A * depends on the state vector (u i , P ) T at time n+1 so that the state equation and the adjoint problem must be solved simultaneously to obtain the sub-optimal shear stresses. The following algorithm is used a each time step of the flow solver to obtain the optimized boundary conditions:

1. Start with the state vector (u i , P ) T at time n, the adjoint vector (η i , π) T and control parameter φ at sub-iteration n k , 2. Use Eq. ( 19) with boundary condition (20) to compute the adjoint vector at sub-iteration n k+1 . For this purpose, the operator for the adjoint velocity, Eq. ( 16), is re-written as:

1 2 η n k+1 i + 1 2 η i + 2β∆t ∂π ∂x i + η j ∂u j ∂x i -u j ∂η i ∂x j - ∂ ∂x j (ν + ν t ) ∂η i ∂x j + ∂η j ∂x i (24)
where only the first term is taken at sub-iteration n k+1 , the others being considered at subiteration n k . A Poisson equation is solved for the adjoint pressure at each sub-iteration to enforce the divergence-free condition on η n k+1 i .

3. Compute the gradient of the cost function at sub-iteration n k+1 by using Eq. ( 22), 4. Compute the control parameter φ at sub-iteration n k+1 by using Eq. ( 23), 5. Compute the L 2 -norm of the difference φ n k+1 -φ n k .

6. If it is more than 0.1% of the L 2 -norm of φ n k , then go back to step 1. Otherwise, use φ u and φ w as approximate boundary conditions to compute the state vector at time n + 1.

5 Results from optimal simulations

Statistics

Several LES's have been performed to test the optimal strategy described in the previous section. In this section, the grid is uniform in all directions with 32x33x32 cells and the domain size is (2πh, 2h, 2πh/3) where h is the channel half height. The Reynolds number based on the friction velocity, u τ , and h is Re τ = 4000. In wall units, the grid spacing is ∆x + ≈ 785, ∆y + ≈ 242 and ∆z + ≈ 262. Since a staggered mesh system is used, the first u velocity point is located at y + ≈ 121. At this coarse resolution, most of near-wall turbulent structures cannot be captured so that an effective wall model is necessary to compensate. The computation was run for a sufficiently long : sub-optimal computation of τ w 12 only; : < u > + = 2.41 ln y + + 5.2 period to be statistically independent of the initial condition and then statistics were accumulated over a time period of order 20 h/u τ . Several numerical tests were performed to fix the coefficients µ and α. It was found that the value µ = 5 × 10 3 ensured good convergence of the steepest descent algorithm, while α = 4 × 10 -5 ensured that the whole algorithm is stable.

Figure 3 shows the mean velocity profile from the LES in which the optimal procedure of Section 4 was used to obtain the approximate boundary condition. The mean value of τ w 12 was either provided by the optimal procedure itself or re-computed so that the mean velocity at the first point coincides with the logarithmic law < u + >= 2.41 log y + + 5.2. The results are very similar in the two cases. The overall agreement is much better than with the shifted model, [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] Eq. ( 1). The profile from Figure 2b is shown in Figure 3 for easier comparison. An artificial boundary layer still develops between the second and the third grid point but its amplitude is much smaller than with the shifted model. The deficit in the log-law intercept is of order 0.25 compared to 1.35 in Figure 2b.

Note that the optimized wall stress boundary conditions produce a mean velocity profile that is nearly exact for the first two grid points. The small error in the channel interior is believed to be due to the sub-optimal formulation in which the wall stresses are optimized only over each time step. However, the result in Figure 3 shows a clear improvement in comparison with the shifted model. An additional computation was performed where only τ w 12 was optimized (u 3 was set to zero at the wall). Figure 3 shows that in this case the mean velocity profile is intermediate between the two previous computations. The conclusion is that both τ w 12 and τ w 32 must be optimized (or modeled carefully) if a velocity profile close to the target one is sought. In the following, only the computation where both shear stress boundary conditions are optimized is considered.

Shear stress structure

From Eq. ( 1), the wall shear stresses in the LES with the shifted model should be perfectly correlated with the velocity components at the first plane, shifted in the upstream direction by the amount ∆ s . This is shown in Figure 4 which displays typical iso-lines of τ w 12 , u, τ w 32 and w. The upstream shift is only a fraction of the cell size, viz. ∆ s ≈ 0.67∆x, and is hardly visible in the figure. Both τ w 12 and u are characterized by elongated structures in the streamwise direction while τ w 32 and w are characterized by more rounded structures in the x -z plane. The reference data from the sub-optimal computation are plotted in Figure 5. The perfect correlation between the wall stresses and the velocity components at the first off-wall plane no longer holds. The streamwise velocity u still displays elongated structures in x as well as τ w 32 . In contrast, τ w 12 and w are less elongated. Eight fields from the sub-optimal LES have been used to compute correlation coefficients between the velocity gradients, the velocity components at the first off-wall plane and the wall stresses. The best results (those with the correlation greater than 0.3) are reported in Figure 6. The position of the perfect correlations between τ w 12 and u, τ w 32 and w which are assumed in the shifted model, [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] Eq. ( 1), is also reported. Note that in the sub-optimal calculation, τ w 12 is almost not correlated with the streamwise velocity downstream shifted by ∆ s . The correlation is better when an upstream shifted version of u is used instead. On the other hand the maximum of correlation between τ w 32 and w is located close to the assumed downstream shift ∆ s . The best correlation for τ w 12 is found with the spanwise derivative of w, whereas τ w 32 best correlates with du/dz. Other non negligible correlation coefficients are found between τ w 12 and du/dx, and τ w 32 and dv/dz. No clear picture is available yet to explain these correlations. Finally, note the fairly good and somewhat expected negative correlation between τ w 12 and v. This correlation supports the underlying idea in the ejection model of Piomelli et al. [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] A more systematic way of exploiting correlations between the approximate boundary condition and the velocity field close to the wall is presented in Section 6

Discussion

The reduced deficit in the intercept with the sub-optimal computation is associated with a better representation of the gradient of the mean velocity profile within the first few grid points. This is better seen in Figure 7 which displays the dimensionless shear, κy + d < u > + /dy + , as a function of the distance to the wall. Theoretically, this quantity should be equal to unity in the case where the sub-optimal strategy is used since the value of κ is taken to be κ = 1/2.41, consistently with the target velocity profile < u > + = 2.41 log y + + 5.2. For the case where the shifted model, [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] Eq. ( 1), is used, this quantity is not expected to be exactly unity since the value picked for κ is somewhat arbitrary. However this quantity is constant if the mean profile follows a logarithmic behavior. In practice, κy + d < u > + /dy + is not found to be equal to unity when computed from the exact log profile < u > + = 2.41 log y + + 5.2 because of the large errors associated with the second-order finite differences on the coarse grid considered. This is shown in Figure 7 which also demonstrates that the non-dimensionalized mean gradient from the sub-optimal computation follows closely (to within 8%) its value from the exact log law. The mean gradient from the computation with the shifted model 9 is 25% below its expected value, which explains the deficit in the intercept observed in Figure 3.

The mean velocity gradient appears in the mean momentum equation which reads: This equation reduces to:

< u ′ v ′ >= y + ν ∂ < u > ∂y + ν t ∂u ∂y (25)
∂ < u > ∂y ≈ < u ′ v ′ > -y ν+ < ν t > (26)
under the usually well verified assumption ν t ∂u ∂y ≈< ν t > ∂u ∂y . Eq. ( 26) reveals that for a given amount of eddy viscosity, the mean velocity gradient is directly related to the difference between the total (-y) and resolved (-< u ′ v ′ >) shear stress, viz. < u ′ v ′ > -y, the larger the difference the larger the gradient. This is confirmed in Figure 8 which shows that the artificial condition provided by the sub-optimal strategy leads to an equilibrium in which the resolved part of the stress is smaller. It is also worth noting that the quantity < u ′ v ′ > -y is quite different in the two computations but that the absolute value of the resolved stress are both close to unity (| < u ′ v ′ > | max ≈ 0.9 for the computation with the shifted model of Piomelli et al., [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] | < u ′ v ′ > | max ≈ 0.86 for the sub-optimal case). Since there is a large difference in the velocity gradient (it is 30% higher in the sub-optimal case), the production of kinetic energy -< u ′ v ′ > d < u > /dy is necessarily greater in the sub-optimal case, as shown in Figure 9a. Figure 9b shows that this is also true for the dissipation rate of kinetic energy, a direct consequence of the fact that production and dissipation balance reasonably well in this high-Reynolds number channel flow. Note that the total dissipation is very well approximated by the subgrid scale (SGS) dissipation ǫ sgs =< 2ν t S ij S ij > that has been reported in Figure 9b. The SGS dissipation can further be approximated (to within a few percents) by ǫ approx = 2 < ν t >< S ij S ij >. The form of the approximate dissipation ǫ approx indicates that the sub-optimized boundary condition can act on the flow field and generate an equilibrium with higher dissipation in the near wall region by either increasing the mean eddy-viscosity or increasing the velocity fluctuations. The mean eddy-viscosity profiles in the two computations are found to be : sub-optimal computation; • : shows the position where the correlation between τ w 12 and u, and τ w 32 and w in the shifted model of Piomelli et al. [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] is assumed to be unity. : not exact differentiation applied to < u > + = 2.41 ln y + + 5.2. : sub-optimal computation; + : Shifted model of Piomelli et al.; [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] : Resolved LES computation. [START_REF] Kravchenko | Zonal embedded grids for numerical simulations of wall-bounded turbulent flows[END_REF] almost identical (not shown). The turbulent fluctuations are higher in the sub-optimal case than in the shifted model 9 computation -see Figure 10. The agreement with the fully-resolved LES data of Kravchenko et al. [START_REF] Kravchenko | Zonal embedded grids for numerical simulations of wall-bounded turbulent flows[END_REF] turns out to be better with the sub-optimized boundary condition regarding the normal and spanwise direction, worse regarding the streamwise direction. In any case, the reasoning given above indicates that the increase in the turbulent fluctuations when the sub-optimal strategy is used -Figure 10-is consistent with the previous findings on the mean velocity gradient -Figure 7-, resolved shear stress -Figure 8-, turbulence production and dissipation -Figure 9. Similar results have been obtained for two other Reynolds numbers, namely Re τ = 640 and 20000 (not shown). For the first case, the mean velocity profile from the sub-optimal computation is slightly closer to the logarithmic law than that obtained from the shifted model. [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] For the case Re τ = 20000, the deficit in the intercept is the same than for the case Re τ = 4000, viz. 0.25. Also the profiles of turbulent fluctuations do not change if plotted in outer coordinates y/h. This means that the asymptotic regime in terms of Reynolds number has been reached: results of the same quality as above can be obtained for arbitrarily large Reynolds numbers with the same coarse grid.

A practical wall model via Linear Stochastic Estimation

While the sub-optimal control strategy for generating wall stresses could be used as a wall model for coarse-grid LES, its cost is approximately 20 times greater than of an LES on the same grid compared to an explicit wall stress model such as (1). Furthermore, a target mean velocity profile must be provided to define the objective function. It may be possible to lower the cost of control strategy, but we have not investigated that possibility here. The real strength of the optimal control strategy is that it yields wall stress boundary conditions that work for coarse-grid, high Reynolds number LES. Thus, a reference data set can be generated against which new wall models can be compared. The reference data can even be used to derive new wall models. Such an approach is described in this section.

The most desirable wall stress model would be similar to (1) in that it would be an explicit relation between the wall stresses and the velocity field. One could, for instance, require the wall stress model to be the best possible mean square estimate of the sub-optimal wall stress as a function of the local velocity field. This is the conditional average of the wall stress given the local velocity field (a conditional average is necessary because the wall stress may have a stochastic, or unpredictable, component with respect to the local velocities). It is denoted by τ w i2 (x, z)|E , where E is a vector of events. In the present study, E will be a vector containing the local velocity field, but it could easily contain pressure, velocity gradients, quadratic products, or any other quantities which might characterize the wall stresses. The conditional average embodies so much statistical information that it is unlikely that it could be found exactly, but it can be approximated by its linear stochastic estimate (LSE), given by (see [START_REF] Adrian | Appoximation of turbulent conditional averages by stochastic estimation[END_REF] for instance):

τ w i2 (x, z)|E ≈ τ w i2 (x, z) = L ij E j i = 1, 3, j = 1, 2, 3, . . . , N, (27)
where N is the number of events being considered, and L ij are estimation coefficients relating τ w i2 to E j . By the statistical orthogonality principle, [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] the mean square error between τ w i2 and τ w i2 is minimized when the event data are uncorrelated with the error e i = τ w i2 -τ w i2 :

e i E k = (τ w i2 -τ w i2 )E k = 0. ( 28
)
Substituting (27), the estimation coefficients L ij are governed by: The matrix E j E k is invertible provided the events are linearly independent. Thus, to obtain the LSE, the correlations τ w i2 E k and E j E k must be found and the events must be selected that best characterize the wall stress.

τ w i2 E k = L ij E j E k . (29)
Though the technique employed here is essentially the same as that of Bagwell, et al., [START_REF] Bagwell | Improved approximation of wall shear stress boundary conditions for large eddy simulation[END_REF] the results are different since the reference data used here is already known to work well for a coarse grid LES, whereas Bagwell's reference data comes from a direct numerical simulation at low Reynolds number. As will be shown below, an event field consisting of the nearby velocities is sufficient to yield wall models of the form (27) that produce wall stresses that have high correlation with the optimal wall stresses. Moreover, the new wall models, when used in an LES, will be shown to reproduce the results of the sub-optimal control LES.

LSE predictions

To implement the LSE, eight velocity fields, well separated in time, and their sub-optimal wall stresses from the Re τ = 4000 simulation discussed in section 4.1 were used. In implementing the LSE, it was found that better results were obtained by estimating only the fluctuating part of the wall stresses from the fluctuating part of the velocity field. To this end, the instantaneous plane averaged wall stress was subtracted from the wall stresses and the average over the first off-wall plane of the horizontal velocity was subtracted from its velocity component for each sample. The LSE can be written as a convolution sum in the wall-parallel, homogeneous directions. The exact

L 11 ijk [u m-i,j,n-k -u(•, y 1 , •) ] + nz+1 k=-nz ny j=1 nx-1 i=-nx L 13 ijk w m-i,j,n-k , (30)
τ w ′ 32 | m,n = nz-1 k=-nz ny j=1 nx+1 i=-nx L 31 ijk [u m-i,j,n-k -u(•, y 1 , •) ] + nz k=-nz ny j=1 nx i=-nx L 33 ijk w m-i,j,n-k , (31)
where the parameters n x , n y , n z determine the number of velocity points used in the convolution sums to estimate the wall stress at each wall location (denoted by the subscript pair m, n). Note that the summation relating τ w ′ 12 and the spanwise velocity w has different indices than the summation for u in (30). This is due to the staggered grid. Also note that the wall-normal velocity does not appear in the LSE wall model. This is because the wall-normal velocity is linearly dependent on the wall-parallel velocities and thus cannot be used as an independent event in the LSE. The coefficients L 11 , L 13 , L 31 , and L 33 , are determined by requiring that the error be orthogonal to the events (velocities) as in equation (28). For an a priori comparison of the wall stress fluctuations predicted by the LSE to those of the optimal strategy one possibility is to compute the correlation coefficient:

C i2 = τ w ′ i2 τ w ′ i2 ( τ w ′ i2 ) 2 1/2 (τ w ′ i2 ) 2 1/2 , (32)
where i = 1, 3 and the average denoted by, • , is taken over all of the samples. Another frequently used quantity for comparison is the relative mean square error:

R i2 = (τ w ′ i2 -τ w ′ i2 ) 2 (τ w ′ i2 ) 2 , (33)
where, again, i = 1, 3.

Table 1 shows the correlation coefficients and relative mean square error for several choices of the parameters n x , n y , and n z . The first noteworthy observation is that the correlations between the sub-optimal wall stress fluctuations and those predicted by the LSE (30), (31) increase significantly when velocity data is included from the first two wall-parallel planes (compare n y = 1 and n y = 2). The estimates with n x = 15 and n z = 15 use nearly all of the velocity data in each wall-parallel plane, but note that this results in little improvement over the n x = 4, n z = 4 results. Figures 11-14 show contour plots of the LSE coefficients for the case n x = 2, n y = 2, n z = 2. Each figure contains two plots, the top plot corresponding to the first plane of velocity data and the bottom plot to the second. The plots indicate several symmetries in the LSE coefficients. For example L 11 is symmetric with respect to reflection in the spanwise direction whereas L 13 is anti-symmetric. These symmetries could be exploited to reduce the number of free coefficients in a wall model, but we have not explored that here.

Bagwell, et al. [START_REF] Bagwell | Improved approximation of wall shear stress boundary conditions for large eddy simulation[END_REF] used entire planes of velocity data to form the LSE, but as results in Table 1 suggest, a small number of data points may be sufficient. These a priori results suggest that the LSE can reproduce the sub-optimal wall stress fluctuations reasonably well, but even correlations in excess of 80% do not guarantee that the LSE can perform well as a wall model. This can only be checked via an actual simulation. n z = 2. Simulations with LSE models based on smaller stencils of velocity data did not work as well, whereas LSE models based on larger stencils reproduced the sub-optimal results. The turbulent fluctuations are also in perfect agreement with those of the sub-optimal wall stress simulation; see Figure 16. These results are remarkable in that the wall stress model (34) reproduces the results of the sub-optimal simulation at a cost only slightly higher than that of a simulation with no wall model. However, the model given by (34) will not be of much use if new coefficients need to be derived for different Reynolds numbers or for different grid resolutions. Fortunately, we have found that simulations with (34) perform well over a large range of Reynolds numbers using the same coefficients derived from the Re τ = 4000 sub-optimal data. Figure 17 shows the results of using (34) at Reynolds numbers: Re τ = 640, 4000, and 20000. The simulations were conducted on the same uniform grid as the Re τ = 4000 simulation with the log law (36) used to determine the mean wall stress. In the Re τ = 640 LES, the log law was used to relate the plane-averaged streamwise velocity from the second plane to the mean wall stress because the first plane is in the buffer region at y + ≈ 20.

Similarly, the same coefficients were used with (34) in a simulation with a refined grid at Re τ = 20000. The number of cells was doubled in each direction to 64 × 65 × 64. The resulting mean flow profile is shown in Figure 18. The log region intercept is still slightly under-predicted, but the mean flow now exhibits a wake-like behavior near the channel center as has been observed in high Reynolds number channel simulations in which the near-wall region is resolved. : < u > + = 2.41 ln y + + 5.2

Conclusion

A sub-optimal control strategy has been successfully applied to a coarse grid LES of a channel flow in order to generate boundary conditions for the wall shear stresses to be used in high Reynolds number LES. The mean velocity profile from the sub-optimal calculations is in very good agreement with the log law, even if the resolution near the wall is very coarse. The method is robust enough to reach the asymptotic range of Reynolds numbers where the same quality of result is obtained for arbitrarily large ret with the same coarse grid. A set of fields from the optimally controlled simulation has been used to compute the convolution coefficients of a linear stochastic model relating the boundary condition to the velocity field in the first few off-wall planes. The results from the explicit model obtained in this fashion are close to those from the sub-optimal computation. The stochastic model can be used in practical LES calculations.

Figure 1 :

 1 Figure 1: Staggered grid system used in the study.

Figure 2 :

 2 Figure 2: Mean velocity profile from LES with coarse grid and existing wall stress models. (a): Re τ = 640 with shifted model (1); (b): Re τ = 4000 + with shifted model, and with the two-layer model. The grid is uniform with 32x33x32 points.: < u > + = 2.41 ln y + + 5.2

Figure 3 :

 3 Figure 3: Mean velocity profiles for Re τ = 4000.: sub-optimal computation;

Figure 4 :

 4 Figure 4: From top to bottom: Contours of τ w 12 , u, τ w 32 and w from the LES with the shifted model of Piomelli et al. 9 at Re τ = 4000. The velocity components are from the first off-wall plane.

Figure 5 :

 5 Figure 5: From top to bottom: Contours of τ w12 , u, τ w 32 and w from the LES with the sub-optimal strategy at Re τ = 4000. The velocity components are from the first off-wall plane.

Figure 6 :

 6 Figure 6: Best correlation coefficients for (a) τ w 12 , (b) τ w 32 .: sub-optimal computation; • : shows the position where the correlation between τ w 12 and u, and τ w 32 and w in the shifted model of Piomelli et al. 9 is assumed to be unity.

Figure 7 :

 7 Figure 7: Non-dimensionalized mean velocity gradient κy + d < u > + /dy + for Re τ = 4000. : sub-optimal computation; + : Shifted model of Piomelli et al.; 9: not exact differentiation applied to < u > + = 2.41 ln y + + 5.2.

Figure 9 :Figure 10 :

 910 Figure 9: Turbulence (a) production (b) subgrid scale dissipation from the sub-optimal computation scaled by the same quantity from the LES with the shifted model of Piomelli et al. 9 for Re τ = 4000.

Figure 11 :

 11 Figure 11: Contours of LSE coefficient L 11 for n x = 2, n y = 2, and n z = 2, in the xz-plane. The top plot is for the first plane (j = 1) and the bottom plot is for the second plane.

Figure 12 :

 12 Figure 12: Contours of LSE coefficient L 13 . Refer to figure 11 for details.

Figure 13 :Figure 14 :

 1314 Figure 13: Contours of LSE coefficient L 31 . Refer to figure 11 for details.

Figure 17 :

 17 Figure 17: Mean velocity profiles from LES's with the LSE model on uniform 32x33x32 grid: : Re τ = 640; : Re τ = 4000; : Re τ = 20000; : < u > + = 2.41 ln y + + 5.2

28

 28

Figure 18 :

 18 Figure 18: Mean velocity profiles from LES's for Re τ = 20000 and different uniform grids. : 64x65x64 grid points and LSE model; : 32x33x32 grid points and LSE model; + : 64x65x64 grid points and shifted model of Piomelli et al.; 9 : < u > + = 2.41 ln y + + 5.2

Table 1 :

 1 Correlation coefficents and relative mean square errors for the LSE wall model fluctuations versus the sub-optimal wall stress fluctuations for different numbers of velocity data use in the estimates. The fifth line, set apart by horizontal lines, is for the coefficients used as an actual wall model in section 6.2 n x n y n z C 12 C 32 R 12 R 32

	1	1	1	.45 .59 .78 .65
	4	1	4	.48 .64 .77 .59
	15 1 15 .51 .66 .74 .57
	1	2	1	.73 .82 .46 .32
	2	2	2	.76 .89 .42 .21
	4	2	4	.78 .91 .40 .18
	15 2 15 .79 .91 .37 .18
	2	3	2	.79 .90 .42 .21
	4	3	4	.81 .92 .34 .15

Acknowledgments This research was partly supported by AFOSR Grant F49620-7-1-0210.

: sub-optimal computation; + : Shifted model of Piomelli et al. : sub-optimal computation; + : Shifted model of Piomelli et al.; [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF] : LSE model; : < u > + = 2.41 ln y + + 5.2

LSE wall model results

The LSE was used to model the fluctuating portion of the wall stresses in terms of the velocity field and must be combined with a model for the mean wall stress to be used as an LES wall model. The actual model used is:

where τ w ′ 12 and τ w ′ 32 are given by ( 30) and (31), respectively. The mean wall stress τ w 12 is found by assuming that the plane-averaged streamwise velocity at the first off-wall grid location and the mean wall stress are related by the logarithmic law of the wall:

By using the LSE wall model (34) in the same flow discussed in section 4.1, with Re τ = 4000 and the same resolution, it was found that the wall models based on the one plane LSE, n y = 1 resulted in mean velocity profiles (not shown here) that were not as good as the profiles from the sub-optimal simulation. However, the wall model based on the two plane LSE was able to reproduce nearly exactly the results of the sub-optimal simulation. Shown in Figure 15 is the mean velocity profile from a simulation with an LSE wall model of the form (34) with n x = 2, n y = 2, and : sub-optimal computation; : LSE model.

List of Tables