Engineering_tools.py

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