Stagnation Flow - Part IV

My last post led us up to the point where we had to integrate...
$$\begin{align} F_{YYY} &= (F_Y)^2 - FF_{YY} - 1\\ F(0) &= 0\\ F_Y(0) &= 0\\ \lim_{Y\to\infty}F_{Y} &= 1\end{align}$$

... and we had to pick \(F_{YY}(0)\) to enforce the last boundary condition. I now sketch a solution suitable for graphing, leaving the detailed solution for you to examine in ZFoil’s source code. The code here is illustrative of ZFoil’s techniques, but should not be taken as "good enough" in terms of numerical accuracy. First, we need an integration routine - this one is only second order, ZFoil’s is much better.

> (define (integrate h g0 gp gpp)     (+ g0 (* h gp)(* 0.5 (sqr h) gpp)))

Now, we apply this to our specific problem.

> (define (integrate-F h F FP FPP)     (define FPPP (- (sqr FP)(* F FPP) 1.0))     (values (integrate h F FP FPP)             (integrate h FP FPP FPPP)             (integrate h FPP FPPP 0.0)))

Next we define a thunk that returns values based on an initial guess at \(F_{YY}(0)\).

> (define (create-soln h FPP0)     (let[(Y  0.0)          (F  0.0)          (FP 0.0)          (FPP FPP0)]       (λ()         (set! Y (+ Y h))         (set!-values (F FP FPP)(integrate-F h F FP FPP))         (values Y F FP FPP))))

Next we try to bracket a good value for \(F_{YY}(0)\).

> (define infty 5.0)

> (define h 0.05)

> (define (FPP->FP FPP)     (let[(step (create-soln h FPP))]                 (let loop [(Y 0.0)(FP 0.0)]                   (if (< Y infty)                       (let-values([(Y F FP FPP)(step)])                         (loop Y FP))                       FP))))

> (define FPP- 1.2)

> (define FP- (FPP->FP FPP-))

> (printf "Too small: F'' = ~a yields F' = ~a~n" FPP- FP-)

Too small: F'' = 1.2 yields F' = 0.208671690207321

> (define FPP+ 1.3)

> (define FP+ (FPP->FP FPP+))

> (printf "Too big:  F'' = ~a yields F' = ~a~n" FPP+ FP+)

Too big:  F'' = 1.3 yields F' = 2.4006733797625706

FP- is too small, and FP+ is too big. We use interpolation to find a value that’s just right.

> (define FPP= (+ FPP- (/ (- 1.0 FP-)(/(- FP+ FP-)(- FPP+ FPP-)))))

> (printf "Just right? F'' = ~a~n" FPP=)

Just right? F'' = 1.2361007162340845

FP= should be just right. We will now see what the velocity profile of stagnation flow looks like by plotting FP versus Y.

> (let[(right-soln? (create-soln h FPP=))]     (plot (points (for/list [(i (in-range 100))]                     (let-values ([(Y F FP FPP)(right-soln?)])                       (vector FP Y))))))

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

Well, we’re close, and my choice of \(\infty\) seems justified. Let’s narrow the search a bit.

> (set! FPP- (- FPP= 0.01))

> (set! FP-  (FPP->FP FPP-))

> (set! FPP+ (+ FPP= 0.01))

> (set! FP+ (FPP->FP FPP+))

> (set! FPP= (+ FPP- (/ (- 1.0 FP-)(/(- FP+ FP-)(- FPP+ FPP-)))))

> (printf "Just right! F'' = ~a~n" FPP=)

Just right! F'' = 1.24051740187513

> (let[(right-soln (create-soln h FPP=))]     (plot (points (for/list [(i (in-range 100))]                     (let-values ([(Y F FP FPP)(right-soln)])                       (vector FP Y))))))

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

This solution seems to be asymptotic to F’=1, not to merely approach it and run away, and so we accept it as good enough (considering the step size). We now possess a single solution to incompressible viscous stagnation flow, that can be scaled to a particular situation with values for velocity gradient and kinematic viscosity. In the next installment, we explore how this solution to the full Navier-Stokes equations shows the way to simpler boundary layer equations.

Stagnation Flow - Part III

