Turbulent–turbulent transient concept in pulsating flows

Abstract The turbulence behaviour of current-dominated pulsating flows has been investigated. Direct numerical simulations have been carried out for Stokes lengths over a range of $l_s^+=5\unicode{x2013}26$, and amplitudes spanning 90 % of the current-dominated regime, about a mean flow of $\overline {Re}=6275$. The results show that the turbulence response in intermediate and low-frequency pulsations is governed by a multistage turbulent–turbulent transition process, which bears a strong similarity to the multistage response of non-periodic acceleration. During the early acceleration period, the flow enters a pretransition stage, in which a new laminar perturbation boundary layer forms at the wall, and the streamwise velocity streaks are stretched. If the low-speed streaks destabilise prior to the deceleration period, then the flow enters a transition stage in which the perturbation boundary layer undergoes a bypass-like transition process. A unique feature of pulsating flows is the ongoing mechanism of turbulence decay, which initiates during the deceleration period and constitutes the main transient turbulence mechanism for much of the cycle. For high-frequency pulsations, the perturbation boundary layer fails to reach the pretransition stage prior to the deceleration period. Instead, the flow alternates between two inertial stages which are characterised by two layers of amplified viscous force; one growing at the wall, and one detached and moving towards the core.


Introduction
Unsteady flows arise throughout many engineering systems, as well as the natural environment (e.g.nuclear power plants, hydraulic turbines, blood flow in large arteries and sediment transport under sea waves).Unsteady flows can be categorised as non-periodic and periodic flows.Non-periodic unsteady flows undergo a temporal acceleration or deceleration from one bulk velocity to another.Periodic unsteady flows are subject to a continuous velocity change represented by a sinusoidal temporal variation.Periodic flow can be further subdivided into oscillating flow, for which the mean bulk velocity of the base flow is zero, and pulsating flow, for which the bulk velocity of the base flow is non-zero.In both cases, the flow oscillates with a frequency of ω, and an amplitude of either A b or A uc , which denote the amplitudes of the bulk velocity and centreline velocity, respectively.In pulsating flow, A b and A uc are scaled by the corresponding time-averaged mean velocity values.

Turbulent-turbulent transient flows
A key feature governing accelerating flows is the delay in the response of the turbulence flow field following an increase of the bulk flow rate.This delay was first observed in the experimental study of Maruyama, Kuribayashi & Mizushina (1976), which identified that the generation of new turbulence at the wall, and its propagation into the core of the flow, was the dominant mechanism of the turbulence response.He & Jackson (2000) found that the propagation of new turbulence into the core of the flow is preceded by two delays in the turbulence response.Firstly, there is a delay in the production of new turbulence at the wall, such that turbulence initially remains 'frozen' throughout the domain.Secondly, there is a further delay between the production of new turbulent energy in the near-wall region, and the transfer of that energy to the wall-normal and spanwise velocity components.The later experiments by Greenblatt & Moss (2004) found that propagation of new turbulence into the outer region eventually leads to the reignition of turbulence production away from the wall, as a secondary mechanism of turbulence growth in this region.
A crucial breakthrough was made when He & Seddighi (2013) utilised high-resolution direct numerical simulation to visualise and quantify the evolution of turbulence structures during step-up acceleration in channel flow.They identified that the evolution of the turbulence structures closely mimicked a bypass transition process, which occurs when laminar-turbulent transition in a developing laminar boundary layer is triggered by instabilities originating from turbulence in the free stream (Jacobs & Durbin 2001;Wu & Moin 2009).They concluded that the generation of new turbulence following a step-up acceleration originates from bypass transition in a newly formed laminar boundary layer at the wall, and that the delay in the turbulent response is due to the time taken for the boundary layer to destabilise.From these observations, He & Seddighi (2013) developed the turbulent-turbulent transient concept for accelerating flows in which the flow moves through three stages of evolution.The pretransition stage (i) is characterised by the growth of a laminar perturbation flow, advection of a frozen turbulence field and elongation of the existing streamwise velocity streaks at the wall.The transition stage (ii) is characterised by the destabilisation of the low-speed velocity streaks, and an explosion of turbulence in the perturbation boundary layer.The fully turbulent stage (iii) is characterised the propagation of new turbulence into the core of the flow.
The multistage process of the turbulent-turbulent transient concept has since been confirmed for temporally accelerating pipe flows (He, Seddighi & He 2016;Guerrero, Lambert & Chin 2021), ramp-up, spatially accelerating channel flows (Falcone & He 2022), and ramp-up, temporally accelerating channel flow (Seddighi et al. 2014).The more gradual, ramp-up acceleration rates resulted in an increase in the time period of the pretransition stage, and a delay in the onset of transition.The experiments of Mathur et al. (2018a) utilised two-component particle image velocimetry and flow visualisation, to confirm the growth of a perturbation boundary layer, and the dominant bypass transition process, which was only previously observed in numerical studies.Further numerical studies have explored the influence of Reynolds number and acceleration rate on the turbulence response.He & Seddighi (2015) observed that the pretransition stage of step-up accelerating channel flows was present even in weakly accelerated flow.However, as the strength of the acceleration was decreased, the strength of the streamwise streaks in the pretransition stage, and the strength of the subsequent turbulent spots, would diminish until the typical bypass transition behaviour was no longer observable.After observing similar trends in ramp-up channel flow, Jung & Kim (2017) suggested that bypass transition may only occur if the additional impulse from the increased flow rate greatly exceeds the additional impulse from the increased shear stress.On the other end of the scale, Mathur, Seddighi & He (2018b) observed that increasing the ratio between the final and initial Reynolds numbers amplified the length of the streaks during the pretransition stage, and increased the strength of the turbulent spots at the onset of transition.He & Seddighi (2015) found that, despite the differences in the transition process, all accelerating boundary layers showed a strong similarity in their perturbation velocity fields.During the pretransition stage, these profiles conformed to the theoretical Stokes solution for a step-up accelerated laminar flow, which confirmed the initial laminar behaviour of the perturbation boundary layer.The validity of the laminar Stokes solution has since been demonstrated for a wide range of accelerating flows with linear ramp-up acceleration rates (Sundstrom & Cervantes 2017) and arbitrary time-dependant acceleration rates (Mathur et al. 2018a).Sundstrom & Cervantes (2018b) further identified that, similar to accelerating flows, linear ramp-down decelerating flows produce a perturbation velocity field which corresponds to the laminar Stokes solutions for a period immediately following the onset of deceleration.
A recent study by Guerrero et al. (2021) identified an additional stage of the flow response which precedes the pretransition stage, referred to as the 'inertial' stage.During this stage, the viscous force (VF) in the near-wall region rapidly increases, starting in the viscous sublayer.These amplified forces serve as a momentum sink in order to preserve the no-slip condition at the wall.Then in the following 'pretransition' stage, the enhanced VF rapidly decreases, whilst the turbulence field remains frozen.Guerrero, Lambert & Chin (2023) found that the 'inertial' stage of a ramp-down decelerating flow is characterised by a rapid increase in the magnitude of VF, as in acceleration.However, the sign of VF changes during deceleration, to become positive, with VF serving as a source of momentum.

