Demo entry 6327588

codododo

   

Submitted by anonymous on Nov 24, 2016 at 01:01
Language: Python. Code size: 2.8 kB.

#coding uft-8

import numpy as np
import math
import matplotlib.pyplot as plt

oper = np.genfromtxt('operating_point.csv',delimiter = ',')
T = oper[0,1] + 273.15 #Kelvin
P = oper[1,1]*10**3#Pa
manometer_rho = 800#kg/m3
manometer_angle = oper[3,1]*math.pi/180
gSI = 9.81#m/s2
air_rho = P/(287*T)#kg/m3
divsize = 0.0254 #m

#Importing datasets
paths = ['notrip_dataset.csv','trip_dataset.csv']
desc = ['smooth cylinder', 'cylinder with tripwire']

#vardef
Cd =Cl = Cdstar = [0,0,0]

Cl = [0,0,0]
for i in range(len(paths)):
    DataSet = np.genfromtxt(paths[i],delimiter = ',')
    #Stripping headers,atm conditions and pitot-static measuremrents
    atm = DataSet[-1,2]*divsize
    PStat = DataSet[29:31,:]*divsize
    DataSet = DataSet[1:29,:]
    angles = DataSet[1:,1]*math.pi/180.0

    #computing the pressure coeficient, Cp
    hs = DataSet[1:,2] * divsize
    Cp = (hs - np.ones(np.size(hs))*PStat[0,2])/(PStat[1,2]-PStat[0,2])
    #Plotting with standard french style
    plt.figure(1)
    if i==0:
        plt.plot(angles,Cp,'k+')
    else:
        plt.plot(angles,Cp,'kx')
    plt.figure(2)
    if i==0:
        plt.plot(angles,Cp,'k^', markerfacecolor='none')
    else:
        plt.plot(angles,Cp,'ks', markerfacecolor='none')

    ##CD coeficients

    Cd[i] = np.trapz(np.multiply(Cp,np.cos(angles)), x=angles)
    Cl[i] = 0

    #from Roshko 1960
    q = (PStat[1,2]-PStat[0,2])*manometer_rho*math.cos(manometer_angle)*gSI
    V = np.sqrt(np.abs(2*(Cp*q + (PStat[0,2] - PStat[1,2])*manometer_rho
        *math.cos(manometer_angle)*gSI/air_rho)))
    Vstar = V * (1 + 0.25*Cd[i]*(0.102/0.46) + 0.82*(0.102/0.46)**2)
    Cdstar[i] = Cd[i]*(1 - 0.50*Cd[i]*(0.102/0.46) - 2.50*(0.102/0.46)**2)
    Cpstar = np.multiply(np.power(np.divide(V,Vstar), np.ones(np.size(Cp))*2),
                         Cp-1)+1
    if i==0:
        plt.plot(angles,Cpstar,'k+')
    else:
        plt.plot(angles,Cpstar,'kx')

#finalazing the graph
plt.figure(1)
plt.grid(True)
plt.title('''Pressure around the cylinder for theoretical,
              smooth and tripwired cylinder''')
plt.xlabel(r'Angle around the cylinder [$-$]')
plt.ylabel(r'Pressure Coeficient [$-$]')

#adding standard theoretical distribution
theory_Cp_funk = np.vectorize(lambda x: 1 - 4*math.sin(x)**2)
theory_angle = np.arange(np.min(angles), np.max(angles), step=0.01)
theory_Cp = theory_Cp_funk(theory_angle)

Cd[2] = np.trapz(np.multiply(theory_Cp_funk(theory_angle),
                    np.cos(theory_angle)), x=theory_angle)

plt.plot(theory_angle,theory_Cp,'k--')


#legend
plt.figure(1)
plt.legend(['smooth','tripwire','theoretical'])
plt.figure(2)
plt.grid(True)
plt.legend(['smooth uncorrected','tripwire uncorrected','smooth corrected',
            'tripwire corrected'])

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).