Блог

 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, color='black', alpha=0.3,linewidth=1)
  

Ваша ошибка находится в этой строке кода cmc=-(np.pi/4)*(a1-a2) , для ff = 2 ваша функция sectional() вернет cmc, который является константой. Я думаю, что в вашей формуле должен отсутствовать алфавит, чтобы вычислить коэффициент момента для разных углов.