Automatic modelling of airfoil data points

JID:AESCTE AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.1 (1-9)

Aerospace Science and Technology ••• (••••) •••–•••

1

Contents lists available at ScienceDirect

2

67 68

3

Aerospace Science and Technology

4 5

69 70 71

6

72

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13 14 15 16

Automatic modelling of airfoil data points

77

F. Pérez-Arribas ∗ , I. Castañeda-Sabadell

79

78 80 81

Universidad Politécnica de Madrid. (UPM), Avenida Arco de la Victoria 4, 28040 Madrid, Spain

82

17 18

83

a r t i c l e

i n f o

84

a b s t r a c t

19 20 21 22 23 24 25 26 27 28 29

85

Article history: Received 23 February 2016 Received in revised form 14 June 2016 Accepted 21 June 2016 Available online xxxx Keywords: Airfoils B-splines Parametric design Least squared fitting

30 31

This paper presents a computer-based technique to construct B-spline parametric models from a large set of airfoil data points, with a reduced number of parameters involved in the geometric representation of the airfoil profile. The proposed method uses different techniques related with the B-spline properties adapted to the geometry of an airfoil (a thin section with great changes of curvature) and produces a B-spline curve that is close to the data points maintaining a maximum tolerance distance. This curve can be used for calculations and is expected to provide a good framework for aerodynamic or hydrodynamic optimization, based on its reduced number of geometric parameters and on its calculation time, when compared with other methodologies. The method stresses the fitting of the airfoil’s leading edge, which has a significant impact on the properties of the airfoil. B-spline curves and surfaces are used in this method because they are widely used in CAD-CAM software products and can be easily exported to other programs. © 2016 Published by Elsevier Masson SAS.

86 87 88 89 90 91 92 93 94 95 96 97

32

98

33

99

34 35

100

1. Introduction

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

Geometric representation and parameterization are crucial when dealing with optimization methods. A wrong parameterization, with a large number of parameters, will produce unpractical solutions due to the computational time dedicated in optimization methodologies. As a result, robust and efficient parameterizations are required. In this paper, B-spline curves are used for the geometrical representation of airfoil shapes, because they are widely used in CAD-CAM software products and can be easily exported to other programs. But the main reason to work with a B-spline curves definition is that they can be used to reduce the number of design parameters as much as possible, while maintaining enough freedom and flexibility to represent a large number of airfoil data points due to their mathematical properties: considering how these curves are defined (see the Appendix for more details), notice that all the equation systems that are used in this paper are linear ones, Eqs. (2), (8), (9), (10), and this has the advantage that they no dot need an initial guess for the control points position. Rational curves such as NURBS need the use of numerical weights that multiply the control points. If these kind of curves are used, non-linear sets of equations have to be solved because the position of the control points is not known at the beginning.

60 61 62 63 64 65 66

*

Corresponding author. E-mail address: [email protected] (F. Pérez-Arribas).

http://dx.doi.org/10.1016/j.ast.2016.06.016 1270-9638/© 2016 Published by Elsevier Masson SAS.

Non-linear techniques would need initial guess and will require more computational time than linear ones. B-spline curves are also used by all CAD software products. These are the main reasons why B-spline curves are used. The objective of this paper is to produce a B-spline approximating curve, which fits a large number of airfoil data points, see Fig. 1(a), under a given tolerance with a reduced number of control points. As seen on this figure, the B-spline does not contain the data points, but is close enough below the tolerance crossing all the circles. In the case of Fig. 1, the radius of the circles or tolerance is 0.001 m for an MH43 airfoil of chord one. This curve will reproduce well the curvature field of the airfoil, when compared with interpolating methods, see Fig. 2, while maintaining the geometric characteristics of the original data points. Several authors demonstrated that curvature is a principal factor in determining the quality of a geometrical representation, [1,2] and this indicator has been used traditionally in CAD methodologies. As mentioned, this is a method that reduces the design parameters involved in the geometric representation of a 2D airfoil and simplify further optimization procedures, or even 3D surfaces constructed using the 2D curves to create wings, keels or foils. The final B-spline curve can be used as an initial guest for an aerodynamic or hydrodynamic optimization procedure, since it will adequately reproduce the original data points (hundreds) with about 7 to 12 control points as in Fig. 1(b). These control points can be used as design variables of the optimization process where the objective function could be the minimization of the difference be-

101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

2

AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.2 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

Fig. 3. Curvature field of the example, interpolating curves of the 3rd and 5th degrees.

20 21 22 23

Fig. 1. Example of the method: airfoil data points, B-spline fitting and leading edge details.

24 25 26 27 28 29 30 31 32 33 34 35 36

Fig. 2. Curvature field of the example at the leading edge, 5th degree.

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

tween the airfoil pressure distribution and a prescribed pressure distribution, or a different target, [3,4]. It is not the goal of this paper to create an optimised design from a performance point of view, but to enable the parametric definition of an airfoil with the use of a B-spline curve. Modelling the geometry of an airfoil is different from modelling other industrial objects; for example, an accurate representation of the shape of the airfoils is required, and this needs working with a large amount of information. A discrete set of data points, usually more than 100 points per airfoil is employed. If standard techniques are used, e.g. interpolation of the data points, very complex surfaces with a large number of control points are obtained and will likely present poor smoothness in the resulting surface computed from an interpolation of the data sets. In the case of an interpolating B-spline, the number of control points would be equal to the number of data points (see Eq. (2), considering N = np), so the same number of variables as the number of data points are needed to obtain the airfoil representation, and this means usually more than 100 variables. Another disadvantage is the aspect of the curvature graph that interpolating curves of large number of data points present: the curvature graph will oscillate, specially near the leading edge. As an example, Fig. 3 shows the curvature field of the same MH43 profile of Fig. 2, using interpolating B-spline curves of different degree. The curvature field of the proposed method is more uniform and faired. A constraint of the present problem is that the leading edge of the airfoil must be accurately reproduced by the curve because most of the aerodynamic characteristics are induced by this region.

