## GOX-kerosene sim
#@ Author Juha Nieminen

#import sys
#sys.path.insert(0, '/Users/juhanieminen/Documents/adamrocket')

import RocketComponents as rc
from physical_constants import poise, inches, Runiv, gallons, lbm, \
    gearth, atm, psi, lbf
from numpy import pi, linspace, cos, radians, sqrt, exp, log, array, full, ceil
from scipy import optimize as opt
import matplotlib.pyplot as plt
from matplotlib import collections as mc
import Flows1D as flows


#DESIGN VARIABLES____________________________________________________________________________________

# nominal parameters
Preg_N2          = 1300*psi         # regulated N2 outlet pressure [Pa]

mdot_fuel_nom   = 0.2               # This is only for cooling jacket pressure drop purposes [kg/s]
Pdrop_jacket_nom= 1*psi             # Cooling jacket pressure drop at mdot_nominal [Pa]
OF_nom          = 2.25              # Oxidizer-to-fuel ratio. This has only effect on initial guesses during solving

# Pressurant tank dimensions
Vprestank       = 0.053             # N2 pressurant tank volume [m3]

# Propellant tank dimensions
Vfueltank       = 4*gallons         # fuel tank volume [m3]
Voxtank         = 4*0.053           # ox tank volume [m3]

# Tubing
d_presfuel_tube = 1.0*inches        # pressurant tank -> fuel tank tube diameter [m]
L_presfuel_tube = 0.5               # pressurant tank -> fuel tank tube length [m]

d_oxtube        = 0.87*inches       # ox tank -> manifold tube diameter [m]
L_oxtube        = 2.4               # ox tank -> manifold tube length [m]
d_fueltube      = 0.87*inches       # fuel tank -> manifold tube diameter [m]
L_fueltube      = 3.0               # fuel tank -> manifold tube length [m]

roughness       = 0.005             # epsilon/diameter, dimensionless

# Valves
Cv_ox_check     = 4.7               # oxidizer check valve flow coefficient, dimensionless
Pcrack_ox_check = 10*psi            # oxidizer check valve opening pressure [Pa]
Cv_pres_check   = 1.8               # nitrogen check valve flow coefficient, dimensionless
Pcrack_pres_check  = 0.33*psi       # nitrogen check valve opening pressure [Pa]
Cv_pres_valve   = 8.8               # nitrogen solenoid valve flow coefficient, dimensionless
Cv_ox_valve     = 8.8               # oxidizer solenoid valve flow coefficient, dimensionless 
Cv_fuel_valve   = 8.8               # fuel solenoid valve flow coefficient, dimensionless

# Injector

cd_oxInjector   = 0.767                                                 # orifice discharge coefficient
diameter_oxInjectorHoles = 2.54e-3 #number xx drill                     # ox orifice diameter [m]
#length_oxHole   = 0.005                                                # ox orifice length [m]
numOxInjectorHoles = 24                                                 # number of ox orifices in the injector
area_oxInjector = numOxInjectorHoles*pi*diameter_oxInjectorHoles**2/4   # total ox flow area [m2]

cd_fuelInjector = 0.767                                                 # orifice discharge coefficient
diameter_fuelInjectorHoles = 0.508e-3 #number xx drill                      # fuel orifice diameter [m]
numFuelHoles    = 59                                                    # number of fuel orifices in the injector
area_fuelInjector = numFuelHoles*pi*diameter_fuelInjectorHoles**2/4     # total fuel flow area [m2]


# Define initial/nominal conditions in the chamber (obtained from CEA code assuming OFratio = 2.25)
TfireInit       = 293                                                   # initial flame temperature [K]
Pfire           = 1*atm                                                 # initial chamber pressure [Pa]
gammaFireInit   = 1.148                                                 # dimensionless
ga              = gammaFireInit
mbarFireInit    = 21.87                              # combustion products' initial molecular mass [kg/kmol]
RfireInit       = Runiv/mbarFireInit                 # combustion products' initial specific gas constant [J/kgK]
Pambient        = atm                                                   # ambient pressure [Pa]

# Nozzle and chamber

d_nozzleThroat  = 1.0*inches                                            # throat diameter [m]
A_nozzleThroat  = pi*d_nozzleThroat**2/4                                # throat area [m2] 
area_ratio      = 7.46                                                  # nozzle exit-to-throat area ratio
A_nozzleExit    = area_ratio*A_nozzleThroat                             # nozzle exit area [m2]
d_nozzleExit    = sqrt(4*A_nozzleExit/pi)                               # nozzle exit diameter [m]
            
