Stochastic optimisation based control of boundary layer transition

Control Engineering Practice 10 (2002) 243–260

Stochastic optimisation based control of boundary layer transition W. MacCormacka, O.R. Tuttyb,*, E. Rogersc, P.A. Nelsona a

Institute of Sound and Vibration Research, University of Southampton, Highfield, Southampton SO17 1BJ, UK b School of Engineering Sciences, University of Southampton, Highfield, Southampton SO17 1BJ, UK c Department of Electronics and Computer Science, University of Southampton, Highfield, Southampton SO17 1BJ, UK Received 18 May 2000; accepted 27 August 2001

Abstract Suction is one of the most promising techniques for delaying the transition from laminar to turbulent flow and hence reducing the drag force acting on an aircraft. However, in order to achieve an overall reduction in energy consumption and thereby operating costs, it is necessary to apply the suction in an efficient if not optimal manner. An investigation is conducted into the use of distributed surface suction with multiple suction panels where the suction distribution is optimised to satisfy a desired objective. Cost functions based on minimising the suction effort while maintaining transition in a fixed position and minimising the total energy consumption of the system are employed. Evolutionary methods, genetic algorithms (GA) and simulated annealing (SA), are used to perform the optimisation. These methods were successful in cases in which gradient based search techniques fail. Physically sensible suction distributions were produced. In general, better behaviour was found with the GA than with SA, which was prone to premature convergence. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Optimisation; Genetic algorithms; Simulated annealing; Controller design; Flow control; Transition control

1. Introduction It is well established that suction applied on the surface of a body can delay transition from laminar to turbulent flow. Since drag is significantly lower with laminar flow, suction applied to e.g. the fin, wings or nacelles of an aircraft has the potential to reduce drag. However, the ultimate objective, minimisation of the total net drag, will be achieved only if the energy saved in the drag is greater than that required to drive the suction system. This paper presents recent results on the optimisation of suction systems for hybrid laminar flow control. The basic approach used in the research programme from which this paper derives consists, in effect, of monitoring the state of the flow together with the automatic control of suction applied through the surface of the body. This programme involves both algorithm development and verification and wind tunnel based experiments. In the case of the latter, an essential *Corresponding author. Tel.: +44-0-2380-592-163; fax: +44-02380-593-058. E-mail address: [email protected] (O.R. Tutty).

element has been the development of techniques for monitoring the position of transition in a boundary layer. Experimentally, pressure measurements are used to determine the position of transition in a boundary layer, and a feedback control loop is then used to maintain transition at a desired location with minimum power consumption by controlling the suction flow rates on suction panels. Previous work has demonstrated that this overall strategy is technically feasible both in simulation studies (of which computational fluid dynamics is an essential part) (Tutty, Hackenberg, & Nelson 2000a, b) and by wind tunnel based experiments (Nelson, Wright, & Rioual, 1997; Wright & Nelson, 1999). This work employed gradient based controllers, which used a local linear model constructed by standard linear system identification tools around selected operating points. For the case of flow over a flat plate these controllers performed well. In somewhat more complicated cases, however, they did not, e.g. when a non-zero pressure gradient is applied to a flat plate (Tutty et al., 2000b). In this paper, it is shown that the gradient based approach fails because it can be incompatible with the basic physics of the flow. This leads to the introduction

0967-0661/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 7 - 0 6 6 1 ( 0 1 ) 0 0 1 4 0 - X

244

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

of a new cost function for which no analytic (one step) solution algorithms exist. Instead, iterative search techniques must be used. Here the use of stochastic optimisation techniques is considered. In particular, genetic algorithms (GA) and simulated annealing (SA) are used and the performance of the resulting control/ optimisation schemes analysed.

2. The basic problem Consider a flat plate with a sharp leading edge aligned with a uniform incoming flow. Over most of the flow field, the effect of the viscosity of the fluid is negligible, and the flow will remain uniform. However, near the plate, the velocity must decrease from the free stream velocity to zero at the surface of the plate, and in this region, the ‘boundary layer’, viscous effects play a significant role. The width of the boundary layer is small compared to the length of the plate, but it increases along the plate, as x1=2 where x is the distance along the plate from the leading edge, with the gradient of the velocity decreasing as x
1=2 : Suppose now that small disturbances are introduced into the flow in the boundary layer. Near the leading edge, for xox0 ; where x0 E6  104 n=UN

ð1Þ

and UN and n are the free stream velocity and kinematic viscosity of the fluid, respectively, these disturbances decay in amplitude, and the flow is stable. However, for x > x0 the disturbances will grow, and the flow is unstable. These disturbances, known as Tollmien– Schlichting waves, initially grow relatively slowly, but further downstream will provoke much more vigorous disturbances which rapidly lead to transition and a fully turbulent flow. One of the primary characteristics of turbulent flow is enhanced mixing due to the rotational character of the small scale motions in the flow. In particular, in a turbulent boundary layer, the turbulence transports momentum from (the higher mean velocity) outer flow towards the wall and increases the velocity near the wall compared to that found in a laminar (non-turbulent) flow in similar conditions. This increases the velocity gradient at the wall, and since the shear stress (the frictional force) is proportional to the velocity gradient, the shear stress increases drastically, by at least an order of magnitude, at transition. The scenario presented above is for a flat plate, but will apply to any body which is basically aligned with the flow. Consequently, if the flow remains laminar, the drag on the body due to the skin friction will be much less than that if the flow is turbulent. Hence, there is great interest in the aerospace industry in delaying transition, or modifying the turbulent flow once transi-

tion has occurred, to decrease drag, and hence the fuel consumption of civil aircraft. This paper is concerned with practical methods of delaying transition. There are a number of ways by which this can be achieved. For example, heating the surface will heat the fluid adjacent to it, and increase its kinematic viscosity, and enhance the stability of the flow: from (1) it can be seen that increasing n will move the point that the flow becomes unstable downstream. However, in the aerospace industry, there is a great deal of interest in using suction to delay transition. It has long been known that, in theory, withdrawing very small amounts of fluid from the boundary layer through the wall can greatly enhance the stability characteristics of the flow. In practice, suction flow rates are at least an order of magnitude greater than those suggested by the basic theory for a flat plate. However, they are still small, and suction is regarded as one of the most viable means for practical drag reduction. This paper considers distributed suction using discrete suction panels of the type commonly used in the aerospace field, with essentially uniform suction on each panel, but with the ability to vary the suction flow rate on each panel independently to optimise the effort required to achieve a desired objective. Different objectives and hence cost functions may be specified in different situations. For an aerofoil the aim will not be to simply delay transition as much as possible, but to minimise the energy consumption of the system taking into account a number of factors. On a nacelle, however, the design requirement is to fix transition at a particular point so that the benefits of laminar flow are obtained ahead of this point but the flow at the trailing edge of the nacelle is turbulent. Turbulence is desirable at the trailing edge of a nacelle as the contribution to the total drag from the pressure rather than the skin friction is substantially decreased in a body with a large blunt trailing edge if the flow is turbulent in this region. Both of these design requirements will be considered below. Note that the transition scenario outlined above, through the amplification of linear Tollmien–Schlichting waves is not the only way that transition can occur. In particular, in a ‘noisy’ environment, a ‘bypass’ mechanism may operate where this linear stage may be absent completely, and transition may occur much more rapidly. However, in an environment with a low free stream turbulence level, there should be a relatively long region in which this linear mechanism is important, and where surface suction as considered below can have a significant effect on the flow.

3. Experimental The experimental arrangement is shown in Fig. 1. A flat plate with, if required, a smooth two-dimensional

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

hump on the wall opposite to generate a pressure gradient, and two suction panels, one upstream and one downstream of the peak of the hump, is positioned in a wind tunnel with mean flow speeds upto 22 m=s: The suction panels are made of laser drilled titanium sheets, which provides a surface which is sufficiently smooth so as not to provoke transition through surface roughness, but a porosity (hole area to surface area) of around 0.5%. Suction flow rates of up to 2 cm=s are possible, where the suction flow rate is defined as the volume flow per unit surface area, and hence the flow in the holes is of the order of the suction flow rate divided by the porosity. The humps are circular arcs 0:36 m long with a radius of 0:555 m; giving a blockage of 10.3%. The suction panels are located so that, according to hot wire

245

anemometry experiments, without suction, transition occurs on the second panel. Each suction panel is connected to a pump which is controlled by a PC via an inverter. Downstream of the two panels an array of flush-mounted microphones is installed. The signals of these microphones are high pass filtered at 800 Hz in order to remove the background noise due to the wind tunnel fan. The microphone signals are sampled at a sampling rate of 4 kHz for a given period, then the rms value is calculated. These rms values are then normalised by pre-measured rms values of the signals for a fully turbulent boundary layer, hence yielding a value close to 0 for laminar flow and a value close to 1 for turbulent flow. The location of the laminar–turbulent transition is determined by comparison of the rms pressures acquired by the microphones to four reference pressures. The resulting error signal can be described by eðkÞ ¼

4 X

½pref ðxm Þ
pðxm ; kÞ

m¼1

¼

4 X

pref ðxm Þ


m¼1

4 X

pðxm ; kÞ ¼ r
yðkÞ;

ð2Þ

m¼1

where k is the iteration index, r is the sum of the desired normalised pressures, and yðkÞ is the sum of the measured normalised pressures. The desired normalised pressures are produced when transition is in the specified location. For the experiments, r and yðkÞ are referred to as the desired plant output and the actual plant output, respectively. These two values are also directly related to the desired transition location d and the actual transition location xT ; such that (2) is equivalent to eðkÞp½d
xT ðkÞ:

ð3Þ

As illustrated in Fig. 2, an appropriate choice of the reference pressures enables the specification of the

Fig. 1. Experimental set-up.

Normalised rms pressure Transition at desired location 1 Transition upstream of desired location

Transition downstream of desired location

x

1

x

2

x

3

x

4

Microphone positions

Fig. 2. Monitoring the transition position using surface mounted microphones.

x

246

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

boundary-layer equations: u*

@u* @u* du* e @2 u* þ v* ¼ u* e þ ; @x* @y* dx* @y* 2

ð4Þ

@u* @*v þ ¼0 @x* @y*

ð5Þ

with the following boundary and initial conditions: y* ¼ 0: y-N: *

u* ¼ 0;

v* ¼ v*wall ;

uð * x; * yÞ* u* e ðxÞ; * 1 * ¼ u* e0 ðxÞ * þ pffiffiffiffi u* e ðxÞ p R

x* ¼ x* 0 : Blasius profile;

Fig. 3. Experimental relationship between transition position and suction flow rates for a two-panel-flat plate case.

location of the transition region with respect to the microphone positions. Fig. 3 shows, for a constant mean flow speed, the position of transition as a function of the suction flow rates on the two panels. The contour lines indicate the position of transition. Fig. 3 also illustrates the purpose of the controller which will regulate the two suction flow rates until the shaded point is reached. This point corresponds to zero error with the minimum possible suction effort.

4. Theoretical modelling When simulating viscous laminar flow, it is possible to solve the full governing equations, the Navier–Stokes equations, numerically. However, with the type of high Reynolds numbers considered here, it is more efficient to apply Prandtl’s theory which assumes that the viscous part of the flow is restricted to a thin boundary layer, and solve a reduced set of equations which contain all the essential physics of the problem. Carter and Wornom (1975) suggested a method of coupling the viscous and inviscid parts of the flow field, where both are calculated from a combination of the axial velocity at the edge of the boundary layer and the displacement thickness. Veldman (1981) developed a numerical algorithm to apply this strategy in which the viscous part of the flow is described by the normalised

ð6Þ ð7Þ Z

x* N x* 0

% @ðu* e d* Þ=@x* * dx; ð8Þ x*
x*

ð9Þ

where u* and v* are the axial and the transverse velocities, respectively, u* e is the axial velocity at the edge of the boundary layer, and x* and y* are the streamwise and normal coordinates, respectively. The variables in (4)–(9) have been non-dimensionalised using the scalings u v pffiffiffiffi x y pffiffiffiffi u* ¼ ; v* ¼ R; R; x* ¼ ; y* ¼ UN UN L L where UN is the free stream velocity, L is the length of the plate, and R ¼ rUN L=m is the free stream Reynolds number, with r the density and m the viscosity of the fluid. The first part of the edge velocity u* e0 in (8) comes from a first order potential flow calculation for the flow past the body and the second part is a correction due to the viscous displacement of the streamlines, where the % displacement thickness is denoted by d* : The computational region is covered by a grid which is equidistant in axial direction and stretched in the vertical direction, so that the majority of the points are placed close to the wall. Eqs. (4) and (5) are rewritten using the axial velocity and the stream function as independent variables. All derivatives are then approximated with second order accuracy, using centred differences for derivatives in the ydirection and three-point-backward differences for derivatives in the x-direction. The finite difference equations are solved at each streamwise position using Newton’s method, and the edge velocity is determined by performing Gauss–Seidel sweeps over the flow field. The effect of the interaction of the boundary layer with the outer inviscid flow is most pronounced at the ends of the suction panels, where the boundary condition at the wall is discontinuous, and the solution with interaction is much smoother than that without interaction. It is assumed that the transition from laminar to turbulent flow is preceded by the linear amplification of small two-dimensional travelling waves, the so-called

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

Tollmien–Schlichting waves, rather than through a bypass mechanism. Linear hydrodynamic stability theory describes the stability of the flow to this kind of disturbances. A detailed description of the linear stability theory is given, for example, by Mack (1984); therefore, only the main features are discussed in this paper. Linear stability theory is based on the study of small perturbations superimposed on a basic laminar flow. The purpose of this theory is to determine the conditions under which these perturbations are amplified or damped. Therefore, the velocity components and the pressure are decomposed into an undisturbed mean part and a small disturbance. This decomposition is introduced into the unsteady two-dimensional Navier– Stokes equations. The resulting equations are simplified by making the approximation that the products of disturbance quantities are negligible. This results in a set of linear disturbance equations. The pressure is then eliminated by cross-differentiation. In the resulting equation, the disturbance is expressed as a normal mode: jðx; y; tÞ ¼ FðyÞeiðax
otÞ ;

ð10Þ

where j is the (dimensionless) disturbance stream function, a is the wave number and o the wave frequency of the disturbance, and F is the (complex) disturbance amplitude. This leads to the following disturbance velocities: u# ¼ F0 eiðax
otÞ ;

ð11Þ

v# ¼
iaFeiðax
otÞ ;

ð12Þ

where the prime stands for differentiation with respect to y: Introducing these velocities into the disturbance equation yields   o 00 u
F
a2 F
u00 F a i ðF0000
2a2 F00 þ a4 FÞ ¼ 0; ð13Þ þ a Re the well-known Orr–Sommerfeld equation. The Reynolds number Re is now based on the local boundarylayer thickness, and u is the streamwise component of the mean-flow velocity, obtained by the solution of the boundary-layer equations. Note that there is an implicit rescaling here between the variables used for the flow calculation and the stability calculation. Eq. (13) presents an eigenvalue problem with homogeneous boundary conditions (u# and v# must vanish at the wall and in the free stream). For a given Reynolds number and real frequency, the solution yields an eigenfunction in form of the complex disturbance amplitude and as an eigenvalue, the complex wave number a ¼ ar þ iai : The imaginary part of the wave number, ai gives the spatial growth rate of the disturbance, which is positive for damped waves and negative for amplified waves. The speed of propagation

247

of the disturbance is given by ar =o: Using linear stability theory, the onset of instability can be predicted, but instability is only a very early precursor of the ultimate transition to turbulent flow. The eN -method is a standard method of relating the point of transition to the amplification of the disturbance waves. The total amplification rate is calculated by integrating the spatial amplification factors ai at each single real frequency o:

Z x A ln ai dx; ð14Þ ¼
A0 x0 where A is the wave amplitude and the index 0 refers to the streamwise position where this wave becomes unstable. The eN -method is based on the envelope of these curves of total amplification, where N is the maximum amplification factor at each point.  A N ¼ max ln : ð15Þ A0 The position of transition is then predicted as the first point at which N exceeds a critical value NT which can vary according to the operating conditions of the flow. Malik (1990) gives an extensive overview over the values of NT found in various wind tunnel and flight experiments. For a low free stream turbulence level, NT ¼ 10 proves to be a reasonable approximation. However, Mack (1977) shows that this critical amplification factor is strongly dependent on the turbulence level. He suggested the following relation: NT ¼
8:43
2:4 lne Tu

ð16Þ

where the turbulence level is calculated as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u02 þ v02 þ w02 ; Tu ¼ 3 with ðu0 ; v0 ; w0 Þ being the dimensionless turbulent disturbance velocities, and the overline represents averaging with respect to time. In the present case, the value of NT is determined by comparison with the experiments of Nelson et al. (1997). This gives NT ¼ 4:3; corresponding to a free stream turbulence level of approximately 0.5%. The calculations reported below use air at standard conditions, with kinematic viscosity of 1:5  10
5 m2 =s; and a free stream velocity of 20 m=s: All lengths will be in metres and velocities in m/s, and wherever convenient these units will be dropped. Given methods for calculating the flow and the transition position, the suction flow rates can be adjusted to achieve the desired result. The cost functions that were used in this study will be discussed below. Assuming, however, that a suitable cost function has been specified, the optimisation was performed using a GA (Goldberg, 1989; Holland, 1975) or SA (Metropolis, Rosenbluth, Teller, & Teller, 1953), both in standard

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

248

form. The GA has binary encoding, roulette wheel reproduction, single point crossover with a probability of 0.9, and mutation using a probability of 0.05. The SA uses Metropolis’ algorithm to update its current model, i.e. a new model mnew is generated by perturbing the current model (changing one of the suction flow rates). If fðmnew Þpfðmcur Þ then the perturbed model is accepted and replaces the current model mcur : If, however, fðmnew Þ > fðmcur Þ the new model is not always rejected, but is accepted with a probability

fðmnew Þ
fðmcur Þ Paccept ¼ exp
; ð17Þ T where T is the temperature at the current iteration. The temperature schedule used is T ¼ T0 ð0:99Þk ;

ð18Þ

where T0 is a constant and k is the iteration number. With (17) and (18), the probability of accepting a model with a higher cost decreases as the iterations continue.

5. Cost functions Consider first the design requirement for a nacelle, where the transition position must be located at a prespecified position while the suction effort is minimised. This can be formulated as the optimisation problem minimise fs ¼ K1

N X

v2i ;

ð19Þ

i¼1

subject to xT
d ¼ 0;

ð20Þ

where there are N suction panels, vi is the (uniform) suction velocity on panel i; xT is the transition position, d is the desired transition location, and K1 is a constant used to provide a convenient scaling for fs : Here, fs is used to estimate the effort made by the suction system. It is not, however, a true measure of the energy used by the pumps in the suction system (in the original experiments, the effort for each pump was Oðvpi Þ where p was between 2 and 3), but should vary monotonically with it, and therefore minimising fs will minimise the effort. It also has the advantage of providing a quadratic cost function. Note that the panels are weighted equally, which in effect assumes that all the panels have the same length. A weighting allowing unequal lengths could easily be included in (19). The problem defined by (19) and (20) forms a constrained optimisation/optimal control problem. In Tutty et al. (2000a, b), this problem was solved using a gradient based controller (optimiser) which worked well for simple cases, but did not work well in more complex cases, and in certain cases, could not converge. In order to use (19) and (20) with the GA and SA, it is necessary to define a cost function incorporating both (19) and

(20). One possible form is f ¼ fs þ K2 jxT
djn ;

ð21Þ

where K2 is a positive constant and n is a positive integer. K1 and K2 are not independent, as either could be scaled out of f without, in principle, affecting the results. However, K1 is retained to force fs to be Oð1Þ: Since the suction velocities will be Oð10
2 Þ or less, K1 ¼ 104 will be used throughout, and K2 will be used to adjust the balance between the suction cost and the constraint in the cost function. Values of one and two were used for n in (21). With n ¼ 1; the distance of transition from the desired value is penalised linearly. With n ¼ 2; the weighting on the penalty increases with distance, and close to the desired position there will be relatively little gain in approaching further, and, in fact, the total cost may increase if the suction cost increases. The idea that the transition position should lie within a region rather than a fixed position can be incorporated into the cost function by setting K2 to zero if jxT
djpd where d is a specified tolerance. Where a tolerance is used, d is taken as 0:025 m: As an alternative to the constrained optimisation problem defined by (19), (20), suppose now that the aim is to minimise the power consumption of the system as a whole, based on the design requirements for a wing or an aerofoil. A model has been developed (Rioual, Nelson, Hackenberg, & Tutty, 1996) which takes into account various factors giving a nonlinear cost in the form WT ¼ Wwd þ Wpump þ Wml
Wthr ;

ð22Þ

where WT is the total power consumption of the system, Wwd is the power consumed in overcoming the wake drag, Wpump is the power consumed by the suction system, Wml is the power due to the momentum loss from withdrawing fluid from the main stream, and Wthr is the power from the thrust gained by expelling this fluid back into the flow in an optimum manner. The wake drag term is given by r 3 Wwd ¼ UN yc ; ð23Þ Zp where Zp ðo1Þ is a coefficient reflecting the efficiency of the propulsive system, and yc the momentum thickness at the trailing edge of the plate, obtained from yc ¼ 0:037ðc
x0 ÞR
1=5 c
x0

ð24Þ

and x0 ¼ xT
36:9xT Rx
3=8 : T

ð25Þ

Rc
x0 and RxT are Reynolds numbers based on their respective lengths, and c is the chord (length) of the plate.

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

The momentum loss is m ’ 2 Wml ¼ UN ; Zp

249

0.007 0.006

ð26Þ

0.005

where m ’ is thePtotal mass P flow rate through the suction rvi li ; where li and mi are the system (m ’ ¼ mi ¼ length and mass flux for panel i). The power consumption of the pumping system is 0 1 !2 P X m 1 @ Zs ’ i Dpi 2 2A ; m UN
m Wpump ¼ ’ ’ i vi þ 2Zs Zp Zs r ð27Þ where Zs ðo1Þ is a coefficient reflecting the efficiency of the suction system, and Dpi is the pressure drop across panel i; given (Bieler & Priest, 1992) by 1 vi 2 t vi Dpi ¼ rC ð28Þ þ32mK 2 : 2 d P P In (28), C and K are empirical constants, d is the typical diameter of a hole in the panel, t is the panel thickness, P the panel porosity, and vi =P the velocity in the holes. Finally, the thrust gained by expelling the sucked fluid back into the flow is Z 2 Wthr ¼ m : ð29Þ ’ 2s UN Zp The values of the various parameters in (23)–(29) used in the calculation reported below are Zp ¼ Zs ¼ 0:9; t ¼ 0:001 m; d ¼ 5  10
5 m; C ¼ 1:6 and K ¼ 0:09; where the last two are taken from Bieler and Priest (1992). Note that the wake drag term, Wwd ; directly relates the power consumption and the predicted transition position, and for the cases considered below is the dominant term in the expression.

6. ResultsFtwo-panel-flat plate As a test case for the GA and SA methods, a twopanel-flat case will be used. In general terms, this case models the original experiments, as illustrated in Fig. 1 but without the lining used to generate a pressure gradient, and can be used for comparison with the gradient projection method (Tutty et al., 2000b). The suction panels are placed at 0:28pxp0:48 and 0:58pxp0:78: For this configuration, the natural transition position, i.e. that with no suction, is at xE0:93; downstream of the trailing edge of the second panel. This configuration has the advantage that the relationship between the suction flow rates and the transition position can be mapped completely, and hence the solution can be determined for any given case. Fig. 4 shows contours of constant transition position against suction rates. This figure has the same general form as that shown in Fig. 2 from the experiments. The

0.004

Panel 1

0.003 0.002 0.001 0

0.001

0.002 0.003 Panel 2

0.004

0 0.005

Fig. 4. Transition position against suction flow rates for the two panel flat plate configuration. The suction flow rates are in m/s and the contours are for constant xT from 1:0 to 1:5 m in steps of 0:1; and 1:55 m; going away from the origin. Also shown is the contour for the optimum value of the suction cost fs when xT ¼ 1:4: The optimum is given by the point closest to the origin where the lines are tangential.

contour for xT ¼ 1:55 has a different shape than for 1:5 or less, with no significant effect from panel 1. This occurs because, once sufficient suction has been applied on a panel, the boundary layer becomes very thin, and there is little to be gained by further suction. Once this level has been reached, which occurs for xT E1:55; panel 1 cannot play an effective role. In fact, for this particular case, with high flow rates on panel 2, transition can be delayed until xT E1:66 with suction on panel 2 alone. To perform a GA optimisation, it is necessary to specify the bounds and the discretisation for the suction flow rates. From Fig. 4, bounds could be chosen which restrict the search area to a small area around the optimum, e.g. for d ¼ 1:4; 0:001pv1 p0:002 and 0:002pv2 p0:003; which should lead to an efficient search. However, the results from this could be deceptive, since for more complicated cases, such information would not be available. Hence, bounds are taken which allow the full range of possible transition delay, by setting 0pvi p0:015: Also, as a first trail, a coarse 16  16 discretisation is taken, with a population size of 10 and a desired transition position of 1.4. The initial population, which is given in Table 1, was set to test the GA, and consists of the corner points of the grid, a point close to the 1.4 contour but not close to the point with the minimum cost, and the other points chosen randomly. Also given are the calculated transition position xT ; and the total and suction costs when d ¼ 1:4; K2 ¼ 10 and n ¼ 1 in (21). With these parameters, the optimum on the grid has ðv1 ; v2 Þ ¼ ð0:003; 0:002Þ with xT ¼ 1:402; f ¼ 0:155 and fs ¼ 0:130: This point is not one of the initial population. The results from this run are shown in Fig. 5. Fig. 5a shows the suction flow rates used for calculations at each iteration, Fig. 5b the corresponding transition position, Fig. 5c the total cost, and Fig. 5d a comparison between the suction and total costs. Also shown are the

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

250

Table 1 Initial population for the GA for two-panel-flat plate with the crude discretisation. Also given are the calculated transition position and the suction and total costs with K2 ¼ 10; n ¼ 1; and d ¼ 1:4 in (21) v2

xT

f

fs

0.0 0.0 0.015 0.015 0.001 0.014 0.007 0.002 0.009 0.012

0.0 0.015 0.0 0.015 0.003 0.001 0.005 0.010 0.013 0.010

0.93 1.66 1.30 1.65 1.45 1.41 1.56 1.62 1.64 1.62

4.0685 4.853 3.216 6.986 0.625 2.105 2.315 3.261 4.910 4.603

0.0 2.25 2.25 4.50 0.10 1.97 0.74 1.04 2.50 2.44

(a) 0.014 Suction Flow Rates (m/s)

v1

0.016

0.012 0.01 0.008 0.006 0.004 0.002 0 -0.002 0

5

10

15 Iterations

20

25

1.7

(b)

Transition Position (m)

1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 0

5

10

15

20

25

Iterations 7

(c)

6

Cost

5 4 3 2 1 0 0

5

10

15

20

25

Iterations 0.8

(d)

0.7 0.6

Cost

current best (least cost) model. There are not always 10 models shown in this figure as the GA does not always recalculate models which have already been tested, and because, as expected, there is increasing duplication of suction combinations as the iterations proceed. In fact, this run terminated after iteration 21 as it had converged on the best cost solution. The optimum solution was obtained on iteration 5, after a total of 38 function evaluations, including the 10 for the initial population (iteration 0), with a relatively large drop in total cost at this point. There was, however, an increase in the suction cost, but this was small compared to the decrease in cost due to moving closer to the desired transition position than in the previous best cost solution which had ðv1 ; v2 Þ ¼ ð0:001; 0:003Þ and xT ¼ 1:452: Runs with K2 ¼ 100 and 1 were performed to assess the effect of adjusting the balance between the two parts of the cost function. With K2 ¼ 100; the optimum solution is the same as with K2 ¼ 10; and again this was obtained on iteration 5, although 44 function evaluations were required, and the run had not terminated by iteration 39 when the specified number of function evaluations (200) was exceeded. However, with K2 ¼ 1 the balance between the terms has been changed sufficiently that the optimum moves to ðv1 ; v2 Þ ¼ ð0:002; 0:002Þ with xT ¼ 1:370; i.e. the weighting placed on the second term in (21) has been decreased to the extent that a comparatively large error in the transition position can be tolerated to accept a drop in the suction cost from 0.13 to 0.08. The optimum in this case was found on iteration 7, taking 43 function evaluations, and there was still a significant amount of diversity in the population by iteration 50. This set of three runs was repeated allowing the tolerance in the transition position. There was no significant difference in the results. As mentioned above, taking n ¼ 2 in (21) changes the balance between the terms in (21) as the desired transition position is approached. Consequently, the

0.5 0.4 0.3 0.2 0.1 0 0

5

10

15

20

25

Iterations

Fig. 5. GA for two panels with the coarse grid: (a) the sampled suction flow rates (þ ¼ v1 ;  ¼ v2 ), (b) the corresponding transition positions, (c) the total cost f ¼ fs þ 10jxT
dj and (d) the total cost and the suction cost fs : The lines show the best case attained so far. In (a) the solid line is for panel 1 and the dashed line for panel 2, and in (d) the solid line is total cost and the dashed line the suction cost.

optimum may also change. Without the tolerance in the transition position, the optimum with K2 ¼ 100 is still ðv1 ; v2 Þ ¼ ð0:003; 0:002Þ; but with K2 ¼ 10 it is ðv1 ; v2 Þ ¼ ð0:002; 0:002Þ; and with K2 ¼ 1 it is ðv1 ; v2 Þ ¼ ð0:001; 0:002Þ which has xT ¼ 1:316: Hence, using n ¼

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

2 rather than n ¼ 1 allows the transition to move further away from the desired position. The GA found the optimum in all three cases, but took longer than before, taking until iteration 14, 37 and 14, respectively. Again, allowing the tolerance made no difference to the results. A series of runs were performed using a fine 216  216 discretisation. With the pseudo-random number generator used in this study, this discretisation gives in effect a continuous distribution. Also, rather than calculating the transition position directly, the values from the 16  16 case were used with a polynomial interpolation procedure to obtain off-grid values. This was much quicker, and enabled a much more comprehensive investigation of the available options than would have been possible if the full calculation procedure has been used. However, a number of runs using the full procedure were performed and the results compared with those using the interpolation. There were no significant differences. Using this procedure the optimum with d ¼ 1:4 is ðv1 ; v2 Þ ¼ ð0:00125; 0:00246Þ with fs ¼ 0:0761: Above a population of 10 was used with a string length of 32 for the binary representation of the suction flow rates. Population sizes of 10, 20 and 40 were used to investigate the effect on the behaviour of the GA. The algorithm was set to terminate on the iteration in which the 600th function evaluation was performed. This implies a different number of iterations with each population size, but will give a better idea of the performance of the algorithm than fixing the number of iterations. To allow a direct comparison with the results using the coarse discretisation, the first 10 of the population were taken from Table 1 and the rest were generated randomly.

251

A full set of fine discretisation GA runs were performed using d ¼ 1:4 and all combinations of K2 ¼ 1; 10 and 100, n ¼ 1 and 2, populations of 10, 20 and 40, and with and without the tolerance on the transition position. Table 2 gives the final values from these runs for no tolerance, and Table 3 the values with the tolerance. Without the tolerance on xT ; the best final results were obtained with K2 ¼ 10; n ¼ 1; and a population of 20. This solution has ðv1 ; v2 Þ ¼ ð0:00127; 0:00245Þ with xT ¼ 1:400 and fs ¼ 0:0761 (Table 2), i.e. in effect the optimum. The results for this run are shown in Fig. 6. This figure is similar to Fig. 5 for the coarse discretisation, particularly in the way it jumps to the region of the optimum cost early in the run, but it also shows a gradual refinement of the solution as the run proceeds further. Changing the population size to 10 gave final values of ðv1 ; v2 Þ ¼ ð0:00142; 0:00238Þ with xT ¼ 1:400 and fs ¼ 0:0774; slightly away from the optimum but still good, particularly considering that these values were obtained on iteration 14. Using a population of 40 gave ðv1 ; v2 Þ ¼ ð0:00292; 0:00318Þ with xT ¼ 1:398 and fs ¼ 0:1017; which is much further away from the optimum. Holding the other parameters the same, but taking K2 ¼ 100; the algorithm had not approached convergence with any of the population sizes by the time the runs terminated, e.g. with 20 the final values were ðv1 ; v2 Þ ¼ ð0:00773; 0:00127Þ with xT ¼ 1:399 and fs ¼ 0:6131: With much longer runs, the optimum was approached, but clearly with the fine discretisation, too much weight is being placed on maintaining the transition position rather than minimising the suction cost to allow the algorithm to function efficiently.

Table 2 Final suction flow rates and costs from the GA for the two-panel-flat plate case with the fine discretisation, d ¼ 1:4 and no tolerance on the transition position K2

n

pop

v1

v2

xT

f

fs

1 1 1 10 10 10 100 100 100 1 1 1 10 10 10 100 100 100

1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2

10 20 40 10 20 40 10 20 40 10 20 40 10 20 40 10 20 40

0.00156 0.00137 0.00148 0.00142 0.00127 0.00292 0.00657 0.00773 0.00515 0.00104 0.00104 0.00087 0.00117 0.00123 0.00115 0.00144 0.00146 0.00141

0.00233 0.00238 0.00234 0.00238 0.00245 0.00318 0.00138 0.00127 0.00155 0.00169 0.00164 0.00157 0.00238 0.00234 0.00240 0.00236 0.00234 0.00237

1.400 1.397 1.397 1.400 1.400 1.398 1.400 1.399 1.399 1.272 1.273 1.247 1.383 1.382 1.384 1.398 1.396 1.398

0.0787 0.0788 0.0794 0.0774 0.0772 0.1242 0.4748 0.7423 0.3829 0.0540 0.0540 0.0558 0.0732 0.0732 0.0732 0.0769 0.0777 0.0767

0.0787 0.0755 0.0768 0.0770 0.0761 0.1017 0.4502 0.6131 0.2897 0.0375 0.0378 0.0324 0.0704 0.0699 0.0795 0.0766 0.0760 0.0763

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

252

Table 3 Final suction flow rates and costs from the GA for the two-panel-flat plate case with the fine discretisation, d ¼ 1:4 and with the tolerance on the transition position K2

n

pop

v1

v2

xT

f

fs

1 1 1 10 10 10 100 100 100 1 1 1 10 10 10 100 100 100

1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2

10 20 40 10 20 40 10 20 40 10 20 40 10 20 40 10 20 40

0.00141 0.00130 0.00166 0.00103 0.00096 0.00178 0.00381 0.00750 0.00399 0.00104 0.00104 0.00087 0.00112 0.00114 0.00109 0.00104 0.00121 0.00147

0.00223 0.00253 0.00229 0.00240 0.00245 0.00236 0.00159 0.00112 0.00308 0.00163 0.00164 0.00157 0.00235 0.00234 0.00240 0.00240 0.00231 0.00220

1.377 1.413 1.398 1.375 1.375 1.380 1.375 1.377 1.396 1.272 1.273 1.247 1.375 1.375 1.380 1.375 1.376 1.377

0.0694 0.0809 0.0798 0.0682 0.0690 0.0693 0.1709 0.5754 0.0965 0.0540 0.0540 0.0558 0.0677 0.0678 0.0694 0.0681 0.0680 0.0701

0.0694 0.0809 0.0798 0.0682 0.0690 0.0693 0.1709 0.5754 0.0965 0.0375 0.0378 0.0324 0.0677 0.0678 0.0694 0.0681 0.0680 0.0701

Taking K2 ¼ 1 produced much better results than with K2 ¼ 100: The best was with a population size of 20 which gave ðv1 ; v2 Þ ¼ ð0:00137; 0:00238Þ with xT ¼ 1:397 and fs ¼ 0:0755: This has a slightly lower suction cost than with K2 ¼ 10; but is further from the desired transition position, and has a slightly higher total cost (0.0788 against 0.0772). Taking n ¼ 2 allows in general a lower cost solution by choosing xT od and hence lower suction cost (Table 2). With K2 ¼ 1 there is a significant mismatch between the predicted and desired transition position, but with K2 ¼ 100 xT is forced back towards d: Given the discretisation and the shape of the solution surface (Fig. 4), not surprisingly allowing the tolerance on the transition position while maintaining the other parameters, produced lower cost solutions by moving the transition towards its lower limit of 1.375, producing essentially the same results as not allowing the tolerance but taking d ¼ 1:375 (Table 3). From these results, it appears that for this two-panel problem at least, when using the GA the best choice of parameters for the cost function (21) are K2 ¼ 10; n ¼ 1 and not allowing a tolerance on the transition position. Turning now to the SA, as well as the parameters in the cost function there are other factors which need to be considered. The value of T0 in the temperature schedule (18) must be specified. Note that this parameter is not independent of K1 and K2 : Instead of changing T0 ; it could be fixed and both K1 and K2 varied to produce exactly the same results as fixing K1 and K2 and varying T0 : However, for consistency with the GA, K1 is fixed at 104 and T0 is allowed to vary. Also, a discrete or continuous (within the limits of the pseudo-random

number generator) distribution could be used for the suction flow rates. If the value of T0 is too large, it will take too long for the algorithm to converge on a solution, while if T0 is too small then there may be insufficient variation and the algorithm may lock onto a local minimum which may be far from the best solution. A number of tests were performed and T0 ¼ 0:1 was chosen as an appropriate value. Consider first the coarse 16  16 discretisation. Fig. 7 shows results using the cost function which worked best with the GA with K2 ¼ 10; n ¼ 1 and no tolerance on xT : The optimum is found eventually, but this takes over 160 iterations, and it is only found once a higher cost solution is accepted at about iteration 120. The algorithm has difficulty moving away from the local optimum on this grid of ðv1 ; v2 Þ ¼ ð0:007; 0:001Þ which was found early in the run. If v1 is held at 0:007 or v2 at 0:001; then ð0:007; 0:001Þ is in fact the best solution. Hence a higher cost model must be accepted before the optimum can be found. Consider now the fine (continuous) discretisation, results for which are shown in Fig. 8. Clearly, although the desired transition position has been found, the solution after 600 iterations is far from optimal. Taking K2 ¼ 1 produced results that were no better, and with K2 ¼ 100; results that were much worse. A large number of other options were tried for the SA, including other values of T0 ; and taking n ¼ 2 in (21). The results were never better than those reported above, with either a worse solution or the algorithm taking much longer to approach the optimum. Convergence to a suboptimal solution was a common

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260 0.016

0.016 0.014

Suction Flow Rates (m/s)

Suction Flow Rates (m/s)

(a)

(a)

0.014 0.012 0.01 0.008 0.006 0.004 0.002

0.012 0.01 0.008 0.006 0.004 0.002

0

0

-0.002 0

5

10

15 20 Iterations

25

30

1.7

0

35

20

40

60

80

Transition Position (m)

1.5 1.4 1.3 1.2 1.1 1

140

160

180

200

(b)

1.5 1.4 1.3 1.2 1.1 1 0.9

0.9 0

5

10

15

20

25

30

0

35

20

40

60

80

7

(c)

6 5

Cost

4 3 2 1 0 0

5

10

15

20

100

120

140

160

180

200

Iterations

Iterations

Cost

100 120 Iterations

1.6

(b)

1.6

Transition Position (m)

253

25

30

35

Iterations

5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

(c)

0

20

40

60

80

100

120

140

160

180

200

Iterations

0.8

Fig. 7. SA for two panels with the coarse grid. (a) Suction flow rates. The solid lines are panel 1 and the dashed panel 2. The thin lines are the current model and the thick the best attained so far (these are the same for panel 2). (b) Transition position. (c) The total cost (solid=þ) f ¼ fs þ 10jxT
dj and the suction cost fs (dashed=). For (b) and (c), the symbols show the current model and the lines the best attained so far.

(d)

0.7

Cost

0.6 0.5 0.4 0.3 0.2 0.1 0 0

5

10

15

20

25

30

35

Iterations

Fig. 6. GA for two panels with the fine grid: (a) the sampled suction flow rates (þ ¼ v1 ;  ¼ v2 ), (b) the corresponding transition positions, (c) the total cost f ¼ fs þ 10jxT
dj and (d) the total cost and the suction cost fs : The lines show the best case attained so far. In (a) the solid line is for panel 1 and the dashed line for panel 2, and in (d) the solid line is total cost and the dashed line the suction cost.

occurrence. Hence it is concluded, that for this twopanel problem, the GA is more reliable and efficient than the SA, although the cost function must be specified so that there is an appropriate balance between the suction cost and the penalty on the mismatch in the transition position.

7. ResultsFsix-panel-flat plate Consider now a flat plate with six suction panels 0:2 m long, spaced 0:1 m apart, placed at 0:3ipxp0:3i þ 0:2; i ¼ 1; y; 6; so that the first panel starts at 0:3 m and the last finishes at 2:0 m: The natural transition position is still at 0:93 m; but with this configuration transition can be delayed to xE3:1 m: It was found in Tutty et al. (2000b) that for d ¼ 1:7; 2:0 and 2:6; a gradient projection (GP) based controller/optimiser converges to and maintains a solution in less than 30 iterations. However, with d ¼ 2:3 more iterations were required and the GP algorithm has trouble maintaining the solution due to the fact that within the tolerance of

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

254 0.016

(a) Suction Flow Rates (m/s)

0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 0

100

200

300 Iterations

400

500

600

Transition Position (m)

1.6

(b)

1.5 1.4 1.3 1.2 1.1 1 0.9 0

100

200

300

400

500

600

Cost

Iterations 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

(c)

0

100

200

300

400

500

600

Iterations

Fig. 8. SA for two panels with the fine grid. (a) Suction flow rates. The solid lines are panel 1 and the dashed panel 2. The thin lines are the current model and the thick the best attained so far (these are the same for panel two). (b) Transition position. (c) The total cost (solid=þ) f ¼ fs þ 10jxT
dj and the suction cost fs (dashed=). For (b) and (c), the symbols show the current model and the lines the best attained so far.

the algorithm the solution was not well-defined. Even when the algorithm converged to a solution in terms of the transition position and suction cost, there were relatively large changes in the individual suction flow rates, with a trade off of effort between different panels with no significant effect on the position of transition or the suction cost. For the GA an upper limit must be set on the suction flow rates. The same value as for the two-panel case above could be used but this would be wasteful in that much of the time the algorithm would be searching in regions that are of no interest. A simple suction distribution, which is of practical interest, is to set the same suction flow rate on all panels. With this distribution, it is found that although the transition

position can be moved to its maximum value with sufficiently large suction, with suction as low as vi ¼ 0:003; xT E ¼ 2:9; and 3:0 with vi ¼ 0:005: Accordingly, the suction flow rates were limited to 0pvi p0:005: The GA was run with 216 choices on each panel, and a population of 20. The first six members of the initial population were given by vi ¼ ðj
1Þ  0:001; j ¼ 1; y; 6; and the rest were generated randomly. Calculations were performed using (21) with n ¼ 1 and no tolerance on xT ; and K2 ¼ 0:1; 1; 10 and 100, with the runs terminated on the iteration the 1200th function evaluation was performed. The best result was with K2 ¼ 1; as shown in Fig. 9. The final values were obtained on iteration 31, after approximately 600 function evaluations. They give xT ¼ 2:30 and a suction cost of fs ¼ 0:197; which compares well with that from the ‘best’ case of 0.190 obtained with the GP algorithm (Tutty et al., 2000a, b). Note that due to the noise added to the suction flow rates in the GP algorithm, this difference is within the numerical uncertainty in the results. In this case, the preferred value of K2 is smaller than for the two-panel case. This arises from the change in the range of possible values of the two terms in (21) due to the change in the parameters, rather than a change in the balance between the terms. With maximum suction the two-panel case has fs ¼ 4:5 and K2 jxt
djE2:6 (K2 ¼ 10 and xT
dE0:26), a ratio of 0.58, whereas with the six-panel case fs ¼ 1:5 and K2 jxt
djE0:7 (K2 ¼ 1 and xT
d ¼ E0:7), a ratio of 0.47. Hence, the balance between the terms is roughly same in both cases. The results for the SA again using K2 ¼ 1; n ¼ 1 and no tolerance are shown in Fig. 10. The algorithm quickly converges on a solution but clearly this is much worse than that obtained by the GA (Fig. 9). The individual suction flow rates for the GA and SA at the end of the runs are given in Table 4. Also included are values from GP algorithm and the suction cost. Although there is little difference in the total cost for the GA and GP, there are relatively large differences between individual suction flow rates, consistent with the large variability found in previous work (Tutty et al., 2000a, b).

8. ResultsFtwo panels with pressure gradient A flat plate, which is commonly used as a test case in problems of this kind, causes minimum disturbance to the external flow field, and usually it is assumed that there is no significant change in the pressure in the fluid. However, for most bodies in an external flow, this does not hold, and the deviation of the flow around the body causes a significant change in the pressure. In the experiments, an external pressure gradient was generated by placing a lining on the wall of the tunnel to

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260 0.0026

0.002

(a)

0.0018 0.0016 0.0014 0.0012

0.004 0.0035 0.003

0.002 0.0015 0.001

0.0008

0.0005 0

0.0006 0

10

20

30 40 Iterations

50

60

Transition Position (m)

Transition Position (m)

(b)

3

0

70

3.5

2.5 2 1.5 1 10

20

30 40 Iterations

50

60

400

600 Iterations

800

1000

1200

(b)

0

0

200

2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8

0.5

200

400

70

600

800

1000

1200

Iterations 1.4

2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

(c)

1.2

(c)

1

Cost

Cost

(a)

0.0025

0.001

0.8 0.6 0.4 0.2 0 0

0

Cost

1 2 3 4 4 5

0.0045

Suction Flow Rates (m/s)

0.0022 Suction Flow Rates (m/s)

0.005

1 2 3 4 5 6

0.0024

255

10

20

30 40 Iterations

50

60

0.28 0.27 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0.19

200

400

70

600

800

1000

1200

Iterations

Fig. 10. SA for six panels with the fine grid: (a) suction flow rates, (b) transition position and (c) the total cost (solid=þ) and the suction cost fs (dashed=). For (b) and (c) the symbols show the current model and the lines the best so far.

(d)

This experimental situation with the lining is modelled by calculating the inviscid flow past a wall with a hump using potential theory, and using the velocity at the wall as ue0 ðxÞ in (8). The hump used is 0

10

20

30

40

50

60

70

Iterations

Fig. 9. GA for six panels with the fine grid: (a) the suction flow rates, (b) the transition positions, (c) the total cost f ¼ fs þ jxT
dj and (d) the total cost (solid) and the suction cost fs (dashed). The lines show the best case attained so far. In (a) only the best is shown for clarity.

yðxÞ ¼ hf ðzÞ with f ðzÞ ¼

(

ð30Þ

1
3z2 þ 2jzj3

if jzjp1;

0

if jzj > 1

ð31Þ

and first accelerate, then decelerate the flow. In general, accelerating a flow improves its hydrodynamic stability while decelerating it tends to destabilise the flow. Hence a flow with a forced pressure gradient can have quite different stability characteristics than one past a flat plate with zero pressure gradient.

z ¼ 2ðx
xh Þ=w0 ;

ð32Þ

where h ¼ 0:005; w0 ¼ 1 and xh ¼ 1: Suppose now that there are two panels at 0:56pxp0:92 and 1:12pxp1:48; so that there is one on the upstream and one on the downstream side of the hump. Contours of the transition position versus the suction flow rates

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

256

Table 4 Individual suction flow rates and suction cost for the different algorithms for the six-panel-flat plate case Panel

GA

SA

GP

1 2 3 4 5 6

0.00099 0.00155 0.00211 0.00168 0.00211 0.00214

0.00160 0.00044 0.00055 0.00382 0.00301 0.00220

0.00038 0.00180 0.00222 0.00296 0.00184 0.00174

fs

0.197

0.315

0.190

transition position, one at 1.7 and greater which is above the discontinuity in Fig. 11 and the other significantly less than 1.7 where the variation in the suction flow rates has placed some of the population below the discontinuity. A similar calculation was performed with the SA, and the results are shown in Fig. 13. This run finished with slightly worse values than the GA: ðv1 ; v2 Þ ¼ ð0:00178; 0:00314Þ with xT ¼ 1:701; f ¼ 0:0131 and fs ¼ 0:130:

9. ResultsFflat plate with the power balance 0.005 0.004 0.003 Panel 1 0.002 0.001

0

0.002

0.004 Panel 2

0 0.006

Fig. 11. Transition position against suction flow rates for the twopanel configuration with a pressure gradient. Contours are shown for xT ¼ 1:3–2.4 in steps of 0.1, and three contours for constant suction cost.

plus some of contours of constant cost are shown in Fig. 11. The transition position–suction flow rate surface no longer has the simple form found with the zero pressure gradient flat plate (Fig. 4). With (31) the contours of xT for xT E1:5–2.1 terminate/originate at the discontinuity in Fig. 11. Further, for much of, although not all, xT in the range of 1.5–2.0, the solution of the fundamental problem (19), (20) is at the discontinuity. Not surprisingly, the GP based algorithm used in previous work (Tutty et al., 2000a, b) failed when confronted with this problem, as would any solution scheme which relies on local smoothness of the data. The GA was run using a population of 20, the fine (216  216 ) discretisation, 0pvi p0:006; and n ¼ 1 and K2 ¼ 1 in (21) with no tolerance on the transition position. The results with d ¼ 1:7; where the optimum is at the discontinuity, are shown in Fig. 12. The initial population was random apart from the first member which had vi ¼ 0: The values at the end of the run were ðv1 ; v2 Þ ¼ ð0:00139; 0:00330Þ with xT ¼ 1:701; f ¼ 0:129 and fs ¼ 0:128: This is on the discontinuity and close to the optimum. From Fig. 12, it can be seen that again a good solution had been found early in the run. Also, the effects of the discontinuity can be seen in Fig. 12b, where there are two distinct bands in the calculated

Consider a flat plate 3:5 m long with 12 suction panels each 0:1 m long running from x ¼ 0:1 to 1:3 with no gaps between the panels, and 0pvi p0:015: The GA was run allowing 16 flow rates on each suction panel, i.e. discrete changes of 0.001 in the vi ; with a population of 20, and vi ¼ 0 as one member of the initial population with the others chosen randomly. The run was terminated once 1200 function evaluations had been completed. The results are shown in Fig. 14. The most expensive solution is that with no suction and transition at its natural position (xE0:93). In the power balance, the largest term by far is the wake drag which depends directly on the transition position. Almost any nontrivial distribution of suction will move transition a significant distance downstream, and therefore decrease the power significantly. As a result, there is a good solution in the initial population, with WT ¼ 35:5 compared with 52.3 with no suction. There is a gradual improvement throughout the run, which finishes with the suction distribution shown in Fig. 14c which has XT ¼ 2:29 and WT ¼ 33:6: This suction distribution may be a good one, but it is not the best possible one. In fact, a better distribution can be generated by dropping the suction on panel five from v5 ¼ 0:004 as shown in Fig. 14c to v5 ¼ 0:002 which has little effect on xT but drops the cost slightly to 33.5. However, if v5 is decreased further to 0.001, transition occurs at xT ¼ 1:07; and consequently the cost is high with WT ¼ 51:3: Examination of the details of the cost function shows that the optimisation procedure essentially seeks to push transition as far downstream as possible with the minimum suction effort, and this is reflected in the suction distribution shown in Fig. 14. Basically, there is a small amount of suction ahead of the natural transition position, which delays transition, and sufficient suction on the last few panels to almost completely suck away the boundary layer and therefore push transition well downstream. Although a realistic suction distribution has been produced, this configuration is not physically sensible as most of the panels are not used. Hence consider instead 12 panels starting at x ¼ 0:8; just ahead of the natural

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260 0.006

257

0.006

(a)

(a) 0.005

Suction Flow Rates (m/s)

Suction Flow Rates (m/s)

0.005 0.004 0.003 0.002 0.001

0.004

0.003

0.002

0.001

0

0

-0.001 0

5

10

15 20 Iterations

25

30

0

35

Transition Position (m)

Transition Position (m)

2.6

(b)

2.4 2.2 2 1.8 1.6 1.4 1.2 1 0

5

10

15

20

25

30

100

200

500

600

(b)

0

35

100

200

300

400

500

600

Iterations 0.8

1.4

(c)

0.7

(c)

1.2

0.6

Cost

1

Cost

400

2 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1

Iterations

0.8

0.5 0.4

0.6

0.3

0.4

0.2

0.2

0.1 0

0 0

5

10

15

20

25

30

35

Iterations

(d)

0.35

0

100

200

300

400

500

600

Iterations

Fig. 13. SA for two panels with a pressure gradient and the fine grid. (a) Suction flow rates. The solid lines are panel 1 and the dashed panel 2. The thin lines are the current model and the thick the best attained so far. (b) Transition position. (c) The total cost f ¼ fs þ jxT
dj (solid=þ) and the suction cost (dashed=). The symbols show the current model and the lines the best attained so far.

0.4 0.3 0.25

Cost

300 Iterations

0.2 0.15 0.1 0.05 0 0

5

10

15 20 Iterations

25

30

35

Fig. 12. GA for two panels with a pressure gradient and the fine grid: (a) the sampled suction flow rates (þ ¼ v1 ;  ¼ v2 ), (b) the corresponding transition positions, (c) the total cost f ¼ fs þ jxT
dj and (d) the total cost and the suction cost fs : The lines show the best case attained so far. In (a) the solid line is for panel 1 and the dashed line for panel 2, and in (d) the solid line is the total cost and the dashed line the suction cost.

transition position, and ending at x ¼ 2: A run with the same parameters as for Fig. 14 was performed, and the results are shown in Fig. 15. Again the initial population

contains a good solution with WT ¼ 23:0; which drops to 19.8 by the end of the run. In the final distribution shown in Fig. 15c, the suction applied on the first panel restabilises the boundary layer to an extent, with some suction on the central panels which delays transition and high suction on the last two panels which forces transition to x ¼ 3:13: Again, this can be interpreted as a strategy to force transition far downstream with the minimum suction effort. The SA was also run with these two configurations. The results for the case with the panels near the leading edge are shown in Fig. 16. This figure shows the familiar pattern of rapid change near the start of the run with a more gradual change afterwards. The final values, which were obtained on iteration 739, are slightly better than

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

258

3.5

2.2

Transition Position (m)

Transition Position (m)

2.4

(a)

2 1.8 1.6 1.4 1.2 1

(a)

3 2.5 2 1.5 1

0.8

0.5 0

10

20

30

40

50

60

70

0

10

20

30

54 52 50 48 46 44 42 40 38 36 34 32

50

60

70

55

(b)

(b)

50 45 40 35 30 25 20 15

0

10

20

30

40

50

60

70

0

10

20

30

Iterations

40

50

60

70

Iterations

0.016

0.016

0.014

(c)

0.012 0.01 0.008 0.006 0.004 0.002

(c)

0.014

Suction Flow Rate

Suction Flow Rate

40

Iterations

Cost

Cost

Iterations

0.012 0.01 0.008 0.006 0.004 0.002

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

x

Fig. 14. GA for 12 panels from x ¼ 0:1 to 1:3 with the power balance and 16 possible suction flow rates on each panel: (a) the transition position, (b) the cost and (c) the final suction distribution.

that found by the GA, with WT ¼ 33:5 and xT ¼ 2:29: Again, there is a small amount of suction ahead of the natural transition position and maximum suction on the last panel. Results for the SA with the panels further downstream can be seen in Fig. 17. The suction distribution has a pattern which is broadly similar to that for the GA (Fig. 15). The cost and transition position are also similar (WT ¼ 19:6 and xT ¼ 3:15). There is however a noticeable difference in the individual suction flow rates. In general, it appears that suction effort can be traded off between panels with little effect on the gross results.

10. Conclusions Simulations have been performed of boundary layer transition control using distributed suction with cost functions based on minimising the suction effort while holding transition in a fixed position, and minimising the total energy consumption of the system. GA and SA

0 0.8

1

1.2

1.4

1.6

1.8

2

x

Fig. 15. GA for 12 panels from x ¼ 0:8 to 2:0 with the power balance and 16 possible suction flow rates on each panel: (a) the transition position, (b) the cost and (c) the final suction distribution.

techniques were used to solve the non-linear optimisation problems. It was shown that both approaches are feasible, and can be used where gradient based methods fail. For the fixed transition/minimum suction case, it was necessary to define the constants in the cost function so that there was a suitable trade off between reducing the suction and penalising the error in the transition position. In general, the GA approach was better behaved in that it consistently found a good solution if enough iterations were performed. In comparison, the SA algorithm was prone to premature convergence, although this might be avoided by suitable adjustment of the temperature schedule so that there was sufficient variation in the early stages of the procedure. With the minimum energy criterion, physically realistic suction distributions were produced with a large number of suction panels. The results suggest that a reasonable general strategy is to place a judicious amount of suction ahead of and around the natural transition position to restabilise the boundary layer, with small amounts of suction downstream of this to further delay transition, and high suction towards the

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260 3.5

2.2

Transition Position (m)

Transition Position (m)

2.4

(a)

2 1.8 1.6 1.4 1.2 1

3

(a)

2.5 2 1.5 1 0.5

0.8 0

200

400

600

800

1000

0

1200

200

400

600

800

1000

1200

Iterations 55

54 52 50 48 46 44 42 40 38 36 34 32

(b)

50

(b)

45

Cost

Cost

Iterations

40 35 30 25 20 15

0

200

400

600

800

1000

0

1200

200

400

Iterations

600

800

1000

1200

Iterations 0.016

0.016 0.014

0.014

(c)

Suction Flow Rate

Suction Flow Rate

259

0.012 0.01 0.008 0.006 0.004

(c)

0.012 0.01 0.008 0.006 0.004 0.002

0.002 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

x

0 0.8

1

1.2

1.4

1.6

1.8

2

x

Fig. 16. SA for 12 panels from x ¼ 0:1 to 1:3 with the power balance and 16 possible suction flow rates on each panel: (a) the transition position, (b) the cost and (c) the final suction distribution.

Fig. 17. SA for 12 panels from x ¼ 0:8 to 2:0 with the power balance and 16 possible suction flow rates on each panel: (a) the transition position, (b) the cost and (c) the final suction distribution.

end of the panels to push transition as far downstream as possible, with the minimum suction effort. It is clear that an efficient suction system will require a good distribution of panels as well as a good distribution of suction on the panels. Work is in progress in this area on design optimisation methods which can choose the number and position of the panels as well as the suction distribution.

Carter, J. E., & Wornom, S. F. (1975). Solutions for incompressible separated boundary layers including viscous–inviscid interaction. Aerodynamic Analysis Requiring Advanced Computers, NASA SP347, 125. Goldberg, D. E. (1989). Genetic algorithms in search, optimization and machine learning. Reading, MA: Addison-Wesley. Holland, J. H. (1975). Adaption in natural and artificial systems. Ann Arbor: The University of Michigan Press. Mack, L. M. (1984). Boundary-layer linear stability theory. AGARDReport 709, 3.1. Mack, L. M. (1977). Transition prediction and linear stability theory. AGARD CP-224, 1.1. Malik, M. R. (1990). Stability theory for laminar flow control design. Progress in Aeronautics and Astronautics, 123, 3. Metropolis, N., Rosenbluth, A. W., Teller, M. N., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 1, 1087–1092. Nelson, P. A., Wright, M. C. M., & Rioual, J.-L. (1997). Automatic control of laminar boundary-layer transition. AIAA Journal, 35, 85–90. Rioual, J.-L., Nelson, P. A., Hackenberg, P., & Tutty, O. R. (1996). Optimum drag balance for boundary-layer suction. Journal of Aircraft, 33, 435–438. Tutty, O. R., Hackenberg, P., & Nelson, P. A. (2000a). Numerical optimisation of the suction distribution for laminar flow control. AIAA Journal, 38, 370–372.

Acknowledgements This work was supported by EPSRC through a ROPA grant.

References Bieler, H., & Priest, J. (1992). HLFC for commercial aircraftFfirst ELFIN test results. Proceedings of the first European forum on laminar flow technology (pp. 558–594). Hamburg, Germany.

260

W. MacCormack et al. / Control Engineering Practice 10 (2002) 243–260

Tutty, O. R., Hackenberg, P., & Nelson, P. A. (2000b). Gradient projection based control and optimisation of the suction distribution for laminar flow control. Proceedings of the Institute of Mechanical Engineering Part I: Systems and Control Engineering, 214, 347–359.

Veldman, A. E. P. (1981). New, quasi-simultaneous method to calculate interacting boundary layers. AIAA Journal, 19, 79. Wright, M. C., & Nelson, P. A. (1999). Four-channel suction distribution optimization experiments for laminar flow control. AIAA Journal, 38, 39–43.