Wall-modeled large-eddy simulation of three-dimensional turbulent boundary layer in a bent square duct

We conduct wall-modeled LES (WMLES) of a pressure-driven three-dimensional turbulent boundary layer (3DTBL) developing on the floor of a bent square duct to investigate the predictive capability of three widely used wall models, namely, a simple equilibrium stress model, an integral nonequilibrium model, and a PDE nonequilibrium model. The numerical results are compared with the experiment of Schwarz and Bradshaw (J. Fluid Mech. (1994), vol. 272, pp. 183-210). While the wall-stress magnitudes predicted by the three wall models are comparable, the PDE nonequilibrium wall model produces a substantially more accurate prediction of the wall-stress direction, followed by the integral nonequilibrium wall model. The wall-stress direction from the wall models is shown to have separable contributions from the equilibrium stress part and the integrated nonequilibrium effects, where how the latter is modeled differs among the wall models. The triangular plot of the wall-model solution reveals different capabilities of the wall models in representing variation of flow direction along the wall-normal direction. On the contrary, the outer LES solution is unaffected by the type of wall model used, resulting in nearly identical predictions of the mean and turbulent statistics in the outer region for all the wall models. This is explained by the vorticity dynamics and the inviscid skewing mechanism of generating the mean three-dimensionality. Finally, the LES solution in the outer layer is used to study the anisotropy of turbulence. In contrast to the canonical 2D wall turbulence, the Reynolds stress anisotropy exhibit strong non-monotonic behavior with increasing wall distance.


1. Introduction
The capability to predict high Reynolds number turbulent flows is essential for many natural and engineering flows such as external aerodynamics of wind turbines and aircraft wings, flow over the hull of marine vehicles, atmospheric boundary-layer flow over complex landscapes and cityscapes, to name a few. However, due to extreme disparity of scales present in high Reynolds number wall-bounded turbulent flows, any attempt to simulate these flows directly on a computational grid without resorting to modeling of some sort results in prohibitively large computational cost. To resolve all the scales of the near-wall turbulent motions using direct numerical simulation (DNS), the required number of grid points scales as ( 37/14 ), where is the characteristic Reynolds number. Wall-resolved large-eddy simulation (WRLES), which resolves only the large (stress-carrying) eddies, reduces the grid point requirement to ( 13/7 ) (Choi & Moin 2012). However, this level of computational cost is still unaffordable when the Reynolds number is high. As a cost-effective alternative to the above two approaches, wall-modeled large-eddy simulation (WMLES) resolves only the energetic eddies in the outer portion of the boundary layer, while the momentum transport within the unresolved near-wall region is accounted for by augmenting the wall-flux through a wall model. Thus, the no-slip condition at the wall is replaced by the Neumann boundary condition supplied by the wall model in the form of the wall shear stress. Note that the wall shear stress is computed by the wall model without ever resolving the near-wall scales. Therefore, the grid point requirement for WMLES reduces to ( ) (Choi & Moin 2012), making the simulations of high Reynolds number flows feasible.
To date, several wall models have been proposed, most of which are based on some form of the law-of-the-wall or solving a set of simplified or full Reynolds-averaged Navier-Stokes (RANS) equations. Deardorff (1970) and Schumann (1975) were the first to recognize the need for wall modeling to perform LES of high Reynolds number plane channels and annuli to overcome lack of computing resources in the 1970s. Grötzbach (1987) later improved the model by removing the necessity of a priori knowledge of the mean wall shear stress. The geometry of the near-wall eddies was incorporated in the work of Piomelli et al. (1989) to account for the inclination of the vortical structures in the streamwise/wall-normal plane. Wang & Moin (2002) proposed an ordinary differential equations (ODE) based wall model derived from the equilibrium assumption (Degraaff & Eaton 2000), which later on was extended to compressible flows (Kawai & Larsson 2012;Bodart & Larsson 2011). The ODE equilibrium wall model excludes the nonequilibrium effects such as pressure gradient, and considers the wall-normal diffusion only. Nonequilibrium wall models based on full 3D RANS equations were investigated by Balaras et al. (1996), Wang & Moin (2002), Cabot & Moin (2000), Kawai & Larsson (2013), and Park & Moin (2014, 2016a. Yang et al. (2015) introduced the integral nonequilibrium wall model based on the integrated boundary layer equations and assumed velocity profiles, which can be considered as a compromise between the aforementioned two classes of wall models. Several efforts have also been directed toward formulating wall models which are not based on RANS. Bose & Moin (2014) and Bae et al. (2019) proposed a differential filter-based wall model which introduced a slipvelocity applied in the form of Robin boundary condition at the wall. Chung & Pullin (2009) proposed a virtual wall model with a slip velocity boundary condition specified on the lifted virtual wall. Gao et al. (2019) extended this virtual wall model in a generalized curvilinear coordinate. Advances on WMLES were reviewed by Piomelli & Balaras (2002), and more recently by Larsson et al. (2016) and Bose & Park (2018). With the development of novel wall models and increase in the computing capacity, WMLES is becoming an indispensable tool for predictive but affordable scale-resolving simulation of practical engineering flows at high Reynolds numbers. Recent applications to external aerodynamics applications include 3 simulation of a wing-body junction flow (Lozano-Durán et al. 2021) and flow over a realistic aircraft model in landing configuration deploying high-lift devices (Goc et al. 2021).
Although WMLES is now gaining popularity as a high-fidelity tool balancing the computational cost and the accuracy, with the potential to be used for design and optimization in practical engineering applications because of its reasonable turnaround times, comprehensive benchmark studies on the comparison of different wall models are lacking. Park (2017) compared the performance of ODE equilibrium wall model and PDE nonequilibrium wall model in a separating and reattaching flow over the NASA wall-mounted hump. Lozano-Durán et al. (2020) tested three RANS based wall models (ODE equilibrium wall model, intergral nonequilibrium wall model and PDE nonequilibrium wall model) in a threedimensional transient channel flow. For the latter study, it is worth noting that the three wall models were not tested using the same LES code. The lack of a like-for-like comparison of different wall models, especially in flows with nonequilibrium effects such as mean flow three-dimensionality and pressure gradient, warrants a systematic study of various wall models under the identical settings of the same solver and the same flow conditions. This will facilitate a clear assessment of the differences in the performance of different wall models, both in terms of the accuracy and the computational cost. Foregoing in view, in the present work, we test three wall models in a three-dimensional turbulent boundary layer (3DTBL) flow: an ODE equilibrium wall model (EQWM), an integral nonequilibrium wall model (integral NEQWM), and a PDE nonequilibrium wall model (PDE NEQWM). These three models respectively represent increasing model complexity with correspondingly increasing physical fidelity for predicting 3DTBLs. The equilibrium wall model assumes that the velocity profile is unidirectional and neglects all nonequilibrium effects, while the latter two are capable of representing skewed velocity profiles and incorporate some or all nonequilibrium effects, albeit in an averaged sense.
Before describing the 3DTBL in more detail, a few remarks are in order regarding the suitability of the current choice of 3DTBL flow to conduct the comparative study of wall models with different physical fidelity. Historically, much of the research on wall turbulence has focused on statistically two-dimensional (2D) equilibrium turbulence in simple geometries (e.g., channel, pipe and flat plate). Different wall models perform equally well therein, making it hard to justify the use of more complex models. Furthermore, many practical flows of interest, such as those found on the swept wings of aircraft, wing/body juncture, bow/stern regions of ships and turbomachinery, are strongly affected by the mean flow three-dimensionality. Such 3DTBLs challenge the validity of the theories and models established from the canonical 2D wall turbulence and thus provide a good stage to exhibit the distinctive capabilities of different wall models. Therefore, the current study of turbulent boundary layer with mean-flow three-dimensionality is well suited to test different wall models, and to explain the physical origins of the differences in the results of these models.
The 3DTBLs can be classified as pressure-driven (also termed skew-induced (Bradshaw 1987) or inviscid-induced (Lozano-Durán et al. 2020)), or shear-driven (also termed viscousinduced (Lozano-Durán et al. 2020)) ones, according to the mechanisms by which the mean three dimensionality is introduced into the flow. For the pressure-driven 3DTBLs, the crossflow is induced by the imposition of spanwise pressure gradient. More specifically, the mean three-dimensioanlity is produced by reorienting (tilting) the existing mean spanwise vorticity to generate non-zero streamwise vorticity. This process is often referred to as "inviscid skewing" due to its quasi-inviscid nature, and streamwise variation in the imposed spanwise pressure gradient often facilitates this vorticity tilting (Coleman et al. 2000). Examples of this type of 3DTBLs include flows in a square duct with a bend (Schwarz & Bradshaw 1994;Flack & Johnston 1994), in an S-shaped duct (Bruns et al. 1999), over wingbody junctures (Rumsey 2018), over swept wings (Bradshaw & Pontikos 1985), and over prolate spheroids (Chesnakas & Simpson 1994). For the shear-driven 3DTBLs, the crossflow is induced by the viscous diffusion of mean spanwise shear from the wall. Examples of this class include flows within a spinning cylinder (Bissonnette & Mellor 1974;Lohmann 1976;Driver & Hebbar 1988), over a rotating disk (Littell & Eaton 1993), over turbomachinery and in Ekman layers. In the present work, we are interested in the skew-induced cases which are more prevalent in external hydrodynamics or aerodynamics applications.
Over the past decades, studies on 3DTBL have unraveled its distinctive features which set it apart from the canonical 2D wall turbulence. First, the mean-flow direction in 3DTBL varies along the wall-normal direction, resulting in a skewed velocity profile. The law-of-thewall, which is the characteristic of the canonical 2D wall turbulence, is therefore challenged in 3DTBL. Second, the Reynolds shear stress vector is not aligned with the mean velocity gradient vector in 3DTBLs. Thus the Reynolds shear stress response in 3DTBL can lag behind or lead that predicted by the isotropic eddy viscosity models which assume perfect alignment of the two. Third, a reduction in the structure parameter (the ratio of the total Reynolds shear stress magnitude to twice the turbulent kinetic energy) is often observed in 3DTBLs, whereas this parameter is nearly constant (roughly 0.15) in the outer layer of 2DTBL. The aforementioned features of 3DTBL pose a fundamental challenge to the validity of the underlying assumptions in many turbulence models (including wall models) that are based on 2DTBL, and therefore bring into question the reliability of these models when applied to practical flows.
The numerical studies of 3DTBLs using direct numerical simulation (DNS) and largeeddy simulation (LES) have mostly focused on deformed 2D wall turbulence. These studies include channel flow subject to sudden crossflow pressure gradients (Lozano-Durán et al. 2020;Sendstad 1992), channel flow with spanwise wall motions, channel flow subject to mean strains (Coleman et al. 2000), TBL over an idealized infinite swept wing generated by a transpiration profile (Coleman et al. 2019), TBL subject to streamwise-varying pressure gradient (Bentaleb & Leschziner 2013), and TBL on a flat plate with a time-dependent freestream velocity vector (Spalart 1989). These numerical studies are limited to relatively low Reynolds number and idealized 3DTBLs due to the large computational cost. The present study focuses on a realistic, spatially-developing, pressure-driven 3DTBL over the floor of a duct with a bend (Schwarz & Bradshaw 1994), which is at a considerably higher Reynolds number than the past studies but still providing a good balance between the physical realism, the tractability of the underlying 3DTBL mechanisms, and the computational cost of the simulations.
The major objective of this work is to evaluate the three aforementioned wall models in the 3DTBL. Another objective is to understand the characteristics of the skew-induced 3DTBL, especially compared with the viscous-induced 3DTBL. The paper is organised as follows. The computational details including the flow configuration, boundary conditions and wallmodel formulations are discussed in section 2. The flow statistics obtained from WMLES are presented in section 3. Based on these results, the performance of different wall models are compared and the characteristics of this pressure-driven 3DTBL are discussed based on the anisotropy invariant map and the Johnston triangular plot. The effects of different nonequilibrium terms in wall models in terms of predicting near-wall flow direction are also quantified. Finally, conclusions are given in section 4.
Focus on Fluids articles must not exceed this page length  Bradshaw (1994)). The measurement locations in the experiment are marked as numbers 0-21 along the duct centerline. Two coordinate systems are employed. ( , , ) is a fixed coordinate system with the origin located at the inlet. ( , , ) is a curvilinear coordinate system aligned with the local duct centerline (measurements in mm).

