In [243]:
import numpy
from matplotlib import pyplot
%matplotlib inline
N_resolution=[32,64,128]
In [244]:
from VortexPanel import Panel,solve_gamma,solve_gamma_kutta,plot_flow,make_jukowski,make_circle
In [245]:
def c_p(gamma): return 1-gamma**2
In [246]:
def C_P(theta): return 1-4*(numpy.sin(theta))**2
In [45]:
pyplot.figure(figsize=(8,6))
for m in range(3):
circle = make_circle(N_resolution[m])
solve_gamma(circle)
coeff = numpy.zeros(N_resolution[m])
for i,p_i in enumerate(circle): coeff[i]=c_p(p_i.gamma)
pyplot.plot(numpy.linspace(0,N_resolution[m],N_resolution[m])/N_resolution[m]*numpy.pi*2,coeff)
pyplot.plot(numpy.linspace(0,numpy.pi,360),C_P(numpy.linspace(0,numpy.pi,360)))
pyplot.legend(["N=32","N=64","N=128","Analytical"],loc=4)
pyplot.title('$C_p$ of a circular cylinder under different resolutions',size=16)
pyplot.xlabel(r'$\theta /rad$',size=20)
pyplot.ylabel('$C_p$',size=20)
Out[45]:
In [110]:
from BoundaryLayer import ddx_delta, heun, g_pohl, g_1, df_0, march
In [111]:
nu=1e-5
In [112]:
def half_c_f(lam, nu, delta, u_e):
Re_d = delta*u_e/nu
return df_0(lam)/Re_d
In [113]:
def taw_w(lam, nu, delta, u_e):
if u_e == 0:
return 0
else:
return half_c_f(lam, nu, delta, u_e)*u_e**2
In [114]:
def C_F_calc(tau, sx, N):
#return numpy.sum(tau[:iSep+1]*sx*numpy.pi/(N-1))
return numpy.trapz(tau[:iSep+1]*sx, dx = numpy.pi/(N-1))
In [257]:
In [259]:
for N in N_resolution:
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta, lam, iSep = march(s,u_e,du_e,nu) # solve!
taw = numpy.full_like(delta, 0)
taw = [taw_w(lam[i], nu, delta[i], u_e[i]) for i in range(N)]
sx = numpy.sin(s[0:iSep+1])
print ('When N = ' + '%i' %N)
C_F_circle = 2*C_F_calc(taw, sx, N)/numpy.pi
print ('Circle frictional coefficient = ' + '%0.2e' %C_F_circle)
C_F_flat = 1.33 * numpy.sqrt(nu/numpy.pi)
print("Flate plate: "+'%0.2e' %C_F_flat)
| Resolution | Circle | Flat plat |
|---|---|---|
| 32 | 4.04e-03 | 2.37e-3 |
| 64 | 4.06e-03 | 2.37e-3 |
| 128 | 4.07e-03 | 2.37e-3 |
In [40]:
pyplot.figure(figsize=(8,6))
pyplot.plot(s,taw*s_x)
pyplot.scatter(s[iSep+1], taw[iSep+1]*s_x[iSep+1],s=50,c="r")
pyplot.xlabel('$s$',size=20)
pyplot.ylabel(r'$\tau_w s_x$', size=20)
Out[40]:
In [207]:
def C_P_sep(s,iSep):
c1=C_P(s[0:iSep+1])
c2=C_P(s[iSep])*numpy.ones(s.size-iSep-1)
return numpy.concatenate((c1,c2),axis=0)
In [254]:
N=128
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta, lam, iSep = march(s,u_e,du_e,nu) # solve!
C_P_s=C_P_sep(s,iSep)
In [255]:
pyplot.figure(figsize=(8,6))
pyplot.plot(s,C_P_sep(s,iSep))
pyplot.xlabel('$s$',size=20)
pyplot.ylabel('$C_p$', size=20)
pyplot.scatter(s[iSep+1],C_P_s[iSep+1],s=50,c="r")
Out[255]:
In [256]:
def C_D(N,sy,C_P_s):
return numpy.trapz(C_P_s*sy, dx=numpy.pi/(N-1))
In [261]:
for N in N_resolution:
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
s_y=numpy.cos(s)
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta, lam, iSep = march(s,u_e,du_e,nu) # solve!
C_P_s=C_P_sep(s,iSep)
print('When N='+'%i'%N+', Drag coefficient= '+'%0.2e'%C_D(N,s_y,C_P_s))
In [264]:
def lift(panels):
c = panels[0].x[0]-panels[len(panels)/2].x[0] # length scale
return -4./c*numpy.sum([p.gamma*p.S for p in panels])
In [265]:
#calculate C_L in a range of angle of attack
def C_L_AA(foil,alpha_max,M):
C_L=numpy.zeros(M)
for i in range(M):
alpha=alpha_max*i*numpy.pi/180./(M-1)
solve_gamma_kutta(foil,alpha) # solve for gamma
C_L[i]=lift(foil) # print the lift
return C_L
In [266]:
def CLA(rr,alpha): return 2.*numpy.pi*(1+4./numpy.sqrt(3)/3.*rr)*numpy.sin(alpha)
In [267]:
#init
N=[32,64,128]
M=100
alpha_max=10.
#plot
pyplot.figure(figsize=(8,6))
for p in range(3):
Na=N[p]
foil1 = make_jukowski(Na,dx=0.15,dy=0,dr=0) # make foil
coff1=C_L_AA(foil1,10,M)
pyplot.plot(numpy.linspace(1,M,M)/M*alpha_max,coff1)
M=101
rr=0.15
coff2=numpy.zeros(M)
for q in range(M):
alpha=alpha_max*q*numpy.pi/180./(M-1)
coff2[q]=CLA(rr,alpha)
pyplot.plot(numpy.linspace(1,M,M)/M*alpha_max,coff2)
pyplot.xlabel('Angle of attack /deg')
pyplot.ylabel('$C_L$')
pyplot.legend(["N=32","N=64","N=128","Analytical"])
Out[267]:
In [268]:
# Pohlhausen Boundary Layer class
class Pohlhausen:
def __init__(self,panels,nu):
self.u_e = [abs(p.gamma) for p in panels] # tangential velocity
self.s = numpy.empty_like(self.u_e) # initialize distance array
self.s[0] = panels[0].S
for i in range(len(self.s)-1): # fill distance array
self.s[i+1] = self.s[i]+panels[i].S+panels[i+1].S
ds = numpy.gradient(self.s)
self.du_e = numpy.gradient(self.u_e,ds) # compute velocity gradient
self.nu = nu # kinematic viscosity
self.xc = [p.xc for p in panels] # x and ...
self.yc = [p.yc for p in panels] # y locations
def march(self):
# march down the boundary layer until separation
from BoundaryLayer import march
self.delta,self.lam,self.iSep = march(self.s,self.u_e,self.du_e,self.nu)
# interpolate values at the separation point
def sep_interp(y): return numpy.interp( # interpolate function
12,-self.lam[self.iSep:self.iSep+2],y[self.iSep:self.iSep+2])
self.s_sep = sep_interp(self.s)
self.u_e_sep = sep_interp(self.u_e)
self.x_sep = sep_interp(self.xc)
self.y_sep = sep_interp(self.yc)
self.delta_sep = sep_interp(self.delta)
In [269]:
# split panels into two sections based on the flow velocity
def split_panels(panels):
# positive velocity defines `top` BL
top = [p for p in panels if p.gamma<=0]
# negative defines the `bottom`
bottom = [p for p in panels if p.gamma>=0]
# reverse array so panel[0] is stagnation
bottom = bottom[::-1]
return top,bottom
In [270]:
# plot panels with labels
def plot_segment(panels):
pyplot.figure(figsize=(10,2))
pyplot.axis([-1.2,1.2,-.3,.3])
for i,p_i in enumerate(panels):
p_i.plot()
In [271]:
def solve_boundary_layers(panels,alpha=0,nu=1e-5):
# split the panels
top_panels,bottom_panels = split_panels(panels)
# Set up and solve the top boundary layer
top = Pohlhausen(top_panels,nu)
top.march()
# Set up and solve the bottom boundary layer
bottom = Pohlhausen(bottom_panels,nu)
bottom.march()
return top,bottom
In [272]:
def predict_jukowski_separation(t_c,alpha=0,N=128):
# set dx to gets the correct t/c
foil = make_jukowski(N,dx=t_c-0.019)
# find and print t/c
x0 = foil[N/2].xc
c = foil[0].xc-x0
t = 2.*numpy.max([p.yc for p in foil])
# solve potential flow and boundary layer evolution
solve_gamma_kutta(foil,alpha)
top,bottom = solve_boundary_layers(foil,alpha)
return top.x_sep, top.y_sep, top.iSep,bottom.x_sep,bottom.y_sep, bottom.iSep
In [278]:
N=128
M=101
rr=0.15
tx = numpy.zeros(M)
ty = numpy.zeros(M)
bx = numpy.zeros(M)
by = numpy.zeros(M)
ti = numpy.zeros(M)
bi = numpy.zeros(M)
for q in range(M):
alpha=alpha_max*q*numpy.pi/180./(M-1)
tx[q], ty[q], ti[q], bx[q], by[q], bi[q]=predict_jukowski_separation(0.15,alpha)
In [280]:
pyplot.figure(figsize=(8,6))
pyplot.plot(numpy.linspace(0,numpy.pi,M),bi)
pyplot.plot(numpy.linspace(0,numpy.pi,M),ti)
pyplot.xlabel('$s$',size=20)
pyplot.ylabel('$C_p$', size=20)
Out[280]:
In [274]:
foil = make_jukowski(N,dx=0.15,dy=0,dr=0) # make foil
top,bottom = split_panels(foil)
In [276]:
plot_segment(top)
In [275]:
plot_segment(top)
for i in range(M):
if i%20 == 0:
pyplot.scatter(tx[i], ty[i], s=100, c='r')
In [234]:
# define the influence of panel_j on panel_i
def influence(panel_i,panel_j):
u,v = panel_j.velocity(panel_i.xc,panel_i.yc,gamma=1)
return u*panel_i.sx+v*panel_i.sy
In [235]:
# construct the linear system
def construct_A_b(panels,alpha=0):
# construct matrix
N_panels = len(panels)
A = numpy.empty((N_panels, N_panels), dtype=float) # empty matrix
numpy.fill_diagonal(A, 0.5) # fill diagonal with 1/2
for i, p_i in enumerate(panels):
for j, p_j in enumerate(panels):
if i != j: # off-diagonals
A[i,j] = influence(p_i,p_j) # find influence
# computes the RHS
b = [-numpy.cos(alpha)*p.sx-numpy.sin(alpha)*p.sy for p in panels]
return [A,b]
In [293]:
# determine the vortex strength on a set of panels
def solve_gamma(panels,alpha=0):
A,b = construct_A_b(panels,alpha) # construct linear system
gamma = numpy.linalg.solve(A, b) # solve for gamma!
for i,p_i in enumerate(panels):
p_i.gamma = -gamma[i] # update panels
In [294]:
def make_wedge(t_c, N=128):
t=0
N_panels=N/3
x_ends=1.5*(numpy.append(numpy.linspace(-1,0,N_panels+1), [numpy.linspace(0,0,N_panels+1), numpy.linspace(0,-1,N_panels+1)])+0.5)
y_ends=1.5*numpy.append(numpy.linspace(0,-t_c/2,N_panels+1),[numpy.linspace(-t_c/2,t_c/2,N_panels+1), numpy.linspace(t_c/2,0,N_panels+1)])
wedge = numpy.empty(N_panels*3, dtype=object)
for i in range(N_panels*3+2):
if(x_ends[i]!=x_ends[i+1])|(y_ends[i]!=y_ends[i+1]):
wedge[t] = Panel(x_ends[i], y_ends[i], x_ends[i+1], y_ends[i+1])
t+=1
return wedge
In [296]:
wedge=make_wedge(1.,64)
solve_gamma(wedge)
plot_flow(wedge)
In [47]:
In [ ]: