Two-scale solution for tripped turbulent boundary layer

Abstract

Recent findings on the Reynolds-number-dependent behaviour of near-wall turbulence in terms of the ‘foot-printing’ of outer large-scale structures call for a new modelling development. A two-scale framework was proposed to couple a local fine-mesh solution with a global coarse-mesh solution by He (Intl J. Numer. Meth. Fluids, vol. 86, 2018, pp. 655–677). The methodology was implemented and demonstrated by Chen & He (J. Fluid Mech, vol. 933, 2022, p. A47) for a canonical turbulent channel flow, where the mesh-count scaling with Reynolds number is potentially reduced from $O(R{e^2})$ for a conventional wall-resolved large-eddy simulation (WRLES) to $O(R{e^1})$. The present work extends the two-scale method to turbulent boundary layers. A two-dimensional roughness element is used to trip a turbulent boundary layer. It is observed that large-scale disturbances originating at the trip have a much shorter lifetime and weaker foot-printing signatures on near-wall flow compared to those long streaky coherent structures in well-developed wall-bounded turbulent flows. Modal analyses show that the impact of trip-induced large scales can be adequately captured by a locally embedded fine-mesh block. For the tripped turbulent boundary layer, a Chebyshev block-spectral mapping is adopted to propagate source terms from the local fine-mesh blocks to the global coarse-mesh domain, driving to a target solution for the upscaled equations. The computed mean statistics and energy spectra are in good agreement with corresponding experimental data, WRLES and direct numerical simulation (DNS) results. The overall mesh count–$Re$ scaling is estimated to reduce from $O(R{e^{1.8}})$ for the full wall-resolved LES to $O(R{e^{0.9}})$ for the present two-scale solution.

1. Introduction

1.1. General background

Numerical simulations of wall-bounded turbulent flows are challenging due to the wide range of length scales involved, from very large eddies away from the wall down to very small Kolmogorov scales in the near-wall region, where the required resolution is the highest attributed to the spatial and temporal micro-scales (Moin 1997; Jimenez & Moser 2000). The near-wall region is thus the focus regarding computational costs for wall-bounded turbulent flows. A higher Reynolds number will result in a more disparate scale difference. As a result, the mesh count needed to resolve all scales of turbulence, as in a direct numerical simulation (DNS), scales with the Reynolds number as $O(R{e^3})$ (Jimenez 2003). For a large-eddy simulation (LES) which filters out the small scales to be modelled, the required mesh count still scales up and increases with the Reynolds number as $O(R{e^2})$ (Mizuno & Jiménez 2013), and thus is still too expensive for many practical engineering applications.

A prevalent framework to reduce the overall computational cost for simulations takes advantage of the near-wall ‘universal’ behaviour where turbulence behaves in a self-sustained ‘autonomous’ manner. The underlying wisdom promotes a hybrid approach to couple the near-wall modelled Reynold-averaged Navier–Stokes (RANS) part with the scale-resolving outer flow (Spalart et al. 1997; Cabot & Moin 2000; Piomelli & Balaras 2002; Spalart 2009; Larsson et al. 2015; Bose & Park 2018). However, the significant findings recently made in understanding wall-bounded turbulent flows challenge the existing modelling approaches. The near-wall turbulence characteristics are revealed not to be completely ‘universal’, but instead are Reynolds-number-dependent with the ‘footprints’ of large-scale coherent structures of the outer flow region, which are also ‘modulated’ in the near-wall region (Hutchins & Marusic 2007; Jimenez 2013; Lee & Moser 2015).

A distinctive alternative approach the authors adopt is the two-scale methodology (He 2018; Chen & He 2022). The method is to couple a locally embedded near-wall fine-mesh DNS block (or a small number of blocks) with a global coarser mesh domain. The influence of large-scale structures on the local fine-mesh block is captured with a scale-dependent coarse–fine mesh interface treatment originated by He (2018), so that the local fine-mesh region can receive the large-scale signals from the global coarse-mesh domain. The global under-resolved near-wall coarse-mesh region outside the local fine-mesh block is governed by the augmented flow governing equations with additional source terms. The source terms are generated by upscaling the target space–time-averaged solution within the locally embedded fine-mesh block. The two-scale methodology has been implemented, verified and demonstrated in the canonical channel flow by Chen & He (2022) where the mesh-count scaling with Reynolds number is potentially reduced from $O(R{e^2})$ for the full wall-resolved LES to $O(R{e^1})$ for the two-scale solutions.

