2.1 Coordinate system

The idealised geometry used by 2½D methods to model swept-tapered wings is illustrated in Fig. 2. The wing is assumed to be of constant section and twist, with leading and trailing edges meeting at a point O. All radii r joining the origin O with a point on the wing surface are generators of that surface. This conical symmetry is assumed to extend to the external inviscid flowfield such that isobars of surface pressure are coincident with the wing generators.

Figure 2. Schematic of a swept-tapered wing illustrating the line-of-flight section AB and the undeveloped computational arc AC for the boundary layer calculation.

Both the magnitude and the direction of the inviscid flow at the boundary layer edge are assumed to be independent of spanwise position. This allows the variation of the viscous flowfield with spanwise position to be reduced to a simple length scale dependence: in laminar boundary layer theory, the boundary layer thickness can be shown to vary with r ^{½(Reference Kaups and Cebeci10,Reference Horton and Stock11)} , while for turbulent boundary layers, the integral length scale varies more or less as r ^{(Reference Ashill and Smith12)}. To simplify reduction of the governing equations to infinite-yawed-wing form as the local taper radius r tends to infinity, a taper curvature κ_T is used.

2.2 Governing equations

Although based upon a well-established theory, the integral boundary layer equations within Callisto have been developed and refined over a number of years, starting with the entrainment method of Head^{(Reference Head1)}, which added to the classical von Kármán equation a transport equation for the boundary layer shape factor, H. Green^{(Reference Green13)} extended Head's method to compressible flows by means of a transformed shape factor H , defined in Equation (6), and then introduced the effects of turbulence history by adding an equation to account for the lag in the response of the turbulence to changes in the strain field^{(Reference Green, Weeks and Brooman3)}. At the same time, Smith^{(Reference Smith14)} extended the basic approach to three-dimensional flows by employing the streamline analogy, whereby turbulence characteristics in a 2D boundary layer are assumed to be replicated by the flow resolved the direction of the external streamline in a 3D boundary layer: closure of the increasing number of momentum thicknesses is achieved by assuming functional dependence of those length scales upon the streamwise momentum thickness, θ₁₁ , the transformed shape factor H , and the twist β of the limiting streamline (at the wall) relative to the direction of the inviscid flow. Smith's fully 3D approach was adapted for the special case of swept-tapered wings described above by Ashill and Smith^{(Reference Ashill and Smith12)}, to provide a low-cost approach to assessing the effect of sweep and taper on turbulent boundary layers on transonic wings. The advantages of lower-cost, approximate methods are still evident today where they can be enablers for CPU-intensive optimization strategies.

The governing equations in the Callisto code are derived from the system of three integral equations for momentum, resolved parallel and normal to the inviscid flow streamline, and for entrainment, as presented by Ashill and Smith^{(Reference Ashill and Smith12)}. The latter system of equations is presented in the appendix, as are the principal modifications to the notation and to the choice of independent variables adopted for the Callisto method. The resulting system, expressed in the conical-surface (x,y) coordinate system illustrated in Fig. 2, can be written:

(1) $$\begin{equation} \left( {\begin{array}{@{}*{3}{c}@{}} {{A_\theta }}&\ {{A_\beta }}&\ {{A_{\bar{H}}}}\\ {{N_\theta }}&{{N_\beta }}&{{N_{\bar{H}}}}\\ {{E_\theta }}&{{E_\beta }}&{{E_{\bar{H}}}} \end{array}} \right)\frac{d}{{dx}}\left( {\begin{array}{@{}*{1}{c}@{}} {{\theta _{11}}}\\ {\tan \beta }\\ {\bar{H}} \end{array}} \right) = \;\left( {\begin{array}{@{}*{1}{c}@{}} {{\lambda _1}}\\ {{\lambda _2}}\\ {{\lambda _3}} \end{array}} \right) \end{equation}$$

and

(2) $$\begin{equation} \frac{{\partial {\theta _{11}}}}{{\partial y}} = - \frac{{{\theta _{11}}}}{r}\quad \left( {{\rm{for\ turbulent\ flows}}} \right)\\ \end{equation}$$

(3) $$\begin{equation} \frac{{\partial \beta }}{{\partial y}}\; = \;\frac{{\partial \bar{H}}}{{\partial y}}\; = \;0 \end{equation}$$