Periodic flows
Periodic flows are traditionally characterised by the formation of the Stokes boundary layer adjacent to the wall.The thickness of this layer, l s , represents the perturbation distance of the oscillating response from the wall.This distance is directly related to the pulsation frequency, ω, such that l s = √ 2ν/ω.Ramaprian & Tu (1983) showed that the oscillations impose an influence on the turbulence flow field up to a distance l t which lies far beyond the Stokes boundary layer.In the present study, l t is referred to as the turbulent penetration depth in accordance with He & Jackson (2009), though in previous studies it is also commonly referred to as the turbulent Stokes length/number.Ramaprian & Tu (1983) constructed a groundbreaking model for characterising pulsating flows, which decomposed the behaviour of the oscillating response into distinct regimes, each defined by the magnitude of the turbulent penetration depth in relation to the channel width or pipe diameter.Scotti & Piomelli (2001) later refined and clarified the limits and characteristic behaviour corresponding to each regime, and proposed a direct relationship between the inner-scaling of l s and l t , dependant solely on the von Kármán constant.
The high-frequency regime is defined as the frequency range for which the Stokes layer lies entirely within the viscous sublayer.In theory, although the turbulent penetration length extends into the buffer layer, the Stokes layer will be uninfluenced by the turbulence modifications occurring beyond the viscous sublayer, and the velocity field will exhibit quasilaminar behaviour.The experimental studies of Tardu, Binder & Blackwelder (1994) and He & Jackson (2009), which utilised fixed amplitudes of A uc = 0.64 and A uc = 0.2, respectively, both identified an upper limit of l + s < 10 for the high-frequency regime.It is important to note that many observations in current-dominated pulsating flows, including the high-frequency regime, are derived from studies in which the amplitude of the pulsation does not exceed A uc = 0.7 ∼ 0.8.The high-resolution numerical simulations of Manna, Vacca & Verzicco (2015, 2012) investigated the influence of amplitude on the turbulent response of wave-dominated pulsating flows.Their direct numerical simulations included one case with an amplitude of A uc = 1.0, and hence, lying on the limit of the current-dominated regime.Although their fixed Stokes length of l + s = 3.1 lay within the high-frequency regime, all Reynolds stress components showed significant variation far beyond the theoretical turbulent penetration depth.Early experimental studies identified an additional 'very high-frequency regime' (l + s < 7) in which the oscillating velocity field would have lower amplitudes than the theoretical quasilaminar prediction (Mao & Hanratty 1986;Tardu & Binder 1993).This was attributed to a resonance with turbulent ejections, as the driving frequency approached the frequency of the turbulent bursting process.The numerical investigations of Papadopoulos & Vouros (2016) confirmed that for l + s < 6.8, velocity field and second-order turbulent statistics have a negligible dependence on the amplitude of pulsation, for amplitudes up to A b = 0.63, but reported no significant interaction between the turbulent bursting process and the oscillations.For similar amplitudes of A b = 0.64, Cheng et al. (2020) proposed that such a regime does exist, but can be characterised by an independence of the Reynolds shear stress cospectra, and is bounded by its highest spectral frequency (l + s ≤ 2.4).The intermediate-frequency regime is defined when the turbulent penetration length extends beyond the buffer layer but does not reach close to the centreline; l + t < 0.5Re τ (Scotti & Piomelli 2001).The turbulence in the core of the flow, beyond the turbulent penetration length, remains effectively frozen whilst the velocity field oscillates as a 'plug' flow.The low-frequency regime is defined when the turbulent penetration length covers the flow domain (l + t O(Re t )).Tardu et al. (1994) found that such pulsations produced the strongest growth of the streamwise Reynolds stress within the turbulence in the buffer layer, which incidentally is the region which provides the greatest contribution to the generation of new turbulent energy.He & Jackson (2009) used two-component laser doppler anemometry to compare the streamwise and wall-normal Reynolds stress components of low-amplitude pulsating pipe flows in the intermediate regime.They observed that the response of the wall-normal Reynolds stress lagged behind the response of the streamwise Reynolds stress during the cycle.The magnitude of this delay was greatest in the buffer layer and decreased when moving away from the wall.This behaviour established a similarity with the multistage response observed in their previous study of accelerating flows (He & Jackson 2000), in which the delay was attributed to the time taken for newly generated turbulent energy to transfer between the velocity components, and its subsequent propagation into the core of the flow.Such a similarity was reinforced by Sundstrom & Cervantes (2018c), who compared experimental results for pulsating and ramp-up pipe flows with acceleration periods beginning at comparable values of Re τ , and confirmed that the skin friction and Reynolds stress components followed similar trends in their response to the acceleration.Scotti & Piomelli (2001) found that pulsations of A uc ≈ 0.7 in the intermediate and low-frequency regimes significantly modified the near-wall turbulent structures, which underwent multiple stages of evolution throughout the cycle.This phenomenon will be discussed further in § 1.3.

Transition in periodic flows
Numerous experimental and numerical studies have observed the phenomenon of laminar-turbulent transition in purely oscillating flows.When the flow alternates between laminar and turbulent states, the transition has been observed to originate from the growth and elongation of low-and high-speed velocity streaks (Vittori & Verzicco 1998;Costamagna, Vittori & Blondeaux 2003).The experiments of Akhavan, Kamm & Shapiro (1991) observed a suppression of turbulent energy production at the wall immediately following the start of acceleration.The onset of flow transition was accompanied by a sudden rapid growth of turbulent energy shortly before the end of the acceleration period.Ozdemir, Hsu & Balachandar (2014) performed direct numerical simulations to observe the receptivity of an initially laminar oscillating flow to initial perturbations, for a wide range of Reynolds numbers.At first, the perturbations induced the growth of two-dimensional spanwise vortical rollers during acceleration.For oscillating flows in the 'intermittently turbulent' regime (A uc √ ωv ≥ 600), these rollers would destabilise prior to the end of acceleration, giving way to three-dimensional instabilities.In such cases, this led to an explosive growth of turbulence which continued during the early deceleration period.Ebadi et al. (2019) explored the temporal evolution of the momentum balance in intermittently turbulent oscillating flows.Once the onset of transition occurred, the turbulent inertia (TI) grew rapidly, in conjunction with a rapid increase in the rate of energy transfer to the wall-normal and spanwise turbulent motions.This evolution would bear some similarity to momentum balance evolution identified by Guerrero et al. (2021) for non-periodic acceleration.Furthermore, the early decelerating regime of Ebadi et al. (2019) was characterised by a shift in VF from a momentum sink to a momentum source, which is also observed in ramp-down decelerating flows (Guerrero et al. 2023).
To date there is a significant lack of high-resolution numerical data for pulsating flows in the intermediate and low-frequency regimes.Where such data exist (Scotti & Piomelli 2001;Weng, Boij & Hanifi 2016), a very limited range of amplitudes is considered.For low-amplitude pulsation (A uc = 0.1), Weng et al. (2016) found that the flow remained in a fully turbulent state throughout the cycle across a range of regimes, for 7 ≤ l + s ≤ 44.7.From their relatively high amplitude (A uc ∼ 0.7) large-eddy simulations, Scotti & Piomelli (2001) were able to visually investigate the occurrence and behaviour of laminar-turbulent transition during the pulsation cycle.As the flow accelerated within the low-frequency regime (l + s = 35), new streamwise velocity streaks formed and continued to grow in the near-wall region.These elongated streaks eventually destabilised, leading to the formation of turbulent spots concentrated around the low-speed streaks, which is consistent with a bypass-transition process.Besides this early investigation, numerical studies which have observed similar phenomena are rare.This is in part due to a strong focus on the high-frequency regime.

Present study
Within the past decade, significant advancements have been made in understanding the evolution of non-periodic accelerating and decelerating flows, stemming from the discovery the turbulent-turbulent transient concept of He & Seddighi (2013).The experimental studies of He & Jackson (2009) and Sundstrom & Cervantes (2018c) have opened the door to a new way of thinking about unsteady flows, in which periodic and non-periodic applications are not truly dissimilar and disconnected.Instead, the full pulsation period may be reframed as a decomposition of individual non-periodic acceleration and deceleration periods which occur successively.
In the present study, the turbulent-turbulent transient concept is investigated for a wide range of frequencies and amplitudes for current-dominated pulsating channel flow.Three values of Stokes length are selected to represent the high-frequency (l + s = 5), intermediate-frequency (l + s = 16) and low-frequency regimes (l + s = 26), whilst an additional Stokes length of l + s = 10 represents the boundary between the high-and intermediate-frequency regimes.For each value of l + s , a set of three amplitudes are computed; A b = 0.1, A b = 0.5 and A b = 1.0, spanning 90 % of the current-dominated regime.The lowest amplitude has previously been explored by Weng et al. (2016) for a wide range of frequencies and is mainly included here to serve as a reference case.The highest amplitude of A b = 1.0 lies on the upper limit of the current-dominated regime, which has not been previously explored for the range of frequencies considered here.

Methodology
A series of pulsating flows were simulated using an in-house code CHAPSim (Seddighi 2011;He & Seddighi 2013;Wang & He 2015).Direct numerical simulations were performed to solve the momentum and continuity equations for an incompressible flow, (2.1) where x 1 , x 2 , x 3 = x, y, z and u 1 , u 2 , u 3 = u, v, w denote the coordinates and velocity components in the streamwise, wall-normal and spanwise directions, respectively.All length and velocity values are normalised using the channel half-height, δ, and the time-averaged bulk velocity, U b , respectively.The Reynolds number is defined as Re = U b δ/ν, or Re = U b δ/ν (time-averaged).A fully explicit, low-storage, third-order Runge-Kutta scheme is used for the temporal discretisation of the nonlinear and viscous terms.The Poisson equation is solved using a fast Fourier transform.An additional time-varying source term has been added to the streamwise Navier-Stokes equation to generate a periodic pulsation flow rate.The bulk velocity follows the waveform in the following (figure 1b), where ω = 2π/T denotes the driving frequency, and the amplitude, A b , is scaled by U b : Table 1 displays the configuration of the pulsating flow cases which were computed for the present study.All cases have a time-averaged, bulk Reynolds number of Re = 6275.The superscript '+' denotes inner-scaled variables based on the kinematic viscosity, ν, and a reference friction velocity, u τ .For the present study, the reference value of u τ is taken as the friction velocity of a smooth-wall, steady-state channel flow simulation with a bulk  Reynolds number of Re = 6275.The approximation for the turbulent penetration depth, l + t , is taken from Scotti & Piomelli (2001), where κ ≈ 0.41 denotes the von Kármán constant: (2.4) With the exception of case L26A10, the domain in each case has streamwise, wall-normal and spanwise dimensions of L x = 4πδ, L y = 2δ and L z = 2πδ, respectively 982 A20-7 (figure 1a), and a mesh distribution of N x × N y × N z = 512 × 256 × 512.For case L26A10, which combined the highest value of the pulsation amplitude A b = 1.0 and lowest driving frequency ω + = 0.0032, the elongation of the streamwise velocity streaks is significantly amplified during certain periods of the cycle, to the extent that the streamwise length of the domain had to be increased to L x = 8πδ.The number of cells in the streamwise direction was increased to N x = 1024.All other properties of the domain geometry and mesh distribution remain unchanged.Adequacy of the computational domain is justified by making sure that the two-point correlation of turbulent fluctuation decays to approximately zero within the domain half-length.The final mesh has streamwise and spanwise resolutions of x + = 8.83 and z + = 4.42, a wall-normal resolution of y + = 0.50 at the wall and a wall-normal resolution of y + = 5.51 at the channel centreline.Further details and discussion of the suitability of the spatial resolution and channel geometry are given in the Appendix.The time step was allowed to vary, and controlled by a Courant-Friedrichs-Lewy condition of CFL ≤ 1, and maximum limit of t + ≤ 0.124.

