Non-uniqueness of Falkner-Skan solutions with adverse pressure gradient
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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
non-unique alpha vs beta
Laminar separation (negative wall shear) with adverse pressure gradient:
non-unique velocity profiles