A fast algorithm for inverse airfoil design using a transpiration model and an improved vortex panel method

TECHNICAL PAPERS

A fast algorithm for inverse airfoil design using a transpiration model and an improved vortex panel method

Denis R. Petrucci^I; Nelson Manzanares Filho^II

^[e-mail]. Universidade Federal do Recôncavo da Bahia. Centro de Ciências Exatas e Tecnológicas. 44380-000 Cruz das Almas, BA. Brazil

^IIEmeritus Member, ABCM. [e-mail]. Federal University of Itajubá - UNIFEI. Institute of Mechanical Engineering. 37500-903 Itajubá, MG. Brazil

ABSTRACT

A fast algorithm for inverse airfoil design using an efficient panel method for potential flow calculation is presented. The method employs linear vortex distributions on the panels and a consistent procedure for imposing the Kutta condition, thus eliminating the spurious aerodynamic loading that usually appears at a cusped trailing edge. The algorithm searches the airfoil ordinates attending to a given surface velocity distribution with fixed abscissas. It begins with a guessed starting shape and successively modifies it by an iterative process, such that the normal velocity vanishes and the calculated velocity distribution gradually approaches the required one. Each iteration is performed in two main steps: 1) the flow calculation step; 2) the geometrical marching step, where the calculated velocity distribution is compared with the required one and a transpiration model is applied to modify the current shape towards another one more close to the target shape. The geometrical marching is conducted by varying the panel slopes as a function of the normal velocity excess induced by the difference between the required and calculated velocities. A scheme is applied in order to close the body shape. Various test cases were carried out and are presented for the efficiency and robustness validation of the proposed inverse algorithm.

Keywords: inverse method, panel method, airfoils, vortex distributions

Introduction

Inverse methods are useful tools for aerodynamic shape design by virtue of their great flexibility. Airfoil shapes are constructed for attending each specific design situation, thus overcoming the relying on standard families of wing sections like the NACA series (Abbott and Doenhoff, 1959).

The inverse methods for aerodynamic shape design can be implemented in various manners. They also compete with other design methods. There is no general accepted classification for airfoil design procedures. For instance, Yiu (1994) classifies them in four broad categories: (1) inverse methods (including iterative correction and nonlinear system approaches); (2) iterative modification methods (including here the very important optimization approach); (3) transformed plane method (including conformal mapping techniques); (4) special methods (including panel methods for incompressible potential flows). There is some superposition and mixing in this classification: for example, pure inverse methods can actually be implemented by means of both conformal mapping and integral equation formulations (for whose solution the panel method represents a classical numerical approach).

The pure inverse methods search a unique body shape solution (if it exists) for attending a specified velocity or pressure distribution on the body contour. On the other hand, optimization methods normally search a best solution regarding to a specified global objective (like minimum drag or maximum lift) and taking into account some restrictions upon the final geometry or flow conditions. Generally, the computational cost of the inverse approach tends to be a small fraction of the optimization approach cost. Thus, the inverse approach is indicated for the earlier stages of the aerodynamic design when some prospective shapes are delineated for later refinements by means of the optimization approach.

For ideal flows, pure inverse methods can be efficiently implemented by means of conformal mappings, streamline coordinate-based transformations or boundary integral formulations, particularly the panel method. In the conformal mapping method, the body shape is searched in a transformed plane where it usually approaches a circle by a suitable construction of parametric mapping functions. This procedure allows a rigorous geometric control during the iterative process. Selig and Maughmer (1992a, b) proposed an inverse airfoil design method based on the multi-point approach due to Eppler (1990). It is possible to adapt the basic procedure for including geometric restrictions such as maximum thickness and camber as well as boundary layer growth criteria upon the required velocity distributions. The procedure was also extended by Selig (1994) to the inverse design of cascade of airfoils.

The great advantage of conformal mapping methods is their ability to express conditions of closeness and geometrical uniqueness in a natural manner. However, they also have some drawbacks: They are not generally extensible to the three-dimensional case, rely on rather sophisticated mathematical techniques and require a prior knowledge of a suitable mapping function in order to handle each particular problem.