Given the homogeneity in both streamwise and spanwise directions as the basic limitation of the canonical channel flow, it will be naturally of interest to examine and demonstrate the feasibility of the two-scale method for a spatially developing flow. The present work is therefore focused on extending the two-scale method to a turbulent boundary layer. There are two related major issues to be addressed. First, how do we start a turbulent boundary layer (TBL) with clearly and consistently applied inflow conditioning? Second, how can we properly implement the two-scale method in the spatially developing boundary layer to achieve a substantial computational performance gain through enacting and taking advantage of the local fine-mesh solutions?

1.2. Inflow conditioning for turbulent boundary layer

The pre-requisite to numerically studying the TBL flow is to have a properly started turbulent flow regime. A popular way to avoid a much-too-long transition route from laminar to turbulent (Sillero, Jimenez & Moser 2013) is to recycle and rescale the downstream characteristics to generate synthetic turbulent inflow conditions (Lund, Wu & Squires 1998). In so doing, the aimed Reynolds number can be reached efficiently with the inflow turbulent profile rescaled from downstream. However, these forms of non-physical approximation introduce a certain level of arbitrariness. The artificial inflow does not reveal the true history of turbulence development, instead, it depends on the future evolution from which it itself should start (Schlatter & Orlu 2010). Other synthetic inflow generation methods include those built on theoretical properties of a well-developed wall-bounded turbulent flow (Kraichnan 1970) and empirical information (Jarrin et al. 2006; Adler et al. 2018). After all, these methods all facilitate efficient calculations in truncated domains, whilst they lack physical coherence and usually require a long distance for the artificial history effects to be sufficiently diminished.

Schlatter & Orlu (2012) introduced an artificial tripping method with a function constructed to exert a volume forcing to trigger a rapid transition. By tuning the representative parameters, a fully developed boundary layer may be generated seemingly almost immediately after the trip, with ‘optimal configurations’ (Hutchins 2012). On the one hand, the optimal tripping in numerical simulations, unlike real-life experiments, takes advantage of the synthetic nature to minimize the spatial length for developing to a fully equilibrium turbulence boundary layer. On the other hand, the optimised tripping function is based on a ‘collective history’ through the tuning procedure instead of a true development history from which the flow is physically transitioned. Therefore, the artificial tripping may not easily correspond to specific set-ups, and thus lack general transferability and comparability with experimental studies. In addition, there is little chance for the optimal tripping to be linked to many practical situations where the corresponding boundary layer flows cannot be idealised as a fully developed boundary layer in an equilibrium state. In fact, an under-developed non-equilibrium TBL may well affect a considerable part of the flow field with a non-negligible impact on aerodynamics performance (e.g. Wheeler, Dickens & Miller 2018). It is recognised that generating effective artificial inflow conditions is needed for studying a fully developed TBL at a high Reynolds number in a truncated domain. The emphasis in the present work, however, is on exploring an easily implementable physical tripping method to cover also a non-equilibrium part of the TBL leading to a fully developed equilibrium one with reasonable controllability and comparability to physical experiments.

It is noted that in many experimental studies and practical applications, a forced bypass transition has been widely explored and relevant experience has been established based on carefully designed trips for the TBL development. The influence of the different tripping devices including wires, grits, and pins on the downstream turbulent regime is scrutinised by Erm & Joubert (1991). Some rather simple tripping devices, for instance two-dimensional (2-D) simple step roughness elements, have proved to be quite successful (Klebanoff & Tidstrom 1972). It is worth noticing that in contrast to many previous experimental studies, very few efforts (e.g. Boudet, Monier & Gao 2015) have attempted to resolve properly physical trips in numerical simulations for TBL inflow generation. For wall-resolved LES, it is recognised that the mesh resolution in the near-wall region should approach that of a DNS (Moin 1997; Jimenez & Moser 2000). In this context, one wonders if and how much extra computational resource would be required to fully resolve a local tripping element. We would thus like to, first of all, examine how fine the local mesh should be to properly resolve a tripping element, particularly in comparison to a typical mesh required for a scale-resolving turbulent boundary layer solution at a similar Reynolds number.

Another relevant aspect of inflow conditioning is the physical behaviour and impact of those tripping-induced large-scale disturbances. Related issues are where in a boundary layer these large-scale disturbances reside, how they evolve streamwise and if they influence the near-wall turbulent flow development. These issues are particularly of interest in the light of the ‘foot-printing’ of large-scale coherent structures in outer flow (log-law region) on the near-wall region for well-developed wall-bounded turbulent flows. An overly disturbed boundary layer with a trip of large size can lead to rather persistent large structures in an outer flow region much further downstream of a TBL, as experimentally observed by Marusic et al. (2015). The sizing of the tripping element thus deserves attention. In general, it would be of interest to see how those tripping-induced large-scale disturbances behave in comparison to those outer flow coherent structures in a well-developed TBL. More specifically, would the tripping-induced large-scale disturbances ‘footprint’ on the corresponding near-wall turbulence?

1.3. Two-scale method with locally embedded fine-mesh blocks

The long-standing wisdom of the autonomous behaviour of the ‘universal’ inner layer has promoted an interest to explore the possibility to reduce the mesh count for scale-resolving turbulent flow simulations by adopting local fine-mesh blocks. Previous efforts include those placing a small locally truncated domain in the near-wall region, often called the ‘minimum flow unit’ (MFU), as adopted by Pascarelli, Piomelli & Candler (2000), Tang & Akhavan (2016), Sandham, Johnstone & Jacobs (2017), Carney, Engquist & Moser (2020).

It must be pointed out however that all previous MFU methods are strictly limited by the spatial periodic condition applied for the MFU domain in the two wall-parallel directions. The periodicity of the MFU unit length is profoundly incompatible with the prevalent ‘foot-printing’ of large scales from the outer flow on the near-wall flow. The large-scale structures are inevitably truncated by the domain lengths of the MFU. In relation to the common limit of the previous MFU methods, the distinctive scale-dependent interface treatment in the original two-scale framework (He 2018) becomes particularly relevant. As illustrated by Chen & He (2022) for a canonical channel flow, the scale-dependent interface treatment enables the large scales to directly pass through the interface, so that the ‘foot-printing’ can now be captured in the local near-wall fine-mesh block. The well-resolved local fine-mesh domain subject to well-captured ‘footprints’ then provides a much more suitable base on which the source terms can be generated to correct the global under-resolved coarse-mesh domain.

The locally embedded fine-mesh solution, as a key part of the two-scale framework, should naturally lend itself to a tripped turbulent boundary layer. In terms of the locality of the fine-mesh blocks in relation to the coarse-mesh domain, considerations need to be given to the two parts of the domain: first for the tripping itself and second for the tripped turbulent boundary layer. The roughness-induced bypass transition results from the forced generation of flow instabilities (Wu, Christensen & Pantano 2019; Kadivar, Tormey & McGranaghan 2021). The process involves the formation of vortical structures which will eventually break down and evolve into hairpin-type smaller-scale vortices in the turbulent regime (Brinkerhoff & Yaras 2011). The early stage of instability involves the growth of flow separation-associated vortical structures over the tripping element. Therefore, it is necessary to locally resolve the associated flow instabilities at different frequencies. A local fine-mesh block can be conveniently embedded around the tripping element to cover the transitional part, while the bulk flow can be efficiently resolved on a coarser base mesh. For a 2-D tripping element, the embedded fine-mesh block may potentially be truncated in the spanwise direction. The sensitivity of the predicted tripping to the span size of the fine-mesh block should thus be examined with the present two-scale method.

In contrast to a tripped boundary layer, for a canonical channel flow, the local solution can be directly mapped to the global coarse-mesh due to the global homogeneity in the two wall-parallel directions for the time-averaged flow, as illustrated previously by Chen & He (2022). Now for a tripped 2-D TBL, the global homogeneity still lies in the spanwise direction, while in the streamwise direction, only a ‘local’ homogeneity exists due to the smooth variation of the time-averaged flow. The local source terms could be analogously generated from the space–time averaged fine-mesh solution. These source terms are amendable to a spectral propagation to accommodate the global streamwise inhomogeneity. Chebyshev and Fourier spectral methods have been widely used in solving differential equations and many other applications of data interpolations (Boyd 2000). The block-spectral (BS) method is applied to inhomogeneous micro-structured surfaces (He 2018; Kapsis et al. 2020), based on Fourier spectral mapping. The Chebyshev method built on the polynomial group is preferred in the present work for the source-term spectral mapping without requiring any special periodicity.

1.4. Overall work scope and structure of the paper

As an extension of Chen & He (2022), the overall objective of the present work is to develop, validate and demonstrate a new approach to simulating TBL with equivalent accuracy to wall-resolved LES but at a much lower computational cost, manifested in terms of the mesh count–$Re$ scaling. The rest of the paper is organised as follows. First, the tripping to obtain properly initiated TBL based on a simple 2-D roughness element is introduced and analysed in § 2. Then, the implementation of the two-scale method for the evolving turbulent boundary layer flow is described in § 3. The section also includes the validation case studies and the mesh count–$Re$ scaling for the two-scale method, compared to fully wall-resolved LES and DNS. Finally, the summary and concluding remarks will be presented.

2. Physical tripping for scale-resolving turbulent boundary layer simulations

Fully wall-resolved LES (WRLES) calculations on the physically tripped TBL are studied in this section. The numerical method and case set-up are first introduced, followed by the full LES validation against DNS databases (Wu & Moin 2009; Schlatter & Orlu 2010, 2012) and the well-established experimental correlations (Smits, Matheson & Joubert 1983; Monkewitz et al. 2007; Chauhan, Monkewitz & Nagib 2009), and the investigation on the local resolution requirement for tripping. Thereafter, the streamwise development of the TBL towards an equilibrium state and the impact of residual large-scale disturbances of the tripping are examined.

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 2009; Lee & Moser 2015), 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 (1959), 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. (2013), 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 2011) 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 (1997) and Jimenez & Moser (2000), 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 (2009). The indicators of the pressure gradient magnitude of Wu & Moin (2009) 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 (2009) and an order of magnitude lower than those of Spalart & Watmuff (1993) 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. (1983). 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. (1983). (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 1997; Jimenez & Moser 2000). 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 (2010), 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. (2009).

Figure 3. Indicative parameters $\delta /{\delta ^\ast }$ and $\delta /\theta $ compared with the correlated curve of Chauhan et al. (2009) (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. (2007) 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 2010).

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. (2007). 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 (2010). 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 (2012) as solid black lines; Wu & Moin (2009) 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 1988; Jeong & Hussain 1995) 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. (2014). The roll-up of vortices becomes spanwise non-uniform and breaks down very soon downstream, similar to that observed by Brinkerhoff & Yaras (2011). 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 2007) 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 2013), 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 (2012) 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. (2009) (solid grey lines). The results of the tuned ‘optimal’ tripping of Schlatter & Orlu (2012) (‘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. 2009) 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. 2018) 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. (2018) 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. 2018). 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 2015), 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. (2015).

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 2015). 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 1967) 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 2022).

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. 1998) 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 (2022) 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 2022). (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 (2010).

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 2022), 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 2010). 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. 2017) 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 2014; Tissot et al. 2014). 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 (2016), 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 2022). 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. Two-scale block-spectral solutions

Having initiated a TBL with the physical tripping as described in the last section, we now turn to solving the tripped TBL efficiently and accurately. To this end, we resort to the two-scale BS method with locally embedded fine-mesh DNS blocks in the near-wall region. In this section, the basics of the two-scale methodology are first introduced briefly, including the dual meshing, the two-scale formulation together with the source-term propagation method. Then the spanwise sizing of the embedded fine-mesh block for the tripping will be examined. Thereafter, the two-scale BS methodology is applied to a boundary layer flow. The validity of the two-scale method for a tripped boundary layer is assessed by comparing the present solutions with the corresponding DNS and full LES results. Finally, the mesh count–$Re$ scaling for the two-scale TBL solutions will be estimated and compared to those of DNS and fully wall-resolved LES.

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 2011). 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 2018) 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 2018):

(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 (2018, 2021) and Chen & He (2022).

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 (2018) 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 2000). 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 (2022) 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 2015) 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. (2013). 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 (2022).

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 (2022). 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 1974). 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 2012).

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 (2012) (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 (2022). 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. 1983) and the DNS data (Schlatter & Orlu 2010).

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. (1983). 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 2010) (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 (2012), 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 (2022) 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. 2013). 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 (2012) 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. (2014). 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).

4. Summary and conclusions

