Quadrature for a Mononomial

In my last post, I mentioned that my next step was to write a lift and moment routine, and that I’d been having difficulties with the interface between Windows 7 and Racket’s Scribble. Sadly, all my time (okay, all the time I’ve devoted to programming) has been devoted to Windows vs Scribble. And the loser is...(drum roll, please) ME! One of the neat things about Scribble documents are their ability to invoke live Racket commands and show their results via the interaction command. Like a good citizen in this unfriendly online world, interaction is properly suspicious of what it perceives to be foreign code. Unfortunately, my use of submodules to implement unit testing set off Racket’s paranoia. Thanks to the friendly Racket community , I eventually found call-with-trusted-sandbox-configuration, and now I can go back to scribbling and programming, at least until I fall down the next rabbit hole.

Okay, enough "mea culpa" about why it’s taken me so long to write a simple integration routine. ZFoil will use a simple enhancement of the trapezoidal rule for all its integration needs. The trapezoidal rule basically assumes that a function f is linear between points a & b, and so can estimate the integral of f by using the formula for the area of a trapezoid:

Note that this still works perfectly well when f crosses the y axis.

Now ZFoil predicts velocity pretty well, as I’ve shown. What’s more, it can predict stagnation points (where velocity touches zero, basically near the leading edge, and at the trailing edge), pretty well. A lot of the numbers we care about depend not on velocity, but of some power of velocity. Pressure (and therefore lift and moment) depends on the square of velocity, and in laminar flow, the fifth power of velocity takes on importance. Consider the plot of \(x\), \(x^2\), and a naive linear fit of \(x^2\) below:

As you can see, there’s a big difference between the area below the lines \(x^2\) and the linear fit (the green horizontal line at the top of the plot). Since ZFoil naturally spits out the sign of the velocity, we should take advantage of that information!

So, we have a function \(v(x)\). We assume \(v\) is linear between \(x_0\) and \(x_1\). We are trying to find \[\int_{x_0}^{x_1} v^n dx\] Let’s introduce the linear transformation \(T(x)=t\), such that when \(x = x_0, t =0\), and when \(x = x_1, t =1\). \[\int_{x_0}^{x_1} v^n dx = \int_0^1 v^n \frac{dx}{dt} dt = (x1 - x0)\int_0^1 v^n dt \] That last integral is easily evaluated: \[\int_0^1 v^n dt = \frac{1}{n+1}(v_1^{n+1}-v_0^{n+1})\] ...which is easily translated into a Racket function:

(define (integrate v0 v1 n)

(/(-(expt v1 (add1 n))

(expt v0 (add1 n)))

(add1 n)))

or, if you C refugees prefer, to be just a bit faster,

(define(integrate v0 v1 n)

(let[(n+1(add1 n))]

(/(-(expt v1 n+1)

(expt v0 n+1))

n+1)))

To me, that’s another plus for Racket, its flexible, and therefore expressive, variable names.

Ok, I’ve laid out my next baby step. Now I just have to take it, and incorporate unit tests, and see if I’ve truly killed the blogging bug I mentioned at the beginning of this post.

It Works!

I've been having issues with Scriblogify and/or Windows 7, so  I can't present the pretty code of my last few posts, but I wanted to show you that I've been making real progress.
Here is my program's analysis of the NACA 4412 at 10 degrees angle of attack:

And here is some real data, extracted from work done at Baylor University:

Qualitatively, I'm on to something. The difference in peak pressure is easily explained -- I've made no attempt yet to compensate for boundary layer effects. And here's the kicker -- my algorithm's execution time scales linearly with increasing number of points -- twice the accuracy costs only twice the computer time!
The source currently posted at https://github.com/SlowThought/SlowFlight is up to date, and sufficient to reproduce the first figure above, with an installation of Racket. I need to write some documentation, a lift & moment routine, and a boundary layer routine, before I can justify pushing a stand-alone executable, but what I've got is what I aimed for so long ago -- a method that doesn't require vortex panels, that doesn't require singularities within the solution space -- solutions where you're allowed to look anywhere within the space!

Testing a Good Error Function - Part II

In my last post, I wrote "Although I am writing a more objective test, this graph was sufficient to give me confidence in my code." Silly me. I shouldn’t have published before developing my objective test. I was right about the usefulness of graphs, however.

(define my-foil(joukowsky 0.25 -0.02+0.01i 101))

generates a perfectly reasonable airfoil...

...but when I analyze it...

(define my-geometry(make-geometry my-foil))

(display(plot-foil(geometry-C my-geometry)))

... the airfoil  fails to map to a circle. The problem seems to lie in choosing the correct root of \(J^{-1}\), but in fact it’s a little more subtle than that. The Joukowsky maps the exterior of a circle \(|c|>a\) onto the entire \(n\) plane, but it also maps the interior onto the same plane(\(|c|<a\)). When trying to perform the inverse mapping, my airfoil wanders across the boundary of \(a\), weird things happen.
So, back to my old approach, picking a point on \(|c|=a\), and seeing if it is close enough to a point \(n\). The heavy lifting on this is done, in that one of the subroutines to my Powell optimizer is pretty good at picking such points, and in that the code that let me find this bug is already incorporated into my test base.
It’s interesting that the size of my code base is shrinking (I get to strip out all the \(J^{-1}\) stuff), but my test base keeps growing. The older I get, the more I value testing discipline.

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

Thanks, Ryan!

Ryan Culpepper has developed a nifty tool for posting Scribble documents to blogs. Not only does this simplify my life as a writer, it allows you, the reader, to see Racket as I see it. Consider my last post’s example, as rendered by Ryan’s tool:

> (require"../private/geometry.rkt")

> (define my-geometry(make-geometry #(1 0.5+0.05i 0+0.1i 0.5+0.05i 1)))

Warning: (j-1 n a) is not yet implemented! Warning: (optimize-c0 c0 n a) is not yet implemented! Warning: (dt/dc ? ? ?) is not yet implemented!

> (geometry? my-geometry)

#t

> (geometry-chord my-geometry)

1.004987562112089

> (geometry-N my-geometry)

'#(0.5024937810560445+6.938893903907228e-018i    0.0+0.0i    -0.5024937810560445-6.938893903907228e-018i    0.0+0.0i    0.5024937810560445+6.938893903907228e-018i)

First up, see how pretty it is? This depiction is much closer to what I see in the Racket environment, although it still doesn't show the syntax coloring. At least I achieve it by writing Racket, writing in my favorite language. All those parentheses aren’t so bad, when you consider HTML’s tags are just the same but <tag>longer</tag>.
And second, my first bug is fixed, as the result of geometry-N is all real (or nearly so, I can accept imaginary components on the order of 1e-18), as it should be. Future topic: incorporation of testing into the development process, or how to document past mistakes as future requirements.

First Code

It's just a sketch of the geometry analysis portion, really, but I've demonstrated that I've got my programming tools reestablished and github working.
Per Racket practice, the ugly code goes in the "Private" directory, and the pretty public face, including contracts enforcing usage and hiding some implementation, will be in the root directory.
Here's a sample interaction, demonstrating the syntax for newcomers to things Lispish:

Welcome to DrRacket, version 5.2.1 [3m].
Language: racket/base; memory limit: 128 MB.
Warning: (test-make-geometry) is not yet implemented!
> (define my-geometry(make-geometry #(1 0.5+0.05i 0+0.1i 0.5+0.05i 1)))
Warning: (j-1 n a) is not yet implemented!
Warning: (optimize-c0 c0 n a) is not yet implemented!
Warning: (dt/dc ? ? ?) is not yet implemented!
> (geometry? my-geometry)
#t
> (geometry-chord my-geometry)
1.004987562112089
> (geometry-N my-geometry)
'#(0.4925434091539446-0.09950371902099892i 0.0+0.0i -0.4925434091539446+0.09950371902099892i 0.0+0.0i 0.4925434091539446-0.09950371902099892i)
> 

We feed make-geometry a vector describing a flat plate set at an angle to the x-axis. We get expected warnings because this is a work in progress. We explore the returned structure. We discover our first bug! geometry-N is supposed to return a vector rotated so as to be parallel to the x axis. I apparently rotated the wrong way.

The Long and Winding Road

The little project I intend to document in this blog has been moving along in fits and starts since my undergraduate days, some (shudder) thirty years ago. At that time I was introduced to panel methods as a means of analyzing airfoils. Modeling an elegantly curved surface as a series of straight lines, and then modelling smooth airflow as a series of tiny tornadoes just struck me as ugly!
On the other hand, I was enchanted with the use of complex numbers to transform simple shapes, easy to analyze, into complex ones. Complex transforms have long been used in the design of airfoils (mapping a cylinder to a shape to achieve a certain velocity distribution), but not in the analysis of airfoils (mapping an airfoil to a circle in order to determine the velocity distribution).
I have gone down many blind alleys over the years, but started making real progress about 3 years ago. My math was solid, my code was solid, and then my laptop crashed.
I had intended to reorganize, anyway. Since I've given up becoming a rich aviation entrepreneur, I've decided to do my coding (https://github.com/SlowThought/SlowFlight) and thinking (here) in public, in the hope that somebody might find my musings useful.