Some limitations of conformal mapping techniques can be overcome by using streamline coordinate-based transformations. In this context, Barron (1990) has proposed an elegant non-iterative method for airfoil design in potential incompressible flows. The method is based on a previously developed formulation employing a von Mises coordinate transformation (Barron, 1989). This approach was extended to general ideal steady flows (Latypov, 1993; Yiu, 1994), but it is also restricted to two-dimensional configurations.

Boundary integral formulations represent a flexible alternative in developing inverse airfoil design procedures, in especial when they are numerically implemented by means of panel methods. A procedure of this type was proposed by Shigemi (1985) for single and multi-element airfoils. In this procedure, straight panels with linear vortex distributions are employed and the Neumann boundary condition is applied in control points. In this case the final vortex distribution is known a priori since it is equivalent to the required contour velocity. The unknown airfoil ordinates are determined by applying the Newton-Raphson method for solving the system of nonlinear algebraic equations that arises. In order to assure the airfoil closeness and also to fix the body edges, Shigemi (1985) applied a least square technique.

Petrucci et al. (1998) proposed an inverse methodology with a modified version of the Hess-Smith method (Hess and Smith, 1967) for the flow calculation step. The method employs straight panels with uniform source distributions of unknown intensities and a vortex distribution of sinusoidal shape with an unknown maximum intensity at the leading edge and a vanishing intensity at the trailing edge. This vortex distribution alleviates spurious aerodynamic loads that usually occur at a cusped trailing edge. The source intensities and the maximum vortex intensity are calculated by applying the Neumann boundary condition at the control points and the Kutta condition at the trailing edge. The geometric marching step is conducted by altering the panel slopes according to a procedure suggested by Murugesan and Railly (1969). Due to the low order of the singularity distributions and also the poor performance of the source distributions in controlling curvature effects, Petrucci (1998) has verified some convergence difficulties in the case of airfoils with very cusped trailing edge and high camber. It was concluded that a higher order panel method would be really necessary in order to overcome these difficulties.

In this paper, one proposes an improved algorithm for inverse airfoil design with basis on the panel method. In the flow calculation step one employs a linear vortex panel method with the Neumann boundary condition and a new scheme for applying the Kutta condition consistently in the case of a cusped trailing edge (Petrucci et al., 2001). In the geometrical marching step, one applies the transpiration model of Muregesan and Railly (1969) with modifications for fixing the trailing edge point and the airfoil abscissas during the iterative process and also for assuring the body closeness.

The effectiveness and robustness validation of the proposed inverse algorithm is carried out by means of some case studies for which analytical results are available in the literature (circular cylinder and Joukowski airfoils).

Nomenclature

a = circle radius

C_D = drag coefficient

Cp = pressure coefficient

ft = accelerating factor

j_st = index of point closest to leading edge stagnation point

m = number of panels

me = eccentricity magnitude

l = chord length

s = natural coordinate along the body contour

S = normalized natural coordinate along the body contour

sl = total countour length

W_a = calculated (analyzed) velocity

Wc_n = calculated normal velocity

W_n = effective normal velocity

W_r = required velocity

W_t = tangential velocity

W_bf = trailing edge velocity

W_nbf = normal velocity at trailing edge

W_tbf = tangential velocity at trailing edge

Wny_j= y-projection of the effective normal velocity of the j-th panel

= conjugated complex velocity

= conjugated complex velocity of the onset uniform flow at infinity

x = abscissa

X = normalized abscissa

y = ordinate of actual body shape

Y = ordinate of new body shape

= effective new ordinate

z = complex coordinate

Greek Symbols

a = angle of attack between the onset flow direction and the chord axis

b* = angle of camber

c_j = angle between the cj-th panel and the x direction

dy = ordinate variation

DS = panel length

Dy = ordinate difference of current body shape

DY = ordinate difference of new body shape

e_st = stopping tolerance

g = vortex intensity

g_fict = fictitious vortex intensities

_fict = fictitious vortex intensities with acceleration

l = filter computed as the mean value of the normal velocity modulus

q = panel angle

z(s) = integration point on the body contour

Subscripts

bf = relative to trailing edge

j = relative to panel j

n = relative to normal velocity

t = relative to tangential velocity

st = relative to stagnation point

Algorithm Description