Flow configuration
The reference configuration for the present study is the experimental setup of Schwarz & Bradshaw (1994). While numerous experimental studies have been reported on 3DTBLs (as discussed in Johnston & Flack (1996)), our choice of the reference experiment was motivated primarily by the following aspects of Schwarz & Bradshaw (1994) which we found to be favorable to the goals of this study: 1) the highest Reynolds number among the pressuredriven 3DTBLs experiments reported in Johnston & Flack (1996), 2) availability of the mean velocity and full Reynolds stress profiles, and of 3) the skin friction (magnitude and orientation) and pressure distribution along the wall. However, some remarks are also in order regarding limitations of the experiment. First, direct wall-stress data is not available, instead, a fit to the log law near + ≈ 100 was used for indirect stress measurement. Second, the description of the zero-pressure gradient region far upstream of the bend region for the purpose of CFD inflow generation is incomplete, therefore requiring an iterative procedure in the inflow generation to match the reported statistics at the first streamwise measurement station.
The experiment featured a spatially developing incompressible turbulent boundary layer, growing along the floor of a square duct with a 30 • bend (see Fig. 1). It should be noted that the boundary layer on the floor was very thin compared to the duct height, with 99 / ranging between 0.026 and 0.07 throughout the test section, where is the width (or height) of the square duct. The flow was far from being fully developed, and the secondary flow near the corner regions was expected to have negligible influence on the centerline region, which was the primary region of interest in the experiment.
Following Schwarz & Bradshaw (1994), two coordinate systems are employed here to facilitate the presentation of the results: ( , , ) denotes the global Cartesian coordinate system; ( , , ) denotes a curvilinear coordinate system aligned with the local duct centerline. = are the wall-normal coordinates (distance from the floor of the duct). In the experiment, the boundary layer on the floor was tripped using a trip wire at the duct inlet located at = 0, thus ensuring a turbulent boundary layer over the entire floor of the test section. Boundary layers on the other three walls of the duct were not tripped (Schwarz, private communication, 2019). Reynolds number is moderately high, with ranging between 4100 and 8500 (or roughly ranging between 1500 and 3900). The flow along the centerline upstream of the bend was reported to exhibit typical characteristics of the canonical 2D zero pressure gradient (ZPG) flat-plate boundary layer. Mean flow threedimensionality was generated in the bend region approximately between = 1626 mm and = 2224 mm due to the cross-stream pressure gradient induced by the bend. The surface streamlines were deflected by up to 22 degrees relative to the local duct centerline. Downstream of the bend, the 3DTBL gradually returned to a 2DTBL owing to the vanished spanwise pressure gradient. The experimental study focused on the boundary layer along the local centerline where the streamwise pressure gradient was found to be small.
The computational domain is identical to the test section in the experiment, which consisted of a square duct ( × = 0.762m × 0.762m) with a total curved length of = 3.748m, as shown in Fig. 1. Two grid resolutions are considered in the present study: a coarse mesh with 8 million control volumes and a fine mesh with 38 million control volumes. Figure 2 shows the near-wall grid spacing distributions along the duct centerline in the two meshes. The computational meshes are designed to maintain adequate wall-modeled LES grid resolution in the test section such that the local boundary layer contains approximately 16∼23 and 32∼45 cells across its thickness in the coarse and fine computational meshes, respectively. Local grid adaptations were applied in the near-wall region with the effect that the grid resolution transitions from the coarser isotropic-cell region in the freestream (Δ = 0.008 m at / > 0.1) toward the finer near-wall region on the duct floor through anisotropic grid refinements. This resulted in wall-parallel grid resolutions (Δ = Δ ) of 4 mm and 2 mm for the coarse and fine meshes, respectively. In the region upstream of the bend ( / 0.43) where the boundary-layer was thin but grew fast, the wall-normal grid spacing were varied with to keep the number of boundary-layer resolving cells approximately constant, resulting in Δ = 0.86 mm∼ 2.2 mm and Δ = 0.43 mm ∼ 1.1 mm for the coarse and fine meshes, respectively. At / 0.43, Δ was fixed at their maximum values aforementioned. Compared to Cho et al. (2021) where WMLES of the same geometry using isotropic voronoi cells was reported, the present study using anisotropic hexahedral cells deploys roughly the same wall-normal resolutions and about twice coarser wall-parallel resolutions in and downstream of the bend. Total cell counts are significantly reduced as a result, while maintaining higher numbers of cells across the thickness of the local boundary layer. The grid-resolution transition zones are located sufficiently away from the shear layer on the duct floor, so that the solution therein is not affected by the accuracy degrade associated with abrupt changes in grid resolution.

Inflow characterization and boundary conditions
Setting the appropriate boundary conditions in the simulation, particularly for the reproduction of flow characteristics upstream of the bend region where a typical equilibrium 2DTBL is expected, is crucial before attempting to compare the simulation results with the experimental results at any downstream location. However, the experiment reports flow statistics at the 22 locations shown in Fig. 1 along the duct centerline, with the first measurement location being far downstream of the test section inlet (at = 826 mm). In the absence of this critical flow information at = 0 mm, where the boundary layer on the floor was tripped in the experiment, we resort to a synthetic turbulence generation based on a digital filter approach (Klein et al. 2003) for approximating the inflow boundary condition, rather than trying to replicate the trip-wire transition in the experiment. This approach requires iterative guesses on the length of the development region (if any) to be appended upstream of the nominal trip location in the experiment ( = 0 mm), and the state of the inflow to be prescribed at the new inlet location. It should be noted that the goal here is to reproduce the 2DTBL upstream of the bend reasonably well, which then acts as the inflow for the 3DTBL within the bend, rather than to exactly match the flow conditions at the test section inlet. After iterating on several inflow conditions and the inlet location, we found that prescribing a flat-plate turbulent boundary layer at = 2560 (Schlatter et al. 2010) at the nominal inlet ( = 0 mm) reproduces the boundary layer statistics well at the first measurement location (station 0: = 826 mm). As shown in Fig. 3, the simulation agrees well with the experiment in terms of the distributions of the boundary layer and momentum thicknesses. In Fig. 4, the mean velocity profile at the first measurement station (station 0: = 826 mm) is also shown to be reproduced well in the present calculations.
The prescription of boundary conditions on the rest of the boundaries is relatively straightforward. A subsonic Navier-Stokes characteristic boundary condition (Poinsot & Lele 1992) is imposed at the outlet of the duct. No attempt was made to resolve the boundary layers on the two side walls and the top wall which were not tripped in the experiment. The no-slip boundary condition is applied to each of these walls. The wall model is applied to the bottom wall, and the wall stress calculated from the wall-model solution is used as the Neumann boundary condition on this wall. All walls are assumed to be thermally adiabatic.