The principal objective of the present work is to explore an efficient and accurate methodology for scale-resolving simulations of turbulent boundary layer flows. Previously, Chen & He (2022) demonstrated a two-scale method for the canonical channel flow where the mesh-count scaling with Reynolds number is potentially reduced from $O(R{e^2})$ of the full WRLES to $O(R{e^1})$. The present work extends the methodology to spatially evolving boundary layers. Two main issues of interest are, first, how to start a properly initiated TBL, and second, how to implement the local embedded two-scale method for a streamwise inhomogeneous TBL.

We resort to a physical tripping method as widely adopted experimentally to initiate a TBL. For a simple 2-D step element to trip the TBL, extensive computational analyses are carried out, first to identify the adequate mesh resolution around the trip and second to examine the evolution and impact of the large-scale disturbances associated with the tripping. We validate the full LES with $R{e_\theta }$ being approximately 1000 at the exit, of which the results match well with well-established DNS results (Wu & Moin 2009; Schlatter & Orlu 2010, 2012) and the experimental correlation (Smits et al. 1983; Monkewitz et al. 2007; Chauhan et al. 2009). It is found that a fine DNS mesh $(\Delta {x^ + } \approx 10,\;\Delta {z^ + } \approx 7)$ should be sufficiently fine for resolving the trip without any extra local mesh refinement. A coarser mesh $(\Delta {x^ + } \approx 20,\;\Delta {z^ + } \approx 15)$ would require local refinement for an adequate resolution, corresponding to a minimal 2 % extra mesh count. The mesh sensitivities with and without local refinement indicate that the mesh resolution required for resolving a trip properly should satisfy a dual requirement of having both a fine enough first wall-normal mesh spacing well within the local viscous sub-layer (e.g. $\Delta {n^ + } < 2$) on the trip element surface, and having a fine enough mesh for the inner boundary layer region around the trip, approaching a DNS mesh for the near-wall region.

Particular attention is paid to the post-tripping non-equilibrium development phase of a turbulent boundary layer before reaching a fully developed equilibrium state. Excessive large-scale disturbances residing in the lower log-region are associated with the tripping, clearly observable from spectral and modal representations. These triping-induced large scales in the outer flow are shown to survive a longer distance downstream compared with the inner near-wall counterpart. Both modal analyses using the EMD (Huang et al. 1998) and the DMD (Schmid 2010) indicate that initial tripping-induced large-scale disturbances do not have noticeable ‘footprints’ on the near-wall flow. Moreover, the DMD modal analyses illustrate clearly that the streamwise decay of the trip-induced disturbances is scale-dependent, with larger scale disturbances decaying seemingly faster than shorter ones. Overall, not only do these outer-flow large scales have no marked ‘footprinting’ on the near-wall region, but they also seem to be more difficult to survive along the streamwise direction than those with shorter-length scales. This is in clear contrast to well-established wall-bounded turbulence at a high $Re$ where long coherent streaks appear persistently in the outer flow region with a marked ‘footprinting’ on the near-wall part (Marusic, Mathis & Hutchins 2010).

When implementing the two-scale method for the tripped boundary layer flow, a locally embedded fine-mesh block is first set up and sized for adequately covering and resolving the trip element and its downstream proximity where the tripping-related large-scale vortical disturbances are initiated. The tripping effects on the downstream turbulent regime are seemingly insensitive to the spanwise size of the frontal fine-mesh block. Characteristics of large-scale unsteadiness around the trip appear to be similar regardless of the spanwise size of the embedded block though detailed spectra over the trip show some differences. These differences manifested in the local spectra around the trip seem to vanish quickly further downstream with very limited history effects. Therefore, the frontal fine-mesh block with a short span comparable to those embedded blocks for the tripped TBL is deemed adequate.

For a tripped turbulent boundary layer, the block-spectral method with a Chebyshev mapping is shown to be effective. The space–time averaged fine-mesh solution is taken as a target, leading to source terms for the upscaled equations in the corresponding coarse-mesh region. The source terms generated from the locally embedded near-wall fine-mesh blocks are effectively propagated to the global coarse-mesh near-wall region by the Chebyshev spectral mapping. The present two-scale method and implementation are validated and demonstrated by the calculated mean statistics and energy spectra in good agreement with full LES and DNS results. Finally, the scaling of computational mesh count with respect to the Reynolds number (based on the boundary layer thickness $\delta $) is estimated. It is highlighted that the overall mesh count–$Re$ scaling may be potentially reduced from $O(Re_\delta ^{1.8})$ for the full wall-resolved LES to $O(Re_\delta ^{0.9})$ for the present two-scale method.