Dchamber        = 0.08                                                  # chamber diameter [m]
Achamber        = pi*Dchamber**2/4                                      # chamber cross sectional area [m2]
Lchamber        = 0.14                                                  # chamber length [m]
Vchamber        = Achamber*Lchamber                                     # chamber volume [m3] 
Lstar           = Vchamber/A_nozzleThroat                               # chamber characteristic length [m]
Mc_nom          = flows.getIsentropicMs(A_nozzleThroat, Achamber, gammaFireInit)[0] # nominal chamber Mach number

print("throat diameter is", '%.1f'%(d_nozzleThroat*1000), 'mm')
print("exit diameter is", '%.1f'%(d_nozzleExit*1000), 'mm')
print("chamber volume is", '%.5f'%Vchamber, "m3")
print("chamber Lstar is", '%.2f'%Lstar, "m")
print("chamber Mach_nom is", '%.2f'%Mc_nom)



# INITIAL CONDITIONS____________________________________________________________________________________________

#Define initial conditions in the tanks

TfuelPresStart  = 293                   # Fuel pressurant (=nitrogen) temp [K]
FFfueltankStart = 0.9                   # Fuel tank fill fraction (Vfuel/Vtank)
PfuelPrestankStart  = 2640*psi - Preg_N2*Vfueltank*(1-FFfueltankStart)/Vprestank     # Fuel pressurant tank pressure once fueltank has been pressurized [Pa]

ToxStart        = 293                   # Oxidizer (GOX) temp [K] 
PoxtankStart    = 1600*psi              # Oxidizer tank pressure [Pa]
TfuelStart      = 293                   # Fuel temp [K]
PfueltankStart  = Preg_N2 -1*psi       # Fuel tank pressure [Pa] (-10psi helps convergence on first timestep)

# initialize propellants
nitrogen        = rc.NitrogenFluid()
GOX             = rc.GOXFluid()
kerosene        = rc.Kerosene()

#initialize nozzle and chamber
nozzle          = rc.ConvergingDivergingNozzle(A_nozzleExit, A_nozzleThroat)
mdot_init_noz   = nozzle.getmdot(gammaFireInit, GOX.R, Pfire, TfireInit, atm)
chamber         = rc.GOXKeroCombustionChamber(nozzle, Vchamber, TfireInit, ga, mbarFireInit, Pfire, atm, mdot_init_noz)

#initialize injector orifices
ox_orifice      = rc.GasOrifice(area_oxInjector, cd_oxInjector, GOX.gamma, GOX.R)
fuel_orifice    = rc.LiquidOrifice(area_fuelInjector, cd_fuelInjector )

#initialize pressurant tanks
fuelprestank    = rc.IdealgasTank(nitrogen, Vprestank, TfuelPresStart, PfuelPrestankStart)

#initialize propellant tanks
oxtank          = rc.IdealgasTank(GOX, Voxtank, ToxStart, PoxtankStart)
fueltank        = rc.LiquidPropellantTank(nitrogen, kerosene, Vfueltank, TfuelStart, TfuelPresStart,\
                    PfueltankStart, FFfueltankStart, Preg_N2)

#initialize pressure regulators
N2_regu         = rc.PressureRegulator(Preg_N2, nitrogen)

#initialize solenoids
fuelSole        = rc.IncompressibleFlowSolenoid( Cv_fuel_valve)
oxSole          = rc.CompressibleFlowSolenoid( Cv_ox_valve, GOX)
presSole        = rc.CompressibleFlowSolenoid( Cv_pres_valve, nitrogen)

#initialize check valves
ox_check        = rc.CompressibleFlowCheckValve( Cv_ox_check, Pcrack_ox_check, GOX)
pres_check      = rc.CompressibleFlowCheckValve( Cv_pres_check, Pcrack_pres_check, nitrogen)

#initialize tubing
ox_tube         = rc.RoughStraightCylindricalTube(d_oxtube, L_oxtube, roughness, True)
fuel_tube       = rc.RoughStraightCylindricalTube(d_fueltube, L_fueltube, roughness, True)
presfuel_tube   = rc.RoughStraightCylindricalTube(d_presfuel_tube, L_presfuel_tube, roughness, True)

#initialize cooling jacket
jacket          = rc.CoolingJacket(mdot_fuel_nom, Pdrop_jacket_nom)

#initialize arrays for various data time histories
T_chamber       = [chamber.T]                       # combustion chamber temperature [K]
Pchamber        = [chamber.get_P_inlet()]           # combustion chamber pressure [Pa]
Pexit           = [nozzle.getPe(Pchamber[0], gammaFireInit, Pambient)] # nozzle exit pressure [Pa]
Mexit           = [nozzle.getMe(Pchamber[0], gammaFireInit, Pambient)] # nozzle exit Mach number
cmass           = [chamber.m]                       # resident propellant mass in combustion chamber [kg]
mdot_nozzle     = [nozzle.getmdot(gammaFireInit, RfireInit, chamber.get_P_inlet(), chamber.T, chamber.Pa)] # mass flow out of the nozzle [kg/s]
Poxtank         = [oxtank.getPtank()]               # ox tank pressure [Pa]
Toxtank         = [oxtank.getTtank()]               # ox tank temperature [K]
mox             = [oxtank.getM()]                   # oxidizer mass in tank [kg]

Pfueltank       = [fueltank.getPtank()]             # fuel tank pressure [Pa]
Tfueltank       = [fueltank.getTpres()]             # pressurant temperature in fuel tank[K]
mPresFueltank   = [fueltank.getMpres()]             # pressurant mass in fuel tank [kg]
mfuel           = [fueltank.getMprop()]             # fuel mass in tank [kg]
FFfueltank      = [fueltank.getFF()]                # fuel tank fill fraction defined as Vfuel/(Vfueltank)

TfuelPres       = [fuelprestank.getTtank()]         # temperature in fuel pressurant tank [K]
PfuelPres       = [fuelprestank.getPtank()]         # pressure in fuel pressurant tank [Pa]
mfuelPres       = [fuelprestank.getM()]             # pressurant mass in fuel pressurant tank [Pa]

time            = [0]                               # time array [s]

mdot_ox         = [0]                               # ox mass flow out of the tank [kg/s]
P1ox            = [0]                               # ox tank presssure [Pa]
P2ox            = [0]                               # ox check valve outlet pressure [Pa]      
P3ox            = [0]                               # ox flow solenoid outlet pressure [Pa]
P4ox            = [0]                               # ox injector inlet pressure [Pa]

T1ox            = [0]                               # ox tank temp [K]
T2ox            = [0]                               # ox check valve outlet temp [K]
T3ox            = [0]                               # ox flow solenoid outlet temp [K]
T4ox            = [0]                               # ox injector inlet temp [K]

mdot_fuel       = [0]                               # fuel mass flow out of the tank [kg/s]
rooFuel         = fueltank.propellant.density       # fuel density, assumed constant [kg/m3]
P1fuel          = [0]                               # fuel tank presssure [Pa]
P2fuel          = [0]                               # fuel solenoid outlet pressure [Pa]      
P3fuel          = [0]                               # fuel cooling jacket inlet pressure [Pa]
P4fuel          = [0]                               # fuel injector inlet pressure [Pa]

mdot_fuel_pres  = [0]                               # fuel pressurant mass flow rate [kg/s]
P3pres          = [0]                               # pressurant pressure at check valve outlet [kg/s]
P4pres          = [0]                               # pressurant pressure at solenoid valve outlet [kg/s]

mTotal          = [0]                               # propellant mass in the system [kg]
mprs            = [mfuelPres[0]+mPresFueltank[0]]   # pressurant mass in the system [kg]
OFratio         = [0]                               # oxidizer to fuel mass flow ratio 
Isp             = [0]                               # specific impulse [s]
Thrust          = [nozzle.getThrust(chamber.get_P_inlet(), Pambient, gammaFireInit) ] # rocket thrust [N]

#SIMULATE_______________________________________________________________________________________________________
# using orifices as follows: ejecting GOX from manifold to chamber, fuel liq-to-liq from manifold to chamber

print("")
print("STARTING SIM...")
print("")
print("mOxStart is", '%.2f'%mox[0], "kg")
print("mKerostart is", mfuel[0], "kg")
print("mN2start in N2 tank is", '%.2f'%mfuelPres[0], "kg")
print("mN2start in fuel tank is", '%.2f'%(fueltank.getMpres()), "kg")

# The first step is to solve oxidizer and fuel mass flow rates from the tank to combustion chamber.

# definitions:

# P1ox              = GOX tank pressure
# P2ox              = check valve outlet pressure
# P3ox              = ox valve outlet pressure
# P4ox              = injector inlet, pressure 
# (P1ox-P2ox)       = ox check valve pressure drop,             eq 1                 
# (P2ox-P3ox)       = ox flow solenoid pressure drop,           eq 2            
# (P3ox-P4ox)       = ox tubing pressure drop,                  eq 3
# (P4ox-Pchamber)   = ox injector pressure drop,                eq 4

# P1pres            = Nitrogen tank pressure
# P2pres            = Regulation pressure
# P3pres            = Check valve outlet pressure
# P4pres            = Nitrogen solenoid outlet
# P5pres            = Nitrogen tubing outlet = fuel tank pressure
# (P2pres-P3pres)   = Nitrogen check valve pressure drop
# (P3pres-P4pres)   = Nitrogen solenoid valve pressure drop
# (P4pres-P5pres)   = Nitrogen tubing pressure drop

# P1fuel            = fuel tank pressure
# P2fuel            = fuel valve outlet pressure
# P3fuel            = cooling jacket inlet pressure
# P4fuel            = injector inlet pressure
# (P1fuel-P2fuel)   = fuel valve pressure drop,                 eq1
# (P2fuel-P3fuel)   = fuel tubing pressure drop,                eq2
# (P3fuel-P4fuel)   = cooling jacket pressure drop,             eq3
# (P4fuel-Pchamber) = injector pressure drop,                   eq4

# In the case of oxidizer, P1 and Pchamber are known, so one must solve for P2, P3, and P4. Fourth unknown is the mass flow rate. The four equations are check valve/solenoid/tubing/injector pressure drops. These equations are defined in oxfunks method below, and underlying physics are in RocketComponents.py under their respective classes.

# With pressurant, P2 (regulation pressure) and P5 (fuel tank pressure) are known, so one must solve for P3 and P4. The third unknown is pressurant mass flow rate. Equations to be solved are pressure drops over the check valve, solenoid valve, and the tubing.

# With fuel P1 and Pchamber are known, so one must solve for P2, P3, and P4. Fourth unknown is mass flow rate.

# fsolve requires sensible initial guesses for all unknowns. They are established by guessing the mass flow rate, because all other pressures trickle down from that.


timestep_small      = 1e-5      # seconds, used during initial transient
timestep_nom        = 1e-4      # seconds, used after 0.01 seconds of simulation time
t_transient         = 0.01      # seconds, estimated time of initial transient

t_simulation        = 3    # seconds

if t_simulation <= t_transient:
    simsteps    = int(ceil(t_simulation/timestep_small))
else:
    simsteps    = int(ceil( t_transient/timestep_small + (t_simulation-t_transient)/timestep_nom ))

print("Sim time is", t_simulation, "s, number of simsteps is", simsteps)

i=0
for i in range(0, simsteps):
    if time[i] < t_transient:
        timestep    = timestep_small                   # use shorter timestep during initial transient
    else: timestep    = timestep_nom                   # proceed with nominal timestep
