1 Introduction

The linear stability theory is concerned with whether a laminar base flow changes its state when perturbations with infinitesimal amplitude are added. If the state changes, the base flow is unstable. If the base flow returns to the original state, it is stable. To date, several successful numerical methods have been developed successfully for the linear stability analysis, such as the modal stability analysis, the non-modal stability analysis, and the so-called direct optimal growth analysis.

Traditionally, the linear stability of a flow has been analyzed using the modal approach. In this approach, the asymptotic behavior of small perturbations to a steady or time-periodic base flow is usually considered. This behavior is determined by the eigenvalues of the linear operator derived from the linearized Navier-Stokes (N-S) equations, because the eigenvalues characterize the asymptotic evolution of the perturbations. For some flows, especially for those with instability driven by thermal or centrifugal forces, e.g., Rayleigh-Bénard convection and Taylor-Couette flow, the predictions of modal analysis agree well with laboratory experiments^[1-2].

Despite remarkable accomplishments of the modal analysis, many questions are left unanswered. For example, discrepancies exist between the computed critical Reynolds numbers and the observed ones in many wall-bounded shear flows, and the theoretically predicted structures are not always observed in unforced experiments^[3]. The drawback of this modal approach is associated with the non-orthogonality of the eigenmodes of the linearized flow system. As a consequence, non-modal analysis theories^[4] are needed for many problems. These theories, referred to here as the transient growth theories, emphasize the linear nature of the non-modal amplification mechanism and are based on the observation that an initial perturbation that is not a pure eigenmode may undergo transient growth, even though all eigenmodes decay monotonically. The transient growth may lead to nonlinear instability or otherwise change the path of instability, leading to, e.g., bypass transition to turbulence. Among this form of initial perturbations, the one which yields the largest amplification is referred to be "optimal".

If the evolution of the perturbations is governed by the Orr-Sommerfeld and Squire (OSS) equations^[5], it is possible to directly evaluate the eigenvalues of the operator matrix. However, it is difficult, and in some cases impossible, to build the operator matrix for general complex base flows. To address the problem, a direct optimal growth method (also called the matrix-free method) for both modal and non-modal instability was recently presented by Barkley et al.^[6]. Because this method is suitable for the stability analysis of flows with arbitrary complexity, it has been applied in stability analysis of various flows, e.g., flow over a backward-facing step^[7], stenotic flows^[8-10], flow past a circular cylinder^[11-12], vortex pair systems^[13], flow through a sudden expansion in a circular pipe^[14], and flow over a turbine blade^[15]. More recently, Mao et al.^[16] demonstrated that both optimization and eigenvalue (called the direct optimal growth method by Barkley) approaches converge to the same outcome, and they applied both methods to calculate both optimal initial and boundary condition problems. In the direct optimal growth method, a Krylov subspace was constructed to approximate the eigenmodes of large matrix by explicitly solving the linearized N-S equations and their adjoints by iterations. This method, however, depends on the specific boundary conditions of the flow. Thus, the linearized adjoint N-S equations should be determined if it is applied to base flows with other boundary conditions.

In this paper, we improve the direct optimal growth method and apply it to the plane Poiseuille flow in a channel. A detailed derivation of the linearized adjoint N-S equations of the channel flow with specific boundary conditions is presented. A direct numerical simulation (DNS) code based on the spectral method is then used to solve the linearized N-S equations and the corresponding adjoint ones. The Krylov subspace algorithm is refined by using a re-orthogonalization Arnoldi technique to improve the orthogonality of the orthogonal basis of the Krylov subspace. The improved algorithm is then combined with the DNS method to investigate the temporally developing global instability in such flows. To validate the algorithm, we examine the one-dimensional plane Poiseuille flow and the two-dimensional plane channel flow perturbed by low speed streaks. A mechanism that induces transient growth of perturbations in the channel flow is discussed.

2 Numerical method 2.1 Transient growth of infinitesimal perturbations