The non-trivial coefficients A_m and λʹ_n are defined in the appendix. The independent variables in the above equations relate to the 3D boundary layer profiles as follows:

(4) $$\begin{equation} {\theta _{11}} = \int_0^\delta {\frac{{\rho u}}{{{\rho _e}{U_e}}}\left( {1 - \frac{u}{{{U_e}}}} \right)dz} ,\\ \end{equation}$$

(5) $$\begin{equation} \tan \beta = {\left. {\frac{{{\raise0.7ex\hbox{${\partial v}$} \! / \!\lower0.7ex\hbox{${\partial z}$}}}}{{{\raise0.7ex\hbox{${\partial u}$} \! / \!\lower0.7ex\hbox{${\partial z}$}}}}} \right|_{z = 0}},\\ \end{equation}$$

(6) $$\begin{equation} \bar{H} = \frac{1}{{{\theta _{11}}}}\int_0^\delta {\frac{\rho }{{{\rho _e}}}\left( {1 - \frac{u}{{{U_e}}}} \right)dz} \end{equation}$$

For completeness, we note that the system of equations above can be solved only for attached flows and that, for flows near to separation or mildly separated, the inviscid velocity gradient becomes an independent variable (rather than a constant factor in the coefficients λ ₁, λ ₂ and λ ₃) governed by an additional equation for the wall-normal velocity at the edge of the boundary layer. The principles of the approach are well described by Lock and Williams^{(Reference Lock and Williams4)}; however for the purposes of this paper, the basic three-equation system is sufficient to calculate the development of a turbulent boundary layer near the AL. We also note here that the entrainment equation can be augmented by Green's lag equation but again, as it is not central to the behaviour of the three-equation system at the AL, it will not be discussed in this paper. One further complication is that the governing equations in Callisto are transformed into a non-orthogonal coordinate system, (ξ,η) in Fig. 2. However, the additional terms do not alter the numerical challenge at the leading edge and so again are not introduced in this paper. Rearranging Equation (1) above, and substituting the expressions listed in the appendix, we obtain the first-order, non-linear differential system:

(7) $$\begin{equation*} \frac{d}{{dx}}\left( {\begin{array}{@{}*{1}{c}@{}} {{\theta _{11}}}\\ {\tan \beta }\\ {\bar{H}} \end{array}} \right) = {\left( {\begin{array}{@{}*{3}{c}@{}} {\cos \psi - \sigma \,{f_2}}&{ - \sin \psi \,{f_2}{\theta _{11}}}&{ - \sigma \,\frac{{d{f_2}}}{{d\bar{H}}}{\theta _{11}}}\\[6pt] {\tan \beta \left( {\cos \psi \,{f_1} - \sigma \,{f_4}} \right)}&\quad {\left( {\cos \psi \,{f_1} - 2\sigma \,{f_4}} \right)\,{\theta _{11}}}&\quad {\tan \beta \left( {\cos \psi \frac{{d{f_1}}}{{d\bar{H}}} - \sigma \,\frac{{d{f_4}}}{{d\bar{H}}}} \right){\theta _{11}}}\\[6pt] {\cos \psi \,{H_1} + \sigma \,{f_3}}&{\sin \psi \,{f_3}{\theta _{11}}}&{\left( {\cos \psi \frac{{d{H_1}}}{{d\bar{H}}} + \sigma \,\frac{{d{f_3}}}{{d\bar{H}}}} \right){\theta _{11}}} \end{array}} \right)^{ - 1}}\;\left( {\begin{array}{@{}*{1}{c}@{}} {{\lambda _1}}\\ {{\lambda _2}}\\ {{\lambda _3}} \end{array}} \right), \end{equation*}$$

where we have introduced the shorthand

(8) $$\begin{equation} \sigma = \sin \psi \tan \beta \end{equation}$$

This type of initial-value problem is well suited to solution by a Runge-Kutta scheme: the variant due to England^{(Reference England15)}, with a facility to reduce the marching step size if the estimated truncation error exceeds a pre-defined threshold, is used in Callisto.

2.3 Conditions for singularity of the viscous derivatives matrix

There are two conditions where the system of Equation (7) can clearly be seen to become singular. The first is most clearly evident in case of 2D flow where tanβ ≡ 0, cosψ ≡ 1 and sinψ ≡ 0 for all x (so that σ ≡ 0), when the normal momentum equation can be de-coupled from the rest of the system, leaving the following equation:

(9) $$\begin{equation} \frac{d}{{dx}}\left( {\begin{array}{@{}*{1}{c}@{}} {{\theta _{11}}}\\ {\bar{H}} \end{array}} \right) = {\left( {\begin{array}{@{}*{2}{c}@{}} 1&0\\ {{H_1}}&\quad{\frac{{d{H_1}}}{{d\bar{H}}}{\theta _{11}}} \end{array}} \right)^{ - 1}}\;\left( {\begin{array}{@{}*{1}{c}@{}} {{\lambda _1}}\\ {{\lambda _3}} \end{array}} \right) \end{equation}$$

The singularity occurs when the parameter $\frac{{d{H_1}}}{{d\bar{H}}}$ in Equation (9) becomes zero (at separation). This condition is well known and is dealt with by introducing a transpiration equation to the system and by changing the dependent variables, as discussed above.

The second singular case is AL flow, as illustrated in Fig. 3 below: locally (at A-A) both cosψ = 0 and tanβ = 0 so that σ = 0, leaving all but two terms in the matrix (7) equal to zero.

Figure 3. The inviscid flow vectors in the vicinity of the AL.

This AL singularity was tackled by Cumpsty and Head^{(Reference Cumpsty and Head16)} and later by Smith^{(Reference Smith17)} by differentiating the spanwise momentum equation with respect to the x coordinate to obtain the following non-linear simultaneous equations for the flow gradients at the AL:

(10) $$\begin{equation} \frac{1}{{{V_1}}}\frac{{d{U_1}}}{{dx}}{\theta _{11}} - {f_2}\,{\theta _{11}}\frac{{d\beta }}{{dx}} = \frac{{{c_f}}}{2}, \end{equation}$$

(11) $$\begin{equation} 2{f_4}{\left( {{\theta _{11}}\frac{{d\beta }}{{dx}}} \right)^2} + \left( {\frac{{{c_f}}}{2} - \frac{{3{f_1}}}{{{V_1}}}\frac{{d{U_1}}}{{dx}}{\theta _{11}}} \right)\left( {{\theta _{11}}\frac{{d\beta }}{{dx}}} \right) + \left( {H + 1} \right)\;{\left( {\frac{1}{{{V_1}}}\frac{{d{U_1}}}{{dx}}{\theta _{11}}} \right)^2} = 0, \end{equation}$$

(12) $$\begin{equation} \frac{{{H_1}}}{{{V_1}}}\frac{{d{U_1}}}{{dx}}{\theta _{11}} + {f_3}{\theta _{11}}\frac{{d\beta }}{{dx}} = {c_E} \end{equation}$$

Note that U ₁ and V ₁ (see appendix) are re-introduced in the above equations, rather than ψ, for clarity. These equations can be solved iteratively but straightforwardly to obtain the momentum thickness and rate-of-change of limiting streamline angle at the AL (the rate of change of U ₁ is a boundary condition given by the inviscid flowfield). The AL solution provides the initial values for Equation (7) which, in principle, can be solved at any point downstream of the AL, given the finite rate of change of both U ₁ (and therefore ψ) and β at the AL. It turns out that, depending upon the resolution of the computational domain very near the AL, Equation (7) is ill-posed until a short distance downstream of the AL – the experience of the authors has been that the system of integral equations cannot be robustly tackled with a Runge-Kutta scheme until the angle ψ has decreased to about 80°. Similar difficulties were reported by Thompson and McDonald^{(Reference Thompson and McDonald18)} and later by Smith^{(Reference Smith17)}, who adopted a numerical fix similar to the one employed by Callisto and described below. The issue may not have been encountered more widely because the problem region accounts for less than 1% of the chord length and would not impact upon less finely resolved computations. However, standards of resolution for the analysis of leading-edge flows have tightened in recent years, owing to the need to resolve the development of crossflow instability on laminar leading edges: adoption of the same numerical domain for the analysis of swept wings with turbulent AL invariably throws up this numerical feature.

2.4 Leading Edge Approximation