#while True:
    print("i=", i) 
    P1ox            = Poxtank[i]
    P1fuel          = Pfueltank[i]
    Pchamb          = Pchamber[i]
    mu_ox           = GOX.getViscosity(P1ox, Toxtank[i])
    roo_ox          = GOX.getDensity(P1ox, Toxtank[i])
    Tox             = Toxtank[i]
    Tpres           = TfuelPres[i]
    mu_fuel         = kerosene.mu
    mu_N2_fuel      = nitrogen.getViscosity(Preg_N2, TfuelPres[i])
    roo_N2_fuel     = nitrogen.getDensity(Preg_N2, TfuelPres[i])

    if i==0:    # First guesses. Based on choked flow at ox injector (multiplied by 0.7 to adjust for better convergence)
        
        
        mdot_injector_choked    = ox_orifice.getMdot(P1ox, Pfire, Tox)
        '''
        mdot_checkvalve_choked  = ox_check.getMdot(P1ox, Pfire, GOX.roo_std, roo_ox, Tox)
        
        if mdot_injector_choked >= mdot_checkvalve_choked:      #check valve is choking
            print("check valve is initially choking")
            mdot_ox_guess   = mdot_checkvalve_choked
            print("mdot_ox_guess is", mdot_ox_guess)
            P4ox_guess      = ox_orifice.getUpstreamPressure(Pchamb, Tox, mdot_ox_guess)
            P3ox_guess      = P4ox_guess + ox_tube.getPressureDrop(mdot_ox_guess, mu_ox, roo_ox)
            P2ox_guess      = P3ox_guess + oxSole.getPressureDrop(mdot_ox_guess, P2ox_guess, roo_ox)
       
        else:
        ''' 
        mdot_ox_guess   = mdot_injector_choked *0.7
        P2ox_guess      = P1ox - ox_check.getPressureDrop(mdot_ox_guess, P1ox, GOX.roo_std, roo_ox, Tox)
        P3ox_guess      = P2ox_guess - oxSole.getPressureDrop(mdot_ox_guess, P2ox_guess, roo_ox)
        P4ox_guess      = P3ox_guess - ox_tube.getPressureDrop(mdot_ox_guess, mu_ox, roo_ox)
        
        print("mdot_ox_guess is", mdot_ox_guess)
        #print("P2ox is", P2ox/psi, "psi")
        #print("P3ox_guess is", P3ox_guess/psi, "psi")
        #print("P4ox_guess is", P4ox_guess/psi, "psi")
        #print("P5ox_guess is", P5ox_guess/psi, "psi")
        #print("P_chamber is", Pchamber[i]/psi, "psi")
        
        mdot_fuel_guess = mdot_ox_guess/OF_nom
        P2fuel_guess    = P1fuel - fuelSole.getPressureDrop(mdot_fuel_guess, rooFuel)
        P3fuel_guess    = P2fuel_guess - fuel_tube.getPressureDrop(mdot_fuel_guess, mu_fuel, rooFuel)
        P4fuel_guess    = P3fuel_guess - jacket.getPressureDrop(mdot_fuel_guess)
        
        mdot_pres_guess = mdot_fuel_guess*roo_N2_fuel/rooFuel #volumetric flowrates of fuel and pressurant are the same
        P3pres_guess    = Preg_N2 - pres_check.getPressureDrop(mdot_pres_guess, Preg_N2, nitrogen.roo_std, roo_N2_fuel, Tpres)
        P4pres_guess    = P3pres_guess - presSole.getPressureDrop(mdot_pres_guess, P3pres_guess, roo_N2_fuel)
        P5pres_guess    = P4pres_guess - presfuel_tube.getPressureDrop(mdot_pres_guess, mu_N2_fuel, roo_N2_fuel)
        
        #print("mdot_pres_guess is is", mdot_pres_guess, "kg/s")
        #print("P3pres_guess is is", P3pres_guess/psi, "psi")
        #print("P4pres_guess is is", P4pres_guess/psi, "psi")
        #print("P5pres_guess is is", P5pres_guess/psi, "psi")
        
        
        #print("mdot_fuel_guess is", mdot_fuel_guess)
        #print("P2fuel is", P2fuel/psi, "psi")
        #print("P3fuel_guess is is", P3fuel_guess/psi, "psi")
        #print("P4fuel_guess is is", P4fuel_guess/psi, "psi")
        #print("P5fuel_guess is is", P5fuel_guess/psi, "psi")
        #print("P_chamber is", Pchamber[i]/psi, "psi")
      
    else :      # guesses for further steps. Use values from previous timestep
        mdot_ox_guess   = mdot_ox[i-1] #ox_orifice.getMdot(Preg_ox, Pchamb, Tox)
        #P3ox_guess      = P2ox - oxSole.getPressureDrop(mdot_ox_guess, P2ox,roo_ox)
        #P4ox_guess      = P3ox_guess - ox_tube.getPressureDrop(mdot_ox_guess, mu_ox, roo_ox)
        P2ox_guess      = P2ox[i-1]
        P3ox_guess      = P3ox[i-1]
        P4ox_guess      = P4ox[i-1]
        #print("mdot_ox_guess is", mdot_ox_guess)
        #print("P2ox_guess is", P2ox_guess/psi, "psi")
        #print("P3ox_guess is", P3ox_guess/psi, "psi")
        #print("P4ox_guess is", P4ox_guess/psi, "psi")
        #print("P_chamber is", Pchamber[i]/psi, "psi")
        
        mdot_fuel_guess  = mdot_fuel[i-1] #mdot_ox_guess/OF_nom*1
        P2fuel_guess    = P2fuel[i-1]
        P3fuel_guess    = P3fuel[i-1]
        P4fuel_guess    = P4fuel[i-1]
        #print("P2fuel is", P2fuel/psi, "psi")
        #print("P3fuel_guess is is", P3fuel_guess/psi, "psi")
        #print("P4fuel_guess is is", P4fuel_guess/psi, "psi")
        #print("P_chamber is", Pchamber[i]/psi, "psi")
        
        mdot_pres_guess = mdot_fuel_pres[i-1]
        P3pres_guess    = P3pres[i-1]
        P4pres_guess    = P4pres[i-1]
        
       
    initial_ox_guesses  = [P2ox_guess, P3ox_guess, P4ox_guess, mdot_ox_guess]