Data reduction
The time-dependant flow statistics are calculated using ensemble averaging.The pulsation cycle is divided into N samp = 32 evenly spaced phases, beginning at φ = 0. Multiple consecutive cycles, N t , are computed for each case.At each phase, the flow fields are spatially averaged in the streamwise and spanwise directions, and temporally averaged over multiple cycles.The method for calculating the phase-averaged statistics is (2.5) In order to minimise the size of the output flow field data, the fluctuating components of the flow statistics are not calculated or stored prior to the ensemble-averaging procedure.The fluctuating components are calculated from the phase-averaged values, (2.6) The time-averaged values of the flow statistics are calculated from a temporal-averaging of the phase-averaged values, (2.7) A total of 50 cycles are averaged for the highest frequency cases of l + s = 5, whilst a total of 30 cycles are averaged for all cases of l + s = 10 and l + s = 16.For l + s = 26, the number of cycles are varied based on the size of the domain.A total of 20 cycles are averaged for amplitudes of A b = 0.1 and A b = 0.5.For A b = 1.0,where the number of cells in the streamwise direction is doubled, the number of cycles is reduced to 10.

Validity and scaling
It is first necessary to address the choice in scaling parameters for the phase-averaged flow field.In non-periodic unsteady flows, the reference value for the friction velocity is typically taken from the initial flow field, u τ 0 .Sundstrom & Cervantes (2018c) demonstrated that pulsating flows of moderate amplitudes are directly comparable to Figure 2. Phase-variance of the skin friction coefficient for all cases in table 1.
non-periodic accelerating flows, by taking the value of u τ at φ = 0 as the reference value of the initial flow field.A similar approach would be desirable in the present study, however, such a reference value cannot be consistently applied across the wide range of flow configurations outlined in table 1. Figure 2 displays the phase-averaged skin friction coefficient, C f = τ w /0.5ρU b 2 .In addition to the wide variation in u τ at φ = 0, there is clear flow reversal at the wall for all values of l + s at the highest driving amplitude.Flow reversal also occurs in the high-frequency regime for a lower amplitude of A b = 0.5.In such cases, the phase-averaged friction velocity vanishes towards zero as u τ changes sign.Hence, in the present study, the reference value of u τ = 0.0572U b is taken from a steady-state smooth value channel flow at Re = 6275.
In order to assess the validity of the present method, the results are compared with two previous numerical studies: Weng et al. (2016) for pulsating channel flows of Stokes lengths l + s = 10 and l + s = 25.8, for a fixed amplitude of A uc = 0.1, where A uc denotes the relative amplitude of the channel centreline velocity; and Manna et al. (2012) for a high-frequency pulsating pipe flow of l + s = 3.1 and A uc = 1.0.Figure 3 displays the time-averaged Reynolds stress and streamwise velocity profiles for all cases in table 1.At low amplitudes, it would be expected that the pulsating motion would have a negligible influence on the first-and second-order velocity statistics.As seen in figure 3, the profiles for all four Reynolds stress components for A b = 0.1 collapse onto a single profile across the full width of the channel, which shows strong agreement with the results of Weng et al. (2016), and the steady-state channel flow.As the amplitude grows, the time-averaged profiles begin to diverge from the steady-state flow.The streamwise, wall-normal and spanwise Reynolds stress profiles of case L05A10 show good agreement with the results of Manna et al. (2012) up to y + = 30.
Finally, a discrete Fourier transform is applied to the phase-averaged velocity field to decompose the oscillating response into a number of sinusoidal waveforms, where the amplitude and phase-lead for a given Fourier mode of variable ψ are denoted by ũ and Φ u for the fundamental mode, and ũm and Φ um for the various harmonic modes (m = 1, 2, 3, . ..). Figure 4 displays the amplitude of the first harmonic mode, and its phase-lead relative to the centreline velocity, of the streamwise velocity profiles for all cases in table 1.For all cases of l + s = 5, the harmonic modes of the streamwise velocity are negligible throughout the domain.The phase-lead and amplitude of the dominant fundamental mode show no dependence on the driving amplitude, even at the upper limit of the current-dominated regime.At the limit of the high-frequency regime (l + s = 10) the fundamental mode remains overwhelmingly dominant, though its phase-lead and amplitude begin to show a weak dependence on A b in the near-wall region.In the intermediate frequency regime (l + s = 16) the phase-lead is strongly dependant on the driving amplitude and grows as A b is increased.At this frequency, the amplitudes of the harmonic modes relative to the fundamental mode become significant, and for A b = 1.0, the amplitude of the fundamental mode accounts for only 75 % of the sum of the amplitudes of the first five modes.For l + s = 26, the influence of A b on the phase-shift is much weaker compared with that of the intermediate-frequency regime.The two lowest amplitudes collapse towards a common profile and show strong agreement with the low amplitude case of Weng et al. (2016) close to the wall, although they exhibit a slightly lower value outside of the viscous sublayer.For the highest amplitude of A b = 1.0, the relative amplitude of the fundamental mode deviates significantly from the prior cases, with a significant reduction within y + < 30, and a significant increase above y + = 30.At this amplitude and frequency, the contribution of the fundamental mode accounts for less than 60 % of the first five modes.

