问题描述
import numpy as np
import math
import matplotlib.pyplot as plt
%matplotlib inline
#Class to calculate CAmber and Thickness distribution
class naca: # defining class naca
def __init__(self,airfoil):
self.airfoil = airfoil
print ('NACA Airfoil is',airfoil,'\n')
self.__maxcAmber=int(airfoil[0:1])/100
print ('Maximum CAmber of airfoil as a %age of chord is',self.__maxcAmber *100)
self.__pos=int(airfoil[1:2])/10
print ('Position of maximum cAmber from leading edge as a %age of chord is',self.__pos*100)
self.__thick=int(airfoil[2:4])/100
print ('Maximum Thickness of airfoil as a %age of chord is',self.__thick * 100,'\n')
#Function to calculate CAmber
def cAmber(self,x): # defining cAmber function
if x < self.__pos and x>=0:
z=(self.__maxcAmber/(self.__pos**2))*((2*self.__pos*x)-x**2)
return z
elif x>=self.__pos and x<=1:
z=(self.__maxcAmber/((1-self.__pos)**2))*(1-(2*self.__pos)+(2*self.__pos*x)-x**2)
return z
else:
raise Exception ('Entered value of x/c should be between 0 to 1')
#Function to calculate Thickness
def thickness (self,x): # defining thickness function
if 0<=x<=1:
t= 5*self.__thick*(0.2969*math.sqrt(x) - 0.1260*x - 0.3516*x**2 + 0.2843*x**3 - 0.1036*x**4)
return t
else:
raise Exception ('Entered value of x/c should be between 0 to 1')
#Class to calculate x-derivative of CAmber and Thickness distribution function
class derivative(naca): # parent class naca mentioned in the bracket
def __init__(self,airfoil):
super().__init__(airfoil) # calling the prevIoUs class attributes
# Finding x-derivative of cAmber
def d_cAmber(self,xa):
h=0.01
if xa==0:
dx=(super().cAmber(xa+h)-super().cAmber(xa))/h # forward differencing
return dx
elif xa==1:
dx=(super().cAmber(xa)-super().cAmber(xa-h))/h # backward differencing
return dx
else:
dx=(super().cAmber(xa+h)-super().cAmber(xa-h))/(2*h) # central differencing
return dx
# Finding x-derivative of thickness
def d_thick(self,xa):
h=0.01
if xa==0:
dt=(super().thickness(xa+h)-super().thickness(xa))/h # forward differencing
return dt
elif xa==1:
dt=(super().thickness(xa)-super().thickness(xa-h))/h # backward differencing
return dt
else:
dt=(super().thickness(xa+h)-super().thickness(xa-h))/(2*h) # central differencing
return dt
#Class to calculate flow around an arbitary airfoil using thin-airfoil theory
class thinairfoiltheory(derivative): #parent class mentioned in the bracket
def __init__(self,airfoil):
super().__init__(airfoil) #calling the prevIoUs class attributes
def sectional(self,alphad,f):
n=15 ### it is not taking value more than 15. n is actually related to the number of division used in intergration
##check how to give n=100.For n>15 it is giving error???
#we have to increase thye value of n for accuracy
alphar=alphad*(np.pi/180)
thetal=0
thetah=np.pi
h=(np.pi-0)/n
th=0
sum1=0
sum2=0
sum3=0
for i in range(1,n): #composite trapezoidal rule used for integration
th=th+h
x=0.5*(1-np.cos(th))
sum1=sum1+super().d_cAmber(x)
sum2=sum2+(super().d_cAmber(x)*np.cos(th))
sum3=sum3+(super().d_cAmber(x)*np.cos(2*th))
xl=0.5*(1-np.cos(thetal))
xh=0.5*(1-np.cos(thetah))
finalsum1=(h/2)*(super().d_cAmber(xl)+super().d_cAmber(xh))+(h*sum1)
a0=(1/np.pi)*finalsum1
finalsum2=(h/2)*(((super().d_cAmber(xl))*np.cos(thetal))+((super().d_cAmber(xh))*np.cos(thetah)))+(h*sum2)
a1=(2/np.pi)*finalsum2
finalsum3=(h/2)*(((super().d_cAmber(xl))*np.cos(2*thetal))+((super().d_cAmber(xh))*np.cos(2*thetah)))+(h*sum3)
a2=(2/np.pi)*finalsum3
cl=2*(np.pi*(alphar-a0)+(np.pi/2)*a1)
cmc=-(np.pi/4)*(a1-a2)
#print("sectional lift coefficient",cl)
#print("quatercord point moment",cmc)
if f==1:
return cl
elif f==2:
return cmc
def plot(self,ff):
deg=np.linspace(0,30,50)
liftcoff=self.sectional(deg,1)
momcoff=self.sectional(deg,2)
if ff==1: #graph for lift coefficient
plt.plot(deg,liftcoff)
plt.rcParams['figure.figsize'] = (10,6) #This fixes the size of plot displayed below
plt.rcParams.update({'font.size': 15}) #This fixes the font size of parameters in the plot below
plt.xlabel('Angle of attack in degress') #This labels the x-axis
plt.ylabel('Liff coeff') #This labels the y-axis
plt.grid(b=True,color='black',alpha=0.3,linewidth=1)
elif ff==2: #graph for moment coefficient
##the moment coeff graph is not plotting.I cannot figure out the problem???
plt.plot(deg,momcoff)
plt.rcParams['figure.figsize'] = (10,6) #This fixes the size of plot displayed
plt.rcParams.update({'font.size': 15}) #This fixes the font size of parameters
plt.xlabel('Angle of attack in degress') #This labels the x-axis
plt.ylabel('mom coeff') #This labels the y-axis
plt.grid(b=True,linewidth=1)
解决方法
您的错误是在代码cmc=-(np.pi/4)*(a1-a2)
的这一行上,对于ff = 2,您的section()函数将返回 cmc ,这是一个常量。我认为您的公式必须缺少 alphar 才能计算出不同角度的力矩系数。