45.sagews

%python

'''Application of vortex lattice method to a planar symmetric wing
Check the reference : pp 271-276, Bertin and Smith, Aerodynamics for Engineers (2nd edition)
'''

# Author : Pankaj Pandey

try:
    import numpy
except ImportError, ierr:
    print '''This program requires 'numpy' and optionally 'pylab' (for plotting) python modules to run,
    get them from http://numpy.scipy.org/ and http://matplotlib.sourceforge.net/ respectively'''
    raise ierr

class VLMPlanarWing(object):
    '''class representing a wing'''
    def __init__(self, rootChord, tipChord, sweep, length, N=20):
        '''sweep should be in degrees'''
        self.rootChord = rootChord
        self.tipChord = tipChord
        self.sweep = sweep * numpy.pi / 180
        self.b = length
        self.N = int(N)
        self.genCoordinates()
        
    def getChord(self, y):
        return self.rootChord - (self.rootChord-self.tipChord)*y*2/self.b
    
    def genCoordinates(self):
        '''To be called after changing any geometric parameter'''
        N = self.N
        self.vortexPoints = numpy.empty((N+1,2))
        self.vortexPoints[:,1] = numpy.arange(N+1) * self.b / (2*N)
        self.vortexPoints[:,0] = self.vortexPoints[:,1] *numpy.tan(self.sweep) + (self.rootChord - (self.rootChord-self.tipChord)*numpy.linspace(0.0,1.0,N+1)) / 4.0
        self.ctrlPoints = numpy.empty((N,2))
        self.ctrlPoints[:,1] = numpy.arange(0.5,N+0.5,1.0) * self.b / (2*N)
        self.ctrlPoints[:,0] = (self.vortexPoints[:-1,0]+self.vortexPoints[1:,0])/2.0 + (self.rootChord - (self.rootChord-self.tipChord)*(numpy.arange(0.5,N+0.5,1.0)) /N ) / 2.0
        
        self.calcAMatrix()
        
    def calcAMatrix(self):
        '''Not to be externally called. Automatically called by genCoordinates()
        Calculates the coefficient matrix from 'Bertin and Smith' book's formulae, without the 4pi factor'''
        dy = self.b / 2 / self.N
        dx = dy * numpy.tan(self.sweep) + (self.tipChord-self.rootChord) / 4.0 / self.N
        dxr = dx
        dyr = - dy
        N = self.N
        A = numpy.empty((N,N))
        PI = numpy.pi
        hypot = numpy.hypot
        for m in xrange(N):
            for n in xrange(N):
                dx1, dy1 = self.ctrlPoints[m] - self.vortexPoints[n]
                dx2, dy2 = self.ctrlPoints[m] - self.vortexPoints[n+1]
                A[m,n] = ( ((dx*dx1+dy*dy1)/hypot(dx1,dy1) - (dx*dx2+dy*dy2)/hypot(dx2,dy2)) / (dx1*dy2-dx2*dy1)
                            - (1+dx1/hypot(dx1,dy1))/dy1
                            + (1+dx2/hypot(dx2,dy2))/dy2 )
                # for the other symmetric part of the wing
                dy1r = self.ctrlPoints[m,1] + self.vortexPoints[n,1]
                dy2r = self.ctrlPoints[m,1] + self.vortexPoints[n+1,1]
                # -ve sign because in the symmetric part (reverse direction) our (n+1)th point is to left of (n)th point 
                A[m,n] -= ( ((dxr*dx1+dyr*dy1r)/hypot(dx1,dy1r) - (dxr*dx2+dyr*dy2r)/hypot(dx2,dy2r)) / (dx1*dy2r-dx2*dy1r)
                            - (1+dx1/hypot(dx1,dy1r))/dy1r
                            + (1+dx2/hypot(dx2,dy2r))/dy2r )
        self.A = A
        return self.A
    
    def solve(self, U_inf=1.0, AoA=180/numpy.pi, rho_inf=1.0):
        '''AoA in degrees'''
        AoA *= numpy.pi / 180.0
        rhs = - U_inf * AoA * numpy.ones((self.N,)) * numpy.pi * 4
        circulation = numpy.linalg.solve(self.A, rhs)
        sectionalLift = rho_inf * U_inf * circulation
        lift = numpy.sum(sectionalLift) * self.b / self.N
        C_l = lift * 4 / rho_inf / U_inf / U_inf / (self.rootChord + self.tipChord) / self.b
        C_l_alpha = C_l / AoA
        self.C_L_alpha = C_l_alpha
        sectional_chord = self.rootChord - (self.rootChord - self.tipChord) * numpy.linspace(0,1.0,self.N+1)
        sectional_control_chord = (sectional_chord[1:] + sectional_chord[:-1]) / 2.0
        sectional_C_l =  sectionalLift * 2 / rho_inf / U_inf / U_inf / sectional_control_chord
        self.sectional_C_l = sectional_C_l
        return circulation, lift, sectionalLift, sectional_C_l, C_l, C_l_alpha
    
    def plot(self):
        '''Use only after calling solve()'''
        try:
            import pylab
        except ImportError:
            print 'Matplotlib not installed. Plotting is disabled'
            return
        pylab.figure()
        pylab.plot(self.ctrlPoints[:,1], self.sectional_C_l, '.-', label='$Sectional\ C_l$ \n $C_{L\\alpha}=%.4f$' % self.C_L_alpha)
        pylab.title('Vortex Lattice Method for planar wing')
        pylab.xlabel('wing semi-span')
        pylab.ylabel('$C_l$')
        pylab.legend(loc='best')
        pylab.savefig('plot.png')

def printdatadict(datadict):
    for key,value in datadict.iteritems():
        print '%25s : %s' %(key,value)

def printdata(data):
    keys = ('Circulation', 'Lift', 'Sectional Lift', 'Sectional C_l', 'C_l', 'C_l_alpha(/radian)')
    print ''
    for i in xrange(len(keys)):
        print '%-19s: %s' %(keys[i],data[i])

def get_result(inputdata):
    wing = VLMPlanarWing(*[float(i) for i in inputdata[:5]])
    data = wing.solve(*[float(i) for i in inputdata[5:]])
    print 'C_L_alpha = ', wing.C_L_alpha
    wing.plot()

@interact
def _(rootChord=0.2, tipChord=0.2, sweep=45.0, length=1.0, N=5, U_inf=1.0, AoA=57.295779513082323, rho_inf=1.0):
    args = rootChord, tipChord, sweep, length, N, U_inf, AoA, rho_inf
    get_result(args)

# input data for the example problem in the 'Bertin and Smith' Book
inputdata = [0.2, 0.2, 45.0, 1.0, 4, 1.0 , 180/numpy.pi, 1.0]
get_result(inputdata)

# Expected approximate results from the example problem in 'Bertin and Smith' book
circulationE = numpy.array([ 0.34281059,  0.36052917,  0.35701059,  0.28261768])
liftE = 0.343313
sectionalLiftE = numpy.array([ 0.34281059,  0.36052917,  0.35701059,  0.28261768])
ClE = 3.43313
Cl_alphaE = 3.43313

C_L_alpha = 3.44422418771

# rootChord tipChord sweep(degrees) length N U_inf AoA(degrees) rho_inf
inputdata = 0.2, 0.1, 45.0, 1.0, 4, 1.0 , 180/numpy.pi, 1.0
get_result(inputdata)

C_L_alpha = 3.85471901591