Numerical solution of Falkner-Skan equation¶

In [44]:
import numpy               as np
import matplotlib.pyplot   as plt
from   scipy.integrate import solve_bvp
from   IPython.display import display, Markdown

Parameters¶

In [45]:
infty = 10.  
beta  = 0.
m     = 101  # initial resolution
tol   = 1e-4

Residual of the equivalent system of 1st-order ODEs¶

In [46]:
def fun(x,y):
    f   = y[0,:]
    u   = y[1,:]
    tau = y[2,:]

    taup = -f*tau - beta*(1-u**2)

    yp = np.array([u, tau, taup])
    return yp

Boundary conditions¶

In [47]:
bc = lambda ya,yb: np.array([ya[0], ya[1], yb[1]-1])

Initial guess¶

In [48]:
eta = np.linspace(0,infty,m)

tau = (1        - eta   /2 ) * (eta<=2)
u   = (eta      - eta**2/4 ) * (eta<=2) + 1.  * (eta>2)
f   = (eta**2/2 - eta**3/12) * (eta<=2) + 4/3 * (eta>2)

y   = np.array([f,u,tau])

Numerical solution¶

In [49]:
res = solve_bvp(fun, bc, eta, y, tol=tol, verbose=2)

   Iteration    Max residual  Max BC residual  Total nodes    Nodes added  
       1          2.65e-06       1.57e-21          101             0       
Solved in 1 iterations, number of nodes 101. 
Maximum relative residual: 2.65e-06 
Maximum boundary residual: 1.57e-21

Plot velocity profile¶

In [50]:
fig,ax = plt.subplots()

# Interpolation
etamax  = 4
Nplot   = 200
etaplot = np.linspace(0,etamax,Nplot)
yplot   = res.sol(etaplot)
uplot   = yplot[1]
ax.plot(uplot, etaplot, '-')

# Data points
ax.plot(res.y[1], res.x, 'o')

# Labels and fig props
ax.set_xlabel(r'$f^\prime$')
ax.set_ylabel(r'$\eta$')
ax.set_xlim(0,1)
ax.set_ylim(0,etamax)
Out[50]:
(0.0, 4.0)
No description has been provided for this image

Properties¶

In [51]:
display(Markdown(r'$f_w^{\prime\prime\prime} = %g $'%(res.y[2,0])))
display(Markdown(r'$\beta_1 = %g $'%(res.x[-1]-res.y[0,-1])))

$f_w^{\prime\prime\prime} = 0.4696 $

$\beta_1 = 1.21678 $

For $\beta = 0$, compare the result to Table 6.1 in Schlichting & Gersten (2017), p. 158.

In [ ]: