- Email: [email protected]
Transition prediction in infinite swept wings using Navier–Stokes computations and linear stability theory
Computers and Structures 82 (2004) 1551–1560 www.elsevier.com/locate/compstruc
Transition prediction in infinite swept wings using Navier–Stokes computations and linear stability theory J.M.M. Sousa *, L.M.G. Silva Instituto Superior Tecnico/Technical University of Lisbon, Mechanical Engineering Department, Avenue Rovisco Pais, 1049-001 Lisboa, Portugal Received 10 March 2004; accepted 12 March 2004 Available online 1 June 2004
Abstract The coupling of Navier–Stokes computations with eN method for transition prediction in infinite swept wings is demonstrated. In the present paper, results are presented illustrating the successful coupling of these methods for a 45 swept wing operating at four different Reynolds numbers. Moreover, the main issues involved in the generation of adapted meshes and in the computation of the mean flow by a Navier–Stokes method are documented. In all studied cases, the subsequent linear stability analysis of the mean flow has shown that the transition process was dominated by cross-flow instability. The transition prediction results have been validated against experimental transition data from the literature. 2004 Elsevier Ltd. All rights reserved. Keywords: Transition prediction; Infinite swept wing; Linear stability theory; eN method; Navier–Stokes; Adapted mesh
1. Introduction It seems nowadays inevitable that, in a near future, Navier–Stokes computations will supersede the combined use of inviscid flow methods and three-dimensional boundary layer methods in the design of laminar wings. On the other hand, the designer must be able to accurately predict the location of laminar-to-turbulent transition irrespectively of the method employed to compute the flow. This is of paramount importance both for laminar wing design and performance assessment of turbulent wings, so that free transition measurements in wind tunnels may be correctly extrapolated to the larger Reynolds numbers characterizing free flight conditions. Since a few decades, the eN method [1,2] has been applied with considerable success for transition pre-
*
Corresponding author. Tel.: +351-21-841-7320; fax: +35121-849-5241. E-mail address: [email protected] (J.M.M. Sousa).
diction in boundary layer flows computed with viscousinviscid codes. This method was initially developed for two-dimensional flows but later it was extended to more complex problems as well [3]. The classical eN method makes use of linear stability theory of parallel flows (the so-called ‘‘local’’ approach) and correlates the amplification of small disturbances with measured transition locations. As a consequence, it is in fact a semi-empirical criterion because the limiting value of the N-factor, i.e. the value of this parameter at transition, must be determined from experiments. Unfortunately, this value depends on the experimental conditions [3], namely the surface roughness of the model and the flow quality in a specific wind tunnel facility. Essentially driven by efforts towards the development of laminar flow technology for commercial aircraft (see [4] for a review), many measurement campaigns both in wind tunnels and free flight have been undertaken since the 1970s, leading to the calibration of the eN method (see, e.g. [5,6]). However, the large scatter often observed in transition data
0045-7949/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2004.03.051
1552
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
suggests that some old experiments should be revisited in the light of recent developments in transition prediction. Despite the fact that the eN method lacks important features of the physics of transition [7], namely showing inability to account for the ingestion of the instabilities entering the boundary layer and generating the unstable waves (the ‘‘receptivity’’ problem [8]), it is the only practical method presently available for industrial applications [9,10]. Recently, the feasibility of coupling this transition prediction tool with a Navier–Stokes method for airfoil computations has been investigated by Stock and Haase [11], thus unveiling new possibilities for the application of the eN method. In their study, a detailed description of the necessary steps required to produce Navier–Stokes solutions compatible with a linear stability analysis was presented. Subsequently, the successful coupling of these two methods, again for airfoil computations, has been demonstrated as well [10,12]. Benefiting of the knowledge harvested from the aforementioned efforts dedicated to two-dimensional flows, the extension of the methodology to infinite swept wings emerged as the following goal. In this paper, the eN method has been used in conjunction with Navier–Stokes computations for transition prediction in infinite swept wing flow. The successful coupling of these two methods relied on the generation of second-order Navier–Stokes solutions of the mean flow employing adapted meshes. Validation of the implemented methodology was provided from transition experiments carried out in the early 1990s at the Arizona State University Unsteady Wind Tunnel [13]. In the investigated test cases, cross-flow instability was the dominant mechanism in transition. Using the coupled methods, the strong effect of the Reynolds number in this type of instability was properly predicted.
2. Elliptic mesh generation As a consequence of the accuracy requirements for the computation of boundary layers using a Navier– Stokes method and subsequent stability analysis, the mesh generation procedure is of fundamental relevance. This assertion was demonstrated for the case of an airfoil by Stock and Haase [11] and naturally it applies to the case of an infinite swept wing as well. The formulation employed in the Navier–Stokes method (see Section 3) is based on a boundary-fitted curvilinear coordinate system. In the present mesh generation procedure, the coordinate values inside de numerical domain are computed from the curvilinear coordinate values specified at the boundaries. Hence, the problem consists in the determination of the values for the curvilinear coordinates xm ðy i Þ inside the physical
domain, given the respective values defining the complex boundary. A permutation between the dependent variables xm and the independent variables y i is carried out first in order to perform the Navier–Stokes computations in a transformed, uniformly spaced, rectangular plane (the computational domain). As a result, after the coordinate transformation from the physical domain to the computational domain, the numerical mesh generation problem is reduced to the calculation of the values of the Cartesian coordinates y i ðxm Þ inside the latter domain. Again, the inputs to this computation are the values specified at the boundaries. The formulation above constitutes an elliptic problem. The algorithm used for the numerical mesh generation in a plane surface is based on the solution of a system of Poisson equations. As shown, e.g. by Thompson et al. [14], these equations may be written as follows: r2 xm ¼ gmm Pm ;
m ¼ 1; 2;
ð1Þ
m
where x denotes the curvilinear coordinate associated to the direction m, gmm corresponds to the contravariant metric tensor and Pm is a control function. Through the above-mentioned permutation of variables, the set of equations to be solved in the transformed plane may be obtained: XX m
l
gml
X o2 y i oy i mm þ g P ¼ 0; m oxm oxl oxm m m
l
i ¼ 1; 2;
ð2Þ
ox . The specification of boundary values where gml ¼ ox oy j oy j does not present any difficulties because the boundary given by the upper and lower surfaces of the wing and by an external border collocated at a prescribed distance defines always a simply connected domain. In addition, the use of suitable control functions provide the desired control of the number of mesh points (and mesh stretching) in the wall-normal direction where large gradients exist. There exist many possibilities in the definition of the foregoing control functions. However, in the present method, the use of exponential functions (see again [14]) has shown to be particularly efficient. Both the definition of control functions and boundary data have been controlled by the Navier–Stokes solutions in order to obtain the required mesh resolution (adapted mesh) inside the viscous layer, near the leading and trailing edges, and in the wake. The mesh requirements for airfoil calculations using a particular Navier–Stokes method have been reported in previous investigations. The experience obtained from these two-dimensional analyses [10–12] was transposed here to infinite swept wing problems, namely the use of a diagnostic function to evaluate the thickness of the viscous layer in Navier–Stokes computations. In the present case, at least 40 mesh points are placed inside the viscous layer in the adapted
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
meshes following an iterative procedure. Very few iterations are required in the adaptation of the numerical mesh to the Navier–Stokes solution, provided that a good approximation is available from previous computations. The same definition of the diagnostic function given in the aforementioned airfoil studies was employed here. However, in contrast with those investigations, a constant number of mesh points inside the viscous layer was not enforced. The set of equations given by (2) was discretized using second-order accurate central, finite differences except near the boundaries where sometimes forward- or backward-differences had to be used. The discretization procedure produces a set of algebraic equations, which was solved by the application of the successive overrelaxation algorithm. Typically, a value of 1.8 has been used for the relaxation parameter. On the other hand, the convergence criterion had to be adjustable so that the quality of the mesh in the vicinity of the wall on the region where the viscous layer was thinner could be controlled, as a consequence of the massive clustering of mesh points demanded in this area. The numerical meshes employed in the Navier– Stokes computations of the flow around an infinite swept wing described in the remainder of this paper have been generated using the above-described methodology. It must be mentioned that elliptic mesh generation strongly facilitates the implementation of local mesh refinements without introducing unwanted discontinuities. This is a key requirement for the generation of boundary layer profiles to be subsequently subjected to a stability analysis. Fig. 1 shows a typical mesh generated for the Navier– Stokes computations constructed with 512 · 128 mesh nodes. As the angle of attack was not changed in the computations and the Reynolds number range was not very wide this mesh suffered only minor adjustments from case to case.
1553
3. Finite volume Navier–Stokes method Although all test cases considered here correspond to incompressible flows, the Navier–Stokes calculation procedure used throughout the present investigation consisted on a variant of the pressure-correction method for the computation of flows at all speeds by Kobayashi and Pereira [15]. In fact, all methods employed herein have been developed for compressible flows and, therefore, major difficulties for the extension of the present methodology to such cases are not anticipated. The aforementioned Navier–Stokes method is briefly described below. The steady continuum flow of a non-reacting, thermally conductive, Newtonian fluid was considered here. The mathematical models expressing transport of mass, momentum and energy are summarized in the Navier– Stokes equations, which in a curvilinear, oblique coordinate system and for the Cartesian components of vectors and tensors are as follows: 1 g1=2 1 g1=2 1 g1=2
o i m qv Ci ¼ 0; oxm
ð3Þ
i o h i j ðqv v þ sij ÞCjm ¼ 0; m ox
ð4Þ
o i 1 o ðqv h þ qi ÞCim ¼ 1=2 m vi dik skj Cjm ; m ox g ox
ð5Þ
respectively for continuity, motion and energy. In Eqs. m (3)–(5), Cim ¼ g1=2 ox , g is the determinant of the metric oy i i tensor, v represents the velocity field and h is the enthalpy. The summation convention is used in the expressions above, whereby summation is implied when an index is repeated for diagonal positions. In addition, for a Newtonian fluid, the stress tensor sij can be written as:
Fig. 1. Numerical mesh generated for the NLF(2)-0415 geometry (partial view).
1554
sij ¼
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
i 2 l ov m kj ovj m ki þ ; p þ lH dij 1=2 C d C d k k oxm oxm 3 g ð6Þ
where p is the static pressure, l is the dynamic viscosity of the fluid and H denotes the velocity divergence. In turn, a Fourier law defines the heat flux vector: qi ¼
k oT m ji C d ; g1=2 oxm j
ð7Þ
where k is the thermal conductivity and T is the static temperature. Equations of continuity, motion and energy can individually be re-written in the form of a general transport equation, as follows: o i ðqv / þ ci ÞDmi ¼ g1=2 s/ ; oxm
ð8Þ
where / is a general transported scalar with a diffusivity coefficient j, Dmi denotes the transpose of the cofactor of oy i in the Jacobian of the coordinate transformation oxm y i ðxm Þ given in the previous section, s/ is a linearized source term and the diffusive flux of / is defined by: m Dj o/ ji ci ¼ j 1=2 d : ð9Þ g oxm The closure of the governing equations is achieved by further introducing the equation of state. In here, the perfect gas relation is considered: p ¼ qRT ;
R ¼ R0 =M;
fore and aft in Eqs. (3)–(9) run from one to three. However, all derivatives along the third (spanwise) dimension of space are neglected. As a result, the equation for the computation of spanwise motion, obtained from Eq. (8) for i ¼ 3 with the aforementioned assumption, consists in the simple convection–diffusion of the third scalar component by the two remaining components of the velocity field. The discrete form of Eq. (8) is obtained by using the finite volume method based on a collocated, oblique mesh. Detailed descriptions of the corresponding discretization procedures may be found profusely in the literature [16–18]. Briefly, it can be said that the integration of Eq. (8) over the numerical domain and subsequent application of the divergence theorem produces a flux balance at each control volume. The resulting total flux of the general transported scalar (component) / is formed by a convective and a diffusive part. Whereas the diffusion term was always approximated here by central differences, the evaluation of the convective flux of / was initially performed using first-order upwinding for robustness. At later stages of convergence, higher-order accuracy was obtained in the present
Table 1 Selected test cases and experimentally determined transition locations [13]
ð10Þ
where R0 stands for the universal gas constant and M is the molecular weight of the gas. Under the infinite swept wing approach, all vector and tensor component indices
Reynolds number: Re 106
Experimental transition location: x=c
1.92 2.37 2.73 3.73
0.78 0.58 0.45 0.30
1.0 0.8
Model experiments [13] 0.6
Upper surface - Lower end Upper surface - Upper end
-Cp
0.4 0.2 0.0
Numerical Predictions Re = 1,920,000 Re = 2,370,000 Re = 2,730,000 Re = 3,730,000
-0.2 -0.4 -0.6 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
x /c
Fig. 2. Final pressure distributions obtained from Navier–Stokes computations.
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
method by implementing a second-order reconstruction of the variable within each control volume making use of the ‘‘minmod’’ limiter. Again, this is a key issue because the stability analysis (see Section 4) is very sensitive to the effect of numerical viscosity on the mean flow computed by the Navier–Stokes method. Finally, it must be mentioned that the solution procedure applied to the nonlinear system of equations resulting from the discretization process is based on a semi-implicit pressurecorrection technique. See again Kobayashi and Pereira [15] for full details. The present method has been used before for the computation of flows governed by the Reynolds-averaged Navier–Stokes equations [19]. The extension of the standard k e turbulence model to low-Reynolds number flows by Coakley and Huang [20] was used in this case in order to describe the near wall region downstream of the transition location. Wall damping effects were incorporated in the model by means of turbulent Reynolds number dependent coefficients. The procedure to determine the transition location is discussed in the next section.
10
The eN method is used in the present work for transition prediction based on the linear, local stability analysis of the boundary layer (see, e.g. [21]). The basic procedure for this analysis is described by Malik [22] and therefore only a brief explanation of the method is given here. According to the conventional normal mode assumption used to derive the linear stability equations in a coordinate system ðx; y; zÞ linked to the wing, the eigensolutions take the form: ^ exp ½ið ax þ bz xtÞ; Eðx; y; z; tÞ ¼ EðyÞ ð11Þ pffiffiffiffiffiffiffi ^ where i ¼ 1 and EðyÞ is a disturbance profile. With the mean flow available from previous computations (see Section 3), the stability problem consists in the determination of six unknowns, namely ar , ai , br , bi , xr and xi . These are the streamwise (x-direction) wavenumber and growth rate, spanwise (z-direction) wavenumber and growth rate, and wave frequency and temporal growth rate, respectively. In the present
0 Hz - 0.005 c 0 Hz - 0.006 c 0 Hz - 0.007 c 0 Hz - 0.008 c 0 Hz - 0.009 c 50 Hz - 0.007 c 50 Hz - 0.008 c 100 Hz - 0.005 c 100 Hz - 0.006 c 100 Hz - 0.007 c 100 Hz - 0.008 c 100 Hz - 0.009 c 150 Hz - 0.007 c 150 Hz - 0.008 c 150 Hz - 0.009 c
8.5 8
7 6.4 6
N-factor
4. Transition prediction method
Re = 1,920,000
9
5
4
3
2
1
0 0.0
0.1
0.2
0.3
1555
0.4
0.5
0.6
0.7
0.8
x/c Fig. 3. Final N -factor curves computed for Re ¼ 1; 920; 000.
1556
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
to Eq. (12), which constitutes a block-tridiagonal system of equations. For the present case, matrix Bij is invertible and a generalized matrix eigenvalue problem, solvable by LR methods, can be formulated. An inverse Rayleigh iteration procedure, together with a block LU factorization to invert the block-tridiagonal system, is further used to perform the local search (or improvement) of the eigenvalues. Linear stability theory provides information regarding the evolution of unstable disturbances in the boundary layer but, alone it is unable to predict the occurrence of transition. This task is accomplished by the use of the semi-empirical eN method where the stability characteristics of the flow are analyzed in terms of local, disturbance growth rates, which in turn are integrated to yield the N -factor. Thus, upon solving the eigenvalue problem and introducing Gaster’s transformation, the amplification rate for a specific disturbance in the form given by Eq. (11) can be obtained from the following expression:
method, the temporal formulation of the problem is preferred to the spatial analysis because the latter is more expensive from the computational point of view. Hence, a and b are real numbers and x is a complex number that is found by solving the stability equations. Temporal growth rates xi are afterwards converted into spatial growth rates ai by the straightforward application of a generalized Gaster’s transformation [23]. The (compressible) stability eigenvalue problem is solved here using a matrix finite-difference method to obtain local solutions when a good guess for the eigenvalue is available. When this is not the case, a globally convergent eigenvalue procedure is used instead. The formulation of the stability equations in matrix form leads to the expression: Aij Ej ¼ xBij Ej ;
ð12Þ
where x naturally stands for the eigenvalue and Ej is the discrete representation of the eigenfunction. The eigenvalue is calculated by the determinant condition applied
10
Re = 2,370,000 0 Hz - 0.004 c 0 Hz - 0.005 c 0 Hz - 0.006 c 0 Hz - 0.007 c 0 Hz - 0.008 c 100 Hz - 0.005 c 100 Hz - 0.006 c 100 Hz - 0.007 c 100 Hz - 0.008 c 150 Hz - 0.006 c 150 Hz - 0.007 c 200 Hz - 0.005 c 200 Hz - 0.006 c 200 Hz - 0.007 c 200 Hz - 0.008 c
9 8.5 8
7
N-factor
6.4 6
5
4
3
2
1
0 0.0
0.1
0.2
0.3
0.4
0.5
x/c Fig. 4. Final N -factor curves computed for Re ¼ 2; 370; 000.
0.6
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
A ¼ exp A0
Z
xi dxg ; Vg
parameter rðxr ; bÞ. Consequently, in order to compute the N-factor, some integration strategy must be followed. However, it must be noted that different integration strategies usually lead to distinct values of the N-factor as well, as shown e.g. by Pereira and Sousa [24]. Thus, the integration strategy should always display an unequivocal link with the physics of the dominant transition mechanism in the considered flow cases. In the present study, the chosen strategy can be described as following ‘‘waves’’ (or modes) of fixed fre2p quency xr and constant total wavelength k ¼ pffiffiffiffiffiffiffiffiffi . 2 2
ð13Þ
^ where A stands for the amplitude of EðyÞ at any particular value of the wall-normal coordinate y, the subscript 0 denotes the value where amplification started, Vg represents the norm of the group velocity vector with r oxr components ðox ; ob Þ and dxg is an infinitesimal arcoa length along the group velocity direction. Hence, for a wave of fixed frequency, Eq. (13) may be rewritten as: Z st A ¼ exp r ds; ð14Þ A0 s0
a þb
This method was first employed by Srokowski and Orzsag [25] for stationary modes only (physically associated with stationary cross-flow vortices). Using this strategy as well, travelling modes have been accounted for by Hefner and Bushnell [7]. Here, the N -factor was computed by maximizing the total amplification of modes characterized by a fixed frequency xr within the unstable range with respect to k, as given by the following expression:
where r is the characteristic growth rate of the disturbance, which is integrated from the neutral point s0 to the transition point st . Looking at Eq. (14), it seems that the computation of the N -factor is now similar to that in a two-dimensional case. However, the additional degree of freedom appearing in a three-dimensional problem is hidden in the
10
Re = 2,730,000 0 Hz - 0.003 c 0 Hz - 0.004 c 0 Hz - 0.005 c 0 Hz - 0.006 c 0 Hz - 0.007 c 100 Hz - 0.004 c 100 Hz - 0.005 c 100 Hz - 0.006 c 100 Hz - 0.007 c 150 Hz - 0.005 c 150 Hz - 0.006 c 200 Hz - 0.004 c 200 Hz - 0.005 c 200 Hz - 0.006 c 200 Hz - 0.007 c
9 8.5 8
7
N - factor
6.4 6
5
4
3
2
1
0 0.0
0.1
0.2
1557
0.3
0.4
0.5
x/c Fig. 5. Final N -factor curves computed for Re ¼ 2; 730; 000.
0.6
1558
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
10
Re = 3,730,000 0 Hz - 0.002 c 0 Hz - 0.003 c 0 Hz - 0.004 c 0 Hz - 0.005 c 0 Hz - 0.006 c 100 Hz - 0.003 c 100 Hz - 0.004 c 100 Hz - 0.005 c 200 Hz - 0.003 c 200 Hz - 0.004 c 200 Hz - 0.005 c 300 Hz - 0.003 c 300 Hz - 0.004 c 300 Hz - 0.005 c 400 Hz - 0.004 c
9 8.5 8
7
N-factor
6.4 6
5
4
3
2
1
0 0.0
0.1
0.2
0.3
0.4
x/c Fig. 6. Final N -factor curves computed for Re ¼ 3; 730; 000.
Table 2 Computed transition locations obtained for stationary and travelling disturbances
Both stationary and travelling modes are to be considered in this analysis as described in the next section.
Reynolds number: Re 106
Computed transition location: x=c Stationary: N ¼ 6:4
Travelling: N ¼ 8:5
5. Results and discussion
1.92 2.37 2.73 3.73
0.75 0.58 0.49 0.31
0.74 0.58 0.49 0.32
In this section, results from coupling Navier–Stokes computations of the mean flow around an infinite swept wing and eN transition predictions are presented and discussed. The selected test cases are constructed with the airfoil NLF(2)-0415 [26], which was designed as a low-drag wing to be used on a commuter aircraft with unswept wings. In the work of Dagenhart [13], detailed results from wind tunnel experiments of a 45 swept wing were obtained by operating the NLF(2)-0415 wing section at a small negative angle of attack a ¼ 4. At this condition, the experimental model was described as a nearly ideal cross-flow generator because the favorable pressure gradient on the upper surface keeps the Tollmien–Schlicthting instability subcritical back to 71% of chord for a 6 0. Attachment-line and G€ ortler insta-
" N ¼ max ln k
A A0
A A0
¼ exp ðx;kÞ
#
; ðx;kÞ
Z
st
r xr ; k ds;
ð15Þ
s0
where an asterisk denotes fixed values. Eq. (15) indicates that the computed N -factor to be used in the prediction of transition actually results from an envelope of curves.
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
1559
1.00 0.90
Experiments [13] Coupled method
0.80 0.70
x/c
0.60 0.50 0.40 0.30 0.20 0.10 0.00 1.5
Transition locations N = 6.4 for stationary disturbances
2.0
2.5
3.0
3.5
4.0
Re x 10-6
Fig. 7. Comparison between numerical predictions of transition and experiments.
bility mechanisms were not found for the investigated test conditions. A strong dependence of the transition location with the Reynolds number Re was observed in the experiments [13], due to the effect of this parameter in the strength of stationary cross-flow vortices. It was expected that the present numerical procedure would be able to reproduce the experimentally observed behavior. The Reynolds numbers and transition locations on the upper surface of the swept wing obtained from the published experiments (Naphtalene flow visualizations) for the investigated test cases are shown in Table 1. Fig. 2 shows the final pressure distributions (in terms of the non-dimensional pressure coefficient Cp ) obtained from the Navier–Stokes computations coupled with the eN method for transition prediction and covering all test cases indicated in Table 1. The results correspond to computations carried out on adapted numerical meshes. Fast convergence has been always achieved in the coupled procedure applied to the upper surface of the wing only. On the lower surface transition was fixed at the leading edge for simplicity. Experimental values obtained from the work of Dagenhart [13] are also displayed in Fig. 2 for comparison purposes (upper surface only). Excellent agreement between numerical calculations and experimental values for the (upper surface) lower end of the model can be observed in the graph. On the other hand, the pressure gradient on the upper end is slightly underpredicted by the computations. In the determination of the transition location resulting from the coupled procedure, two different approaches have been applied to the obtained envelopes of N -factor curves. The first considered stationary modes only in the linear, local stability analysis. In this case, the best correlation with the experimental transition loca-
tions was obtained and a value of the N -factor at transition of 6.4 applies. The second considered travelling modes in the analysis and therefore the maximization in Eq. (15) was extended to the frequency as well. In this case, a satisfactory correlation with the experimental transition locations was obtained for N-factor values of 8.5. These values of limiting N-factor are acceptable for experiments carried out in a wind tunnel environment. Figs. 3–6 present the final N-factor curves for each test case. The legends in the figures list the range of frequencies and wavelengths (expressed as a fraction of the chord c ¼ 1:83 m) of the disturbances employed in the analysis. Filled symbols have been used to denote stationary modes, whereas open symbols were used for travelling modes. Flow instability is detected very early in all cases ðx=c < 0:05Þ, which is in consonance with the development of cross-flow instability. The corresponding values of wave propagation direction (not shown) corroborated this assertion. A summary of the transition locations obtained by assuming the aforementioned values for the N-factor at transition in all test cases is also shown in Table 2. In addition, the excellent agreement between the results of the coupled procedure for stationary modes and the experimental data of Dagenhart [13] is confirmed by the transition plot portrayed in Fig. 7. Experimental error bars obtained from the same reference are also shown in this figure.
6. Conclusions In the present work, a demonstration of successful coupling between Navier–Stokes computations and eN method for transition prediction using linear, local stability analysis applied to an infinite swept wing has been
1560
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
presented. The procedures for the elliptic generation of adapted numerical meshes and to obtain second-order accurate Navier–Stokes solutions of the mean flow have been described as well. These constitute fundamental requirements for the success of the methodology. Four test cases have been constructed using the NLF(2)-0415 wing section at a sweep angle of 45 and varying the Reynolds number from 1.92 to 3.73 million. In all cases, the transition process was dominated by cross-flow instability. Additionally, the implemented methodology was able to reproduce the strong dependence of the transition location with the Reynolds number. The results have shown that, for the investigated test cases, excellent agreement with published experimental transition data can be obtained from the stability analysis of stationary modes only, by assuming a value of the N-factor at transition of 6.4. If travelling modes are also considered in the analysis, satisfactory agreement with experimental results may also be obtained by taking N ¼ 8:5 at transition instead.
Acknowledgements This research work has been partially supported by the European Commission through project ALTTA within the framework of the GROWTH programme.
References [1] Smith AMO, Gamberoni N. Transition, pressure gradient and stability theory. Report No. ES26388, Douglas Aircraft Co., 1956. [2] Van Ingen JL. A suggested semi-empirical method for the calculation of the boundary layer transition region. Report VTH-74, University of Delft, The Netherlands, 1956. [3] Arnal D. Boundary layer transition: predictions based on linear theory. AGARD Report No. 793. Special course on progress in transition modelling, 1994. [4] Joslin RD. Aircraft laminar flow control. Ann Rev Fluid Mech 1998;30:1–29. [5] Bushnell DM, Malik MR, Harvey WD. Transition prediction in external flows via linear stability theory. In: Zierep J, Oertel H, editors. IUTAM Symposium Transsonicum III. Berlin: Springer-Verlag; 1988. [6] Schrauf G. Large-scale laminar flow tests evaluated with linear stability theory. AIAA Paper No. 2001-2444, 2001. [7] Hefner JN, Bushnell DM. Status of linear boundary-layer stability theory and the eN method, with emphasis on swept-wing applications. NASA Technical Paper No. 1645, 1980.
[8] Morkovin MV. On the many faces of transition. In: Wells CS, editor. Viscous drag reduction. New York: Plenum Press; 1969. [9] Arnal D, Casalis G. Laminar-turbulent transition prediction in three-dimensional flows. Prog Aerospace Sci 2000; 36(2):173–91. [10] Stock HW. Airfoil validation using coupled Navier–Stokes and eN transition prediction methods. J Aircraft 2002; 39(1):51–8. [11] Stock HW, Haase W. A feasibility study of eN transition prediction in Navier–Stokes methods for two-dimensional airfoil computations. AIAA J 1999;37(10):1187–96. [12] Stock HW, Haase W. Navier–Stokes airfoil computations with eN transition prediction including transitional flow regions. AIAA J 2000;38(11):2059–66. [13] Dagenhart JR. Crossflow stability and transition experiments in a swept-wing flow. NASA Technical Memorandum No. 108650, 1992. [14] Thompson JF, Thames FC, Mastin CW. TOMCAT––A code for numerical generation of boundary-fitted curvilinear coordinate systems on fields containing any number of arbitrary two-dimensional bodies. J Comput Phys 1977; 24(3):274–302. [15] Kobayashi MH, Pereira JCF. Characteristic-based pressure correction at all speeds. AIAA J 1996;34(2):272–80. [16] Anderson DA, Tannehill JC, Pletcher RH. Computational fluid mechanics and heat transfer. New York: Hemisphere; 1984. [17] Fletcher CAJ. In: Computational techniques for fluid dynamics, vol. 1/2. Berlin: Springer-Verlag; 1988. [18] Hirsch C. In: Numerical computation of internal and external flows, vol. 1/2. NewYork: John Wiley; 1990. [19] Pereira JCF, Kobayashi MH, Marques NPC. Numerical solution of the turbulent flow over a fence using twoequation models. In: Dervieux A, Braza M, Dussauge JP, editors. Notes on numerical fluid mechanics, vol. 65. Braunschweig: Vieweg-Verlag; 1998. [20] Coakley TJ, Huang PG. Turbulence modelling for highspeed flows. AIAA Paper No. 92-0436, 1992. [21] Mack LM. Boundary-layer linear stability theory. AGARD Report No. 709. Special course on stability and transition of laminar flows, 1984. [22] Malik MR. Finite-difference solution of the compressible stability eigenvalue problem. NASA Contractor Report No. 3584, 1982. [23] Gaster M. A note on the relation between temporally increasing and spatially increasing disturbances in hydrodynamic stability. J Fluid Mech 1962;14:222–4. [24] Pereira JCF, Sousa JMM. Comparison of several N -factor integration strategies in a compressible boundary layer over a swept curved surface. ICAS Paper No. 94-2.4.2, 1994. [25] Srokowski A, Orszag SA. Mass flow requirements for LFC wing design. AIAA Paper No. 77-1222, 1977. [26] Somers DM, Horstmann KH. Design of a medium-speed natural-laminar-flow airfoil for commuter aircraft applications. DFVLR Report IB/29-85/26, 1985.
Transition prediction in infinite swept wings using Navier–Stokes computations and linear stability theory J.M.M. Sousa *, L.M.G. Silva Instituto Superior Tecnico/Technical University of Lisbon, Mechanical Engineering Department, Avenue Rovisco Pais, 1049-001 Lisboa, Portugal Received 10 March 2004; accepted 12 March 2004 Available online 1 June 2004
Abstract The coupling of Navier–Stokes computations with eN method for transition prediction in infinite swept wings is demonstrated. In the present paper, results are presented illustrating the successful coupling of these methods for a 45 swept wing operating at four different Reynolds numbers. Moreover, the main issues involved in the generation of adapted meshes and in the computation of the mean flow by a Navier–Stokes method are documented. In all studied cases, the subsequent linear stability analysis of the mean flow has shown that the transition process was dominated by cross-flow instability. The transition prediction results have been validated against experimental transition data from the literature. 2004 Elsevier Ltd. All rights reserved. Keywords: Transition prediction; Infinite swept wing; Linear stability theory; eN method; Navier–Stokes; Adapted mesh
1. Introduction It seems nowadays inevitable that, in a near future, Navier–Stokes computations will supersede the combined use of inviscid flow methods and three-dimensional boundary layer methods in the design of laminar wings. On the other hand, the designer must be able to accurately predict the location of laminar-to-turbulent transition irrespectively of the method employed to compute the flow. This is of paramount importance both for laminar wing design and performance assessment of turbulent wings, so that free transition measurements in wind tunnels may be correctly extrapolated to the larger Reynolds numbers characterizing free flight conditions. Since a few decades, the eN method [1,2] has been applied with considerable success for transition pre-
*
Corresponding author. Tel.: +351-21-841-7320; fax: +35121-849-5241. E-mail address: [email protected] (J.M.M. Sousa).
diction in boundary layer flows computed with viscousinviscid codes. This method was initially developed for two-dimensional flows but later it was extended to more complex problems as well [3]. The classical eN method makes use of linear stability theory of parallel flows (the so-called ‘‘local’’ approach) and correlates the amplification of small disturbances with measured transition locations. As a consequence, it is in fact a semi-empirical criterion because the limiting value of the N-factor, i.e. the value of this parameter at transition, must be determined from experiments. Unfortunately, this value depends on the experimental conditions [3], namely the surface roughness of the model and the flow quality in a specific wind tunnel facility. Essentially driven by efforts towards the development of laminar flow technology for commercial aircraft (see [4] for a review), many measurement campaigns both in wind tunnels and free flight have been undertaken since the 1970s, leading to the calibration of the eN method (see, e.g. [5,6]). However, the large scatter often observed in transition data
0045-7949/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2004.03.051
1552
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
suggests that some old experiments should be revisited in the light of recent developments in transition prediction. Despite the fact that the eN method lacks important features of the physics of transition [7], namely showing inability to account for the ingestion of the instabilities entering the boundary layer and generating the unstable waves (the ‘‘receptivity’’ problem [8]), it is the only practical method presently available for industrial applications [9,10]. Recently, the feasibility of coupling this transition prediction tool with a Navier–Stokes method for airfoil computations has been investigated by Stock and Haase [11], thus unveiling new possibilities for the application of the eN method. In their study, a detailed description of the necessary steps required to produce Navier–Stokes solutions compatible with a linear stability analysis was presented. Subsequently, the successful coupling of these two methods, again for airfoil computations, has been demonstrated as well [10,12]. Benefiting of the knowledge harvested from the aforementioned efforts dedicated to two-dimensional flows, the extension of the methodology to infinite swept wings emerged as the following goal. In this paper, the eN method has been used in conjunction with Navier–Stokes computations for transition prediction in infinite swept wing flow. The successful coupling of these two methods relied on the generation of second-order Navier–Stokes solutions of the mean flow employing adapted meshes. Validation of the implemented methodology was provided from transition experiments carried out in the early 1990s at the Arizona State University Unsteady Wind Tunnel [13]. In the investigated test cases, cross-flow instability was the dominant mechanism in transition. Using the coupled methods, the strong effect of the Reynolds number in this type of instability was properly predicted.
2. Elliptic mesh generation As a consequence of the accuracy requirements for the computation of boundary layers using a Navier– Stokes method and subsequent stability analysis, the mesh generation procedure is of fundamental relevance. This assertion was demonstrated for the case of an airfoil by Stock and Haase [11] and naturally it applies to the case of an infinite swept wing as well. The formulation employed in the Navier–Stokes method (see Section 3) is based on a boundary-fitted curvilinear coordinate system. In the present mesh generation procedure, the coordinate values inside de numerical domain are computed from the curvilinear coordinate values specified at the boundaries. Hence, the problem consists in the determination of the values for the curvilinear coordinates xm ðy i Þ inside the physical
domain, given the respective values defining the complex boundary. A permutation between the dependent variables xm and the independent variables y i is carried out first in order to perform the Navier–Stokes computations in a transformed, uniformly spaced, rectangular plane (the computational domain). As a result, after the coordinate transformation from the physical domain to the computational domain, the numerical mesh generation problem is reduced to the calculation of the values of the Cartesian coordinates y i ðxm Þ inside the latter domain. Again, the inputs to this computation are the values specified at the boundaries. The formulation above constitutes an elliptic problem. The algorithm used for the numerical mesh generation in a plane surface is based on the solution of a system of Poisson equations. As shown, e.g. by Thompson et al. [14], these equations may be written as follows: r2 xm ¼ gmm Pm ;
m ¼ 1; 2;
ð1Þ
m
where x denotes the curvilinear coordinate associated to the direction m, gmm corresponds to the contravariant metric tensor and Pm is a control function. Through the above-mentioned permutation of variables, the set of equations to be solved in the transformed plane may be obtained: XX m
l
gml
X o2 y i oy i mm þ g P ¼ 0; m oxm oxl oxm m m
l
i ¼ 1; 2;
ð2Þ
ox . The specification of boundary values where gml ¼ ox oy j oy j does not present any difficulties because the boundary given by the upper and lower surfaces of the wing and by an external border collocated at a prescribed distance defines always a simply connected domain. In addition, the use of suitable control functions provide the desired control of the number of mesh points (and mesh stretching) in the wall-normal direction where large gradients exist. There exist many possibilities in the definition of the foregoing control functions. However, in the present method, the use of exponential functions (see again [14]) has shown to be particularly efficient. Both the definition of control functions and boundary data have been controlled by the Navier–Stokes solutions in order to obtain the required mesh resolution (adapted mesh) inside the viscous layer, near the leading and trailing edges, and in the wake. The mesh requirements for airfoil calculations using a particular Navier–Stokes method have been reported in previous investigations. The experience obtained from these two-dimensional analyses [10–12] was transposed here to infinite swept wing problems, namely the use of a diagnostic function to evaluate the thickness of the viscous layer in Navier–Stokes computations. In the present case, at least 40 mesh points are placed inside the viscous layer in the adapted
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
meshes following an iterative procedure. Very few iterations are required in the adaptation of the numerical mesh to the Navier–Stokes solution, provided that a good approximation is available from previous computations. The same definition of the diagnostic function given in the aforementioned airfoil studies was employed here. However, in contrast with those investigations, a constant number of mesh points inside the viscous layer was not enforced. The set of equations given by (2) was discretized using second-order accurate central, finite differences except near the boundaries where sometimes forward- or backward-differences had to be used. The discretization procedure produces a set of algebraic equations, which was solved by the application of the successive overrelaxation algorithm. Typically, a value of 1.8 has been used for the relaxation parameter. On the other hand, the convergence criterion had to be adjustable so that the quality of the mesh in the vicinity of the wall on the region where the viscous layer was thinner could be controlled, as a consequence of the massive clustering of mesh points demanded in this area. The numerical meshes employed in the Navier– Stokes computations of the flow around an infinite swept wing described in the remainder of this paper have been generated using the above-described methodology. It must be mentioned that elliptic mesh generation strongly facilitates the implementation of local mesh refinements without introducing unwanted discontinuities. This is a key requirement for the generation of boundary layer profiles to be subsequently subjected to a stability analysis. Fig. 1 shows a typical mesh generated for the Navier– Stokes computations constructed with 512 · 128 mesh nodes. As the angle of attack was not changed in the computations and the Reynolds number range was not very wide this mesh suffered only minor adjustments from case to case.
1553
3. Finite volume Navier–Stokes method Although all test cases considered here correspond to incompressible flows, the Navier–Stokes calculation procedure used throughout the present investigation consisted on a variant of the pressure-correction method for the computation of flows at all speeds by Kobayashi and Pereira [15]. In fact, all methods employed herein have been developed for compressible flows and, therefore, major difficulties for the extension of the present methodology to such cases are not anticipated. The aforementioned Navier–Stokes method is briefly described below. The steady continuum flow of a non-reacting, thermally conductive, Newtonian fluid was considered here. The mathematical models expressing transport of mass, momentum and energy are summarized in the Navier– Stokes equations, which in a curvilinear, oblique coordinate system and for the Cartesian components of vectors and tensors are as follows: 1 g1=2 1 g1=2 1 g1=2
o i m qv Ci ¼ 0; oxm
ð3Þ
i o h i j ðqv v þ sij ÞCjm ¼ 0; m ox
ð4Þ
o i 1 o ðqv h þ qi ÞCim ¼ 1=2 m vi dik skj Cjm ; m ox g ox
ð5Þ
respectively for continuity, motion and energy. In Eqs. m (3)–(5), Cim ¼ g1=2 ox , g is the determinant of the metric oy i i tensor, v represents the velocity field and h is the enthalpy. The summation convention is used in the expressions above, whereby summation is implied when an index is repeated for diagonal positions. In addition, for a Newtonian fluid, the stress tensor sij can be written as:
Fig. 1. Numerical mesh generated for the NLF(2)-0415 geometry (partial view).
1554
sij ¼
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
i 2 l ov m kj ovj m ki þ ; p þ lH dij 1=2 C d C d k k oxm oxm 3 g ð6Þ
where p is the static pressure, l is the dynamic viscosity of the fluid and H denotes the velocity divergence. In turn, a Fourier law defines the heat flux vector: qi ¼
k oT m ji C d ; g1=2 oxm j
ð7Þ
where k is the thermal conductivity and T is the static temperature. Equations of continuity, motion and energy can individually be re-written in the form of a general transport equation, as follows: o i ðqv / þ ci ÞDmi ¼ g1=2 s/ ; oxm
ð8Þ
where / is a general transported scalar with a diffusivity coefficient j, Dmi denotes the transpose of the cofactor of oy i in the Jacobian of the coordinate transformation oxm y i ðxm Þ given in the previous section, s/ is a linearized source term and the diffusive flux of / is defined by: m Dj o/ ji ci ¼ j 1=2 d : ð9Þ g oxm The closure of the governing equations is achieved by further introducing the equation of state. In here, the perfect gas relation is considered: p ¼ qRT ;
R ¼ R0 =M;
fore and aft in Eqs. (3)–(9) run from one to three. However, all derivatives along the third (spanwise) dimension of space are neglected. As a result, the equation for the computation of spanwise motion, obtained from Eq. (8) for i ¼ 3 with the aforementioned assumption, consists in the simple convection–diffusion of the third scalar component by the two remaining components of the velocity field. The discrete form of Eq. (8) is obtained by using the finite volume method based on a collocated, oblique mesh. Detailed descriptions of the corresponding discretization procedures may be found profusely in the literature [16–18]. Briefly, it can be said that the integration of Eq. (8) over the numerical domain and subsequent application of the divergence theorem produces a flux balance at each control volume. The resulting total flux of the general transported scalar (component) / is formed by a convective and a diffusive part. Whereas the diffusion term was always approximated here by central differences, the evaluation of the convective flux of / was initially performed using first-order upwinding for robustness. At later stages of convergence, higher-order accuracy was obtained in the present
Table 1 Selected test cases and experimentally determined transition locations [13]
ð10Þ
where R0 stands for the universal gas constant and M is the molecular weight of the gas. Under the infinite swept wing approach, all vector and tensor component indices
Reynolds number: Re 106
Experimental transition location: x=c
1.92 2.37 2.73 3.73
0.78 0.58 0.45 0.30
1.0 0.8
Model experiments [13] 0.6
Upper surface - Lower end Upper surface - Upper end
-Cp
0.4 0.2 0.0
Numerical Predictions Re = 1,920,000 Re = 2,370,000 Re = 2,730,000 Re = 3,730,000
-0.2 -0.4 -0.6 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
x /c
Fig. 2. Final pressure distributions obtained from Navier–Stokes computations.
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
method by implementing a second-order reconstruction of the variable within each control volume making use of the ‘‘minmod’’ limiter. Again, this is a key issue because the stability analysis (see Section 4) is very sensitive to the effect of numerical viscosity on the mean flow computed by the Navier–Stokes method. Finally, it must be mentioned that the solution procedure applied to the nonlinear system of equations resulting from the discretization process is based on a semi-implicit pressurecorrection technique. See again Kobayashi and Pereira [15] for full details. The present method has been used before for the computation of flows governed by the Reynolds-averaged Navier–Stokes equations [19]. The extension of the standard k e turbulence model to low-Reynolds number flows by Coakley and Huang [20] was used in this case in order to describe the near wall region downstream of the transition location. Wall damping effects were incorporated in the model by means of turbulent Reynolds number dependent coefficients. The procedure to determine the transition location is discussed in the next section.
10
The eN method is used in the present work for transition prediction based on the linear, local stability analysis of the boundary layer (see, e.g. [21]). The basic procedure for this analysis is described by Malik [22] and therefore only a brief explanation of the method is given here. According to the conventional normal mode assumption used to derive the linear stability equations in a coordinate system ðx; y; zÞ linked to the wing, the eigensolutions take the form: ^ exp ½ið ax þ bz xtÞ; Eðx; y; z; tÞ ¼ EðyÞ ð11Þ pffiffiffiffiffiffiffi ^ where i ¼ 1 and EðyÞ is a disturbance profile. With the mean flow available from previous computations (see Section 3), the stability problem consists in the determination of six unknowns, namely ar , ai , br , bi , xr and xi . These are the streamwise (x-direction) wavenumber and growth rate, spanwise (z-direction) wavenumber and growth rate, and wave frequency and temporal growth rate, respectively. In the present
0 Hz - 0.005 c 0 Hz - 0.006 c 0 Hz - 0.007 c 0 Hz - 0.008 c 0 Hz - 0.009 c 50 Hz - 0.007 c 50 Hz - 0.008 c 100 Hz - 0.005 c 100 Hz - 0.006 c 100 Hz - 0.007 c 100 Hz - 0.008 c 100 Hz - 0.009 c 150 Hz - 0.007 c 150 Hz - 0.008 c 150 Hz - 0.009 c
8.5 8
7 6.4 6
N-factor
4. Transition prediction method
Re = 1,920,000
9
5
4
3
2
1
0 0.0
0.1
0.2
0.3
1555
0.4
0.5
0.6
0.7
0.8
x/c Fig. 3. Final N -factor curves computed for Re ¼ 1; 920; 000.
1556
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
to Eq. (12), which constitutes a block-tridiagonal system of equations. For the present case, matrix Bij is invertible and a generalized matrix eigenvalue problem, solvable by LR methods, can be formulated. An inverse Rayleigh iteration procedure, together with a block LU factorization to invert the block-tridiagonal system, is further used to perform the local search (or improvement) of the eigenvalues. Linear stability theory provides information regarding the evolution of unstable disturbances in the boundary layer but, alone it is unable to predict the occurrence of transition. This task is accomplished by the use of the semi-empirical eN method where the stability characteristics of the flow are analyzed in terms of local, disturbance growth rates, which in turn are integrated to yield the N -factor. Thus, upon solving the eigenvalue problem and introducing Gaster’s transformation, the amplification rate for a specific disturbance in the form given by Eq. (11) can be obtained from the following expression:
method, the temporal formulation of the problem is preferred to the spatial analysis because the latter is more expensive from the computational point of view. Hence, a and b are real numbers and x is a complex number that is found by solving the stability equations. Temporal growth rates xi are afterwards converted into spatial growth rates ai by the straightforward application of a generalized Gaster’s transformation [23]. The (compressible) stability eigenvalue problem is solved here using a matrix finite-difference method to obtain local solutions when a good guess for the eigenvalue is available. When this is not the case, a globally convergent eigenvalue procedure is used instead. The formulation of the stability equations in matrix form leads to the expression: Aij Ej ¼ xBij Ej ;
ð12Þ
where x naturally stands for the eigenvalue and Ej is the discrete representation of the eigenfunction. The eigenvalue is calculated by the determinant condition applied
10
Re = 2,370,000 0 Hz - 0.004 c 0 Hz - 0.005 c 0 Hz - 0.006 c 0 Hz - 0.007 c 0 Hz - 0.008 c 100 Hz - 0.005 c 100 Hz - 0.006 c 100 Hz - 0.007 c 100 Hz - 0.008 c 150 Hz - 0.006 c 150 Hz - 0.007 c 200 Hz - 0.005 c 200 Hz - 0.006 c 200 Hz - 0.007 c 200 Hz - 0.008 c
9 8.5 8
7
N-factor
6.4 6
5
4
3
2
1
0 0.0
0.1
0.2
0.3
0.4
0.5
x/c Fig. 4. Final N -factor curves computed for Re ¼ 2; 370; 000.
0.6
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
A ¼ exp A0
Z
xi dxg ; Vg
parameter rðxr ; bÞ. Consequently, in order to compute the N-factor, some integration strategy must be followed. However, it must be noted that different integration strategies usually lead to distinct values of the N-factor as well, as shown e.g. by Pereira and Sousa [24]. Thus, the integration strategy should always display an unequivocal link with the physics of the dominant transition mechanism in the considered flow cases. In the present study, the chosen strategy can be described as following ‘‘waves’’ (or modes) of fixed fre2p quency xr and constant total wavelength k ¼ pffiffiffiffiffiffiffiffiffi . 2 2
ð13Þ
^ where A stands for the amplitude of EðyÞ at any particular value of the wall-normal coordinate y, the subscript 0 denotes the value where amplification started, Vg represents the norm of the group velocity vector with r oxr components ðox ; ob Þ and dxg is an infinitesimal arcoa length along the group velocity direction. Hence, for a wave of fixed frequency, Eq. (13) may be rewritten as: Z st A ¼ exp r ds; ð14Þ A0 s0
a þb
This method was first employed by Srokowski and Orzsag [25] for stationary modes only (physically associated with stationary cross-flow vortices). Using this strategy as well, travelling modes have been accounted for by Hefner and Bushnell [7]. Here, the N -factor was computed by maximizing the total amplification of modes characterized by a fixed frequency xr within the unstable range with respect to k, as given by the following expression:
where r is the characteristic growth rate of the disturbance, which is integrated from the neutral point s0 to the transition point st . Looking at Eq. (14), it seems that the computation of the N -factor is now similar to that in a two-dimensional case. However, the additional degree of freedom appearing in a three-dimensional problem is hidden in the
10
Re = 2,730,000 0 Hz - 0.003 c 0 Hz - 0.004 c 0 Hz - 0.005 c 0 Hz - 0.006 c 0 Hz - 0.007 c 100 Hz - 0.004 c 100 Hz - 0.005 c 100 Hz - 0.006 c 100 Hz - 0.007 c 150 Hz - 0.005 c 150 Hz - 0.006 c 200 Hz - 0.004 c 200 Hz - 0.005 c 200 Hz - 0.006 c 200 Hz - 0.007 c
9 8.5 8
7
N - factor
6.4 6
5
4
3
2
1
0 0.0
0.1
0.2
1557
0.3
0.4
0.5
x/c Fig. 5. Final N -factor curves computed for Re ¼ 2; 730; 000.
0.6
1558
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
10
Re = 3,730,000 0 Hz - 0.002 c 0 Hz - 0.003 c 0 Hz - 0.004 c 0 Hz - 0.005 c 0 Hz - 0.006 c 100 Hz - 0.003 c 100 Hz - 0.004 c 100 Hz - 0.005 c 200 Hz - 0.003 c 200 Hz - 0.004 c 200 Hz - 0.005 c 300 Hz - 0.003 c 300 Hz - 0.004 c 300 Hz - 0.005 c 400 Hz - 0.004 c
9 8.5 8
7
N-factor
6.4 6
5
4
3
2
1
0 0.0
0.1
0.2
0.3
0.4
x/c Fig. 6. Final N -factor curves computed for Re ¼ 3; 730; 000.
Table 2 Computed transition locations obtained for stationary and travelling disturbances
Both stationary and travelling modes are to be considered in this analysis as described in the next section.
Reynolds number: Re 106
Computed transition location: x=c Stationary: N ¼ 6:4
Travelling: N ¼ 8:5
5. Results and discussion
1.92 2.37 2.73 3.73
0.75 0.58 0.49 0.31
0.74 0.58 0.49 0.32
In this section, results from coupling Navier–Stokes computations of the mean flow around an infinite swept wing and eN transition predictions are presented and discussed. The selected test cases are constructed with the airfoil NLF(2)-0415 [26], which was designed as a low-drag wing to be used on a commuter aircraft with unswept wings. In the work of Dagenhart [13], detailed results from wind tunnel experiments of a 45 swept wing were obtained by operating the NLF(2)-0415 wing section at a small negative angle of attack a ¼ 4. At this condition, the experimental model was described as a nearly ideal cross-flow generator because the favorable pressure gradient on the upper surface keeps the Tollmien–Schlicthting instability subcritical back to 71% of chord for a 6 0. Attachment-line and G€ ortler insta-
" N ¼ max ln k
A A0
A A0
¼ exp ðx;kÞ
#
; ðx;kÞ
Z
st
r xr ; k ds;
ð15Þ
s0
where an asterisk denotes fixed values. Eq. (15) indicates that the computed N -factor to be used in the prediction of transition actually results from an envelope of curves.
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
1559
1.00 0.90
Experiments [13] Coupled method
0.80 0.70
x/c
0.60 0.50 0.40 0.30 0.20 0.10 0.00 1.5
Transition locations N = 6.4 for stationary disturbances
2.0
2.5
3.0
3.5
4.0
Re x 10-6
Fig. 7. Comparison between numerical predictions of transition and experiments.
bility mechanisms were not found for the investigated test conditions. A strong dependence of the transition location with the Reynolds number Re was observed in the experiments [13], due to the effect of this parameter in the strength of stationary cross-flow vortices. It was expected that the present numerical procedure would be able to reproduce the experimentally observed behavior. The Reynolds numbers and transition locations on the upper surface of the swept wing obtained from the published experiments (Naphtalene flow visualizations) for the investigated test cases are shown in Table 1. Fig. 2 shows the final pressure distributions (in terms of the non-dimensional pressure coefficient Cp ) obtained from the Navier–Stokes computations coupled with the eN method for transition prediction and covering all test cases indicated in Table 1. The results correspond to computations carried out on adapted numerical meshes. Fast convergence has been always achieved in the coupled procedure applied to the upper surface of the wing only. On the lower surface transition was fixed at the leading edge for simplicity. Experimental values obtained from the work of Dagenhart [13] are also displayed in Fig. 2 for comparison purposes (upper surface only). Excellent agreement between numerical calculations and experimental values for the (upper surface) lower end of the model can be observed in the graph. On the other hand, the pressure gradient on the upper end is slightly underpredicted by the computations. In the determination of the transition location resulting from the coupled procedure, two different approaches have been applied to the obtained envelopes of N -factor curves. The first considered stationary modes only in the linear, local stability analysis. In this case, the best correlation with the experimental transition loca-
tions was obtained and a value of the N -factor at transition of 6.4 applies. The second considered travelling modes in the analysis and therefore the maximization in Eq. (15) was extended to the frequency as well. In this case, a satisfactory correlation with the experimental transition locations was obtained for N-factor values of 8.5. These values of limiting N-factor are acceptable for experiments carried out in a wind tunnel environment. Figs. 3–6 present the final N-factor curves for each test case. The legends in the figures list the range of frequencies and wavelengths (expressed as a fraction of the chord c ¼ 1:83 m) of the disturbances employed in the analysis. Filled symbols have been used to denote stationary modes, whereas open symbols were used for travelling modes. Flow instability is detected very early in all cases ðx=c < 0:05Þ, which is in consonance with the development of cross-flow instability. The corresponding values of wave propagation direction (not shown) corroborated this assertion. A summary of the transition locations obtained by assuming the aforementioned values for the N-factor at transition in all test cases is also shown in Table 2. In addition, the excellent agreement between the results of the coupled procedure for stationary modes and the experimental data of Dagenhart [13] is confirmed by the transition plot portrayed in Fig. 7. Experimental error bars obtained from the same reference are also shown in this figure.
6. Conclusions In the present work, a demonstration of successful coupling between Navier–Stokes computations and eN method for transition prediction using linear, local stability analysis applied to an infinite swept wing has been
1560
J.M.M. Sousa, L.M.G. Silva / Computers and Structures 82 (2004) 1551–1560
presented. The procedures for the elliptic generation of adapted numerical meshes and to obtain second-order accurate Navier–Stokes solutions of the mean flow have been described as well. These constitute fundamental requirements for the success of the methodology. Four test cases have been constructed using the NLF(2)-0415 wing section at a sweep angle of 45 and varying the Reynolds number from 1.92 to 3.73 million. In all cases, the transition process was dominated by cross-flow instability. Additionally, the implemented methodology was able to reproduce the strong dependence of the transition location with the Reynolds number. The results have shown that, for the investigated test cases, excellent agreement with published experimental transition data can be obtained from the stability analysis of stationary modes only, by assuming a value of the N-factor at transition of 6.4. If travelling modes are also considered in the analysis, satisfactory agreement with experimental results may also be obtained by taking N ¼ 8:5 at transition instead.
Acknowledgements This research work has been partially supported by the European Commission through project ALTTA within the framework of the GROWTH programme.
References [1] Smith AMO, Gamberoni N. Transition, pressure gradient and stability theory. Report No. ES26388, Douglas Aircraft Co., 1956. [2] Van Ingen JL. A suggested semi-empirical method for the calculation of the boundary layer transition region. Report VTH-74, University of Delft, The Netherlands, 1956. [3] Arnal D. Boundary layer transition: predictions based on linear theory. AGARD Report No. 793. Special course on progress in transition modelling, 1994. [4] Joslin RD. Aircraft laminar flow control. Ann Rev Fluid Mech 1998;30:1–29. [5] Bushnell DM, Malik MR, Harvey WD. Transition prediction in external flows via linear stability theory. In: Zierep J, Oertel H, editors. IUTAM Symposium Transsonicum III. Berlin: Springer-Verlag; 1988. [6] Schrauf G. Large-scale laminar flow tests evaluated with linear stability theory. AIAA Paper No. 2001-2444, 2001. [7] Hefner JN, Bushnell DM. Status of linear boundary-layer stability theory and the eN method, with emphasis on swept-wing applications. NASA Technical Paper No. 1645, 1980.
[8] Morkovin MV. On the many faces of transition. In: Wells CS, editor. Viscous drag reduction. New York: Plenum Press; 1969. [9] Arnal D, Casalis G. Laminar-turbulent transition prediction in three-dimensional flows. Prog Aerospace Sci 2000; 36(2):173–91. [10] Stock HW. Airfoil validation using coupled Navier–Stokes and eN transition prediction methods. J Aircraft 2002; 39(1):51–8. [11] Stock HW, Haase W. A feasibility study of eN transition prediction in Navier–Stokes methods for two-dimensional airfoil computations. AIAA J 1999;37(10):1187–96. [12] Stock HW, Haase W. Navier–Stokes airfoil computations with eN transition prediction including transitional flow regions. AIAA J 2000;38(11):2059–66. [13] Dagenhart JR. Crossflow stability and transition experiments in a swept-wing flow. NASA Technical Memorandum No. 108650, 1992. [14] Thompson JF, Thames FC, Mastin CW. TOMCAT––A code for numerical generation of boundary-fitted curvilinear coordinate systems on fields containing any number of arbitrary two-dimensional bodies. J Comput Phys 1977; 24(3):274–302. [15] Kobayashi MH, Pereira JCF. Characteristic-based pressure correction at all speeds. AIAA J 1996;34(2):272–80. [16] Anderson DA, Tannehill JC, Pletcher RH. Computational fluid mechanics and heat transfer. New York: Hemisphere; 1984. [17] Fletcher CAJ. In: Computational techniques for fluid dynamics, vol. 1/2. Berlin: Springer-Verlag; 1988. [18] Hirsch C. In: Numerical computation of internal and external flows, vol. 1/2. NewYork: John Wiley; 1990. [19] Pereira JCF, Kobayashi MH, Marques NPC. Numerical solution of the turbulent flow over a fence using twoequation models. In: Dervieux A, Braza M, Dussauge JP, editors. Notes on numerical fluid mechanics, vol. 65. Braunschweig: Vieweg-Verlag; 1998. [20] Coakley TJ, Huang PG. Turbulence modelling for highspeed flows. AIAA Paper No. 92-0436, 1992. [21] Mack LM. Boundary-layer linear stability theory. AGARD Report No. 709. Special course on stability and transition of laminar flows, 1984. [22] Malik MR. Finite-difference solution of the compressible stability eigenvalue problem. NASA Contractor Report No. 3584, 1982. [23] Gaster M. A note on the relation between temporally increasing and spatially increasing disturbances in hydrodynamic stability. J Fluid Mech 1962;14:222–4. [24] Pereira JCF, Sousa JMM. Comparison of several N -factor integration strategies in a compressible boundary layer over a swept curved surface. ICAS Paper No. 94-2.4.2, 1994. [25] Srokowski A, Orszag SA. Mass flow requirements for LFC wing design. AIAA Paper No. 77-1222, 1977. [26] Somers DM, Horstmann KH. Design of a medium-speed natural-laminar-flow airfoil for commuter aircraft applications. DFVLR Report IB/29-85/26, 1985.