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.
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.
No comments:
Post a Comment