1 Introduction

There are two critical Reynolds numbers for channel, pipe, and boundary layer flows. One is the critical Reynolds number for the existence of small amplitude unstable disturbance waves, which is denoted by Re_cr1. The other is the critical transition Reynolds number, which is denoted by Re_cr2.

For incompressible plane channel flows, Re_cr1 =U_CL^*h^*/υ^*=5 772 ^[1], where U_CL^* is the laminar flow velocity at the channel center, h^* is the half channel width, and υ^* is the kinematic viscosity. It is first obtained by the linear stability theory (LST), and then confirmed by experiments. For incompressible pipe flows, the LST does not yield a finite value for Re_cr1, and it is widely believed to be infinity. Experiments do not yield a definite value, which depends on the environment and the initial state of the experiment. For Re_cr2, no theoretical result is available, but experiments yield a value around 2 300 ^[2-3] with the reference velocity as the mean turbulent flow velocity and the length scale as the pipe diameter.

For the incompressible boundary layer on a flat plate, Re_cr1 =U^*δ^*/υ^* = 519 ^[4], where U^* is the velocity at the edge of the boundary layer, and δ^* is the displacement thickness. The result is also confirmed by experiments. For Re_cr2, experiments show that, it does exist (see Fig. 1) ^[5]. However, no exact value of Re_cr2 has ever been determined. For compressible boundary layers, the situation is even complicated, and the critical Reynolds number depends not only on the Reynolds number but also on the Mach number, the temperature condition on the wall, etc. Therefore, there is no unique critical transition Reynolds number.

The determination of the critical transition Reynolds number Re_cr2 for compressible boundary layers is important for some engineering problems, e.g., the scramjet engine of hypersonic vehicles. The boundary layer at the first section of an inlet is expected to be turbulent, otherwise the performance of the engine may be deteriorated. However, natural transition usually does not take place early enough. Therefore, a certain artificial measure has to be used to trigger the transition. If the Reynolds number at the expected transition location is lower than Re_cr2, transition cannot be triggered as expected, no matter what measure is taken.

So far, there has been no theoretical method for this purpose, and experimental methods cannot solve such problems. Therefore, numerical simulation method is the only possible candidate.

In this paper, we will first explore how numerical method can be used to determine the critical transition Reynolds number for plane channel flows. Then, the method will be extended to determine the critical Reynolds number for the transition of supersonic/hypersonic boundary layer flows in the subsequent papers. First, the direct numerical simulation (DNS) method is carried out to obtain a turbulent flow. Obviously, the corresponding Reynolds number is higher than Re_cr1. Then, we take the flow field to obtain the initial flow field for the DNS method, with which the Reynolds number is lowered step by step until the turbulent flow can no longer be maintained. The obtained smallest Reynolds number for the turbulent flow is just the critical transition Reynolds number.

For experiments, obviously, one can start either from a laminar flow, i.e., with a small Reynolds number, and then gradually increases the Reynolds number until the flow becomes turbulent, or start from a turbulent flow, and then gradually decreases the Reynolds number until the flow becomes laminar. There is no difficulty in using the first approach. However, in the second approach, when the Reynolds number approaches and crosses the critical transition Reynolds number, the turbulent flow will take a long time or distance to become laminar, which may cause a certain ambiguity in the determination of the critical Reynolds number. Therefore, in experiments, people prefer to start with a Reynolds number smaller than the critical Reynolds number. So far, there has been no report that the two different approaches yield different critical transition Reynolds numbers.

When the DNS method is used for the determination of the critical transition Reynolds number, one prefers to start from a turbulent flow. This is because that, if one starts from a laminar flow, the corresponding Reynolds number must be smaller than the critical transition Reynolds number, which will results in the introduction of a certain initial disturbance. However, since all small amplitude disturbances are stable according to the stability theory, once the DNS starts, the initial disturbances may die out before the transition is triggered. Thus, one has to try many different initial disturbances. If transition does take place, the corresponding Reynolds number has already been larger than the critical Reynolds number. Therefore, one has to lower the Reynolds number a little bit, and repeats the same process until the critical transition Reynolds number is found. Obviously, such a method is time consuming. Therefore, in our method, we start from a turbulent flow.

2 Governing equations and numerical method

The DNS is done under a Cartesian coordinate system. The origin is located at the channel center, and x^*, y^*, and z^* represent the coordinates in the stream-wise, wall-normal, and span-wise directions, respectively. The reference velocity U_CL^* is the laminar velocity at the channel center, and the length scale is the half channel width h^*. In the non-dimensional equation, the Reynolds number is defined by Re_L =U_CL^*h^*/υ^*. The superscript * represents the dimensional quantities in this paper. The non-dimensional disturbance equations are

(1)

where the subscript 0 implies the quantities of the laminar basic flow U₀ =1-y², the superscript' implies the disturbance quantities, Π' = p'+(u^'2+v^'2+w^'2) /2, and the nonlinear terms are

(2)

The pseudo-spectrum method is used in the stream-wise and span-wise directions. The functions are expanded as the Fourier series as follows:

(3)

where φ represents the disturbance vector, is the spectral component, and α and β are the fundamental wave numbers in the stream-wise and span-wise directions, respectively. The uniform grid and periodic boundary condition are employed in the stream-wise and span-wise directions. In the wall-normal direction, the following finite difference scheme is used:

(4)

which is of the 4th-order accuracy.

Non-uniform grids are applied in the wall-normal direction, i.e.,

(5)

The non-slip boundary condition is used on the top and bottom walls.

Then, Eq. (1) can be rewritten as follows:

(6)

where L₀ (φ) and L(φ) represent the linear terms, and N(φ) represents the nonlinear term. The time advancing scheme is

(7)

The mesh for our computation is 128×128×128 in the x-, y-, and z-directions, respectively. In the turbulent flow in Stage Ⅰ, there are 21 points below y⁺=30, and the grid point closest to the wall is at y⁺ =0.07. In the following stages, because the viscosity length υ^*/U₊^* increases when the computation Reynolds number Re_L decreases, the mesh becomes even denser in the wall-normal direction.

3 Method for the critical Reynolds number

As stated above, when the critical Reynolds number is solved, we start from a turbulent flow, and then decrease the Reynolds number gradually until the turbulent flow can no longer be maintained. To obtain the initial turbulent flow, we start from a flow for which the Reynolds number is higher than the critical Reynolds number for the instability of the laminar flow. By introducing three unstable Tollmien-Schlichting (T-S) waves as the initial disturbances, transition can readily be triggered to obtain a turbulent flow. Then, we reduce the Reynolds number step by step. At each step, we keep the DNS going until the flow becomes statistically steady. We assure that a turbulent flow can exist under the given Reynolds number only when a statistically steady state is reached.

To lower the Reynolds number, physically, we can either reduce the velocity U_CL^* of the corresponding laminar flow at the channel center with the channel half width h^* or increase the viscosity υ^*. Obviously, in experiments, the most natural choice is to change the reference velocity, because otherwise, one has to change the working fluid or the experimental facility.