Considering a general incompressible flow U in a three-dimensional domain Ω, the linearized N-S equations for infinitesimal perturbations can be expressed as

(1)

Here, u and p are the infinitesimal perturbation velocity and kinematic pressure, respectively, υ is the fluid viscosity defined by υ = Q/Re, where Q is the constant-volume-flow rate, and Re is the bulk Reynolds number, and U is the base flow.

The linear evolution of a perturbation from t= 0 to a later time t under Eq. (1) can be expressed concisely as the action of a linear evolution operator A(t) on the initial perturbation u(x, 0),

(2)

Usually, the modal analysis focuses on the asymptotic behavior of the perturbations, which is characterized by the eigenvalues of the system and the evolution of the eigenmodes. However, owing to the non-orthogonality of the eigenmodes of the linearized N-S equations, the dynamics of interest may not be in the form of an exponential function of time multiplying a fixed modal shape. Therefore, the eigenvalue problem in modal analysis is not directly relevant^{[2, 5]}. Rather, a transient growth method is preferred to quantify such dynamics, which is concerned with the maximum energy growth for all possible initial perturbations over a finite time interval.

Typically, the total kinetic energy E of a perturbation field over the full flow domain is chosen to quantify the size of the perturbations^[5], derived from the L₂ inner product of the perturbed velocity u,

(3)

where is the standard L₂ inner product in the space. Transient growth is described by the growth of the energy norm of the perturbations over a given time interval and can be quantitatively measured by the ratio of final energy at the time t to initial energy at the time zero. Setting the norm of the initial perturbations to unit, i.e., the transient energy growth over the interval t is

(4)

where A^*(t) is the adjoint evolution operator of A(t).

The goal of transient growth analysis is to find the maximum energy growth and the corresponding initial perturbations, i.e., the optimal growth and the optimal perturbations. From Eq. (4), it is obvious that seeking the optimal growth is equivalent to finding the leading eigenvalue of operator A^*(t)A(t), and the corresponding eigenmode is the optimal perturbation. Let λ_j and v_j denote an eigenvalue and the corresponding normalized eigenmode of A^*(t)A(t), respectively. Then, we have

(5)

Thus, the maximum energy growth at the time t, denoted as G(t), can be derived as

(6)

2.2 Linearized adjoint N-S equations

A^*(t) in Eq. (4) represents the evolution operator for the linear adjoint equations of the adjoint perturbations with the same boundary conditions as the linearized N-S equations. In this section, we derive the adjoint N-S equations with the specific boundary conditions for the channel flow.

The adjoint N-S equations can be divided into three parts: the advection term, the viscous and pressure term, and the time derivative term. Let H be an operator representing one of these three terms. Then, H and its adjoint operator H^* must satisfy the following relation based on the L₂ inner product in the domain of space Ω and time [0, τ]:

(7)

for arbitrary functions u and u^* with homogeneous boundary conditions on the walls and periodic boundary conditions in the streamwise and spanwise directions.

For the advection term, we define an advection operator N as

(8a)

or

(8b)

According to the definition of inner product in Eq. (3), one obtains

(9)

The boundary conditions in the channel flow presented in this study are such that the perturbations on the no-slip wall are zero and are periodic in the streamwise (x) and spanwise (z) directions. Thus, the surface integral in Eq. (9) is zero, and the adjoint advection operator N^* can be expressed according to Eq. (7),

(10a)

or

(10b)

The viscous and pressure terms in Eq. (1), as well as the continuity equation, are all linear, and hence can be treated together. These terms can be written as

(11)

According to the definition of inner product in Eq. (3), we have

(12)

Note that and in Eq. (12) by the divergence theorem and boundary conditions. Thus, Eq. (12) can be rewritten as follows:

(13)

From Eqs. (7) and (13), we can see that the adjoint operator of the coupled linear operator is the same as itself, that is, this operator is self-adjoint.

The adjoint operator of the time derivative term in Eq. (1) can be deduced as follows:

(14)

