483

Multi-dimensional aircraft surface pressure interpolation using radial basis functions T C S Rendall and C B Allen∗ Department of Aerospace Engineering, University of Bristol, Bristol, UK The manuscript was received on 23 July 2007 and was accepted after revision for publication on 22 January 2008. DOI: 10.1243/09544100JAERO263

Abstract: Multi-dimensional interpolation via radial basis functions is applied to the problem of using aircraft surface pressure data obtained both computationally and experimentally to obtain pressure distribution predictions through parameter space. In the most complicated cases, the data may be a function of spatial position, Mach number, Reynolds number, and angle of attack as well as other more intricate variables such as control surface deflections. Amalgamation of computational fluid dynamics and wind tunnel data for load prediction is currently a timeconsuming task, especially given the large number of load cases that need to be evaluated to achieve aircraft certification, so that an efficient tool for making rapid estimates based on all the information available would be of great use. The approach, using radial basis functions, is tested on a combination of simple computational and experimental results and found to offer great flexibility, while still being capable of reproducing relatively detailed features of the pressure distribution. Keywords: radial basis functions, multi-dimensional interpolation, datafusion, experimental aerodynamics data processing, compact functions, aerodynamic data interpolation, computational fluid dynamics, wind tunnel measurements

1

INTRODUCTION

During the design and development phase, aerodynamic calculations are required to cover the full-flight envelope of an aircraft. In the early part of the design process, less sophisticated methods such as vortex lattice, panel, or semi-empirical methods are used. As the design iteration process proceeds towards a fixed configuration, wind tunnel (WT) experiments, and computational fluid dynamics (CFD) simulations are used more extensively. Traditionally, WT experiments were the main data source, but it has now become common for these to be supplemented through the use of CFD simulations [1]. However, for the foreseeable future, the aerodynamic calculation process

∗ Corresponding

author:

Department

of

Aerospace

Engi-

neering, University of Bristol, Queens Building, University Walk, Bristol, BS8 1TR, UK. email: [email protected]; [email protected] JAERO263 © IMechE 2008

has to rely on multiple data sources. The main priority for aerodynamics engineers is to get high-quality data (regardless of the source) for the complete flight envelope at a given timescale and cost. The data analysis procedures for the two sources (WT and CFD) are markedly different because the CFD simulations provide a vast array of relatively dense data, whereas WT experiments provide smaller sets of sparsely arranged data. Also by the nature of experimental testing these points must be decided on before the experiment begins and may therefore not be placed at the most suitable locations for the final pressure distribution. Owing to the differences in these data, engineers have devised dedicated sets of tools appropriate for analysis of each data source. The authors believe significant gains can be achieved if the multiple data sources can be combined to produce one continuous data set, which can then be used as the final data set encompassing the full-flight envelope. The quantity of interest is normally the surface pressure and this depends on factors such as the Reynolds number, angle of attack, Mach number, Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

484

T C S Rendall and C B Allen

Fig. 1 Flow diagram showing the main components of a multi-dimensional and multi-source interpolation tool

and of course the three spatial dimensions. Hence, the aim of this work is to explore the applicability of a multi-dimensional interpolation tool for merging data arriving from many different sources (Fig. 1), together with its ability to interpolate this data onto any required point. Methods capable of solving such problems come with numerous parameters. In general, these parameters can be used to introduce intelligence into the process. For example, higher weighting can be attributed to CFD data for points that are known to be reliable, or higher weighting can be attributed to WT results for cases when CFD data are known to be inaccurate. Potentially, the method-specific parameters can also be used to locally change the interpolation characteristics when special flow characteristics are encountered (shocks, flow separation, etc). Examples of the use of artificial neural networks to interpolate multi-dimensional pressure data on an aircraft wing can be found in Hsu [2], Cao et al.[3], and Allen and Dibley [4]. In the present study, an interpolation technique based on radial basis functions (RBFs) is used. The useful properties of this method are its ability to operate in a space of any dimension, as well as giving smooth interpolants based on a number of variable tuning parameters. Section 2 is devoted to the theoretical description of RBF’s, whereas the RBF interpolation method has been evaluated using various test cases in sections 3, 4, and 5, and these are followed by concluding remarks in section 6.

coefficient, Cp , which depends on position as well as flow parameters such as Mach number, Reynolds number, and angle of attack. The method works by expressing the function in terms of the Euclidean distance to all of the points in the interpolation space where the function is specified, which is done via the basis functions. In order to recover the initial data, a linear system needs to be solved at the beginning of the process but after this, evaluation only consists of evaluating the basis function values from the distances. An evaluation may be carried out for any point in the parameter space. If f (x) is the function to be approximated (in this case Cp ), given as fi at N discrete points x i in a space of any number of dimensions, and s(x) an interpolation using the basis function φ, then the solution begins with the form of the required interpolation [5, 6], defined by s(x) =

2.1

INTERPOLATION SCHEME Radial basis functions

The goal here is to approximate a discrete function in any number of dimensions and from multiple sets of data, so that this interpolation may then be evaluated at other points in the parameter space to make predictions. In this work, the function is the surface pressure Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

βi φ(x − x i ) + p(x)

(1)

i=1

It is normal to append the polynomial p(x) to allow recovery of polynomials of that order and improve the quality of the interpolation. The polynomial is written as p = γ 0 + γx x + γ y y + γ z z + γ α α + γ M M

(2)

Higher-order terms could also be included at this point, as well as potentially non-linear ones. For example, it might be speculated that √a superior choice for the Mach function could be 1/ 1 − M 2 , based on the Prandtl–Glauert rule. However, for simplicity in this work linear terms have been used. The coefficients βi are found by requiring exact recovery (in the case where no smoothing is used) of the original function; this means that none of the original data are lost through the interpolation process. When the polynomial term is included, the system is completed by the additional requirement for all polynomials q(x) with degree less than or equal to that of p(x) (sometimes referred to as a ‘side condition’) i=N 

2

i=N 

βi q(x) = 0

(3)

i=1

This condition originates from a dot product ensuring that the polynomial space is orthogonal to the space of basis function coeffcients [5]. For easy computational treatment that maximizes the use of matrix multiplication, the problem is written in the following fashion. Vectors containing all of the data fi and interpolation coefficients βi and γ must JAERO263 © IMechE 2008

Multi-dimensional aircraft surface pressure interpolation

first be defined

radius, R the basis function argument is redefined as ⎞ γ0 ⎜ γx ⎟ ⎜ ⎟ ⎜ γy ⎟ ⎜ ⎟ ⎜ γz ⎟ ⎜ ⎟ ⎜ ⎟ a = ⎜ γα ⎟ ⎜ γM ⎟ ⎜ ⎟ ⎜ βs 1 ⎟ ⎜ ⎟ ⎜ .. ⎟ ⎝ . ⎠ βsN ⎛

⎛ ⎞ 0 ⎜ 0⎟ ⎜ ⎟ ⎜ 0⎟ ⎜ ⎟ ⎟ f =⎜ ⎜ 0⎟ ⎜ 0⎟ ⎜ ⎟ ⎝ 0⎠ f

⎞ f s1 ⎜ ⎟ f  = ⎝ ... ⎠ fsN ⎛

xscaled =

(4)

0 ⎜0 ⎜ ⎜0 ⎜ ⎜0 ⎜ ⎜0 ⎜ ⎜ C = ⎜0 ⎜ ⎜1 ⎜ ⎜ ⎜ .. ⎜. ⎜ ⎝ 1

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

x s 1 ys 1 zs 1 α s 1 Ms 1 .. .

.. .

.. .

.. .

.. .

1 1 xs2 xs1 y s1 ys 2 z s1 zs 2 α s1 αs 2 M s 1 Ms 2 φs 1 s 1 φ −S1 s1 s2 .. .. . .

··· ··· ··· ··· ··· ···

1 xs 1 s N ys 1 s N zs 1 s N αs1 sN Ms1 sN

· · · φs1 sN ..

.

x s N ys N zs N α s N M s N φs N s 1 φs N s 2 · · ·

.. . φ sN sN −SsN

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (5)

with φs1 s2 = φ(x s1 − x s2 )

(6)

where site s1 is the centre and s2 is the point for which the evaluation is being made. Linear terms are included in all the possible dimensions, with x, y, z, α, M given here as an example. The condition for the recovery of the data at the centres may now be written as f = Ca

(7)

Compact functions that decay to zero at a certain distance from their centre are used as this gives them a local character. The distance to this zero-boundary is termed the support radius and it is normally scaled to some desired value. A small support radius leads to a better condition number for the problem, but a less plausible interpolation between data, whereas a large support radius gives good interpolation but with a worse condition number. To include the support JAERO263 © IMechE 2008

x R

(8)

The basis function chosen for this work is Wendland’s C0, which is given for xscaled  1 as [5] φ(xscaled ) = (1 − xscaled )2

There is no distinction between CFD and WT data in f ; this only arises when it is necessary to select the correct smoothing parameter. The interpolation matrix, containing all basis function, polynomial and smoothing terms Ssi is also required ⎛

485

(9)

and for xscaled > 1 as φ = 0. Although all dimensions involved in an interpolation may have quite different scales, data are normalised to be in the range 0–1 in all directions. Hence, a certain support radius still has intuitive meaning in all directions, while all support radii quoted are given in these normalized units. For all calculations, the numerical value of the support radius has been kept at 1; this means that within any patch (section 2.2) all points have influence. This support radius and the C0 function are retained throughout for consistency but examples showing the influence of these parameters on interpolation may be found in Wendland [5] and Buhmann [6]. A smaller support radius localizes the interpolation and improves the condition number of the system, while the function changes the condition number but also the smoothness. Given that some smoothing will already be needed to accomplish the unification of the different data sets, a particularly smooth function is not required and the C0 function was therefore chosen to give the lowest condition number of the available compact functions of this type. The definition of the norm has a significant impact on the interpolation. Typically the Euclidean norm is used, so that

d=D

 x =  xd2 (10) d=1

where D is the number of dimensions. This implies, working in three dimensions, that a RBF will have a constant value on the surface of a sphere (or hypersphere in dimensions higher than three) positioned on a given centre. One possible modification is to include weighting coefficients in each of the dimensions

d=D

 kd xd2 (11) x =  d=1

Here, the coefficients may be regarded as a measure of how much a given axis direction should be considered insignificant for building the interpolation. This is because if the coefficient in front of any component is large then large distances in that dimension increase the norm more than the same distance in Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

486

T C S Rendall and C B Allen

other dimensions. Although the further away a point is the less influence it has on the interpolation (for a decaying basis function), this will bias the interpolation to concentrate on certain points. In effect the original hypersphere of influence has become a higher-dimensional ellipsoid, where the above coefficients are proportional to the inverse square of the size in that axis direction. Thus, a large kz corresponds to a higher-dimensional ellipsoid that does not extend far in the z-direction, thus ignoring points far away in z. Changing all coefficients simultaneously is equivalent to changing the support radius. Since there will always be significant variation betweenWT and CFD results on certain parts of the aircraft, typically in regions dominated by features such as shocks and largely separated flows, it is necessary to prevent overfitting (undesirable oscillations caused by fitting every point exactly) by applying smoothing. This is done, following Wahba [7] and Carr et al.[8], by subtracting from the diagonal entries of the basis function in matrix C and means that the interpolation minimizes a function that is a sum of the mean square error plus a function that decreases with improving smoothness of the interpolation, as demonstrated by Wahba [7]. Such smoothing may be applied differently to any subset of the data or even individual points, so that it is possible to smooth the CFD data more than the WT data, or vice versa. When this is combined with compact functions and norm weighting it is clear that a very flexible framework exists for defining interpolants.