Flow solver and SGS / near-wall modeling
The simulations were performed with CharLES, an unstructured cell-centered finite-volume compressible LES solver developed at Cascade Technologies, Inc. The solver employs an explicit third-order Runge-Kutta (RK3) scheme for time advancement and a second-order central scheme for spatial discretization. More details regarding the flow solver can be found in Khalighi et al. (2011) and Park & Moin (2016a). The Vreman model (Vreman 2004) is used to close the SGS stress and heat flux.
In WMLES, LES equations are solved on a coarse mesh, where the stress-carrying eddies in the near-wall region are mostly unresolved. The LES mesh alone cannot represent the sharp velocity gradients and the momentum transport near the wall. This causes SGS models to produce insufficient levels of modeled stresses. Wall modeling aims to compensate for such numerical and modeling errors in the underresolved near-wall region of LES, by augmenting the total stresses directly through the imposition of the modeled stress boundary condition at the wall in lieu of the no-slip condition. In the present work, wall models solve simplified, vertically integrated, or full RANS equations on a separate near wall mesh. The grid for wall models have fine resolution in the wall-normal direction (with the exception of the integral nonequilibrium wall model, which does not require a wall-normal grid), but the wall-parallel grid resolution (if any) is identical to or coarser than the LES grid. All wall models in this study are driven by the LES states imposed at their top boundaries, which are taken at a specified matching height in the LES grid. At each time step, wall stress and heat flux obtained from the wall model are used as the Neumann wall boundary condition for the LES. Figure 5 shows a schematic of the wall modeled LES procedure employed in the present work.
In the present work, the location at which the wall-models take input from the LES (matching height, denoted as ℎ ( )) is fixed across different grid resolutions by setting it to the centroids of the third off-wall cells or the top faces of the fifth off-wall cells in the coarse and fine meshes, respectively, corresponding to 10 3 ℎ / = 6 ∼ 14. This has effect of fixing the wall-modeled regions in LES during grid refinements, so that improvement in LES prediction is not associated with the change in wall-modeling details, but it is attributed The three wall models considered in the present study are: an equilibrium stress model (EQWM) in the form of ordinary differential equations (ODE), an integral nonequilibrium wall model (integral NEQWM) that solves the vertically-integrated Navier-Stokes equations, and a PDE nonequilibrium wall model (PDE NEQWM) that retains the complexity of the full Navier-Stokes equations. All three wall models parameterize the unresolved turbulence in the wall-model domain in a statistical sense using simple RANS models based on the mixing-length formulation. Note that the EQWM and PDE NEQWM have previously been implemented in CharLES, and they were tested extensively through various studies (Park & Moin 2014, 2016aPark 2017;Bodart & Larsson 2012). The integral NEQWM was recently integrated into CharLES, the implementation aspects of which will be discussed in a future article (Hayat & Park 2021). A brief description of each of these models is given below.
The EQWM (Kawai & Larsson 2012; Bodart & Larsson 2011) solves the simplified boundary layer equations which account only for the wall-normal diffusion.
where is the local wall-normal coordinate, | | is the wall-parallel velocity magnitude, is the temperature, is the molecular viscosity, is the molecular thermal conductivity, and and are the turbulent eddy viscosity and conductivity, respectively. The velocity vector is assumed to be aligned with the LES velocity at the matching height. Owing to this intrinsic assumption, the equilibrium wall model is incapable of predicting skewed velocity profiles within the wall-modeled domain. Also, due to the unidirectionality and the condition + > 0, the EQWM can represent monotonic velocity profiles only, and it cannot predict velocity profiles with sign changes in the slope as found in the near-wall regions of separated shear layers. The wall-model eddy viscosity is based on the following mixing-length formula, On the other hand, the PDE NEQWM (Park & Moin 2014, 2016a) solves the full 3D unsteady RANS equations, (2.5) where is the density and is the velocity component, is pressure and = /[ ( −1)] + /2 is the total energy. The stress tensor and heat flux are given by, = 2( + ) and = −( + ) . For the RANS closure, a novel mixing-length model is used, which dynamically accounts for the resolved Reynolds stresses carried by the wall model (Park & Moin 2014). The wall-model mesh for the PDE NEQWM has the same wall-parallel grid content as in the coarse LES mesh, but it is refined in the wall-normal direction to resolve the viscous sublayer.
Lastly, the integral NEQWM formulation solves a similar set of equations as the PDE NEQWM, albeit in a wall-normal integrated form. Currently, this formulation is limited to incompressible flows, and therefore the energy equation is not solved. For the sake of brevity, only the 2D formulation (the wall-normal and one wall-parallel velocity components) is presented below, and the reader is referred to Yang et al. (2015) for the details of full 3D formulation. The vertically integrated momentum equation is given by, (2.7) where and represent the local wall-parallel and wall-normal coordinates, ℎ is the matching height, is the time-filtered velocity from the LES solution at the matching location. = =0 and ℎ = ( + ) =ℎ are the shear stresses at the wall and at the matching location, respectively. The integral terms are evaluated by assuming an analytical composite profile for the velocity within the wall model, which has the form: where the unknown parameters , , and are determined from the solution of equation (2.7) along with suitable matching and boundary conditions. For the full 3D formulation consisting of two wall-parallel velocity components, like the one employed for the present study, the composite profiles have a total of 11 unknown parameters. This approach attempts to model the effects of pressure gradient and advection through the last term in Eq. (2.9) representing linear departure from the log law.
It is worth mentioning here that in the original 3D formulation in Yang et al. (2015), the assumed form of the viscous-sublayer velocity profiles in the two wall-parallel directions (Eq. (C5) in Yang et al. (2015)) resulted in inconsistent asymptotic behavior of velocity near the wall. This made the wall-stress predictions of the wall model highly sensitive to the choice of the local / coordinates. In our current integral NEQWM formulation, we modify the assumed viscous-sublayer profile to ensure consistent near-wall asymptotic behavior as given by the Taylor series expansion. The details of this modified formulation along with its implementation in an unstructured solver are presented in Hayat & Park (2021). The MATLAB implementation for the EQWM and the integral NEQWM are available on GitHub at https://github.com/imranhayat29/Wall-Models-for-LES.
A remark is in order regarding the overall cost of simulations with different wall models. The computational costs of the three wall models were compared by running the simulations on the fine LES mesh with 256 CPU cores for three convective flow-through times. When the cost of the simulation without any wall model (no-slip LES) is taken as the unity, the simulation costs are 1.27, 1.2, and 2.2 with the EQWM, the integral NEQWM, and the PDE NEQWM, respectively. The higher cost with the PDE NEQWM is due to the inversion of a large linear system required as a part of implicit time advancement.

