3.1. Grid convergence

Figure 6 shows the mean-flow statistics and the Reynolds stresses at $x'=1875$ mm (the eighth measurement station in figure 1), as well as the centreline skin-friction coefficient $C_{f}$ distribution obtained from the EQWM LES with the coarse and the fine grids described in § 2.1. The skin-friction coefficient is defined as

(3.1a,b)\begin{equation} C_{f}=\frac{\tau_{w}}{\frac{1}{2}\rho U^2_{e}}, \quad \frac{U_e}{U_{ref}}=(1-C_{p})^{1/2}, \end{equation}

where $\tau _w = \sqrt {\tau _{w,x}^2+\tau _{w,z}^2}$ is the magnitude of the mean wall shear stress vector and $U_e$ is the local free stream velocity. The pressure coefficient is defined as

(3.2)\begin{equation} C_{p} = \frac{p-p_{ref}}{\frac{1}{2}\rho U^2_{ref}}. \end{equation}

Following the experiment, $p_{ref}$ is the static pressure at $x'=0$ mm and $U_{ref}$ is the free stream velocity at $x'=826$ mm (defined at the spanwise centreline). 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 figure 6(a), both the mean velocity and the mean-flow direction (defined by the angle between the mean velocity vector and the free stream velocity vector) 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 towards the grid convergence of the results for the refinement level used in this study. Figure 6(b) shows the skin-friction distribution along the centreline of the duct. Between the first and the last measurement stations, we observe a reasonably converged $C_{f}$, with the fine-grid calculation producing slight improvement in $C_{f}$. In figure 6(c), 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 6.

Grid convergence of the EQWM LES. (a) Blue, mean-flow direction versus distance from the wall; red, mean velocity magnitude profile. (b) Centreline distribution of the skin-friction coefficient and (c) Reynolds-stress profiles. Red, $\overline {u'u'}$; blue, $-\overline {u'v'}$; green, $\overline {v'w'}$. Solid lines, coarse grid; dashed line, fine grid. Profiles in panels (a) and (c) are at $x'=1875$ mm (station 8). Black vertical dashed lines in panel (b) denote the start and end of the bend region.

3.2. Instantaneous flow field

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 $x'=826$ mm (figure 4) following the typical law of the wall observed in a 2-D 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 within the bend region (at $z<0$) and the original boundary layer developed from the upstream section is turning rapidly towards the convex side (at $z>0$). This overall flow behaviour is visually in fair agreement with the surface oil visualization in the experiment as shown in the inset in figure 7.

Figure 7.

Visualization of the near-wall vortical structures using the isosurfaces of $Q$ ($Q = Q_{rms}$, where $Q$ is the second invariant of the velocity gradient tensor), coloured by the distance from the floor of the duct. Flow is from the bottom left to the top right. Surface oil visualization from the experiment (Schwarz & Bradshaw Reference Schwarz and Bradshaw1994) is shown in the inset (figure reprinted with permission from Schwarz & Bradshaw (Reference Schwarz and Bradshaw1994).

3.3. 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-modelling details. Figure 8 shows the distribution of the static wall-pressure coefficient on the floor of the duct. The wall-pressure probing lines are parallel to the duct centreline, 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 centreline. However, 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 $C_p$ in the experiment remains to be slightly negative near the outlet, $C_p$ in the simulation naturally vanishes to its upstream zero value as the flow relaxes back to its 2-D ZPG state. Note that the experiment reported only the centreline distribution in this region and that extending the duct further downstream in the simulation did not change the trend.

Figure 8.

Variation of the wall-pressure coefficient from coarse EQWM simulation (results from the other wall models and resolutions are almost identical). Symbols are from the experiment. Colours denote different spanwise locations corresponding to the inset figure. Here, $z':0$ mm, ${\pm }127$ mm, ${\pm }254$ mm, ${\pm }368$ mm.

Figure 9 shows the distribution of the skin-friction coefficient along the duct centreline. The centreline mean-flow is expected to agree well with the canonical ZPG 2DTBL in this region. A deviation of the skin friction from the ZPG 2DTBL near the inlet is the artefact 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 ($\delta _{in}$), through which coherence-lacking artificial structures mature into fully developed turbulence (e.g. Patterson, Balin & Jansen (Reference Patterson, Balin and Jansen2021), Sandberg (Reference Sandberg2012), Larsson (Reference Larsson2021) where $Re_\tau = 400 \sim 500$). The present high-Reynolds-number case simulated with very coarse meshes produced longer development lengths ($30\unicode{x2013}40 \delta _{in}$). 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 Reference Fernholz and Finley1996) and a wall-resolved LES of ZPG 2DTBL (Eitel-Amor, Orlu & Schlatter Reference Eitel-Amor, Orlu and Schlatter2014). In figure 9(b), slight overprediction of the skin friction from WMLES is observed throughout the duct. A similar trend was reported by Cho et al. (Reference Cho, Lozano-Durán, Moin and Park2021), 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 Reference Patel1965) for a 2DTBL. Patel (Reference Patel1965) reported that errors as large as 6 % could occur when a 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 free stream direction is a measure of the mean-flow three-dimensionality. As shown in figure 10(a), the cross-flow 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^{\circ }$ occurring at $x'/L=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 figure 10(a) can be thought of as an approximation of the near-wall total flow turning angle. Related to this, figure 10(b) visually highlights predictions of the 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 centreline; however, the deviation is not predicted evenly across the different wall models, consistent with our observation in the flow turning angles in figure 10(a). The surface flow from the PDE NEQWM turns much more rapidly than that from the EQWM. In figure 11, profiles of the mean-velocity magnitude are compared between the different wall models and the experiment, at several locations along the centreline, 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 modelling, a significant improvement is achieved in the predicted mean-velocity profiles. In line with the predictions of the skin-friction coefficient in figure 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.

Figure 9.

Centreline distribution of the skin-friction coefficient ($C_f$). (a) $C_f$ versus $Re_{\theta }$ upstream of the bend. (b) $C_f$ versus axial location. Squares, experiment; red dash–dotted line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM. In panel (a): black dotted line, ZPGFPBL empirical correlation ((9) in Fernholz & Finley Reference Fernholz and Finley1996); black dashed line, WRLES of ZPGFPBL (Eitel-Amor et al. Reference Eitel-Amor, Orlu and Schlatter2014).

Figure 10.

(a) Centreline distribution of the surface flow turning angles with respect to the free stream ($\gamma _w = \tan ^{-1}(\tau _{w,z}/\tau _{w,x})$ is the wall shear stress direction, $\gamma _{\infty } = \tan ^{-1}(W_e/U_e)$ is the free stream direction). (b) Streamlines of wall shear stress. Squares, experiment; red line, equilibrium wall model; blue line, PDE non-equilibrium wall model; green line, integral non-equilibrium wall model. Solid and dashed lines are for coarse and find grids, respectively.

Figure 11.

Profiles of the mean-velocity magnitude at five measurement locations (stations 4, 8, 12, 16 and 20, from left to right). Station 4 is upstream of the bend; station 8 is within the bend; stations 12, 16 and 20 are downstream of the bend. Red dash–dotted line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM; magenta dotted line, no-slip LES; black circle, experiment. Profiles are shifted along the abscissa by 1.

To complete this picture, figure 12 shows how the flow direction changes along the wall-normal direction. 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 cross-flow away from the wall. A difference of approximately $3^{\circ }$ 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 figure 10b) is not felt by the LES solutions away from the wall.

Figure 12.

Mean-flow direction relative to the local free stream as a function of wall distance. Red dash–dotted line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM; magenta dotted line, no-slip LES; black solid line, experiment. Symbols are used to differentiate stations along the duct floor centreline only (square, station 0; triangle, station 6; diamond, station 10). Lines with the same symbols denote the results at the same stations.

3.4. 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 centreline 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 modelling is also observed.

Figure 13.

Reynolds normal stress profiles at the same five measurement stations as the mean velocity profiles in figure 11. Red dash–dotted line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM; magenta dotted line, no-slip LES; black circle, experiment. Profiles are shifted along the abscissa by multiples of 10, 2 and 4 for $\overline {u'u'}$, $\overline {v'v'}$ and $\overline {w'w'}$, respectively.

Figure 14.

Reynolds shear stress profiles at the same five locations as the mean velocity profiles in figure 11. Red dash–dotted line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM; magenta dotted line, no-slip LES; black circle, experiment. Profiles are shifted along the abscissa by multiples of 1.5, 2 and 1 for $-\overline {u'v'}$, $-\overline {u'w'}$ and $\overline {v'w'}$, respectively.

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 as

(3.3a,b)\begin{equation} \gamma_{\tau}=\tan^{-1}\left(\frac{\overline{v'w'}}{\overline{u'v'}}\right), \quad \gamma_{g}=\tan^{-1}\left(\frac{\partial W/\partial y}{\partial U/\partial y}\right). \end{equation}