86 87 88

According [17], pressures in the subsonic flow around an airfoil are governed by surface curvature, so in this case it will be important to reproduce well the curvature especially at the leading edge. This part of the surface presents high curvature values when compared with the rest of the shape of the blades. If for example a rudder is modelled with an inadequate definition of the leading edge of the sections that make up its body, additional drag to the one that has been predicted will be present in the final design. An iterative least squares (ILS) fitting process [5] with an original selection of the parameterisation based on the Haussdorf metric (the minimum of the maximum Euclidean distances) is described in Section 3 and also considers accurate modelling of the leading edge, see Fig. 1(c) and (d). The proposed method allows a reduction in the number of control points of the curve without reducing accuracy. The main reason for this reduction is the use of a non-uniform knot vector, as will be explained in Section 4 of the paper, with the knots located near the leading edge. This way, the control points are automatically located near this high curvature area, so this part of the airfoil will be well reproduced under a tolerance. Another key factor is the iterative parameterization as will be described in Section 3. Validation is produced in Section 7, where the presented method is compared with different ones available at the literature, showing a reduction in the number of control points for airfoils with equivalent tolerance. Data reduction in an object representation speeds up most of the downstream processes, improves the fairing of the resulting surface and decreases storage requirements in following design stages. Regarding related work about the modelling of airfoils, most of the references are based on interpolation of certain points amongst the large set of data points, trying to reduce the number of parameters. Reference [6] proposed a NURBS method that interpolates the data points instead of approximating them. They modified the data points to produce more points in the leading edge, whereas the presented method does not alter the data points. The same authors in [7] produced a new NURBS method based on interpolation resulting in 26 control points. Both references are related with aerodynamic optimization of wings. An interesting reference is [8], which focused on aerodynamic airfoil optimisation and used B-spline curves for its definition, obtaining curves for different cases with an error of about 0.002 and 14 control points, when applied to airfoils of chord one. The examples in that paper are compared with real airfoil shapes in

89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.3 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

Section 7. That work, based on optimization techniques that need an initial guess curve, began with a given definition of the control points of a B-spline that models the 2D airfoil data points and optimises the position of the control points until a given tolerance is reached. The presented method in this paper does not require an initial B-spline for modelling the airfoil points. Not all the references use B-spline curves. Reference [9] use a specific parameterization of curves called Class Shape Transformation (CST), that in a similar way with Bézier curves are based on the Bernstein Polynomials, but they do not depend on control points but in a set of 10 scalar parameters. This methodology produces errors of about 0.01 for the tested airfoils of chord one but can not be exported to other software unless a standard exchange export format (maybe NURBS) is used. An Appendix with the notation for this paper and a brief discussion of B-splines is at the end of this paper and will help to understand the mathematical properties of these curves. As mentioned, a technique based on the Least Squared (LS) fitting is used in this paper. In the case of curve adjustment, NURBS methods can use control points and/or weight factors (B-spline curves use unit weights), [10,11]. The present paper uses a progressive and iterative least squared (ILS) method that is efficient and intuitive for data fitting and requires low computational time (∼5 sec.). The present method also uses a specific technique to select the B-spline knot vector for an airfoil. The knot vector sets the end points of the different pieces or Bézier curves, which make up a B-spline curve. The knot vector consists in a set of parametric values usually considered between 0 and 1, see the Appendix for more details. If you use a CAD software, a uniform knot vector is used when you define a B-spline curve by indicating the curve’s degree and click the control points, and most of the fitting methods use this kind of uniform vector. The main originalities of this paper are the use of the nonuniform knot vector that depends on the shape of the airfoil, as described in Section 4, together with the use of the iterative parameterization described in the following section. The effect of knot vector on the shape of a B-spline curve has been studied for many years, [12,13]. This last reference used genetic algorithms. Most of the methods that optimize the knot vector are complicated and very time consuming due to the optimization process that they use. This paper presents an efficient method to set the knot vector for an airfoil, based on a specific tolerance and optimization techniques are not used. The method outline follows; First, the presented method reads the discrete data definition of the airfoil. Then the definition of a knot vector that adapts well to the airfoil shape is defined, producing a non-uniform knot vector (Section 4). After this, a discrete 2D definition of the airfoils is defined, and an iterative curve fitting scheme is applied to the discrete data, resulting in a B-spline curve with its control points (Section 6). This fitting considers a fairing parameter in order to obtain realistic curves, since the problem is strongly non-linear and can produce unfaired solutions. Finally, the method is tested considering different airfoils shapes from the literature in Section 7, and compared with some validation results. The Appendix about B-splines completes this work.

57 58

2. Discrete Airfoil definition

59 60 61 62 63 64 65 66

Airfoils, at their most basic conceptual level, can be defined as a thickness distribution, often with a camber deflection. For many airfoils the geometry is mathematically derived, and the airfoil data is put into an ordinate file. When the user goes to design a wing’s airfoil in a surface modelling package, the starting point is that ordinate file and then, the creation of the best possible representation of the airfoil so that the wing will fly in a predictable

3

67 68 69 70 71 72 73 74 75 76 77

Fig. 4. Obtaining the parameters t j for the data points fitting.

78

and expected manner. When starting an optimization process, the coordinates contained in this file will produce an initial guess for the mathematical techniques that are used. So, the most usual way is to use a file that starts in the airfoils trailing edge, and then goes in order through the upper part of the airfoil to the lower part, crossing the leading edge near the origin of the file. These files of data points can be obtained from software or from airfoil data bases and they can contain a higher density of points in the leading edge, in order to accurately reproduce this important part of the airfoil. Notice that these data are subject to measurement errors and noise. The airfoil data points are defined as Qj = Q(x j , y j ), j = 0, . . . , np where np usually is over 100 points. The origin O of the airfoils can be placed well within the chord midpoint, at the point of maximum thickness or usually at the leading edge.

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

3. Fitting data points with B-splines using the least squares technique

98 99 100

Once the np + 1 discrete data points Qj have been obtained, the B-spline fitting can be conducted. A large number of data points do not suggest the use of an interpolating B-spline. Instead, an approximating curve c(t) is required. B-spline are parametric curves and if they are normalized, the interval range will be: t ∈ [0, 1]. See the Appendix for more details. This curve does not pass through the data points exactly but close enough to the points to capture their inherent shape. This is the well-known least squares (LS) approximation, which is discussed in [5]. In a general fitting problem, np + 1 data points Q0 , . . . , Qnp are approximated by a B-spline of the pth degree with N + 1 control points P0 , . . . , PN , N < np that are unknown and are obtained as the final results of the calculations, see Fig. 4. The general LS problem is described by an over-determined set of np + 1 equations with N + 1 unknown variables:

101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119

p

p

120

p

B 0 (t 0 ) · P 0 + B 1 (t 0 ) · P 1 + · · · + B N (t 0 ) · P N = Q 0 p B 0 (t 1 ) ·

p B 1 (t 1 ) ·

p B N (t 1 ) ·

P0 + P1 + · · · + PN = Q1 .. .. .. .. . . . . p p p B 0 (tnp ) · P 0 + B 1 (tnp ) · P 1 + · · · + B N (tnp ) · P N = Q np p

121 122

(1)

B i corresponds to the ith basis function of a pth degree B-spline that is calculated using de Boor’s algorithm, Eq. (12) in the Appendix, and considering a uniform knot vector. In Eq. (1), t j ( j = 0, n) represents the parameters associated with the data points. The graphical meaning of this t j is depicted in Fig. 4

123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.4 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

4

Matrix forms are convenient to solve the problem:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

⎤  p p p B 0 (t 0 ) B 1 (t 0 ) ... B N (t 0 )  P0 ⎥  P ⎢ p p p ⎢ B 0 (t 1 ) B 1 (t 1 ) ... B N (t 1 ) ⎥  1 ⎥· . ⎢ .. .. .. .. ⎥  . ⎢ ⎦  . ⎣ . . . .  PN p p p B 0 (tnp ) B 1 (tnp ) ... B N (tnp ) ⎡

 ⎡ ⎤  Q0   ⎢ Q1 ⎥  ⎢ ⎥  = ⎢ .. ⎥  ⎣ . ⎦   Q np

(2)

[M ] · [ P ] = [ Q ] ⇒ [M ] · [M ] · [ P ] = [M ] · [ Q ] T

T

This system of (np + 1) by (N + 1) equations is solved by multiplying both sides of Eq. (2) by [ M ] T , which creates a determined (N + 1) by (N + 1) linear system. This type of system can be poorly conditioned, especially if a large number of control points are used. A conventional technique should not be used to solve this ill-conditioned system. Instead, a single-value decomposition of [ M ] T · [ M ] and a later back-substitution process are performed. The solutions of this system are the control points of the best B-spline fitting. Up to now, a standard LS fitting has been described. Approaching this problem with a standard parameterisation such as centripetal or chord-length parameterisation is correct but does not consider the effect of the distance of the data points to the Bspline. In this method, parameterisation based on a minimum distance metric is adopted (Haussdorf metric). The process is iterative and is described by the following three steps:

26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

1. The method begins with a centripetal parameterisation of the Qi points, and system of Eq. (2) is solved. This produces a starting curve of the iterative process only for the first loop, so no initial guess is needed. 2. For each Qj , the minimum distance (Euclidean) to the B-spline is calculated, see Fig. 4, which leads to a solution to Eq. (3).







Q j − c i (t j ) · c i (t j ) = 0

(3)

Once the solution has been found, this t i value is the parameter associated with the point Qi when solving system of Eq. (2). 3. After obtaining the t j ( j = 1, . . . , np) values, the distance d j = (Qj − c(t j )) is computed, which is the Euclidean distance between Qi and the B-spline. This distance is used to verify the shape requirements. If the maximum distance d j ( j = 1, n) is above a given tolerance, steps 2 and 3 are repeated until an acceptable maximum distance is obtained. The quality of the obtained curve is measured using the tolerance constraint, and the shape of the B-spline is amended using parameterisation of Eq. (3).

48 49 50 51 52 53 54

Hitherto, the iterative least squares (ILS) of a set of points has been described. This technique is used first to determine the best knot vector for the B-spline that fits the airfoil data points, and once this knot vector has been obtained, a new stage of the ILS is used to obtain the final B-spline. These two steps are described in the following sections.

55 56

4. Obtaining the degree and the knot vector

57 58 59 60 61 62 63 64 65 66