Given a required velocity distribution at the airfoil contour and a starting body shape, the inverse design algorithm searches the required body shape by an iterative process. Each iteration is subdivided into two main steps: 1) the flow calculation step; 2) the geometrical marching step. At the end of each iteration one finds a new body shape that is presumably more close to the target shape. The iterative process is repeated until the differences between the ordinates of the current and last iterations stay within a prescribed tolerance. In what follows, these iterative steps are described in more detail. A flow chart of the algorithm is shown in Fig. 1.

Flow Calculation Step

In this work, the corresponding flow field for each iterative body shape is computed by means of an efficient panel method based on linear vortex distributions. The flow is assumed to be two-dimensional, steady, potential and incompressible. The body contour is approximated by m straight panels with z₁, z₂,..., z_m, z_m+1 denoting the panel nodes and z₁= z_m+1 representing the trailing edge (Fig. 2). In each panel, the central point is designated as a control point for applying the Neumann boundary condition (of zero normal velocity). Linear vortex distributions are placed on the panels such that the vortex intensity associated with the cj-th node is g_j. This value is equivalent to the required velocity at the cj-th node.

The conjugated complex velocity in a generic point z = x + iy of the complex plane, i = (- 1)^1/2, is expressed by means of a Cauchy integral on the boundary contour C:

where represents the conjugated complex velocity of the onset uniform flow at infinity, a being the angle of attack between the onset flow direction and the chord axis; z(s) represents an integration point on the body contour. The integral is calculated by summation of each sub-integral value corresponding to the contribution of the panels joined by the node z_j (Fig. 2a). Assuming straight panels with linear vortex distributions, this contribution can be calculated analytically, resulting the following expression for j = 2, ..., m (c _j is the angle between the ccj-th panel and the x direction)

The panels joined by the trailing edge z₁ = z_m+1 receive a special treatment (Fig. 2b). A regularization condition is imposed to make the trailing edge velocity finite. Thus the equation system resulting from application of the Neumann boundary condition in each panel becomes determined (m x m). However, the regularization condition is not sufficient for satisfying the Kutta condition consistently in the case of a very cusped trailing edge. Thus one further requires the total velocity at the trailing edge to be equal to the vortex intensity there, as it should theoretically be, i. e., . In this way, the system of equations becomes over-determined (m + 1 equations with m unknowns g ₁, ..., g _m) being solved by least squares (Petrucci et al, 2001).

The tests made up to now have demonstrated the versatility and good precision of this procedure for airfoils with very cusped trailing edges, even with a relatively low number of panels. Figures 3 and 4 show some coefficient pressure results for Joukowski airfoils with and without camber, respectively, using m = 16 panels only. The abscissas are normalized with respect to the airfoil chord. The airfoils were generated by conformal mappings from the circle with the following parameters: angle of attack, a, angle of camber, b*, eccentricity ratio, a/me, a being the radius of the circle and me the eccentricity magnitude (Karamcheti, 1980). The results are compared with the low order panel method of Hess-Smith modified by Petrucci (1998).

One can see that the low order panel method is not able to represent adequately the pressure distribution with this low number of panels, mainly for the cambered airfoil (Fig. 4). On the other hand, the present method leads to satisfactory results in both cases. It is important to mention that increasing the number of panels m — as in fact becomes necessary in applications — will produce results that converge rapidly to the analytical results. Drag coefficients C_D(m) were computed to be compared with the zero analytical value. Some results are C_D(16) = 0.01077, C_D(32) = 0.00117, C_D(64) = 0.00011, C_D(128) = 0.00000 for the symmetrical airfoil and C_D(16) = 0.01896, C_D(32) = 0.00415, C_D(64) = 0.00097, C_D(128) = 0.00009 for the cambered airfoil. For the Hess & Smith method with m = 128 one has C_D(128) = 0.00143 for the symmetrical airfoil and C_D(128) = 0.00114 for the cambered airfoil.

Geometrical Marching Step

At the end of the flow calculation step one has obtained a new vortex distribution. This could be compared with the required vortex (velocity) distribution in order to decide if one has attained convergence or not. However, before verifying the convergence, one proceeds directly with a geometrical marching step by constructing another shape that is more close to the target shape.