Equation (14) shows that the adjoint of the time derivative term is

According to Eqs. (10), (13), and (14), the adjoint of the linearized N-S equations for infinitesimal perturbations are written as

(15)

2.3 Algorithm of direct optimal growth

Equation (5) indicates that the main purpose of the transient growth is searching for the leading eigenvalue (λ_max) of the combined operator A^*(t)A^*(t) and the corresponding eigenvector. In the direct optimal growth method, a Krylov subspace is first constructed by applying the operator M to the initial vector u₀ repeatedly, where the operator M represents A^*(t)A^*(t). The Krylov subspace can be obtained by repeatedly integrating the linearized forward and subsequently adjoint N-S equations starting from u₀. Subsequently, a standard QR decomposition is employed to obtain the orthogonal basis of the Krylov subspace, which is used to obtain the eigenmode information of the linearized system. However, it is sometimes ill-suited to construct the orthogonal basis of the Krylov subspace by the standard QR decomposition. Sometimes, the application of the Krylov subspace on specific flows is ill-conditioned, and the orthogonality of the orthogonal basis of the Krylov subspace deteriorates in the construction process of the Krylov space. In practice, it is the orthogonal basis of the Krylov subspace, not the Krylov subspace itself, that is used. Thus, it is preferred to use the classical Arnoldi technique^[17] to improve orthogonality, which does not construct the Krylov subspace but directly calculates its orthogonal basis by the standard QR decomposition. In addition, to further improve orthogonality of the orthogonal basis, a re-orthogonalization technique is employed on the basis of the classical Arnoldi method. The performance of this technique will be discussed in detail in Section 4.

We now give the algorithm of direct optimal growth analysis based on the re-orthogonalization Arnoldi technique, by which we compute the leading eigenvalue (i.e., the maximum growth) and the corresponding eigenvector (i.e., the optimal initial perturbation) of the linearized channel flow.

Define a set of orthogonal basis of a k-dimensional Krylov subspace as

(16)

(ⅰ) Initialization: Provide a positive integer k_max as the maximum dimension of Q_k in Eq. (16), a residual norm tolerance n_tol (which is the measurement of convergence in iterations), and an initial vector q₀ with unit norm. Set k=0.

(ⅱ) Integrate linearized forward and adjoint N-S equations (Eqs. (1) and (15)), taking qk as the initial perturbation:

(ⅲ) Perform orthogonalization procedure of classical Arnoldi method by using a modified Gram-Schimdt approach:

(ⅳ) Perform the re-orthogonalization procedure to improve the orthogonality of orthogonal basis of the Krylov subspace,

(ⅴ) Calculate the leading eigenvalue λ_max and the corresponding eigenvector v_max of H, which consists of h_ik from Step (ⅳ). In other words, h_ik is the element of H. Compute the Ritz eigenvector Q_kv_max. Calculate h_{k+1, k} and the residential norm r_n as

where λ_max^(k) and λ_max^(k-1) are the leading eigenvalues in kth and (k-1) th iterations, respectively.

In Step (ⅱ), we modify the DNS code based on the standard Fourier-Chebyshev spectral method^[18] to integrate the linearized forward and adjoint N-S equations. In this simulation, the Chebyshev-τ method and the no-slip condition are used in the non-homogeneous wall-normal direction (y-direction), while the Fourier expansion is used in the homogeneous directions, i.e., the streamwise (x) and spanwise (z) directions. A Chebyshev-τ influence-matrix method, including a τ-correction step, is employed for the viscous and pressure term to ensure that the computed solutions satisfy both the incompressibility constraint and the momentum equation. The aliasing errors in the x-and z-directions are removed by truncation according to the 3/2-rule. Time advancement is carried out by using a semi-implicit backward-difference scheme with the third-order accuracy. This numerical method has been well tested in our previous studies^[19-20].