A generic B-spline of the pth degree and with N P control points, is composed of N P − p Bézier curves of the same pth degree, with p + 1 control points each. The joining points of the Bézier pieces are called knots. The mathematical continuity C between two Bézier curves at their joining knot depends on the degree of the curves, and is defined as C = p − 1. This means that if p = 2 ⇒ C = 1 and only the first derivative is continuous between enjoyment Bézier pieces. The curvature is a function of the curve’s first and second derivatives, so the minimum degree that produces

a continuous curvature is p = 3 when C = 2 between pieces. Increasing the B-spline degree increases the number of derivatives that are continuous at its knots, and therefore the continuity of the curvature at these points. The knot vector also plays a role in this continuity. See [16] for more details on Bézier and B-spline curves. For example, the curvature graph for a p = 3 curve as depicted in Fig. 3, presents jumps at the knots while the graph is more continuous for p = 5, the lower curve in that figure. So, the lowest degree that is normally used in interpolation or approximation techniques is 3. Increasing the degree does not guarantee good results in the curvature as shown in that figure, specially if a large number of control points is used. That is where the selection of the knots produces a major effect. When the designer is using a CAD software to produce B-spline curves and surfaces, apart from selecting the position of the control points, the only parameter that can be set is the curve’s degree. The knots, that are mathematically computed in the called knot vector, are taken uniformly distributed. That means that if a B-spline is formed of two pieces, independently from the degree, the knot value in the parametric space where the B-spline c(t) is computed is placed at t = 0.5. If the B-spline is composed of 3 pieces its knots will be placed at t = 1/3, and t = 2/3, and so on. The reason for this is that a uniform knot vector, see the Appendix for more details, produces good results in most of the engineering applications, since the user can add more control points where needed to increase the detail of specific areas of the curve. In the case of interpolating B-splines in a CAD software, the user clicks points on the curve (knots) and in this case the software may ask the user to use uniform knots, or compute knots based on distance between points or in the Sqrt[distance between points]. In the proposed method, knot vector is selected in a nonuniform as described in this section. The ILS technique of Section 3 could be used, increasing the number of control points N + 1 until obtaining a B-spline under a desired tolerance, with the limit N + 1 = np + 1 when an interpolating curve would be obtained. As mentioned in the introduction, a large number of control points will produce complex curves and surfaces constructed based on these curves, and the large number of parameters (control points) makes optimization procedures more difficult and longer. The present method produces fitting curves under a given tolerance by playing with the knot vector, that is selected in a nonuniform way following this steps:

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112

1. For a given degree p, the data points are first fitted with independent Bézier curves joined between them in C = 0 continuity (the end control point of one curve is the first one of the following one) that are obtained with the ILS technique of Sec. 3. The use of C0 is the most simple one. Authors also test the use of Bézier pieces in C1 (First derivatives: tangent) and even C2 (second derivatives: curvature) in a generic fitting problem, but in the case of airfoils with a high curved leading edge, the used of C1 and C2 worked very bad, producing unrealistic solutions without an airfoil shape. C1 and C2 are obtained considering the position of the Bézier control points, and these conditions can be

113 114 115 116 117 118 119 120 121 122 123 124 125

So we decided to use C0 pieces. The use of C0 is used to detect when the curvature is increased, that is near the leading edge, and when it is reduced, after passing the leading edge. This way knots will be placed near this area, and as a consequence also the control points. The use of periodic knots will produce a closed B-spline, and the trailing edge, sharp, will not be well reproduced.

126 127 128 129 130 131 132

JID:AESCTE AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.5 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

5

1

67

2

68

3

69

Fig. 5. Example of knots distribution, p = 3.

4

70

5

71

6

72

Fig. 7. Example of knots distribution, p = 5.

7

73

8

74

9

75

10

76

11

77

12

78

Fig. 8. Example of knots distribution, p = 4.

13

79

14 15

Fig. 6. Detail of knot distribution at the leading edge, p = 3.

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

2. Starting in the trailing edge Q0 , data points are introduced sequentially one by one maintaining the number of control points of the Bézier curves equal to p + 1, first solving the system of Eq. (2) and then obtaining the maximum distances according Eq. (3). 3. When a new data point Qk is introduced, the maximum distance is computed and checked with the tolerance value. If the maximum distance is bellow the tolerance, the following point Qk + 1 is introduced, and equations (2) and Eq. (3) are again computed considering the new extra data point. If the maximum distance is over the tolerance, that means that a new Bézier piece should start, as described in Fig. 5. In this figure, the large black dots mark the star or the end of the mentioned Bézier pieces. The figure shows the different p = 3 Bézier pieces for a AG03 airfoil of chord one with 100 data points. Every piece will have p + 1 = 4 control points. In this case the airfoil can be modelled in 4 pieces with a tolerance tol = 0.001 m. Notice how for the last piece, the 4 control points are mostly aligned, since this airfoil has a flat aft part. It seems clear how the leading edge will need more pieces than the rest of the airfoil parts. 4. The following Bézier piece will have the first control point equal to the last one of the previous piece, and the data points are now introduced starting from Qk + 1 and finishing when the maximum distance of the new subset of data points computed with Eq. (3) overpasses the tolerance. Fig. 6 shows a detail of the leading edge of the previous AG03. This part of the airfoil will need 2 pieces to fit the points with the tolerance, the second piece fits from points 42 to 50, the third piece from 50 to 56. In this figure, the radius of the circles has been increased to 5 mm, so the numbers inside are clearer. Notice how the Bézier curves do not present a continuity between them, as is clear at point 50 and 55. 5. Once the B different Bézier pieces have been obtained with the previous steps, the knots ki will be obtained as sum of lengths, divided by the total length:

i

1 Lj

55

ki =  B

56

1

57 58 59 60 61 62 63 64 65 66

Lj

i = 1, B − 1

(4)

Where L j is the length of the Bezier piece of index j, and B is the number of pieces that models the airfoil data points. In the case of Fig. 5 B = 4, and the 3 knots will have the parametric values of 0.433, 0.516, and 0.556. If a uniform knot vector is employed, the parametric values would have been 0.25, 0.5 and 0.75. Increasing the degree, produces more control points in the curves and this has the effect of diminishing the number of pieces that are necessary to model the data points. Fig. 7 shows the knots for a DFVLR R-4 unit chord airfoil with 200

data points, fitted with p = 5, tol = 0.001 m. Notice how the leading edge is modelled with a single piece, and in this case the trailing edge concave part will need another one. In this case p = 5 means 6 control points per curve, and B = 3 with the knots at 0.484 and 0.929. A uniform knot distribution would place the knots at 0.333 and 0.666. It can be seen how most of airfoils can be modelled with a low number of Bézier pieces, and that the leading edge will need one or two of these pieces to be reproduced accurately. Fig. 8 shows a MH 121 unit chord airfoil of 100 data points, fitted with p = 4, that means 5 control points per curve. The leading edge will need 2 pieces, with the knots parameters 0.320, 0.505 and 0.793. Although the number of pieces B can be quite similar for all the airfoils, the knots values depends on the airfoil geometric characteristics.

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

5. Subdivision of the knot vector close to the leading edge

96

Up to now, the position of the knots has been calculated in accordance with the previous sections. This determines the mathematical knot vector (see Appendix) that is necessary to compute a B-spline according to Eq. (5).

98

[0, . . . , 0, ki =1, B −1 , 1, . . . , 1]    

103

p +1

97

(5)

p +1

Because of the properties of B-spline curves, the size of the knot vector (knots), the degree p and the number of control points N P are related as knots = N P + p + 1. So, the number of control points is calculated directly once the knot position has been obtained following Section 4. For the airfoil of Fig. 5, N P = 11 − 4 = 7 control points, and in the case of the example of Fig. 7, N P = 14 − 6 = 8 control points. Nevertheless, this number of control points and knot vector do not guarantee that the tolerance is fulfilled when an ILS technique of the np + 1 airfoil data points is computed using Section 3, as depicted in Fig. 9 where the fitting of the example of Fig. 5 has been computed with the previously mentioned 7 control points. The approximation of the leading edge is quite poor, and the maximum distances to the data points are much larger than the tolerance (red circles of 0.001 m). This is solved in the present method by subdividing some Bézier pieces in a desired number P of equal parts, starting from the second piece, and finishing in the next to the last one. As mentioned and depicted in Fig. 5, Fig. 7 and Fig. 8, the leading edge is covered starting from the second piece until the next to last one. That means that extra knots are added in the knot vector close to the leading edge, and this way new control points will appear near this area and so, a better fitting will be obtained there. This effect is showed in Fig. 10 where the same example of Fig. 5 is fitted subdividing the 2nd and 3rd pieces in two. The knots are now 0.433, 0.474, 0.516, 0.536 and 0.556 in the parametric space, and the ILS technique of Section 3 can be computed with 9 control

99 100 101 102 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

6

AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.6 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

Fig. 9. B-spline fitting, no subdivision.

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

Fig. 10. B-spline fitting with subdivision.

31

97

32 33 34 35 36 37 38 39

Fig. 11. Fitting without minimizing derivatives.

40 41 42 43 44 45 46 47

points instead of 7. These two new control points together with the new knots placed near the leading edge, have a great effect over the fore area that now matches the tolerance of 0.001 m as shown in Fig. 10, in comparison with Fig. 9. Subdivision in 2, 3 and up to 4 parts according to the mentioned scheme has produced good results in the fitting of all the airfoils that authors have tested.

48 49

6. Final adjustment of the airfoil data points

The improvement of the ILS presented in this paper is the addition of some kind of fairing of the control polygon. The problem is that the ILS only minimizes the distance between the points and the B-spline, and does not consider the shape of the resulting curve. This means that the result of Eq. (2) will be a curve that may squirm or curl producing big disturbances in the curvature, yet it is as close to the data points as a result of the ILS technique. Fairing curves is a very subjective matter and there are plenty of criteria. In this paper a global fairing criteria is applied by using an integral formula, following Ref. [14]. The fairing will try to minimize function

54 55 56 57 58 59 60 61 62 63 64 65 66

100 101 102 103 104 105 106 107 108 109 111

 1  2 d  du f =  c ( u )  du 2 

112

(6)

0

113 114 115 116

51 53

99

110

50 52

98

Some final improvements for the ILS method of Section 3 have been developed in this work. The reason is that the system of equations (1) together with the way that the parameters of the data points are obtained by iteration, can be unstable and produce strange results, specially if the number of data points is large and with high order curves. If one wants to place the control points of a B-spline manually to resemble the airfoil data points, the control point will mimic the shape of the airfoil, and this can not be the case when using the LS technique for the afore mentioned reasons. This is depicted in Fig. 11 where the fitting is made for an Eppler 335 airfoil of chord one, with a B-spline of degree p = 5, and subdividing the pieces in P = 3 parts as explained in Section 5. This produces a 10 control points B-spline that can not fit the data points with a good tolerance, unless the control polygon is deformed as shown.