On the basis that the rates of change of both ψ and β are small for some distance downstream of the AL, a numerical fix was devised, involving a simple extrapolation of the AL similarity solution, used by Cumpsty and Head^{(Reference Cumpsty and Head16)} and Smith^{(Reference Smith17)} over the region where the integral equations were ill-posed, as defined by the criterion ψ > 80°. Such a sacrifice of accuracy for robustness was justified on the basis of the very short region concerned and the strongly favourable pressure gradient over this region, leading to the belief that the downstream development of the boundary layer integral quantities would not be significantly affected, hence providing acceptable accuracy for the calculation of the profile drag. This assumption appears to have been justified, but a fuller discussion is left until later in the paper, along with the demonstration of the effect of this numerical fix on the solution of the boundary layer equations. At this point we simply observe that there were no experimental data available to validate this approximation – hence the experimental study reported below.

3.1 Experimental set-up and procedures

The experimental AL was generated on NACA0050 aerofoil profile, at zero angle of incidence, swept back by 60°. The model was made of wood with a normal-to-LE chord length of 466mm and a LE radius of curvature of 114mm. It was mounted between the floor and ceiling of the test section of the T2 wind tunnel in the Handley Page Laboratory at City, University of London, as shown in Fig. 4.

Figure 4. Schematic of the swept-wing model mounted in the T2 wind tunnel.

The T2 tunnel has a speed range of 4–55m/s. The test section measures 1120mm wide by 800mm high by 1.2m long: the model extended into the tunnel diffuser by 150mm. The blockage ratio introduced by such a large model is mitigated to an extent by the high sweep angle. Because the present investigation is focused on the flow very near the swept leading edge, the necessary numerical adjustments to the accompanying calculations can be achieved by an ‘effective sweep’ correction which infers the true magnitude of the velocity component parallel to the attachment line based upon the measured value of the pressure coefficient at the attachment line. The reader is referred to Ref. Reference Gowree, Atkin and Grupppetta9 for a comparison of the chordwise pressure distributions obtained from parallel rows of pressure taps around the model near the root (bottom), mid-span and tip sections. These pressure distributions and subsequent flow visualisations suggest that the flow separated at about 80% chord over most of the span of the model.

A surface-mounted boundary layer traverse probe with micro-displacement capability was employed to capture the velocity profile, with a resolution of 2.5μm per step achievable within an accuracy of ±0.09μm. This traverse, and the optical system used to calibrate probe displacement from the surface, is described in detail by Gowree et al^{(Reference Gowree, Atkin and Grupppetta9)} for measurements on the attachment line (almost exactly on the model centreline). The right-hand image in Fig. 5 shows the additions made to the traverse, still fixed to the centreline line of the model, to allow the hot-wire probe to survey the flow away from the attachment line. The left-hand image illustrates how the axis of the optical alignment system was inclined to enable the same technique as described in Ref. Reference Gowree, Atkin and Grupppetta9 to be used for determining distance between probe and model surface.

Figure 5. Arrangement used to apply the techniques used by Gowree et al^{(Reference Gowree, Atkin and Grupppetta9)} to probe the flow away from the leading edge of the model.

A single normal (SN) hot-wire probe, Dantec 55P15, was used to capture the single velocity component at the AL and a Single Yawed (SY) probe, 55P12, for the measurement of the two in-plane velocity components downstream of the AL. The Constant Temperature Anemometry (CTA) technique was used and the hot-wire probes were connected to the DISA-55M10 CTA Standard Bridge (M-Unit) module. The M-Unit was in turn interfaced with a National Instruments (NI)-DAQ card with built-in A/D converter and installed in a PC for data acquisition using NI LabVIEW. The hot-wire output signal, sampled for 10 seconds at 10 kHz, was pre-filtered through a low-pass filter rated at 4.8kHz prior to recording. King's law was applied for the reduction of hot-wire output voltage. The calibration of the SN probe and measurement of a single velocity component at the AL was straightforward, but more challenging for the SY probe used to measure the two velocity components immediately downstream of the AL. In this case, Bradshaw's method described in Ref. Reference Gowree19 was employed and a yaw calibration was required due to the directional sensitivity of the SY probe. The estimated error for the hot-wire measurements was less than 5%^{(Reference Gowree19)} but the excellent agreement between the measured velocity profiles along a laminar attachment line and the theoretical solution for swept Hiemenz flow (Fig 17 of Ref. Reference Gowree, Atkin and Grupppetta9) suggests that the actual error, including probe interference, was marginal.

Preston's^{(Reference Preston20)} technique was employed for the measurement of local surface shear stress. At the AL, the flow resolves into a single spanwise velocity component, similar to the streamwise flow along a flat plate, thus the method was expected to yield reasonable accuracy. This technique has thus far been restricted to 2D flows where the skin friction acts in the same direction as the velocity component; however, for the flow downstream of the AL, an attempt was made to extend this technique to the 3D boundary layer under the highly curved streamline at the LE of a swept wing. The surface shear stress measurement was made by aligning the pitot tube in the direction of the local external streamline, obtained from the velocity components measured using the CTA technique at the edge of the boundary layer. The angle between the external and the surface streamlines was small (approximately ±5^° as discussed in Section 4.2), and therefore this approach was felt to add negligible additional error: the uncertainty analysis conducted by Gowree^{(Reference Gowree19)} yielded the same confidence level as for the 2D case of purely spanwise flow along the attachment line.

5.1 Modification to leading edge model

As discussed previously, in Callisto a numerical approximation to the governing equations is imposed near the AL, specifically while ψ > 80°. For the current experiment this zone ends somewhere in the region 0.005 < x/c < 0.01, as shown in Fig. 11. As discussed at the end of Section 4.2, in this region the limiting streamline angle appears to be small in magnitude, |β| < 3°. This experimental observation suggests that the governing equations might remain ill-posed immediately downstream of the AL not simply because cosψ ≈ 0, but because tanβ is also small in the same region: the slow development of cosψ and tanβ means that the system of Equations (7) remains dominated by just two terms out of nine close to the attachment line (see discussion in Section 2.3). This results in numerical error and poor robustness when marching the boundary layer equations using a Runge-Kutta scheme. A new, more subtle numerical approximation to the governing equations was therefore developed: to extend the assumption – by Cumpsty and Head^{(Reference Cumpsty and Head16)} and Smith^{(Reference Cumpsty and Head16)} – of constant, finite ∂β/∂x at the AL to the flow immediately downstream. This allows the normal momentum equation to be de-coupled, allowing the remaining two equations to be solved as a coupled system, as shown below:

(18) $$\begin{eqnarray} \frac{d}{{dx}}\left( {\begin{array}{@{}*{1}{c}@{}} {{\theta _{11}}}\\ {\bar{H}} \end{array}} \right) = {\left( {\begin{array}{@{}*{2}{c}@{}} {\cos \psi - \sigma \,{f_2}}&{ - \sigma \,\frac{{d{f_2}}}{{d\bar{H}}}{\theta _{11}}}\\ {\cos \psi \,{H_1} + \sigma \,{f_3}}&\quad{\left( {\cos \psi \frac{{d{H_1}}}{{d\bar{H}}} - \sigma \,\frac{{d{f_2}}}{{d\bar{H}}}} \right){\theta _{11}}} \end{array}} \right)^{ - 1}}\!\left( {\begin{array}{@{}*{1}{c}@{}} {{\lambda _1} + \sin \psi \,{f_2}{\theta _{11}}{{\left. {\frac{{d\beta }}{{dx}}} \right|}_{AL}}}\\ {{\lambda _3} - \sin \psi \,{f_3}{\theta _{11}}{{\left. {\frac{{d\beta }}{{dx}}} \right|}_{AL}}} \end{array}} \right)\!,\nonumber\\ \end{eqnarray}$$

(19) $$\begin{equation} \frac{{d\beta }}{{dx}} = {\left. {\frac{{d\beta }}{{dx}}} \right|_{AL}} \end{equation}$$

The experimental results actually show that ∂β/∂x does vary downstream of the attachment line, so Equation (20) can be thought of as the first term of an expansion for ∂β/∂x. In this respect we are primarily looking for a more robust numerical treatment rather than for a model of the highest accuracy. Moreover no attempt has been made (nor could be, given the sparseness of the crossflow profile data obtained) to improve the Mager^{(Reference Preston25)} crossflow profile modelling (see Appendix) which – as is well-known – cannot cater for the type of crossover profiles seen in the experimental results.

The revised solution process is then as follows: the first task is to solve the simultaneous system described by Equations (10) through (12) for θ ₁₁, H (and therefore H ) and ∂β / ∂x at the AL; immediately downstream of the AL, Equations (18) and (19) are used to march the solution until tan β recovers to a threshold value, whereupon the full 3D system of Equations (7) can be restored. This threshold value is determined by numerical robustness and is so small that the switch from Equation (18) to the full 3D form – Equation (7) – is hardly detectable in the results.