Assignment

Circular Cylinder


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]:
<matplotlib.text.Text at 0xbeae7b8>

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)


When N = 32
Circle frictional coefficient = 4.04e-03
Flate plate: 2.37e-03
When N = 64
Circle frictional coefficient = 4.06e-03
Flate plate: 2.37e-03
When N = 128
Circle frictional coefficient = 4.07e-03
Flate plate: 2.37e-03
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]:
<matplotlib.text.Text at 0xbeafeb8>

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]:
<matplotlib.collections.PathCollection at 0x169b1b38>

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))


When N=32, Drag coefficient= 2.42e+00
When N=64, Drag coefficient= 2.38e+00
When N=128, Drag coefficient= 2.36e+00

Jukowski foil


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]:
<matplotlib.legend.Legend at 0x17541ba8>

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]:
<matplotlib.text.Text at 0x178fee10>

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')


Wedge


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 [ ]: