Non-uniqueness of Falkner-Skan solutions with adverse pressure gradient

integrate_fs2.py

import scipy as sp
from scipy.integrate import ode 
from scipy.optimize import fsolve 
from matplotlib import rc
rc('text', usetex=True)
tick_size = 'large'
label_size = 'xx-large'
rc('xtick', labelsize=8)
rc('ytick', labelsize=8)
rc('legend', fontsize=12)
rc('axes', labelsize=12)
fig_width_pt = 469.75502 
inches_per_pt = 1.0/72.27
golden_mean = (sp.sqrt(5.0)-1.0)/2.0
fig_width = fig_width_pt *inches_per_pt 
fig_height = fig_width * golden_mean
fig_size = [fig_width, fig_height] 
rc('figure', figsize=fig_size)

from matplotlib import pylab as plt

from fs import *

def int_fs(beta_0, beta, alpha, eta_inf):
    args = sp.array([beta_0, beta, eta_inf])    
    y0 = sp.array([0.0, 0.0, alpha])
    r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-14, atol=1e-14, method='adams', with_jacobian=True, order=12)
    t0 = 0.0 
    t1 = 1.0 
    nt = 200
    dt = t1 / float(nt)
    r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
    for i in xrange(nt):
        r.integrate(r.t+dt)
    return(r)

def int_tanh_fs(beta_0, beta, alpha, eta_inf):
    args = sp.array([beta_0, beta, eta_inf])    
    y0 = sp.array([0.0, 0.0, alpha])
    r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-14, atol=1e-14, method='adams', with_jacobian=True, order=12)
    t0 = 0.0 
    t1 = 1.0 
    # instead of using a uniform point distribution, use a tanh
    # mapping, which clusters the points near the wall where the
    # variation is largest
    nt = 64  
    t = sp.arctanh(sp.linspace(t0,t1,nt,endpoint=False)) 
    t = t/t.max() 
    dt = sp.diff(t) 
    r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
    for i in xrange(nt-1):
        r.integrate(r.t+dt[i])
    return(r)

def objective_func(x, beta_0, beta):
    alpha = x[0] 
    eta_inf = x[1]
    r = int_fs(beta_0, beta, alpha, eta_inf)
    return(sp.array([(r.y[1]-1.0).real, (r.y[2]).real]))

def objective_func2(x, beta_0, beta):
    alpha = x[0] 
    eta_inf = x[1]
    r = int_tanh_fs(beta_0, beta, alpha, eta_inf)
    return(sp.array([(r.y[1]-1.0).real, (r.y[2]).real]))

def extrap(y):
    return(2*y[1]-y[0])
def extrap2(y):
    return(-(4*y[1]-y[0]-5*y[2])/2.0)

eta_inf = 5.0  
alpha = 1e-2
beta_0 = 1.0 
x = sp.array([alpha, eta_inf])
nbeta = 71  
beta = sp.linspace(-0.199010445, 0.199010445, nbeta)  
foo = sp.zeros((nbeta,2))
for i in xrange(nbeta):
    args = (beta_0, beta[i]) 
    foo[i] = fsolve(objective_func2, x, args=args, xtol=1e-12)
    if(i==0):
        x[0] = alpha # foo[i,0]
    elif(i==1): 
        x[0] = extrap(foo[i-1:i+1,0]) 
    else:
        x[0] = extrap2(foo[i-2:i+1,0]) 
    x[1] = eta_inf 
goo = sp.zeros((nbeta,2))
alpha2 = -3e-1
eta_inf2 = 6.0  
for i in xrange(nbeta):
    args = (beta_0, beta[i]) 
    goo[i] = fsolve(objective_func2, x, args=args, xtol=1e-12)
    if(i==0):
        x[0] = alpha2 # foo[i,0]
        x[1] = eta_inf2 
    elif(i==1): 
        x[0] = extrap(goo[i-1:i+1,0]) 
        x[1] = eta_inf2 
    else:
        x[0] = extrap2(goo[i-2:i+1,0]) 
        x[1] = goo[i,1] 