Figure 15(a) clearly shows that the Reynolds shear stress vector lags behind the mean velocity gradient vector within the bend where mean-flow three-dimensionality is strongest. Downstream of the bend where mean-flow three-dimensionality declines, the difference between the two vectors also decreases (figure 15b). 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 $y/\delta _{99} \approx 0.7$ within the bend. Downstream of the bend, the Reynolds shear stress vector even starts to lead to a mean velocity gradient vector above $y/\delta _{99} \approx 0.8$, and this leads to increases with the wall distance. These behaviours 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 $y/\delta _{99} \le 0.7$, a reasonable agreement between the simulations and the experiment is observed, and results from the different wall models not showing a notable difference.

Figure 15.

Directions (relative to the local free stream) of the mean velocity gradient vector and the Reynolds shear stress vector: (a) station 10; (b) station 17. Red, EQWM; blue, PDE NEQWM; green, integral NEQWM; black, experiment. Circle, angle between the Reynolds shear stress vector $\gamma _{\tau }$ and the local free stream; cross, angle between the mean velocity gradient vector $\gamma _{g}$ and the local free stream.

3.5. 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 analyse the Reynolds-stress anisotropy. Following Pope (Reference Pope2000), the normalized anisotropy tensor is defined as

(3.4)\begin{equation} {\mathsf{b}}_{ij} = \frac{\langle u_{i}u_{j} \rangle}{\langle u_{k}u_{k} \rangle} - \frac{1}{3}\delta_{ij}. \end{equation}

The anisotropy tensor has zero trace and thus has two independent principal invariants. It is convenient to define two variables $\xi$ and $\eta$, corresponding to the two invariants, as

(3.5)$$\begin{gather} 6\eta^2=-2I_{2}={\mathsf{b}}^2_{ii}={\mathsf{b}}_{ij}{\mathsf{b}}_{ji}, \end{gather}$$

(3.6)$$\begin{gather}6\xi^3=3I_{3}={\mathsf{b}}^3_{ii}={\mathsf{b}}_{ij}{\mathsf{b}}_{jk}{\mathsf{b}}_{ki}. \end{gather}$$

The state of anisotropy can be characterized by the above two variables $\xi$ and $\eta$. All realizable Reynolds-stress states must be located within a triangular region in the $\xi$–$\eta$ plane, which is known as the Lumley triangle. The boundary of the Lumley triangle corresponds to some special states of the Reynolds-stress tensor: the origin corresponds to the isotropic turbulence; the left corner point corresponds to the two-component (2C) axisymmetric state; the right corner point corresponds to the one-component (1C) state; the left straight line corresponds to the axisymmetric contraction and the right straight line corresponds to the axisymmetric expansion; the top curve represents the two component (2C) turbulence.

The evolution of the Reynolds stress anisotropies through the bend is shown in figure 16. Comparisons to a statistically 2-D channel flow at $Re_\tau = 2003$ (Hoyas & Jimenez Reference Hoyas and Jimenez2006) and a transient statistically 3-D channel flow at $Re_\tau = 546$ (Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2020) can be also made from these data shown in figure 16(d). Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) studied a transient three-dimensional channel flow, where an initially statistically one-dimensional flow in a doubly periodic channel evolves to a new state after a sudden imposition of a spanwise pressure gradient. This flow models a shear-driven 3DTBL, and it can be used for comparison to the pressure-driven case discussed in the current work. Right upstream of the bend, the wall-normal distribution of the anisotropies away from the wall shows some similarity to that in the 2-D channel (figure 16a), exhibiting 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 towards 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 2-D channel. Although station 18 is considerably downstream of the bend region, the anisotropies are still seen to further depart from its 2-D behaviour. This is consistent with the observation in figure 15 that the Reynolds stresses respond more slowly than the mean to the imposed three-dimensionality. Also, in figure 16(d), note the similarity of the anisotropy distributions in the duct and the shear-driven 3DTBL from the transient channel flow, although departure from the 2-D behaviour is much stronger in the duct case.

Figure 16.

Lumley triangle of the WMLES results along the wall-normal direction: (a) station 4; (b) station 10; (c) station 12; (d) station 18. Coloured dots, fine grid PDE NEQWM result (all the wall models give almost identical results, and thus only PDE-NEQWM is shown here); black dash–dotted line, canonical 2-D channel flow at $Re_{\tau }=2003$ (Hoyas & Jimenez Reference Hoyas and Jimenez2006); red solid line, shear-driven 3DTBL from transient channel at $Re_{\tau }=546$ at $t^+=192$ (Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2020). Colour bar denotes the wall distance normalized by the local boundary layer thickness.

