(X, Y, Z)	incoming-flow coordinate system;	Re,	Reynolds number;
(x, y, z)	chordwise coordinate system;	Rec,	chord-length Reynolds number;
(xs, ys, z)	body-fitted coordinate system;	α,	eigenvalue;
(xt, ys, zt)	crossflow coordinate system;	β,	spanwise wavenumber;
u0	streamwise velocity;	ω,	frequency;
v0	normal velocity;	σ,	growth rate;
w0	spanwise velocity;	ã,	complex streamwise wavenumber;
wt	crossflow velocity;	Ã	, complex amplitude;
U X	incoming-flow-direction velocity;	A,	amplitude.

1 Introduction

Hydrodynamic stability and laminar-turbulent transition are important topics in fluid mechanics and aerospace engineering ^[1]. For civil aircraft, the most common configuration is the swept wing ^[2]. Maintaining a larger laminar flow area on the wing can greatly reduce friction and heat transfer. The design of laminar wing requires a transition prediction tool of high efficiency and high precision. Currently, the most widely used prediction tool is the semi-empirical e ^N method ^[3]. The e ^N method can be used either with a local stability approach such as the linear stability theory (LST) or a streamwise marching stability approach such as the linear parabolized stability equation (LPSE).

The LST is a local eigenvalue problem based on the local parallel flow hypothesis ^[4]. The advantage of the LST is that it is a local problem. Therefore, it can be conveniently applied to the stability analyses and transition predictions for boundary layers. The disadvantage of the LST is that it ignores the influence of non-parallelism on the instability. In weakly non-parallel boundary layers ^[5], the growth rate predicted by the LST is almost the same as that of the direct numerical simulation (DNS). However, there is a considerable discrepancy between them in their predictions regarding strongly non-parallel boundary layers ^[6].

The LPSE is a parabolic partial differential equation, and can be solved by a marching procedure in the streamwise direction. It considers the non-parallelism of the base flow and the shape function ^[7-8]. The LPSE can predict the growth rates accurately, even for strongly non-parallel boundary layers ^[9-10], and the computational burden is much smaller than that of the DNS. However, compared with the LST (local approach), the LPSE is inflexible and inefficient when it is applied to e ^N calculations. To obtain the envelope of the N factor, the LPSE must integrate all unstable disturbances with different frequencies and spanwise wavenumbers, whereas for the LST, the maximum N factor can be obtained simply by integrating the local maximum growth rate ^[11].

The expansion of the parabolized stability equation (EPSE) is obtained by the Taylor expansion of the LPSE in the streamwise direction. The EPSE considers the non-parallelism of the boundary layers, similar to the LPSE. The EPSE together with the homogeneous boundary conditions forms a local eigenvalue problem, similar to the LST. Therefore, the EPSE can predict growth rates more accurately than the LST, and can predict e ^N more efficiently and flexibly than the LPSE. The EPSE was first proposed by Bertolotti and Herbert ^[12], who established a form of the first-order EPSE (EPSE1) for two-dimensional (2D) incompressible flows to provide more accurate inlet conditions for the LPSE. They used the EPSE to study the Tollmien-Schlichting (T-S) instability of the incompressible flat boundary layer, and obtained results closer to those obtained from the LPSE than the LST. Yu et al. ^[13] presented a form of the EPSE1 for three-dimensional (3D) compressible flows, and used it to study the stability of the supersonic flat boundary layer and mode exchange. The obtained results were in agreement with those obtained from the LPSE and DNS. Huang and Wu ^[14] derived a similar equation by the Taylor expansion of the disturbance equation. They referred to their equation as the extended eigenvalue equation (EEV), and extended the equation to higher orders. The first-order EEV (EEV1) has the same form as the EPSE1. They studied the stability of the incompressible flat boundary layer, and found that the results of the second-order EEV (EEV2) were closer to the results of the DNS than the EEV1.

