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 Reference Moin1997; Jimenez & Moser Reference Jimenez and Moser2000). 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 Reference Jimenez2003). 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 Reference Mizuno and Jiménez2013), 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. Reference Spalart, Jou, Strelets and Allmaras1997; Cabot & Moin Reference Cabot and Moin2000; Piomelli & Balaras Reference Piomelli and Balaras2002; Spalart Reference Spalart2009; Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2015; Bose & Park Reference Bose and Park2018). 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 Reference Hutchins and Marusic2007; Jimenez Reference Jimenez2013; Lee & Moser Reference Lee and Moser2015).
A distinctive alternative approach the authors adopt is the two-scale methodology (He Reference He2018; Chen & He Reference Chen and He2022). 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 (Reference He2018), 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 (Reference Chen and He2022) 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 Reference Sillero, Jimenez and Moser2013) is to recycle and rescale the downstream characteristics to generate synthetic turbulent inflow conditions (Lund, Wu & Squires Reference Lund, Wu and Squires1998). 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 Reference Schlatter and Orlu2010). Other synthetic inflow generation methods include those built on theoretical properties of a well-developed wall-bounded turbulent flow (Kraichnan Reference Kraichnan1970) and empirical information (Jarrin et al. Reference Jarrin, Benhamadouche, Laurence and Prosser2006; Adler et al. Reference Adler, Gonzalez, Stack and Gaitonde2018). 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 (Reference Schlatter and Orlu2012) 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 Reference Hutchins2012). 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 Reference Wheeler, Dickens and Miller2018). 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 (Reference Erm and Joubert1991). Some rather simple tripping devices, for instance two-dimensional (2-D) simple step roughness elements, have proved to be quite successful (Klebanoff & Tidstrom Reference Klebanoff and Tidstrom1972). It is worth noticing that in contrast to many previous experimental studies, very few efforts (e.g. Boudet, Monier & Gao Reference Boudet, Monier and Gao2015) 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 Reference Moin1997; Jimenez & Moser Reference Jimenez and Moser2000). 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. (Reference Marusic, Chauhan, Kulandaivelu and Hutchins2015). 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 (Reference Pascarelli, Piomelli and Candler2000), Tang & Akhavan (Reference Tang and Akhavan2016), Sandham, Johnstone & Jacobs (Reference Sandham, Johnstone and Jacobs2017), Carney, Engquist & Moser (Reference Carney, Engquist and Moser2020).
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 Reference He2018) becomes particularly relevant. As illustrated by Chen & He (Reference Chen and He2022) 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 Reference Wu, Christensen and Pantano2019; Kadivar, Tormey & McGranaghan Reference Kadivar, Tormey and McGranaghan2021). 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 Reference Brinkerhoff and Yaras2011). 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 (Reference Chen and He2022). 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 Reference Boyd2000). The block-spectral (BS) method is applied to inhomogeneous micro-structured surfaces (He Reference He2018; Kapsis et al. Reference Kapsis, He, Li, Valero, Wells, Krishnababu, Gupta, Kapat and Schaenzer2020), 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 (Reference Chen and He2022), 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 Reference Wu and Moin2009; Schlatter & Orlu Reference Schlatter and Orlu2010, Reference Schlatter and Orlu2012) and the well-established experimental correlations (Smits, Matheson & Joubert Reference Smits, Matheson and Joubert1983; Monkewitz et al. Reference Monkewitz2007; Chauhan, Monkewitz & Nagib Reference Chauhan, Monkewitz and Nagib2009), 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
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.
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.
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.
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.
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 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 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$.
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 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.
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
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$).
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).
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).
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.
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
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:
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
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$,
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.
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.
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.
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. 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.
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:
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
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}$:
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:
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
The source term $\boldsymbol{S}{\boldsymbol{T}_{\boldsymbol{st}}}(\boldsymbol{x})$ is time-invariant and consists of two parts (He Reference He2018):
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:
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:
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
The sampling points on the transformed interval [−1, 1] are based on the Chebyshev–Gaus–Lobatto points as
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
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
where the constructor matrix $\boldsymbol{\varLambda }$ is a linear system of $N\; + \; 1$ dimensions (Trefethen Reference Trefethen2000). The constructor matrix can be expressed as
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.
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.
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.
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.
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.
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).
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.
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$.
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.
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.
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).
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.
The total mesh count required for the whole TBL $({x_0} \le x \le {L_x})$ is
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
The three terms on the right-hand side of (3.15) corresponding to the three domains (figure 31b) are estimated as
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
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
As ${f_i}$ for subdomain i is only a function of the local multiplier ${C_i}$ independent of Reynolds number, we will then have
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.
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 (Reference Chen and He2022) 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 Reference Wu and Moin2009; Schlatter & Orlu Reference Schlatter and Orlu2010, Reference Schlatter and Orlu2012) and the experimental correlation (Smits et al. Reference Smits, Matheson and Joubert1983; Monkewitz et al. Reference Monkewitz2007; Chauhan et al. Reference Chauhan, Monkewitz and Nagib2009). 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. Reference Huang1998) and the DMD (Schmid Reference Schmid2010) 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 Reference Marusic, Mathis and Hutchins2010).
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.
Acknowledgements
The authors would also like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC).
Funding
The work is supported in part by the EU Horizon 2020 grant, Marie Skłodowska-Curie Actions (Control of Turbulence Fricion Force, CTFF). C. Chen is on a scholarship from the China Scholarship Council (CSC).
Declaration of interests
The authors report no conflicts of interest.
Author contributions
C. Chen carried out the method implementation, case study set-up, simulations and data-processing. L. He instigated the framework methodology and case study design. Both authors contributed equally to analysing and interpreting results, and in writing the paper.