Flow visualisations
This section explores the evolution of the near-wall turbulent structures during the pulsation cycle for three cases (L26A05, L26A10 and L10A10), to provide an overview of the flow behaviour and turbulent mechanisms in play.The turbulent vortical structures are visualised through isosurfaces of λ 2 , which denote the second eigenvalue of the symmetric component of the velocity gradient tensor.The near-wall velocity streaks are visualised through positive and negative values of the streamwise fluctuating velocity component, u = u − u .The mean velocity reference value of u is taken from the phase-averaged value at the corresponding wall-normal location at each point.Each visualisation comprises of eight flow fields which are evenly spaced by a phase interval of 0.25π.Two additional flow fields are included at φ = 0.875π and φ = 1.125π, in order to provide greater clarity to the transitional behaviour of case L26A10.
At the start of the acceleration period in case L26A05 (figure 5), the near-wall flow is still undergoing a process of turbulence decay, which was initiated during the preceding deceleration phase.As the flow accelerates, the remaining vortical structures continue to gradually degrade, whilst the streamwise velocity fluctuations are further suppressed.At φ = 0.25π, the rate of decay in the vortical structures has slowed substantially, leaving the remaining vortical structures, which remain loosely dispersed over the surface, to convect downstream.At the same time, the remaining weak streamwise velocity streaks begin to grow and elongate.As the flow passes the midpoint of the acceleration period, new vortical structures begin to form in clusters that are mostly concentrated around the low-speed velocity streaks.The clusters of vortical structures, which indicate the presence of highly concentrated turbulent spots, continue to grow in size until they begin to merge.The evolution of the vortical structures between φ = 0.25π and φ = 0.875π closely resembles the pretransition stage of non-periodic acceleration (He & Seddighi 2015).By comparing figure 5 to the ramp-up flow case of Seddighi et al. (2014) these similarities become even more striking.During the remainder of the acceleration period and the early deceleration period there is no significant change in the density or distribution of the vortical structures at the wall.However, by φ = 1.5π the density of the vortical structures is visibly reduced as turbulence begins to decay towards the weakly turbulent flow field seen at φ = 0.
As for case L26A05, the flow field at the start of acceleration in case L26A10 (figure 6) is still in a stage of turbulence decay.However, in case L26A10 the vortical structures continue to diminish until φ = 0.25π, when the flow in the near-wall region comes close to resembling a fully laminar state.The flow does not immediately enter a stage of velocity streak growth, and there is a delay in the formation of low-and high-speed velocity streaks.The flow remains frozen in a weakly turbulent (almost laminar-like) state until the first streamwise velocity streaks begin to form anew at φ = 0.375π.The flow then follows the same pattern seen in figure 5 and Seddighi et al. (2014), and thus begins a process of bypass-like transition.However, there are key differences in the evolution of the flow in figure 6 when compared with the lower amplitude case in figure 5. Firstly, and most notably, there is a delay in the formation of the turbulent spots around the low-speed streaks.These spots do not begin to form until φ = 0.75π, at which time the turbulent spots in case L26A05 were already beginning to merge.Due to this delay, there is insufficient time for the bypass-like transition process to complete, and the flow reaches the start of the deceleration period before the turbulent spots have started to merge.Secondly, the turbulent spots in case L26A10 are much more refined, and consist of a greater number of vortical structures which are packed together with greater density.Finally, the streamwise velocity streaks in case L26A10 undergo far greater elongation, with a maximum length that is approximately three times larger than that seen in case L26A05.Figure 6 shows that the flow in case L26A10 begins the deceleration period in the early stages of bypass-like transition.At this point, the developing turbulent spots are limited to a few discrete locations and vortical structures are absent over the majority of the surface.It is natural to expect that such an incomplete transition process will fail to reach a fully turbulent state under a continuously increasing deceleration rate.However, during the initial deceleration period the turbulent spots continue to grow and start to merge.By φ = 1.25π the vortical structures cover the majority of the surface.Throughout the second half of the deceleration period, these vortical structures diminish along with the streamwise velocity streaks, as the turbulence decays towards the weakly turbulent state at φ = 0.The persistence of the transition process can be explained by considering the response of an equivalent non-periodic flow.When a fully developed turbulent flow is subject to a ramp-down deceleration, there is a delay in the response of the turbulence in the early stages of deceleration (Guerrero et al. 2023).During this delay, the Reynolds shear stress only undergoes a very gradual change in the near-wall region, such that turbulence remains effectively frozen in its initial state.This delayed turbulence response can also be confirmed for case L26A05, as can be seen in figure 5, and as will be further discussed in more detail in § 3.3.These findings suggest that the concept of a frozen turbulence response can be expanded in the context of pulsating flow to include both the existing turbulence structures, and the temporally evolving mechanism of bypass-like transition.
The similarity between free-stream turbulence induced bypass transition and the behaviour observed in figures 5 and 6 can be confirmed through the assessment of turbulent kinetic energy variations during streak development.Figure 7 displays the evolution of the wall-normal maximum values of the turbulent kinetic energy and its three velocity components; u u + max , v v + max and w w + max , throughout the cycle.In the early stages of free-stream turbulence bypass transition, the elongation of the streaks is reflected in the linear growth of u u + max and turbulent kinetic energy (2 TKE + max ) with streamwise distance (Matsubara & Alfredsson 2001;Fransson, Matsubara & Alfredsson 2005), or time in the case of the pretransition stage of temporal acceleration (He & Seddighi 2013).Furthermore, since this early growth of u u + max is not caused by amplified turbulence generation, the maximum wall-normal and spanwise velocity motions should be relatively unaffected prior to the destabilisation of the streaks.From figures 5 and 6 it is clear that in cases L26A05 and L26A10 the initial growth of turbulent energy is attributed almost entirely to the amplification of the streamwise turbulent motions, whilst the gradual decay of v v + max and w w + max adversely impacts the turbulent energy growth.Hence, it can be concluded that the growth of u u + max is an artificial amplification which is attributed primarily to streak elongation, as characteristic of bypass transition.In figure 7(c,d) similar behaviour is observed for cases L16A05 and L16A10, such that together, these four cases represent the four instances of a bypass-like transition observed in the present study.Reducing the Stokes length further (as shown in figure 7e) extends the artificial growth of turbulent kinetic energy far into the deceleration period.The flow then deviates from bypass-like transition behaviour, since the decay of u u + max counteracts the natural growth of v v + max and w w + max as the streaks destabilise.When the Stokes length is reduced to l + s = 10, for an amplitude of A b = 1.0 (case L10A10, figure 8), the flow undergoes the early stages of bypass-like transition, but fails to reach a fully turbulent state before the end of the deceleration period.Therefore, the critical limit for which a current-dominated flow can undergo a full turbulent-turbulent transition process can be assumed to lie within the range of l + s = 10 to l + s = 16.At the start of the acceleration period, turbulence remains strong at the wall, with vortical structures covering the surface, though with a greater density around the elongated low-speed velocity streaks.At this point, the vortical structures and velocity streaks are convecting in the reverse direction, due to a reversal of the flow at the wall during the late deceleration period.As the flow accelerates, the turbulence starts to decay and the vortical structures are diminished, whilst the streamwise velocity streaks break up.New streamwise velocity streaks begin to form towards the end of the acceleration period, though by φ = π there is no sign of turbulent spots, and a significant amount of 'old' turbulence still lingers in the flow.At φ = 1.5π a clustering of new vortical structures can be observed around the low-speed velocity streaks, indicating the growth of turbulent spots consistent with bypass transition.Throughout the second half of the deceleration period these clusters grow very slowly.This growth continues even as the near-wall flow is reversed, and the early stages of merging can be observed.However, this is as far as the process gets, and growth of turbulent spots is heavily suppressed throughout the remainder of the deceleration period.

Streamwise velocity streaks
The properties of the near-wall velocity streaks can be quantified through a two-point correlation of the streamwise velocity field (3.1a,b).The streamwise correlation R 11x quantifies the length of the streaks, whilst the spanwise correlation R 11z quantifies the  spanwise spacing between the streaks, which is assumed to be equal to twice the spanwise distance at which R 11z reaches its minimum value, (R z11 ) min .The magnitude of (R z11 ) min indicates the strength of the streaks relative to the surrounding flow in the corresponding x-z plane (Mathur et al. 2018b;Falcone & He 2022), Case L26A05 starts the acceleration period with weak streaks with no clearly defined or uniform spacing.However, figure 9(a) confirms that that by φ = 0.625π, strong velocity streaks have formed, with a consistent spacing of 2z + = 160.As the streaks are stretched, they maintain this spacing, whilst growing in strength.Following their destabilisation and resultant breakup, the strength of the streaks degrades, and their spacing shrinks, before plateauing around 2z + = 80.This strength and spacing remains effectively frozen during early deceleration.After φ = 1.25π the streaks begin to decay and lose their form, until a clear minimum of R z11 is no longer visible, though they never fully dissipate.
During the initial acceleration period in case L26A10 (figure 9b), R 11z remains positive for all spacings of 2z + < 400 throughout the domain.The velocity streaks that emerged during the preceding turbulent-turbulent transition process have almost completely decayed, and by φ = 0 their remaining decay has a minimal impact on the weak turbulence field.At φ = 0.375π, new velocity streaks are beginning to take shape, which converge to a clearly defined spacing of 2z + ≈ 160 by φ = 0.5π (figure 9b).This approximate spacing is maintained throughout the remaining life cycle of the streaks, right up until their disintegration.The streaks elongate and grow in strength throughout the remainder of the acceleration period, even after φ = 0.75π, when the magnitude of the minimum value of R 11z falls as the streaks begin to destabilise.The strength of the streaks prior to transition is significantly greater than that in case L26A05.However, at φ = 1.25π,R z11 has once again become positive for 2z + < 400 throughout the domain, indicating a complete breakdown of the streaks by this point.This behaviour in R 11z bears a strong similarity to that observed in the step-up accelerating flows of Mathur et al. (2018b), in which the Reynolds number was increased by a large factor of 6.5.When this factor was Case L10A10 (figure 9c) starts off in a similar manner to case L26A05.During acceleration, the existing streaks become much weaker as the existing turbulence continues to decay following the preceding cycle.At φ = 0.875π a clearly defined spacing can be observed, as new streaks begin to develop.By the end of the acceleration period the spacing between the streaks has shrunk to 2z + ≈ 100, and the strength of the streaks has increased substantially.Throughout the first half of the deceleration period there is minimal change in both the spacing and the strength of the streaks (figure 9c).However, small sinuosities start to form in the low-speed streaks.Such sinusoidal deformations have been identified as the initial instability from which streak breakdown typically occurs in bypass transition (Schlatter et al. 2008).The development of the sinuosities is slow, relative to the phase of the cycle, such that the elongated velocity streaks remain mostly intact by φ = 1.5π.In terms of inner-scaled time, these sinuosities grow at comparable rates between different strokes lengths for A b = 1.0.Sinuosity development for l + s = 10 spans a period of t + = 78.Consider that the snapshots of case L26A10 in figure 10 are spaced by t + = 130.8 apart.After the appearance of the first sinuosity, it takes t + = 130.8 for the instability to develop into a local breakdown of the streak, and at least an additional t + = 130.8 for the turbulent spot to grow sufficiently to begin merging.Additional time will be required before new turbulence growth will propagate into the core (y + = 106.4),where the turbulence field has shown only minor changes between φ = 0.75π and φ = π.By comparison, time between the first sinuosity and the first turbulent spot in a flow of l + s = 16 is t + = 100, with a further t + = 75.0for merging to occur.Although these processes occur at different phases in the cycle, and hence, at different bulk velocities, the growth rate of the sinous instabilities appears to be unaffected by the local bulk velocity, at least in the late acceleration period and early deceleration period.