Swept-wing boundary layers are typical 3D boundary layers with high non-parallelism degree. There are four types of instabilities for swept-wing flows ^[15-16], i.e., attachment line, streamwise, centrifugal, and crossflow. The crossflow instability is considered to be the most "dangerous" instability for swept-wing boundary layers. Because there are inflectional points in the crossflow profiles, the crossflow instability is a type of inflection instabilities ^[17]. It has two instability modes, i.e., stationary crossflow vortices and unsteady traveling waves ^[18]. Although Bippes and Nitschke-Kowsky ^[19] have shown that the transition criteria based on the LST might be limited, transitions are commonly induced by the stationary crossflow vortices in the low-turbulence environments such as flight. The main tools for the study of the crossflow instability are the LST, PSE, and DNS. The LST can provide the stability characteristics, but the obtained growth rates are less accurate than the LPSE in common cases. Haynes and Reed ^[20] studied the crossflow instability of the NLF(2) -0415 incompressible swept-wing, and found that the N factor predicted by the LPSE was larger than that predicted by the LST by approximately 20%. Although the PSE and DNS can predict more accurate growth rates, it is difficult to use them to obtain the transition predictions for engineering purposes because of their computational complexity. The EPSE can overcome the shortcomings of the LST and LPSE. One can therefore expect that the EPSE would be a more accurate tool than the LST and a more efficient tool than the LPSE in swept-wing transition predictions.

In this paper, the general forms of the nth-order EPSE (EPSE n) (n=1, 2, …) are formulated. The EPSE1 and EPSE2 are applied to the crossflow instability analysis and the transition prediction of an infinite-span swept wing. We take the NLF(2) -0415 airfoil as the research model, and select the airliner cruise parameters as the calculation conditions, i.e., the Mach number 0.8 and the gas parameters at 10 km altitude. The results obtained from the EPSE n (n=1, 2) are compared with those obtained from the LST, LPSE, and DNS.

2 Formulation 2.1 LPSE

To formulate the EPSE, we first express the LPSE. We consider a wall with curvature, and apply the body-fitted coordinates (x_s, y_s, z) to it. x_s and z are tangent to the wall and perpendicular to each other, x_s is the streamwise direction, y_s is normal to the wall, and z is the spanwise direction. The streamwise Lame coefficient is H=1+ Ky, where K is the wall curvature along the streamwise direction.

The instantaneous quantity φ=(ρ, u, v, w, T) can be divided into the base flow quantity φ₀ and the disturbance quantity φ', i.e.,

(1)

Substitute Eq. (1) into the full Naiver-Stokes (N-S) equations, and subtracte the N-S equations φ₀ that are satisfied. Then, we obtain the disturbance equation. The disturbances are assumed to be small. Therefore, the equation can be linearized, and be written as follows:

(2)

where Γ, M, B, C, D, V _xx, V _xy, V _xz, V _yy, V _yz, and V _zz are the coefficient matrices ^[21], the elements of which are functions of the base flow, the curvature, the Mach number, and the Reynolds number. We express φ' as follows:

(3)

where , α, β, and ω are the shape function, the streamwise wavenumber, the spanwise wavenumber, and the frequency, respectively. Substituting Eq. (3) into Eq. (2) and ignoring and , we obtain the LPSE as follows:

(4)

where

For convenience, we rewrite Eq. (4) as follows:

(5)

where

2.2 General form of EPSE n

We set x_s= x_s0+, where is a small parameter. Let α be a constant in this neighborhood. The Taylor expansions of , and at x_s0 are

(6)

where

Substituting Eq. (6) into Eq. (5) yields

(7)

Since 1 are linearly independent, the coefficients of like powers of must be equal to zero.

(8)

We can rewrite these equations in a matrix equation form as follows:

(9)

which is the general form of the EPSE n. The boundary conditions are

(10)

where y = 0, and y = ∞. Equation (9), together with the homogeneous boundary conditions (10), constitutes the eigenvalue problem.

2.3 Specific forms of EPSE1 and EPSE2

The form of the EPSE1 is

(11)

which can be rewritten as

(12)

where

Similarly, the form of the EPSE2 is

(12)

where

Equations (12) and (13), together with their homogeneous boundary conditions, constitute the eigenvalue problems.

2.4 Growth rate determination

The complex streamwise wavenumber ã is composed of the growth rate - ã _i and the streamwise wavenumber ã _r. The definition is

(14)

where à is the complex amplitude of the disturbance. In general, the maximum of the streamwise perturbation velocity u' at y_s= y_sm is selected as the complex amplitude, i.e.,

(15)

where α is the eigenvalue.

In the LST, with its parallel flow hypothesis, is independent of x_s. Therefore, the complex streamwise wavenumber ã is equal to the eigenvalue α of the O-S equation. In the EPSE, the complex streamwise wavenumber is composed of two parts, i.e., the eigenvalue and the streamwise variation ratio of the shape function. The relationship between ã and α is