Note that the magnitude of h_{k+1, k} is used to judge the convergence of the iterations in Ref. [6]. However, the convergence of h_{k+1, k} is not necessary in this study, since we are focusing on the optimal growth. Therefore, we take as the measure of convergence in our computations, and discuss the validity in Section 4. In this study, we set n_tol=1×10^-4. Usually, convergence can be obtained in about 6 to 13 iterations.

3 Results and discussion

In this paper, we study the transient growth of perturbations for two base flows, that is, a one-dimensional plane Poiseuille flow and a two-dimensional Poiseuille flow with a straight streamwise low-speed streak. Actually, the algorithm presented in this work is applicable to any complex base flow in a channel, i.e., there is no limitation on the form of U in Eqs. (1) and (15).

3.1 Transient growth of plane Poiseuille flow

In this section, we choose a one-dimensional plane Poiseuille flow solution as the base flow U, whose streamwise velocity profile is given by

(17)

where y is the coordinate in the direction normal to the channel walls, and y∈[-1, 1].

The reason to choose the Poiseuille-type base flow in this work is twofold. Firstly, we intend to validate the algorithm by the OSS method. Secondly, the effects of the spanwise width L_z and the Reynolds number (Re) of the channel flow on the transient growth of perturbations are investigated. The channel, i.e., the computational domain, has a fixed streamwise length of L_x=π, and a fixed wall-normal height of L_y=2 (ranged from -1 to 1), and a variable spanwise width L_z. Table 1 gives the computational cases with different values of L_z and Re.

Figure 1 gives the profiles of transient growth (optimal growth) G(t) of perturbations at different L_z for Re=2 670. One can see that, in all cases, G(t) increases with time until it reaches the maximum value (defined as G_max), and then decreases. Besides, Fig. 1(a) shows that the maximum growth G_max in each case increases with the spanwise width L_z for L_z between π/3 and π. Beyond L_z >π, G_max no longer increases, as shown in Fig. 1(b). Figure 2 gives the variations of G(t) for various spanwise widths at Re=1 000. The profiles of G(t) and corresponding G_max for Re=1 000 show similar behaviors with those at Re=2 670.

Figure 3 further depicts the variations of the maximum growth G_max with the spanwise width for both Re. It is clear that G_max increases at small L_z until it reaches its largest value at L_z=π and then fluctuates and eventually tends to a constant. In addition, G_max shows similar behaviors for both Re, except that the value of the maximum growth for Re=1 000 is smaller than that for Re=2 670 at the same spanwise width. Figure 4 plots t_max vs. G_max for the two Reynolds numbers, where t_max is defined as the time when G(t) reaches G_max. One can see a good linear relationship between t_max and G_max. This relationship implies that the variation of t_max along with L_z is similar to G_max vs. L_z in Fig. 3.

To understand the dependence of G_max on L_z, we choose three typical spanwise widths of the computational domain to inspect the optimal streamwise velocity perturbation u at t_max for Re=1 000, as shown in Fig. 5. Note that the structures of the optimal perturbation field are x-independent. Therefore, only structures in the yz-section are shown. It can be seen that the alternating low speed and high speed regions in the z-direction are observed. In addition, such alternating regions occur both in the upper and bottom parts of the channel with half phase-shift in the z-direction. The high and low speed regions actually represent the high-speed and low-speed streaks in the optimal flow field. The occurrence of these streaks implies that the most dangerous state to destabilize the flow is the streaky structure which has been observed in transitional and developed boundary layer flows^[21]. Due to periodicity in the z-direction, the streaks always appear in pairs (one low-speed streak and one high-speed streak). Smaller spanwise width leads to narrower spacing between the low and high speed streaks (see Fig. 5(a)), whereas for the spanwise width L_z=π, the spacing between the streaks is optimal and the appropriate shape of the streaks is shown (see Fig. 5(b)). When the spanwise width of the computational box is further extended to 2π (see Fig. 5(c)), the shape and the spacing of the streaks remain the same as those in the case with the spanwise width π. These results indicate that π is the critical size to attain the maximum growth of the perturbations, as also shown in Fig. 3.