3.6. 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 cross-flow so developed is the ‘Johnston triangular plot’ (Johnston Reference Johnston1960), which is the triangular plot of $U_s$ against $U_n$, where $U_s$ and $U_n$ are along and normal to the local free stream 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 (Hawthorne Reference Hawthorne1951; Squire & Winter Reference Squire and Winter1951; Bradshaw Reference Bradshaw1987):

(3.7)\begin{equation} \frac{U_{n}}{U_{e}}=2\gamma_{e}\left(1-\frac{U_{s}}{U_{e}}\right), \end{equation}

where $\gamma _{e}$ is the free stream 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 towards the right end of the triangular plot. In figure 17, we present the mean velocity for the current bent duct flow and a temporally developing shear-driven 3-D channel flow from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2020) 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 free stream turning angle with respect to the upstream flow, and the free stream slope in the triangular plot (figure 17a) therefore increases towards the downstream direction.

Figure 17.

Johnston triangular plots. (a) WMLES of the bent square duct (all the wall models give almost identical results, and thus only PDE-NEQWM is shown here): square, station 0; circle, station 4; triangle, station 6; cross, station 8; diamond, station 10; star, station 12. (b) DNS of the shear-driven 3DTBL from the transient channel flow at $Re_{\tau }=546$ at $t^+=192$ (Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2020). Red straight line, the SWH formula (3.7). Colour bar denotes the wall distance normalized by the local boundary layer thickness.

The wall model solutions (defined only close to the wall) are also visualized along with the outer LES solution in figure 18. It gives us a vivid picture of the different wall models’ capabilities to depict skewed mean-velocity profiles. In the triangular plot, the wall model solutions are from the origin to the third 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. However, 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 cross-flow developing stage (figure 18a), 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 cross-flow decaying stage (figure 18b), the difference among three wall model solutions gradually decreases. All three wall model solutions become almost unidirectional close to the end of the duct.

Figure 18.

Johnston triangular plot of LES and wall-model solutions. (a) Cross-flow developing stage: square, station 0; circle, station 4; triangle, station 6; cross, station 8; diamond, station 10; star, station 12. (b) Cross-flow decaying stage: square, station 12; circle, station 14; triangle, station 16; cross, station 18; diamond, station 20; star, station 21. Red solid straight line, relation given by the SWH formula (3.7); red dash–dotted line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM. Lines with symbols represent the WMLES solution similar to figure 17(a). Colour bar denotes the wall distance normalized by the local boundary layer thickness.

3.7. Quantification of the non-equilibrium contributions to the wall shear stress direction

In this section, the non-equilibrium effects neglected in the EQWM are analysed through the full RANS equations used in the non-equilibrium wall models. This analysis highlights the importance of accounting for the non-equilibrium effects in wall modelling 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 using the fully resolved flow fields, as attempted by Hickel et al. (Reference Hickel, Touber, Bodart and Larsson2012). 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 recast as

(3.8)$$\begin{gather} \frac{\partial}{\partial y}\left[(\nu+\nu_t)\frac{\partial \langle u\rangle }{\partial y}\right] = S_x, \end{gather}$$

(3.9)$$\begin{gather}\frac{\partial}{\partial y}\left[(\nu+\nu_t)\frac{\partial \langle w\rangle }{\partial y}\right] = S_z, \end{gather}$$

where $S_x$ and $S_z$ are the non-equilibrium source terms comprising the following terms ($S = A + P - D$):

(3.10)$$\begin{gather} \textrm{advection } A_x = \frac{\partial \langle u\rangle^2}{\partial x} + \frac{\partial \langle u\rangle \langle v\rangle }{\partial y} + \frac{\partial \langle u\rangle \langle w\rangle }{\partial z} +\frac{\partial \langle u'u'\rangle}{\partial x} + \frac{\partial \langle u'v'\rangle}{\partial y} + \frac{\partial \langle u'w'\rangle}{\partial z}, \end{gather}$$

(3.11)$$\begin{gather}\textrm{advection } A_z = \frac{\partial \langle w\rangle^2}{\partial z} + \frac{\partial \langle v\rangle \langle w\rangle }{\partial y} + \frac{\partial \langle u\rangle \langle w\rangle }{\partial x} +\frac{\partial \langle w'w'\rangle}{\partial z} + \frac{\partial \langle v'w'\rangle}{\partial y} + \frac{\partial \langle u'w'\rangle}{\partial x}, \end{gather}$$

(3.12)$$\begin{gather}\textrm{pressure gradient } P_x = \frac{1}{\rho}\frac{\partial \langle P\rangle}{\partial x}, \end{gather}$$

(3.13)$$\begin{gather}\textrm{pressure gradient } P_z = \frac{1}{\rho}\frac{\partial \langle P\rangle}{\partial z}, \end{gather}$$

(3.14)$$\begin{gather}\textrm{lateral diffusion } D_x = \frac{\partial}{\partial x}\left[(\nu+\nu_t)\frac{\partial \langle u\rangle }{\partial x}\right] + \frac{\partial}{\partial z}\left[(\nu+\nu_t)\frac{\partial \langle u\rangle }{\partial z}\right], \end{gather}$$

(3.15)$$\begin{gather}\textrm{lateral diffusion } D_z = \frac{\partial}{\partial x}\left[(\nu+\nu_t)\frac{\partial \langle w\rangle }{\partial x}\right] + \frac{\partial}{\partial z}\left[(\nu+\nu_t)\frac{\partial \langle w\rangle }{\partial z}\right]. \end{gather}$$

By integrating (3.8) and (3.9) twice, the following expressions of the wall shear stress for the PDE NEQWM are obtained:

(3.16)$$\begin{gather} \tau_{w,x} = \frac{\displaystyle U_{LES}-\int_{0}^{h_{wm}}\frac{\displaystyle \int_{0}^{y}S_x\, {\rm d} y}{\nu+\nu_t}{{\rm d} y}}{\displaystyle \int_{0}^{h_{wm}}\frac{1}{\nu+\nu_t}{{\rm d} y}}, \end{gather}$$

(3.17)$$\begin{gather}\tau_{w,z} = \frac{\displaystyle W_{LES}-\int_{0}^{h_{wm}}\frac{\displaystyle \int_{0}^{y}S_z\, {\rm d} y}{\nu+\nu_t}{{\rm d} y}}{\displaystyle \int_{0}^{h_{wm}}\frac{1}{\nu+\nu_t}{{\rm d} y}}, \end{gather}$$

where $U_{LES}$ and $W_{LES}$ are the LES velocity components at the matching location. A similar expression has been reported by Wang & Moin (Reference Wang and Moin2002). As shown in the previous sections, $U_{LES}$ and $W_{LES}$ 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

(3.18)\begin{equation} \frac{\tau_{w,z}}{\tau_{w,x}} = \frac{W_{LES}-I_z}{U_{LES}-I_x}, \end{equation}

where $I_z = \int _{0}^{h_{wm}}({\int _{0}^{y}S_z\, {\rm d} y}/({\nu +\nu _t}))\,{{\rm d} y}$ and $I_x = \int _{0}^{h_{wm}}({\int _{0}^{y}S_x\, {\rm d} y}/({\nu +\nu _t}))\,{{\rm d} y}$. When all the non-equilibrium effects are neglected (i.e. letting $S_x=S_z=0$), this relation reduces to the wall shear stress direction predicted by the EQWM which assumes unidirectional flow (${\tau _{w,z}}/{\tau _{w,x}} = {W_{LES}}/{U_{LES}}$). The fidelity with which the constitutive terms of $I_x$ and $I_z$ are modelled is, therefore, crucial to the performance of wall models in predicting the surface flow direction.

To separate the non-equilibrium contributions from the flow direction predicted by the EQWM, we can first reorganize (3.18) as

(3.19)\begin{equation} \frac{\tau_{w,z}}{\tau_{w,x}} = \frac{W_{LES}}{U_{LES}}\left(\frac{1-I_z/W_{LES}}{1-I_x/U_{LES}}\right). \end{equation}

For the present flow, it is shown in figure 19 that $I_x$ is relatively small compared to $U_{LES}$, which permits the use of a truncated Taylor series expansion of ${1}/({1-I_x/U_{LES}})$,

(3.20)\begin{equation} \frac{\tau_{w,z}}{\tau_{w,x}} = \frac{W_{LES}}{U_{LES}}\left(1-\frac{I_z}{W_{LES}}\right) \left[1+\frac{I_x}{U_{LES}}+O\left(\frac{I_x}{U_{LES}}\right)^2\right]. \end{equation}

Figure 19.

Centreline distribution of the non-equilibrium contributions (from PDE NEQWM). Red solid line, $U_{LES}/U_{ref}$; red dashed line, $W_{LES}/U_{ref}$; black solid line, $-I_x/U_{ref}$; black dashed line, $-I_z/U_{ref}$.

Equation (3.20) can be further expanded as

(3.21)\begin{equation} \frac{\tau_{w,z}}{\tau_{w,x}} = \frac{W_{LES}}{U_{LES}}\left(1-\frac{I_z}{W_{LES}}+ \frac{I_x}{U_{LES}}+\cdots\right). \end{equation}

The terms $I_x/U_{LES}$ and $I_z/W_{LES}$ can be viewed as representing the corrective contributions from the non-equilibrium effects to the surface flow direction predicted by the EQWM. In essence, (3.21) enables the surface flow angle to be decomposed into distinct contributions originating from the equilibrium and the non-equilibrium terms. Although not shown here for brevity, this truncated relation was found to provide almost identical description as the actual wall flow direction ($\tau _{w,z}/\tau _{w,x}$). Note that (3.21) applies to the integral NEQWM as well. However, the terms $I_x$ and $I_z$ therein are largely modelled 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 (3.21) for the two non-equilibrium wall models in figure 20 (all non-equilibrium terms are zero for EQWM and they are not plotted.) Several interesting observations are made. First, the two non-equilibrium wall models show that the total non-equilibrium 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 non-equilibrium models properly sense the region where non-equilibrium effects are important and attempt to model them therein. Second, the two non-equilibrium wall models produce comparable distributions of $I_x/U_{LES}$, implying that the axial ($x$) contents of the non-equilibrium effects are modelled almost identically by the two wall models. Third, the difference between the wall models appears to be concentrated in $-I_z/W_{LES}$, i.e. in the way the models sense and model the cross-flow ($z$) component of the non-equilibrium effects. The integral NEQWM underpredicts $-I_z/W_{LES}$ 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 comprises advection ($A_z$), pressure gradient ($P_z$), lateral diffusion ($D_z$) and $W_{LES}$. Among these terms, $P_z$ and $W_{LES}$ are imposed largely from the LES solution (which are seen to be identical from the simulations using the two wall models), and $D$ is seen to be negligible in its magnitude. Therefore, we conjecture that the difference of $-I_z/W_{LES}$ in the integral wall model originates largely from its assumed velocity profile, which is directly used in computing the advection term $A_{z}$. It should be noted that the lines representing $-I_z/W_{LES}$ in figure 20 show rapid changes in the upstream region ($0.2\leq x'/L \leq 0.4$) because the spanwise velocity $W_{LES}$ is almost zero therein. This makes these terms ill-behaved (division by 0) in this region with an extreme data range. It should also be noted that these terms will be multiplied by $W_{LES}/U_{LES}$ in the end, as in (3.21). When $W_{LES}$ is almost zero, the flow angle will just be zero. Therefore, the focus should be on the bend and the downstream regions.

Figure 20.

Centreline distribution of the non-equilibrium contribution to the flow direction. Solid line, $I_x/U_{LES}$; dashed line, $-I_z/W_{LES}$; dotted line, total non-equilibrium correction. Blue, NEQWM; green, integral NEQWM. Black vertical dashed lines, start and end of the bend region.

Furthermore, the individual contributions from different non-equilibrium effects are also analysed. Starting from (3.18) with $I_x = I_z =0$ (corresponding to the EQWM prediction), we examine the change in the surface flow turning angle by systematically including the contributions from various non-equilibrium effects (advection, pressure gradient and lateral diffusion) to $I_x$ and $I_z$. Note that this is done in a post-processing manner using the solution of the PDE NEQWM. The results are shown in figure 21. When all the non-equilibrium 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 free stream, while the advection tends to make the flow deviate less from the local free stream. 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.

Figure 21.

Centreline distribution of the surface flow turning angles with respect to the free stream ($\gamma _w$ is the wall shear stress direction, $\gamma _{\infty }$ is the free stream direction): (a) including all non-equilibrium effects; (b) including diffusion only; (c) including pressure gradient only; and (d) including advection only. Red dash–dotted line, EQWM; blue dashed line, PDE NEQWM; blue solid line with circles, reconstruction with non-equilibrium contributions; black squares, experiment; black vertical dashed line, start and end of the bend region.