2.2

Partition of unity

Partition of unity (PoU) methods [9, 10, 5, 11–13] are frequently used in many branches of engineering and science. These are a method of assembling many small local solutions in order to find the global solution to a much larger problem. The name ‘partition of unity’ stems from the fact that this is done using a weighting function that sums to unity throughout the domain. The global solution consists only of a weighted combination of the local solutions. Motivation for a PoU approach in RBF interpolation is straightforward. As the number of centres, N , increases the time taken to solve for the unknown coefficients grows as N α , where a value for α of around three would be typical (but depends on the exact method chosen). This is because the matrix contains N 2 elements and to perform the inversion it is necessary to operate on these order N times. Some speed savings might be made through the use of compact functions and a procedure designed for sparse matrices but a large interpolation might still not be attainable. Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

The PoU approach begins by separating the large set of centres, , into K patches k for k = 1, K . An interpolation is then found for each of these patches (an operation, which is naturally perfectly parallel) and these are blended together via the weighting function w. If the local interpolation is sk and the global interpolation is sg then a partition of unity is formed such that sg =

k=K 

wk (x)sk

(12)

k=1

Subject to the summing constraint k=K 

wk (x) = 1

(13)

k=1

The weighting function normally takes the form of a Shepard interpolant [14] ψk (x) wk (x) = n=K n=1 ψn (x)

(14)

where the function ψ may be defined as desired, providing it is smooth. A fair choice is to use any of Wendland’s compact functions [5]. The inherent problem in this approach is how to separate the centres into ‘sensible’ patches, while a further complication is that these patches must overlap to some extent so that the weighting function may vary smoothly across the overlap to prevent any discontinuities in the global interpolation. Here, the patches are found by choosing a small set of well-separated points and growing a patch outwards from each of these until they overlap suitably. A data tree (for example a bd or kd tree [5]) is built and used to find the nearest neighbours to form these patches efficiently with the ANN [15] package. Indeed, it is the construction of suitable patches that is likely to be the limiting factor for the largest interpolations, but if such calculations are desired then it would be possible to presort the data when it was first generated. A degree of regularity will always exist in CFD and WT results and this may be exploited to make the patching process more efficient, but at present maximum generality is maintained. In order to find a suitable set of points from which to grow the patches outwards a minimum separation distance algorithm is used. This simply selects a subset of the data, which are all separated from one another by at least some specified distance and gives an even distribution of points suitable for defining the patches. The point at the centre of each patch (centre point) is very useful as it fixes a location about which to define the patch weighting function. Similarly, the definition of the overlap is straightforward to define in terms of JAERO263 © IMechE 2008

Multi-dimensional aircraft surface pressure interpolation

487

the distance to the centre point. Although any point beyond the extent of the available data will have a sum of patch weights equal to zero, there is also a safeguard against interpolating unreasonably beyond the extent of the data. Once the patches are found the interpolation is built easily. A data file is written for each patch, which may be worked on independently of all other patches in a process that is perfectly parallel. Evaluating the interpolation is also speeded by the partition as a list of patch sizes and locations may be kept from when the interpolation was built, so that it is only necessary to refer to this list for each query point to determine which patches are needed. The appropriate patches for that query point are read, the basis functions calculated, weightings read, and the evaluation of that point completed. Typically no more than five patches will be needed per query, each containing a few hundred points, meaning that evaluation of each query point is fast. Again, evaluation is an easily parallelized process. 3

SPATIAL TEST CASES