(16)

(17)

3 Swept-wing flow

We consider the compressible flow through an infinite-span swept wing. An NLF(2) -0415 airfoil with the 1.29 m chord length c, the -4 degree attack angle θ, and the 45 degree swept angle Λ is selected. A schematic diagram of the NLF(2) -0415 wing and coordinates is shown in Fig. 1. (X, Y, Z), (x, y, z), (x_s, y_s, z), and (x_t, y_s, z_t) denote the incoming-flow Cartesian coordinate system, the chordwise Cartesian coordinate system, the body-fitted orthogonal curvilinear coordinate system, and the crossflow orthogonal curvilinear coordinate system, respectively.

The freestream Mach number is 0.8, the gas parameters correspond to an altitude of 10 km, the freestream temperature is 223.3 K, and the chord-length Reynolds number Re_c is 8.793×10⁶. In this paper, the quantities with a superscript * are dimensional unless otherwise stated. The velocities, temperature, and density are non-dimensionalized by the freestreams U_∞^*, T_∞^*, and ρ_∞^*, respectively. The length is normalized by c×10^-3. The Reynolds number Re=8 793 is obtained by Re= Re_c×10^-3.

The base flow is computed by the DNS ^[22]. The second-order NND scheme is used for the discretization of the convective terms after flux splitting, the second-order center difference scheme is used for the viscous terms, and the LU-SGS scheme is used for the time terms. The boundary conditions are given as follows: at the outer boundary, the far-field boundary conditions are implemented; at the wall, the no-slip and adiabatic conditions are used; and in the spanwise direction, the periodic conditions are imposed at the boundaries.

3.1 Computational domain and grid

The DNS is computed in the chordwise Cartesian coordinate system (x, y, z). The outer computational domain is approximately 20 times larger than the chord length. The structured grid is used with 400×260×3 grid points in the chordwise, wall-normal, and spanwise directions, respectively. In the wall-normal and chordwise directions, the grid stretches sufficiently to resolve the physically relevant regions, and the size of the minimum wall-normal grid is 0.02 mm. The grid is uniform in the spanwise direction. The computational domain and the grid distribution in the xy-plane are shown in Fig. 2.

3.2 Code verification

To verify the code, we compute a case of Reibert's experiment, in which the chord Reynolds number is Re_c=1.6×10⁶× cos45°. Figure 3(a) shows a comparison of the pressure coefficient obtained by the DNS with that obtained by the experiment. In Figs. 3(b) and 3(c), the distributions of the incoming-flow-direction velocities U _X along the wall-normal direction at x^*/ c=0.1 and x^*/ c=0.4 are compared between the DNS and the experiment. Both the pressure coefficient and the incoming-flow-direction velocities agree with the results of the experiment, suggesting that the code is reliable.

3.3 Results

Let (u₀, v₀, w₀) denote the velocity components in the coordinate system (x_s, y_s, z). Figures 4(a)-4(c) show the distributions of u₀, v₀, and w₀, respectively, along the wall-normal direction at different chordwise positions.

