1 Introduction

The transition prediction in a boundary layer is an important problem in aircraft design. It has received much attention in the past decades due to its key role in drag reduction and heat protection. In the natural transition scenario of a boundary layer, the linear growth of an infinitesimal disturbance is a fundamental stage. Therefore, the investigation of the linear evolution of a small disturbance plays a significant role in transition predictions.

The e ^N method, proposed in aircraft industry, has been widely used to predict the transition onset of the boundary layers from laminar to turbulent flow. Although its usefulness has been well validated in two-dimensional and axially symmetric boundary layers, the extension of this method to general three-dimensional boundary layers remains to be solved. The e ^N method is a semi-empirical method based on the linear stability theory (LST). In this method, the LST analysis ^[1] is used to calculate the logarithmic amplification ratio of the amplitude versus its initial value for each wave. The envelope of these ratios, i.e., the N factor which represents the amplitude evolution of the most amplified disturbance, is used to predict the transition according to the N criterion which has been established by transition experiments ^[2-3].

For a two-dimensional flow, the frequency and the spanwise wavenumber of a disturbance remain unchanged when it propagates downstream. Then, the wavenumber and the growth rate in the streamwise direction can be solved according to the dispersion relationship. After that, the amplitude logarithmic amplification ratio of each wave with different frequencies and/or spanwise wavenumbers and the N factor can be calculated by integrating the growth rate along the x-direction. For a three-dimensional swept wing boundary layer with infinite spanwise length, the flow is uniform in the spanwise direction. Therefore, the spanwise wavenumber of the disturbance is invariant when it propagates downstream, and the same transition prediction method can be used as in two-dimensional boundary layers ^[4-5].

However, for general three-dimensional boundary layers, the propagation direction of a disturbance is not always evident. Moreover, the flow is nonuniform in both the streamwise direction and the spanwise direction. Therefore, both the streamwise complex wavenumber and the spanwise complex wavenumber vary slowly when the disturbance propagates downstream, i.e., there are two unknowns for the equation given by the dispersion relationship. The solution is not unique unless the supplement conditions are specified. In contrast to the situation for two-dimensional flows, there are two additional problems that need to be dealt with when the e ^N method is extended to general three-dimensional boundary layers. The first problem is how to choose an appropriate path along which the growth rate is integrated to obtain the N factor. It seems to have reached a consensus. Extensive theoretical and computational studies on this subject ^[6-9] indicate that the group velocity line of the disturbance is a reasonable and appropriate integration path. The second problem is how to establish supplement conditions so as to make the equation given by the dispersion relationship well-posed. On this subject, Nayfeh ^[6] and Cebeci and Stewartson ^[7] theoretically analyzed the stability of three-dimensional boundary layers. With the group velocity concepts which led to the requirement , Cebeci and Stewartson ^[7] actually established the saddle point method (SPM), which could determine the far-field asymptotic behavior of complex wavenumbers. It seems that the SPM has been the only popular approach to extend the e ^N method to general three-dimensional boundary layers.

Although the e ^N method based on the SPM, called the e ^N-SPM in the present paper, has been widely used to predict the transition in general three-dimensional boundary layers ^[10-12], this method has two inherent defects. First, from the theoretical point of view, the SPM is only suitable for evaluating the far-field asymptotic evolution of a wave packet with a fixed frequency. However, for practical applications of the e ^N method, the near-field local evolution of a single wave should always be followed. In the e ^N-SPM method, one needs to, at each location along the propagation path, find out the wave that not only satisfies the condition of the saddle point but also has the largest growth rate ^[7]. However, the waves determined by the e ^N-SPM method on the propagation path do not coincide with the actual evolution process of a disturbance. Second, the SPM has an implicit assumption that the dispersion relationship is invariant during the evolution of a disturbance. However, due to the nonuniform flow in general three-dimensional boundary layers, the dispersion relationship gradually varies along both the streamwise direction and the spanwise direction, which aggravates the first defect in turn. Due to the above two inherent defects, the amplitude of a disturbance predicted by the e ^N-SPM, sometimes, obviously deviates from its actual one.

Two kinds of supplement conditions used to close the equation given by the dispersion relationship are proposed in the present paper by introducing the compatibility relationship (CR) and the general theory of ray tracing (RT) which have been used to describe the near-field local evolution of acoustic waves ^[13-14]. Then, two new e ^N methods based on the CR and RT, i.e., the e ^N-CR method and the e ^N-RT method in this paper, are established to predict the evolution of the disturbances in general three-dimensional boundary layers. The CR provides a constraint condition between the two wavenumber components of an acoustic wave. In fact, Chang ^[9] has used the CR in the linear parabolized stability equation (LPSE) approach to predict the linear evolution of the disturbances in three-dimensional boundary layers. The RT describes the evolution of the wavenumbers along and only along the group velocity line in anisotropic systems ^[15]. Fortunately, the group velocity line is just the integration path along which the N factor can be obtained. To the best knowledge of the authors, the RT has not been used in the subject of extending the e ^N method to general three-dimensional boundary layers.

The rest of the present paper is organized as follows. In Section 2, we formulate the dispersion relationship for a general three-dimensional boundary layer, and theoretically analyze the relation among the imaginary parts of the wave parameters. How to establish the two new e ^N methods is also given in this section. The results and comparisons with the direct numerical simulation (DNS) and the e ^N-SPM are presented in Section 3. A summary and some concluding remarks are given in Section 4.

2 Physical problem and methodology 2.1 Dispersion relationship for disturbances in general three-dimensional boundary layers

For two-dimensional boundary layers, such as flat plate boundary layers and axially symmetric cone boundary layers, the flow field varies slowly in the streamwise direction. With the LST, one can adopt a quasi-parallel hypothesis in the streamwise direction, and then derive an Orr-Sommerfeld equation. For the general three-dimensional boundary layer considered in the present study, the flow field gradually varies in both the streamwise direction and the spanwise direction. Therefore, the quasi-parallel hypothesis needs to be adopted in both the streamwise direction and the spanwise direction. A disturbance in the general three-dimensional boundary layer can be presented as the following form of travelling waves:

(1)

where x, y, and z represent the coordinates in the streamwise, wall-normal, and spanwise directions, respectively, and t denotes the time. φ'=( ρ', u', v', w', T') stands for the disturbance in the boundary layer, ρ' is the density fluctuation, T' is the temperature fluctuation, and u', v', and w' represent the velocity fluctuations in the streamwise, wall-normal, and spanwise directions, respectively. A₀ is the amplitude of u' at the initial location ( x₀, z₀), is the shape function normalized by A₀, and c.c. represents the conjugate counterpart. The complex θ is concerned with the wavenumbers and growth rates of the disturbance. α and β are complex. Their real parts ( α_r and β_r) represent the streamwise wavenumber and the spanwise wavenumber, respectively, and the negative imaginary parts (- α_i and - β_i) represent the growth rates in the streamwise direction and the spanwise direction, respectively. ω_r is the frequency of the disturbance, and ω_i represents the temporal growth rate. From Eq. (1), we have

(2)

This implies that the three wave parameters ω, α, and β are not mutually independent, but connect with each other through the complex phase θ. The three parameters all vary slowly with respect to x and z.

Substitute Eqs. (1) and (2) into the linearized disturbance equations in the three-dimensional boundary layer. Then, we can derive the familiar Orr-Sommerfeld equation for a general three-dimensional boundary layer flow. Combined with the homogeneous boundary conditions, we can establish the dispersion relationship for the wave parameters as follows:

(3)

The flow field is nonuniform, and thus the dispersion relationship ω explicitly depends on x and z. For a fixed location in the xz-plane, if two of the three parameters ω, α, and β are given, we can get the third one by solving the dispersion relationship.

Since the spatial mode, where ω is real, is unstable, the local amplitude of a disturbance can be written as follows:

(4)

where N, i.e., the N factor, represents the logarithmic amplification ratio of the local amplitude versus the initial one. The key problem for the prediction of the evolution of the disturbance is to obtain the N factor with Eq. (4).

Now, we consider the evolution of the wave parameters of the disturbance in the three-dimensional boundary layer. Apparently, the wave parameters are invariant with respect to time, which yields

From Eq. (2), we have

(5a)

(5b)

(5c)

It suggests that the complex frequency does not change with x and z. In other words, its frequency remains unchanged when the disturbance evolves in the steady boundary layer.

For two-or three-dimensional flows being uniform in the spanwise direction, leads to . Consequently, the frequency and spanwise wavenumber do not change with x. For an unstable wave with the specific frequency ω and the complex spanwise wavenumber β, the complex streamwise wavenumber α can be solved through Eq. (3) at any location downstream. Then, the N factor can be obtained according to Eq. (4).

However, for general three-dimensional boundary layers, , which means that , i.e., β keeps changing when the disturbance propagates downstream. With the two unknowns α and β, Eq. (3) cannot yield the unique solution. Therefore, additional relationships between α and β are needed to make Eq. (3) well-posed.

2.2 Conservation law of imaginary parts for wave parameters

The CR and RT are originally utilized in conservative systems, e.g., water waves and acoustic waves, in which the wave parameters are real. However, the boundary layers belong to nonconservative systems, in which the wave parameters are complex. In this paper, the CR and RT are only used to describe the behavior of the real parts of the wave parameters. Therefore, we must additionally consider the imaginary parts of the wave parameters alone.

In practical general three-dimensional boundary layers, the imaginary parts of the wave parameters are much less than the real parts. This is a prerequisite not only in this subsection but also in Subsection 2.4.

Considering the dispersion relationship (3), we define the group velocity in the general three-dimensional boundary layer as follows:

(6)

The group velocity u_g is generally complex, but its imaginary part is extremely small in practice. Therefore, we define the magnitude of the group velocity as follows:

The direction l is presented by the direction angle φ=arctan( g _z/ g _x), where

At a specific location ( x, z), the total differential of Eq. (3) is

(7)

which specifies the variation of ω caused by the variations of α and β. Keep the real parts of α and β constant, and change their imaginary parts. Then, the imaginary part of Eq. (7) can be derived as follows:

(8)

For a fixed wave vector ( α_r, β_r), by integrating Eq. (8) and neglecting the high order small quantities introduced by the tiny changes of g _x and g _z, we have

(9)

where

in which φ is the direction angle of the group velocity. Equation (9) indicates that a conservation law exists among the imaginary parts of the wave parameters, which is called the conservation law of the imaginary parts for wave parameters (CLIP). In fact, the left side of Eq. (9) is precisely the growth rate of a disturbance along the group velocity direction. The CLIP suggests that the growth rate along the group velocity direction only depends on the real parts of the wave parameters, and is independent of the imaginary parts. For an unstable wave of the spatial mode, the frequency is real. Therefore, the growth rate along the group velocity direction is

(10)

According to the CLIP, the value of β_i does not affect the growth rate σ_g. For simplicity, β_i is set to be zero when the N factor is calculated by Eqs. (3) and (4). Consequently, the first supplement condition that we establish is

(11)

and the N factor can be acquired by the new integration as follows:

(12)

where the integration path is along the group velocity line. For the specific wave vector ( α_r, β_r), the group velocity direction, i.e., the directional angle φ, can be easily obtained. Next, we will establish the relationship between α_r and β_r to determine the evolutions of α_r and β_r.

2.3 Compatibility relation application in general three-dimensional boundary layers

Equation (5a) describes the relation between α and β in the flow field. Here, we just use their real parts, i.e.,

(13)

as the compatibility relationship of the wave vector. The CR actually provides a constraint condition for the variation of the wavenumbers caused by the nonuniform flow field.

Let Eqs. (11) and (13) be the supplement conditions for the dispersion relationship (3). Then, a new e ^N method based on the CR, which we call the e ^N-CR method, is established, with which we can predict the evolution of the disturbance in the general three-dimensional boundary layer. The computing strategy is as follows:

(ⅰ) Specify the frequency ω that we concern, let β_i=0 in the total process, and initialize the spanwise distribution of β_r at an initial location upstream.

(ⅱ) With known ω and β_r, α can be solved by Eq. (3) along the spanwise curve. Meanwhile, the partial derivative and hence the partial derivative can be obtained.

(ⅲ) Integrate the partial derivative from the present location to the next location downstream to obtain the β_r at the next location.

(ⅳ) March to the next location and repeat Steps (ⅱ) and (ⅲ) until α and β are solved in the entire computation zone.

(ⅴ) Integrate the growth rate along the group velocity line according to Eq. (12), and acquire the N factor.

2.4 Extension of the ray tracing theory to general three-dimensional boundary layers

Take the partial derivatives of the dispersion relationship (3) with respect to x and z. Then, we have

(14)

Substitute Eq. (5) into Eq. (14). Then, we have

(15a)

(15b)

Since the imaginary parts of u_g, α, and β are all small quantities, the real parts of Eq. (15a) and (15b) can be derived as follows:

(16a)

(16b)

Here, we neglect some second-order small quantities caused by the products of the imaginary parts. Equations (16a) and (16b) indicate that the evolution of the (streamwise/spanwise) wavenumber along the group velocity direction is closely connected with the inhomogeneity of the flow. This inhomogeneity is reflected by the partial derivatives of the dispersion relation.

Since Eq. (16a) is not independent of Eq. (16b), only one of them shall be utilized. For the transition prediction in the general three-dimensional boundary layer, Eq. (16b) is generally chosen because that the streamwise direction is usually close to the group velocity direction. From Eq. (16b), we can obtain the directional derivative of β_r along the group velocity direction, i.e.,

(17)

Combining Eqs. (11) and (17) with Eqs. (3) and (12), we can establish another new e ^N method based on the RT, i.e., the e ^N-RT method. Then, we can predict the evolution of the disturbance along the group velocity line in the general three-dimensional boundary layer. The computing strategy is as follows:

(ⅰ) Specify the frequency ω that we concern, let β_i=0 in the total process, and initialize the β_r at an initial location ( x₀, z₀).

(ⅱ) With known ω and β_r, α can be solved from Eq. (3). With the specific α and β, solve U_g and . Then, we can obtain the directional derivative

(ⅲ) Integrate the directional derivative along the group velocity line. Then, β_r can be obtained at the next location.

(ⅳ) March to the next location, and repeat Steps (ⅱ) and (ⅲ) until α and β are solved along the hole group velocity line.

(ⅴ) Integrate the growth rate along the group velocity line according to Eq. (12), and acquire the N factor.

3 Results and discussion

In this section, the reliability of the proposed e ^N methods is verified and validated by performing the DNS in a real general three-dimensional (rather than a quasi three-dimensional) boundary layer. The evolution of the disturbances is calculated by the DNS and the proposed e ^N methods, and the computing results are compared. Meanwhile, the results are also compared with those obtained by the e ^N-SPM method which is recently the most popular transition prediction method for general three-dimensional boundary layers.

3.1 Model and computation domain

An aircraft model consisting of half of a cone and a flat plate (see Fig. 1) is investigated here. The flight altitude is 30 km, and the corresponding atmospheric parameters are considered as the reference scales. The reference length is 2.06 mm, and the flight Mach number is 10. The physical variables and parameters are non-dimensionalized by the corresponding reference quantities. The corresponding Reynolds number is 7 762. Hereinafter, all the physical variables and parameters presented in the results are non-dimensional.

The quadrilateral computational domain is chosen on the side face of the cone, which is outlined by the thick solid lines in Fig. 1. In Fig. 1, the thin lines indicate the group velocity direction of the disturbance with the spanwise wavenumber β=0. In fact, for disturbances with different spanwise wavenumbers, the group velocity directions are almost identical. They all align with the potential flow direction. Since the disturbance propagates along the group velocity direction, the computational domain is chosen so that the side boundaries are nearly parallel to the group velocity lines. The inlet boundary is chosen behind the locations where the disturbance begins to amplify. In the following calculations, the x-axis is set to be aligned with the potential flow direction, the z-axis is perpendicular to the x-axis, and both the x-and z-axes are along the surface.

3.2 Disturbance evolution in general three-dimensional boundary layers

In order to calculate the evolution of a single wave, the disturbance is introduced at the inlet of the computation domain as follows:

(18)

where β_r is independent of z at the inlet. The distribution of the initial amplitude A₀( z) is illustrated in Fig. 2. The initial amplitude vanishes on the side boundaries so that the zero boundary conditions can be specified in the DNS calculation. The wave packet just propagates along the group velocity direction, and it will not extend its spanwise extent obviously. Because the spanwise boundaries of the DNS domain are almost parallel to the group velocity, it is appropriate to specify the zero boundary conditions at the spanwise boundaries in the DNS. The initial amplitude is uniform in most parts of the spanwise range (see Fig. 2), which makes the calculation of the DNS consistent with the situation of a single wave to some extent. In fact, the initial amplitude A₀( z) does contain some components with specific spanwise wave numbers, but the spanwise wave numbers of the main components are concentrated in an extremely small neighborhood around zero. Therefore, the DNS calculation can well describe the evolution of a single wave. This has been validated in the early work by Liu ^[16]. represents the local eigenfunction obtained by the LST. Once introduced into the domain and calculated by the DNS, such disturbances will experience a slight modulation process near the inlet.

The disturbance evolution with the frequency ω=1.2 and the spanwise wavenumber β_r=0 is calculated by the DNS. The peak value of A₀ is 1× 10^-5. Therefore, the nonlinear interaction is absent in the DNS calculation. Four theoretical methods are used to calculate the evolution of the disturbance along the group velocity line. Two of them are the proposed methods in this paper, i.e., the e ^N-CR method and the e ^N-RT method. The third one is the e ^N-SPM method, and the last one, termed as the e ^N-fixed β method, is an e ^N method in which the spanwise wavenumber remains unchanged. The results obtained by the four methods are compared with those obtained by the DNS. Figure 3(a) illustrates the changes of the spanwise wavenumber predicted by the different methods. The comparison of the amplitude magnitude is presented in Fig. 3(b). Since the disturbance experiences a slight modulation process, the disturbance amplitude is normalized at the location x=28 for comparison. As shown in Fig. 3(a), the DNS result suggests that the spanwise wavenumber increases with the propagation. Therefore, it is improper to keep β unchanged during the calculation. Meanwhile, the amplitude predicted by the e ^N-fixed β method is larger than that obtained by the DNS. The comparison also shows that both the e ^N-CR method and the e ^N-RT method can well predict the evolution of the spanwise wavenumber. Therefore, the amplitudes predicted by them also essentially agree with that predicted by the DNS (see Fig. 3(b)). However, as shown in Fig. 3(a), the spanwise wavenumber determined by the e ^N-SPM method severely deviates from the result determined by the DNS. This is attributed to the fact that the e ^N-SPM method considers the most unstable wave at each location along the propagation path of the disturbance. Therefore, this method extremely overestimates the amplitude of the disturbance compared with the DNS (see Fig. 3(b)).

As shown in Fig. 3(a), the result obtained by the e ^N-SPM method suggests that the spanwise wavenumber of the most unstable wave is about -0.5 at the inlet. Therefore, similar calculations are carried out for the introduced disturbance with the same frequency ω=1.2 but a different spanwise wavenumber β=-0.5. Figures 4(a) and 4(b) present the evolutions of the spanwise wavenumber and the disturbance amplitude, respectively. In Fig. 4(b), the disturbance amplitude is normalized at the location x=28. Figure 4(a) suggests that, both the e ^N-CR method and the e ^N-RT method can well predict the change of the spanwise wavenumber, and thus the amplitude evolutions predicted by them are most approximate to that obtained by the DNS. As can be seen in Fig. 4(a), the difference between the spanwise wavenumbers predicted by the e ^N-SPM method and the DNS is small in the upstream. Therefore, the difference between the amplitudes predicted by the e ^N-SPM method and the DNS slightly decreases (see Fig. 4(b)). However, overall, the e ^N-SPM method still overestimates the amplitude of the disturbance. In this case, it is different from the situation in Fig. 3(b) that the e ^N-fixed β method underestimates the disturbance amplitude.

Another disturbance with the spanwise wavenumber β=-0.9 at the inlet is also calculated in a similar way. The disturbance has the same frequency, i.e., ω=1.2. Figures 5(a) and 5(b) show the evolutions of the spanwise wavenumber and the disturbance amplitude. The disturbance amplitude is normalized at the location x=20. As suggested in the figures, both the e ^N-CR method and the e ^N-RT method well predict the evolution of the spanwise wavenumber, and thus well predict the evolution of the disturbance amplitude. However, the wavenumber determined by the e ^N-SPM method severely deviates from the result given by the DNS, and hence it extremely overestimates the disturbance amplitude compared with the DNS. In this case, the amplitude predicted by the e ^N-fixed β method is much less than that predicted by the DNS, and the unstable region predicted by the e ^N-fixed β method is also much smaller than that predicted by the DNS.

The results presented in Figs. 3-5 indicate that both the e ^N-fixed β method and the e ^N-SPM method cannot correctly describe the change of the spanwise wavenumber, and hence neither of them can correctly predict the evolution of the disturbance amplitude. The e ^N-fixed β method may overestimate or underestimate the disturbance amplitude depending on the specific spanwise wavenumber. The e ^N-SPM method always overestimates the amplitude of the disturbance. In the calculation of the e ^N-SPM method, it picks out the spanwise wavenumbers which have the largest growth rates at each location along the propagation path of the disturbance, and the N factor is calculated by integrating these largest growth rates. The result predicted by the e ^N-SPM method is not the evolution of any disturbance. Because the disturbances with different spanwise wavenumbers have almost an identical group velocity direction, the result given by the e ^N-SPM method only depends on the base flow, and has nothing to do with the initial spanwise wavenumber. Therefore, the e ^N-SPM method provides the identical prediction in Figs. 3-5.

By contrast, both the e ^N-CR method and the e ^N-RT method can well describe the evolution of the spanwise wavenumber for the disturbances with different spanwise wavenumbers. Therefore, these two proposed methods can well predict the amplification of the disturbance amplitude. This is attributed to two facts. First, different from the SPM suitable for evaluating the far-field asymptotic evolution of a wave packet, what the CR and/or RT describe is the near-field local evolution of a single wave. In the e ^N method, one just needs to follow the local evolution of a single wave. Second, the CR and RT do not need the assumption that the dispersion is invariant, i.e., the flow field is uniform, which is implicitly contained in the SPM. In fact, the CR provides a constraint condition for the variation of the wavenumbers caused by the nonuniform flow field, and the RT describes the relationship between the evolutions of the wavenumbers and the inhomogeneity of the flow field.

In the prediction for the disturbance evolution by the e ^N-CR method, the spanwise wavenumber β_r needs to be specified along the inlet, and the calculation must be performed in the whole domain. For the e ^N-RT method, we just give the initial spanwise wavenumber β_r at the point that we concern. Accordingly, the calculation just needs to be performed along a group velocity line. Therefore, it should be more convenient to use the e ^N-RT method for the transition prediction in engineering.

4 Conclusions

The e ^N method has made great achievements in the transition prediction for two-dimensional and axially symmetric boundary layers. In the present study, in order to extend the e ^N method to a general three-dimensional boundary layer, based on the fact that the imaginary parts of the wave parameters and the dispersion relationship are small quantities, the CLIP is deduced, and provides a supplement condition with the simple form of β_i=0. The CR and RT are used in the general three-dimensional boundary layer. The CR describes the relationship between α and β in the entire domain, and the RT describes the evolutions of β and/or α along the group velocity direction. Combining the dispersion relationship with the CR and RT, we propose two new e ^N methods, i.e., the e ^N-CR method and the e ^N-RT method, to predict the evolution of the disturbance in the general three-dimensional boundary layer. Their reliabilities and superiorities are validated by a DNS in the general three-dimensional hypersonic boundary layer. Compared with the DNS results, both the e ^N-CR method and the e ^N-RT method can well predict the evolutions of the spanwise wavenumber and the amplitude of the disturbance, indicating that the proposed e ^N methods have great potentials in improving the transition prediction accuracy in general three-dimensional boundary layers. Most notably, this is the first study, to our knowledge, to find out the CLIP and apply the RT to the transition prediction in general three-dimensional boundary layers.

Acknowledgements The authors are grateful to Professor Xuesong WU of Imperial College London for valuable discussion and suggestions.