1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
|
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from math import *
import xlrd
import scipy.integrate as integr
g=9.81
R=8.31 #J/K/mol
gamma = 1.40 #gaz parfait diatomique, coeff Laplace
m0=768789.3 #masse initiale en kg
r=5.425 #diamètre coiffe en m
Cx = 0.6 #supposé constant car vitesses élevées atteintes rapidement (graphe de reynolds, mach...)
workbook = xlrd.open_workbook("mettre le nom du fichier")
worksheet = workbook.sheet_by_index(0)
tps = worksheet.col_values(0)
pousse = worksheet.col_values(1)
masse = worksheet.col_values(2)
debit = worksheet.col_values(3)
l=len(tps)
tps, pousse, masse, debit = [], [],[], []
for i in range (l-1):
tps.append(float(worksheet.col_values(0)[i]))
pousse.append(float(worksheet.col_values(1)[i]))
masse.append(float(worksheet.col_values(2)[i]))
debit.append(float(worksheet.col_values(3)[i]))
tps1=np.zeros(l-2)
for i in range (l-2) :
tps1[i]=tps[i]
debitm, pas, v_ejec=[],[],[]
#calcul du débit massique
for i in range(l-2):
debitm.append(abs((masse[i]-masse[i+1])/(tps[i+1]-tps[i])))
#pas variable pour equa diff
pas=np.zeros(l)
for i in range (l-2):
pas[i]=(tps[i+1]-tps[i])
#température en deg, pression en pascal et masse volumique de l'air en fonction de l'altitude z
a=6.5*0.01 #en K/m
T0=(15+273.15) #en kelvin
p0= 101325 #en pascal
M=28.966*0.01 #en kg/mol
def temperature(z):
return -a*z + T0
def pression(z): #hypothèse d'un gradient de température constant
return p0*((1-z*a/T0)**(M*g/(R*a)))
def masse_volumique(z):
return M*pression(z)/(R*temperature(z))
def vitesse_son(z):
return sqrt(gamma*R*temperature(z)/masse_volumique(z))
#rk4fusee eau.xlsx
alt,vit=[0],[0]
z2=0
y,x =0,0.6371009*(10**(-19))
w=2*np.pi/(23*3600+59*60+4) #rad/s#vitesse de rotation de la terre(465,1 m/s )
for i in range (l-3):
if z2<100000:
k1 = -g + (pousse[i] - 0.5*Cx*maitre_couple(tps[i])*masse_volumique(z2)*v**2)/masse[i] + x*w**2
else:
k1 = -g + pousse[i]/masse[i] + x*w**2
if z2< 100000:
k2 = -g + ((pousse[i]+pousse[i+1])/2 - 0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*((y+(pas[i]*k1/2))**2))/(masse[i]+masse[i+1])/2 + x*w**2
else:
k2 = -g +(pousse[i]+pousse[i+1])/2 + x*w**2
if z2< 100000:
k3 = -g + ((pousse[i]+pousse[i+1])/2 -0.5*Cx*(np.pi*(5.425/2)**2 + 2*np.pi*(3.15/2)**2)*masse_volumique(z2)*(y+1/2*pas[i]*k2)**2)/(masse[i]+masse[i+1])/2 + x*w**2
else:
k3 = -g + (pousse[i]+pousse[i+1])/2 + x*w**2
if z2<100000:
k4 = -g + ((pousse[i+1]+(y+pas[i]*k3)*debitm[i+1])/masse[i+1]) + x*w**2
else:
k4 = -g + (poussee[i+1]/masse[i+1]) + x*w**2
y+=pas[i]*(k1/6 + k2/3 + k3/3 + k4/6)
z2+=y*pas[i]
x+=y*pas[i]
alt.append(z2)
vit.append(y) |
Partager