Let w_t denote the crossflow velocity and (ρ w_t')' denote the generalized second-order derivative with y_s of the crossflow. Figure 5 shows the distributions of w_t along the wall-normal direction at different chordwise positions. Figure 6 displays the distributions of (ρ w_t')' along the wall-normal direction at different chordwise positions. There are inflection points, i.e., (ρ w_t')'=0 in Fig. 6, indicating that the crossflow instability can be produced.

4 Non-parallelism effect on instability

To show the non-parallelism effect on the crossflow instability calculated by the DNS, Figs. 7(a) and 7(b) give the amplitudes of the stationary crossflow vortex with A₀=1×10^-4 and β=1.35 and the amplitudes of the unsteady travelling wave with A₀=1×10^-4, β=1.35, and ω=0.1. The results are compared with those predicted by the LST and LPSE. Both figures show that, when the amplitudes predicted by the DNS are less than 0.1, the results of the LPSE agree well with those of the DNS, suggesting that the perturbation is in the linear stage and the LPSE can predict the crossflow instability accurately. However, the amplitudes predicted by the LST are much less than those obtained by the DNS and LPSE, suggesting that the non-parallelism increases the severity of the crossflow instability.

Figures 8(a)-8(b) show the growth rates corresponding to Figs. 7(a) and 7(b), respectively. Both figures show that the difference between the results of the LST and the results of the LPSE is distinct when x^*/ c < 0.3, indicating that the non-parallelism degree is strong in this region whereas small when x^*/ c>0.3, suggesting that the non-parallelism degree is weak in this region.

5 EPSE results 5.1 Neutral curve

Figures 9(a) and 9(b) show the comparisons of the neutral curves predicted by the LST, EPSE1, and EPSE2 at x^*/ c=0.05 and x^*/ c=0.2, respectively. At x^*/ c=0.05, where the non-parallelism degree is strong, the results of the LST are significantly different from those of the EPSE1 and EPSE2, indicating that the non-parallelism effect on the neutral curve can be represented by the EPSE n (n=1, 2). The non-parallelism enlarges the instability range for β>1.5, while narrows the instability range for β < 1.5. At x^*/ c=0.2, where the non-parallelism degree is weaker than that at x^*/ c=0.05, the neutral curves predicted by the LST, EPSE1, and EPSE2 are slightly different except for the range ω < -0.1. Both figures show that the difference between the neutral curves of the EPSE1 and EPSE2 is small, despite the strong non-parallelism at the location x^*/ c=0.05. Overall, the non-parallelism extends the unstable regions.

Figures 10(a)-10(c) show the comparisons of the neutral curves predicted by the EPSE1, EPSE2, and LST for three different frequencies. When ω=0, the upper and lower branches of the neutral curves predicted by the EPSE1 and EPSE2 shift upward relative to those of the LST when x^*/ c < 0.2. The closer the position is to the leading edge, the greater the magnitude of the shift is. When x^*/ c>0.2, the results of the three methods are consistent. For both ω=0.1 and ω=0.2, the lower branches of the three results are consistent, whereas the upper branches predicted by the LST are lower than those of the EPSE1 and EPSE2. The difference between the neutral curves of the EPSE1 and the EPSE2 is small for all these frequencies.

5.2 Growth rates

Figures 11(a)-11(c) and Figs. 12(a)-12(c) display the variations of the growth rates with β for three different frequencies at x^*/ c=0.05 and x^*/ c=0.2, respectively, and the results of the LST, EPSE1, and EPSE2 are compared. At the position x^*/ c=0.05, where the non-parallelism degree is strong, the growth rates predicted by the LST are less than those predicted by the EPSE1 and EPSE2 when β>1.5, and the results of the EPSE1 and EPSE2 are slightly different in the vicinity of the maximum growth rate. At the position x^*/ c=0.2, where the non-parallelism degree is much weaker than that at x^*/ c=0.05, the growth rates predicted by the LST are much closer to those predicted by the EPSE1 and EPSE2 than those shown in Figs. 11(a)-11(c). For ω=0 at x^*/ c=0.05, the maximum growth rates predicted by the LST, EPSE1, and EPSE2 are approximately 0.07, 0.082, and 0.084, respectively, while the spanwise wavenumbers corresponding to the maximum growth rates predicted by them are almost the same, and are approximately equal to 2.5. The non-parallelism degree is found to cause the maximum growth rate to increase, and has a very little effect on the spanwise wavenumber corresponding to the maximum growth rate. The same conclusions are drawn for the other frequencies. The maximum N factor is usually obtained by integrating the maximum growth rate. Therefore, we can expect that the maximum N factor predicted by the LST is much less than those predicted by the EPSE1 and EPSE2 in the e ^N prediction.

5.3 Single-wave linear evolution

Figures 13(a)-13(d) show the evolutions of the amplitudes and the growth rates for four different single waves, where the curves with the sign A denote the amplitudes, and the curves with the sign σ denote the growth rates. The results of the EPSE1 and EPSE2 are compared with those of the LST and LPSE. All figures show that the differences among the growth rates predicted by the LST, EPSE1, EPSE2, and LPSE are small in the positions away from the leading edge, where the non-parallelism degree is weak, whereas near the leading edge, where the non-parallelism degree is strong, the growth rate predicted by the LST for each single wave is significantly different from those of the EPSE1, EPSE2, and LPSE because the LST ignores the non-parallelism effect. The results of the EPSE1 and LPSE for each single wave are found to have some differences, suggesting that the non-parallelism degree is so strong that the first-order approximation of the LPSE is not adequate for the single-wave evolution. Compared with the growth rates and amplitudes of the EPSE1, the results predicted by the EPSE2 are closer to those of the LPSE. This may be because that the EPSE2 is a higher-order approximation to the LPSE than the EPSE1. From Section 4 and this subsection, the four methods can be ranked in the order of accuracy from high to low as LPSE, EPSE2, EPSE1, and LST.

5.4 e ^N calculations

The e ^N method is one of the most popular transition prediction tools in engineering. Because the LPSE is an initial value problem, the envelope of the N factor can be achieved only by integrating the growth rates of all unstable waves when the LPSE is applied to the e ^N prediction. The expression is

(18)

For the LST, the maximum N factor can also be obtained by Eq. (18), but it is generally achieved by integrating the local maximum growth rate, i.e.,

(19)

Because the EPSE is a local eigenvalue problem, as the case for the LST, Eq. (19) can also be used to calculate the maximum N factor for the EPSE.

Compared with Eq. (18), Eq. (19) is simpler and more efficient. For the same calculation tool, such as the LST, Eq. (19) is a more conservative integration strategy, and obtains a maximum factor value higher than that of Eq. (18). Therefore, the engineering designs by use of Eq. (19) are safer and more conservative.

Figure 14(a) shows a comparison of the maximum N factors predicted by the LST, EPSE1, EPSE2, and LPSE for stationary vortices. Figure 14(b) displays a comparison of the maximum N factors predicted by the LST, EPSE1, EPSE2, and LPSE for all unstable disturbances. Equation (18) is used to calculate the envelope of the N factor for the LPSE, whereas Eq. (19) is used to obtain the maximum N factors for the LST and EPSE n (n=1, 2). It can be observed that the LPSE must calculate a large number of N so as to obtain the envelope of the N factor. Therefore, the LPSE is an inflexible and inefficient method for the calculation of the maximum N factor. Although Eq. (19) is more conservative, the maximum N factors predicted by the LST for both the stationary vortices and all unstable disturbances are found to be significantly less than those predicted by the LPSE because of the non-parallelism effect. The maximum N factors for both the stationary vortices and all unstable disturbances predicted by the EPSE1 are almost identical to those of the LPSE in the condition N < 10, whereas they are slightly larger than those of the LPSE when N>10. Although the single-wave growth rates predicted by the EPSE2 are almost the same as those of the LPSE, the maximum N factor calculated by the EPSE2 is slightly larger than that of the LPSE because Eq. (19) is more conservative than Eq. (18).

6 Conclusions

The crossflow instability is the main instability mechanism for swept-wing boundary layers. It occurs at the leading edge of a wing, where the non-parallelism degree is so strong that it has a non-negligible effect on the stability of the boundary layer. Therefore, the LST, which is based on the parallel assumption, cannot predict the crossflow instability in swept-wing boundary-layer flows accurately. However, when the LPSE is used to calculate e ^N, all unstable disturbances with different frequencies and spanwise wavenumbers must be integrated. As a result, the LPSE is an inflexible and inefficient method for the calculation of e ^N. The EPSE can effectively avoid these shortcomings. Therefore, it can be used to study the crossflow instability of swept-wing boundary layers better.

In this paper, we derive the nth-order form of the EPSE firstly, and then use the EPSE1 and EPSE2 to study the swept-wing crossflow instability. This dissertation adopts the EPSE1 and EPSE2 to predict the evolutions of the single-wave amplitudes and the growth rates along the streamwise direction. Compared with the amplitudes and the growth rates predicted by the LST, the results of the EPSE1 and EPSE2 agree well with those of the LPSE. The results of the EPSE2 are almost the same as those of the LPSE. The neutral curves and the variation curves of the growth rates with the spanwise wavenumbers obtained from the EPSE1, EPSE2, and LST are compared. The results show that the stronger the non-parallelism degree is, the greater their differences are, indicating that the LST cannot predict the neutral curves and growth rates accurately. Finally, we study the e ^N calculations, and find that the LPSE is a more time-consuming method for obtaining the envelope of the N factor than the other three methods. The maximum N factors of both the stationary vortices and all unstable disturbances predicted by the EPSE1 are almost identical to those of the LPSE.

Our work shows that, in contrast with the LST and LPSE, the EPSE1 has a better e ^N predictive effect on the calculation accuracy and computational efficiency in the swept-wing boundary layers. The EPSE1 is also expected to be advantageous when it is used to analyze other boundary layers, for which the non-parallelism degree is relatively strong, e.g., complicated 3D boundary layers.