Here one applies a modified version of the transpiration model due to Murugesan and Railly (1969). In this model, a vortex-only boundary method is employed and the transpiration effects are represented by normal velocities induced by fictitious vortex representing the difference between current and target velocities. The introduced modifications were first presented by Petrucci (1998). They do not significantly alter the basic transpiration model and were mainly proposed for maintaining the trailing edge fixed during the iterative process and also for assuring a closed body shape at the end of the geometrical marching step. In the original model, the leading stagnation point of the iterated shape must be found and a re-paneling of the intrados and extrados is required in any iteration, starting from the stagnation point towards the trailing edge. Further, the body closeness issue is not apparently addressed by Muregesan ans Raily (1969). Another modification is also applied for fixing the body abscissas between iterations and so facilitating the global convergence. With these modifications, the angle of attach of the onset flow do not need to be given since the iterative process is able to set the correct body attitude relative to the onset flow at convergence by itself. This issue can be very useful for design purposes. Also, a re-paneling of the body surface becomes unnecessary.

Basically, the geometrical marching step is as follows: firstly, one calculates a fictitious vortex distribution in terms of the difference between the required and calculated velocities (vortex intensities) at the panel nodes. Then one calculates a normal velocity excess (transpiration) induced by this fictitious vortex distribution. Finally, one conducts a suitable variation of the panel slopes trying to annul this normal velocity excess. A straightforward scheme is further applied in order to assure the body closeness. In what follows, the whole procedure will be described in more detail.

Calculation of the Fictitious Vortex Distribution

The intensity of the fictitious vortex distribution depends upon the difference between the calculated (analyzed) velocity W_a(s) and the required (target) velocity W_r(s). However, one must be careful about the velocity signal relative to the adopted path orientation around the body. Here the path starts and ends at the trailing edge, and its orientation is such that the interior of the body remains at the right side (Fig. 5). Normally one has a stagnation point next to the leading edge. Thus, the velocity is negative on the bottom portion of the body (starting from the trailing edge up to the stagnation point) and positive on the top portion (following from the stagnation point back to the trailing edge.)

Initially, one finds the nodal point j_st closest to the stagnation point corresponding either to the calculated velocity or to the required velocity, that one which first occurs. That is, we take the first point where these velocities have different signals or at least one of them is null (with exception of the trailing edge) and use this point to divide the boundary into two branches:

Then one computes the fictitious vortex intensities at the panel nodes by the following expressions:

Acceleration of the Iterative Process

Some providence can be taken at this point in order to accelerate the global iterative process. First, it was verified that a constant accelerating factor ft can be applied on the fictitious vortex distribution without altering its basics characteristics. Tests have indicated that this factor must lie in a conservative range, 1 < ft < 3.5. Values much above 3.5 may produce catastrophic oscillations in the iterative process while retarding factors ft <1 have not needed to be applied. One writes:

Calculation of the Induced Normal Velocities

The fictitious vortex distribution will induce a normal velocity distribution on the body contour which is calculated using the current influence matrix or preferably the starting influence matrix (this issue will be later discussed in detail). It was observed that this normal velocity may suffer severe variations especially near the leading edge where the fictitious vortices may also vary intensely. In producing a new shape, these variations would cause the appearance of spurious bumps or concavities at the leading edge region which impair the convergence process. For avoiding these drawbacks, one applies an automatic filter l to the calculated normal velocities Wc_n(s) by bounding their excessive values and maintaining those that are sufficiently mild. The filter is automatically computed as the mean value of the normal velocity modulus:

where sl represents the total contour length. The filter acts as an upper bound to the modulus of the effective normal velocities W_n(s) such that

Generation of a New Geometrical Shape

Using the effective normal velocity in each control point, a scheme was conceived in order to alter the panel slopes, i. e., the body ordinates, without varying the body abscissas. For the current body shape the ordinate difference between the nodal points of cj-th panel is Dy_j = y_j+1 – y_j. For the new body shape the corresponding ordinate difference DY_j is calculated by adding a suitable ordinate variation dy_j as follows

According to Fig. 6a, a relation should exist between the ordinate variation, the required velocity and the normal velocity calculated in the current iteration. Since one seeks for a new body shape with zero normal velocity on it, the dashed line on bottom of Fig. 6a represents an approximation to the new tangential direction. But this direction is also geometrically represented by the new panel slope which is related to the current panel slope and the ordinate variation. Thus, a similarity exists between the cinematic and geometric triangles in Fig. 6a, from which the required ordinate variation follows

Starting from the first panel, a cumulative summation of ordinate variations calculated in Eq. (11) is performed for the nodal points j = 2, ..., m+1. Thus, the new ordinates becomes

Procedure for Closing the New Body Shape

The application of Eq. (12) does not assure that the resulting shape is closed as required. The ordinate of the last node m+1 do not necessarily coincides with that of the first node, m =1. To overcome this drawback, one divides the difference Y₁ – Y_m+1 by the number of panels and adds this remaining value to the nodal points cumulatively, starting with the 2^nd node and ending with the node m+1. This scheme assures the body closeness (Fig. 6b). Thus the effective new ordinates becomes

Stopping Criteria for the Iterative Process

At least two kinds of stopping criterion can be applied. One is based on the difference between the calculated and required velocities at the end of the flow calculation step. Normally such type of criterion is not adequate for the present algorithm, since the velocity distribution may exhibit a very large variation at the leading edge region for very small geometrical variations. The same may occur with a criterion based on the velocity difference between successive iterations.

A more adequate and realistic criterion can be formulated using the geometric marching itself. In this paper one defines the stopping criterion using the mean quadratic error of the differences between the ordinates of the current iteration and those of the last iteration. The iterative process stops when this value becomes smaller than a prescribed tolerance, e_st. At convergence, the calculated velocities approach the required ones, the fictitious vortices and induced normal velocities approach zero and the current shape approaches the target shape.

Use of the Starting Influence Matrix during the Geometrical Marching Step; Further Discussion about the Computations

Some tests have raised another important convergence issue: in calculating the normal velocities induced by the fictitious vortex distribution, convergence improvements can be attained if one uses the influence matrix of the starting shape instead of the current iterated shape. At first glance, this may sound strange because the final and starting shapes may be very different from each other and the influence matrix is strongly connected to the shape of the airfoil.

For explaining this situation, it is important first to note that the current matrix is always used in the flow calculation step. The normal velocities induced by the actual vortex distribution must counterbalance the normal effect of the uniform onset flow. Roughly speaking, being [A] the current influence matrix (of the current shape) and (W¥ _n) the vector of onset flow normal velocities, the corresponding vector (g ) of actual vortex intensities (equal to tangential velocities) is calculated by solving the system [A](g ) = - (W¥ _n) in the flow calculation step. Due to the Kutta condition treatment, the original system of equations is modified and it is solved by least squares, as already discussed before. But in synthesis, the current influence matrix must somehow be retained in order to obtain a well calculated flow field that indeed represents the current airfoil shape.

On the other hand, in the geometrical marching step, a vector of fictitious vortex intensities is calculated, as previously explained. Now, the actual uniform onset flow is not present and the vector (W_n) of induced normal velocities is calculated by a simple matrix multiplication: (W_n) = [A] . Here, it is easy to see that the use of the current matrix is not strictly necessary, provided the algorithm is convergent. Indeed, in this case both vectors and (W_n) will converge to zero and the matrix [A] will eventually become irrelevant during the iterative procedure.

The reason why the use of the starting matrix instead of the current matrix has been shown advantageous is more subtle. One hypothesis is that the successive shape alterations would impair the conditioning properties of the influence matrix with resulting error propagation. Some tests were made for showing that this hypothesis is really appealing. For these tests, one has chosen as the target shape a Joukowski airfoil with conformal mapping parameters b * = 12º , a/me = 4.5, a = 4º, and as the starting shape an ellipse with aspect ratio equal to 0.1. The stopping tolerance was taken as e_st = 10- ⁴. Tests were made for m = 24 panels and m = 50 panels, with the accelerating factor ft set equal 2.2.

