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.