To understand the special role of L_z=π, we consider a wave-like solution in the z-direction where β is the spanwise wave number. We have

(18)

where n is a positive integer. Thus, if L_z < π, we have we have β>2n. Hence, β>2 since n≥1. As was documented in Ref.[5], G_max reaches its maximum at β=2. Therefore, G_max cannot be the maximum if L_z < π. Moreover, β=2n when L_z=π. Thus, β=2 and G_max reaches its maximum, when n=1 according to Ref. [5]. This argument is consistent with Fig. 5(b), which shows that the optimal perturbation contains only one period. Hence, n=1. Figure 5(c) is always consistent with β=2, corresponding to L_z=2π and n=2.

The results of Fig. 3 also show that the critical size π is independent of the Reynolds number. Thus, we can explore the variations of G_max/Re² along with β, as shown in Fig. 6. It is obvious that the curves overlap for both Re=2 670 and Re=1 000, and G_max/Re² reaches the peak when β=2. In other words, the maximum growths reach the peak at β=2 regardless of variations of Re.

3.2 Transient growth of two-dimensional flow

The results in Subsection 3.1 show that the most dangerous state to destabilize the flow is the streaky structure. Thus, we further study the transient growth of a two-dimensional Poiseuille flow with an initial straight low-speed streak in the streamwise direction. The two-dimensional base flow U(y, z) is adopted based on the form proposed by Schoppa and Hussain^[22] and can be represented as follows:

(19)

where U₀(y) is the Poiseuille flow as shown in Eq. (17), Δu is the streak velocity defect (Δu=0.8), β is the spanwise wave number (β=1), and Δh represents the most active height of the streak (Δh=0.35). Based on Eq. (19), we can obtain the two-dimensional base flow for the Poiseuille flow with a low-speed streak located at the bottom wall of the channel, as shown in Fig. 7.

For the purpose of comparison, we set L_z=π and Re=2 670 in the two-dimensional base flow, corresponding to one of the cases in the one-dimensional plane Poiseuille flow calculation (see Subsection 3.1). Figure 8 shows the comparison of the growth factor between the one-and two-dimensional base flows at the given L_z and Re. It can be seen that G(t) in the two-dimensional case grows more quickly than that in the one-dimensional case. G_max(=714.8) in the two-dimensional case is larger than that (G_max=502.4) in the one-dimensional case. Also, t_max(=105) in the two-dimensional case is smaller than that (t_max=125) in the one-dimensional case. The results imply that the two-dimensional base flow where the initial low-speed streak is present is more susceptible to perturbations and may break down and transit into the turbulent state more easily.

Figure 9 gives the distribution of optimal streamwise velocity perturbation at t_max=105 for the two-dimensional base flow. Due to the presence of the initial low-speed streak at the bottom of the channel, the optimal perturbations show a high-speed region and two low-speed regions at the lower half channel. The high-speed and low-speed regions are observed on the strongest shearing region in the base flow introduced by the low-speak streak. The upper half channel still presents a pair of streaks, i.e., alternating high-and low-speed regions, which is similar to the optimal perturbation for the one-dimensional base flow.

The strength of the streak is a critical factor controlling the instability of the base flow. We thus further study the effects of streak strength on the transient growth of the perturbations in the two-dimensional case. A base flow vortex line inclination angel θ on the streak flank is used to represent the strength of the streak. As proposed by Schoppa and Hussain^[22], θ is given locally by θ=arctan(|Ω_y|/|Ω_z|), where Ω_y and Ω_z denote the normal and spanwise vorticities, respectively. We select the angle at y=-0.8 (corresponded to a distance 0.2 from the bottom wall of the channel) as the measure of the streak strength. Specifically, we use θ_0.2 to denote the angle. For the base flow given by Eq. (19), its expression is given by θ=arctan(Ω_y|_max/(dU₀/ dy))_y=-0.8, where Ω_y|_max=βΔu/2, and dU₀/dy is the mean flow spanwise vorticity.

