Cea integration

.gitignore

@@ -5,3 +5,4 @@ env
 *.pyc
 *.pyc
 *.pyc
+*.png

README.md

@@ -1,3 +1,42 @@
+# To Do List
+small tasks:
+make everything numpy arrays
+make flow seperation calculator
+make volume calculator take take convergine 
+
+cooling calculator to do:
+make "geometry" take input "contour" (see if radius or diameter is needed)
+
+major tasks
+1. bartz equasion needs a complete overhall to test accuracy, viscosity is a constant for some reason
+2. make engines class to run multiple engines
+3. make prop feed class that has tank stats and stuff
+4. make rocket class that takes engine and prop feed classes as inharitence
+5. make engines class somhow able to take rocket classes too and iterate
+6. look into making rocketCEA run faster
+7. make engine able to run with lower contour resolution
+8. look into making solveMach run faster
+9. code thrust level so it can do a linear sweep that is O(n) rather than the current O(nlog(n)) method
+
+major things to do:
+make rocket simulator with rocketpy
+complete overhall of equations
+check output curves to other caluculators to verify data accuracy
+look into making CEA run faster somehow(either by writing my own CEA wrapper of modifying one of the existing ones)
+code thrust level so it can do a linear sweep that is O(n) rather than the current O(nlog(n)) method
+code a UI that is easy to distribute, intuitive to use, and has full functionallity
+add 1d regen cooling channel caluclator(in progress by emmett)
+add film cooling calculation capabilities
+calculate flow seperation
+calculate combustion instability modes
+
+eventually
+1. get ansyst o run through python
+2. get rocket trajectory and aero simulator
+3. cooling channel calculator
+4. film cooling calculator
+5. injegrate injector class
+
 # Rocket Configuration Visualizer
 
 <!-- Add information here explaining the visualizer -->
@@ -10,10 +49,31 @@ This code sizes a rocket engine contour based on design specifications and therm
 2. Pull a clone of this repository
 3. Create a python virtual enviorment with the following command in the root directory of the Rocket Configuration Visualizer. `python -m venv venv/`
 4. Finally install all the dependencys once in your virtual enviorment of python with `pip3 install -r requirements.txt` (do this in the Terminal of VS Code, not cmd)
+5. install rocketcea library. this is different per operating system but this  goes over how to do it(https://rocketcea.readthedocs.io/en/latest/quickstart.html) 
+note: this library can be rather difficult to install correctly on windows
+sometimes the interpreter won't be able to find rocketcea when using a virtual enviernment; if that happens you can either not use a venv or use the following code to manually locate rocketcea
+# import sys
+# sys.path.insert(0, r"C:\Users\<insert_user>\AppData\Local\Programs\Python\Python310\Lib\site-packages")
+the specific file location and names will be defferent depenging on the version and location of your python files
+
+## Mac install guide
+1. Install python (It can be probably be the newest one but there could be issues becuase the rocketcea website says it is compatibly with 3.10 and below) you can check python version in terminal with "python3 --version" or maybe "python --version"
+2. download github desktop (or prefered github tool)
+3. pull repository from "https://github.com/LiquidPropulsionGroup/RocketConfigurationVisualizer/tree/CEA-integration" make sure it is the CEA integration branch
+4. download vscode
+5. download python extension in vscode
+6. select python interpreter to use with ctl-g and then typing ">Python: select interpreter" (make sure the selected interpreter and the one being called by pip3 are the same)
+7. update operating system
+8. download xcode through app store
+9. (in terminal)download homebrew
+10. (in terminal)download gcc (make sure it is linked properly), can check for gfrotran with "which gfortran"
+11. (vscode command line)install all the dependencys with `pip3 install -r requirements.txt` NOTE: this is not in a virtual enviernment beucase it causes issues with rocketcea library accessing gfortran
+12. install rocketcea seperatly with 'pip3 install rocketcea', this fails in many cases so this is the website for refference (https://rocketcea.readthedocs.io/en/latest/quickstart.html)
+
 
 ## Basic Overview
 
-Chemistry class is instantiated per every imported chemistry data point from CEA. Rocket class generates a unique rocket contour and computes flow and heat properties upon instantiation. Mulitple Rocket objects represent different rocket sizes. There are inline comments where possible in the code. Further documentation coming soon.
+The engine class takes the input veriables creates the optimal engine contour for the given constraints and throttle range. To create this contour chemistryCEA is used. this class takes the raw data output, which is from the nasaCEA code that is run through the rocketcea library, and parces the data into veriables. durring this optimization of the contour, thrustLevel instances are created for the max and min thrusts. when thrustLevel is instantiated basic veriables are calculated. from there we go back to rocket where the contour is made into an array and this array is then given back to the thrust level instances in the heatCalcs method, which calculates the heat flux at every point int eh contour array using the bartz approximation and other simple operations. at this point all of the calculations have been done and you can view the data using the variablesDisplay and graphDisplay methods in engine. 
 
 ## Conventions
 
@@ -33,3 +93,93 @@ my_contourPoints is a function that generates x and y positions for critical poi
 
 This function interpolates the contour for numerical calculations. It uses lambda functions: the contour consists of straight lines and circles per Sutton. Nozzle is conical, with 80% bell optimization capability to be added in the future.
 
+## Math Review
+
+### inputs
+
+ox  # full propellant selection is availible at https://rocketcea.readthedocs.io/en/latest/propellants.html
+fuel
+customFuel = [
+    "Isopropanol70", 
+    """fuel C3H8O-2propanol C 3 H 8 O 1    wt%=70.0
+h,cal=-65133.0     t(k)=298.15   rho=0.786
+fuel water H 2.0 O 1.0  wt%=30.0
+h,cal=-68308.0  t(k)=298.15 rho,g/cc = 0.9998"""
+]
+pMaxCham    #max thrust chamber pressure in bar
+Mr  #propellant mixture ratio
+pAmbient
+pMinExitRatio
+#set veriables
+mdotMax  #max thrust mass flow rate
+filmCoolingPercent
+Lstar
+Dcham
+conv_angle
+div_angle
+wall_temp
+fuel_delta_t
+fuel_cp
+r1
+r2
+r3
+step
+nozzle_type
+### outputs
+
+
+
+### CEA and parser
+the thermochemistry part of this code is done through the CEArocket package which is a wraper for the NASA CEA FORTRAN code. the chemistry CEA
+#### CEA inputes
+fuel
+pCham
+Mr
+pAmbient or ea
+frozen or equilibrium
+
+#### CEA outputs
+
+
+### chamber nozzle contour
+points a,b,c,d,o,n,e from left to right where the exit is to the right are used to define the nozzle contour
+r1, r2, r3 are the radii ratios compared to the throat for the throat and chamber curvatures
+o stands for origin
+e for exit
+the nozzle can be either connical or a bell nozzle:
+for connical a 15 degree strait line is drawn from point n to the exit radius
+
+#### 80% bell nozzle calcs
+for the bell nozzle a parabola is drawn with the x distance being defined as 80% of the length of the respective connical nozzle. the paper below has detailes for the equations we used for making this.
+"Design and analysis of contour bell nozzle and comparison with dual bell nozzle"
+https://core.ac.uk/download/pdf/154060575.pdf
+
+### mach area
+
+( 2 / (self.thr.gam + 1) * ( 1 + (self.thr.gam - 1)/2 * mach**2 ))**((self.thr.gam+1)/(2*(self.thr.gam-1))) - mach * area_ratio
+https://web.mit.edu/16.unified/www/SPRING/propulsion/notes/node104.html
+IMPORTANT NOTE: area ratio might be wrong, might need to be inversed
+
+### bartz
+(0.026 / math.pow(d_throat, 0.2) * math.pow((p_chamber / c_star), 0.8) * math.pow((d_throat / d),
+            1.8) * c_p * math.pow(visc, 0.2) * math.pow((t_gas / t_boundary), (0.8 - 0.2 * 0.6)))
+
+h_g_arr[1,i] = self.bartz(self.thr.d, self.cham.p, self.Cstar, self.contour[1,i]*2, self.cham.cp*1000, 0.85452e-4, self.temp_arr[1,i], self.wall_temp)
+
+https://arc.aiaa.org/doi/pdf/10.2514/8.12572
+
+### adiabatic pressure, temperature, density
+
+t_stag = self.cham.t * (1 + ((gam-1)/2 * self.cham.mach**2))
+myreturn = t_stag * (1 + ((gam-1)/2 * mach**2))**(-1)
+
+p_stag = self.cham.p * (1 + ((gam-1)/2 * self.cham.mach**2))**(gam/(gam-1))
+myreturn = p_stag * (1 + ((gam-1)/2 * mach**2))**(-gam/(gam-1))
+
+d_stag = self.cham.rho * (1 + ((gam-1)/2 * self.cham.mach**2))**(1/(gam-1))
+myreturn = d_stag * (1 + ((gam-1)/2 * mach**2))**(-1/(gam-1))
+
+
+### other heat calcs
+
+### wall thermal calcs

big_engine.py

@@ -0,0 +1,46 @@
+import math
+from src_cea.engine import Engine
+
+#test engine 
+title = 'Big Engine'
+ox = 'LOX'        # full propellant selection is availible at https://rocketcea.readthedocs.io/en/latest/propellants.html
+fuel =  'RP1'# 'Isopropanol70''Ethanol''CH4'
+customFuel = [
+    "Isopropanol70", 
+    """fuel C3H8O-2propanol C 3 H 8 O 1    wt%=70.0
+h,cal=-65133.0     t(k)=298.15   rho=0.786
+fuel water H 2.0 O 1.0  wt%=30.0
+h,cal=-68308.0  t(k)=298.15 rho,g/cc = 0.9998"""
+]
+#print('{}'.format(customFuel[1]))
+pMaxCham = 125     #max thrust chamber pressure in bar
+Mr = 2.1 # propellant mixture ratio
+pAmbient = 1.01325 #bar
+#pMaxCham = 25*14.5     #max thrust chamber pressure in bar
+#pAmbient = 1.01325*14.5 #bar
+pMinExitRatio = [] #trottle exit pressure
+#set veriables
+mdotMax = 5        #max thrust mass flow rate
+filmCoolingPercent = 0.0
+Lstar = 1.02
+Dcham = 3.75* 0.0254 #in meters
+conv_angle = math.pi / 4 # rad, 45deg
+div_angle = math.pi / 12  # rad, 15deg
+wall_temp = 650 # K
+fuel_delta_t = 300 # K
+fuel_cp = 2010 # J/KgK
+r1 = 1
+r2 = 1
+r3 = 0.4
+step = 1e-3 #array resolution
+nozzle_type = 'bell80' #'conical'
+doContours = True
+
+test = Engine(title, fuel, ox, nozzle_type, Mr, pMaxCham, mdotMax, Lstar, Dcham, wall_temp, r1, r2, r3, conv_angle, fuel_delta_t, pMinExitRatio = pMinExitRatio, filmCoolingPercent = filmCoolingPercent, div_angle = div_angle, contourStep = step, customFuel = customFuel, frozen = 1, pAmbient = pAmbient, doContours = doContours)
+test.variablesDisplay(minthrust = False)
+#test.debugAndRawVariablesDisplay()
+test.graphDisplay(minthrust = False)
+#test.runTime()
+#test.fullCEAOut()
+#test.testRunCEAOut()
+#print(test.cea.get_Densities(Pc=pMaxCham, MR=Mr, eps = test.max.exit.aeat))

cea_test.py

@@ -1,7 +1,171 @@
 from rocketcea.cea_obj import CEA_Obj
+import rocketcea.cea_obj as cea_obj
+import rocketcea.py_cea as py_cea
+import time
+from CEA_Wrap import Fuel, Oxidizer, RocketProblem
+import platform
+
+#: RocketCEA wraps the NASA FORTRAN CEA code to calculate Isp, cstar, and Tcomb
+#: This object wraps the English unit version of CEA_Obj to enable desired user units.
+#: Same as CEA_Obj with standard units except, input and output units can be specified.
+#:  parameter             default             options
+#: isp_units            = 'sec',         # N-s/kg, m/s, km/s
+#: cstar_units          = 'ft/sec',      # m/s
+#: pressure_units       = 'psia',        # MPa, KPa, Pa, Bar, Atm, Torrs
+#: temperature_units    = 'degR',        # K, C, F
+#: sonic_velocity_units = 'ft/sec',      # m/s
+#: enthalpy_units       = 'BTU/lbm',     # J/g, kJ/kg, J/kg, kcal/kg, cal/g
+#: density_units        = 'lbm/cuft',    # g/cc, sg, kg/m^3
+#: specific_heat_units  = 'BTU/lbm degR' # kJ/kg-K, cal/g-C, J/kg-K (# note: cal/g K == BTU/lbm degR)
+#: viscosity_units      = 'millipoise'   # lbf-sec/sqin, lbf-sec/sqft, lbm/ft-sec, poise, centipoise
+#: thermal_cond_units   = 'mcal/cm-K-s'  # millical/cm-degK-sec, BTU/hr-ft-degF, BTU/s-in-degF, cal/s-cm-degC, W/cm-degC
+#: fac_CR, Contraction Ratio of finite area combustor (None=infinite)
+#: if make_debug_prints is True, print debugging info to terminal.
+#__init__(propName='', oxName='', fuelName='', useFastLookup=0, makeOutput=0, isp_units='sec', cstar_units='ft/sec', pressure_units='psia', temperature_units='degR', sonic_velocity_units='ft/sec', enthalpy_units='BTU/lbm', density_units='lbm/cuft', specific_heat_units='BTU/lbm degR', viscosity_units='millipoise', thermal_cond_units='mcal/cm-K-s', fac_CR=None, make_debug_prints=False)
 
 C = CEA_Obj(oxName='LOX', fuelName='RP1')
 
-full_cea = C.get_full_cea_output(Pc=100, MR=1.0, eps=None, subar=None, PcOvPe=100, show_transport=1, pc_units='bar', output='KJ', short_output=1 )
-print('this is the cea string:------------------------------------------------------------------------------------------------------------------')
-print(full_cea)
+#full_cea = C.get_full_cea_output(Pc=15.2, MR=2.3, eps=None, PcOvPe=15.2/0.9, pc_units='bar', output='KJ', short_output=1 )
+#print('this is the cea string:------------------------------------------------------------------------------------------------------------------')
+#print(full_cea)
+
+#cea_obj.print_py_cea_vars()
+#cea_obj.clearCache()                       #not working
+#print(cea_obj.getCacheDict())
+print(cea_obj.get_rocketcea_data_dir())    #not working
+
+#testing setup cards
+"""
+setupCards(Pc=100.0, MR=1.0, eps=40.0, subar=None, PcOvPe=None, frozen=0, ERphi=None, ERr=None, frozenAtThroat=0, short_output=0, show_transport=0, pc_units='psia', output='calories', show_mass_frac=False)
+#: Pc = combustion end pressure (psia)
+#: eps = Nozzle Expansion Area Ratio
+#: if PcOvPe has a value, use it instead of eps to run case
+#: ERphi = Equivalence ratios in terms of fuel-to-oxidant weight ratios.
+#: ERr = Chemical equivalence ratios in terms of valences.
+#: pc_units = 'psia', 'bar', 'atm', 'mmh'(mm of mercury)
+#: frozen flag (0=equilibrium, 1=frozen)
+#: frozenAtThroat flag, 0=frozen in chamber, 1=frozen at throat
+"""
+#run cea
+st = time.time()
+C.setupCards(Pc=15.2, MR=2.3, eps=None, PcOvPe=15.2/0.9, frozen=0, pc_units='bar', output='KJ', short_output=1)
+et = time.time()
+rt = et-st
+print(f'rocketcea run time = {rt}s')
+
+if platform.system() == "Windows":
+    h2 = Fuel("H2(L)", temp=20) # Liquid H2, Temp in Kelvin
+    lox = Oxidizer("O2(L)", temp=90)
+    # Rocket at 2000psi and supersonic area ratio of 5
+    st = time.time()
+    problem1 = RocketProblem(pressure=2000, materials=[h2, lox], phi=1, sup=5)
+    results = problem1.run()
+    et = time.time()
+    rt = et-st
+    print(f'CEA_Wrap run time = {rt}s')
+else:
+    print("windows not detected. skipping CEA_Wrap")
+
+"""
+py_cea.rockt.vaci[ self.i_thrt ] =  0.0 # Vacuum Isp at throat
+py_cea.rockt.vaci[ self.i_exit ] =  0.0 # Vacuum Isp at exit
+py_cea.rockt.cstr = 0.0                 # cstar
+py_cea.prtout.ttt[ self.i_chm ] = 0.0   # chamber temperature
+py_cea.rockt.app[ self.i_thrt ] = 0.0   # Pc/Pt
+py_cea.rockt.app[ self.i_exit ] = 0.0   # Pc/Pe
+py_cea.rockt.aeat[ self.i_exit ] = 0.0  # exit area/throat area
+py_cea.rockt.vmoc[ self.i_exit ]  = 0.0 # MACH, mach number
+py_cea.miscr.eqrat = 0.0                # equivalence ratio
+
+for i in range(3):
+    py_cea.rockt.sonvel[i] = 0.0    # SON, sonic velocity, speed of sound
+    py_cea.prtout.hsum[i] = 0.0     # enthalpies
+    py_cea.prtout.wm[i] = 0.0       # MolWt, molecular weight
+    py_cea.prtout.gammas[i] = 0.0   # gamma, heat capacity ratio
+    py_cea.prtout.vlm[i] = 0.0      # densities
+    py_cea.prtout.cpr[i] = 0.0      # heat capacities Cp
+
+    py_cea.trpts.vis[i] = 0.0    # viscosity
+    py_cea.trpts.cpeql[i] = 0.0  # equilibrium specific heat
+    py_cea.trpts.coneql[i] = 0.0 # equilibrium thermal conductivity
+    py_cea.trpts.preql[i] = 0.0  # equilibrium prandtl number
+    py_cea.trpts.cpfro[i] = 0.0  # frozen specific heat
+    py_cea.trpts.confro[i] = 0.0 # frozen thermal conductivity
+    py_cea.trpts.prfro[i] = 0.0  # frozen prandtl number
+"""
+
+# ['P,', 'p'],
+# ['T,', 't'],
+# ['Cp,', 'cp'],
+# ['GAMMAs', 'gam'],
+# ['SON', 'son'],
+# ['MACH', 'mach'],
+# ['Ae/At', 'aeat'],
+# ['CSTAR,', 'cstar'],
+# ['Ivac,', 'ivac'],
+# ['Isp,', 'isp'],
+# ['RHO,', 'rho'],
+# ['VISC,MILLIPOISE', 'mu'],
+# ['PRANDTL', 'Pr']
+
+"""
+# ------------ Finite Area Combustor String ------
+
+#print('In setupCards, fac_CR =', self.fac_CR)
+if self.fac_CR is not None:
+    N += 1
+    fac_str = ' fac ' + make_inp_str(" ac/at=", self.fac_CR )
+
+    set_py_cea_line(N, fac_str )
+
+    self.i_injface = 0
+    self.i_chm = 1
+    self.i_thrt = 2
+    self.i_exit = 3
+else:
+    self.i_injface = 0
+    self.i_chm = 0
+    self.i_thrt = 1
+    self.i_exit = 2
+"""
+
+#print variables-----------------------------------------------------------------------------------------------------
+#self.i_injface = injector, self.i_chm = chamber, self.i_thrt = thrust, self.i_exit = exit
+# using these instead of the indicies manually allows for finite area combustor support
+temperatures = py_cea.prtout.ttt[C.i_chm] 
+cp = py_cea.prtout.cpr
+gammas = py_cea.prtout.gammas
+son = py_cea.rockt.sonvel
+mach = py_cea.rockt.vmoc
+aeat = py_cea.rockt.aeat
+Cstar = py_cea.rockt.cstr
+ivac = py_cea.rockt.vaci
+#isp = ... def estimate_Ambient_Isp(self, Pc=100.0, MR=1.0, eps=40.0, Pamb=14.7, frozen=0, frozenAtThroat=0):
+rho = py_cea.prtout.vlm
+mu = py_cea.trpts.vis
+h = py_cea.prtout.hsum
+m = py_cea.prtout.wm
+"""
+        if self.fac_CR is not None:
+            sonicList = list(py_cea.rockt.sonvel[1:4])
+        else:
+            sonicList = list(py_cea.rockt.sonvel[:3])
+"""
+
+print(f'temperatures = {temperatures}')
+print(f'cp = {cp}')
+print(f'gammas = {gammas}')
+print(f'son = {son}')
+print(f'mach = {mach}')
+print(f'aeat = {aeat}')
+print(f'Cstar = {Cstar}')
+print(f'ivac = {ivac}')
+#print(f'isp = {isp}')
+print(f'rho = {rho}')
+print(f'mu = {mu}')
+print(f'h = {h}')
+print(f'm = {m}')
+
+
+
+