How Complex Numbers Make Things Simple - Part II

How Complex Numbers Make Things Simple - Part II (in which we complicate things just a teensy bit).
Last time, we computed the inviscid, incompressible flow over a flat plat, arriving at the surprisingly simple set of relations...
\[V ={{sin(\theta - \alpha) + sin \alpha }\over{sin \theta}}\] \[f = 2 cos \theta\]
... where \(\theta\)(theta) ranges from 0 to \(2\pi\)(a full circle in radians), \(\alpha\) is the angle of attack of the flat plate, \(V\) is the velocity on the flat plate (with \(V_\infty = 1\)), and \(f\) is the position on the flat plate. Values of \(\theta > \pi\) represent the lower surface of the flat plate.
These functions (correctly) predict that the velocity at the leading edge will be infinite. Correct, given the assumptions, but a problem when dealing with numerical integration schemes such as I’m using to predict lift coefficients.
Let’s back up just a little bit – the expression for \(f\) above is a combination and simplification of the more general relations...
\[f = w + {1\over w}\] \[w = e^{i\theta}\]
The definition for f is an instance of the Joukowski transformation and the definition for w describes the unit circle.
We are going to use Racket to investigate how a slightly different circle looks under the Joukowski transform.

> (require plot)

; Racket's plot library deals in 2D vectors, not complex

> (define (complex->vector z)(vector (real-part z)(imag-part z)))

; Unit circle

> (define (w1 theta)     (exp (* 0+1i theta)))

; Slightly offset & expanded circle

> (define (w2 theta)     (- (* 1.01 (exp (* 0+1i theta)))        0.01))

; Warping space, what geek doesn't want to do that?

> (define (f w)(+ w (/ w)))

; Plot the transformation of our slightly different circles

> (plot (list          (parametric (λ (theta)                        (complex->vector (f (w1 theta))))                      0 (* 2 pi)                      #:color "red")          (parametric (λ (theta)                        (complex->vector (f (w2 theta))))                      0 (* 2 pi)                      #:color "blue"))         #:y-min -0.5         #:y-max  0.5)

(object:2d-plot-snip% ...)

I have intentionally distorted the image, quadrupling the apparent thickness (note the limits on the y scale), to make clear that the new curve is a typical airfoil shape. Even at this scale, it’s obvious that it’s a very thin, symmetrical airfoil, not far from the flat plate we started with, but the blunt leading edge is enough to prevent infinite velocities.
The curve generated by the second circle is called a Joukowski airfoil. It’s not a good airfoil, but it is easily analyzed. We will now construct a complex potential function that describes the flow about the second circle.
\[\Phi =w\:e^{-i\alpha} + {{{1.01^2}e^{i\alpha}} \over {w+0.01}} + 2.02i\:sin(\alpha)\: ln(w+0.01)\] \[{d\Phi\over dw}=e^{-i\alpha}- {{1.01^2}e^{i\alpha} \over (w+0.01)^2}+{2.02i\:sin(\alpha)\over {w+0.01}}\]
This is a slight tweak of the function I created in my last posting. We will first plot in w space to confirm that I got the circle I intended, and then plot in f space to see how close our airfoil is to a flat plate.

> (define alpha 0.2)

> (define (Phi w)     (+ (* (exp (* 0-1i alpha))w)        (* (sqr 1.01)(exp (* 0+1i alpha))(/(+ w 0.01)))        (* 0+2.02i (sin alpha)(log (+ w 0.01)))))

> (plot    (contours     (λ(x y)       (imag-part (Phi (make-rectangular x y))))     -3.0 3.0 -3.0 3.0     #:levels '(-2.5 -2.0 -1.5 -1.0 -0.5 -0.05 0 0.05 0.5 1.0 1.5 2.0 2.5)))

(object:2d-plot-snip% ...)

> (define (w f)     (let* [(radical (sqrt (- (sqr f) 4)))            (w1 (+ f radical))]       (if (>= (magnitude w1) 2.0)           (/ w1 2.0)           (/ (- f radical) 2.0))))

> (plot    (contours     (λ(x y)       (imag-part (Phi (w(make-rectangular x y)))))     -3.0 3.0 -3.0 3.0     #:levels '(-2.5 -2.0 -1.5 -1.0 -0.5 -0.05 0.0 0.05 0.5 1.0 1.5 2.0 2.5)))

(object:2d-plot-snip% ...)

Great, it looks as if I’ve got things right. Now, did I fix the pesky problem of infinite velocity at the leading edge?

> (define (dPhi/dw w)     (+ (exp (* 0-1i alpha))        (* -1 (sqr 1.01)(exp (* 0+1i alpha))(/(sqr(+ w 0.01))))        (* 0+2.02i (sin alpha)(/(+ w 0.01)))))

> (define (df/dw w)     (- 1 (/(sqr w))))

> (plot    (parametric (λ(theta)                  (let[(w(w2 theta))]                    (vector (real-part (f w))                            (magnitude (/(dPhi/dw w)(df/dw w))))))                0 (* 2 pi)                #:x-min -3                #:x-max  3                #:y-min -0.1                #:y-max 4.9))

(object:2d-plot-snip% ...)

Perfect! A big, but finite, velocity spike at the leading edge. Now I’ve got something solid to test zfoil against.