We compute the transient growths of perturbations for the two-dimensional flow with low-speed streaks with different strengths. The strength of the streak is changed by varying Δu and Δh=0.35 fixed. Figure 10 gives the variations of G_max and t_max as a function of the strength of the low-speed streak in the base flow, in which θ_0.2=0 corresponds to the case of one-dimensional base flow without streak. It can be observed that, with the increase of θ_0.2, G_max rapidly rises in a quadratic form, whereas t_max reduces linearly. These results indicate that a stronger low-speed streak leads to larger G_max in earlier time. Therefore, compared with the pure one-dimensional Poisueille flow, this flow, with the presence of the low-speed streak, is more efficient to excite perturbation growth and to promote the subsequent transition to turbulence.

3.3 Mechanism that induces transient growth

Due to similar behaviors of the maximum growth of perturbations at different Reynolds numbers in the one-dimensional case (see Figs. 3 and 6), we choose the case with Re=2 670 and L_z=π to further analyze the mechanism inducing the transient growth of the perturbations.

Figure 11 shows the evolution of the component-wise root mean square with time for both one-and two-dimensional base flows. It can be seen that both cases are similar to each other. The perturbations in the wall-normal and spanwise components (v_rms and w_rms) are somewhat larger than or approach to those in the streamwise component (u_rms) initially, while the streamwise component rapidly increases and reaches the maximum at t=125 and t=105 for the one-dimensional and two-dimensional cases, respectively (see vertical lines in Fig. 11). These results indicate transfer of energy from the mean base flows to u_rms due to transverse perturbation v_rms and w_rms during the time evolution as a consequence of the lift-up mechanism^[4].

Figure 12 shows the distributions of streamwise vortices and velocity vectors of the perturbations in the yz-section at time t=1 and t=125 at the optimal initial condition for the one-dimensional base flow. At an early stage (see Fig. 12(a)), there is a pair of strong positive and negative streanwise vortices alternating in the z-direction, due to larger v_rms and w_rms at t=1 in Fig. 11(a). The drastic upward and downward perturbation velocities at z=π(0) and respectively, induce the perturbative motion in the wall-normal direction which becomes the source of the lift-up mechanism. As time passes (see Fig. 12(b)), perturbation velocities in the yz-section and the corresponding streamwise vorticity decrease. The energy is extracted and transferred to the streamwise direction which leads to stronger streamwise perturbation velocity at t=125, as shown in Fig. 11(a).

The results of transient growth for the two-dimensional case also manifest the effects of the lift-up mechanism, as shown in Fig. 13. Although the distributions of streamwise vortices and perturbation velocity at the initial time (see Fig. 13(a)) in this case are different from those in the one-dimensional case (see Fig. 12(a)), they are similar with each other at the optimal time (see Figs. 12(b) and 13(b)) except for higher magnitude of streamwise vorticities and larger velocity vectors for the two-dimensional case, which suggests a stronger lift-up mechanism.

4 Orthogonality of eigenvectors and validation of algorithm

As discussed in Subsection 2.3, we use the re-orthogonalization Arnoldi technique to improve the orthogonality of the algorithm. Figure 14 checks the orthogonality of the orthogonal basis of the Krylov subspace for both one-and two-dimensional cases. Here, the orthogonality measure used in the figure is defined as in which I is the unit matrix, and Q is the orthogonal basis which is derived by the standard QR decomposition of the Krylov subspace (Method Ⅰ), the classical Arnoldi technique (Method Ⅱ), and the re-orthogonalization Arnoldi technique (MethodⅢ), respectively.

The results of Fig. 14 indicate that the orthogonality of the basis calculated by Method Ⅰ is much worse than that by the other two methods (Methods Ⅱ andⅢ), deteriorating quickly with iterations for both base flows. Figure 14(b) shows that the differences of the orthogonality of the three methods are even more significant in the two-dimensional base flow, and MethodⅢ, i.e., the re-orthogonalization Arnoldi method, preserves the orthogonality much better. The results suggest that this method will have a much better performance than the others in complex flow fields.

Furthermore, we test the validity of the algorithm. Recall that for the plane Poiseuille flow, the perturbation growth reaches its maximum at β=2 (see Fig. 6), regardless of the Reynolds number. This conclusion is consistent with the result obtained at Re=1 000 by Reddy and Henningson^[23] and that at Re=5 000 by Butler and Farrell^[4], both computed from the OSS equations for the Poiseuille flow. To further validate the algorithm presented in this work, we compute the optimal growth at L_z=π (corresponded to β=2) and Re=2 670 for both one-and two-dimensional base flow cases in two other ways. One is to use the OSS equation with β=2 given. The other is to cross-check the growth of the optimal initial perturbation over time^[6] with the spanwise width L_z=π. For the cross-check, the perturbation growth can be expressed as

(20)

where u(0) is the optimal initial perturbation, and u(t) is the perturbation at time t, which is carried out by linearized forward DNS codes and taking u(0) as initial perturbation.

Figure 15 depicts the validity of the algorithm in the present study by comparing the G(t) curves with those obtained by the OSS equations and by the cross-check method. The one-dimensional flow (see Fig. 15(a)) shows that the profile of G(t) computed by the algorithm agrees well with that by the OSS equations. However, the cross-check method slightly under-predicts G(t). This is probably due to the mesh resolution of DNS. The maximum values of perturbation growth (G_max) for our algorithm, the OSS method, and the cross-check method are 502.4, 502.6, and 483.6, respectively. Therefore, the maximum relative error of our result is only 0.04% compared with the OSS method. Moreover, Fig. 15(b) shows the comparison of G(t) with that calculated from the cross-check method for the two-dimensional flow. It shows about 3% of maximum relative error. Therefore, the application of the algorithm in the complex flow field is feasible. Note that since the OSS equations are not applicable to calculate the transient growth of perturbations for complex flow fields, only the cross-check method is used to validate the result of the two-dimensional flow.

5 Conclusions

In this study, we improve the direct optimal growth method, which is proposed by Barkley et al.^[6], and develop an algorithm that computes the transient growth behavior of perturbations in channel flows. In this algorithm, a re-orthogonalization Arnoldi technique is adopted to improve the orthogonality of orthogonal basis of the Krylov subspace. The linearized adjoint N-S equations are deduced to meet the specific boundary conditions of channel flow. In addition, a convergence criterion based on the residential norm of leading eigenvalue is proposed to obtain the solutions.

The algorithm is first applied to the one-dimensional plane Poiseuille flow in a channel box. The effects of spanwise width of the channel and Reynolds number on transient growth of perturbations are discussed. The largest optimal growth of perturbations is found to be at the spanwise width L_z=π for both Re=1 000 and Re=2 670 cases. The width corresponds to the spanwise wave number β=2, a value that has been observed in the literature. Next, the algorithm is applied to a two-dimensional base flow, namely, a plane Poiseuille flow with a low-speed streak in the streamwise direction. The results show that the transient growth of perturbations in the two-dimensional flow develops more quickly and is stronger than that in the one-dimensional flow. The effects of the strength of the low-speed streaks on transient growth are discussed in the two-dimensional case, using a vortex line inclination angle of the base flow as a measure of the strength. The discussion shows that the strength of the streaks is a critical factor to excite the perturbation growth. The optimal flow field shows that the structures of high-and low-speed streaks and streamwise vortices induce the transient growth of perturbations via the lift-up mechanism for both one-and two-dimensional base flows.

The re-orthogonalization Arnoldi technique in the algorithm is tested and proved to be more accurate in more complex base flows in this study. The algorithm is validated by the results of the OSS equations and of the cross-check method for both one-and two-dimensional Poiseuille-type base flows. However, it is not limited to such flows and is applicable to any complex base flows in a channel.