However, at the end of each step, we have a statistically steady flow field. If we take the so obtained flow field as the initial flow field for the DNS of the next step, since the reference velocity has changed, it has to be rescaled as follows to obtain the velocity component under the new reference velocity. The velocity component u^' is first converted to its dimensional value by being multiplied by the reference velocity U_CL1^*, through which we can obtain u^'U_CL1^*. In the next computational step, the reference velocity is U_CL2^*, and its non-dimensional value is u^'U_CL1^*/U_CL2^*. The mean flow is treated similarly. Therefore, the total velocity is (u^'+1-y²) U_CL1^*/U_CL2^*. Since in the new computational step, we still use the non-dimensional disturbance equations, the initial disturbance velocity is (u^'+1-y²) U_CL1^*/U_CL2^*-(1-y²). The other components are treated similarly.

After the rescaling, the flow field does not correspond to a statistically steady flow field under the new Reynolds number. Therefore, it will take another long time to reach a new statistically steady state.

However, if we do not rescale the original flow field obtained in the previous stage as stated above, we will keep the non-dimensional flow field unchanged and use it as the initial condition for the next computational stage, for which the Reynolds number is different from the previous stage. Then, actually, although it is not a real flow field under the new reference velocity, it is actually close to a statistically steady state of the turbulent flow under the new Reynolds number. The reasons are as follows:

(ⅰ) For turbulent channel flows, the non-dimensional mean flow profile bears a high degree of similarity when the Reynolds number changes, i.e., they are close to each other in the non-dimensional form.

(ⅱ) The distribution of the averaged turbulent quantities, such as the turbulent energy and the Reynolds stresses, also bears certain similarity in their non-dimensional form. This is true even for supersonic boundary layers ^[6].

Therefore, the non-dimensional mean flow profile and turbulent quantities obtained in the end of the previous stage can be used directly as the initial condition in the next stage, for which the Reynolds number is different. Since the non-dimensional basic flow profile is unchanged when the Reynolds number changes, the mean flow modification caused by the disturbance velocity is also unchanged. In such a way, the transient period for reaching a new statistically steady state should be very short, which will be confirmed in the next section.

4 Results and discussion

As discussed in the first section, we start with a turbulent flow, and then decrease the Reynolds number gradually until the turbulent flow can no longer be maintained.

To obtain a turbulent flow, we start with a laminar flow with the Reynolds number Re_L =8 000, which is appreciably larger than Re_cr1, so that unstable waves for triggering the transition can readily be found. The process of lowering the Reynolds number is shown in Table 1. In Stage Ⅰ, a turbulent flow is obtained. Then, the Reynolds number is reduced step by step, while the flow keeps being turbulent until in Stage Ⅵ. If the Reynolds number is reduced further as in Stage Ⅶ, even slightly, the flow will decay to be laminar.

The unstable waves for triggering the transition at the beginning of Stage Ⅰ consist of three large amplitude T-S waves, the parameters of which are the same as those of Case 3 in Ref. [7]. Figures 2 and 3 show the disturbance energy e₂ and the mean velocity gradient on the wall du/dy in Stage Ⅰ. The quick growth of e₂ and du/dy implies that the transition process takes place immediately, while the subsequent decay of e₂ and du/dy implies that the transition process is ended at approximately t=45. However, it is not statistically steady yet, because du/dy is still decaying at the end of Stage Ⅰ, and its value for a statistically steady state is 2.

The transition in this case is a bypass transition, because the amplitudes of the imposed initial disturbances are quite large (see Ref. [7]). As shown in Ref. [5], during the breakdown process of the laminar-turbulent transition, the stability characteristics of the mean flow profile undergo a profound change, i.e., the unstable zone encircled by the neutral curve in the frequency-wavenumber space for the T-S waves becomes larger and larger, so does the maximum amplification rate of the T-S waves. Therefore, a positive feedback mechanism for producing more unstable waves is formed, and the flow becomes more and more chaotic, accompanied by the quick increases in the turbulent energy and wall friction coefficient. Actually, at the end of the transition process, the positive feedback mechanism quickly disappears, but the turbulent energy and wall friction coefficient are larger than those for a fully developed turbulent flow. In our case, it is also true. The neutral curves of the mean flow profiles during the transition process are plotted in Fig. 4. It is shown that the unstable zone enlarges very quickly from t=0 to t=1, keeps much larger than that of the laminar flow until t=35, and then shrinks quickly. At t=40, it is very small, though it is still a bit larger than that of the laminar flow.

Figure 5 shows the mean velocity profile at the end of Stage Ⅰ. The profile looks like a turbulent profile, but it is not the final statistically equilibrium profile. Figure 6 shows the mean velocity distribution u⁺ near the wall. It reflects the typical characteristics of a turbulent flow.

Stage Ⅱ starts from t=230, and Re_L in the computational program is reduced drastically to 2 590. At t=3 500, the flow remains to be turbulent, and becomes statistically equilibrium. The consumed time is quite long. In Stage Ⅲ, Re_L is further decreased to 2 500, and the flow remains to be turbulent and reaches its statistical equilibrium state quite soon.

As stated in the previous section, we expect that the transient period for reaching a new statistically steady state can be very short by use of our method with the initial condition. This is true. As shown in Figs. 7 and 8, the disturbance energy and velocity gradient at the wall remain almost unchanged under the sudden change of the Reynolds number at the beginning of Stage Ⅲ, which shows that they approach the new statistically steady state almost immediately.

Moreover, at the beginning of Stage Ⅱ, the flow field obtained in Stage Ⅰ is still far from being statistically steady. In other words, it does not correspond to a fully developed turbulent flow. Therefore, the above similarity argument does not hold in this case, and the turbulent energy and mean velocity gradient at the wall do have an apparent surge. It is a long period for the flow to become fully developed (see Figs. 7 and 8).

At the end of Stages Ⅱ and Ⅲ, the flow profiles correspond to fully developed turbulent profiles (see Figs. 9 and 10). The almost exact coincidence of the two curves is a manifestation of the similarity property as mentioned above.

In Stage Ⅳ, since the Reynolds number is not too far from the critical Reynolds number, the Re_L in our computation is decreased to several different values. For each of them, the flow field is calculated from t=7 000 up to t=9 740. The flows with Re_L ≥ 2 000 remain to be turbulent, while the others, i.e., cases with Re_L < 2 000, eventually become laminar (see Figs. 11-14). From Fig. 11, we can see the evolution of the respective disturbance energy of the flow. If it decays to zero, the flow becomes laminar, while if it tends to a finite value, the flow remains to be turbulent. The mean velocity gradient at the wall, the mean velocity profile, and the mean velocity distribution in the wall region are shown, respectively, in Figs. 12, 13, and 14, from which we can judge if the corresponding flow is laminar or turbulent in a similar way.

Since the Reynolds number 2 000 is even closer to the critical Reynolds number, in Stages Ⅴ to Ⅶ, we search the critical Reynolds number more carefully. Since the basic procedures are similar to those for Stage Ⅳ, the details are omitted here. In Stage Ⅶ, the flow is laminar. Therefore, the critical Reynolds number must be between 1 700 and 1 650.

In experiments, two Reynolds numbers are often used, i.e., Re and Re_b, for which Re=U_c^*h^*/υ^* and Re_b =1.5U_m^*h^*/υ^*, where U_c^* is the turbulent mean velocity at the channel center, and U_m^* is the mean velocity of the turbulent flow profile. Table 2 shows the corresponding Re and Re_b at the end of each stage. They are 1 081 and 1 233, respectively, at the end of Stage Ⅵ. Therefore, the critical transition Reynolds number for the plane channel flow determined in this paper is nearly Re_cr2 =1 081 or R_bcr2 =1 233, which is very close to those obtained by the experiments ^[8-11].

We have checked if our computational domain is large enough to ensure the correctness of our DNS in Table 3. The results show that, for turbulent flows, the stream-wise and span-wise domains should be large enough compared with the average length and average spacing of low speed streaks in the wall region of the turbulent flow ^[12]. We also checked the results by doubling the stream-wise and span-wise extents of our computational domain in each computation stage. The results are satisfactory.

5 Conclusions

One can find the critical transition Reynolds number by the DNS method. The way we choose is to start with a turbulent flow, the Reynolds number of which is higher than the critical transition Reynolds number. Then, the value of the Reynolds number is decreased step by step until the flow can no longer remain turbulent. The lowest value for maintaining a turbulent flow must be close to the critical transition Reynolds number.

For plane channel flows, the obtained critical Reynolds number is Re_cr2=1 081, which is based on the mean velocity at the center of the mean turbulent flow profile, or Re_bcr2 =1 233, which is based on 1.5 times of the mean velocity of the turbulent flow profile. The results agree well with the experimental data, implying that the searching method employed in this paper is appropriate.

Acknowledgements The author is grateful to Professor Heng ZHOU in Tianjin University for his constant encouragement and helpful suggestions.