Turbulent statistics
During turbulent-turbulent transition in non-periodic acceleration, the temporal evolution of the Reynolds stress components follow a multistage process which correlates to the changes of near-wall turbulent structures (He & Seddighi 2015).These stresses are important in accurately defining both the temporal limits and characteristic behaviours of the distinct stages of the turbulent-turbulent transition process.To understand the flow of energy between the components, we also assess the temporal evolution of the terms in the turbulent energy budget for the streamwise Reynolds stress component (3.2); the production P 11 , viscous dissipation e 11 , pressure strain Π 11 , pressure diffusion Γ 11 , turbulent transport T 11 and viscous diffusion D 11 .The wall-normal and spanwise components do not gain energy directly from newly generated turbulence.Their primary source comes from the transfer of energy from the streamwise component, which is governed by the pressure strain.Hence, P 11 , which acts at the primary source of newly generated turbulent energy, and Π 11 , serve as the most critical budget terms which govern The phase-variance of the production, pressure strain, viscous dissipation, turbulent transport and viscous diffusion streamwise budget terms are displayed alongside the Reynolds stress components for cases L26A05 (figure 11), L26A10 (figure 12) and L16A10 (figure 13).The pressure-diffusion term Γ 11 is negligible for the streamwise component, and hence, is excluded from this analysis.The evolution of Reynolds stress components and streamwise budgets terms (3.2) in case L16A05 closely resembles that of case L16A10, and so the corresponding plots for the former case have been excluded for clarity.For these pulsating flows, we have decomposed the turbulent-turbulent transition cycle in five distinct stages: the (I) residual decay; (II) pretransition; (III) transition; (IV) post-transition; (V) turbulence decay stages.
In case L26A05 (figure 11), immediately following acceleration there is a short delay before the elongation of velocity streaks which signals the beginning of a pretransition stage.During this delay, v rms and − u v show minimal change at y + = 12.0 and y + = 5.4, whilst they gradually decay farther away from the wall.Meanwhile w rms shows a continuous decay throughout the domain.These same trends are also observed in varying degrees for other intermediate and low-frequency cases (figures 12 and 13).At first glance,  it is tempting to attribute this behaviour of decay to the growth of a new laminar boundary layer, and the displacement of existing turbulent structures away from the wall.However, this cannot be the case, since in non-periodic accelerating flows the perturbation flow field remains independent of the base flow, and exerts no significant change in v rms and w rms , prior to the onset of transition (Seddighi et al. 2014;He & Seddighi 2015;Guerrero et al. 2021).Instead, this response represents a continuation of ongoing turbulence decay which started within the preceding deceleration period, as will be further discussed later Figure 12.Phase-variance of the Reynolds stress components and streamwise energy budget terms for case L26A10.
in this section.We refer to the period of this delay as the 'residual decay' stage, which is unique to periodic acceleration in pulsating flows.The pretransition stage is marked by a gradual growth of u rms .As in the case of non-periodic accelerating flow, the growth of u rms stems primarily from the elongation of high-and low-speed streamwise velocity streaks and not from the generation of . Phase-variance of the Reynolds stress components and streamwise energy budget terms for case L16A10.
new turbulence.However, this elongation does produce a mild response in the streamwise production term P 11 .This occurs because the streamwise production term depends on the Reynolds shear stress, see (3.2), which in turn relies on the strength of both the streamwise and wall-normal turbulent motions.Since v rms remains effectively frozen during pretransition, the streak elongation causes production to grow gradually in line with u rms .Case L26A10 (figure 12) presents an exception to this behaviour, as the near-total decay of v rms close to the wall ensures that the Reynolds shear stress, and hence streamwise production, remains suppressed despite the growth of u rms .In contrast to non-periodic acceleration, the wall-normal and spanwise Reynolds stresses do not remain frozen during pretransition.Whilst v rms shows minimal change at the wall, its strength continues to diminish within the core, w rms falls continuously throughout the domain.Clearly, the ongoing process of turbulence decay, and its propagation into the core, from preceding stages (I and V), remains active.Since the pressure strain is heavily suppressed prior to destabilisation, this decay process initially constitutes the primary mechanism governing the evolution of v rms and w rms .
In the present study, the transition stage is the most sharply defined out of all stages in the periodic turbulent-turbulent transition process.Its beginning is marked by the rapid shift in the growth rate in the streamwise production term.This is followed by a rapid growth of Π 11 , which indicates that the critical turbulent mechanism of redistributing energy to the wall-normal and transverse motions becomes highly active.The elongation of velocity streaks is no longer the primary mechanism for extracting turbulent energy from the flow, and the evolution of P 11 now represents the amplification of true turbulent energy production.Comparison between figures 5 and 6 makes it clear that this new energy originates from destabilisation and breakup of low-speed velocity streaks, and the growth of highly concentrated turbulent spots.In all cases, this process causes a 'violent explosion' of the Reynolds stresses at the wall, at least in the case of v rms , w rms and − u v .So far, this behaviour closely mimics the turbulent-turbulent bypass transition of a perturbation boundary layer for non-periodic acceleration (He & Seddighi 2013).In each case, the transition stage is marked by an accelerated growth rate in u rms , spurred on by the enhanced production of new turbulent energy.However, at the end of transition, case L26A05 (figure 5) is the only case to display an overshoot of u rms in the buffer layer, as is characteristic of the transition in non-periodic acceleration (Jung & Kim 2017).The reason for this relates to the deceleration period and will be elaborated on further in the discussion of stage V.
In non-periodic acceleration, completion of transition is followed by a final stage of relaxation, in which near-wall turbulence settles into a state of equilibrium, whilst the core region waits to settle through the slow process of turbulent propagation from the wall.In all cases where the transition stage extends into the deceleration period (figure 12 and 13) there is insufficient time for such relaxation, as the growing deceleration rate quickly ensures that turbulence decay takes over as the main transient turbulence process.However, in case L26A05, figure 11 shows that if the transition stage ends whilst the flow is still accelerating, then an intermediate stage emerges.By φ = 0.875π the growth rates of P 11 and Π 11 are beginning to slow down.As a result, the growth rates of all Reynolds stress components in the near-wall region fall substantially, particularly in the case of u rms , v rms and w rms .These three components then settle onto a relatively flat peak around their maximums.Meanwhile the strength of Reynolds stress components away from the wall (y + = 106.4)continues to grow, as new turbulence at the wall continues to propagate into the core.Whilst this stage shares many characteristics with the 'core relaxation' stage of non-periodic acceleration (Guerrero et al. 2021), its lifespan is only a fraction of the time scales needed for propagation of enough turbulence to bring the core into equilibrium.By φ = 1.25π, the Reynolds stress components close to the wall have already begun to decay whilst those at y + = 106.4are continuing to grow.Furthermore, at y + = 5.4 and y + = 12.4 the production and viscous dissipation (figure 11e,g), which reached their peak at the end of the transition stage, are becoming gradually weaker throughout the latter half of stage IV.With these differences in mind, this stage of the pulsating cycle could be reclassified in more general terms as a 'post-transition' stage, in which non-periodic and periodic acceleration share a common characteristic of a continuous growth of Reynolds stress components in the core, and a partially shared characteristic regarding the reduced growth rate of these stresses in the near-wall region.Their only variation lies in the severity of growth rate reduction, which does not completely reduce to zero in periodic acceleration.
Over a majority of the deceleration period all Reynolds stress components show a continuous reduction in strength which slows down as the flow approaches its initial state at φ = 2π(0).In each case, the process begins at the wall, first with u rms , and then, after a short delay, continuing with v rms and w rms , which occurs simultaneously.Each response is precipitated by the reduction of its respective primary source term (P 11 and Π 11 ).Once again, there is a delayed response in the core, as turbulence generated during the preceding transition or post-transition stage is initially still propagating towards the core region.For case L26A05, turbulent shear stress begins to decay at φ = π, whilst for the remaining three cases its decay occurs simultaneously with that of v rms and w rms .During the pretransition stage u rms and − u v were shown to not be reliable indicators of true turbulence generation.With this in mind, we opt to define a starting point for the turbulence decay stage as the point at which v rms and w rms begin to decay.
It is important to note that, in the absence of a distinct post-transition stage (figures 12 and 13), any potential overshoot in u rms towards the end of transition is masked by the sudden dominance of the decay process in the turbulence decay stage.By comparison, this overshoot is clearly visible in figure 11, such that u rms peaks at the wall, before briefly settling towards a pseudoequilibrium value.This distinction is important to keep in mind, as the overshoot of u rms is a typical feature of turbulent-turbulent transition, and has been offered as a reliable criterion for its detection in non-periodic acceleration (Jung & Kim 2017).The present results confirm that this criterion can be valid for certain pulsating flows, but should be applied selectively based upon the phase within which the transition stage concludes.

Perturbation flow field
The turbulent-turbulent transient process is initiated by the growth and destabilisation of a temporally developing perturbation boundary layer at the wall.He & Seddighi (2015) showed that the perturbation boundary layer can be obtained from a decomposition of an unsteady flow into two flow fields: a perturbation flow and a base flow.The base flow is time-independent, and is usually taken as either the initial flow field in the case of non-periodic flow, or the time-averaged flow field in the case of periodic flow.The perturbation flow is time-dependant and represents the difference between the flow field at a point in time and the base flow.In practical terms, the perturbation flow represents the change in the flow field which results from the unsteady forcing.Although the present investigation concerns a pulsating flow, the perturbation flow is defined as though the flow was a non-periodic accelerating flow.The base flow is taken as the flow field at the point of minimum bulk velocity (φ = 0) for both the acceleration and deceleration periods of the pulsation cycle.The perturbation velocity for the present case is given as Two theoretical laminar flow scenarios are considered: Stokes first problem for a flow accelerating from rest; and Stokes second problem for an oscillating flow.At first glance, the former scenario may be assumed to be irrelevant in the context of pulsating flows.However, as shown, Stokes first problem provides a close representation of laminar perturbation development during the pretransition stage of certain low-frequency, pulsating flows.The solution to Stokes second problem for a laminar flow in a channel is given as (3.4) For an accelerating or decelerating flow with an arbitrary, time-varying acceleration rate, the solution to Stokes first problem is given by the extended Stokes solution, (Schlichting & Gersten 2017) where 'erf' represents the error function of the form erf(ψ) = (2/ √ π) ψ 0 e −γ 2 dγ .The cosh function in (3.4) can be solved through a series expansion, although the expansion is slow to converge for all but very low frequencies (ω → 0).For very high frequencies (ω → ∞) the cosh function can be expressed in the asymptotic formula cosh(ψ) → 0. This simplification produces the quasilaminar solution for Stokes second problem which is valid in the high-frequency regime of pulsating turbulent flows (Manna et al. 2012;Papadopoulos & Vouros 2016).In its current form the error function in (3.5) must be solved either through numerical integration or approximation.The error function can be eliminated by integrating (3.5), however, the resulting solution becomes highly unstable as ξ → t.Hence, in the present study, (3.5) is approximated through numerical integration using dξ = 3.125 × 10 −6 (2π/ω).
By taking t = 0 at φ = 0, the perturbation velocity in each Stokes solution can be expressed in terms of a perturbation function, ζ , as follows: Extending (3.5) to account for consecutive acceleration and deceleration periods, and rearranging (3.4), produces two self-similar perturbation functions for a pulsating flow, corresponding to the quasilaminar (3.7) regime of an oscillating flow, and the non-periodic accelerating/decelerating solution (3.8), referred to here as the extended laminar solution, where η s = y/l s , and g(y, χ) represents an integral function of the form (3.9) The perturbation velocity for all cases of l + s = 10, l + s = 16 and l + s = 26 is displayed in figure 14, along with the quasilaminar solution, u ∧ ql , and extended laminar solution, 982 A20-25 u ∧ ext , for each specified frequency.The corresponding velocity profiles for l + s = 5 are not shown.Within the high-frequency regime, l + s = 5, the perturbation velocity profile corresponds very closely with the quasilaminar solution throughout the majority of the cycle for all values of A b considered.This would be expected as for l + s < 10 the laminar Stokes boundary layer remains confined within the viscous sublayer (Scotti & Piomelli 2001;Papadopoulos & Vouros 2016).
Firstly, consider the relationship between the quasilaminar and extended laminar Stokes solutions.Initially, for all cases the quasilaminar solution strongly overshoots the extended laminar solution close to the wall.Given that the streamwise velocity in the near-wall region has a phase-lead relative to the centreline velocity, the value of u − u 0 will grow at a faster rate than u c − u c0 during the early acceleration phase.Sundstrom & Cervantes (2017) found that the small value of u c − u c0 at the beginning of acceleration would magnify errors in the calculation of u ∧ .From figure 14, it would appear that these small values also magnify the effect of the phase-lead in pulsating flows.As the flow accelerates, the contribution of the phase-lead as a percentage of u ∧ ql diminishes, and the overshoot of the quasilaminar solution gradually shrinks as its profile moves closer to that of the extended laminar solution.The two solutions maintain excellent agreement with each other throughout the remainder of the cycle, although, due to the phase-lead in the quasilaminar solution, a complete collapse is never achieved.As the flow returns to the beginning of the cycle, both solutions predict a reversal of the perturbation velocity field at the wall.The close similarity that is maintained by (3.7) and (3.8) during this drastic change in flow behaviour shows that the proposed similarity between periodic and non-periodic flows is supported by the underlying theoretical models.
Out of all values of l + s , the progression of the perturbation velocity for l + s = 26 (figure 14c) shows the strongest dependence on A b .At the smallest amplitude of A b = 0.1 the perturbation velocity does not collapse towards either Stokes solution.The perturbation velocity gradually increases in the near-wall region, until it starts to overshoot the extended laminar solution at φ = 0.5π.Similarly, the perturbation velocity starts to undershoot the extended laminar solution away from the wall, such that the profile resembles that of a fully turbulent perturbation flow throughout the remainder of the acceleration period and the first half of the deceleration period.Out of the three amplitudes considered, the progression of u ∧ for A b = 0.5 (case L26A05) shows the strongest resemblance to the progression of an accelerating flow undergoing turbulent-turbulent transition (He & Seddighi 2015).Following the initial overshoot, u ∧ converges towards the extended laminar solution close to the wall.By φ = 0.125π, u ∧ has achieved a complete collapse within the region of y + ≤ 10, and by φ = 0.25π this region has grown to y + ≤ 25.The laminar behaviour of u ∧ within the newly formed perturbation boundary layer remains virtually unchanged between φ = 0.25π and φ = 0.5π , i.e. the pretransition stage (II), which can be explained by the frozen state of turbulent energy generation confirmed in figure 11( f ).After φ = 0.5π, u ∧ begins to grow within the region of y + ≤ 25, such that it now overshoots the extended laminar solution at the wall.At the same time, the profile of u ∧ undershoots the extended laminar solution, such that the profile of u ∧ resembles that of a turbulent perturbation boundary layer (He & Seddighi 2015;Mathur et al. 2018a).The growth of u ∧ coincides with the destabilisation of the near-wall velocity streaks (figure 5), and the rapid growth of turbulent energy production at the wall (figure 11e,f ).By the time that the streamwise Reynolds stress overshoots its pseudoequilibrium value at φ = 0.875π (figure 11a), the profile of u ∧ has shifted to overlap with the u ∧ profile   of the case A b = 0.1.The u ∧ profile for A b = 0.5 maintains a strong agreement with that for A b = 0.1 throughout the post-transition stage (IV), as the remaining Reynolds stresses, v rms , w rms and − u v , and the turbulent energy production P 11 settle around their respective peaks (figure 11).For the highest driving amplitude of A b = 1.0 (case L26A10), the profile of u ∧ also converges towards the extended laminar solution, although, whilst the overshoot in u ∧ falls significantly during the pretransition stage, it does not achieve a complete near-wall collapse as seen for A b = 0.5.Compared with A b = 0.5, there are significant delays in the departure of u ∧ from the extended laminar solution.The reasons for this delay are evident when considering the near-negligible values of v rms and − u v at the wall prior to φ = 0.8125π (figure 12b,d).The breakdown of the near-wall velocity streaks in the perturbation boundary layer (figures 5 and 6) depends on existing turbulence within the channel to induce the initial perturbations.Hence, in case L26A10 the pretransition stage (II) of the cycle is prolonged, during which time the perturbation boundary layer remains in a near-laminar state.However, following the explosive growth of the turbulence production term P 11 , the u ∧ profile immediately follows the evolution seen for A b = 0.5; to diverge from the extended laminar solution and move towards the profile for A b = 0.1.The fact that the transition stage (III) extends into the decelerating period appears to not significantly disrupt the evolution of u ∧ during this process.Given that case A b = 0.1 shows no significant distortion of its boundary layer structures throughout the cycle (see § 3.5), its profile for u ∧ can be assumed to represent the reference turbulent solution for the perturbation velocity field at higher amplitudes.Hence, the transition stage of the turbulent-turbulent transition process for l + s = 26 is characterised by the transition of the perturbation velocity profile between the theoretical extended laminar solution in (3.8) to the numerical turbulent solution for a reference low amplitude flow.
For cases L16A05 and L16A10, the perturbation velocity initially overshoots the extended laminar solution, similar to cases L26SA05 and L26A10.As the flow accelerates, the profiles of u ∧ settle into the general form of the extended laminar solution, but maintain a persistent overshoot prior to the onset of transition.During this period, the production of new turbulent energy, as indicated by −Π 11 in figure 13( f ), remains relatively static (figure 13), indicating that the perturbation flow field remains in a laminar state.The rapid growth of −Π 11 during the transition stage, φ = 0.875π and φ = 0.9375 for cases L16A05 and L16A10, respectively, correlates to a sudden growth of u ∧ at the wall, as seen for l + s = 26.Similarly, above y + ≈ 30 u ∧ starts to rapidly decay.The relative changes in u ∧ following φ = 0.875π mimics the transition-based response seen for the lower frequency cases of l + s = 26.However, given the initial overshoot in u ∧ , comparing this response with the extended laminar solution does not serve as a reliable criterion for identifying the laminar, transitional and turbulent states of the perturbation boundary layer, as it does for l + s = 26.
Since the perturbation velocity demonstrates a self-similar distribution for pulsating flows, it follows that the skin friction coefficient for the perturbation flow field will also demonstrate self-similarity.He & Seddighi (2015) introduced a perturbation form of the skin friction coefficient for accelerating flows, C ∧ f , in which the perturbation shear stress is scaled by the initial friction velocity and the final perturbation bulk velocity.An equivalent definition of C ∧ f can be derived for pulsating flow, as follows, where 2A b represents the perturbation bulk velocity at the end of the acceleration period, and u τ is the reference friction velocity for a steady state flow at Re = 6275: (3.10) The solution for C ∧ f in Stokes second problem for ω → ∞ is obtained directly from the derivation of (3.7), as follows: The quasisteady solution of C ∧ f is taken from the second-order approximation in Tardu et al. (1994): (3.12) The derivation of the extended laminar solution for Stokes first problem, (3.5) and (3.9), is unstable and relies heavily on the precision of the numerical integration of g(y, χ), due to the vanishing of the denominator as ξ → χ .Therefore, the solution of C ∧ f is approximated by assuming a linear gradient of u ∧ between the static wall, and an arbitrary point close to the wall:  f (qs) .In all cases, there is a notable similarity between the trend of the quasilaminar and extended laminar solutions for C ∧ f in terms of both the phase of the profile, and the relative amplitudes between their peaks.The extended laminar solution underpredicts the quasilaminar solution throughout the cycle.However, the relative magnitude of the discrepancy between these two solutions at different points in the cycle remains consistent for varying values of the Stokes length.The magnitude of | grows as the flow accelerates, reaching a peak at the same point at which C ∧ f (ql) reaches its maximum.Beyond this point |C ∧ f (ql) − C ∧ f (ext) | shrinks and the profiles converge, such that there is good agreement between the two solutions for a majority of the deceleration period.
The high-frequency regime l + s = 5 (not shown) shows a strong agreement with the quasilaminar solution, with the amplitude exerting no influence on C ∧ f .This agreement is also observed for l + s = 10 (figure 15a-c), particularly at the highest amplitude of A b = 1.0; however, as the amplitude is reduced, C ∧ f starts to undershoot C ∧ f (ql) as C ∧ f approaches its peak.This undershoot at l + s ≈ 10 is well known, with previous studies consistently observing a ratio between the shear stress amplitude of the oscillating flow field to that of the high-frequency Stokes solution of less than 1 (Weng et al. 2016;Sundstrom & Cervantes 2018a).However, the observation that this undershoot exists independently of 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 0.5π 1.0π 1.5π 2.0π 0 the driving amplitude is not universal to the whole current-dominated regime, particularly towards its upper limit, as shown in figure 15(c).
For l + s = 16 and l + s = 26, the behaviour of C ∧ f during the cycle is heavily influenced by the pulsation amplitude.At the lowest amplitude of A b = 0.1 (figure 15d,g), C ∧ f follows a simple sinusoidal temporal profile, which approaches the quasisteady approximation as the frequency is reduced (l + s = 16 → 26).For cases L16A05, L16A10, L26A05 and L26A10, the behaviour of C ∧ f varies considerably between the pretransition (II), transition (III) and turbulence decay (V) stages identified in § 3.3.Firstly, it must be noted that in all four cases the value of C ∧ f does not undershoot the extended laminar solution at any point in the cycle.Furthermore, C ∧ f does not overshoot the quasilaminar solution prior to 982 A20-30 entering the transition stage (III).Hence, it would appear that in a pulsating flow exhibiting a laminar-like response, C ∧ f remains confined to a region which is bounded by these two laminar solutions.In cases L26A05 and L26A10, the acceleration of u ∧ at the wall during the transition stage is reflected in the growth rate of C ∧ f .During this stage, C ∧ f overshoots the quasilaminar solution, and rises towards the quasisteady approximation.This rapid growth of C ∧ f towards the quasisteady approximation is consistent with the observations of Sundstrom & Cervantes (2018c) in experiments of pulsating pipe flow under similar conditions to case L26A05.However, the present results further show that delaying the transition stage until late in the accelerating period, as seen in case L26A10 (figure 15i), does not seem to have a significant impact on the response of C ∧ f .The delayed growth of C ∧ f simply leads to the intersection with the quasisteady approximation taking place at a later point in the cycle.By comparing figures 15(h) and 15(i) with figures 11 and 12, respectively, it can be seen that the shift of C ∧ f is reflective of the rapid generation of new turbulent energy at the wall, which can persist beyond the end of the accelerating period.The transition stages of cases L16A05 and L16A10 (figure 15e,f ) are also characterised by a departure of C ∧ f from the bounded laminar region, and its convergence towards the quasisteady approximation.However, unlike cases L26A05 and L26A10, this departure is not the result of an accelerated growth of C ∧ f , but instead a reduction in its rate of decay.Despite the rapid growth of turbulent energy generation (figure 13f ), C ∧ f effectively plateaus.Eventually it reaches close proximity to the quasisteady approximation, after which point, the rate of decay in C ∧ f recovers, such that C ∧ f temporally follows the trend of the quasisteady approximation.Comparing figures 15( f ) and 13(d,f ) shows that the sudden change in C ∧ f correlates closely to the point at which the Reynolds stress reaches its peak, which is subsequently followed by the rapid reduction of turbulent energy generation.Whilst the evolution of C ∧ f during the transition stage may vary for l + s = 16 and l + s = 26, the underlying mechanisms governing its evolution are consistent.It can be concluded that, as for non-periodic accelerating flow (He & Seddighi 2015;Sundstrom & Cervantes 2017, 2018a), the evolution of C ∧ f provides a strong criterion for identifying the presence of turbulent-turbulent behaviour in pulsating flows.

Momentum balance and the inertial regime
The temporal rate of change in momentum is governed by three forces: (i) VF, represented by the wall-normal gradient of the local shear stress ∂ 2 u + /∂y + 2 ; (ii) TI, represented by the wall-normal gradient of Reynolds shear stress −∂ u v + /∂y + ; (iii) pressure force (PF), represented by the streamwise pressure gradient − ∂p/∂x /ρ.Together they form the mean momentum balance, A fully developed channel flow can be divided into four distinct layers, each of which is characterised by a specific distribution of the forces in (3.14) (Wei et al. 2005;Cheng et al. 2020).Inside the classical viscous sublayer (y + ≤ 5) there exists a secondary 'inner viscous balance' layer (Layer 1) within which VF and the PF are in balance, and dominate over TI.Within the 'stress-gradient balance' layer (Layer 2) VF and TI are balanced but with opposite signs such that they effectively counteract each other, whilst PF is close to zero.The 'meso viscous balance' layer (Layer 3) contains the point at which TI crosses the zero axis, such that TI transitions from a source term to a sink term.The final layer is the inertial balance layer (Layer 4) where TI and PF are balanced whilst VF vanishes to zero.Note that this interpretation represents a steady-state channel flow where ∂ u /∂t ≈ 0, whilst in the present case |∂ u /∂t| > 0 due to the driving force.
Figures 16 and 17 show, respectively, the temporal variation of VF and TI during the pulsation cycle for all cases in table 1.In figure 17, only positive values of TI are shown.Our interest in exploring (3.14) for analysing pulsating flows comes from three considerations.Firstly, the emergence a four-layer structure constitutes a robust criterion for assessing the progress of laminar-turbulent (Ebadi et al. 2019) or turbulent-turbulent (Guerrero et al. 2021) transition in unsteady flows.Secondly, TI is assumed to be the primary momentum source which supplies the overshoot in the perturbation velocity field (Sundstrom & Cervantes 2017) (see § 3.4).Finally, for a short period of time immediately after non-periodic acceleration or deceleration, the shear stress and VF constitute the only significant dynamic transient response, whilst the turbulent field remains frozen (Guerrero et al. 2021(Guerrero et al. , 2023)).This final point is critical for defining the evolution of high-frequency pulsating flows in the framework of the turbulent-turbulent transient concept, as discussed below.
We first consider the evolution of VF for the highest pulsating frequency, l + s = 5.For A b ≥ 0.5 (figure 16b,c), the transient response is dominated by a pair of alternating viscous layers which form at the wall and expand towards the core.One layer grows as the flow begins to accelerate, and contains a negative VF which acts as an amplified momentum sink.Within the second layer, which grows as the flow beings to decelerate, the sign of the amplified VF is reversed, such that it now serves as a momentum source.The purpose of each layer is to maintain the no-slip condition at the wall by counteracting the rapid change of momentum in the flow brought about by the driving force.Taken in isolation, these individual sink and source layers closely mimic the flow response in the inertial stages of non-periodic acceleration and deceleration, respectively (Guerrero et al. 2021(Guerrero et al. , 2023)).For clarity, we refer to the individual stages in the acceleration and deceleration periods as the 'inertial sink' and 'inertial source' stages, respectively.In non-periodic acceleration, the 'inertial sink' stage serves as a precursor to the pretransition stage.However, in case L05A05 (figure 16b) and case L05A10 (figure 16c), the pretransition stage is never reached, as this is replaced by the 'inertial sink' stage.Since these stages occur in immediate proximity, the existing viscous layer is forced away from the wall as the new viscous layer forms.The detached layer is then pushed into the core of the flow until it disperses.Whilst the presence of the existing layer does not appear to hamper the growth of the new source or sink layer, it does lead to some distortion of TI (figure 17b,c), which typically remains frozen during the inertial stage of non-periodic flows (Guerrero et al. 2021(Guerrero et al. , 2023)).During the inertial source stage, the turbulent inertia shows significant growth during the deceleration period within a region which slowly moves away from the wall.This overshoot in TI occurs when the boundary between the developing momentum source layer and the detached momentum sink layer passes through the buffer layer of the mean flow field.This suggests that the temporal changes to TI in high-frequency pulsating flow stems primarily from the presence of the detached momentum source layer which is not present in non-periodic unsteady flows.
From figure 16 it can be seen that the strength of VF in the viscous layers becomes weaker with decreasing amplitude and increasing Stokes length.More specifically, figure 16 suggests that the existence of these viscous source/sink layers is dependant on  the acceleration and deceleration rates of the flow.Guerrero et al. (2021) observed the existence of an inertial stage for transient forcing periods of up to t + ≈ 208 (rescaled based on Re τ in the present study).For l + s = 5, the equivalent forcing period of similar magnitude is t + = 39, which explains the characteristic inertial behaviour even for low amplitudes.However, since the acceleration rate is not independent of the time period   ).The mechanisms of turbulent production and redistribution are frozen.However, the process of decay from the preceding deceleration period remains active.The wall-normal and shear Reynolds stress components become gradually weaker away from the wall, whilst the spanwise component continues to decay throughout the domain.It is important to note that the necessity of defining a distinct residual decay stage is dependant on the relative position of φ = 0 in the cycle.In the turbulent-turbulent transient concept, where φ = 0 lies at the point of minimum bulk velocity, the residual decay stage is distinguished by the formation of a new perturbation boundary layer, the growth of which is taken to originate at φ = 0.When the position of φ = 0 is shifted to an arbitrary point in the cycle, this stage could justifiably be merged with the final stage of the preceding cycle, in order to simplify the respective model.
As for non-periodic accelerating flows, the pretransition stage (II) of pulsating flows is characterised by the stretching and elongation of high-and low-speed streamwise velocity streaks, whilst the production and redistribution of new turbulent energy remains relatively frozen.The perturbation flow field follows a self-similar laminar distribution, and remains confined within a region bounded by the periodic quasilaminar Stokes solution and non-periodic extended laminar Stokes solution.Reynolds stress components farther away from the wall continue to fall, as the propagation of old turbulence from the preceding residual decay and turbulence decay stages remains active, and serves as the dominant transient turbulent process during this stage.
The transition stage (III) of a pulsating flow follows a very similar development to the transition stage in a non-periodic accelerating flow.The elongated velocity streaks in the perturbation boundary layer begin to destabilise, triggering a bypass-like transition process.Turbulent spots begin to form around low-speed streamwise velocity streaks.The spots grow and merge together until the near-wall region is populated by new vortical structures.In pulsating flows, the transition stage is not limited to the acceleration period of the cycle.Turbulent spots which form late in the acceleration period are able to continue growing and merging even as the flow starts to decelerate.During pulsation, the frozen turbulence phenomenon, in which the effects of the deceleration on the turbulence flow field are initially delayed, extends to include ongoing turbulence transient processes, such as bypass transition, in addition to the existing turbulence in the flow.During the transition stage of the low-frequency regime, C f grows rapidly such that it collapses towards the quasisteady solution, whilst for the intermediate frequency regime C f plateaus.
The post-transition stage (IV) only occurs when the transition stage concludes before the end of the acceleration period.The Reynolds stress components settle around a pseudoequilibrium value in the near-wall region, whilst there is a continuing growth of these components in the core of the flow, as turbulence gradually propagates away from the wall.The flow remains in the post-transition stage during the early deceleration period, due to the delay in the turbulence response to the deceleration.
The turbulence decay stage (V) is characterised by a gradual decay of the Reynolds stresses, as the streamwise velocity streaks break up and the vortical structures diminish throughout the domain.The skin friction of the perturbation flow reduces in accordance with the theoretical solution for a quasisteady flow field.If the laminar-turbulent transition process concludes within the first half of the deceleration period, the flow immediately moves from the transition regime into the turbulence decay stage.In such cases, any potential overshoot of the maximum value of the streamwise Reynolds stress, which serves as a strong indicator of the occurrence of turbulent-turbulent bypass transition, is masked, as no time is given for the turbulent statistics to settle towards any pseudoequilibrium value.
In the high-frequency regime (l + s = 5), the vortical structures, streamlines and turbulent statistics all show a very weak response throughout the cycle.However, significant dynamic behaviour has been identified for the VF and TI terms of the mean momentum balance.During the acceleration period, a new layer of highly amplified VF grows at the wall, and acts as a momentum sink, in order to maintain the non-slip condition at the wall.Similarly, a second layer grows during the deceleration period, though this time with VF acting as a momentum source.The growth of such layers mimics the initial response, or inertial stages, observed in non-periodic acceleration and deceleration, which are precursors for either the onset of bypass transition or the suppression of turbulence in the flow.However, the acceleration period in the current high-frequency regime is insufficient for the boundary to interact with the buffer layer.The flow alternates between two consecutive stages, referred to here as the inertial sink stage during the acceleration period and the inertial source stage during the deceleration period.With this in mind, the relative inactivity of the turbulence in the high-frequency regime can be attributed to the failure of the flow to leave the inertial sink regime before the end of the acceleration period.One critical difference in the inertial stages of periodic flows, compared with non-periodic flows, is the existence an amplified VF layer from the preceding stage, which detaches from the wall and moves into the flow as the new layer grows.Therefore, TI is not completely inert, as its strength and distribution is influenced by the presence of the boundary between two highly viscous but opposing regions passing through the buffer layer.distance of x + 1 < 0.5L + x throughout the cycle.Furthermore, R 11z falls to zero within a spanwise distance of z + 1 < 0.5L + z throughout the cycle.Figure 19 indicates that the width and length of each domain is sufficient to capture at least two streak lengths in the streamwise direction, and more than two streak widths in the spanwise direction, including

Figure 1 .
Figure 1.Numerical set-up of the code for (a) spatial configuration of the channel domain and (b) cosine waveform of the driving force: (solid) A b = 0.1; (dash) A b = 0.5; (dot-dash) A b = 1.0.

Figure 3 .Figure 4 .
Figure 3. Wall-normal profiles of the time-averaged streamwise velocity and Reynolds stresses.All cases in table 1 are shown.The results of Weng et al. (2016) for a case of l + s = 10 at A uc = 0.1 (circles), and Manna et al. (2012) for a case of l + s = 3.1 at A uc = 1.0 (triangles) are shown for comparison.
Seddighi et al. (2014) identified such differences in the growth of turbulent spots and velocity streaks when comparing ramp-up and step-up acceleration in non-periodic accelerating flows at an equivalent Reynolds number.During the acceleration, an amplitude of A b = 0.5 (case L26A05) more closely mimics a ramp-up flow acceleration, whilst an amplitude of A b = 1.0 (case L26A10) more closely mimics a step-up acceleration.As in the ramp-up acceleration ofSeddighi et al. (2014), for an amplitude of A b = 0.5 the flow completes transition before the maximum bulk velocity is reached, whilst this is not true for the higher amplitude of A b = 1.0.

Figure 7 .
Figure 7. Phase-variance of the maximum Reynolds stress components and turbulent kinetic energy for cases (a) L26A05, (b) L26A10, (c) L16A05, (d) L16A10 and (e) L10A10.Arrows indicate the respective axis for the intersecting curves.Vertical blue lines indicate the initiation of turbulence growth/transition.The appropriate location of each line was determined from the analysis in § 3.3.

Figure 10 .
Figure 10.Growth of a secondary sinuous instability of the streamwise velocity streaks at y + = 12.4, and its propagation at y + = 106.4,for case L26A10.

Figure 11 .
Figure 11.Phase-variance of the Reynolds stress components and streamwise energy budget terms for case L26A05.
Figure 15 displays the transient behaviour of the perturbation skin friction coefficient during the pulsation cycle.For cases L16A05, L16A10, L26A05 and L26A10 the distinct stages of the transient-transient process, which are identified in § 3.3, are indicated in their respective plots.Each value of l + s is compared with the quasilaminar (Stokes) solution, C ∧ f (ql) , extended laminar (Stokes) solution, C ∧ f (ext) , and the quasisteady (turbulent) approximation, C ∧f (qs)  .In all cases, there is a notable similarity between the trend of the quasilaminar and extended laminar solutions for C ∧ f in terms of both the phase of the profile, and the relative amplitudes between their peaks.The extended laminar solution underpredicts the quasilaminar solution throughout the cycle.However, the relative magnitude of the discrepancy between these two solutions at different points in the cycle remains consistent for varying values of the Stokes length.The magnitude of|C ∧ f (ql) − C ∧ f (ext) | correlates with the value of C ∧ f (ql) , such that |C ∧ f (ql) − C ∧ f (ext)| grows as the flow accelerates, reaching a peak at the same point at whichC ∧ f (ql) reaches its maximum.Beyond this point |C ∧ f (ql) − C ∧ f (ext)| shrinks and the profiles converge, such that there is good agreement between the two solutions for a majority of the deceleration period.The high-frequency regime l + s = 5 (not shown) shows a strong agreement with the quasilaminar solution, with the amplitude exerting no influence on C ∧ f .This agreement is also observed for l + s = 10 (figure15a-c), particularly at the highest amplitude of A b = 1.0; however, as the amplitude is reduced, C ∧ f starts to undershoot C ∧ f (ql) as C ∧ f approaches its peak.This undershoot at l + s ≈ 10 is well known, with previous studies consistently observing a ratio between the shear stress amplitude of the oscillating flow field to that of the high-frequency Stokes solution of less than 1(Weng et al. 2016;Sundstrom & Cervantes 2018a).However, the observation that this undershoot exists independently of