2.1. Baseline method and case set-up

A nominally canonical zero pressure gradient (ZPG) TBL flow is simulated by using the open-source CFD solver OpenFOAM. The unsteady incompressible flow equations in a flow domain $\varOmega = [0,\; {L_x}] \times [0,\; {L_y}] \times [0,\; {L_z}]$ are

(2.1a)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,\end{gather}

(2.1b)\begin{gather}\partial \boldsymbol{u}/\partial t + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} - \nu {\nabla ^2}\boldsymbol{u} + (1/\rho )\boldsymbol{\nabla }p = 0,\end{gather}

where $\boldsymbol{u}$ is the velocity vector, $\rho $ is the constant density in the incompressible flow and $\nu $ is the kinematic viscosity. The pressure gradient $\boldsymbol{\nabla }p$ should be nominally zero both in the laminar and turbulent regimes. The domain lengths ${L_x}$, ${L_z}$ and ${L_y}$ are those in streamwise x, spanwise z and wall-normal y directions, respectively.

The fluid domain is illustrated in figure 1. The corresponding boundary conditions are as follows. The inlet velocity profile is set as constant with freestream velocity ${U_\infty }$. The pressure gradient is set as zero at the outlet. The periodic boundary conditions are applied to the spanwise directions. The upper boundary is a slip wall mimicking the ideally infinite wall-normal height. The lower boundary is the no-slip wall. The domain length ${L_x}$ is designed to reach the targeted nominal Reynolds number $R{e_\theta } = 1000$, where $R{e_\theta }$ is the Reynolds number based on the freestream velocity ${U_\infty }$, the kinematic viscosity $\nu $ and the momentum thickness $\theta $. The domain height is taken to contain approximately seven boundary layer thicknesses $\delta $ at the outlet. The spanwise domain size ${L_z}\sim 0.5{L_y}$ should be large enough to contain the largest scale spanwise structures. A preliminary sensitivity study with doubled ${L_z}$ has been carried out and shows negligible differences. The overall domain covers $R{e_x}$ from 0 to $7 \times {10^5}$, which is based on the freestream velocity ${U_\infty }$, the kinematic viscosity $\nu $ and the streamwise location x.

Figure 1. Schematic of the computational domain (the black cone indicating the trip location, with a close-up view of local mesh).

The flow equations are discretised in a finite volume scheme, with a second-order central difference scheme in space. The time advancement is achieved by a second-order Crank–Nicolson scheme. A constant time-step is taken in keeping the maximum Courant numbers below 0.5. The pressure-implicit splitting operators (PISO) algorithm is used for solving the incompressible flow field. In the present calculations, no explicit sub-grid model is involved. All scales are directly resolved without any turbulence models. The simulations are therefore effectively implicit LESs where the sub-grid scales are dampened only by numerical dissipations. The simulations are run for a sufficiently long time to reach a statistically steady state when the mean statistics differences are confined within 1 %. All the statistics to be reported are time-averaged for at least ten eddy turnover times $\delta /{u_\tau }$ (Wu & Moin Reference Wu and Moin2009; Lee & Moser Reference Lee and Moser2015), where $\delta $ here is the boundary layer thickness at the exit and ${u_\tau }$ is the wall shear velocity there. The statistics are spatially averaged along the homogeneous spanwise direction.

The tripping element covers the entire span $({L_z})$. The element height is set to be ${\sim} 0.78{\delta ^\ast }$, as suggested in the experimental study by Dryden (Reference Dryden1959), where ${\delta ^\ast }$ is the displacement thickness at the tripping location based on the Blasius laminar boundary layer solution. The streamwise length of the tripping element is the size of the outer mesh spacing and roughly equals the element height. The distance from the inlet (the leading edge) to the trip is approximately 300 streamwise trip element lengths. The distance from the trip to the outlet depends on the Reynolds number, and for the nominal case is approximately 500 trip element lengths.

The inner region (coloured in blue in figure 1) and outer region (coloured in green in figure 1) meshes are connected by a wall-parallel interface at a wall-normal distance ${y_s}$, which is nominally determined by $y_s^ +\approx 3{({\delta ^ + })^{0.5}}$, following Marusic et al. (Reference Marusic, Monty, Hultmark and Smits2013), as an estimation on the lower bound of the log-law layer. The superscript ‘$+ $’ denotes values normalised by shear velocity ${u_\tau }$ and viscosity $\nu $. The inner–outer mesh interface is located roughly at the beginning of the log region in the y direction, which is around $y_s^ +\approx 70$ at $R{e_\theta } \approx 1000$. Given that the inner and outer meshes have different node distributions, there should be hanging nodes at the mesh interface. The interface is treated as the non-conforming arbitrary mesh interface (AMI) patches based on the Galerkin projection (Farrell & Maddison Reference Farrell and Maddison2011) for conservation and instantaneous flow field interpolation.

The mesh resolutions of different cases are shown in table 1. Note that the wall-normal mesh resolution is kept the same in all cases, gradually stretching from $\Delta y_w^ +{\sim} 0.5$ at the lower solid wall to the upper boundary. There are approximately 10 mesh points below ${y^ + } = 5$, 55 points below $y_s^ + $ and 80 points located within the boundary layer at $R{e_\theta } \approx 1000$. Following Moin (Reference Moin1997) and Jimenez & Moser (Reference Jimenez and Moser2000), the near-wall mesh of WRLES should be as fine as that of DNS. A uniformly distributed mesh is used in both spanwise and streamwise directions as the baseline mesh.

Table 1. Details of mesh resolution of the cases simulated.

Three cases with different inner mesh resolutions, denoted as ‘Base-x40z30’ (coarse LES), ‘Base-x20z15’ (fine LES) and ‘Base-x10z7’ (DNS), are set up with quadrupole refinement in both x and z directions. Note that resolutions in wall units (inner and outer mesh resolutions in table 1) are calculated based on wall shear velocity ${u_\tau }$ at $R{e_\theta } \approx 1000$. Thereafter, the investigation on the sensitivity of local mesh refinement around the trip is carried out for the three cases with a locally refined mesh around the trip: ‘RFtrip-x40z30’, ‘RFtrip-x20z15’ and ‘RFtrip-x10z7’. The trip is refined as indicated in figure 1 with the mesh grids clustered around the trip for RFtrip-x20z15 as an example. The resolutions at the trip itself are calculated using local time-averaged ${u_\tau }$ at each corresponding trip surface based on the first mesh cell size in the wall normal directions ${\boldsymbol{n}_x}$ and ${\boldsymbol{n}_y}$ (figure 1).

2.2. Full WRLES validation and mesh-resolution requirement for tripping

For the cases considered, the streamwise pressure gradient should be nominally zero in a zero-pressure-gradient (ZPG) TBL flow, but this is practically almost impossible to achieve either experimentally or numerically as discussed by Wu & Moin (Reference Wu and Moin2009). The indicators of the pressure gradient magnitude of Wu & Moin (Reference Wu and Moin2009) are used here to evaluate the credibility of the simulated ZPG TBL: $\beta = ({\delta ^\ast }/{\tau _w})(\mathrm{\partial }\bar{p}/\mathrm{\partial }x)$; $- 20{P^ + } ={-} 20(\nu /\rho u_\tau ^3)(\mathrm{\partial }\bar{p}/\mathrm{\partial }x)$, where $\bar{p}$ is the time-averaged pressure, ${\tau _w}$ denotes the wall shear stress and $\rho $ is the constant fluid density in the incompressible flow. The averaged values in the fully turbulent regime are compared in table 2. The results are of the same order of magnitude as those of Wu & Moin (Reference Wu and Moin2009) and an order of magnitude lower than those of Spalart & Watmuff (Reference Spalart and Watmuff1993) for assurance.

Table 2. Averaged streamwise pressure gradient at approximately $y\sim {L_y}/8$.

The calculated time-averaged friction coefficients ${C_f}$ as a function of the momentum thickness based Reynolds number $R{e_\theta }$ are shown in figure 2. All cases overlap in the laminar regime. In the turbulent regime, as the inner mesh is quadropoly refined in the streamwise and spanwise directions, a clear mesh-convergence can be observed. It can also be observed that the transition begins at $R{e_\theta } \approx 250\unicode{x2013} 300$ and the flow becomes turbulent at $R{e_\theta } \approx 450$ after the trip. The local refinement at the trip slightly improves the results downstream close to the trip, shown by the comparison between Base-x20z15 (figure 2a) and RFTrip-x20z15 (figure 2b). The results match well with the fitted curve based on the experimental data of Smits et al. (Reference Smits, Matheson and Joubert1983). Note that the results in the transition regions (covering the tripping and downstream recirculation regions) are not shown here as a calculated $R{e_\theta }$ for this region with reverse flows would be hardly meaningful.

Figure 2. Skin friction coefficients as functions of $R{e_\theta }$ in the streamwise direction compared with Blasius laminar solution and the correlated curve for the turbulent part (solid grey lines) by Smits et al. (Reference Smits, Matheson and Joubert1983). (a) For the base meshes: Base-x10z7, open brown circles; Base-x20z15, open blue squares; Base-x40z30, open red triangles. (b) For the trip-refined meshes: RFTrip-x10z7, solid brown circles; RFtrip-x20z15, solid blue squares; RFtrip-x40z30, solid red triangles. Note that not all data are shown for clarity.

Based on the results in figure 2, we can make some more general observations regarding the mesh resolution required for the trip. Note first that the DNS base inner mesh (Base-x10z7) by itself without any local refinement seems to be capable of resolving the roughness element well. Perhaps more relevantly, the corresponding first wall-normal mesh spacings for this mesh $(\Delta n_x^ += 2,\;\Delta n_z^ += 1.4)$ turn out to be all within the local viscous sublayer (see table 1). However, it needs to be stressed that the local wall normal mesh spacing within the sublayer should not be taken as a sufficient condition. This point is highlighted by the two cases of the coarse LES meshes with the local refinement (RFtrip-x20z15, RFtrip-x40z30). Even when the local refinement leads to an adequate first wall-normal mesh spacing at the trip ($\Delta n_x^ += 1.7,\;\Delta n_z^ += 1.3$, table 1), the overall resolution for the tripping may still be insufficient, as shown in figure 2(b). It follows therefore that the mesh resolution required for resolving a trip properly should satisfy a dual requirement: (a) a fine enough first wall-normal mesh spacing well within the local viscous sub-layer (e.g. $\Delta {n^ + } < 2$) on the trip element surfaces; and (b) a fine enough mesh for the inner boundary layer region around the trip (e.g. a DNS mesh for the near-wall region).

The present finding as discussed above is consistent with the established consensus that for a near-wall flow, the mesh resolution of LES should approach that of DNS (Moin Reference Moin1997; Jimenez & Moser Reference Jimenez and Moser2000). Indeed, only would the two DNS inner meshes with and without the local trip refinement (Base-x10z7 and RFtrip-x10z7) satisfy the dual requirement for being sufficiently fine both on the wall surfaces of the trip and in the inner boundary layer field around the trip. It is thus not surprising that for the DNS base inner mesh, the extra local refinement around the trip (RFtrip-x10z7) should not be necessarily needed.

The finer LES inner mesh (Base-x20z15) seems to be a borderline case close to the DNS resolution for the overall inner flow region. In this case, the local refinement around the trip (RFtrip-x20z15’) with only approximately 2 % extra mesh count does lead to a more identifiable improvement for the tripping solution, as seen by comparing the corresponding results between figures 2(a) and 2(b). The present results suggest that if the inner region mesh is sufficiently fine for a normal TBL, resolving a physical trip element should not require a significant extra mesh resolution locally.

More sensitive indicators, the ratio between the boundary layer thickness and the displacement or momentum thickness, $\delta /{\delta ^\ast }$ and $\delta /\theta $, as used by Schlatter & Orlu (Reference Schlatter and Orlu2010), are also presented for a closer examination of the local mesh refinement around the trip. For the two cases Base-x20z15 and RFtrip-x20z15, the indicators converge at approximately $R{e_\theta } \approx 700$, as shown in figure 3. There is little difference at a higher Reynolds number further downstream and both cases asymptotically converge to the correlated experimental data of Chauhan et al. (Reference Chauhan, Monkewitz and Nagib2009).

Figure 3. Indicative parameters $\delta /{\delta ^\ast }$ and $\delta /\theta $ compared with the correlated curve of Chauhan et al. (Reference Chauhan, Monkewitz and Nagib2009) (solid grey lines). The present results are shown for RFTrip-x20z7 (solid black lines) and Base-x20z7 (dashed red lines).

Figure 4 shows the shape factor, ${H_{12}} = {\delta ^\ast }/\theta $, as the function of $R{e_\theta }$ for the boundary layer development. The results from the case RFtrip-x10z7 match up with the fitted curve of TBL development by Monkewitz et al. (Reference Monkewitz2007) in the turbulent regime with the shape factor ${H_{12}}$ asymptotically approaching ${H_{12}} \approx 1.4$. Figure 5 shows that the equivalent shear Reynolds number $({\delta ^ + })$ based on local boundary layer thickness and shear velocity for case RFTrip-x10z7, which matches well in the turbulent regime with the fitted curve on TBL development (Schlatter & Orlu Reference Schlatter and Orlu2010).

Figure 4. Shape factor ${H_{12}}$ as functions of $R{e_\theta }$ compared with fitted curve as a solid black line from Monkewitz et al. (Reference Monkewitz2007). RFTrip-x10z7, solid brown circles. Note that not all data are shown for clarity.

Figure 5. Boundary layer thickness ${\delta ^ + }$ as functions of $R{e_\theta }$ compared with the fitted curve as a solid black line of Schlatter & Orlu (Reference Schlatter and Orlu2010). RFTrip-x10z7, solid brown circles. Note that not all data are shown for clarity.

Figure 6 presents the time-mean velocity profiles, the root mean squared (rms) velocity fluctuations and the Reynolds stresses (ensemble averaged values are denoted as ‘$\langle \cdot \rangle $’) as a function of wall-normal distances ${y^ + }$. Results at two Reynolds numbers are compared with the published DNS results. The present results of RFtrip-x10z7 match with the DNS data well at $R{e_\theta } \approx 1000$, better than that at $R{e_\theta } \approx 670$ where the flow presumably may not yet be fully developed. Interestingly, the difference between the two DNS datasets at $R{e_\theta } \approx 670$ is also non-trivial, probably due to their different ways of initiating TBL. In addition, we see that outer flows are quite well resolved by the coarser mesh (table 1) with the Reynolds stresses matching with the DNS data at $R{e_\theta } \approx 1000$.

Figure 6. Mean statistics with respect to wall-normal distances ${y^ + }$. (a,b) Mean velocity profiles. (c,d) Streamwise velocity fluctuations. (e, f) Reynolds shear stresses. Results are at (a,c,e) $R{e_\theta } \approx 670$ and (b,d, f) $R{e_\theta } \approx 960$. The DNS data are shown at the corresponding $R{e_\theta }$: Schlatter & Orlu (Reference Schlatter and Orlu2012) as solid black lines; Wu & Moin (Reference Wu and Moin2009) as dashed green lines. The present results are shown as solid brown circles. The grey lines mark the inner–outer mesh interfaces. Note that not all data are shown for clarity.

2.3. Development of tripped boundary layer towards an equilibrium state

An instantaneous view of vortical flow structures from the trip to the developed turbulent boundary layer state downstream is illustrated in figure 7. The structures are visualised in terms of the Q-criteria (Hunt, Wary & Moin Reference Hunt, Wary and Moin1988; Jeong & Hussain Reference Jeong and Hussain1995) for case RFtrip-x10z7. The laminar boundary layer is buffeted by the 2-D step with vortices shed downstream from the top surface, similar to a roughness-induced bypass transition, e.g. Rao et al. (Reference Rao, Jefferson-Loveday, Tucker and Lardeau2014). The roll-up of vortices becomes spanwise non-uniform and breaks down very soon downstream, similar to that observed by Brinkerhoff & Yaras (Reference Brinkerhoff and Yaras2011). The instability develops downstream as the spanwise vortices break down into three-dimensional turbulent structures. From $R{e_\theta } \approx 350$ to $R{e_\theta } \approx 420$, when the new boundary layer is formed after the reattachment, the hairpin vortices (Adrian Reference Adrian2007) start to emerge. Note that a natural transition triggered by Tollmien–Schlichting (TS) waves follows the classical route from initial linear primary disturbances to secondary instabilities (Sayadi, Hamman & Moin Reference Sayadi, Hamman and Moin2013), where the formation of linear $\varLambda $-shape vortices should be spatially long. These features are obviously bypassed in the present situation.

Figure 7. Iso-surfaces of the Q-criterion in the inner region with refined mesh. The coloured contour levels indicate the normalised velocity $U/{U_\infty }$.

Figure 8 shows the results of the present calculation (RFtrip-x10z7), the TS wave-associated natural transition (TS), and the artificial-trip-induced transition (OP-KTH) by Schlatter & Orlu (Reference Schlatter and Orlu2012) in terms of the indicative parameters $\delta /{\delta ^\ast }$ and $\delta /\theta $. An obviously shorter development length of the present calculation can be observed compared with the natural transition (TS). The synthetic volume force tripping with ‘optimally’ tuned parameters (OP-KTH) enforces a fully developed turbulent regime resulting in a nearly straight line shortly after the tripping, distinctively different from the other results.

Figure 8. Indicators $\delta /{\delta ^\ast }$ and $\delta /\theta $ compared with the experimentally correlated curve of Chauhan et al. (Reference Chauhan, Monkewitz and Nagib2009) (solid grey lines). The results of the tuned ‘optimal’ tripping of Schlatter & Orlu (Reference Schlatter and Orlu2012) (‘OP-KTH’) are shown as dash-dotted black lines; the natural TS-wave transitions (‘TS’) are shown as dashed blue lines; the present calculations (RFTrip-x10z7) are shown as solid brown circles.

It may be more informative that the results shown in figure 8 are categorised into two groups in terms of how each solution approaches asymptotically to the eventual almost constant values (indicative of an equilibrium state), approximately 8.5 for $\delta /\theta $ and 5.8 for $\delta /{\delta ^\ast }$. The first group consists of those physical transition routes: the natural transition (TS), the tripping experimentally implemented (Chauhan et al. Reference Chauhan, Monkewitz and Nagib2009) or numerically simulated (RFtrip-x10z7). The results of this group first overshoot and then decrease gradually to approach the equilibrium state. The initial overshoot indicates a relatively thicker boundary layer due to the outer flow region being disturbed by large-scale disturbances, which seemingly have a relatively small influence on the inner part with a large velocity deficit more responsible for contributing to the two integral parameters, ${\delta ^\ast }$ and $\theta $. However, in the second group, the ‘optimally’ tripped one (OP-KTH), the two parameters $\delta /{\delta ^\ast }$ and $\delta /\theta $ approach the final state rather differently, both increasing with $R{e_\theta }$ monotonically. The similarity among those physically tripped turbulent flows in approaching the final equilibrium state draws further attention to the characteristics of a post-transition non-equilibrium TBL.

Now, we first use the present results to elaborate on the overall tripping process from the viewpoint of turbulent energy production and dissipation. Then we will examine the characteristics and influence of those tripping-induced large-scale disturbances through modal analyses.

We introduce the dissipation coefficient ${C_d}$ for time-averaged incompressible flow (e.g. Wheeler et al. Reference Wheeler, Dickens and Miller2018) as

(2.2)\begin{equation}{C_d} = \epsilon + Pr = \frac{1}{{U_\infty ^3}}\int_0^\delta {\nu {{\left( {\frac{{\partial {U_i}}}{{\partial {x_i}}}} \right)}^2}\,\textrm{d}y} + \left( { - \frac{1}{{U_\infty^3}}\int_0^\delta {\overline {{u_i}{u_j}} \frac{{\partial {U_i}}}{{\partial {x_j}}}\,\textrm{d}y} } \right),\end{equation}

where the dissipation coefficient ${C_d}$ represents the loss of mean kinetic energy as a combination of the viscous dissipation of the mean flow $(\epsilon )$ and the production of turbulence kinetic energy $(Pr)$ in the spatially evolving boundary layer flow.

The full wall-resolved LES solutions (‘RFTrip-x10z7’ introduced in § 2.2) are processed based on (2.2), as shown in figure 9. The entire route of the boundary layer development may be categorised into four stages: laminar, transitional, non-equilibrium and equilibrium. Wheeler et al. (Reference Wheeler, Dickens and Miller2018) showed that ${C_d}$ is initially very high within the transitional stage, plunges to a normal level for turbulent boundary layers and eventually approaches asymptotically to a value of 0.002. A similar process can be observed in figure 9. Initially, in the laminar regime, there is zero turbulence kinetic energy production. Then it takes a journey into the non-equilibrium turbulent phase after the transition, when the boundary layer becomes fully turbulent yet not fully developed. The non-equilibrium state is attributed to the fact that the turbulence kinetic energy generation rate reacts much faster in transition than the viscous dissipation rate of the mean flow, resulting in an imbalance (Wheeler et al. Reference Wheeler, Dickens and Miller2018). Note that the peak of the viscous dissipation rate $\epsilon $ lags behind the peak of the total dissipation coefficient ${C_d}$ (as shown in the later part of the transitional stage, figure 9). The turbulent boundary layer gradually develops towards an equilibrium state after $R{e_x}/{10^5} \approx 6$ with ${C_d} \approx 0.002$ ($\epsilon \approx 0.001$ and $Pr \approx 0.001$).

Figure 9. Dissipation coefficient ${C_d}$ as functions of $R{e_x}$, categorised in four stages: laminar, transitional, non-equilibrium and equilibrium. The blue solid line depicts the viscous dissipation of the mean flow $\epsilon $. The orange solid line depicts the dissipation coefficient ${C_d} = \epsilon + Pr$. The grey dashed line shows the equilibrium ${C_d} \approx 0.002$.

To examine the post-trip non-equilibrium turbulent region (from $R{e_x}/{10^5} \approx 3.4\unicode{x2013} 3.8$) more closely, we apply the fast Fourier transfer to obtain the one-dimensional pre-multiplied energy spectra for this region, as shown in figure 10(a). Notice that the wall shear Reynolds number $R{e_\tau }$ in the middle of the region is approximately 190. The present results are thus compared with the DNS data at $R{e_\tau } \approx 192$ for the canonical channel flow (Lee & Moser Reference Lee and Moser2015), as shown in figure 10(b). An energy overshoot of the present results (figure 10a) can be clearly seen around the trip height, $y = {\varDelta _{y,trip}}({y^ + } \approx 40)$. This location is in the log-region, so we label it as being of an outer flow region in contrast to the near-wall region below the buffer layer. Large-scale disturbances in the outer flow region of a tripped turbulent boundary layer are similarly observed by Marusic et al. (Reference Marusic, Chauhan, Kulandaivelu and Hutchins2015).

Figure 10. One-dimensional pre-multiplied energy spectra generated from (a) present calculation in the non-equilibrium TBL and (b) DNS channel flow data at approximately the same $R{e_\tau } \approx 190$ (Lee & Moser Reference Lee and Moser2015). The contour levels are from 1.5 to 4 with an interval of 0.5. The white dashed lines mark the two wall-normal locations: ${y^ + } \approx 40$ (roughly $y \approx {\varDelta _{y,trip}}$ locally) and ${y^ + } \approx 13.5$. The solid red box highlights the excessive energy due to the tripping-induced large scales in the outer flow region.

Given the presence of large-scale disturbances in the outer flow region, we now would like to know how persistent downstream they are and whether they have a significant influence on the near-wall region. We will first look at the power spectra density (PSD) at different wall-normal locations in different streamwise stations. The PSD of the flow velocity component is calculated using Welch's method with the Hann window (Welch Reference Welch1967) to reduce finite sampling effects. The dimensionless frequencies are in the form of a Strouhal number $St = fL/U$, where the characteristic length and velocity are the trip height ${\varDelta _{y,trip}}$ and the freestream velocity, respectively.

Figure 11 shows the PSD at different streamwise and wall-normal locations. For the outer flow location taken at the trip height, vortex-shedding features can be observed right after the trip $(R{e_x}/{10^5} \approx 3.1)$, as shown in figure 11(a). The differences in the streamwise development between the outer region $(y = {\varDelta _{y,trip}},{y^ + } \approx 40)$ and the inner region $({y^ + } \approx 13.5)$ are more clearly illustrated by comparing figure 11(b) with figure 11(c). For the outer flow location (figure 11b), the spectrum is characterised by a rather flat low-frequency range, indicating a relatively higher magnitude of low-frequency large-scale disturbances. The residual large-scale disturbances from the non-smoothness of the spectrum around $St \approx 1$, are clearly detectable at $R{e_x}/{10^5} \approx 3.6$, and remain very persistent even at $R{e_x}/{10^5} \approx 7$ (figure 11b). In a clear contrast, at the inner location ${y^ + } \approx 13.5$ (figure 11c), the boundary layer develops into a typical turbulent flow shortly after the trip $(R{e_x}/{10^5} \approx 3.6)$, and it becomes just marginally more developed and smoother further downstream at $R{e_x}/{10^5} \approx 7$ (figure 11c). The excessive energy generated at the trip height seemingly affects the outer flow region more significantly than the inner near-wall region. This is in line with the prevalent observation that the outer region requires a longer evolution distance than the inner one to reach a full equilibrium state (Devenport & Lowe Reference Devenport and Lowe2022).

Figure 11. Power spectral density (PSD) with regard to the non-dimensional frequencies. (a) $y = {\varDelta _{y,trip}}$ right after the trip as a blue line; spectra in panels (b) and (c) are for two wall-normal distances and at two streamwise stations: (b) outer flow ($y = {\varDelta _{y,trip}}$, ${y^ + } \approx 40$); and (c) inner flow ${y^ + } \approx 13.5$. The results right after the trip $(R{e_x}/{10^5} \approx 3.6)$ are shown as red lines; results further downstream $(R{e_x}/{10^5} \approx 7)$ are shown as green lines. All plots are taken in the mid-span of the domain.

2.4. Evolution and impact of large-scale disturbances of tripping

To further examine the streamwise evolution of the trip-induced large-scale disturbances and their potential impact on the near-wall turbulence, we carry out modal analyses. First, the empirical mode decomposition (EMD) method (Huang et al. Reference Huang1998) is implemented to identify if large-scale structures would correlate between the outer and inner regions. The EMD algorithm splits the original signals into a set of intrinsic mode functions (IMFs) based on local characteristic scales without introducing any cut-off wavelengths. The EMD method of a 2-D version is applied to decompose small- and large-scale components in the present work. The last three of eight IMFs are retained as the long-pass filtered large-scale flow field. Essentially, the same procedure of the EMD method as in Chen & He (Reference Chen and He2022) is applied to the snapshot plane cuts of the instantaneous flow field, as shown in figures 12–14.

Figure 12. (a,c) Instantaneous velocities on cut planes at two wall-normal locations: $y_{out}^ +\approx 52$ and $y_{in}^ +\approx 13.5$, respectively. (b,d) Corresponding large-scale structures retained from the regions marked by black dashed boxes in panels (a) and (c), respectively. The contour maximum values in panels (b) and (d) are set as ${\pm} 2.74{u_\tau }$.

Figure 13. (a,c) Instantaneous velocity field from $R{e_\theta } \approx 800$ to 950 at two wall-normal locations: $y_{out}^ +\approx 70$ and $y_{in}^ +\approx 13.5$, respectively. (b,d) Large-scale structures retained from panels (a) and (c), respectively. The contour maximum values in panels (b) and (d) are set as ${\pm} 2.4{u_\tau }$.

Figure 14. Clear ‘foot-printing’ evidence at high $Re$ in a fully developed channel flow (reproduced from Chen & He Reference Chen and He2022). (a,c) Instantaneous velocity field from $R{e_\tau } \approx 2000$ at two wall-normal locations: $y_{out}^ +\approx 3.9\sqrt {R{e_\tau }} \approx 175$ and $y_{in}^ +\approx 13.5$, respectively. (b,d) Large-scale structures retained from panels (a) and (c), respectively. The inner and outer plane locations have the same ${y^ + }$ values as those of Marusic, Mathis & Hutchins (Reference Marusic, Mathis and Hutchins2010).

As seen from figure 12, there seems to be little correlation between the outer and inner large-scale structures filtered from the instantaneous velocity flow fields of the transitional region up to $R{e_\theta } \approx 460$. Reaching the end of the domain at a relatively higher $Re$ (figure 13), large scales start to show up more clearly in the near-wall region (figure 13d). The correlation of the large scales between the inner and outer flows however still appears to be quite low.

In direct contrast to figures 12 and 13, figure 14 illustrates strong footprints as seen in the canonical channel flow at a high $Re$ (Chen & He Reference Chen and He2022), where the fully developed turbulent flow is regarded as in an equilibrium state. The outer and inner large scales appear not only as long streaky structures but also highly correlated in both shapes and magnitudes (figure 14b versus figure 14d).

The EMD analysis above provides a relatively straightforward way to identify and correlate distinctive large length-scales between two snapshots. To gain more systematic insights into those dynamically active temporal and spatial flow structures, we also carry out some analyses using the dynamic mode decomposition (DMD) method (Schmid Reference Schmid2010). A snapshot matrix is constructed first when applying the DMD method as

(2.3a,b)\begin{equation}\boldsymbol{X}_0^{n - 1} = [{\boldsymbol{x}_0},{\boldsymbol{x}_1}, \ldots ,{\boldsymbol{x}_{n - 1}}] \in {\boldsymbol{R}^{m \times (n - 1)}},\quad \boldsymbol{X}_1^n = [{\boldsymbol{x}_1},{\boldsymbol{x}_2}, \ldots ,{\boldsymbol{x}_n}] \in {\boldsymbol{R}^{m \times (n - 1)}},\end{equation}

where in the present case, ${x_i}$ represents the instantaneous streamwise velocity field, consisting of m data points at time $t = i\Delta t$, with a sampling interval set to be $\Delta t$. The matrix $\boldsymbol{X}_0^{n - 1}$ effectively contains snapshots from $i = 0$ to $n - 1$. Assume there is a linear operator between two consecutive snapshots:

(2.4)\begin{equation}\boldsymbol{X}_1^n = \boldsymbol{AX}_0^{n - 1},\end{equation}

where $\boldsymbol{A}$ is the dynamics matrix and can be obtained by $\boldsymbol{A} = \boldsymbol{X}_1^n{(\boldsymbol{X}_0^{n - 1})^{ - 1}}$, where ${[{\cdot} ]^{ - 1}}$ denotes the pseudo-inverse operator. Given that the size of $\boldsymbol{X}$ is usually huge, the singular value decomposition is introduced (Taira et al. Reference Taira2017) to reduce the rank of $\boldsymbol{A}$ to r. The jth eigenvalue and eigenvector of A are ${\lambda _j}$ and ${\boldsymbol{\varPhi }_j}$, respectively. The dominant modes as the eigenvector and the corresponding frequency calculated by the eigenvalue can thereafter be picked up as the representatives of the raw field to reveal the dominant features of the original flow field.

To reconstruct the field, the amplitude of each mode should be first determined as $\boldsymbol{b} = [{b_1},\; {b_2}, \ldots ,{b_r}]$. The ith snapshot at time t is then given by

(2.5)\begin{equation}{\overline{\overline{\boldsymbol{x}}}_i} = \sum\limits_{j = 0}^r {\lambda _j^i{b_j}{\boldsymbol{\varPhi }_j},\; i} = t/\Delta t.\end{equation}

The loss function is introduced for comparison as the difference between the reconstructed field ${(\overline{\overline {\boldsymbol{X}_0^{n - 1}}} )_{DMD}}$ and the original field $\boldsymbol{X}_0^{n - 1}$ accumulated from $t = 0$ to $T = (n - 1)\Delta t$,

(2.6)\begin{equation}\textrm{loss} = ||\boldsymbol{X}_0^{n - 1} - {(\overline{\overline {\boldsymbol{X}_0^{n - 1}}} )_{DMD}}|{|^2}/||\boldsymbol{X}_0^{n - 1}|{|^2}.\end{equation}

For the verification purpose, the DMD method is first applied to the present full WRLES case. The aforementioned procedure is applied to $n = 201$ snapshots extracted from the x–y plane at the middle span of the domain. The streamwise range is chosen to be right downstream of the trip $(R{e_x} = 3.1\unicode{x2013} 3.6 \times {10^5})$. The convergence of the modal reconstruction is shown in figure 15, by the means of the loss of accuracy as a function of the number of modes N retained. The loss gets to less than 5 % when N is larger than 120.

Figure 15. Loss function in relation to the number of modes retained for field reconstruction N.

The reconstructed fields using 5, 25 and 50 modes and the original field are shown in figure 16. The time instant is taken at $t\; = \; 0.5T$, i.e. $t = 0.5(n - 1)\Delta t$. The reconstructed field using 50 modes and the original field are in fairly good agreement.

Figure 16. (a) Original simulated flow in the full-LES case (RFTrip-x10z7), (b–d) Reconstructed field using 5, 25 and 50 modes respectively.

There are various selection criteria to identify the dominant modes (Jovanovic, Schmid & Nichols Reference Jovanovic, Schmid and Nichols2014; Tissot et al. Reference Tissot, Cordier, Benard and Noack2014). The corresponding mode frequency is non-dimensionalised in the form of a Strouhal number $St = fL/U$, where the frequency f is calculated from the imaginary part of the eigenvalue ${\lambda _j}$, together with the trip height as the characteristic length and the freestream velocity. In the present study, the most dominant mode amplitude is selected at the primary vortex shedding frequency $St \approx 0.3$ (figure 11a) as the first mode. The dominant modes are shown in figure 17, with the second and third modes chosen at the second and third harmonics of the primary shedding frequency.

Figure 17. Dominant modes extracted from the mid x–y plane from full WRLES with frequency around (a) $St \approx 0.3$, (b) $St \approx 0.6$ and (c) $St \approx 0.9$. The blue contour colour corresponds to a negative value, while the yellow one indicates a positive value. The red dashed line indicates roughly the location of ${y^ + } \approx 13.5$, calculated with local shear velocity.

The alternating velocity components of the first mode (figure 17a) above the trip height right after the tripping element ($R{e_x}/{10^5} \approx 3.1$, $y/{\varDelta _{y,trip}} \approx 1.5$) appear as the key feature of vortex shedding in the streamwise velocity flow field. They begin to break down at the streamwise location range of $R{e_x}/{10^5} \approx 3.2$ to 3.3. Interestingly, the vortical signatures are seemingly dissected from $y/{\varDelta _{y,trip}} \approx 1$, either keeping convected downstream or propagated towards the near-wall region. There is a clear difference in magnitude between the downstream convected and the near-wall propagated. Thus, the outer-flow structures do not seem to be able to penetrate into the inner layer with a significant effect. Note that in Agostini & Leschziner (Reference Agostini and Leschziner2016), the magnitudes of the outer coherent large scales are identified to be similar to those of their ‘footprints’ in the near-wall region. In the present work, we see no clear evidence of such ‘footprints’ from the outer-flow region.

In relation to the streamwise evolution, we can also see that the large-scale trip-associated disturbances start to decay soon after being generated. The largest scales in the most dominant mode (figure 17a) decay the fastest. The shorter scales are generated later, but also seem to decay at a slower pace, as shown in figures 17(b) and 17(c). This is in clear contrast to the ‘footprinting’ features in a well-developed turbulent flow where large coherent streaky structures persist with a long lifetime.

Overall, the excessive energy observed in the outer flow region (figure 10) is clearly associated with the shed vortices over the tripping element. These vortical flow structures residing in the log-region make the outer flow take longer to settle down to a complete equilibrium state. However, these tripping-induced large-scale residual disturbances share little similarity with the coherent large structures of long meandering streamwise streaks in a fully developed turbulent flow. The tripping-induced large-scale disturbances clearly decay in the streamwise direction, seemingly at a faster rate than those of shorter wavelengths. These large scales also do not show clearly discernible ‘footprints’ on the near-wall region. The characteristics as observed may, to some extent, explain why the near-wall region is usually observed to experience a shorter development distance than the outer flow region before reaching an equilibrium state (Devenport & Lowe Reference Devenport and Lowe2022). The related physical understanding also provides a useful basis for implementing the two-scale method for a tripped turbulent boundary layer, more particularly for sizing a locally embedded fine-mesh block for resolving the trip itself in relation to other fine-mesh blocks embedded in the near-wall region for the rest of the boundary layer.

3.1. Two-scale methodology for turbulent boundary layer flows

The local fine-mesh DNS blocks embedded in a near-wall region, as illustrated in figure 18(a), are generated by subdividing corresponding coarse-mesh cells. At the inlet to the computational domain, a laminar Blasius profile can be specified to reduce the domain size for the laminar part. The distance from the inlet to the trip is kept around 20 tripping element streamwise lengths. No identifiable influence of the trip in the form of upstream propagated pressure waves on the inlet velocity profile is observed, thus the distance is regarded as adequately long. The whole computational domain can be divided into three parts as shown in figure 18(b): the global outer flow region (coarse-mesh) where the large-scale structures are directly resolved on the base coarse-mesh, the local near-wall DNS blocks (fine-mesh) and the global inner region (coarse-mesh). The key working of the present two-scale method is that the solution of the under-resolved global inner coarse-mesh region will be corrected using the time-invariant source terms originated from the local fine-mesh blocks and propagated spatially by a block-spectral mapping.

Figure 18. Illustration of computational domain with the near-wall embedded DNS blocks. (a) Overall view of the configuration (the solid magenta line indicates the trip). (b) Close-up x–y cut plane view.

The interface flux conservation between the coarse- and fine-mesh regions is achieved through the non-conforming arbitrary mesh interface (AMI) (Farrell & Maddison Reference Farrell and Maddison2011). It is applied directly as the interface condition to the top surface of each embedded block. For a pair of side faces of a block embedded in the developed turbulent region (either ${F_{x,a2}}$ and ${F_{x,b2}}$, or ${F_{z,a2}}$ and ${F_{z,b2}}$, as marked in figure 18a), the scale-dependent interface treatment (He Reference He2018) is adopted. A local flow variable is decomposed into a coarse-mesh base value and a fine-mesh perturbation:

(3.1)\begin{equation}\boldsymbol{u}(\boldsymbol{x},t) = {\boldsymbol{u}_C}(\boldsymbol{x},t) + {\boldsymbol{u^{\prime\prime}}_f}(\boldsymbol{x},t),\end{equation}

where ${\boldsymbol{u}_C}(\boldsymbol{x},t)$ is the coarse-mesh resolved, obtained directly via the baseline AMI. For the fine-scale fluctuation part ${\boldsymbol{u^{\prime\prime}}_f}(\boldsymbol{x},t)$, the periodic condition is applied. As an example of the pair faces of the block, ${F_{x,a2}}$ and ${F_{x,b2}}$, we have

(3.2) \begin{gather} \left.{\begin{gathered} {{{\boldsymbol{u^{\prime\prime}}}_f}{{(x,t)}_{x,a2}} = {{\boldsymbol{u^{\prime\prime}}}_f}{{(x,t)}_{x,b2}}\; }\\ {\boldsymbol{u}{{(\boldsymbol{x},t)}_{x,a2}} = {\boldsymbol{u}_C}{{(\boldsymbol{x},t)}_{x,a2}} + {{\boldsymbol{u^{\prime\prime}}}_f}{{(\boldsymbol{x},t)}_{x,b2}}}\\ {\boldsymbol{u}{{(\boldsymbol{x},t)}_{x,b2}} = {\boldsymbol{u}_C}{{(\boldsymbol{x},t)}_{x,b2}} + {{\boldsymbol{u^{\prime\prime}}}_f}{{(\boldsymbol{x},t)}_{x,a2}}} \end{gathered}}\right\}.\end{gather}

The scale-dependent interface treatment can also be applied to the spanwise pairing points of the frontal block (${F_{z,a1}}$ and ${F_{z,b1}}$). The upstream boundary of the frontal block $({F_{x,a1}})$ is consistent with the global region, subject to the Blasius laminar velocity profile. The downstream boundary of the frontal block $({F_{x,b1}})$ is simply dealt with by the AMI treatment.

The principal working of the two-scale methodology consists of two related aspects: (a) generating the source terms from the local fine-mesh blocks and (b) propagating the source terms to the global inner coarse-mesh region (figure 18b). Consider the flow governing equations in a simple form for flow variable vector $\boldsymbol{u}$:

(3.3)\begin{equation}\frac{{\partial \boldsymbol{u}}}{{\partial t}} + R(\boldsymbol{u}) = 0,\end{equation}

which in its original form will be directly solved numerically in the global outer domain and local inner blocks (figure 18b). For the under-resolved global coarse-mesh inner domain in the near-wall region (figure 18b), the upscaling is introduced leading to the augmented equations with extra source terms to drive the coarse-mesh solution toward the target solution equivalently subject to a fine-mesh locally.

The upscaling is facilitated through space–time averaging. For a coarse-mesh cell embedded with ${n_f}$ fine-mesh cells, the space–time averaged flow variable can be simply defined as the local volume-average of time-averaged variables of ${n_f}$ fine-mesh cells:

(3.4)\begin{equation}{\widetilde {\bar{\boldsymbol{u}}}_f}(\boldsymbol{x}) = \frac{1}{{\Delta {v_C}}}\sum\limits_{i = 1}^{{n_f}} {{{[{{\bar{\boldsymbol{u}}}_f}(\boldsymbol{x})\Delta {v_f}]}_i}} ,\end{equation}

where $\Delta {v_C}$ is the volume of the coarse-mesh cell, $\Delta {v_C} = \sum\nolimits_{i = 1}^{{n_f}} {{{(\Delta {v_f})}_i}} $. The overbars ‘–’ and ‘${\sim} $’ denote the time-averaged and the spatial-averaged values, respectively. The upscaled equations in the coarse-mesh domain become

(3.5)\begin{equation}\frac{{\partial {\boldsymbol{u}_{\boldsymbol{C}}}}}{{\partial t}} + R({\boldsymbol{u}_{\boldsymbol{C}}}) = \boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{st}}}(\boldsymbol{x}).\end{equation}

The source term $\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{st}}}(\boldsymbol{x})$ is time-invariant and consists of two parts (He Reference He2018):

(3.6)\begin{equation}\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{st}}} = {(\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{st}}})_f} + {(\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{t}}})_C}.\end{equation}

The first part is generated by simply computing the net flux residuals on the coarse-mesh in a discrete form with the space–time averaged variables from the corresponding fine-mesh solution:

(3.7)\begin{equation}{(\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{st}}})_f} = R({\widetilde {\bar{\boldsymbol{u}}}_f}).\end{equation}

The second part exists due to the nonlinear time-averaging effect of the unsteady solution on the coarse-mesh. It should be accounted for so that the upscaled equations for the targeted space–time-averaged solution in (3.5) can be balanced:

(3.8)\begin{equation}{(\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{t}}})_C} = \overline {R({\boldsymbol{u}_{\boldsymbol{C}}})} - R({\bar{\boldsymbol{u}}_C}).\end{equation}

When a scale-resolving flow solution in the global inner region is statistically converged, the local time-mean coarse-mesh solution ${\bar{\boldsymbol{u}}_C}$ should converge to the target space–time averaged fine-mesh solution ${\widetilde {\bar{\boldsymbol{u}}}_f}$ as intended. For more detailed formulations, the reader is referred to He (Reference He2018, Reference He2021) and Chen & He (Reference Chen and He2022).

Once the locally sampled source terms are generated, they need to be effectively propagated to the global coarse-mesh domain spatially in a wall-parallel direction with inhomogeneity. He (Reference He2018) resorts to constructing and propagating the global spatial variations in the two wall-parallel directions through a Fourier spectral mapping for a periodic domain or a half of the Fourier spectrum for a non-periodic domain. In the present work, the Chebyshev spectral method is used for the non-periodic TBL flow in the streamwise direction. The orthogonal group of Chebyshev polynomials are constructed as

(3.9)\begin{equation}{T_n}(x) = \cos [n \times \textrm{arccos}(x)],\quad n = 1,2, \ldots .\end{equation}

The sampling points on the transformed interval [−1, 1] are based on the Chebyshev–Gaus–Lobatto points as

(3.10)\begin{equation}{x_n} ={-} \cos \left( {n \times \frac{{\rm \pi} }{N}} \right),\quad n = 0,1,2, \ldots ,N.\end{equation}

The corresponding values at sampling points are ${f_n}({x_n}),\;n = 0,1, \ldots ,N$. The globally mapped variable ${f_{map}}$ can be expressed in terms of the basis functions ${T_n}(x)$ with coefficients ${C_n}$ as

(3.11)\begin{equation}{f_{map}}(x) = \sum\limits_{n = 0}^{N - 1} {{C_n}{T_n}(x)} .\end{equation}

Herein, the discrete cosine transform can be used to acquire the coefficients $C = {[{C_1},{C_2}, \ldots ,{C_N}]^T}$ from the locally sampled function values $f = {[{f_1},{f_2}, \ldots ,{f_N}]^T}$ as

(3.12)\begin{equation}\boldsymbol{C} = \boldsymbol{\varLambda f},\end{equation}

where the constructor matrix $\boldsymbol{\varLambda }$ is a linear system of $N\; + \; 1$ dimensions (Trefethen Reference Trefethen2000). The constructor matrix can be expressed as

(3.13) \begin{gather}{\varLambda _{nm}} = \begin{cases} {\dfrac{1}{N}}&{m = 0}\\ {\dfrac{2}{N}\cos \dfrac{{nm {\rm \pi}}}{N}}&{m = 1,\; \ldots ,\; N - 1}\\ {\dfrac{1}{N}\cos (n{{\rm \pi} })}&{m = N} \end{cases}\end{gather}

A flowchart showing the main parts of the present two-scale block spectral method during one time-marching step is given in figure 19, as implemented in the OpenFOAM for incompressible flows.

Figure 19. Flowchart of the two-scale block-spectral solution process as implemented.

3.2. Sizing locally embedded fine-mesh block for tripping

Chen & He (Reference Chen and He2022) discussed the sizing of the embedded block for a turbulent channel flow, chiefly based on the energy spectra generated from the DNS database (Lee & Moser Reference Lee and Moser2015) with the analysis of the length scales of the near-wall ‘universal’ portion as well as in relation to the inner region in terms of the low bound wall-normal distance of the log-law region, based on Marusic et al. (Reference Marusic, Monty, Hultmark and Smits2013). The corresponding block sizes should also be applicable to a fine-mesh block embedded in a fully turbulent part of a boundary layer after the trip. Thus, the streamwise and spanwise block lengths (figure 18a) should accordingly be $l_{x2}^ +\approx 3500$ and $l_{z2}^ +\approx 350$ respectively, where ‘+’ denotes the wall units based on the local shear velocity ${u_\tau }$ and viscosity $\nu $. The wall-normal height of the block is chosen to be $l_{y2}^ +\approx 75$ in the turbulent region, consistent with the baseline WRLES study presented in § 2 and that of Chen & He (Reference Chen and He2022).

For sizing the frontal fine-mesh block for the tripping, extra considerations should be given particularly in the context of the physical characteristics, as identified and discussed in § 2. First, regarding the streamwise length of the fine-mesh block, it is required to cover the key streamwise development phase with essential features of the break-down process of the initial vortices shed from the tripping element. Given the observations made in § 2 on how far downstream the laminar-to-turbulent transition would be largely completed, the streamwise length is chosen to reach $R{e_\theta } \approx 650$ in the present case study.

Second, the wall-normal height of the fine-mesh block for the tripping should cover and resolve the trip element itself as well as the trip-generated large-scale vortical disturbances. Those large-scale vortical disturbances, once generated, would evolve streamwise in the lower part of the log-region which should be covered by the base coarse-mesh. However, the near-trip vortex shedding process, responsible for generating these large-scale disturbances, needs to be resolved locally with the fine-mesh resolution. As such, the wall-normal distance of the frontal fine-mesh block for tripping should be higher than the low bound of the log-region, which is used to determine the fine-mesh block height in a fully turbulent regime, as shown by Chen & He (Reference Chen and He2022). This requirement can be easily met however if we take a value of the log region low-bound at a downstream position (thus a more developed and thicker boundary layer). In the present work, we take the wall-normal distance of the frontal block for the tripping to be the same as that of the second one for the tripped boundary layer, $l_{y1}^ += l_{y2}^ +\approx 75$. This should provide sufficient coverage for resolving the generation of those tripping-associated large-scale vortical disturbances evolving in the outer flow region of the tripped turbulent boundary layer.

Third, the sensitivity to the spanwise length of the fine-mesh tripping block ${l_{z1}}$ is studied through several test cases shown in table 3. The ‘Full-Span’ case has the frontal fine-mesh block span ${l_{z1}}$ covering the full spanwise width ${L_z}$. The ‘Half-Span’ case reduces the frontal block span by half to ${L_z}/2$. The minimum width of ${l_{z1}}$ (‘Mini-Span’) here is chosen to be the same size as that of the second fine-mesh block in the tripped TBL ${l_{z2}}$. The spectra of flow velocity taken from a numerical probe placed at a location just downstream of the trip height show clearly stand-out peaks (figure 20). These peaks correspond to the dominant shedding frequency and its second and third harmonics at the corresponding Strouhal number $St$, approximately 0.3, 0.6 and 0.9 as discussed in § 2. At the primary shedding frequency $(St\; \approx 0.3)$, the energy peak is well captured in all three cases. At higher frequencies with much less energy (less than 1 % of the dominant peak), the Half-Span and Mini-Span are seemingly subject to more broadband disturbances.

Table 3. Frontal block widths of the cases simulated.

Figure 20. Turbulence power spectra density (PSD) with respect to the non-dimensional frequency, Strouhal number $St$. The black, magenta and blue lines are for the Full-Span, the Half-Span and the Mini-Span results, respectively.

How would these detailed differences at the trip between different spanwise block sizes affect the downstream flow? Flow statistics including the mean velocity profiles and the fluctuations are extracted at two locations, as shown in figure 21. The first location is in the rear part of the frontal fine-mesh block at $R{e_\theta } \approx 600$ and the second location is in the rear part of the second embedded fine-mesh block in the fully turbulent region at $R{e_\theta } \approx 800$. Also shown (figure 22) are the energy spectra at these two locations. We see good agreement in the mean statistics and the energy spectra at both locations between the Full-Span, Half-Span and Mini-Span sizes. For the spanwise size range tested, the downstream flow field seems largely insensitive to the spanwise size of the frontal block ${l_{z1}}$. The rest of the two-scale calculations is thus all taken with the Mini-Span for the frontal embedded fine-mesh block.

Figure 21. Mean velocity profiles as functions of wall-normal distance ${y^ + }$ at (a) $R{e_\theta } \approx 600$ and (b) $R{e_\theta } \approx 800$. Profiles of streamwise velocity fluctuations as functions of ${y^ + }$ at (c) $R{e_\theta } \approx 600$ and (d) $R{e_\theta } \approx 800$. Full-Span, solid black lines; Half-Span, solid magenta squares; Mini-Span, open blue circles. The grey lines mark the inner–outer mesh interface. Note that not all data are shown for clarity.

Figure 22. PSD as a function of non-dimensional frequency k: (a) $R{e_\theta } \approx 600$ and (b) $R{e_\theta } \approx 800$. Full-Span, black lines; Half-Span, magenta lines; Mini-Span, blue lines. The energy spectra are taken at ${y^ + }\sim 13.5$ as the peak energy location in figure 21(c,d).

To see more closely how the large-scale disturbances are triggered and evolve downstream in the embedded block, we now apply the DMD modal analysis here, in the same procedure as in § 2. We compare the mid-plane of the local fine-mesh block with a reduced-span to that of the full LES solutions. The streamwise coverage range is from the trip $(x = 0.1{l_{x1}})$ to the end of the frontal block $(x = {l_{x1}})$. The comparison between the two cases is shown in figure 23. The transition process from shed vortical disturbances to the breakdown and formation of a new TBL is well captured in the embedded block, shown in figure 23(b). The detailed features in the vicinity of the trip are slightly different, e.g. a slightly later break-down process for the full LES (figure 23a) compared to that for the locally embedded block (figure 23b). Downstream, however, the differences between the two seem rather insignificant with the rapid decay of the large-scale disturbances towards the end of the block. The overall picture of the streamwise evolution for the trip-induced large scales is in line with the observed insensitivities in the mean statics and energy spectra at $R{e_\theta } \approx 600$, indicating the adequate length of the locally embedded fine-mesh block for capturing and resolving the tripping.

Figure 23. Comparison of large-scale disturbances $(St \approx 0.3)$ on the mid x–y plane from (a) full LES and (b) the local DNS block.

3.3. Two-scale solution for tripped TBL

3.3.1. Validation against baseline full LES

A key variable sensitive to the near-wall mesh resolution is the wall shear stress, as shown in the mesh-dependence study in § 2.2. The accurate calculation of the wall shear stress is of great importance in many spatially developing flows (Launder & Spalding Reference Launder and Spalding1974). The computed friction coefficient ${C_f}$ is shown in figure 24 as the function of the local Reynolds number. The full LES result is also shown as the baseline reference which has been validated against well-established experimental results and DNS databases, as presented in § 2. The direct solution with locally embedded blocks but without the source-term coupling is labelled as the ‘one-scale’ solution, showing a significant discrepancy, because of the under-resolution of the near-wall coarse-mesh. In the present two-scale solution, the mean flow errors are markedly reduced by the source terms propagated through BS mapping.

Figure 24. Friction coefficient ${C_f}$ as the function of Reynolds number $R{e_x}$. The blue solid line is for the present two-scale block-spectral solution in the global coarse region. The solid black squares indicate the full LES results. The red dash line shows the ‘one-scale’ solution in the global coarse-mesh region without the source-term coupling.

The time-mean boundary layer velocity profiles on the global coarse-mesh at the exit of the domain where $R{e_\theta } \approx 1000$ are compared, as shown in figure 25. The under-resolved ‘one-scale’ solution with no source-term coupling gives poor results with large errors (~19 %). The impact of the source term correction is clearly underlined by the two-scale BS solution agreeing well with the DNS results (Schlatter & Orlu Reference Schlatter and Orlu2012).

Figure 25. Mean velocity profile in relation to wall normal distance ${y^ + }$ in the global coarse-mesh region. The present two-scale BS solution (blue squares) is compared with the DNS result from Schlatter & Orlu (Reference Schlatter and Orlu2012) (black solid lines) at the Reynolds number $R{e_\theta } \approx 1000$. The red dash line is for the ‘one-scale’ solution without the source-term coupling.

It should be noted that the source term mapping for the tripping region covered by the first embedded fine-mesh block is different from that for the rest of the TBL. As shown previously in § 2, the trip-induced transition process is subject to large-scale disturbances of large amplitudes highly interactive particularly shortly downstream of the trip element. As a result, the region from $R{e_x}/{10^5} \approx 3$ to 4 (figure 24) experiences a much higher gradient and stronger local history effect in the streamwise direction than the rest of the TBL. Consequently, for the frontal tripping region, the source terms generated in the first fine-mesh block are only propagated in the spanwise direction. They are directly mapped to the full-span coarse-mesh region. The source term components in the x, y and z directions, $S{T_x}$, $S{T_y}$ and $S{T_z}$, are shown in figure 26, where dash lines indicate those source terms generated from the frontal fine-mesh tripping block to be mapped directly to the full span of the corresponding coarse-mesh region.

Figure 26. Distribution of source terms $\boldsymbol{ST} = (S{T_x},\; S{T_y},S{T_z})$ along the streamwise direction with respect to local Reynolds numbers $R{e_x}$ in the viscous sublayer $({y^ + } \approx 6)$. The blue dash lines indicate $S{T_x}$ in the frontal trpping block with a spanwise mapping only. The solid blue lines indicate the spectrally mapped results from the Chebyshev method for the rest of TBL. The solid blue triangles mark the sample points. The other two scalar components $S{T_y}$ and $S{T_z}$ are shown in yellow and green, respectively.

From the streamwise location $R{e_x}/{10^5} \approx 4$ onwards, the streamwise change of the flow development and the wall shear stress is much more gradual as it enters a fully turbulent region. The source term distribution along the streamwise direction also becomes much smoother. The Chebyshev spectral method introduced in § 3.1 is applied here for the source-term propagation. Two sample points for constructing the streamwise Chebyshev spectrum in this case are marked by triangle symbols in figure 26. One is at the rear portion of the frontal fine-mesh block and the other is at the rear part of the second fine-mesh block embedded in the fully turbulent region.

The mean velocity and the fluctuations in rms in the near-wall fine-mesh blocks and adjacent coarse-mesh region are compared with the full LES results in good agreement at two streamwise locations, as shown in figure 27. Figure 27(a,b) are for the results at the rear part of the frontal block at $R{e_x}/{10^5} \approx 4$ with the corresponding $R{e_\theta } \approx 600$. Figure 27(c,d) are for the results at the rear part of the second block at $R{e_x}/{10^5} \approx 5.7$ with the corresponding $R{e_\theta } \approx 800$.

Figure 27. Mean statistics including (a,c) the mean velocity profiles and (b,d) streamwise velocity fluctuations with respect to the wall-normal distance. The present two-scale BS solutions (blue circles) are compared with the full LES solutions (black solid lines) at two Reynolds numbers: (a,b) $R{e_\theta } \approx 600$ and (c,d) $R{e_\theta } \approx 800$. The grey lines mark the inner–outer mesh interfaces. Note that not all data are shown for clarity.

Shown in figure 28 are the energy spectra at these two Reynolds numbers in the fully turbulent regime. The spectra are taken in the near-wall region at ${y^ + } \approx 13.5$. The dimensionless frequency is again calculated based on the boundary layer thickness and wall shear velocity: $k = f\delta /{u_\tau }$. First, note that the results from the present two-scale method overlap with the full-LES spectra. Second, the PSD taken from the embedded region confirms a full coverage of the spectral range in the local fine-mesh block. More remarkably, the fine-mesh spectrum overlaps smoothly with the coarse-mesh PSD at lower frequencies without any spectral gap. This is attributed to the local embedded fine-mesh block receiving low-frequency large-scale signals directly from the global coarse-mesh domain, thanks to the scale-dependent interface method. The global coarse mesh is capable of resolving the large scales, but experiences high numerical dissipations for shorter wavelengths in the higher frequency range, as expected. The corresponding coarse-mesh under-resolution in the near-wall region is corrected by the source terms to drive towards the target time-mean flow. The present observation for a tripped TBL, where a smooth coverage of full turbulence spectrum without a spectral gap or scale separation can be achieved in a locally embedded fine-mesh block, is consistent with that for a canonical channel flow made by Chen & He (Reference Chen and He2022). It underscores the advantage of the scale-dependent interface treatment over the direct periodic condition commonly adopted in previous MFU-based methods. In the previous MFU methods, the spatial periodicity of the unit length would have to artificially truncate the local large-scale disturbances in the near-wall region footprinted by those large-scale coherent structures residing in the outer flow region.

Figure 28. PSD with respect to the dimensionless frequencies. The present two-scale BS solutions (blue lines) are compared with the full LES solutions (green lines) at two Reynolds numbers: (a) $R{e_\theta } \approx 600$ and (b) $R{e_\theta } \approx 800$. Red lines are for the results in the global coarse-mesh inner region. The energy spectra are taken at ${y^ + }\sim 13.5$.

3.3.2. Cases with higher Reynolds numbers

For further validation and demonstration, we have two additional cases simulated at higher Reynolds numbers ($R{e_\theta } \approx 1500$ and 2600) at the end of the TBL. The simulation parameters are shown in table 4. Both cases have the same frontal block embedded as that in § 3.3.1 and the second block embedded close to the exit. All configurations in wall units are normalized by the local shear velocity ${u_\tau }$ at the specific Reynolds numbers $R{e_\theta }$. The fluid domain simulated is $\varOmega \; = \; [0,\; {L_x}] \times [0,\; {L_y}] \times [0,\; {L_z}]$, with the domain used in § 3.3.1 denoted as $\varOmega ^{\prime} = [0,\; {L^{\prime}_x}] \times [0,\; {L^{\prime}_y}] \times [0,\; {L^{\prime}_z}]$ for comparison purposes.

Table 4. Parameters of simulation at higher $Re$.

Validations of the solution accuracy are shown in figures 29 and 30. Figure 29 presents the friction coefficient ${C_f}$ distribution with respect to the Reynolds number $R{e_\theta }$. Figure 30 shows the mean statistics at the highest Reynolds number calculated $(R{e_\theta } \approx 2600)$ with the local DNS block embedded. All two-scale solutions match well with the correlated curves (Smits et al. Reference Smits, Matheson and Joubert1983) and the DNS data (Schlatter & Orlu Reference Schlatter and Orlu2010).

Figure 29. Friction coefficient ${C_f}$ as a function of Reynolds number $R{e_\theta }$. The solid grey lines indicate the correlated curve by Smits et al. (Reference Smits, Matheson and Joubert1983). The solid blue circles are for the present two-scale BS solution. The red dashed lines are for the one-scale solution without the source-term coupling. Both results are taken from the global coarse-mesh region. Note that not all data points are shown for clarity.

Figure 30. Mean statistics at Reynolds number $R{e_\theta } \approx 2540$: (a) mean velocity profiles and (b) streamwise velocity fluctuations with respect to the wall-normal distance. The present two-scale block-spectral solutions (open blue squares) are compared with the DNS results (Schlatter & Orlu Reference Schlatter and Orlu2010) (black solid lines). The grey lines mark the inner–outer mesh interfaces. Note that not all data points are shown for clarity.

3.4. Mesh count– $Re$ scaling

Finally, we examine the mesh count–Reynolds number scaling for the spatially developing TBL to evaluate the potential benefit in terms of computational cost for the present two-scale method. As a laminar-to-turbulent transition (regardless of the way of initiation) should occur at a largely fixed Reynolds number, we shall only consider a fully turbulent part of the boundary layer starting from an upstream station ${x_0}$ to a downstream one $x = {L_x}$, where the Reynolds number based on either the total length or the maximum boundary layer thickness is calculated.

We assume that the whole TBL can be divided into ${N_s}$ subdomains, as shown in figure 31(a). Within each subdomain, we have one locally embedded near-wall fine-mesh block. As such, a subdomain can be approximately viewed as a channel flow with the local mean boundary layer thickness $\delta $ being taken as the half channel height.

Figure 31. Sectioned TBL domain: (a) illustration of the TBL domain with ${N_s}$ sub-domains; (b) a close-up view of the sub-domain i.

The total mesh count required for the whole TBL $({x_0} \le x \le {L_x})$ is

(3.14)\begin{equation}{N_{total}} = \sum\limits_{i = 1}^{{N_s}} {{N_i}(R{e_{{\delta _i}}})} ,\end{equation}

where the total mesh count ${N_{total}}$ is a summation of local mesh counts for ${N_s}$ subdomains (figure 31a). For sub-domain i, its mesh count ${N_i}$ depends on the local Reynolds number $R{e_{{\delta _\textrm{i}}}}$ based on the freestream velocity ${U_\infty }$, the kinematic viscosity $\nu $ and the local boundary layer thickness ${\delta _i}$. For each sub-domain, the boundary layer thickness can be approximately taken as a constant. Thus, a local TBL flow in a sub-domain is equivalent to a channel flow (figure 31b).

Following Choi & Moin (Reference Choi and Moin2012), the local mesh count is estimated as $NLE{S_{local}}\sim Re_\delta ^{11/6}$ for wall-resolved LES and $NDN{S_{local}}\sim Re_\delta ^{11/4}$ for DNS. For the present two-scale method, there should be one embedded fine-mesh DNS block in each subdomain. As indicated above, the mesh count ${N_i}(R{e_\delta })$ for subdomain i can be estimated similarly to a channel flow. Chen & He (Reference Chen and He2022) estimated the mesh count for a channel flow with $\delta $ being the half channel height as

(3.15) \begin{align}\begin{aligned} NLE{S_{channel}} & = NLE{S_{fine\;inner}} + NLE{S_{coarse\;outer\; }}\\ & \quad + NLE{S_{coarse\;inner}}\sim \; {\delta ^ + }. \end{aligned}\end{align}

The three terms on the right-hand side of (3.15) corresponding to the three domains (figure 31b) are estimated as

(3.16) \begin{equation}\left. {\begin{gathered} {NLE{S_{fine\;inner}} = NDN{S_{embed\;block}}\; \sim \; y_s^ + \; \sim \; {\delta^{ + 0.5}}}\\ {NLE{S_{coarse\;outer\; }}\sim {{\int_{{y_s}}^\delta {{\delta^2}} } / {{y^3}\,\textrm{d}y}}\sim {\delta^ + }}\\ {NLE{S_{coarse\;inner}}\; \sim \; {\delta^ + }} \end{gathered}} \right\},\end{equation}

where ${y_s}$ marks the starting point of the log region, being proportional to ${\delta ^ + }^{0.5}$ (Marusic et al. Reference Marusic, Monty, Hultmark and Smits2013). In the present TBL case, we replace $\delta $ denoting the half channel-height for channel flow by the local boundary layer thickness ${\delta _i}$. Choi & Moin (Reference Choi and Moin2012) represented scaling with different Reynolds number definitions through the conversion, $R{e_\delta }\sim {({\delta ^ + })^{12/11}}$. Thus, correspondingly, (3.15) can be rewritten in terms of local TBL thickness for subdomain i as

(3.17)\begin{equation}NLE{S_i}(R{e_{{\delta _i}}})\sim \; Re_{{\delta _i}}^{11/12}.\end{equation}

Furthermore, we assume the local Reynolds number for subdomain i can be approximately correlated to the overall Reynolds number $R{e_\delta }$ based on the downstream boundary layer thickness by $R{e_{{\delta _i}}} = {C_i}R{e_\delta }$, where ${C_i}$ is a multiplier coefficient for subdomain i $(0 < {C_i} < 1)$. The multiplier coefficient should remain roughly constant when the TBL Reynolds number $R{e_\delta }$ varies. Equation (3.14) can then be expressed as

(3.18)\begin{equation}NLE{S_{total}}\sim \; Re_\delta ^{11/12}\sum\limits_{i = 1}^{{N_s}} {{f_i}({C_i})} .\end{equation}

As ${f_i}$ for subdomain i is only a function of the local multiplier ${C_i}$ independent of Reynolds number, we will then have

(3.19)\begin{equation}NLE{S_{total}}\sim \; Re_\delta ^{11/12}.\end{equation}

Therefore, the overall mesh count for TBL can now be reduced from $O(Re_\delta ^{11/6})$ for the WRLES to $(Re_\delta ^{11/12})$ for the two-scale solutions. Figure 32 shows the mesh count–$Re$ scaling for the different solution approaches.

Figure 32. Mesh count–$Re$ scaling. DNS (black dots); wall-resolved LES (red dashed line); the present two-scale method (blue solid line). The open black square and red circle are for the actual ${N_{total}}$, as reviewed in Deck et al. (Reference Deck, Renard, Laraufie and Sagaut2014). Solid red circle is for the full WRLES used in the present case study (§ 2); solid blue squares are for the present two-scale calculations (§ 3).