Experimental and CFD data has been used in the form of chordwise pressure distributions at a range of angles of attack, Mach number, and spanwise positions for a given wing. The coordinate system is orientated such that x corresponds to the chord, y the span, and z the depth. A first test of the method is to explore the interpolations built using a reduced data set and to compare these results with the full original data. For this, a sectional pressure distribution is chosen, from which every second or fourth point is used to produce two interpolations. In this case only the CFD data are taken at the 10 per cent span position. These interpolations are then compared with the original data at the locations where points were removed. The results are shown in Fig. 2, indicating that even using every fourth point the pressure distribution is reproduced well. The second spatial test is to reproduce sectional pressure distributions for a full wing. Here, the CFD data are considered, but with every other spanwise location removed from the input data and defined as a query station so that a comparison may then be made to the original pressure distribution. In terms of a normalized spanwise position η, the sectional pressure distribution was known from η = 0.1 to η = 0.95 in increments of 0.05, allowing the query stations to be defined from η = 0.15 to η = 0.85 in increments of 0.1. Figure 3 gives the results of these interpolations and shows that generally the distribution is recovered well, with the poorest results occurring in the trailing edge region and in some cases the suction-peak near the leading edge. This may be of some concern since JAERO263 © IMechE 2008

Fig. 2

Interpolation quality using reduced sectional datasets

this area contributes significantly to the drag coefficient. Such an effect near the trailing edge may be due to adverse influence from points on the other side of the wing, since although on separate surfaces these points are spatially close to each other. As shown by Fig. 4, the problem may be alleviated by increasing the weighting coefficient in the chordwise direction and reducing it somewhat in the wing depth direction. It does seem sensible that since the greatest variations in Cp are to be found along the chord weighting the interpolation away from this direction is helpful; Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

488

T C S Rendall and C B Allen

Fig. 3

Spanwise sectional pressure distributions, kx/y/z = 1, 1, 1

if certain coordinate directions exhibit more gentle variations then a superior interpolation may be found by focusing on these. Interpolation was carried out with kx,y,z = 10, 1, 0.2.

4 4.1

NON-SPATIAL TEST CASES Angle of attack

Non-spatial test cases were defined for the angle of attack α in a similar way to the spatial test cases. Data were given for 12 angles between −7.0◦ and 16.1◦ at η = 0.399 and M = 0.70. The first case to consider is that of interpolation through a range of angles of Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

attack, from a negative value right up to beyond the stall, conducted here for only the experimental data. Every second angle of attack point was removed and used as a query point for a comparison to the original data. It must be noted that the number of chordwise points are significantly less than that available from CFD results. Figure 5 shows these results and indicates that the method successfully captures the presence of a transition bubble on the lower surface as well as separation at higher angles of attack, but that the region of poorest interpolation remains the trailing edge. Again, this may be due to the influence of points on the other side of the wing. In an attempt to rectify this the biasing from section 3 was included, but with more limited success, as shown in Fig. 6. JAERO263 © IMechE 2008

Multi-dimensional aircraft surface pressure interpolation

JAERO263 © IMechE 2008

Fig. 4

Spanwise sectional pressure distributions, kx/y/z = 10, 1, 0.2

Fig. 5

Pressure distributions with varying angle of attack – no bias

489

Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

490

T C S Rendall and C B Allen

Fig. 6

4.2

Pressure distributions with varying angle of attack – with bias

Mach number and spanwise position

Data to demonstrate a two parameter interpolation, through both Mach number and spanwise position, were produced by running an inviscid upwind code [16] at M = 0.70, 0.75, 0.80, 0.85 with a fixed angle of attack of 0.88◦ for the isolated wing of the multidiscipline optimization (MDO) aircraft [17]. Query points were defined in Mach number as M = 0.725, 0.775, 0.825 and additional CFD runs performed at these conditions to provide comparison data. Alternating spanwise stations were removed from the input data to allow some of these to be used as query points (η = 0.27, 0.50, 0.69, as shown in Fig. 7). Figures 8 and 9 show how the results may be influenced by changing the bias. An optimal combination uses a biasing of kx/y/z/M = 10, 1, 10, 1 and gives the best recovery of the leading edge suction peak and aft shock. It should be noted that this is a fairly harsh test case, as the original mesh used only 1617 surface points and with every other spanwise section removed the data are quite sparsely arranged. Also, although discrepancies exist in the pressure distribution, the sectional force coefficients are captured quite accurately. These bias coefficients were also used on the full data set (all spanwise stations included, as shown in Fig. 7) to demonstrate interpolation purely in Mach number. Figure 10 illustrates that the method does capture the variation in pressure distribution reasonably well. It also proved necessary to scale the Mach number dimension into the range 0–0.2 so that all Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

Fig. 7

MDO wing query stations and input data. Black dots denote input data points used when testing interpolation in both Mach number and spanwise position. Mesh indicates full set of data points, used when interpolating in only Mach number

patches included more than one Mach number, otherwise the points for some patches could lie on a plane in higher-dimensional space and the linear system became singular. JAERO263 © IMechE 2008

Multi-dimensional aircraft surface pressure interpolation

Fig. 8

MDO wing Cp , kx/y/z/M = 1, 1, 1, 1, C0, left to right η = 0.27, 0.50, 0.69. CFD data shown by heavy line, interpolation shown by light line

Fig. 9

MDO wing Cp , kx/y/z/M = 10, 1, 10, 1, C0, left to right η = 0.27, 0.50, 0.69. CFD data shown by heavy line, interpolation shown by light line

JAERO263 © IMechE 2008

491

Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

492

T C S Rendall and C B Allen

Fig. 11

Fig. 10

5

MDO wing Cp , kx/y/z/M = 10, 1, 10, 1, C0, interpolation left, CFD data right

MERGED DATA TEST CASE

Although interpolation based on data obtained exclusively from either a CFD or WT source is in itself useful, a major goal is to be able to combine both. Indeed, data might also be available from semi-empirical sources or even flight tests and this would also need to be incorporated in the long term. Using two data sets for Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

Spanwise loading for different smoothing parameters

an interpolation is therefore an important test case to consider. Sectional lift coefficient distributions, obtained through integration of the chordwise pressure distributions, are compared for results from a mixed data set in Fig. 11. Since the chord was already normalized along the span when the calculations were made, the load distribution deviates from the elliptical form that would normally be observed. The set used to build the interpolation contained every second spanwise CFD point and all experimental points, so that the unused CFD points could again be considered as query points. When mixed data sets are used the most important parameters are the smoothings applied to each of the sets of data as these determine whether or not the interpolation is overfitted or overly smoothed. An overfitted interpolation recovers every input point, meaning that it must zig-zag in an undesirable and meaningless fashion between all of the points, whereas an overly smoothed interpolation fits none of the points accurately and is equally meaningless. These interpolations were constructed for three different values of the smoothing parameters. Evidently, the greater the smoothing applied to any set the less accurately those points are recovered from the interpolation, whereas with no smoothing the result is overfitted because the CFD and WT data will always differ to some extent. It would be usual to lend greater weight to WT results and therefore, the smoothing applied to this data set should normally be small. An example showing the improvement made in a predicted pressure distribution through inclusion of experimental data is given in Fig. 12 (with no JAERO263 © IMechE 2008

Multi-dimensional aircraft surface pressure interpolation

Fig. 12

Effect of merging data sets on a single pressure distribution

smoothing). Originally, the CFD data did not show what appears to be a transition bubble on the lower surface, but this was present in the WT data at η = 0.662, and the influence of this has been carried over into the interpolation. The set used consisted of the CFD from η = 0.1 to η = 0.9 in increments of 0.1 and all wind tunnel data; this allowed a query station to be defined at η = 0.65. 6

CONCLUSIONS AND FUTURE WORK

Surface pressure interpolation in multiple dimensions, using RBFs, has been demonstrated on a variety of simple CFD and WT test cases both for interpolating along dimensions and for between data sets. The conclusion of the spatial tests is that RBFs do produce sound interpolations in these dimensions, but this is understandable owing to the huge amount of spatial information available. In comparison, information regarding either the angle of attack or Mach number is restricted and therefore it is more challenging to interpolate along these dimensions. The poorest results were obtained in the region of the trailing edge, with nonphysical influence from the opposite side of the aerofoil being the most likely cause. Ongoing work includes the following. 1. Consideration of restricting the interpolation to concentrate on physically relevant points; this could include removing influence from points on an opposing surface or allowing for an upwind form of the interpolation to prevent any downstream influence on an upstream supersonic pocket, and could be expected to improve the resolution of shocks. 2. Currently, the smoothing parameter is defined for each set, i.e. the parameter is only a function of that data set, but a more sensible option would be to use a smoothing parameter that is a function of a number of different variables. For example, under JAERO263 © IMechE 2008

493

conditions likely to produce separated flow (blunt edges or high angles of attack) CFD points should be smoothed much more than WT points, and this is currently being investigated. 3. The partition of unity also incurs a limit of the space dimension in which it may be used. If it is assumed that there is a maximum patch size suitable for rapid computational processing Nmp and that the data are evenly spaced with nd points in all dimensions up to a maximum of D then Nmp = ndD

(15)

An interpolation will not be possible, or worthwhile, if a patch contains fewer than two points in any direction, so that the maximum dimension of the interpolation space that may be used is Dlim =

ln(Nmp ) ln(2)

(16)

where Dlim is rounded down to the nearest integer. Even if it is possible to handle 1000 points at a time this limits the method to D = 9, while more practical values of Nmp = 600 and nd = 5 gives Dlim = 3. This restriction is problematic and solutions are under development. 4. Further questions arise over the optimal topology for constructing a partition of unity. It is clear that spheres pack the least tightly of all shapes and this means that the partition is not necessarily as efficient as it could be. A more suitable option would be to use blocks, but this is a more elaborate method and has not yet been implemented. 5. Ultimately, it may be preferable to skip one level of fidelity and go directly interpolation of global parameters, that is Cl = f (η, α, M , Re) Cd = f (η, α, M , Re) Cm = f (η, α, M , Re) (17) Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

494

T C S Rendall and C B Allen

If every section of the wing used 100 points for the Cp distribution then a saving of this order could be made in the number of points used, which is massive. Even more conveniently, the work would be taking place in a space of lower dimensions so that many more points from each dimension would be included in a given patch. Given that interpolation of order 104 points is easily possible with a partition of unity then a very large amount of information could be used. 6. All the points above require interpolations of large data sets and, hence, methods that reduce the number of centres required in the interpolation, for example using greedy methods [8, 18–21], are currently in development by the authors. 7. These reduction methods will be essential for the extension to volume data, which is the next step. This will be extremely valuable, for example, interpolation throughout an entire CFD solution so that the flow variables may be extracted for any spatial coordinate, or for interpolation of an entire volume solution onto a different mesh.

