<<<<<<< HEAD import numpy as np from numpy import linalg as LA def initial_vel(latitude,retrograde=False): """ Initial velocity (m.s-1) depending on the launch point latitude (°) """ T_side = 86164.10 # sidereal day (s) R_T = 6378e3 # earth raduis (m) initial_velocity = (2*np.pi/T_side)*R_T*np.cos(np.deg2rad(np.abs(latitude))) if retrograde: return -initial_velocity else: return initial_velocity def final_target_vel(h, a): """ inputs: h: injection altitude (m) a: orbit semi major axis (m) output: orbit absolute injection velocity (m.s-1) """ mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2) R_T = 6378e3 # earth raduis (m) return np.sqrt(2*mu*((1/(R_T+h))-(1/(2*a)))) def staging_optim(isp1,isp2,eps1,eps2,Vc,Mp): """ Optimise mass loading in the case of two stage rocket inputs: isp1, isp2 : specific impulse of the 1st and 2nd stage (s) eps1, eps2 : structural coefficient of the 1st and 2nd stage (adimentional) Vc : target critical velocity (m.s-1) Mp : payload mass (Kg) outputs: M1 : mass of the 1st stage M2 : mass of the 2nd stage NB: structure mass structural coefficient = -------------------------------- structure mass + propellant mass ctritical velocity = final velocity - initial velocity + losses """ # Notation : # # Mt = total mass # q1 = M1/Mt # q2 = M2/Mt # qp = Mp/Mt # omega = qp + q2 g0 = 9.81 # standard gravity (m.s-2) omega_tab=np.linspace(0.001,0.9,100) # array of omega to test qp_tab=[] for omega in omega_tab: q2 = omega+eps1*(1-omega) q2 = np.power(q2,isp1/isp2)*np.exp(Vc/(g0*isp2)) q2 = 1/q2 q2 = 1-q2 q2 *= omega/(1-eps2) qp = omega-q2 qp_tab.append(qp) qp=max(qp_tab) omega=omega_tab[qp_tab.index(qp)] q2=omega-qp q1=1-q2-qp M1=(Mp/qp)*q1 M2=(Mp/qp)*q2 return M1, M2 def standard_atmosphere(altitude): """ Compute atmospheric properties according to the International Standaard Atmosphere (ISA) model input: altitude (m) outputs: temperature (K) pressure (Pa) density (Kg.m-3) """ assert altitude < 84852 , str(altitude)+' is higher than max altitude of the ISA model' g0 = 9.81 # standard gravity (m.s-2) r = 287 # specific gas constant (J.Kg-1.K-1) # ISA model h0_tab = np.array([0, 11e3, 20e3, 32e3, 47e3, 51e3, 71e3, 84852]) # base altitude (m) a_tab = np.array([-6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2.0e-3]) # lapse rate (K.m-1) # the base properties of each range are tabulated to reduce computational cost T0_tab = np.array([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65]) # base temperature (K) P0_tab = np.array([101325, 22632, 5474.9, 868.02, 110.91, 66.939, 3.9564]) # base pressure (Pa) index = 0 while altitude > h0_tab[index+1]: index += 1 a = a_tab[index] T0 = T0_tab[index] P0 = P0_tab[index] T = T0+a*(altitude-h0_tab[index]) # Toussaint equation if T == T0 : # isohermal region P = P0*np.exp(-g0*(altitude-h0_tab[index])/(r*T)) else : # gradient region P = P0*np.power(T/T0,-g0/(a*r)) rho = P/(r*T) return T,P,rho def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r,tolerance=1e-5): """ Compute the thrust assuming a started nozzle inputs: Pt : total pressure in the combustion chambre (Pa) Tt : total temperature in the combustion chambre (K) Pamb : ambient pressure (Pa) Ath : throat section area (m2) Ae : exit section area (m2) gamma : ratio of specific heats (adimentional) r : specific gas constant (J.Kg-1.K-1) tolerance : outputs: Thrust (N) Specific impulse (s) N.B : Thrust = mass flow rate * exit velocity + (exit pressure - ambiant pressure) * exit area """ assert gamma > 1, 'gamma must be > 1' g0 = 9.81 # standard gravity (m.s-2) A_ratio = lambda M,gamma: (1/M)*np.power((2/(gamma+1))*(1+((gamma-1)/2)*M*M),(gamma+1)/(2*gamma-2)) # compute the nozzle exit section area corresponding to the adapted regime Madapted = np.sqrt((2/(gamma-1))*(np.power(Pamb/Pt,(1-gamma)/gamma)-1)) Aadapted = Ath*A_ratio(Madapted,gamma) if Aadapted > Ae: # the nozzle is under-expanded # compute the exit Mach number iteratively using dichotomy algorithm A_ratio_target = Ae/Ath Me_upper = Madapted Me_lower = 1 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error > 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure elif Aadapted < Ae: # the nozzle is over-expanded Mshock = 1 step = 0.1 A_ratio_target = Ae/Ath Me_upper = 1 Me_lower = 0 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error < 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) while (Pe-Pamb)*1e-5 > tolerance: Mshock_old = Mshock Mshock += step Pe_old = Pe Pt_ratio = np.power(((gamma+1)*Mshock*Mshock)/((gamma-1)*Mshock*Mshock+2),gamma/(gamma-1))*np.power((gamma+1)/(2*gamma*Mshock*Mshock-gamma+1),1/(gamma-1)) Pi = Pt*Pt_ratio Ai = Ath/Pt_ratio A_ratio_target = Ae/Ai Me_upper = 1 Me_lower = 0 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error < 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) Pe = Pi*np.power(Tt/Te,gamma/(1-gamma)) if (Pe > Pamb) and (A_ratio(Mshock,gamma) > Ae/Ath): # the normal shock is outside the divergent and all the flow inside is supersonic A_ratio_target = Ae/Ath Me_upper = Mshock Me_lower = 1 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error > 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure break if Pe < Pamb: Pe = Pe_old Mshock = Mshock_old step /= 2 else: # the nozzle is adapted Me = Madapted Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature Pe = Pamb Ve = Me*np.sqrt(gamma*r*Te) # exit velocity # maximum mass flow rate in the case of a started nozzle Q = np.power(2/(gamma+1),(gamma+1)/(2*(gamma-1)))*np.sqrt(gamma/(Tt*r))*Pt*Ath Thrust = Q*Ve+(Pe-Pamb)*Ae isp = Thrust/(Q*g0) return Thrust,isp def adiabatic_flame(mech,T,P,X): """ Compute adiatabic flame temperature inputs: mech : chemical mechanism in cti format T : initial temperature (K) P : pressure (Pa) X : initial chemical composition """ # work in progress def injector(): """ """ #%% Trajectory tools def solve_traj(): """ """ # ....... #%% Astrodynamics tools def state_to_orbital_elem(r,v): """ Convert state vector to orbital elements inputs: r : coordinate vector (m) v : velocity vector (m.s-1) output: {'a','e','i','capital_omega','small_omega','nu'} where: a : semi_major axis (m) e : eccentricity (adimentional) i : inclination (deg) capital_omega : argument of the ascending node (deg) small_omega : argument of the perigee (deg) nu : true anomaly (deg) NB: - An internal quantity (small) is considered for the case of undefined classical orbital elements. When an element is undefined it is set to zero and other(s) following it are measured according to this convention - The function stops when r and v are aligned - This function is widely inspired by the Matlab function written by Prof. Josep Masdemart from UPC """ orbel={'a':0,'e':0,'i':0,'capital_omega':0,'small_omega':0,'nu':0} mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2) small = 1e-10 h = np.cross(r,v) # angular momentum rm = LA.norm(r) # module of the coordinate vector vm = LA.norm(v) # module of the velocity vector hm = LA.norm(h) # module of the angular momentum vector if hm < small: print('rectilinar or collision trajectory') return n = np.array([-h[1],h[0],0]) nm = LA.norm(n) if nm < small: n = np.array([1,0,0]) nm = 0 else: n = n/nm nm = 1 dotrv = np.dot(r,v) e = ((vm*vm-mu/rm)*r-dotrv*v)/mu em = LA.norm(e) if em < small: e = n em = 0 else: e = e/em orbel['e'] = em orbel['a'] = (hm*hm/mu)/(1-em*em) aux = h[2]/hm if np.abs(aux) > 1: aux = np.sign(aux) # all the angles are first computed in rad orbel['i'] = np.arccos(aux) if nm > small: aux = n[0] if np.abs(aux) > 1: aux = np.sign(aux) orbel['capital_omega'] = np.arccos(aux) if n[1] < 0: orbel['capital_omega'] = 2*np.pi-orbel['capital_omega'] if em > small: aux = np.dot(e,n) if np.abs(aux) > 1: aux = np.sign(aux) orbel['small_omega'] = np.arccos(aux) if e[2] < 0: orbel['small_omega'] = 2*np.pi-orbel['small_omega'] aux = np.dot(e,r)/rm if np.abs(aux) > 1: aux = np.sign(aux) orbel['nu'] = np.arccos(aux) if em > small: if dotrv < 0: orbel['nu'] = 2*np.pi-orbel['nu'] else: ha = np.cross(e,r) if np.dot(h,ha) < 0: orbel['nu'] = 2*np.pi-orbel['nu'] # convert the angles to deg orbel['i'] = np.rad2deg(orbel['i']) orbel['capital_omega'] = np.rad2deg(orbel['capital_omega']) orbel['small_omega'] = np.rad2deg(orbel['small_omega']) orbel['nu'] = np.rad2deg(orbel['nu']) return orbel def propagate_orbit(): """ """ #%% Thermodynamic cycles # the following functions provide the pressure, the mixuture composition # and the initial temperature in the combustion chambre and properties of some # components of the cycle. def pressure_fed(R,L,pi,Ltube): """ input p_fp: pressure of the fed pressure tank p_co: pressure of the oxizider tank p_ca: pressure of the fuel tank d_fp: denisty of the fed pressure fluid d_ox: density of the oxidizer d_fu: density of the fuel rajouter valves et ou et lesquelles ?? définir tank He définir pipe He -> heat exchanger définir pipe heat exchanger -> comb définir pipe heat exchanger -> carb définir tank comb définir tank carb définir pipe comb -> output définir pipe carb -> output """ #vi = 0 t = 0 dt = 1e-1 tank_Hein = [phe,vhe,0,dhe,rhe,0,rshe,lhe,che] tank_He = tank(tank_Hein+[dt]) tank_O2in = [po,vo,do,co,ro,reo] tank_O2in = tank(p_fp,0,0,d_fp,) tank_CH4in = [pch,vch,dch,cch,rch,rech] tank_CH4 = tank(p_fp,0,0,d_fp,) pipeHeHe = pipes(p,v,d,mu,r,l,h) tank_ while t < T tank_He = tank(p_fp,0,0,) t += dt def tank(p,v,d,dc,r,re,rs,l,c,dt): g=9.81 pi=3.1415926 """ input p : pression d'entrée v : vitesse d'entrée d : densité d'entrée dc : densité du contenant r : rayon re : rayon d'entrée rs : rayon de sortie l : longueur c : capacité de remplissage dt : pas de temps output ps : pression sortie vs : vitesse sortie dc : densité cs : capacité remplissage """ if d_c >= 0 pe = p+((d*v^2)/2*(1-(re^2/r^2))) ve = v*(re^2/r^2) pm = pe*(dc/d)+((dc*(ve^2)/2)*(1-(d/dc)^2))+dc*l*g*(1-2*c) vm = ve*(d/dc) else pm = p - (d*g*h) - (8*mu*l*(1-c)*v/(r^2)) vm = v ps = pm+((dc*vm^2)/2*(1-(r^2/rs^2))) vs = vm*(r^2/rs^2) cs = c-((vs*dc*dt)/(pi*r^2*l)) return ps,vs,dc,cs def pipes(p,v,d,mu,r,l,h): g = 9.81 """ input p : pression d'entrée v : vitesse d'entrée d : densité d'entrée mu : viscosité r : rayon l : longueur h : longueur total projeté sur l'axe vertical output ps : pression sortie vs : vitesse sortie """ pe = p + (d*g*h) - (8*mu*l*v/(r^2)) ve = v return ps,vs def heat_exchanger(p,v,d,mu,r,l,L): g = 9.81 """ input p : pression d'entrée v : vitesse d'entrée d : densité d'entrée mu: viscosité dynamique r : rayon l : longueur h : longueur total projeté sur l'axe vertical L : chaleur latente de vaporisation """ pe = p + (d*g*h) - (8*mu*l*v/(r^2)) ve = v #heat_exchanger ps = p - (L*d*S*l) vs = ve return ps,vs def valves(x,p,v,dp): pi = 3.1415926 """ input x: état imposé d'ouverture de la valve p: pression d'entrée v: vitesse d'entrée dp: pressure drop output xs: état actuel d'ouverture de la valve ps: pression de sortie de valves vs: vitesse d'entrée """ ps = p - dp #f(p,v,...) vs = v xs = x return xs,vs,ps def electric_turbo(): """ """ def closed_expander(): """ """ def dual_expander(): """ """ def open_expander(): """ """ def tap_off(): """ """ def gas_generator(): """ """ def ox_rich_staged_combust(): """ """ def full_flow_staged_combust(): """ """ def inv_ratioareasec_mach(rsa,gamma,subsup) """ inputs rsa : ratio des sections des aires gamma : coefficient de Laplace subsup : recherche d'un Mach sub ou supersonic (0|1) Ce dernier input est lier à l'initialisation de la méthode. """ np.float128 gamma=1.4 epsilon=1e-10 if subsup=1 : x=5 else: x=0.2 iter=0 while x > epsilon and iter < 20: iter += 1 y=np.sqrt((2/(gamma+1)*(1+(gamma-1)*(x**2)/2))**((gamma+1)/(gamma-1))/x**2)-rsa xp=x+epsilon yp=np.sqrt((2/(gamma+1)*(1+(gamma-1)*(xp**2)/2))**((gamma+1)/(gamma-1))/xp**2)-rsa x=x-epsilon*(y/(yp-y)) #pour afficher chaque itération, décommenter #print(x,y,iter) return x ======= import numpy as np from numpy import linalg as LA def initial_vel(latitude,retrograde=False): """ Initial velocity (m.s-1) depending on the launch point latitude (°) """ T_side = 86164.10 # sidereal day (s) R_T = 6378e3 # earth raduis (m) initial_velocity = (2*np.pi/T_side)*R_T*np.cos(np.deg2rad(np.abs(latitude))) if retrograde: return -initial_velocity else: return initial_velocity def final_target_vel(h, a): """ inputs: h: injection altitude (m) a: orbit semi major axis (m) output: orbit absolute injection velocity (m.s-1) """ mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2) R_T = 6378e3 # earth raduis (m) return np.sqrt(2*mu*((1/(R_T+h))-(1/(2*a)))) def staging_optim(isp1,isp2,eps1,eps2,Vc,Mp): """ Optimise mass loading in the case of two stage rocket inputs: isp1, isp2 : specific impulse of the 1st and 2nd stage (s) eps1, eps2 : structural coefficient of the 1st and 2nd stage (adimentional) Vc : target critical velocity (m.s-1) Mp : payload mass (Kg) outputs: M1 : mass of the 1st stage M2 : mass of the 2nd stage NB: structure mass structural coefficient = -------------------------------- structure mass + propellant mass ctritical velocity = final velocity - initial velocity + losses """ # Notation : # # Mt = total mass # q1 = M1/Mt # q2 = M2/Mt # qp = Mp/Mt # omega = qp + q2 g0 = 9.81 # standard gravity (m.s-2) omega_tab=np.linspace(0.001,0.9,100) # array of omega to test qp_tab=[] for omega in omega_tab: q2 = omega+eps1*(1-omega) q2 = np.power(q2,isp1/isp2)*np.exp(Vc/(g0*isp2)) q2 = 1/q2 q2 = 1-q2 q2 *= omega/(1-eps2) qp = omega-q2 qp_tab.append(qp) qp=max(qp_tab) omega=omega_tab[qp_tab.index(qp)] q2=omega-qp q1=1-q2-qp M1=(Mp/qp)*q1 M2=(Mp/qp)*q2 return M1, M2 def standard_atmosphere(altitude): """ Compute atmospheric properties according to the International Standaard Atmosphere (ISA) model input: altitude (m) outputs: temperature (K) pressure (Pa) density (Kg.m-3) """ assert altitude < 84852 , str(altitude)+' is higher than max altitude of the ISA model' g0 = 9.81 # standard gravity (m.s-2) r = 287 # specific gas constant (J.Kg-1.K-1) # ISA model h0_tab = np.array([0, 11e3, 20e3, 32e3, 47e3, 51e3, 71e3, 84852]) # base altitude (m) a_tab = np.array([-6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2.0e-3]) # lapse rate (K.m-1) # the base properties of each range are tabulated to reduce computational cost T0_tab = np.array([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65]) # base temperature (K) P0_tab = np.array([101325, 22632, 5474.9, 868.02, 110.91, 66.939, 3.9564]) # base pressure (Pa) index = 0 while altitude > h0_tab[index+1]: index += 1 a = a_tab[index] T0 = T0_tab[index] P0 = P0_tab[index] T = T0+a*(altitude-h0_tab[index]) # Toussaint equation if T == T0 : # isohermal region P = P0*np.exp(-g0*(altitude-h0_tab[index])/(r*T)) else : # gradient region P = P0*np.power(T/T0,-g0/(a*r)) rho = P/(r*T) return T,P,rho def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r,tolerance=1e-5): """ Compute the thrust assuming a started nozzle inputs: Pt : total pressure in the combustion chambre (Pa) Tt : total temperature in the combustion chambre (K) Pamb : ambient pressure (Pa) Ath : throat section area (m2) Ae : exit section area (m2) gamma : ratio of specific heats (adimentional) r : specific gas constant (J.Kg-1.K-1) tolerance : outputs: Thrust (N) Specific impulse (s) N.B : Thrust = mass flow rate * exit velocity + (exit pressure - ambiant pressure) * exit area """ assert gamma > 1, 'gamma must be > 1' g0 = 9.81 # standard gravity (m.s-2) A_ratio = lambda M,gamma: (1/M)*np.power((2/(gamma+1))*(1+((gamma-1)/2)*M*M),(gamma+1)/(2*gamma-2)) # compute the nozzle exit section area corresponding to the adapted regime Madapted = np.sqrt((2/(gamma-1))*(np.power(Pamb/Pt,(1-gamma)/gamma)-1)) Aadapted = Ath*A_ratio(Madapted,gamma) if Aadapted > Ae: # the nozzle is under-expanded # compute the exit Mach number iteratively using dichotomy algorithm A_ratio_target = Ae/Ath Me_upper = Madapted Me_lower = 1 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error > 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure elif Aadapted < Ae: # the nozzle is over-expanded Mshock = 1 step = 0.1 A_ratio_target = Ae/Ath Me_upper = 1 Me_lower = 0 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error < 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) while (Pe-Pamb)*1e-5 > tolerance: Mshock_old = Mshock Mshock += step Pe_old = Pe Pt_ratio = np.power(((gamma+1)*Mshock*Mshock)/((gamma-1)*Mshock*Mshock+2),gamma/(gamma-1))*np.power((gamma+1)/(2*gamma*Mshock*Mshock-gamma+1),1/(gamma-1)) Pi = Pt*Pt_ratio Ai = Ath/Pt_ratio A_ratio_target = Ae/Ai Me_upper = 1 Me_lower = 0 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error < 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) Pe = Pi*np.power(Tt/Te,gamma/(1-gamma)) if (Pe > Pamb) and (A_ratio(Mshock,gamma) > Ae/Ath): # the normal shock is outside the divergent and all the flow inside is supersonic A_ratio_target = Ae/Ath Me_upper = Mshock Me_lower = 1 Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) while np.abs(error) > tolerance: if error > 0: Me_lower = Me else: Me_upper = Me Me = (Me_upper+Me_lower)/2 error = A_ratio_target-A_ratio(Me,gamma) Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure break if Pe < Pamb: Pe = Pe_old Mshock = Mshock_old step /= 2 else: # the nozzle is adapted Me = Madapted Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature Pe = Pamb Ve = Me*np.sqrt(gamma*r*Te) # exit velocity # maximum mass flow rate in the case of a started nozzle Q = np.power(2/(gamma+1),(gamma+1)/(2*(gamma-1)))*np.sqrt(gamma/(Tt*r))*Pt*Ath Thrust = Q*Ve+(Pe-Pamb)*Ae isp = Thrust/(Q*g0) return Thrust,isp def adiabatic_flame(mech,T,P,X): """ Compute adiatabic flame temperature inputs: mech : chemical mechanism in cti format T : initial temperature (K) P : pressure (Pa) X : initial chemical composition """ # work in progress def injector(): """ """ #%% Trajectory tools def solve_traj(): """ """ # ....... #%% Astrodynamics tools def state_to_orbital_elem(r,v): """ Convert state vector to orbital elements inputs: r : coordinate vector (m) v : velocity vector (m.s-1) output: {'a','e','i','capital_omega','small_omega','nu'} where: a : semi_major axis (m) e : eccentricity (adimentional) i : inclination (deg) capital_omega : argument of the ascending node (deg) small_omega : argument of the perigee (deg) nu : true anomaly (deg) NB: - An internal quantity (small) is considered for the case of undefined classical orbital elements. When an element is undefined it is set to zero and other(s) following it are measured according to this convention - The function stops when r and v are aligned - This function is widely inspired by the Matlab function written by Prof. Josep Masdemart from UPC """ orbel={'a':0,'e':0,'i':0,'capital_omega':0,'small_omega':0,'nu':0} mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2) small = 1e-10 h = np.cross(r,v) # angular momentum rm = LA.norm(r) # module of the coordinate vector vm = LA.norm(v) # module of the velocity vector hm = LA.norm(h) # module of the angular momentum vector if hm < small: print('rectilinar or collision trajectory') return n = np.array([-h[1],h[0],0]) nm = LA.norm(n) if nm < small: n = np.array([1,0,0]) nm = 0 else: n = n/nm nm = 1 dotrv = np.dot(r,v) e = ((vm*vm-mu/rm)*r-dotrv*v)/mu em = LA.norm(e) if em < small: e = n em = 0 else: e = e/em orbel['e'] = em orbel['a'] = (hm*hm/mu)/(1-em*em) aux = h[2]/hm if np.abs(aux) > 1: aux = np.sign(aux) # all the angles are first computed in rad orbel['i'] = np.arccos(aux) if nm > small: aux = n[0] if np.abs(aux) > 1: aux = np.sign(aux) orbel['capital_omega'] = np.arccos(aux) if n[1] < 0: orbel['capital_omega'] = 2*np.pi-orbel['capital_omega'] if em > small: aux = np.dot(e,n) if np.abs(aux) > 1: aux = np.sign(aux) orbel['small_omega'] = np.arccos(aux) if e[2] < 0: orbel['small_omega'] = 2*np.pi-orbel['small_omega'] aux = np.dot(e,r)/rm if np.abs(aux) > 1: aux = np.sign(aux) orbel['nu'] = np.arccos(aux) if em > small: if dotrv < 0: orbel['nu'] = 2*np.pi-orbel['nu'] else: ha = np.cross(e,r) if np.dot(h,ha) < 0: orbel['nu'] = 2*np.pi-orbel['nu'] # convert the angles to deg orbel['i'] = np.rad2deg(orbel['i']) orbel['capital_omega'] = np.rad2deg(orbel['capital_omega']) orbel['small_omega'] = np.rad2deg(orbel['small_omega']) orbel['nu'] = np.rad2deg(orbel['nu']) return orbel def propagate_orbit(): """ """ #%% Thermodynamic cycles # the following functions provide the pressure, the mixuture composition # and the initial temperature in the combustion chambre and properties of some # components of the cycle. def pressure_fed(R,L,pi,Ltube): """ input p_fp: pressure of the fed pressure tank p_co: pressure of the oxizider tank p_ca: pressure of the fuel tank d_fp: denisty of the fed pressure fluid d_ox: density of the oxidizer d_fu: density of the fuel rajouter valves et ou et lesquelles ?? définir tank He définir pipe He -> heat exchanger définir pipe heat exchanger -> comb définir pipe heat exchanger -> carb définir tank comb définir tank carb définir pipe comb -> output définir pipe carb -> output """ #vi = 0 t = 0 dt = 1e-1 tank_Hein = [phe,vhe,0,dhe,rhe,0,rshe,lhe,che] tank_He = tank(tank_Hein+[dt]) tank_O2in = [po,vo,do,co,ro,reo] tank_O2in = tank(p_fp,0,0,d_fp,) tank_CH4in = [pch,vch,dch,cch,rch,rech] tank_CH4 = tank(p_fp,0,0,d_fp,) pipeHeHe = pipes(p,v,d,mu,r,l,h) tank_ while t < T tank_He = tank(p_fp,0,0,) t += dt def tank(p,v,d,dc,r,re,rs,l,c,dt): g=9.81 pi=3.1415926 """ input p : pression d'entrée v : vitesse d'entrée d : densité d'entrée dc : densité du contenant r : rayon re : rayon d'entrée rs : rayon de sortie l : longueur c : capacité de remplissage dt : pas de temps output ps : pression sortie vs : vitesse sortie dc : densité cs : capacité remplissage """ if d_c >= 0 pe = p+((d*v^2)/2*(1-(re^2/r^2))) ve = v*(re^2/r^2) pm = pe*(dc/d)+((dc*(ve^2)/2)*(1-(d/dc)^2))+dc*l*g*(1-2*c) vm = ve*(d/dc) else pm = p - (d*g*h) - (8*mu*l*(1-c)*v/(r^2)) vm = v ps = pm+((dc*vm^2)/2*(1-(r^2/rs^2))) vs = vm*(r^2/rs^2) cs = c-((vs*dc*dt)/(pi*r^2*l)) return ps,vs,dc,cs def pipes(p,v,d,mu,r,l,h): g = 9.81 """ input p : pression d'entrée v : vitesse d'entrée d : densité d'entrée mu : viscosité r : rayon l : longueur h : longueur total projeté sur l'axe vertical output ps : pression sortie vs : vitesse sortie """ pe = p + (d*g*h) - (8*mu*l*v/(r^2)) ve = v return ps,vs def heat_exchanger(p,v,d,mu,r,l,L): g = 9.81 """ input p : pression d'entrée v : vitesse d'entrée d : densité d'entrée mu: viscosité dynamique r : rayon l : longueur h : longueur total projeté sur l'axe vertical L : chaleur latente de vaporisation """ pe = p + (d*g*h) - (8*mu*l*v/(r^2)) ve = v #heat_exchanger ps = p - (L*d*S*l) vs = ve return ps,vs def valves(x,p,v,dp): pi = 3.1415926 """ input x: état imposé d'ouverture de la valve p: pression d'entrée v: vitesse d'entrée dp: pressure drop output xs: état actuel d'ouverture de la valve ps: pression de sortie de valves vs: vitesse d'entrée """ ps = p - dp #f(p,v,...) vs = v xs = x return xs,vs,ps def combustion_chamber(dth, gamma): """ dth: Diameter of throat """ # Determiner Acc/At: r_a = (8*(dth)**(-0.6))+1.25 # Déterminer Mcc Mcc = inv_ratioareasec_mach(r_a, gamma) # Déterminer la section de sortie Ath = pi*(dth**2)/4 Acc = r_a*Ath dcc = 2*sqrt(Acc/pi) # Déterminer la température de sortie Tcc = 4000 return Tcc, Pcc, Acc def electric_turbo(): """ """ def closed_expander(): """ """ def dual_expander(): """ """ def open_expander(): """ """ def tap_off(): """ """ def gas_generator(): """ """ def ox_rich_staged_combust(): """ """ def full_flow_staged_combust(): """ """ def inv_ratioareasec_mach(rsa,gamma,subsup) """ inputs rsa : ratio des sections des aires gamma : coefficient de Laplace subsup : recherche d'un Mach sub ou supersonic (0|1) Ce dernier input est lier à l'initialisation de la méthode. """ np.float128 gamma=1.4 epsilon=1e-10 if subsup=1 : x=5 else: x=0.2 iter=0 while x > epsilon and iter < 20: iter += 1 y=np.sqrt((2/(gamma+1)*(1+(gamma-1)*(x**2)/2))**((gamma+1)/(gamma-1))/x**2)-rsa xp=x+epsilon yp=np.sqrt((2/(gamma+1)*(1+(gamma-1)*(xp**2)/2))**((gamma+1)/(gamma-1))/xp**2)-rsa x=x-epsilon*(y/(yp-y)) #pour afficher chaque itération, décommenter #print(x,y,iter) return x >>>>>>> 1c1d15bbcbe958e44641336266bad5cfe11f9079