The problem arises because there is not a direct definition of the B-spline until the ILS method is applied and it is not possible to include Eq. (6) directly into the definition of the ILS fitting. The key is to change the focus from the curve itself to its control polygon, since it is inside the definition of the system of Eq. (2) and is the final result of the ILS fitting. As a result of the properties of the B-Splines, the curve mimics the shape of its control polygon: if the polygon conducts well (does not wiggle or twist), so does the B-Spline. So, instead of applying Eq. (6) directly, the method considers the minimization of the second derivatives included in Eq. (6) changing derivatives by differences of the control points. If the differences of the curve’s control points decrease, so do the derivatives of the B-spline. This can be accomplished setting the conditions of Eq. (7) to the control points for the second differences:

117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.7 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

1 2

P 0 − 2P 1 + P 2 = 0 P 1 − 2P 2 + P 3 = 0

4 5 6 7 8 9 10 11

⎡

1 −2 ⎢0 1 ⎢

1 ··· −2 · · ·

.. .

.. .

13

0

0

0

14

19 20 21 22 23 24 25 26 27 28 29

70 71

0 0 0 0

⎤ ⎡

0 0

⎤

0 P0 ⎢ P1 ⎥ 0⎥ ⎥ ⎢ ⎥

⎡

p B 0 (t 0 ) p B 0 (t 1 )

⎢ ⎢ ⎢ .. ⎢ . ⎢ p ⎢ B 0 (tnp ) ⎢ ⎢ 1 ⎢ 0 ⎢ ⎣ ..

75

0 ⎢0⎥ ⎢ ⎥

76 77

PN

78 79

0

80

(8)

34 36 38 39

−2

.

0

0

···

⎡

0

⎤

··· ··· .. . ···

p B N −3 (t 0 ) p B N −3 (t 1 )

p B N −2 (t 0 ) p B N −2 (t 1 )

p B N −1 (t 0 ) p B N −1 (t 1 )

. . . p B 1 (tnp )

. . .

. . .

B N −2 (tnp )

B N −1 (tnp )

0 0

0 0

0 0

. . .

. . .

0

1

. . . −2

p

p

p B N (t 0 ) p B N (t 1 )

⎤

⎥ ⎥ . ⎥ . ⎥ . ⎥ p B N (tnp ) ⎥ ⎥ ⎥ 0 ⎥ 0 ⎥ . ⎦ . .

1

Q0 ⎢ Q1 ⎥

PN

37

···

. . .. . . . . . p p B 1 (tnp ) B 2 (tnp ) · · ·

. . .

·⎢ . ⎣ ..

35

p B 2 (t 0 ) p B 2 (t 1 )

. . .

P0 ⎢ P1 ⎢

33

p B 1 (t 0 ) p B 1 (t 1 )

1

⎡

32

⎥ ⎤ ⎢ ⎢ .. ⎥ ⎢ . ⎥ ⎥ ⎥ ⎢ Q np ⎥ ⎥ ⎢ ⎥ ⇒ [M E ] · [ P ] = [ Q E ] ⎢ ⎥=⎢ ⎦ ⎢ 0 ⎥ ⎥ ⎢ 0 ⎥ ⎥ ⎢ ⎢ . ⎥ ⎣ .. ⎦

(9)

0

40

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

74

.. ⎥ · ⎢ .. ⎥ = ⎢ .. ⎥ .⎦ ⎣ . ⎦ ⎣.⎦

· · · 0 1 −2 1

1 −2

31

42

73

81 82

So, the LS system of equations of Eq. (2) can be completed to include the minimization of the second derivatives, producing an extended matrix M E of dimensions (np + N ) × ( N + 1) and also an extended independent term Q E with the form of Eq. (9).

30

41

72

⎡ ⎤

⇒ [ S ] · [ P ] = [0]

15

18

69

This set of N − 1 equations with N + 1 unknown control points, can be developed in matrix form easily as in Eq. (8):

⎢ .. ⎣.

17

68

(7)

P N −2 − 2P N −1 + P N = 0

12

16

67

.. .

3

7

To summarize the whole method, matrix formulation (8) contains an implicit fairing of the control polygon, and is included in the definition of the ILS fitting in Eq. (9) since these expressions are a function of the control points. A weighted combination of the derivatives/differences is used, and factor α inside the interval [0, 1] that counts the effect of the differences will also be an input to ILS fitting. A large α produces a softer distribution of the control points, but may reduce the effect of fitting, increasing the distances between the data points and the B-spline. This system is shown in Eq. (10) and is solved as explained in Section 3.



 (1 − α )[ M ] · [ P ] = (1 − α )[ Q E ] α[ S ]

(10)

Increasing α produces a control polygon without wiggles and with a nicer aspect, but the distance between the data points and the B-spline will be higher than with the use of a lower α value. The parameters t i corresponding to the data points that are used in Eq. (9) and Eq. (10) are iteratively calculated based on the minimum distance from the data points to the B-spline, as explained in Section 3. So, in order to obtain the fitting of the airfoil data with the desired tolerance points following the previous sections, the user has only to indicate the degree p, the number of parts P to subdivide the Bézier pieces that gives the non-uniform knot vector, and the value of α . A schematic flowchart of the complete method is presented in Fig. 12.

83 84 85 86

Fig. 12. Flowchart of the proposed method.

87 88

7. Validation examples

89 90

Different researchers have developed several methods for airfoil data fitting, usually as a previous step before optimization procedures. In this Section some airfoils found in the literature are fitted with the proposed methodology comparing the tolerance and the number of control points of the curves. A summary is produced in Table 1. The first two examples correspond to reference [15], where airfoils NASA SC(2)-0712 and NACA 64(1)-212 were fitted with a NURBS technique, considering two different cubic B-spline with N P = 9 control points each, one for the upper and a second one for the lower part. Authors of this reference have obtained a maximum deviation of 0.0019 m and a mean deviation of 0.0007 m for unit chord airfoils. Results of the present method with p = 3, P = 3 subdivisions in the leading edge and α = 0.6, produces a B-spline with N P = 14 control points, that fits the data points with a maximum distance of 0.00097 m and an average distance of 0.00038 m. This curve is presented in Fig. 13, including a curvature graph and a detail of the leading edge. The tolerance used (radius of the circles) is 0.001 m. If the tolerance is increased, the number of control points can be reduced by increasing the degree of the fitting B-spline. For example, p = 5, P = 3 and α = 0.5 produces a B-spline of N P = 10 control points, that fits the data with a tolerance of 0.002 m, and has Max/Med distance of 0.00174/0.00081 m. This is comparable with the results of reference [15], but considering only 10 instead of 2 × 9 control points. The curve is depicted in Fig. 14, together with a comparison of the curvature graphs of p = 3 and p = 5, and with a detail of the airfoil’s leading edge. In the case of the second airfoil of reference [15], equivalent results are obtained considering p = 5, P = 2 and α = 0.5, producing a B-spline with 9 control points that has a Max/med distance with the data points of 0.00153/0.00058 m, that are lower than the values presented in the mentioned reference and also with a smaller number of control points. This is depicted in Fig. 15. More examples of airfoil fitting can be found in reference [8]. The author of that work maintained the number of control points for B-splines of degree p = 5, and tested different airfoils. He considered a MH121 unit chord airfoil and obtained a Maximum distance of 0.0014 m considering 14 control points (Fig. 16). The presented method when applied to the mentioned airfoil using the same degree p = 5, P = 4 and α = 0.4 obtains a Max/med distances of 0.00092/0.00043 m with 11 control points.

91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

8

AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.8 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

1

67

2

68

3

69

4

70

5

71

Fig. 16. Fitting for MH121, p = 5.

6

72

7

73

8

74

9

75

10

76

11

77

12

78

Fig. 17. Fitting for MH43, p = 5.

13

79

14

80

15

Table 1 Comparison of the presented method with other ones (Intel Core I7 Processor, 3.4 GHz).

16 17

81 82 83

18

Airfoil 1m Chord

Reference

NPC

Max. Dist. (mm)

Comp. time (s)

84

19

NASA SC(2)-0712

[15] This [15] This [8] This [8] This

18 14 18 10 14 11 14 7

0.0019 0.00097 0.0019 0.0017 0.0014 0.00092 0.0016 0.00131

N.A. 3.3 N.A. 3.0 N.A. 3.1 N.A. 2.9

85

20

NACA 64(1)-212

21 22

MH121

23 24

MH43

25

paper paper paper paper

26

29

Fig. 13. Fitting for NASA SC(2)-0712, p = 3.

30 31 32 33 34 35 36 37 38

In the same reference, the author tested an MH43 one meter chord airfoil, obtaining with p = 5 and 14 control points, a Maximum distance of 0.0016 m. In the present method, considering the same degree p = 5, P = 4 and α = 0.5, a B-spline of 7 control points is obtained with a Max/med distance of 0.00131/0.00050 m. This airfoil is showed in Fig. 17. The same reference also fits an airfoil of the NACA64 family, quite similar to the one previously presented in Fig. 15, obtaining higher maximum distances than the ones of the present methodology. 8. Conclusions

39 40 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 59

Fig. 14. Fitting for NASA SC(2)-0712, p = 5.

60 61 62 63 64 65 66

88 89 90 91 93 94 95 96 97 98 99 100 101 102 103 104 105 106

41

58

87

92

27 28

86

Fig. 15. Fitting for NACA 64(1)-212, p = 5.

This paper has presented a practical method for producing a Bspline representation of a large number of airfoil data points. The method produces a fitting of the 2D airfoil data points allowing the leading edge, which has a major effect on the airfoil properties, to be well reproduced and with a good curvature distribution improving interpolating methods. The main mathematical reason to use B-spline curves is that they have enough flexibility to solve this particular problem of fitting a large number data points: notice that all the equation systems are linear and no dot need initial guess for the control points. The fitting also uses an iterative selection of the parameterisation, which allows setting to a given tolerance, which is very useful from a construction perspective. The method produced an airfoil with a minimum number of control points that fits the data points below a given allowance. The main reason for this reduction is the use of a non-uniform knot vector, as explained in Section 4, with the knots located near the leading edge. This way, the control points are automatically located near the high curved leading edge, and this part of the airfoil is well reproduced under a tolerance. This reduced number of control points permits a more simple surface definition, which is important in subsequent stages of the design process including numerical CFD calculations and optimization. When fitting a large number of offset data points, the B-spline may produce a curve with a control polygon that wiggles as a result of the solution of the large equation systems, see Eq. (2). In

107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:3702 /FLA

[m5G; v1.180; Prn:28/06/2016; 15:55] P.9 (1-9)

F. Pérez-Arribas, I. Castañeda-Sabadell / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8

the present method this is improved with the fairing conditions of the control polygon described in Section 6. The presented method was tested for airfoils with different characteristics and compared using various methods from the literature. It produced good results for the modelling a wide range or airfoil shapes and improves the maximum deviations from the data points using a lower number of control points than the available literature, as detailed at Section 7.

9 10

to the order of the curve at its ends, where order = degree + 1. In this way, the B-spline interpolates the ends of its control polygon at u = 0 and u = 1 and is tangent at its ends to the first and last segment of its control polygon. This last property simplifies the mathematical definition of the curves used in the method. Note that the derivative s (u ) of a B-spline is a linear combination of the derivatives of the basis functions:

None declared.

13 14

=

Acknowledgements

17 18 19

This work was partly supported by the Spanish Ministerio de Economía y Competitividad through research grant TRA2015-67788 -P, and by UPM Research Project RP160850004. Thanks to Sergio del Castillo Tello for Grasshopper support.

20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

Appendix To introduce the notation for this paper, a brief discussion of B-spline curves follows. A B-spline curve is formed by several pieces of polynomial curves called Bézier pieces, and the entire curve is C 2 (common curvature or second derivatives) at the junctions in the case of cubic (p = 3) B-splines. The curve is defined with a polygon called the control polygon and with an interpolation algorithm that allows for its construction by relating the curve to the control polygon. The interpolation steps are encoded in a family of piecewise polynomial functions B nj (u ) called B-spline functions of the nth degree and are calculated with Cox de Boor’s algorithm of Ref. [16]. A B-spline curve, s(u ) in Eq. (11), is a linear combination of basis functions with N + 1 control points P j as coefficients. Therefore, B-spline curves are parametric, x = X (u ), y = Y (u ), z = Z (u ), and the parameter u is normally considered between [0, 1]. In the plane, P j = ( X j , Y j ), j = 0, . . . , N generate a B-spline s(u ) of pth degree:

40 41 42

s (u ) =

45

=

48

N 

51 52 53 54 55 56



Xj ·

p B j (u ), Y j

·

p B j (u )

(11)

where the basis functions are obtained with Cox de Boor’s algorithm in Eq. (12):



49 50



p

P j · B j (u ) = X (u ), Y (u )

j =0

46 47

N  j =0

43 44

p





B 0j (u ) B nj (u )

= =

1 0

u ∈ [ u j −1 , u j ) u∈ / [ u j −1 , u j ) u − u j −1

u j +n−1 − u j −1 p

B nj −1 (u ) +

u j +n − u u j +n − u j

−1 B nj + (u ) 1

(12)

The basis function B j (u ) depends on the knot vector u j , which, in this work was chosen to be non-uniform with a multiplicity equal

68 69 70 71 72 73 75

P j · B j (u ) = X (u ), Y (u )

76

N 

77

X j · B j (u ), Y j · B j (u ) p

p

78

(13)

j =0

15 16

N  j =0

11 12

67

74

s  (u ) =

Conflict of interest statement

9

When the number of control points (N + 1) is equal to the order of the curve (p + 1), then the B-spline curve is formed from just one piece and can be called a Bézier curve of the pth degree. A variable presented in bold letters indicates that this element is formed from several components, such as points, vectors or B-spline curves, that may have X, Y and Z components. References

79 80 81 82 83 84 85 86 87 88 89 90 91

[1] D. Hoitsma, A note on surface inflections, Comput. Aided Geom. Des. 13 (1) (1996) 89–93. [2] N.M. Patrikalakis, T. Maekawa, Umbilics and lines of curvature, in: Shape Interrogation for Computer Aided Design and Manufacturing, Springer, 2002, pp. 231–264. [3] J.H.S. Fincham, M.I. Friswell, Aerodynamic optimisation of a camber morphing aerofoil, Aerosp. Sci. Technol. 43 (2015) 245–255. [4] M. Robitaille, A. Mosahebi, É. Laurendeau, Design of adaptive transonic laminar  θ t transition model, Aerosp. Sci. Technol. 46 (2015) airfoils using the γ − Re 60–71. [5] L.A. Piegl, W. Tiller, The NURBS Book, Springer, 1997, pp. 410–413. [6] A. Bentamy, F. Guibault, J.Y. Trépanier, Cross-sectional design with curvature constraints, Comput. Aided Des. 37 (2005) 1499–1508. [7] A. Bentamy, F. Guibault, J.Y. Trépanier, Wing shape optimization using a constrained NURBS surface geometrical representation, in: Proceedings of ICAS 2002 Conference, Toronto, 2002, pp. 123.1–123.11. [8] X. Mauclère, Automatic 2D airfoil generation, evaluation and optimisation using MATLAB and XFOIL, Master Thesis, Technical University of Denmark, 2009. [9] N.A. Vua, J.W. Lee, Aerodynamic design optimization of helicopter rotor blades including airfoil shape for forward flight, Aerosp. Sci. Technol. 42 (2015) 106–117. [10] L. Piegl, W. Tiller, The NURBS Book, second edn., Springer, 1997. [11] I. Juhosz, Weight-based shape modification of NURBS curves, Comput. Aided Geom. Des. 16 (5) (1999) 377–383. [12] I. Juhász, M. Hoffmann, Modifying a knot of B-spline curves, Comput. Aided Geom. Des. 20 (2003) 243–245. [13] R. Goldenthal, M. Bercovier, Spline curve approximation and design by optimal control over the knots, Computing 72 (2004) 53–64. [14] H. Nowacki, X. Lü, Fairing composite polynomial curves with constraints, Comput. Aided Geom. Des. 11 (1994) 1–15. [15] Yu Liang, XiaoQuan Cheng, ZhengNeng Li, JinWu Xiang, Multi-objective robust airfoil optimization based on non-uniform rational B-spline (NURBS) representation, Sci. China, Technol. Sci. 53 (10) (2010) 2708–2717. [16] G. Farin, Curves and Surfaces for CAGD, fifth edition, Morgan Kaufmann, San Francisco, 2001. [17] A. Sobester, A.I.J. Forrester, Aircraft Aerodynamic Design. Geometry and Optimization, Wiley, 2015.

92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132