Testing a Good Error Function

First on my to-do list in my last posting was writing a good function for evaluating how closely a particular Joukowski airfoil fit an arbitrary user’s foil. I had done this before (the great crash ), and in that process learned what is now obvious to me: in marching CCW around the points on an airfoil/circle, one must inevitably progress from the first quadrant to the second to the third to the fourth. Obvious, yes? But unless you specifically code it this way, and instead, say, rely on your airfoil coordinates to determine which root of the inverse Joukowski transformation to use, then you will often find yourself dealing with some bizarre results, such as the lower surface of undercambered airfoils being erroneously assigned to the first quadrant. That one took me a long time to figure out, and it was only Neil Toronto’s new plot package that let me figure it out. In fact, those plotting routines are so simple and useful, that they are integral not only to zfoil, the product, but zfoil, the program, including associated test code.

Case in point, my newest function, (airfoil-fit-error N a c0). N is a "normalized" vector of airfoil coordinates, a the parameter defining the Joukowski transformation, and c0 the center of the circle that the transform of which is being compared to N. Applying this function to the NACA 4412 gives us the graph:

Although I am writing a more objective test, this graph was sufficient to give me confidence in my code, and insight into the next step in writing my program, coding an optimization routine to find the best value of c0, the best point to use for the center of a circle that will map to the user’s airfoil. The skew in the contour of the error function suggests that there’s some benefit in estimating the gradient, and as it happens,while getting my master’s degree, I came up with a 2D special case of Powell’s Method that should work quite well here. That’ll be a later post, for sure.
Re "explicitly coding it this way" (quadrant by quadrant), I ended up trying it two ways, via Racket’s "comprehensions" and "old school" with a sequence of named lets, basically unrolling the loop on the quadrant variable.
Method 1, using for/fold:

(define (airfoil-fit-error N a c0)

(let[(R(R a c0))

(q 1)]

(for/fold[(sum 0.0)]

[(n N)]

(cond

[(= q 1)

(cond[(< (real-part n)0)

(set! q 2)])]

[(= q 2)

(cond[(< (imag-part n)0)

(set! q 3)])]

[(= q 3)

(cond[(> (real-part n)0)

(set! q  4)])])

(let*[(c(cond

[(= q 1)(J^-1a a n)]

[(= q 2)(J^-1b a n)]

[(= q 3)(J^-1a a n)]

[(= q 4)(J^-1b a n)]))

(j(J a (+ c0 (make-polar R (angle (- c c0))))))

(k(- n j))]

(+ sum (sqr (magnitude k)))))))

Method 2/0 (0 being the "lost version", as best I can recollect):

(define (airfoil-fit-error N a c0)

(let[(r(radius a c0))]

(let loop1[(error 0.0)

(i 0)]

(let[(n(vector-ref N i))]

(if(>(real-part n)0)

(loop1(+ error

(let*[(c(J^-1a a n))

(j(J a (+ c0 (make-polar r (angle (- c c0))))))]

(sqr (magnitude (- n j)))))

(add1 i))

(let loop2[(error error)

(i i)]

(let[(n(vector-ref N i))]

(if(>(imag-part n)0)

(loop2(+ error

(let*[(c(J^-1b a n))

(j(J a (+ c0 (make-polar r (angle (- c c0))))))]

(sqr (magnitude (- n j)))))

(add1 i))

(let loop3[(error error)

(i i)]

(let[(n(vector-ref N i))]

(if(<(real-part n)0)

(loop3(+ error

(let*[(c(J^-1a a n))

(j(J a (+ c0 (make-polar r (angle (- c c0))))))]

(sqr (magnitude (- n j)))))

(add1 i))

(let loop4[(error error)

(i i)]

(if(< i (vector-length N))

(let[(n(vector-ref N i))]

(loop4(+ error

(let*[(c(J^-1a a n))

(j(J a (+ c0 (make-polar r (angle (- c c0))))))]

(sqr (magnitude (- n j)))))

(add1 i)))

error)))))))))))))

I first went for the comprehension to "modernize" my code, but IMHO it’s less readable. Truly "modern" Racket would probably use local defines to get rid of the indentation associated with the nested lets, but that never bothered me, so screw it.

The Big Idea

My target audience is probably familiar with the Joukowsky transformation:

\[n = c + {{a^2}\over{c}}\]
Note that I have replaced the traditional variable names (ok, avoided picking from among the many traditional variable names) with my own: \(c\) represents "cylinder space", where calculating inviscid flows are easy, and \(n\) represents "normal space", where "user space" has been rotated and translated to place the center of the airfoil on the origin, and the chord in line with the x axis. For the remainder of this article, the above relation will be referred to as \(J(c)\). This transformation is well known for transforming a cylinder of radius \(a\) centered on the origin onto a flat plate extending from \(-2a\) to \(2a\), and for transforming cylinders centered near the origin, and passing through the point \(c=a\), into airfoil like shapes.

Abbott and Van Daenhoff, in "Theory of Wing Sections", demonstrated that the Joukowsky transformation could be extended as a Laurent series to create a transformation that can describe an arbitrary airfoil:
\[T(c) = c + \sum_{i=1}^\infty {{a_i} \over {c^i}}\]
This transform is nifty in that it is analytic, and formulas exist to calculate all those coefficients, but there are all those cooefficients, and the series is not typically at all well behaved (I know, I’ve tried brute force, and when that failed, pondered the significance of large negative powers of small numbers). So calulating \(T\) directly is out, but we know that an anlalytic \(T\) exists, and that is enough. Consider
\[T(c) = J(c) + K(c)\]
\(K\) is my "correcting" transform. Since \(T\) and \(J\) are analytic, so must \(K\) be. Assuming c describes a suitable cylinder in "cylinder space", \(T(c)\) describes our airfoil, \(J(c)\) describes a similar Joukowski airfoil, and \(K(c)\) is the difference between the two.

It’s easy to describe the flow around the cylinder. It’s easy to transform that into the flow about the corresponding Joukowski airfoil. It’s not too hard to calculate \(K\)(we have a set of points from the user corresponding to \(T\), and we have our wholly analytic function \(J\)), and \(K\)’s derivative (via somewhat fancy finite difference methods ). With exact knowledge of \({dJ}\over{dc}\), and good knowledge of \({dK}\over{dc}\), we have a good estimate of \({dT}\over{dc}\), which with our exact knowledge of the flow about a cylinder via the complex potential (\({d\Phi}\over{dc}\)), we have a good estimate of the flow about our arbitrary airfoil:
\[{{d\Phi}\over{dT}} = {{{d\Phi}/{dc}}\over{dT/dc}} = {{{d\Phi}/dc}\over{{d\over{dc}}(J+K)}} = {{{d\Phi}/dc}\over{dJ/dc+dK/dc}}\]
This is the big picture. I am left with two smaller sub problems:

• Find an appropriate circle(minimize the magnitude of K)
• Find the derivative of K (the method described in the link above is still not all I'd hoped for)

The better my choice of a circle, the smaller my K function, and the less important the accuracy of K’s derivative. At the same time, the more accurate the derivative, the less important the choice of circle. Art remains.

Biplanes and Induced Drag

I recently got involved in a discussion of biplanes and induced drag. It reinforced the impression that I had that many (to include formally trained engineers) were perpetuating some misconceptions about induced drag, about the price you have to pay to get lift.
The biggest misconception is the idea that to minimize induced drag, you must maximize aspect ratio, you must favor long, skinny wings. This misconception is reinforced by the commonly cited formula:
\[C_{d_i} = \frac{C_l^2}{\pi e AR}\]
Yes, the larger the aspect ratio, the smaller the induced drag coefficient, but when you're dealing with non-dimensional coefficients, you're distancing yourself from the real world parameters that truly make or break a design.
\[\frac{D_i}{\frac{1}{2}\rho V^2 S} = \frac{(\frac{L}{\frac{1}{2}\rho V^2 S})^2}{\pi e \frac{b}{c}}\]
\[D_i=\frac{L^2}{\frac{1}{2}\rho V^2 bc \pi e \frac{b}{c}}\]
\[D_i=\frac{L^2}{\frac{1}{2}\pi e\rho V^2 b^2}\]
Induced drag is proportional to the square of lift (and therefore the square of weight), and inversely proportional to the square of the span.
\[D_i \varpropto \frac{W^2}{b^2}\]
The chord doesn't enter into it, aspect ratio doesn't matter! Of course, area matters to profile drag, which drives designers to narrow chords, which results in high aspect ratios. But the aspect ratio, per se, isn't the important bit!
So how does this relate to biplanes? Suppose we divide the load between two wings:
\[D_{i_{biplane}} \varpropto 2 \frac{(W/2)^2}{b^2} = \frac{2}{4}\frac{W^2}{b^2} \varpropto \frac{1}{2} D_i\]
That's right, twice the wings, half the induced drag (assuming you don't shrink the span). Some of those flying Venetian blinds from the early days of aviation were not the work of cranks.
A carefully designed biplane can offer an overall drag reduction, but it does tend to have a narrow performance envelope because of the need to balance profile and induced drag, and to carefully manage interference effects.
Examples:

• Rutan Quickie
• Flying Flea
• Early Flying Wing Biplane

A few weeks ago, I saw a page on recent drone work along these lines, but of course as I write this, I can't find it.

Symbols:

• \(C_{d_i}\): Coefficient of induced drag
• \(C_l\): Coefficient of lift
• \(\pi\): A mathematical constant, approximately 3.14159
• \(e\): An efficiency factor, 1.00 for elliptical lift distributions
• \(AR\): Aspect ratio, for a rectangular wing, the span divided by the chord
• \(D_i\): Induced drag
• \(L\): Lift
• \(\rho\): Density
• \(V\): Velocity
• \(S\): Wing area
• \(b\): Wing span
• \(c\): Wing chord
• \(W\): Weight