Most of this post is based on the relevant discussion in Frank White’s "Viscous Fluid Flow" (I’m not using the current edition, and I surely didn’t pay $240 for it back in the ’80s). As this material comes from work done in about 1914, I don’t think I’m violating any copyright laws!
A couple of posts ago I introduced the stream function for inviscid flow at a stagnation point, such as that near the leading edge of an airfoil...
$$\Psi=Bxy$$
...where B is a real constant. Using the results of my last post, we find that the velocity parallel to the surface is...
$$u=Bx$$
Intuitively, we expect friction to retard the flow near the surface, and to have less effect the further from the surface we get. We express this with the following boundary conditions...
$$\begin{align} \lim_{y\to 0} u &=0\\ \lim_{y\to\infty} u &=Bx\end{align}$$
We’d also like to keep the form of the inviscid solution as much as we can (engineers being at least as lazy as anyone else), and so the following modified stream function and boundary conditions are proposed...
$$\begin{align} \Psi&=Bxf(y)\\ u&=Bxf’(y)\\ f’(0)&=0\\ \lim_{y\to\infty}f’(y)&=1\end{align}$$
What the above sort of says (to me) is that the function f(y) is merely a distortion of y, very distorted near the surface, and indistinguishable from y far enough away from the surface.
Is this distortion justified? Yes, if we can prove that introducing it has not violated any laws of physics [conservation of mass (continuity) and conservation of momentum (Navier-Stokes)].
Well, since we’ve retained the stream function, we don’t have to worry about continuity. We do need to look at Navier-Stokes(NS), starting in the "y" direction, assuming steady incompressible flow (note that subscripts denote partial derivatives here; yes, I’m not at all consistent in this regard, and very sloppy about noting what convention I’m using at a given moment... but I’m not charging you $240, either)...
$$\begin{align} -p_y+\mu(v_{xx}+v_{yy})&=\rho(uv_x+vv_y)\\ &\text{and now substitute derivatives of }\Psi\\ &\text{as appropriate}\\ -p_y+\mu[((-Bxf)_x)_{xx}-((Bxf)_x)_{yy}]&=\rho[(Bxf)_y((-Bxf)_x)_x-(Bxf)_x((-Bxf)_x)_y]\\ -p_y+\mu[(-Bf)_{xx}-(Bf)_{yy}&=\rho[-(Bxf_y)(Bf)_x+(Bf)(Bf)_y]\\ -p_y+\mu(0-Bf_{yy})&=\rho B^2[-xf_y(0)+ff_y]\\ -p_y-\mu Bf_{yy}&=\rho B^2 ff_y\\ &\text{and as nothing remaining is a function of x}\\ -p_{xy}&= 0 \end{align}$$
... which implies that \(p_x\) is a function only of x. This lets us pick the y value where we calculate the pressure gradient across the surface, and we pick a value where we know the most, far from the surface where friction can be ignored. In such circumstances, Bernoulli rules...
$$\begin{align} p&=p_0-\frac{1}{2}\rho(u^2+v^2)\\ &(p_0\text{ represents total pressure})\\ p_x&=-\rho(uu_x+vv_x)\\ &\text{and again we substitute our }\Psi \text{ derivatives}\\ p_x&=-\rho[(Bxf_y)(Bf_y)+(-Bf)(-Bf_x)]\\ p_x&=-\rho B^2[x(f_y)^2+f\cdot 0]\\ &f_y \text{ approaches 1 far from the surface.}\\ p_x&=-\rho B^2 x \end{align}$$
Now we are ready to tackle NS in the "x" direction.
$$\begin{align} -p_x+\mu(u_{xx}+u_{yy})&=\rho(uu_x+vu_y)\\ \rho B^2 x+\mu[(Bxf_y)_{xx}+(Bxf_y)_{yy}]&=\rho[(Bxf_y)(Bxf_y)_x-(Bf)(Bxf_y)_y]\\ \rho B^2 x+\mu(0+Bxf_{yyy})&=\rho xB^2[(f_y)^2-ff_{yy}]\\ 1+\frac{\nu}{B}f_{yyy}&=(f_y)^2-ff_{yy}\end{align}$$
It is convenient to scale both y and f by the factor \(\sqrt{\frac{B}{\nu}}\) (a study of why this is convenient is worthwhile), yielding the non-dimensional equation and boundary conditions...
$$\begin{align} 1+F_{YYY}&=(F_Y)^2-FF_{YY}\\ F(0)&=F_Y(0)=0\\ \lim_{Y\to\infty} F_Y &= 1\end{align}$$
Next time, we use Racket to generate a numerical solution to the above, and to explore some of the consequences of that solution.

Stagnation Flow - Part II

In my last post, I mentioned that we’d be working with the stream function \(\Psi\), rather than the complex potential \(\Phi\), but I haven’t really defined it well. You’ll remember that \(\Phi\) is related to the velocity field as follows:
$${{d\Phi}\over{dz}}=\overline{V}=u-iv$$
where u and v are the velocities along the x and y axes, respectively. In plotting streamlines, I’ve noted that the imaginary part of \(\Phi\) is \(\Psi\).
$$\Phi =\phi + i\Psi$$
We shall ignore \(\phi\)(little "fee"), as its properties where friction is a factor are the reason we’re tossing \(\Phi\) aside in the first place! So how does \(\Psi\) relate to velocity?
$$\begin{align} u-iv & = {{d\Phi}\over{dz}}\\ & = {{d\Phi}\over{dx}}\\ & = {{d\phi}\over{dx}}+i{{d\Psi}\over{dx}}\\ -iv & = i{{d\Psi}\over{dx}}\\ v & = -{{d\Psi}\over{dx}}\end{align}$$
Above we let dz = dx. Now we let dz = idy.
$$\begin{align} u-iv & = {{d\Phi}\over{dz}}\\ & = {{d\Phi}\over{i dy}}\\ & = {{d\phi}\over{i dy}}+i{{d\Psi}\over{i dy}}\\ & = -i{{d\phi}\over{i dy}}+{{d\Psi}\over{dy}}\\ u & = {{d\Psi}\over{dy}}\end{align}$$
Now we have velocity in terms of the stream function, which will prove to be useful in my next posting, in which we plug the stream function into Navier-Stokes.

Stagnation Flow - Part I

Now that the inviscid portion of ZFoil is working (although I already have some tweaks I want to try), it’s time to work on boundary layer code. The appropriate place to begin is, of course, the beginning, and the boundary layer begins at the leading edge, at the stagnation point where the flow is brought to a halt. Let us build such a flow using some of the tools I’ve already introduced.
Let us start with a uniform flow:
$$\overline{V}={{1}\over{2}}B$$ $$\Phi={{1}\over{2}}Bz$$
B is a real constant, and the fraction will prove to be convenient. The flow is too simple to even bother plotting, but let us apply the following transform:
$$w^2=z$$ $$\Phi={{1}\over{2}}Bw^2$$ $$\overline{V}=Bw$$

> (define (Phi w) (* 1/2 B w w))

> (define B 1)

> (plot (contours          (λ(x y)            (imag-part (Phi (make-rectangular x y))))         -3.0 3.0 0.0 1.5         #:levels '(-4.0 -2.0 -1.0 -0.5 -0.25 0.25 0.5 1.0 2.0 4.0))         #:height 200)

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

In z space, B is a measure of velocity, but in w space, it is a measure of velocity gradient. The stagnation point is at w = 0. Note that I have plotted the imaginary component of the complex potential \(\Phi\) ("fee"), the stream function \(\Psi\) ("psi", very close to "sigh"). When I introduce friction in my next post, we will be able to retain \(\Psi\), but \(\Phi\) will have to go, as it assumes irrotational flow. I now present an explicit form of the inviscid stream function:
$$\begin{align} \Phi&={{1}\over{2}}Bw^2\\ \Phi&={{1}\over{2}}B(x+iy)^2\\ \Phi&={{1}\over{2}}B[(x^2-y^2)+i(2xy)]\\ \Psi&=Bxy\\ \end{align}$$
Yes, the streamlines are hyperbolas. In my next post, we will deal with how to modify this stream function to account for friction.

A Momentous Occasion

The pitching moment code is done (and the lift code is fixed; my integration code had a bug), and it seems to work! Watch:

> (require "../zfoil.rkt"            "../airfoils.rkt")

> (define my-geo (make-geometry naca-4412))

> (define alphas '(-0.05 0.0 0.05 0.1 0.15))

> (define CMs (for/list [(alpha alphas)]                 (CM alpha my-geo)))

> (plot (points (for/list [(alpha alphas)                            (CM CMs)]                   (vector alpha CM)))         #:x-label "alpha"         #:y-label "Cm"         #:y-max -0.1         #:y-min -0.14)

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

This agrees pretty well with other numerical methods, and more importantly, agrees with actual data available here.

Next, viscous flow and drag. This is to be strictly "cook book" for now. We’ll see if I am bored into motionless, or code it up quickly so I can actually play with my new toy!

Linux Migration Discovers Windows Bugs

Okay, they were my bugs, but Windows lets me get away from them, and Linux won't. If you've examined my code, you know that I have organized it pretty much per Racket convention, with a directory named "Private" trying to hide some of the messy details. Private, with a capital "P".
When I first tried to run my old code on my new Linux setup, I suddenly got file not found errors when clearly, the files could be found. But Racket was looking in "private", with a small "p".
At first, I attributed this to some weird Windows/Linux file system quirk, invoked when I cut and pasted my "Source" directory tree from my old Windows drive to my new Linux drive, and Windows two-tiered file name system was mapped onto Linux' straightforward one.
It was not a Windows/Linux quirk, but a Racket implementation choice and my hidden typo.
Under the hood running in Windows, Racket uses the short file names, the old xxxxxxxx.yyy format, which is case insensitive. I had a nested require statement that used "private" instead of "Private", and the error had been quietly ignored under Windows, but of course not in Linux.
The nested part was especially vexing, because I had required the same file in the module I was working on, without the typo. Thus the error seemed mysterious, and Racket's long, long, long error trace hid the nested attempts to require the same file.
I found the issue in Linux by using the grep command, and what a joy it is to have the Linux tools at my fingertips again (yes, Windows implementations exist, but I always found them a royal pain).
Anyway, since making the migration early this week, things are progressing nicely now that I'm not banging my head against the OS every 20 minutes with another BSOD. The code at GitHub is current, and even runnable! Errors in calculated velocities are on the order of 8% of free stream, and lift coefficients are within 0.1 of the correct ones (which sounds large, but that's only about a 1 degree error in angle of attack; I'm hoping that Cl vs Cm and Cl vs Cd will look pretty good).
Next, pitching moments.

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.

How Complex Numbers Make Things Simple

I need a robust test case that can validate both calculated velocities and lift coefficients. Luckily, a special case of the Joukowsky transformation on which zfoil is based provides an analytic solution of the flow over a flat plate. I’ll use this for my test case, as a means to convey to you the power of complex numbers, and as an excuse to play with Racket’s plot library.
I’ll be taking some baby steps here, both to refresh mysaelf on some of the fine points of complex functions and of fluid mechanics, and to make the workings of zfoil accessible to users whose comfort zone hasn’t pushed much beyond algebra.
Ok, what’s a complex number? It’s the sum of a real number, and an imaginary one. Imaginary numbers are multiples of \(\sqrt{-1}\), usually denoted by i or j. All the properties of complex numbers emerge from this definition, and many of the properties of real numbers extend to the complex.
\[z=x+iy\]
Why complex numbers? Well, from their very definition, they do something reals can’t – handle the square root of negative numbers. They can also handle a number of other situations where real numbers seem to faulter, either by increasing the space to look for solutions (from the number line to the complex plane), or by tweaking the definition of common real functions to accomodate complex numbers. As a result, we say that complex numbers are closed – any operation on a complex number (or its subset, the reals) yields a complex number.
That’s pretty cool, mathematically. Complex numbers clean up some of the messes of real numbers, and don’t create any new messes. But how they do this, introducing an "imaginary" second dimension, means that a complex number can represent a vector, and a complex function can map one vector into another. Since in fluid mechanics, we are often mapping positions to velocities, both vectors, wouldn’t it be convenient if we could use complex numbers to describe fluid flows?
Consider a generic complex function:
\[F(z)=u(x,y)+iv(x,y)\]
Very loosely speaking, if \(F(z)\) is well behaved, if you can write it much as you might write a real function, then \(F\) is said to be analytic, and its constituent real functions \(u\) and \(v\) satisfy the Cauchy-Riemann (CR) equations:
\[u_x - v_y = 0\] \[u_y + v_x = 0\]
The CR equations are easily derived by finding \({dF}\over{dz}\) via the traditional limit process, first by letting \(\Delta z=\Delta x\), and then by letting \(\Delta z=i\Delta y\), and then by supposing that the results ought to be identical.
Now consider a generic description of a velocity field:
\[V(z)=u(x,y)+iv(x,y)\]
Fluid flows are obviously subject to the laws of physics, such as conservation of mass, which for an incompressible fluid in steady motion can be expressed as:
\[u_x + v_y = 0\]
The above is a special case of the continuity equation, derived by converting a Lagrangian description of an infinitesimal fluid element’s mass to an Eulerian description of a density field – less scary than it sounds, maybe another article.
Conservation of momentum, aka Newton’s second law, also must be obeyed. In other words, you have to consider the forces acting on a fluid. The force of friction can be eliminated if:
\[u_y - v_x = 0\]
The above is a statement that the flow is "irrotational" (see fluid kinematics). If the above relation is used to generate second derivatives to substitute into Navier-Stokes (the complete fluid mechanics description of momentum), you will find that, if it holds, then friction is irrelevant.
So the physical laws that velocity must obey (given that the assumptions above hold, the flow is incompressible, and friction can be ignored) differ from the mathematical laws that analytic complex functions must obey only by the sign of the v term, the sign of the imaginary component. This means that the conjugate of the velocity is analytic!
\[\bar{V}=u(x,y)-v(x,y)\]
So any analytic complex function we can write represents the conjugate of a physically possible velocity field! That is very cool!
Let’s use Racket to plot a complex function, \(\bar{V}=1-{1\over{z^2}}\)(or rather, it’s conjugate, which will appeal to our sense of what a velocity field should look like).

> (require plot)

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

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

> (define (vbar z)     (- 1.0 (expt z -2.0)))

> (plot (vector-field          ; Cast vbar into a form suitable for plot          (λ(x y)(complex->vector                  (if (> (+ (sqr x)(sqr y)) 1.0)                      (conjugate                       (vbar (make-rectangular x y)))                      ; Suppress large velocities near origin                      0.0)))          -2.0 2.0 -2.0 2.0))

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

Hey, that almost looks like the flow about a cylinder! Guess what, it is the (incompressible inviscid) flow about a cylinder, simply because we can write it!
Earlier, I referred to Lagrangian & Eulerian descriptions. Most of you were exposed to the Lagrangian point of view in high school or college physics, where the "answer" was often the position of a particle as a function of time, and velocity and acceleration could be found by taking the derivative of that function, a derivative in time.
In fluid mechanics, we are not interested in the current position of a given particle, but the properties of whatever particle is at a point at a given time. This is the Eulerian viewpoint. When you investigate the mathematics that relate the two viewpoints, you find that the temporal derivative in the Lagrangian approach is expressed as a spatial derivative in the Eulerian.
So if the temporal integral of a particle’s velocity is its position, what do you get when you integrate the conjugate of velocity over space?
I’ll let Racket show you:

; Integrate vbar

> (define (Phi z)     (+ z (expt z -1.0)))

; Values of contour lines

> (define cv '(-1.5 -1.2 -0.9 -0.6 -0.3 -0.01 0.01 0.3 0.6 0.9 1.2 1.5))

; Plot imaginary component of Phi

> (plot (list          (contours           (λ(x y)             (let [(norm (+ (sqr x)(sqr y)))]               (if (> norm 1.0)                   (imag-part (Phi (make-rectangular x y)))                   ; Suppress interior flow in hopefully aesthetic way                   1.0)))           -2.0 2.0 -2.0 2.0           #:levels cv)))

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

The function \(\Phi\)(Phi) is called the complex potential, and its imaginary component is called the stream function. As you can see, contours of the stream function in the second plot correspond well with the direction of the velocity vectors in the first plot, and the cylindrical flow we’re investigating is clearly apparent.
Now the really nifty part – what happens if we take the derivative of \(\Phi\) with respect to a different spatial variable?
Let \(w=ze^{i\alpha}\), or \(z=we^{-i\alpha}\). "W" and "Z" space differ by a simple rotation. Our cylindrical flow shouldn’t change, as we’re just changing our viewpoint. As proper lazy programmers, we’ll phrase things to reuse as much code as possible by using composition and the chain rule:
\[\Phi_{1} = \Phi(z(w))\] \[{{d\Phi_1}\over{dw}}={{d\Phi}\over{dz}}{{dz}\over{dw}}\] \[{{d\Phi_1}\over{dw}}={\bar{V}(z(w))}{{dz}\over{dw}}\]

> (define alpha 0.3)

> (define (z w)(* w (exp (* 0-1i alpha))))

> (define dz/dw (exp (* 0-1i alpha)))

> (define (Phi1 w)     (Phi (z w)))

> (plot    (list (contours           (λ(x y)             (let [(norm (+ (sqr x)(sqr y)))]               (if (> norm 1.0)                   (imag-part (Phi1 (make-rectangular x y)))                   1.0)))           #:levels cv)          (vector-field           (λ(x y)(complex->vector                   (if (> (+ (sqr x)(sqr y)) 1.0)                       (conjugate                        (* dz/dw (vbar (z (make-rectangular x y)))))                       0.0)))))    #:x-min -2.0    #:x-max  2.0    #:y-min -2.0    #:y-max  2.0)

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

Notice that using composition with \(\Phi\) is somewhat easier than the chain rule with \(\bar{V}\). This tends to be true generally. It pays to manipulate \(\Phi\) and invoke its derivative only when neccessary.
Ok, in yet another display of foreknowledge, I’m going to add circulation to the picture.
\[\Phi_{2} = \Phi_{1} + {k}\cdot{ln(w)}\]
We want the strength of the circlation to be such that the point \(w = 1\) is a stagnation point.
\[{{d\Phi_{2}}\over{dw}}(1) = 0\] \[{{{d\Phi_1}\over{dw}}(1) + {k\over 1}}=0\] \[k = -{{d\Phi_1}\over{dw}}(1)\] \[k = -{{d\Phi}\over{dz}}(z(1))\cdot{{dz}\over{dw}}\] \[k = -{\bar{V}(z(1))}\cdot{e^{-i\alpha}}\] \[k = -(1-{1\over{z(1)^2}})\cdot{e^{-i\alpha}}\] \[k = -(1-{1\over{(e^{-i\alpha})^2}})\cdot{e^{-i\alpha}}\] \[k = -(e^{-i\alpha}-e^{-i\alpha}\cdot{e^{2i\alpha}})\] \[k = e^{i\alpha}-e^{-i\alpha}\] \[k = 2i sin \alpha\]

> (define (Phi2 w)     (+ (Phi1 w)        (* 0+2i           (sin alpha)           (log w))))

> (plot    (contours           (λ(x y)             (let [(norm (+ (sqr x)(sqr y)))]               (if (> norm 1.0)                   (imag-part (Phi2 (make-rectangular x y)))                   1.0)))           -2.0 2.0 -2.0 2.0           #:levels cv))

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

... which sets us up for my last great feat of foreknowledge. Let \(f = w + {1\over w}\), and \(\Phi_{3}=\Phi_{2}(w(f))\). Unlike the earlier simple rotation, the transformation between "F" and "W" space presents a little challenge. Note that as \(w \rightarrow \infty\), \(f \rightarrow w\). In other words, from far enough away, our flow field is unchanged by the transformation. However, in the near field, the flow is distorted in an interesting way. \[f = w + {1\over w}\] \[fw = w^2 +1\] \[0 = w^2 - fw + 1\] \[w ={{f \pm \sqrt{f^2 - 4}}\over{2}}\] Note that \(w(f)\) is double valued, but we are interested in the solutions lying outside our cylinder.

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

> (define (Phi3 f)     (Phi2 (w f)))

> (plot    (contours     (λ(x y)       (imag-part (Phi3 (make-rectangular x y))))     -3.0 3.0 -3.0 3.0     #:levels cv))

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

Yes, the flow about our cylinder is mapped to the flow about a flat plate (hence "F" space). Racket’s plotting ability again provides backup to my mathematical meandering, and I have an exact description of the flow around a flat plate to match zfoil’s calculations against. The only step left is to confine z to \(e^{i\theta}\), which will result in f being confined to the x axis between -2 and 2... my test flat plate. This results in a neat real-valued function for \(\bar{V}\) (and since we’re dealing with a flat plate, also for \(V\)). \[V ={{sin(\theta - \alpha) + sin \alpha }\over{sin \theta}}\] \[f = 2 cos \theta\] By relating f and V to \(\theta\), we explicitly select whether we’re dealing with the upper or lower surface of our flat plate. Note that near the trailing edge, \(V \rightarrow {{\theta - \alpha + \alpha}\over{\theta}} \rightarrow 1\), but at the leading edge, \(V \rightarrow {{{sin(\pi - \alpha)}+sin {\alpha}}\over{sin{\pi}}}\rightarrow {{2 sin{\alpha}}\over 0} \rightarrow \infty\). Infinity sucks. It can be dealt with, symbolically, and the result is exactly what you’d expect for lift over a flat plate, but no numerical method can compare to its symbolic certainty, and so I do not quite have my perfect test case yet.
Next post, we'll approach the performance of a flat plate with a very thin airfoil with a blunt leading edge, which will keep velocities finite everywhere.