Propagation mechanism of localized wave packet in plane-Poiseuille flow

Localized turbulent structures are revealed recently to be the key features near the onset of turbulence in linearly stable shear flows, e.g. puffs in pipe flow and oblique turbulent stripes or bands in channel flows Eckhardt et al. (2007); Tuckerman et al. (2020). For two-dimensional (2D) plane-Poiseuille flow, the corresponding structure is localized wave packet (LWP) Rozhdestvensky and Simakin (1984); Jimenez (1990); Price et al. (1993), whose relations with finite-amplitude periodic waves were analyzed theoretically and numerically Soibelman and Meiron (1991); Drissi et al. (1999); Mellibovsky and Meseguer (2015). LWP has a strong downstream edge and a slowly decaying upstream edge Jimenez (1990), and the corresponding decay and growth rates were explained in terms of the linear spatial modes Barnett et al. (2017). Similar asymmetry between the upstream and downstream sides was also found for three-dimensional coherent structures in channel flows Zammert and Eckhardt (2016). LWP in linearly unstable channel flows, where the Reynolds numbers are larger than the linear critical value 5772, is studied as well, but a damping filter is used in the simulations to restrain LWP from expanding Teramura and Toh (2016). The crucial questions for LWP during the subcritical transition are the following: What is the localization mechanism? What determines its convection velocity? What is the selection criterion for the dominating frequency? If the wavelength is not uniform in the streamwise direction, what is its selection criterion?

It is postulated that the subcritical transitions of shear flows may fall into the universality class of directed percolation (DP) Pomeau (1986); Sipos and Goldenfeld (2011); Shih et al. (2016); Lemoult et al. (2016); Pomeau (2016); Chantry et al. (2017). For plane-Poiseuille flow, recent experiments defined a critical Reynolds number $Re_{c}$ of 830 based on the DP power law Sano and Tamai (2016), while numerical simulations revealed that the DP power law is retrieved only as $Re$ is above 924 Shimizu and Manneville (2019). When $Re$ is far below these DP thresholds, it has been found numerically and experimentally that the localized turbulent bands can extend obliquely Xiong et al. (2015); Kanazawa et al. (2017); Paranjape (2019); Xiao and Song (2020), and the periodic turbulent bands can sustain in a sparse turbulent state Tao et al. (2018). By applying random initial disturbances, LWP density, the corresponding parameter of turbulence fraction for two dimensional plane-Poiseuille flow, was shown numerically to approach zero as $Re_{c}$ was approached from above, and it was concluded that the subcritical transition was more like a continuous phase transition rather than a first-order one Wang et al. (2015). It is known that the subcritical transition may start at different ${Re_{c}}^{\prime}s$ depending on different initial or upstream disturbances Mullin (2011). Are the random disturbances the most effective perturbation to trigger the transition at the lowest $Re_{c}$? Finding answers to these crucial questions is the motivation of this paper.

The incompressible two-dimensional Navier-Stokes equations are solved for plane-Poiseuille flow with a spectral code Chevalier et al. (2007). The half height of the channel $h$ and 1.5 times of the bulk velocity $U_{m}$ are chosen as the characteristic length and velocity, respectively. The flow rate is kept constant and the Reynolds number is defined as $Re=1.5U_{m}h/\nu$, where $\nu$ is the kinematic viscosity of the fluid. Periodic boundary conditions are used in the streamwise $x$ direction and no-slip conditions are imposed at the parallel walls ($y=\pm 1$). 65 Chebyshev modes in $y$ direction and 512 Fourier modes per 100 length units in $x$ direction are used. For details of the simulation methods, we refer to the previous papers Xiong et al. (2015); Tao et al. (2018).

Following the method proposed for turbulent bands Tao et al. (2018), the center’s coordinate $x_{c}$ and length $l$ of the localized wave packet at a given time are defined as,