REFERENCES 1 Tinoco, E. N., Bogue, D. R., Kao, T.-J., Yu, N. J., Li, P., and Ball, D. N. Progress towards CFD for full flight envelope. Aeronaut. J., 2005, 109(1100), 451–460. 2 Hsu, K. C. The use of artificial neural networks in aircraft loads modelling. PhD Thesis, RMIT University, Australia, December 2005. 3 Cao, X., Sugiyama,Y., and Mitsui,Y. Application of neural networks to load identification. Comput. Struct., 1998, 69, 63–78. 4 Allen, M. J. and Dibley, R. P. Modeling aircraft wing loads from flight test data using neural networks. Technical report NASA/TM-2003-212032, NASA, 2003. 5 Wendland, H. Scattered data approximation, 1st edition, 2005 (Cambridge University Press, Cambridge, UK). 6 Buhmann, M. Radial basis functions, 1st edition, 2005 (Cambridge University Press, Cambridge, UK). 7 Wahba, G. Spline models for observational data, 1st edition, 1990 (Society for Industrial and Applied Mathematics, Philadelphia, USA). 8 Carr, J. C., Beatson, R. K., Cherrie, J. B., Mitchell, T. J., Fright, W. R., McCallum, B. C., and Evans, T. R. Reconstruction and representation of 3D objects with radial basis functions. In Proceedings of the ACM SIGGRAPH Conference on Computer Graphics, 2001, pp. 67–76. 9 Wendland, H. Fast evaluation of radial basis functions: methods based on partition of unity. In Approximation theory X: wavelets, splines, and applications, 2002, pp. 473–483 (Vanderbilt University Press, Nashville, Texas, USA). 10 Ahrem, R., Beckert, A., and Wendland, H. A meshless spatial coupling scheme for large-scal fluid-structure interation problems. Comput. Model. Eng. Sci., 2006, 12, 121–136. Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering

11 Schaback, R. and Wendland, H. Kernel techniques: from machine learning to meshless methods. Acta Numer., 2006, 12, 1–97. 12 Wendland, H. Computational aspects of radial basis function approximation. In Topics in multivariate approximation and interpolation (Ed. K. Jetter), 2005, pp. 1–26 (Elsevier B.V., Amsterdam). 13 Babuska, I. and Melenk, J. M. The partition of unity method. Int. J. Numer. Methods Eng., 1997, 40, 727–758. 14 Shepard, D. A two-dimensional interpolation function for irregularly-spaced data. In Proceedings of the 23rd ACM National Conference, Association for Computing Machinery, New York, 27–29 August 1968, pp. 517–524. 15 Mount, D. M. and Arya, S. ANN: a library for approximate nearest neighbor searching, May–June 2006, Version 1.1.1, release date 4 August 2006. Available from http://www.cs.umd.edu/ mount/ANN/. 16 Allen, C. B. Convergence of steady and unsteady formulations for inviscid hovering rotor solutions. Int. J. Numer. Methods Fluids, 2002, 41, 931–949. 17 Allwright, S. Multi-discipline optimisation in preliminary design of commercial transport aircraft. In Computational methods in applied sciences (Eds J.-A. Desideri, C. Hirsch, P. Le Tallec, E. Onate, M. Pandolfi, J. Periaux, and E. Stein), 1996, pp. 523–526 (ECCOMAS, Wiley). 18 Schaback, R. and Wendland, H. Adaptive greedy techniques for approximate solution of large RBF systems. Numer. Algorithms, 2000, 24, 239–254. 19 Ohtake, Y., Belyaev, A., and Seidel, H. Multi-scale and adaptive CS-RBFs for shape reconstruction from cloud of points. In MINGLE Workshop on Multiresolution in Geometric Modelling, Cambridge, UK, September 2003, pp. 337–348, available from citeseer.ist.psu.edu/ohtake03 multiscale.html. 20 De Marchi, S. On optimal locations for radial basis function interpolation: computational aspects. Rendiconti Del Seminario Matematico, 2003, 63(3), 343–357. 21 De Marchi, S., Schaback, R., and Wendland, H. Nearoptimal data-independent point locations for radial basis function interpolation. Adv. Comput. Math., 2005, 23, 317–330.

APPENDIX Notation a C Cp d D f f f N R s

vector of RBF coefficients interpolation matrix pressure coefficient index for dimension number of dimensions function to be interpolated vector of data, including leading zeros vector of data, excluding leading zeros number of data points support radius interpolation JAERO263 © IMechE 2008

Multi-dimensional aircraft surface pressure interpolation

Si w

smoothing at point i patch weighting function

β ψ

RBF coefficient smooth function used to define patch weighting

JAERO263 © IMechE 2008

495

Subscript g i k K

global centre index patch index number of patches

Proc. IMechE Vol. 222 Part G: J. Aerospace Engineering