RESULTS
Results from the WMLES simulations are discussed in this section. Overall characteristics of the flow are first highlighted from the instantaneous flow field standpoint. Mean and turbulence statistics obtained with different wall models are then assessed against the experimental data. Furthermore, the three-dimensionality of the outer portion of the boundary layer is examined with the aid of Reynolds stress anisotropy and the Johnston triangular plot.
Some remarks are in order concerning the ways the main results are presented in this paper. The primary interest of the experiment was to examine the effect of the mean three dimensionality in the absence of strong streamwise pressure gradient. To this end, the experiment presented the key flow statistics along the floor centerline, where the axial pressure gradient was observed to be nearly zero. It should be noted, however, that the meanflow trajectory near the wall deviates somewhat substantially from the centerline as the flow passes through the bend region, as observed from the instantaneous flow field (Fig. 7) and the near-wall streamlines (Fig. 10( )). This leaves some ambiguity in the interpretation of the statistics presented along the centerline, because any two fluid particles on the centerline (in and after the bend region) would have traveled along different Lagrangian trajectories, experiencing different history effects (most notably, they are subject to different upstream axial/spanwise pressure gradients.) With this limitation in mind, we still choose to show our results along the centerline, as all experimental data (except the wall-pressure) were presented so. Figure 6 shows the mean flow statistics and the Reynolds stresses at = 1875 mm (the eighth measurement station in Fig. 1), as well as the centerline distribution obtained from the EQWM LES with the coarse and the fine grids described in Sec. 2.1. Although only the EQWM results are shown here for brevity, it is noted that the other two wall models exhibited similar grid-convergence characteristics. In Fig. 6( ), both the mean velocity and the mean flow direction (see Sec.3.3 for definition) profiles obtained from the coarse-grid calculation are already in good agreement with the experiment, and the results only improve marginally with the grid refinement. More importantly, this points toward the grid convergence of the results for the refinement level used in this study. Figure 6( ) shows the skin-friction distribution along the centerline of the duct. Between the first and the last measurement stations, we observe a reasonably converged , with the fine-grid calculation producing slight improvement in . In Fig. 6( ), a similar trend is observed for all the Reynolds stress components shown, with the exception of the streamwise component of the Reynolds normal stress, which shows noticeable variation with grid refinement in the near-wall region; however, the Reynolds stresses in the outer portion of the boundary layer have largely converged. Having established reasonable evidence of grid-convergence for most of the flow statistics on the coarse grid, the remainder of this paper will focus largely on discussing the results obtained with the coarse grid unless stated otherwise. Figure 7 visualizes the vortical structures in the floor boundary layer. Near the inlet (approximately within 20 times the inlet boundary-layer thickness from the inlet), structures with less coherence resulting from the synthetic inflow turbulence generation are observed. The floor boundary layer then gradually develops into a coherent fully developed state far upstream of the bend region, which is also verified by the velocity profile at = 826mm ( Fig. 4) following the typical law of the wall observed in 2D turbulent boundary layer. When the flow enters the bend region, a clear contrast of two boundary layers with different origins (blue and red regions) are observed. The boundary layer in the red region is thicker than that in the blue region, showing that a new boundary layer is emerging from the concave sidewall (at < 0) and the original boundary layer developed from the upstream section is turning toward the convex side. This overall flow behavior is visually in fair agreement with the surface oil visualization in the experiment as shown in the inset in Fig. 7.