Figures 7 (a) and (b) show that the condition number of the current matrix increases rapidly and exhibits some oscillations during the first iterations. This behavior is intensified when the number of panels is increased. These results serve to explain why the use of the starting influence matrix helps to accelerate the convergence process: Even more important than the smaller condition number of the starting matrix, the constancy of its spectral properties implies a faster shape evolution in comparison with the use of the current matrix. More specifically, the oscillations in the condition number of the current matrix impart oscillations in the normal velocities induced by the fictitious vortex distribution which in turn generate oscillations in the iterated shape and finally a relative delay in the convergence process. This delay can be clearly observed in Figs. 7 (c) and (d) where it is shown the iterative evolution of the mean value of the normal velocity modulus l (the filter defined in Eq. 7). The "convergence points", corresponding to the iteration for which the stopping tolerance is achieved, are indicated by black circles. When using the starting matrix, the convergence is attained in 18 iterations for m = 24 panels and in 30 iterations for m = 50 panels; on the other hand, when using the current matrix the number of iterations for convergence practically doubles — to 34 and 68, respectively. As now expected from Figs. 7 (a) and (b), it is also clear from Figs. 7 (c) and (d) that the starting matrix leads to smoother convergence behavior in comparison with the current matrix.

 

Some Remarks about the Generality of the Proposed Inverse Algorithm

In this work, the vortex panel method was applied for the whole iteration process. Nevertheless, it is important to note that the flow calculation step is independent of the geometrical marching step. It is conceivable the application of other flow models in the first step, while maintaining the proposed approach for modifying the body geometry in the second step. It is possible to use viscous or compressible flow solvers in the first step as long as one applies a consistent procedure for defining an equivalent velocity at the boundary contour. For example, suppose that a viscous incompressible method based on primitive variables is being applied. It is possible to use the Bernoulli theorem for converting the target and calculated pressure distributions into "equivalent potential velocities distributions" for computing the fictitious vortex distribution required in the geometrical marching step.

Further, the vortex-only panel method here used can be extended to three dimensions using a dipole-only panel method (Morino and Kuo, 1974): The dipole intensity is equivalent to the potential function on the boundary and its tangential derivative is equivalent to the boundary velocity (which in turn is equal to the vortex intensity in two dimensions). The use of vortex-only (2D) or dipole-only (3D) formulations can advantageously replace the source-based formulations used in direct solutions (Hess and Smith, 1967) and inverse design (Malone, 1982) since the addition of vortex or dipoles to the source distribution is always necessary for representing lifting-bodies like airfoils and wings.

Implementation and Validation of the Methodology

Computer Implementation Issues

For implementing the previously described algorithm, a Fortran double precision program was developed using the MS Fortran PS 4.0® compiler. The least square subroutine DLSQRR of the IMSL® library was called for solving the equation system in each flow calculation step. For the condition number calculations presented in Figs. 7 (a) and (b), the subroutine DLFCRG of the same library was called. The running tests were made on an IBM PC-like computer with a 350 MHz Pentium II® processor and 64 Mb RAM. Using 50 panels, the average CPU time per iteration was about 45 ms.

Case Studies

Some case studies were carried out in order to validate the proposed methodology and to evaluate its characteristics of convergence and robustness. The circular cylinder and the Joukowski airfoil were chosen as the target shapes since they have analytical results available in the literature, so being suitable for benchmarking. In all of the cases, the adopted starting shape (initial body) was an ellipse with aspect ratio equal to 0.1, and the stopping tolerance was taken as e_st = 10- ⁴.

Four cases will be shown here. According to target shape, number of panels, m, and adopted acceleration factor, ft, they are: (1) circular cylinder without circulation, m = 24, ft = 3; (2) symmetrical Joukowski airfoil at no incidence (b * = 0º, a/me = 12,5, a = 0º), m = 24, ft = 2.1; (3) cambered Joukowski airfoil at incidence (b * = 12º, a/me = 4.5, a = 4º), m = 24, ft = 2.1; (4) cambered Joukowski airfoil at incidence (b * = 12º, a/me = 4.5, a = 4º), m = 50, ft = 2.1. The nomenclature for conformal mapping parameters of the Joukowski airfoils follows that used by Karamcheti (1980).

Graphs for these cases are presented in Figs. 8 to 15, showing the convergence behavior for shape, velocity distribution W(s) and ordinate distribution y(s). In all of the graphs, the target shape and required (target) distributions are always represented by solid lines. The current iterated shapes are represented by dashed lines. The current velocity and ordinate distributions are indicated by proper keys within the corresponding graph. For the Joukowski airfoils, the ordinates were scaled up for better visualization of the iterative process for shape (the real shapes are relatively thin). The shape coordinates (x, y) and the boundary coordinate s were normalized with respect to the airfoil chord length l and the total contour length sl, respectively (X = x/l, Y = y/l, S = s/sl).

 

Discussion

One observes that a relatively small number of iterations was necessary for attaining geometric convergence. The symmetrical bodies (cases 1 and 2, Figs. 8 to 11) required less iteration than the non-symmetrical body (cases 3 and 4, Fig. 12 to 15). One also observes that the number of panels has a strong effect on the convergence ratio: for the same target airfoil, the use of 24 panels required 18 iterations (case 3, Figs. 12 and 13) while the use of 50 panels required 30 iterations (case 4, Figs. 14 and 15). Obviously, there is a precision improvement as the number of panels increases: This is a reminder about the need of a panel method able to achieve precise results with relatively few panels. Of course, a sufficiently high number of panels must be chosen for achieving an adequate geometric representation. For airfoils, m = 50 may be considered a recommended inferior bound.

Petrucci (1998) have previously employed a source panel method of low order for solving the inverse problem. This author has observed convergence difficulties in treating airfoils with high camber and/or with very cusped trailing edges. Velocity oscillations and geometrical crossings were then observed at the trailing edge region. These difficulties have arisen in solving the inverse problem even when the direct problem was satisfactorily solved with the same number of panels. It was observed that the source distributions made a poor job in controlling the airfoil curvature. It was concluded that both the type and the order of the singularity distributions must be chosen with care in order to adequately solve inverse design problems. In fact, one verifies in the present paper that the use of linear vortex panels and a consistent application of Kutta condition have decisively improved the iterative process.

One observes in cases 3 and 4 (Figs. 12 and 14) that the iterative process is able to adjust the correct angle of attack at convergence. In these cases, the target shape was set at a = 4º while the starting shape was set at a = 0º. The recovering of the correct angle of attack increases the flexibility of the inverse design approach.

Other important aspect is the departure between the starting and target shapes. Although the algorithm has exhibited a good global convergence performance one should be careful about guessing the starting geometry. In case 3, for instance, if one starts from the circular cylinder instead of the thin ellipse, the required number of iterations for convergence will increase from 18 to 235. This is an extreme situation: the algorithm is asked to search a thin and cusped airfoil starting from a very blunt body without an effective trailing edge. A slow-down in the convergence ratio is natural in this case due to the severe geometric adjustments to be performed, mainly at the trailing edge region. But the algorithm after all converges, exhibiting a global convergence capability and robustness.

In all of the cases, the calculation of the normal velocity distribution induced by the fictitious vortices was made by means of the starting influence matrix. This has substantially accelerated the iterative process in all of the cases. The numbers of iterations required for cases 1, 2, 3 and 4 were 8, 6, 18 and 30, respectively. On the other hand, when the influence matrix of the current iteration was used, the numbers of iterations for convergence were 18, 15, 28 and 50, respectively. But the characteristic of global convergence still remains.

With regard to the accelerating factor ft one has verified that it exerts a strong effect on the convergence ratio. In all of the cases the values of ft were chosen in order to produce a convergence ratio next to the best. For no acceleration (ft =1), the number of iterations would be doubled in average.

Concluding Remarks

A fast algorithm for inverse airfoil design with basis on linear vortex panels has been proposed. The presented test results have validated the algorithm in terms of versatility, convergence ratio and good precision even with a relatively low number of panels.

The convergence of the algorithm is fast when the departure between the starting and target shapes is not too large but it may be slowed-down in other cases. Nevertheless, in all of the tests the algorithm has exhibited good characteristics of global convergence. One can take advantage of this aspect in combining the proposed algorithm with methods having strong local convergence behavior, like the Newton-Raphson method.

The proposed algorithm is able to recover the correct attitude of the target airfoil with respect to the onset flow (the angle of attack) by itself. This issue increases the flexibility of the inverse design approach.