#    x[1] = eta_inf2 

sp.savetxt("foo.txt", foo)
sp.savetxt("beta.txt", beta)

r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-14, atol=1e-14, method='adams', with_jacobian=True, order=12)
nt = 64 
u1 = sp.zeros(nt)
u1[0] = 0.0 
y0 = sp.array([0.0,0.0,foo[0,0]])
t0 = 0.0
t1 = 1.0 
t = sp.arctanh(sp.linspace(t0,t1,nt,endpoint=False))
dt = sp.diff(t) 
args = sp.array([beta_0, beta[0], foo[0,1]])    
r.set_initial_value(y0, 0.0).set_f_params(args).set_jac_params(args)
for i in xrange(nt-1):
    r.integrate(r.t+dt[i])
    u1[i+1] = r.y[1].real  
u2 = sp.zeros(nt)
u2[0] = 0.0 
y0 = sp.array([0.0,0.0,foo[int(nbeta/2),0]])
args = sp.array([beta_0, beta[int(nbeta/2)], foo[int(nbeta/2),1]])    
r.set_initial_value(y0, 0.0).set_f_params(args).set_jac_params(args)
for i in xrange(nt-1):
    r.integrate(r.t+dt[i])
    u2[i+1] = r.y[1].real  
u3 = sp.zeros(nt)
u3[0] = 0.0 
y0 = sp.array([0.0,0.0,foo[nbeta-1,0]])
args = sp.array([beta_0, beta[nbeta-1], foo[nbeta-1,1]])    
r.set_initial_value(y0, 0.0).set_f_params(args).set_jac_params(args)
for i in xrange(nt-1):
    r.integrate(r.t+dt[i])
    u3[i+1] = r.y[1].real  
ni = 14
u4 = sp.zeros(nt)
u4[0] = 0.0 
y0 = sp.array([0.0,0.0,goo[ni,0]])
args = sp.array([beta_0, beta[ni], goo[ni,1]])    
r.set_initial_value(y0, 0.0).set_f_params(args).set_jac_params(args)
for i in xrange(nt-1):
    r.integrate(r.t+dt[i])
    u4[i+1] = r.y[1].real  
u5 = sp.zeros(nt)
u5[0] = 0.0 
y0 = sp.array([0.0,0.0,foo[ni,0]])
args = sp.array([beta_0, beta[ni], foo[ni,1]])    
r.set_initial_value(y0, 0.0).set_f_params(args).set_jac_params(args)
for i in xrange(nt-1):
    r.integrate(r.t+dt[i])
    u5[i+1] = r.y[1].real  

plt.figure()
plt.plot(beta, foo[:,1]) 
plt.xlabel(r'$\beta$')
plt.ylabel(r'$\eta_{\infty}$')
plt.savefig("eta_inf.png")
plt.figure()
plt.plot(beta, foo[:,0])
plt.plot(beta, goo[:,0]) 
plt.xlabel(r'$\beta$')
plt.ylabel(r'$\alpha$')
plt.savefig("alpha.png")
plt.figure()
plt.plot(u1, foo[0,1]*t, label=r'$\beta=%g$'%beta[0])
plt.plot(u2, foo[int(nbeta/2),1]*t, label=r'$\beta=%g$'%beta[int(nbeta/2)])
plt.plot(u3, foo[nbeta-1,1]*t, label=r'$\beta=%g$'%beta[nbeta-1])
plt.plot(u4, goo[10,1]*t, label=r'$\beta=%g$'%beta[10])
plt.plot(u5, foo[10,1]*t, label=r'$\beta=%g$'%beta[10])
plt.legend(loc=0)
plt.xlabel(r"$f'$")
plt.ylabel(r'$\eta$')
xmin,xmax,ymin,ymax = plt.axis() 
plt.axis([xmin,xmax,0.0,8.0])
plt.savefig("fs_profiles.png")
plt.show()

non-unique alpha vs beta

Laminar separation (negative wall shear) with adverse pressure gradient:
non-unique velocity profiles