Mean flow statistics
The cross-stream pressure gradient is the source of the mean three dimensionality in the bend region. It acts to deflect the streamlines close to the wall more strongly than those near the free stream. It is therefore important to first establish close agreement in the pressure distribution close to the bend between the simulation and the experiment. Note that in the current case without flow separation, the pressure distribution is determined largely by the wall geometry and the inviscid effect, presumably unaffected by the wall-modeling details. Figure 8 shows the distribution of the static wall-pressure coefficient on the floor of the duct.   (3.1) Following the experiment, is the static pressure at = 0 mm, and is the freestream velocity at = 826 mm (defined at the spanwise centerline). The wall-pressure probing lines are parallel to the duct centerline as shown in the inset figure. The figure shows good agreement between the simulation and the experiment, except in the recovery region downstream of the bend. The axial pressure gradient is almost zero along the centerline. On the other hand, a significant spanwise pressure gradient starts to develop upstream of the bend, reaches a maximum within the bend region, and eventually decays to zero downstream of the bend. The reason for disagreement in the recovery region remains unclear to us. While in the experiment remains to be slightly negative near the outlet, in the simulation naturally vanishes to its upstream zero value as the flow relaxes back to its 2D ZPG state. Note that the experiment reported only the centerline distribution in this region, and that extending the duct further downstream in the simulation did not change the trend. Figure 9 shows the distribution of the skin-friction coefficient along the duct centerline. The skin-friction coefficient is defined as, , is the magnitude of the mean wall-shear stress vector, and is the local free-stream velocity. Figure 9( ) shows the centerline distribution of the skin friction upstream of the bend as a function of the momentum thickness Reynolds number. The centerline mean flow is expected to be agreeing well with the canonical ZPG 2DTBL in this region. A deviation of the skin friction from the ZPG 2DTBL near the inlet is the artifact of the inflow treatment. Note that the synthetic inflow turbulence generation methods when applied to DNS or wall-resolved LES of low Reynolds number are known to produce a development length of 10 ∼ 20 initial boundary layer thicknesses ( ), through which coherence-lacking artificial structures mature into fully-developed turbulence (e.g., Patterson et al. (2021); Sandberg (2012); Larsson (2021) where = 400 ∼ 500). The present high-Reynolds number case simulated with very coarse meshes produced longer development lengths (30∼40 ). The flow was observed to be fully developed from slightly upstream of the first measurement station, after which the WMLES results are in reasonable agreement with the experiment as well as with an empirical correlation (Fernholz & Finley 1996) and a wall-resolved LES of ZPG 2DTBL (Eitel-Amor et al. 2014). In Fig. 9( ), slight overprediction of the skin friction from WMLES is observed throughout the duct. A similar trend was reported by Cho et al. (2021), where the EQWM was used with up to 76 million control volumes in LES. It should be noted that the wall shear stresses were measured indirectly in the experiment using a Preston tube and Patel's calibration (Patel 1965) for a 2DTBL. Patel (1965) reported that errors as large as 6% could occur when Preston tube is used in flows with moderate streamwise pressure gradient.
Next, we examine the mean three-dimensionality of the flow in the duct. The variation of the surface flow direction relative to the freestream direction is a measure of the mean flow three-dimensionality. As shown in Fig. 10( ), the crossflow is almost zero in the upstream, and it grows rapidly as the flow approaches the bend. The resulting turning angle reaches the maximum near the end of the bend and decays gradually thereafter. These observations are consistent with the development of the spanwise pressure gradient. All three wall models predict the general trend in the turning-angle variation correctly; however, the PDE NEQWM gives the most accurate prediction among the three (especially within the bend region), followed by the integral NEQWM, and then the EQWM. The maximum difference between the PDE NEQWM and the EQWM is roughly 5 degrees occurring at / = 0.49 within the bend. Note that the total flow turning is an accumulative effect of the local flow change, and the area under the curve in Fig. 10( )   near-wall flow direction by different wall models in terms of the select surface streamlines calculated from the mean wall shear-stress vector. It can be clearly seen that the flow deviates from the local centerline; however, the deviation is not predicted evenly across the different wall models, consistent with our observation in the flow turning angles in Fig. 10( ). The surface flow from the PDE NEQWM turns much more rapidly than that from the EQWM. This has a direct implication for practical engineering flows involving 3DTBL. For instance, in the skew-induced 3DTBL over the swept wing of an aircraft, the circulation over the wing and the downwash characteristics could potentially be affected by the near-wall change in the flow direction, which may alter the lift. Additionally, errors in the surface-flow directionality can alter three-dimensional separation patterns observed near the wing-body juncture of wings at high-lift configurations (Evans et al. 2020) to increase uncertainty in prediction of maximum lift and the associated drag.
In Fig. 11, profiles of the mean-velocity magnitude are compared between the different wall models and the experiment, at several locations along the centerline, including upstream of, within, and downstream of the bend. It can be seen that the no-slip LES, which does not employ a wall model, gives a very poor prediction of the mean velocity. Here, a higher momentum is imparted to the boundary layer as a consequence of the underpredicted wall shear force. With the introduction of wall modeling, a significant improvement is achieved in the predicted mean-velocity profiles. In line with the predictions of the skin-friction coefficient in Fig. 9, the mean velocity profiles across the three wall models are almost identical. Note that here, the profiles only show the magnitude of the mean velocity, thus lacking information on the three-dimensionality of the mean flow.
To complete this picture, Fig. 12 shows how the flow direction changes along the wallnormal direction. The flow direction characterizes the mean flow three-dimensionality and it is defined by the angle between the mean velocity vector and the freestream velocity vector. The mean flow three-dimensionality is the strongest at the wall and becomes weaker with the increase in distance from the wall, as evident from the diminishing crossflow away from the wall. A difference of approximately 3 degrees is observed between the WMLES and the experimental results. However, the predicted angles from the different wall models are almost identical, indicating that the difference in the wall-model outputs (the wall-shear force direction observed in Fig. 10( )) is not felt by the LES solutions away from the wall.

Reynolds stress
We now turn our attention to the turbulent content of the 3DTBL and its role in distinguishing this flow from the canonical 2DTBL, by examining the Reynolds stress related statistics. Indeed, the Reynolds stresses in the 3DTBL exhibit unique characteristics not seen in the 2DTBL, as we will see shortly. Figure 13 shows the profiles of the Reynolds normal stresses at the same five centerline locations which were chosen to depict the mean velocity profiles. Note that the experiment did not have access to the Reynolds-stress data in the inner layer, therefore missing information on the peak values and their locations. The no-slip LES acutely underpredicts the Reynolds normal stress, pointing towards the grid resolution being insufficient for the no-slip LES to resolve the near wall eddies. All three wall models significantly improve the prediction of the normal stresses, bringing the profiles closer to the experimental results. The predicted Reynolds normal stresses in the wall-parallel directions show remarkable agreement with the experiment, whereas those in the wall-normal direction are underpredicted near the wall. Figure 14 shows the profiles of the Reynolds shear stresses, where substantial improvement with wall modeling is also observed.
An important characteristic of the 3DTBL is that the Reynolds shear stress vectors are not necessarily aligned with the mean velocity gradient vectors, which challenges the fundamental assumption of the commonly used isotropic eddy viscosity model. The directions of these two vectors are characterized by the angles defined below, between the two vectors also decreases ( Fig. 15( )). Furthermore, this lag appears to be a function of the distance from the wall. The experiment shows that the lag decreases with wall distance in the outer layer above / 99 ≈ 0.7 within the bend. Downstream of the bend, the Reynolds shear stress vector even starts to lead mean velocity gradient vector above / 99 ≈ 0.8, and this lead increases with wall distance. These behaviors and the shear-stress angles therein are not captured well in LES. We conjecture that computation of the angles in this region is prone to contamination by numerical error or measurement noise, because both the Reynolds shear stress and the mean velocity gradient values are very small there. At / 99 0.7, a reasonable agreement between the simulations and the experiment is observed, and results from different wall models do not show notable difference.

Lumley triangle: anisotropy invariant map
We have noted that the Reynolds stresses from the simulations agree reasonably well with the experiment. This makes it possible to further investigate the anisotropy of the Reynolds stress in the outer layer of this 3DTBL using the WMLES solution. In this section, using the fine grid prediction, we employ the Lumley triangle to analyze the Reynolds-stress anisotropy.  Following Pope (2000), the normalized anisotropy tensor is defined as, The anisotropy tensor has zero trace and thus has two independent principal invariants. It is convenient to define two variables and , corresponding to the two invariants, as 6 2 = −2 2 = 2 = (3.5) (3.6) The state of anisotropy can be characterized by the above two variables and . All realizable Reynolds-stress states must be located within a triangular region in the -plane, which is known as the Lumley triangle.  -Durán et al. 2020) can be also made from these data shown in Fig. 16( ). Right upstream of the bend, the wall-normal distribution of the anisotropies away from the wall shows some similiarty to that in the 2D channel (Fig. 16( )), exhiting a characteristic S-shape lying close to the axisymmetric-expansion (AE) limit. As the flow passes through the bend, the left cusp of this S curve rapidly dislocates toward inside the triangle, leaving less points close to the AE limit. The non-monotonic decrease of the anisotropies (with increasing wall distance) is also observed vividly, which is only weakly present in the 2D channel. While the station 18 is considerably downstream of the bend region, the anisotropies are still seen further departing from its 2D behavior. This is consistent with the observation in Fig. 15 that the Reynolds stresses respond more slowly than the mean to the imposed three-dimensionality. Also, in Fig. 16( ), note the similarity of the anisotropy distributions in the duct and the shear-driven 3DTBL from the transient channel flow, although departure from the 2D behavior is much stronger in the duct case.

Triangular plot
In the present work, the mean three-dimensionality in the outer layer is created by the inviscid skewing mechanism, where the streamwise vorticity is produced by reorientation of the spanwise vorticity. A popular way of representing the crossflow so developed is the "Johnston triangular plot" (Johnston 1960), which is the triangular plot of against , where and are along and normal to the local freestream direction, respectively. In particular, the outer-layer mean velocity profile of the skew-induced 3DTBLs can be accurately approximated by the Squire-Winter-Hawthorne (SWH) relation (Squire & Winter 1951;Hawthorne 1951;Bradshaw 1987 where is the freestream turning along the streamline. This is a special case of the vorticity transport equation in which viscous terms and Reynolds stresses are neglected. The SWH relation shows up as a straight line with a negative slope toward the right end of the triangular plot. In Fig. 17, we present the mean velocity for the current bent duct flow and a temporally developing shear-driven 3D channel flow from Lozano-Durán et al. (2020) in the triangular plot. It is observed that the mean velocities in the outer layer from the duct flow satisfy the SWH formula well, whereas they deviate from the SWH relation in the shear-driven case, as expected. The slope in the SWH relation represents the freestream turning angle with respect to the upstream flow, and the freestream slope in the triangular plot (Fig. 17( )) therefore increases toward the downstream direction. The wall model solutions (defined only close to the wall) are also visualized along with the outer LES solution in Fig. 18. It gives us a vivid picture of different wall models' capabilities to depict skewed mean-velocity profiles. In the triangular plot, the wall model solutions are from the origin to the 3rd cell LES solutions (coarse grid resolution). The EQWM, due to its unidirectional assumption, cannot describe skewed mean-velocity profiles, and it shows up as a straight line starting from the origin in the triangular plot. On the other hand, the PDE NEQWM and the integral NEQWM are able to represent skewed mean-velocity profiles. They show up as curved lines in the triangular plot, implying that the flow direction changes with the wall distance. During the crossflow developing stage (Fig. 18( )), the difference between the two NEQWM solutions and the EQWM solution gradually grows. Notice that the PDE NEQWM is able to represent a richer variation of the slope (flow direction) along the wall-normal direction than the integral NEQWM. During the crossflow decaying stage

Quantification of the nonequilibrium contributions to the wall shear stress direction
In this section, the nonequilibrium effects neglected in the EQWM are analyzed through the full RANS equations used in the nonequilibrium wall models. This analysis highlights the importance of accounting for the nonequilibrium effects in wall modeling for accurate prediction of the surface flow direction, as well as the subtle difference in how these effects are incorporated in different wall models. Our analysis is based on the solutions of the PDE NEQWM and the integral NEQWM. Ideally, this analysis should be done in a priori sense utilizing the fully resolved flow fields, as attempted in Hickel et al. (2012). This was deemed infeasible in the present case due to the high Reynolds number.
Assuming incompressible flow, the time-averaged momentum equations in the PDE NEQWM can be recasted as where and are the LES velocity components at the matching location. A similar expression can be found in Wang & Moin (2002). As shown in the previous sections, and are almost identical among the simulations with the different wall models. The wall shear stress direction, which is the quantity of interest exhibiting the most significant difference among the three wall models, is then expressed as ). The fidelity with which the constitutive terms of and are modeled is, therefore, crucial to the performance of wall models in predicting the surface flow direction.
To separate the nonequilibrium contributions from the flow direction predicted by the EQWM, we can first reorganize Eqn. For the present flow, it is shown in Fig. 19 that is relatively small compared to , which permits the use of truncated Taylor series expansion of (3.21) The terms / and / can be viewed as representing the corrective contributions from the nonequilibrium effects to the surface flow direction predicted by the EQWM. In essence, Eqn. (3.21) enables the surface flow angle to be decomposed into distinct contributions originating from the equilibrium and the nonequilibrium terms. Although not shown here for brevity, this truncated relation was found to provide almost identical description as the actual wall flow direction ( , / , ). Note that Eqn. (3.21) applies to the integral NEQWM as well. However, the terms and therein are largely modeled using the assumed velocity profile, in contrast to the PDE NEQWM where these terms are largely solved for, using the full RANS equations.
To further compare the wall models, we plot the two leading order terms of Eqn. (3.21) for the two nonequilibrium wall models in Fig. 20 (all nonequilibrium terms are zero for EQWM and they are not plotted.) Several interesting observations are made. First, the two nonequilibrium wall models show that the total nonequilibrium angle corrections are large at the beginning of the bend region, and that they gradually decrease to zero towards the end of the duct. That is, the nonequilibrium models properly sense the region where nonequilibrium effects are important, and attempt to model them therein. Second, the two nonequilibrium wall models produce comparable distributions of / , implying that the axial ( ) contents of the nonequilibrium effects are modeled almost identically by the two wall models. Third, the difference between the wall models appears to be concentrated in − / , i.e., in the way the models sense and model the cross-flow ( ) component of the nonequilibrium effects. The integral NEQWM underpredicts − / throughout the duct compared to the PDE NEQWM. Furthermore, within the bend, the signs of this term are opposite in the two wall models. This term is comprised of advection ( ), pressure gradient ( ), lateral diffusion ( ), and . Among these terms, and are imposed largely from the LES solution (which are seen to be identical from the simulations using the two wall models), and is seen to be negligible in its magnitude. Therefore, we conjecture that the difference of − / in the integral wall model originates largely from its assumed velocity profile, which is directly used in computing the advection term .
Furthermore, the individual contribution from different nonequilibrium effects are also analyzed. Starting from Eqn. (3.18) with = = 0 (corresponding to the EQWM prediction), we examine the change in the surface flow turning angle by systematically including the contributions from various nonequilibrium effects (advection, pressure gradient, and lateral diffusion) to and . Note that this is done in a post-processing manner using the solution of the PDE NEQWM. The results are shown in Fig. 21. When all the nonequilibrium effects are taken into account, the reconstructed surface flow turning-angle agrees with the PDE NEQWM prediction, as expected. The lateral diffusion terms have negligible contributions to the surface flow turning angle. The pressure gradient and advection terms are significant within the bend, where the mean three-dimensionality develops. However, these terms appear to have a competing effect within the bend. The pressure gradient tends to make the flow deviate more from the local freestream, while the advection tends to make the flow deviate less from the local freestream. These two contributions largely cancel out each other, but a subtle balance between the two appears to be crucial in prediction of the surface flow direction.

Conclusion
In the present work, we have studied a spatially-developing pressure-driven 3DTBL over the floor of a square duct with a 30 • bend using WMLES. The major focus of this study is to contrast the performance of three commonly used wall models in a high Reynolds number pressure-driven 3DTBL.
These models represent varying degrees of physical fidelity: the EQWM solves the simplified boundary-layer equations which neglect all nonequilibrium effects and assume the flow is unidirectional in the wall-modeled region; the PDE NEQWM solves unsteady 3D RANS equations which maintain much of the nonequilibrium effects; the integral NEQWM represents a compromise between the two models. Algebraic complexity is kept low thanks to the presumed velocity profile, while the nonequilibrium effect is represented by the linear perturbation to the logarithmic law of the wall.
The mean flow statistics and the Reynolds stresses are predicted reasonably well by the WMLES under the current grid resolutions, which are too coarse for the no-slip LES. The more comprehensive PDE NEQWM does show improvement against the integral NEQWM (which is in turn better than the EQWM) in predicting the direction of mean wall shear stress. The error in the local wall-shear force direction accumulates along the surface streamlines, leading to significant difference of the surface flow at the end of the duct. Budget analyses have been conducted to elucidate precise mechanisms by which the three wall models produce different predictions of the wall-shear stress directions given almost identical inputs. For the present flow, the surface flow direction is shown to have separable contributions from the equilibrium part of the wall models and the integrated nonequilibrium effects (advection and pressure transports). It was shown that difference in how the cross-flow component of the nonequilibrium contribution is modeled leads to different behaviors of the models.
Additionally, the pressure gradient and the advection are shown to have a competing effect in deflecting the surface flow within the bend. Although these terms largely cancel out each other, neglecting any of them produces large errors in the surface flow direction, and a subtle balance between the two appears to be crucial in prediction of the surface flow.
However, such difference in the wall shear stress direction predicted by the wall models appears to be not felt by the (outer) LES solution. The three wall models produce almost the same mean velocity and Reynolds stresses profiles. A possible explanation for this phenomenon comes from the nature of the pressure-driven 3DTBL. In the duct flow under consideration, the mean flow three-dimensionality in the outer layer is largely controlled by the "inviscid skewing" mechanism, which is not affected by the near wall viscous effects. In this class of 3DTBL, the outer-layer mean flow appears to be robustly set up by the inviscid effect, provided that the momentum drain by the wall is specified with reasonable accuracy only (in particular, its magnitude rather than the direction).
The characteristics of the 3DTBL are also analyzed. The anisotropy of turbulence along the wall-normal direction is investigated with the Lumley triangle. Compared to the 2DTBL, a large in-ward pointing sharp corner is presented in the Lumley triangle plot of 3DTBL in the downstream section. This large sharp corner represents a non-monotonic decrease of anisotopy for increasing wall distance. When the mean crossflow is generated by the inviscid effect (reorientation of spanwise vorticity), the relation between spanwise and streamwise velocity in the local freestream coordinate system in the outer part of the boundary layer shows as a straight line predicted by the SWH formula. In the present duct flow, the outer LES results show good agreement with the SWH formula, which shows that the"inviscid skewing" mechanism is the major effect for generating mean three-dimensionality in the outer layer. The triangular plot of the inner wall-model solution reveals that the PDE NEQWM is most capable of representing the change in the flow direction along the wall-normal direction. The EQWM is most restrictive in this sense with its unidirectional flow assumption. The integral NEQWM is in between the two.

Acknowledgement
This research was sponsored by NASA's Transformational Tools and Technologies Project of the Transformative Aeronautics Concepts Program under the Aeronautics Research Mission Directorate (Grant 80NSSC18M0155), and by the Office of Naval Research (ONR) (Grant N000141712310). Computational resources supporting this work were provided by the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center.