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.