Finally, it is important to remark that the geometrical marching step of the proposed algorithm can be used in conjunction with other flow solvers, not necessarily of potential type. By using this approach, it would be possible to accelerate inverse design methods based on viscous or compressible flow codes. Thus the proposed algorithm can be useful in more realistic design situations.

Acknowledgments

The development of this work has been supported by the CNPq — National Council of Scientific and Technological Development of Brazilian Government — by means of a doctoral fellowship.

Paper accepted January, 2007.

Technical Editor: Francisco Ricardo Cunha.

• Abbot, I. H. and von Doenhoff, A. E., 1959, "Theory of Wing Sections", 2nd Edition, Dover Publications, Inc., New York.
• Barron, M. R., 1989, "Computation of incompressible potential flow using von Mises coordinates", Mathematics and Computers in Simulation, Vol. 31, pp. 177-188.
• Barron, M. R., 1990, "A non-iterative technique for design of aerofoils in incompressible potential flow", Communications in Applied Numerical Methods, Vol. 6, pp. 557-564.
• Eppler, R., 1990, "Airfoil Design and Data", Spring Verlag, Berlin-Heidelberg.
• Hess, J.L. and Smith, A.M.O., 1967, "Calculation of potential flow about arbitrary bodies", Progress in Aeronautical Sciences, Vol. 8, pp. 1-138.
• Karamcheti, K., 1980, "Principles of Ideal-Fluid Aerodynamics", R. E. Krieger Publishing Company, Florida.
• Latypov, A.M., 1993, "Numerical Solution of Steady Euler Equations in Streamline-Aligned Orthogonal Coordinates", IMA Preprint Series, No.1182, University of Minnesota.
• Malone, J. B., 1982, "A subsonic panel method for iterative design of complex aircraft configurations", Journal of Aircraft, Vol. 19, No. 10, pp. 820-825.
• Morino, L. and Kuo, C. C., 1974, "Subsonic potential aerodynamics for complex configurations: a general theory", AIAA Journal, Vol.12, No. 2, pp. 191-197.
• Murugesan, K., Railly, J. W., 1969, "Pure design method for aerofoils in cascade", Journal of Mechanical Engineering Science, Vol. 11, No. 5, pp. 454-467.
• Petrucci, D. R., 1998, "Inverse problem for the flow around isolated airfoils and turbomachine cascades (in Portuguese)", M. Sc. Dissertation, EFEI, Itajubá-MG, Brazil.
• Petrucci, D. R., Manzanares Filho, N. and Oliveira, W., 1998, "A numerical technique for solving the inverse problem of potential flow in turbomachine cascades (in Portuguese)", Proceedings of ENCIT 1998 - 7th Brazilian Congress of Thermal Engineering and Sciences, ABCM, Rio de Janeiro-RJ, Brazil, pp. 1305 1310.
• Petrucci, D. R., Manzanares Filho, N. and Ramirez, R. G. C., 2001, "An efficient panel method based on linear vortex distributions for flow analysis in turbomachinery cascades (in Portuguese)", Proceedings of COBEM 2001, 16th Brazilian Congress of Mechanical Engineering, ABCM, Uberlândia, MG, Brazil, Vol. 8, pp. 256-265.
• Selig, M. S., Maughmer, M. D., 1992a, "Multipoint inverse airfoil design method based on conformal mapping", AIAA Journal, Vol. 30, No 5, pp. 1162-1170.
• Selig, M. S., Maughmer, M. D., 1992b, "Generalized multipoint inverse airfoil design", AIAA Journal, Vol. 30, No 11, pp. 2618-2625.
• Selig, M. S., 1994, "Multipoint inverse design of an infinite cascade of airfoils", AIAA Journal, Vol. 24, No 4, pp. 774-782.
• Shigemi, M., 1985, "A solution of an inverse problem for multi-element aerofoils through application of panel method", Trans. Japan Soc. Aero. Space Sci., Vol. 28, No. 80, pp. 97-107.
• Yiu, K. F. C., 1994, "Computational methods for aerodynamic shape design", Mathl. Comput. Modelling, Vol. 20, No. 12, pp. 3-29.

Publication Dates

• Publication in this collection
  07 Apr 2008
• Date of issue
  Dec 2007

History

• Received
  Jan 2007
• Accepted
  Jan 2007