Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Parameter sweep for Falkner-Skan flat plate flow
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 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 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 = 2.35e-2
beta_0 = 1.0
x = sp.array([alpha, eta_inf])
nbeta = 331
beta = sp.linspace(-0.1981, 0.1981, nbeta)
foo = sp.zeros((nbeta,2))
exa = sp.zeros(nbeta)
exe = sp.zeros(nbeta)
for i in xrange(nbeta):
args = (beta_0, beta[i])
foo[i] = fsolve(objective_func, x, args=args, xtol=1e-12)
if(i==0):
x[0] = 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
r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-14, atol=1e-14, method='adams', with_jacobian=True, order=12)
nt = 200
u1 = sp.zeros(nt)
u1[0] = 0.0
y0 = sp.array([0.0,0.0,foo[0,0]])
dt = 1.0 / float(nt)
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)
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]])
dt = 1.0 / float(nt)
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)
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]])
dt = 1.0 / float(nt)
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)
u3[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.xlabel(r'$\beta$')
plt.ylabel(r'$\alpha$')
plt.savefig("alpha.png")
plt.figure()
plt.plot(u1, foo[0,1]*sp.linspace(0,1,u1.shape[0]), label=r'$\beta=%g$'%beta[0])
plt.plot(u2, foo[int(nbeta/2),1]*sp.linspace(0,1,u1.shape[0]), label=r'$\beta=%g$'%beta[int(nbeta/2)])
plt.plot(u3, foo[nbeta-1,1]*sp.linspace(0,1,u1.shape[0]), label=r'$\beta=%g$'%beta[nbeta-1])
plt.legend(loc=0)
plt.xlabel(r"$f'$")
plt.ylabel(r'$\eta$')
plt.axis([0.0,1.0,0.0,10.0])
plt.savefig("fs_profiles.png")
#plt.show()
@jstults

This comment has been minimized.

Copy link
Owner Author

jstults commented Feb 5, 2012

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.