Convert Python to Matlab
lavanda993
1.pdf
#Import standard math functions from math import * import numpy as np import matplotlib.pyplot as plt #Prandtl-Meyer expansion angle function: def nu(M,y): return sqrt((y+1)/(y-1))*atan(sqrt((y-1)*(M**2-1)/(y+1)))- atan(sqrt(M**2-1)) #Prandtl-Meyer Expansion angle for numeric extraction of Mach number def nu2(M,nu,y): return sqrt((y+1)/(y-1))*atan(sqrt((y-1)*(M**2-1)/(y+1)))- atan(sqrt(M**2-1))-nu #Derivative of Prandtl-Meyer expansion angle function def nudif(M, y = 1.25): return sqrt(M**2-1)/(M*(1+(y-1)*M**2/2)) #Stagnation pressure from freestream pressure. def Pstag(Pe,M,y): return Pe*(1+(y-1)/2*M**2)**(y/(y-1)) #Stagnation pressure ratio across a shockwave with P01 calculated from Pstat and M. def Pstagratio(M,P02,P1,y): P01 = Pstag(P1,M,y) f1 = 2/((y+1)*(y*M**2-(y-1)/2)**(1/(y-1))) f2 = (((y+1)/2*M)**2/(1+(y-1)/2*M**2))**(y/(y-1)) return f1*f2-P02/P01 #Stagnation pressure ratio across a shockwave with P01 calculated from Pstat and M. def Pstagshock(M,P1,y):#returns P02 P01 = Pstag(P1,M,y) f1 = 2/((y+1)*(y*M**2-(y-1)/2)**(1/(y-1))) f2 = (((y+1)/2*M)**2/(1+(y-1)/2*M**2))**(y/(y-1)) return f1*f2*P01 #Pressure across an oblique shockwave as a function of upstream pressure, gamma, and Mach. def Pstatshock(P,M,y,theta):
beta = shockangle(M,y,theta)[1] return P*(1+2*y/(y+1)*((M*sin(beta))**2-1)) #Numerical derivative for this particular function with only x (which is M) varying. def centraldiff(func, x, A, Ast, y, h = 0.0000001): return (func(x+h,A,Ast,y)-func(x-h,A,Ast,y))/(2*h) #Newton's method for four-term fuction. def Newts(func, fin,y,A,Ast, err, iMax = 100000): x1 = fin it = 0 ea = np.Infinity while ea > err and it < iMax: x0 = x1 x1 = x0 - func(x0,A,Ast,y)/ (centraldiff(func,x0,A,Ast,y)) #x1 = x0 - func(x0,A,Ast,y)/(RPdif(x0,y)) it +=1 ea = abs((x0-x1)/x1) return x1 #Newton's method def Newt(Mguess,nu,gamma,err,iMax=1000): x1 = Mguess it = 0 ea = 100000 while ea > err and it < iMax: x0 = x1 x1 = x0 - nu2(x0,nu,gamma)/nudif(x0,gamma) it +=1 ea = abs((x0-x1)/x1) return x1 #Pressure across expansion fan def Pexpand(M1, M2, P1, y): top = 1+(y-1)/2*M1**2 bot = 1+(y-1)/2*M2**2 return P1*(top/bot)**(y/(y-1))
#Explicit solution for beta angle. Requires theta input in degrees. def shockangle(M,y,thd): th = thd*pi/180 lt1 = (M**2-1)**2 lt2 = (1+(y-1)/2*M**2) lt3 = (1+(y+1)/2*M**2) l = (lt1-3*lt2*lt3*tan(th)**2)**.5 xt1 = (M**2-1)**3 xt2 = lt2 xt3 = (1+(y-1)/2*M**2+(y+1)/4*M**4) x = (xt1-9*xt2*xt3*tan(th)**2)/l**3 tanb = [] for i in range(2): tt1 = (M**2-1) tt2 = 2*l*cos((4*pi*float(i)+acos(x))/3) tt3 = 3*(1+(y-1)/2*M**2)*tan(th) tanb.append((tt1+tt2)/tt3) return (atan(tanb[0]),atan(tanb[1])) #Mach number after normal shockwave def normshock(M1,y): top = (1+(y-1)/2*M1**2) bot = y*M1**2-(y-1)/2 return (top/bot)**.5 #Loop function for finding freestream Mach number def Maftershock(M1,theta,gamma): th = theta*pi/180 beta = shockangle(M1,gamma,theta)[1] Mn1 = M1*sin(beta) Mn2 = normshock(Mn1,gamma) #if Mn2/sin(beta-th) < 1: # beta = shockangle(M1,gamma,theta)[0] # Mn1 = M1*sin(beta) # Mn2 = normshock(Mn1,gamma) return Mn2/sin(beta-th) #Calculate drag and lift coefficients.
def Cd(Ps, Pinf, delta, alpha, M, gamma): d = delta*pi/180 a = alpha*pi/180 part1 = 1/(2*cos(d))*1/(gamma*M**2/2) part2 = (Ps[0]/Pinf-Ps[3]/Pinf) part3 = (Ps[1]/Pinf-Ps[2]/Pinf) part4 = (Ps[2]/Pinf-Ps[1]/Pinf) Cdrag = part1*(part2*sin(d+a)+part3*sin(d-a))#Might need an abs here, but I kinda doubt it. Clift = part1*(part2*cos(d+a)+part4*cos(d-a)) return (Cdrag,Clift) #Wrapper to solve for the Mach number after expansion fan. def PmMach(M,theta,gamma): vM1 = nu(M,gamma) vM2 = theta*pi/180+vM1#Theta is assumed to be in degrees, so this converts it to radians. M2 = Newt(2,vM2,gamma,0.00001) return M2 #Calculate pressures on the four edges of the diamond air foil. def getPressures(Pinf, delta, alpha, M, gamma): M2 = Maftershock(M,delta+alpha,gamma) P2 = Pstatshock(Pinf,M,gamma,delta+alpha) #if abs(delta-alpha) <= 0.00001: # M3 = M # P3 = Pinf if alpha >= delta:#Remember to get rid of the try things before this is sent #try: M3 = PmMach(M,alpha-delta,gamma) #except ValueError: # print(str(M)+" "+str(alpha)+" "+str(Pinf)) P3 = Pexpand(M,M3,Pinf,gamma) else: M3 = Maftershock(M,delta-alpha,gamma) P3 = Pstatshock(Pinf,M,gamma,delta-alpha) #try:
M4 = PmMach(M2,2*delta,gamma) #except ValueError: # print(str(M)+" "+str(alpha)+" "+str(Pinf)) P4 = Pexpand(M2,M4,P2,gamma) #try: M5 = PmMach(M3,2*delta,gamma) #except ValueError: # print(str(M)+" "+str(alpha)+" "+str(Pinf)) P5 = Pexpand(M3,M5,P3,gamma) return ((P2,P3,P4,P5),(M2,M3,M4,M5)) #Ambient pressure of the atmosphere at altitudes from 10km to 50km in 10km increments. Palt = [2.650E+4,5.529E+3,1.197E+3,2.871E+2,7.977E+1] #Sweep through alpha, M, and altitudes. LDalt = [[[] for j in range(40)] for i in range(5)] Alphalt = [[[] for j in range(40)] for i in range(5)] Pressurealt = [[[] for j in range(40)] for i in range(5)] Machsalt = [[[] for j in range(40)] for i in range(5)] Lalt = [[[] for j in range(40)] for i in range(5)] Dalt = [[[] for j in range(40)] for i in range(5)] delta = 5 for j in range(5): for k in range(40): for i in range(151): Machin = float(k)/4.0+1.5 alpha = float(i)/10.0 try: Press = getPressures(Palt[j],delta,alpha,Machin, 1.4) Pressures = Press[0] Machs = Press[1] except ValueError: break CD = Cd(Pressures,Palt[j],delta,alpha,Machin,1.4) Lalt[j][k].append(CD[1]) Dalt[j][k].append(CD[0])
LDalt[j][k].append(CD[1]/CD[0]) Alphalt[j][k].append(alpha) Pressurealt[j][k].append(Pressures) Machsalt[j][k].append(Machs) LDmax = [[[],[]] for i in range(5)] for i in range(len(LDalt)): for j in range(len(LDalt[i])): LDmax[i][0].append(max(LDalt[i][j])) LDmax[i][1].append(float(j)/4.0+1.5) plt.figure(0) for i in range(len(LDalt)): plt.plot(Alphalt[i][2],LDalt[i][2]) Presstat = [[] for i in range(4)] Machstat = [[] for i in range(4)] for j in range(len(Pressurealt[0][2])): for i in range(4): Machstat[i].append(Machsalt[0][2][j][i]) Presstat[i].append(Pressurealt[0][2][j][i]/Palt[0]) plt.ylabel("L/D") for i in range(4): plt.figure(i) plt.xlabel("Alpha (Degrees)") plt.figure(1) plt.plot(Alphalt[0][2],Presstat[0],Alphalt[0] [2],Presstat[1],Alphalt[0][2],Presstat[2],Alphalt[0] [2],Presstat[3]) plt.ylabel("Ramp Pressure - P/Pinf") plt.legend(['P2/Pinf','P3/Pinf','P4/Pinf','P5/Pinf']) plt.figure(2) plt.plot(Alphalt[0][2],Lalt[0][2],Alphalt[0][2],Dalt[0][2]) plt.ylabel("Lift and Drag Coefficients") plt.legend(['Lift Coefficient','Drag Coefficient']) plt.figure(3) plt.plot(Alphalt[0][2],Machstat[0],Alphalt[0] [2],Machstat[1],Alphalt[0][2],Machstat[2],Alphalt[0] [2],Machstat[3])
plt.ylabel("Mach Number") plt.legend(['M2','M3','M4','M5']) plt.figure(5) plt.plot(LDmax[0][1],LDmax[0][0]) plt.xlabel("Mach Number") plt.ylabel("L/D Max") #plt.show() #Add variable viscosity, density, temperature, and speed of sound. Rhos = [4.135E-1,8.891E-2,1.841E-2,3.995E-3,1.027E-3] mus = [14.58e-6,14.22e-6,14.75e-6,16.01e-6,17.04e-6] cs = [299.5,295.1,301.7,317.2,329.8] Tamb = [223.3,216.6,226.5,250.4,270.6] Cpair = 1003#kJ/kg*K Dfalt = [[[] for j in range(40)] for i in range(5)] LDfalt = [[[] for j in range(40)] for i in range(5)] Res = [[[] for j in range(40)] for i in range(5)] def Tavgturb(T,V,Cp,Rf): return T+V**2/(2*Cp)*(Rf-7/9) def Tavgtu(T,M,y): return T*(1+2/9*((y-1)/2*M**2)) def Tavglam(T,M,y): return T*(1+7/15*((y-1)/2*M**2)) for i in range(len(Dalt)): for j in range(len(Dalt[i])): for k in range(len(Dalt[i][j])): Vinf = (float(j)/4.0+1.5)*cs[i] Tavg = Tavgtu(Tamb[i],Vinf/cs[i],1.4) #Tavg = Tavgturb(Tamb[i],Vinf/cs[i],Cpair,1) Re = (Rhos[i]*Vinf*2*cos(float(k)/10*pi/180)/ cos(delta*pi/180))/mus[i] Res[i][j].append(Re) Cdincomp = 7/(225*Re**(1/7)) Cdcomp = 7/(225*(Re**(1/7)*((Tamb[i]/ Tavg)**(5/2)*((Tavg+120)/(Tamb[i]+120)))**(1/7)))
Dfalt[i][j].append(Dalt[i][j][k]+ 2*Cdcomp)#Cdincomp/((Tamb[i]/Tavg)**(5/2)*((Tavg+120)/(Tamb[i] +120))**(1/7)))#My indexing might be off here. LDfalt[i][j].append(Lalt[i][j][k]/Dfalt[i][j][k]) plt.figure(6)#Something's wrong with the friction drag. It's getting too high too fast. for i in range(len(LDfalt)): plt.plot(Alphalt[i][26],LDfalt[i][26]) plt.xlabel("Alpha (Degrees)") plt.ylabel("L/D") plt.legend(['10 km','20 km','30 km','40 km','50 km']) LDfmax = [[[],[]] for i in range(5)] for i in range(len(LDfalt)): for j in range(len(LDfalt[i])): LDfmax[i][0].append(max(LDfalt[i][j])) LDfmax[i][1].append(float(j)/4.0+1.5) plt.figure(7) for i in range(len(LDfmax)): plt.plot(LDfmax[i][1],LDfmax[i][0]) plt.xlabel("Mach number") plt.ylabel("L/D max") plt.legend(['10 km','20 km','30 km','40 km','50 km']) Dflalt = [[[] for j in range(40)] for i in range(5)] LDflalt = [[[] for j in range(40)] for i in range(5)] Rels = [[[] for j in range(40)] for i in range(5)] for i in range(len(Dalt)): for j in range(len(Dalt[i])): for k in range(len(Dalt[i][j])): Vinf = (float(j)/4.0+1.5)*cs[i] #Tavg = Tavgturb(Tamb[i],Vinf/cs[i],Cpair,1) Re = (Rhos[i]*Vinf*2*cos(float(k)/10*pi/180)/ cos(delta*pi/180))/mus[i] Rels[i][j].append(Re) if Re < 500000: #Laminar Tavg = Tavglam(Tamb[i],Vinf/cs[i],1.4)
Cdincomp = 1.328/(Re**(1/2)) Cdcomp = Cdincomp/((Tamb[i]/ Tavg)**(5/2)*((Tavg+120)/(Tamb[i]+120)))**(1/2) elif Re > 500000 and Re < 10000000: #Transitional Tavg = Tavgtu(Tamb[i],Vinf/cs[i],1.4) Cdincomp = 1/Re**(1/7)*(-3.116/ Re**(5/14)+0.04096) Cdcomp = Cdincomp/((Tamb[i]/ Tavg)**(5/2)*((Tavg+120)/(Tamb[i]+120)))**(1/7) else: #Turbulent Tavg = Tavgtu(Tamb[i],Vinf/cs[i],1.4) Cdincomp = 7/(225*Re**(1/7)) Cdcomp = 7/(225*(Re**(1/7)*((Tamb[i]/ Tavg)**(5/2)*((Tavg+120)/(Tamb[i]+120)))**(1/7))) Dflalt[i][j].append(Dalt[i][j][k]+ 2*Cdcomp)#Cdincomp/((Tamb[i]/Tavg)**(5/2)*((Tavg+120)/(Tamb[i] +120))**(1/7)))#My indexing might be off here. LDflalt[i][j].append(Lalt[i][j][k]/Dflalt[i][j][k]) LDflmax = [[[],[],[]] for i in range(5)] def specmax(l): Ma = 0 idx = 0 for i in range(len(l)): if l[i] > Ma: Ma = l[i] idx = i return (Ma,idx) for i in range(len(LDflalt)): for j in range(len(LDflalt[i])): LDflmax[i][0].append(specmax(LDflalt[i][j])[0]) LDflmax[i][1].append(float(j)/4.0+1.5) LDflmax[i][2].append(Rels[i][j][specmax(LDflalt[i][j]) [1]]) plt.figure(8)
for i in range(len(LDfmax)): plt.plot(LDflmax[i][1],LDflmax[i][0]) plt.xlabel("Mach number") plt.ylabel("L/D max") plt.legend(['10 km','20 km','30 km','40 km','50 km']) plt.figure(9) for i in range(len(LDfmax)): plt.plot(LDflmax[i][2],LDflmax[i][0]) plt.xscale('log') plt.xlabel("Reynolds Number") plt.ylabel("L/D max") plt.legend(['10 km','20 km','30 km','40 km','50 km']) plt.axvline(x=500000, linestyle='--') plt.axvline(x=10000000, linestyle='--') plt.show()
