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.