where $e$ is the kinetic energy of the velocity field after deducting the basic flow. The packet velocity $c_{p}$ can be calculated by tracking $x_{c}$, e.g., $c_{p}=0.69$ is obtained based on the $x_{c}$ data during 2000 time units for an isolated LWP at $Re=2400$. It is found that $l$ is large and increases with $Re$ at moderate Reynolds numbers, e.g., $l$ increases from 16.3 at $Re=2500$ to 33.1 at $Re=4500$. In order to diminish the influences of periodic boundaries and obtain the intrinsic properties of an isolated LWP, a computational domain of at least 3$l$ long is required in simulations.

In a frame moving with the packet velocity $c_{p}$, which is referred to as $S$ frame hereafter, the envelope of LWP looks static and the velocity field is decomposed into three parts,

where the basic flow $U_{0}=1-y^{2}-c_{p}$, $\textbf{U}_{1}=(U_{1},V_{1})$ is the mean-flow modification, the mean flow after deducting $\textbf{U}_{0}$, and $\textbf{u}^{\prime}=(u^{\prime},v^{\prime})$ is the fluctuation velocity. From now on, $x$ and $y$ denote the coordinates in the $S$ frame. As shown in Fig.1(a), two rows of vortices lie near the top and bottom walls respectively, a typical feature of LWP Jimenez (1990). The fluctuation component $\textbf{u}^{\prime}$ travels upstream in a wave form in the $S$ frame [Fig. 1(c)], and is referred to as fluctuation wave (FW) in this paper. The mean flow in $S$ frame is calculated with 200 fields sampled every 10 time units. It is illustrated in Fig. 1(b) that the vorticity field of the mean-flow modification exhibits a prominent structure: a vortex dipole at the downstream end of LWP sandwiched between vortex layers extending upstream near the walls. The vortex dipole, to the best of our knowledge, has not been reported before, and is shown next to be responsible for the localization property, the packet velocity $c_{p}$, and the frequency of LWP.

In order to understand the dynamic behavior of FW, linear stability analyses are carried out for the mean flow ($\textbf{U}_{0}+\textbf{U}_{1}$) corresponding to Fig. 1(b) at different $x$ based on the parallel-flow approximation. The perturbations are assumed in the form of $\sim e^{i[(k_{r}+ik_{i})x-(\omega_{r}+i\omega_{i})t]}$, where $k_{r}$, $k_{i}$, $\omega_{r}$, and $\omega_{i}$ are solved from the local dispersion relation Schmid and Henningson (2001) and shown in Fig. 2. Several interesting features should be noted. First, there is a finite unstable region ($\omega_{i,max}>0$) surrounded by stable regions, i.e. the region between $x=21$ and 26.3 shown in Fig. 2(a), a main part of the vortex dipole region [Fig. 1(b)]. Second, $\omega_{r,max}$, the frequency of the most unstable temporal mode ($k_{r}>0,\ k_{i}=0$) is negative as shown in Fig. 2(b), representing an unstable wave mode traveling upstream in the $S$ frame just as the FW found in simulations. This result is reasonable by considering that the mode’s maximum amplitude lies at the neighborhoods of walls [Fig. 2(f)], where the mean flow moves upstream. These features suggest a localization mechanism for LWP: traveling wave mode amplified in the unstable region decays in the stable regions, forming a localized wave packet.

According to the simulations, FW propagates with a unique global frequency, e.g., $\omega_{g}=-0.345$ as shown by the blue line in Fig. 2(b). It is noted that though the frequency can be measured as $x>30$, the fluctuation velocity is too weak to be recognized as shown in Fig. 1(c) because the downstream side of the present LWP is the far upstream end of another LWP. In order to understand the selection criterion of the global frequency, spatio-temporal stability analyses are carried out based on the mean flow, and the absolute growth rate $\omega_{0,i}$ and absolute frequency $\omega_{0,r}$, where the group velocity is zero, are shown in Fig. 2. $\omega_{0,r}$ satisfying the saddle-point condition Chomaz et al. (1991); Hammond and Redekopp (1997) is computed as -0.39, which is different from $\omega_{g}$, and hence the saddle-point criterion seems not applicable for LWP. The absolute frequency at the upstream boundary of the absolutely unstable region ($\omega_{0,i}>0$) is -0.35, which is labeled by the red point in Fig. 2(b) and almost coincides with $\omega_{g}$, suggesting an marginal stability criterion for the global frequency selection Dee and Langer (1983); Monkewitz and Nguyen (1987).

Based on the fluctuation vorticity recorded during 2000 time units along the bottom wall, the phase velocity and then the wave number of FW are determined with $\omega_{g}$ at each $x$ position, and the spatial growth rate of FW is calculated based on the the envelope of the normal fluctuation velocity $v^{\prime}$ along the centerline. Since the wave following the dipole decays in the upstream direction, we look for the spatial mode with the minimal $-k_{i}$ at each $x$. As shown in Fig. 2(c) and 2(d), the spatial growth rate and wave number of FW acceptably agree with $-k_{i}$ and $k_{r}$ of the spatial mode with $\omega_{g}$, and the flow structure of FW shown in Fig. 2(e) agrees with that of the spatial mode [Fig. 2(f)] as well, indicating that the spatial properties of FW are mainly determined by the local dispersion relation of the mean flow.

The spatial instability of LWP was discussed before based on the basic flow and the frequency obtained in simulations, and constant $k_{r}$ and $k_{i}$ were solved and found to be consistent with the simulations at upstream and downstream tails of LWP Barnett et al. (2017). Such a consistency is an asymptotic case for the present study. As shown in Fig. 2, the spatial growth rate and wave number of FW are not constant but vary in the streamwise direction. However, at the far upstream tail of LWP, where the mean flow modification almost diminishes, $-k_{i}$ and $k_{r}$ are close to the asymptotic values corresponding to the basic flow as shown by the horizontal dashed lines in Fig. 2(c) and 2(d), respectively.

Besides FW, another key feature of LWP is its propagation velocity. As shown in Fig. 1(b), the vorticity field is antisymmetric about the centerline, and hence at the centerline $V_{1}$ and $\partial U_{1}/\partial y$ or $U_{1y}$ are zero. Consequently, the mean flow along the centerline may be described by a steady one-dimensional nonlinear model. Considering that the subcritical transitions occur as $1/Re\sim 10^{-3}$ and the streamwise fluctuation velocity $u^{\prime}$ is very small along the centerline as shown in Fig.2(e), the viscous diffusion of the mean flow modification and the Reynolds stress ($\overline{u^{\prime}v^{\prime}}$ and $\overline{u^{\prime}u^{\prime}}$) are ignored along the centerline. As a result, the gradient of the mean pressure modification $P_{1x}$ should mainly depend on the variation and the derivatives of the centerline mean velocity, i.e. $(1-c_{p}+U_{1})$, $U_{1x}$, $U_{1xx}$, $U_{1xxx}$, …. According to the Bernoulli equation, deceleration corresponds to an adverse pressure gradient, suggesting that $P_{1x}$ remains roughly the opposed phase to $U_{1x}$, an odd derivative of $U_{1}$, as shown in Fig. 3(a). Inspired by the fact that even derivatives will cause phase shift from the odd derivatives for a harmonic wave, only the odd derivatives of $U_{1}$ are considered and the dependency of $P_{1x}$ on $U_{1}$ is simplified to a linear relation:

where $A$ and $B$ are coefficients. As shown in Fig. 3(a), the pressure gradient estimated with the velocity derivatives (red curve) agrees with the numerical data (black curve) along the centerline of the vortex dipole, indicating that Eq. (3) grasps the main relation between $P_{x}$ and $U_{1}$. Note that $A=-0.347$ and $B=0.154$ are used in Eq. (3) to guarantee that the estimated $P_{1x}$ has the same minimum and maximum values as the numerical one.

Substituting Eq. (3) into the mean $x-$momentum equation, the steady centerline model for the mean flow becomes:

This is a KdV-type equation for a steady solitary wave in the moving $S$ frame or equivalently, a solitary wave traveling with a velocity of $c_{p}$ in a motionless frame. Its solution can be solved easily as:

where $x_{1}$ is a constant defining the $x$ coordinate of the maximum $U_{1}$. As shown in Fig. 3(b), the maximum velocity predicted by the centerline model [Eq. (5)] is consistent with the corresponding simulation value. Therefore, when the maximum mean velocity at the centerline of the vortex dipole region is given, $c_{p}$, the convection velocity of the solitary wave and LWP in a motionless frame, is determined by Eq. (5). Interestingly, though $A$ and $B$ used in Fig. 3 are determined with the numerical data for $Re=2400$, the nonlinear centerline model [Eq. (3)-(5)] does not include the Reynolds number explicitly, suggesting that $c_{p}$ does not strongly depend on $Re$. This suggestion agrees qualitatively with the numerical simulations, where the convection velocity of LWP increases slightly with $Re$, e.g. $c_{p}=0.69$, 0.74, and 0.78 for $Re=2400$, 3000, and 4000, respectively.

It should be noted that without the Reynolds stress contributed by the finite-amplitude FW the vortex dipole will decay due to viscous diffusion and dissipation. It is checked numerically that both the finite-amplitude $\textbf{U}_{1}$ and the finite-amplitude FW are necessary to sustain LWP, denoting that there exists a threshold energy and the subcritical transition is a first order transition. LWP becomes longer with the increase of $Re$, and was found to split at $Re=5000$ Jimenez (1990). According to the present simulations, the critical Reynolds number for an isolated LWP to split in a long domain is about 4950, above which LWP splitting will cause more LWPs and the whole domain will be occupied by LWPs eventually, leading to a statistically steady or equilibrium state.

When we decrease $Re$ gradually from the equilibrium state, e.g. $Re=6000$, the perturbation kinetic energy in the whole domain area $S$, $E_{k}=\frac{1}{2S}\int_{S}\mid\textbf{U}_{1}+\textbf{u}^{\prime}\mid^{2}dS$, decreases but does not vanish if $Re$ is larger than 2331, where 9 LWPs are reserved in a domain of 400 units long. In order to examine the $Re$ threshold for sustained LWP, different numbers ($N$) of a sample LWP obtained at $Re=2350$ are evenly spaced and introduced as initial perturbations. As shown in Fig. 4(a), relaminarization occurs when $Re$ becomes lower than a critical value $Re_{c}$, and the same $Re_{c}=2332.5\pm 0.5$ corresponds to different $N^{\prime}s$ [Fig. 4(b)] or LWP densities (the LWP number per unit streamwise length) less than 0.017, indicating that when LWPs are far from each other they behave just like an isolated one. According to the dashed curve shown in Fig. 4(b), the lowest $Re_{c}$ or threshold $Re$ for sustained LWP is $2330.5\pm 0.5$ with a nonzero threshold LWP density about 0.022, confirming again that the present transition is a discontinuous type. The LWPs initially arranged too tightly (e.g. the case of $N=10$ shown in Fig. 4) tend to decay at low Reynolds numbers due to LWP interaction. Considering that random initial disturbances are not as effective as LWP themselves to trigger LWPs and the LWPs caused by random disturbances through transient growth may stay tightly and tend to diminish, using random initial disturbances may lead to less sustained LWPs than localized initial perturbations.

According to the present study, the localized wave packet found in two-dimensional plane-Poiseuille flow represents a symbiosis between the vortex dipole in the mean flow modification and the fluctuation wave: the dipole defines the convection velocity of the whole packet with the solitary wave velocity at its centerline, provides an unstable region to amplify FW and a global frequency for FW, while FW feeds back the Reynolds stress to prevent the vortex dipole from decaying, and travels upstream obeying the local dispersion relation of the mean flow. Therefore, nonlinear effects are necessary to sustain LWP, and the subcritical transition is a first order transition with a nonzero threshold LWP density. In addition, isolated LWP can sustain as $Re<4950$, suggesting that the initial transition stage is characterized by a sparse structure state instead of an equilibrium state, which can be achieved only as $Re>4950$ due to the wave packet split.