Closed RasmusAPedersen closed 2 years ago
In the first run, you are using units of Pa for pressure and degK for output temperature. In the second, you are using units of psia for pressure and degR for output temperature.
The following code will correct the values.
from rocketcea.cea_obj_w_units import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8',
pressure_units='Pa',
cstar_units='m/s',
temperature_units='K',
sonic_velocity_units='m/s',
enthalpy_units='kJ/kg', # Need to be aware that this is kJ
density_units='kg/m^3',
specific_heat_units='J/kg-K',
thermal_cond_units='W/cm-degC') # Need to be aware of cm here
P_comb = 5.31e5 # 5.31 bar
P_comb_psia = 0.000145038*P_comb # == 77.015
print( 'P_comb=%g, P_comb_psia=%g'%(P_comb, P_comb_psia) )
er = 1
print(combustor.get_Temperatures(Pc=P_comb, MR=8.9,eps=er))
print(combustor.get_Temperatures(Pc=P_comb, MR=9.0,eps=er))
print(combustor.get_Temperatures(Pc=P_comb, MR=9.1,eps=er))
print('---')
from rocketcea.cea_obj import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8')
print( [t/1.8 for t in combustor.get_Temperatures(Pc=P_comb_psia, MR=8.9,eps=er)] )
print( [t/1.8 for t in combustor.get_Temperatures(Pc=P_comb_psia, MR=9.0,eps=er)] )
print( [t/1.8 for t in combustor.get_Temperatures(Pc=P_comb_psia, MR=9.1,eps=er)] )
```output is:
P_comb=531000, P_comb_psia=77.0152
[2923.661514085015, 2775.686615971293, 2775.6866159713377]
[2915.793048468526, 2767.0340827573127, 2767.0340827572713]
[2907.9173515455013, 2758.358534637467, 2758.3585346375326]
---
[2923.661662450697, 2775.6867365248213, 2775.686736524833]
[2915.793194752911, 2767.034201047843, 2767.034201047821]
[2907.9174957518553, 2758.358650669505, 2758.3586506694633]
I am aware of the units. I remade the CEA object in imperial units to show that the issue was similar for both unit sets.
You will notice that the outputs of my function calls are identical for MR=9.0 and 9.1
That is the issue
The problem is the use of an integer for area ratio, for future versions I will consider casting the input to float. When I run the following code:
from rocketcea.cea_obj_w_units import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8',
pressure_units='Pa',
cstar_units='m/s',
temperature_units='K',
sonic_velocity_units='m/s',
enthalpy_units='kJ/kg', # Need to be aware that this is kJ
density_units='kg/m^3',
specific_heat_units='J/kg-K',
thermal_cond_units='W/cm-degC') # Need to be aware of cm here
P_comb = 5.31e5 # 5.31 bar
er = 1.0
print(combustor.get_Temperatures(Pc=P_comb, MR=8.9,eps=er))
print(combustor.get_Temperatures(Pc=P_comb, MR=9.0,eps=er))
print(combustor.get_Temperatures(Pc=P_comb, MR=9.1,eps=er))
"""
I get
[2923.661514085015, 2775.686615971293, 2775.6866159713377]
[2915.793048468526, 2767.0340827573127, 2767.0340827572713]
[2907.9173515455013, 2758.358534637467, 2758.3585346375326]
"""
Correcting from
er = 1
to
er = 1.0
Does not seem to fix the issue for me, unfortunately. In fact it seems sort of random when it uses the same value or not.
Tinkered a little more with it. Running this code:
from rocketcea.cea_obj import set_rocketcea_data_dir
set_rocketcea_data_dir( r'D:\temp\RocketCEA' )
from rocketcea.cea_obj_w_units import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8',
pressure_units='Pa',
cstar_units='m/s',
temperature_units='K',
sonic_velocity_units='m/s',
enthalpy_units='kJ/kg', # Need to be aware that this is kJ
density_units='kg/m^3',
specific_heat_units='J/kg-K',
thermal_cond_units='W/cm-degC') # Need to be aware of cm here
P_comb = 5.31e5 # 5.31 bar
er = 1.0
print(combustor.get_Temperatures(Pc=float(P_comb), MR=float(8.9),eps=float(er)))
print(combustor.get_Temperatures(Pc=float(P_comb), MR=float(9.0),eps=float(er)))
print(combustor.get_Temperatures(Pc=float(P_comb), MR=float(9.1),eps=float(er)))
Gives me this:
Setting ROCKETCEA_DATA_DIR = D:\temp\RocketCEA
[2994.1379217115486, 2852.486415593293, 2852.4864153875865]
[2915.777219727355, 2767.0212830485807, 2767.02128304865]
[2915.777219727355, 2767.0212830485807, 2767.02128304865]
Which is still the same issue. If I change MR to 10, it gives me a different return, and it only seems to update its output on whole integers of MR.
Running CEA with no nozzle at all (i.e. an area ratio of 1.0) is a bit unusual. Try putting a small nozzle on. For example a 2:1 nozzle. When I run er = 2.0, I get the same chamber and throat temperatures as er = 1.0.
from rocketcea.cea_obj_w_units import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8',
pressure_units='Pa',
cstar_units='m/s',
temperature_units='K',
sonic_velocity_units='m/s',
enthalpy_units='kJ/kg', # Need to be aware that this is kJ
density_units='kg/m^3',
specific_heat_units='J/kg-K',
thermal_cond_units='W/cm-degC') # Need to be aware of cm here
P_comb = 5.31e5 # 5.31 bar
er = 2.0
print(combustor.get_Temperatures(Pc=P_comb, MR=8.9,eps=er))
print(combustor.get_Temperatures(Pc=P_comb, MR=9.0,eps=er))
print(combustor.get_Temperatures(Pc=P_comb, MR=9.1,eps=er))
"""
with 1.0 I get
[2923.661514085015, 2775.686615971293, 2775.6866159713377]
[2915.793048468526, 2767.0340827573127, 2767.0340827572713]
[2907.9173515455013, 2758.358534637467, 2758.3585346375326]
with 2.0 I get
[2923.661514085015, 2775.686615971293, 2376.9687762950966]
[2915.793048468526, 2767.0340827573127, 2364.698294153291]
[2907.9173515455013, 2758.358534637467, 2352.3556307779886]
"""
The issue remains still.
Is there some kind of printout of my system that I can give you?
OK... try the following code to generate standard CEA output. Since you are running such low pressure and extremely oxidizer rich, it is possible that the FORTRAN in CEA, when driven by RocketCEA, is leaving some internal iteration termination criteria as satisfied when the mixture ratio changes so slightly (that's a total guess on my part)
Each of the runs should give a detailed output with the correct results.
Also, since I can not reliably reproduce your error, tell me what version of python you are using and a little about your Windows environment.
from rocketcea.cea_obj import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8')
P_comb = 5.31e5 # 5.31 bar
P_comb_psia = 0.000145038*P_comb # == 77.015
er = 2.0
print(combustor.get_full_cea_output(Pc=P_comb_psia, MR=8.9,eps=er, short_output=1))
print(combustor.get_full_cea_output(Pc=P_comb_psia, MR=9.0,eps=er, short_output=1))
print(combustor.get_full_cea_output(Pc=P_comb_psia, MR=9.1,eps=er, short_output=1))
I ran the code as you suggested.
from rocketcea.cea_obj import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8')
P_comb = 5.31e5 # 5.31 bar
P_comb_psia = 0.000145038*P_comb # == 77.015
er = 2.0
print(combustor.get_full_cea_output(Pc=P_comb_psia, MR=8.9,eps=er, short_output=1))
print(combustor.get_full_cea_output(Pc=P_comb_psia, MR=9.0,eps=er, short_output=1))
print(combustor.get_full_cea_output(Pc=P_comb_psia, MR=9.1,eps=er, short_output=1))
And the output is below.
It's rather curious that CEA reports the same OF being calculated twice despite it seemingly having understood the input?
*******************************************************************************
NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA, OCTOBER 18, 2002
BY BONNIE MCBRIDE AND SANFORD GORDON
REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996
*******************************************************************************
reac
fuel C3H8(L) C 3 H 8 wt%=100.
h,cal=-30372. t(k)=231.08 rho=0.5808
oxid O2(G) O 2
h,cal=0. t(k)=298.15 wt%=100.
prob case=RocketCEA,
rocket equilibrium p,psia=77.015178, supar=2.000000,
o/f=8.900000,
output calories short transport
end
THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM
COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR
Pinj = 77.0 PSIA
CASE = RocketCEA,
REACTANT WT FRACTION ENERGY TEMP
(SEE NOTE) CAL/MOL K
FUEL C3H8(L) 1.0000000 -30372.000 231.000
OXIDANT O2(G) 1.0000000 0.000 298.000
O/F= 8.00000 %FUEL= 11.111111 R,EQ.RATIO= 0.453543 PHI,EQ.RATIO= 0.453543
CHAMBER THROAT EXIT
Pinf/P 1.0000 1.7253 7.5280
P, ATM 5.2395 3.0369 0.69600
T, K 2994.14 2852.49 2483.59
RHO, G/CC 6.1074-4 3.7603-4 1.0172-4
H, CAL/G -76.531 -186.49 -452.19
U, CAL/G -284.29 -382.07 -617.89
G, CAL/G -7235.75 -7007.01 -6390.65
S, CAL/(G)(K) 2.3911 2.3911 2.3911
M, (1/n) 28.639 28.983 29.785
(dLV/dLP)t -1.02493 -1.02032 -1.00922
(dLV/dLT)p 1.5385 1.4628 1.2449
Cp, CAL/(G)(K) 1.2077 1.1202 0.8221
GAMMAs 1.1249 1.1244 1.1319
SON VEL,M/SEC 988.9 959.2 885.9
MACH NUMBER 0.000 1.000 2.001
TRANSPORT PROPERTIES (GASES ONLY)
CONDUCTIVITY IN UNITS OF MILLICALORIES/(CM)(K)(SEC)
VISC,MILLIPOISE 1.0165 0.98285 0.89347
WITH EQUILIBRIUM REACTIONS
Cp, CAL/(G)(K) 1.2077 1.1202 0.8221
CONDUCTIVITY 2.3679 2.1005 1.3179
PRANDTL NUMBER 0.5185 0.5242 0.5573
WITH FROZEN REACTIONS
Cp, CAL/(G)(K) 0.3921 0.3900 0.3832
CONDUCTIVITY 0.5610 0.5373 0.4745
PRANDTL NUMBER 0.7106 0.7134 0.7216
PERFORMANCE PARAMETERS
Ae/At 1.0000 2.0000
CSTAR, FT/SEC 4828.9 4828.9
CF 0.6517 1.2046
Ivac,LB-SEC/LB 184.8 220.7
Isp, LB-SEC/LB 97.8 180.8
MOLE FRACTIONS
*CO 0.03840 0.03072 0.01276
*CO2 0.17809 0.18837 0.21240
*H 0.00569 0.00429 0.00146
HO2 0.00017 0.00012 0.00004
*H2 0.00711 0.00585 0.00274
H2O 0.24372 0.25426 0.27999
H2O2 0.00001 0.00000 0.00000
*O 0.03140 0.02490 0.01062
*OH 0.06978 0.05961 0.03346
*O2 0.42564 0.43187 0.44654
* THERMODYNAMIC PROPERTIES FITTED TO 20000.K
NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS
*******************************************************************************
NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA, OCTOBER 18, 2002
BY BONNIE MCBRIDE AND SANFORD GORDON
REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996
*******************************************************************************
reac
fuel C3H8(L) C 3 H 8 wt%=100.
h,cal=-30372. t(k)=231.08 rho=0.5808
oxid O2(G) O 2
h,cal=0. t(k)=298.15 wt%=100.
prob case=RocketCEA,
rocket equilibrium p,psia=77.015178, supar=2.000000,
o/f=9.000000,
output calories short transport
end
THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM
COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR
Pinj = 77.0 PSIA
CASE = RocketCEA,
REACTANT WT FRACTION ENERGY TEMP
(SEE NOTE) CAL/MOL K
fuel C3H8(L) 100.0000000 -30372.000 231.000
oxid O2(G) 100.0000000 0.000 298.000
O/F= 9.00000 %FUEL= 10.000000 R,EQ.RATIO= 0.403149 PHI,EQ.RATIO= 0.403149
CHAMBER THROAT EXIT
Pinf/P 1.0000 1.7296 7.6469
P, ATM 5.2395 3.0294 0.68518
T, K 2915.78 2767.02 2364.69
RHO, G/CC 6.3999-4 3.9416-4 1.0674-4
H, CAL/G -68.878 -174.15 -427.75
U, CAL/G -267.14 -360.27 -583.20
G, CAL/G -6875.00 -6633.04 -5947.51
S, CAL/(G)(K) 2.3342 2.3342 2.3342
M, (1/n) 29.225 29.543 30.228
(dLV/dLP)t -1.01876 -1.01443 -1.00498
(dLV/dLT)p 1.4184 1.3411 1.1408
Cp, CAL/(G)(K) 1.0235 0.9279 0.6362
GAMMAs 1.1298 1.1312 1.1488
SON VEL,M/SEC 968.1 938.6 864.4
MACH NUMBER 0.000 1.000 2.005
TRANSPORT PROPERTIES (GASES ONLY)
CONDUCTIVITY IN UNITS OF MILLICALORIES/(CM)(K)(SEC)
VISC,MILLIPOISE 0.99688 0.96114 0.86251
WITH EQUILIBRIUM REACTIONS
Cp, CAL/(G)(K) 1.0235 0.9279 0.6362
CONDUCTIVITY 1.9194 1.6479 0.9228
PRANDTL NUMBER 0.5316 0.5412 0.5946
WITH FROZEN REACTIONS
Cp, CAL/(G)(K) 0.3818 0.3794 0.3714
CONDUCTIVITY 0.5319 0.5077 0.4407
PRANDTL NUMBER 0.7155 0.7183 0.7270
PERFORMANCE PARAMETERS
Ae/At 1.0000 2.0000
CSTAR, FT/SEC 4708.3 4708.3
CF 0.6540 1.2075
Ivac,LB-SEC/LB 180.3 215.0
Isp, LB-SEC/LB 95.7 176.7
MOLE FRACTIONS
*CO 0.02613 0.01954 0.00583
*CO2 0.17270 0.18145 0.19983
*H 0.00364 0.00254 0.00058
HO2 0.00016 0.00011 0.00003
*H2 0.00480 0.00373 0.00131
H2O 0.22907 0.23875 0.26128
H2O2 0.00001 0.00000 0.00000
*O 0.02526 0.01895 0.00609
*OH 0.05866 0.04836 0.02262
*O2 0.47956 0.48657 0.50244
* THERMODYNAMIC PROPERTIES FITTED TO 20000.K
NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS
*******************************************************************************
NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA, OCTOBER 18, 2002
BY BONNIE MCBRIDE AND SANFORD GORDON
REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996
*******************************************************************************
reac
fuel C3H8(L) C 3 H 8 wt%=100.
h,cal=-30372. t(k)=231.08 rho=0.5808
oxid O2(G) O 2
h,cal=0. t(k)=298.15 wt%=100.
prob case=RocketCEA,
rocket equilibrium p,psia=77.015178, supar=2.000000,
o/f=9.100000,
output calories short transport
end
THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM
COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR
Pinj = 77.0 PSIA
CASE = RocketCEA,
REACTANT WT FRACTION ENERGY TEMP
(SEE NOTE) CAL/MOL K
fuel C3H8(L) 100.0000000 -30372.000 231.000
oxid O2(G) 100.0000000 0.000 298.000
O/F= 9.00000 %FUEL= 10.000000 R,EQ.RATIO= 0.403149 PHI,EQ.RATIO= 0.403149
CHAMBER THROAT EXIT
Pinf/P 1.0000 1.7296 7.6469
P, ATM 5.2395 3.0294 0.68518
T, K 2915.78 2767.02 2364.69
RHO, G/CC 6.3999-4 3.9416-4 1.0674-4
H, CAL/G -68.878 -174.15 -427.75
U, CAL/G -267.14 -360.27 -583.20
G, CAL/G -6875.00 -6633.04 -5947.51
S, CAL/(G)(K) 2.3342 2.3342 2.3342
M, (1/n) 29.225 29.543 30.228
(dLV/dLP)t -1.01876 -1.01443 -1.00498
(dLV/dLT)p 1.4184 1.3411 1.1408
Cp, CAL/(G)(K) 1.0235 0.9279 0.6362
GAMMAs 1.1298 1.1312 1.1488
SON VEL,M/SEC 968.1 938.6 864.4
MACH NUMBER 0.000 1.000 2.005
TRANSPORT PROPERTIES (GASES ONLY)
CONDUCTIVITY IN UNITS OF MILLICALORIES/(CM)(K)(SEC)
VISC,MILLIPOISE 0.99688 0.96114 0.86251
WITH EQUILIBRIUM REACTIONS
Cp, CAL/(G)(K) 1.0235 0.9279 0.6362
CONDUCTIVITY 1.9194 1.6479 0.9228
PRANDTL NUMBER 0.5316 0.5412 0.5946
WITH FROZEN REACTIONS
Cp, CAL/(G)(K) 0.3818 0.3794 0.3714
CONDUCTIVITY 0.5319 0.5077 0.4407
PRANDTL NUMBER 0.7155 0.7183 0.7270
PERFORMANCE PARAMETERS
Ae/At 1.0000 2.0000
CSTAR, FT/SEC 4708.3 4708.3
CF 0.6540 1.2075
Ivac,LB-SEC/LB 180.3 215.0
Isp, LB-SEC/LB 95.7 176.7
MOLE FRACTIONS
*CO 0.02613 0.01954 0.00583
*CO2 0.17270 0.18145 0.19983
*H 0.00364 0.00254 0.00058
HO2 0.00016 0.00011 0.00003
*H2 0.00480 0.00373 0.00131
H2O 0.22907 0.23875 0.26128
H2O2 0.00001 0.00000 0.00000
*O 0.02526 0.01895 0.00609
*OH 0.05866 0.04836 0.02262
*O2 0.47956 0.48657 0.50244
* THERMODYNAMIC PROPERTIES FITTED TO 20000.K
NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL OXIDANTS
It looks like the FORTRAN CEA code is not reading the input deck correctly. (i.e. the O/F in the input deck is correct, but CEA runs the integer value of the O/F)
What can you tell me about how you compiled the FORTRAN?
Perhaps to make sure that is the problem for your calls to get_Temperatures, try the following code. I get the smooth curve shown below, you may get stairsteps.
from pylab import *
from rocketcea.cea_obj import CEA_Obj
combustor = CEA_Obj(oxName='GOX', fuelName='C3H8')
er = 2.0
P_comb_psia = 77.015
mrL = [imr/10.0 for imr in range(1,100)]
TcL = [combustor.get_Temperatures(Pc=P_comb_psia, MR=mr,eps=er)[0] for mr in mrL]
plot(mrL, TcL)
title("Chamber Temperature")
show()
I ran the code you pasted, and it outputs this. It's quite clear it only runs a different OF value on whole integers. Something somewhere seems to get casted to an int?
As for how I compiled Fortran, I don't think I'm qualified enough to give a detailed report on that. I'm on a Win10 machine, and I followed your provided instructions here: https://rocketcea.readthedocs.io/en/latest/installgfortran.html
Previously, I had some issues getting rocketcea to work, as discussed in this issue: #18
I'm pretty sure that the FORTRAN compile is the problem. (i.e. RocketCEA creates a correct input deck to the FORTRAN code, but the FORTRAN misinterprets a float as an integer.
The FORTRAN compile is administered by the numpy.f2py code, so it is certainly possible that python version or numpy version could affect it. Anaconda may also have some special handling for FORTRAN.
Because Windows is so tricky to get the compile correct, and I am so unfamiliar with your particular setup, I don't think I can solve the FORTRAN compile from afar.
I have 3 recommendations, none of which are perfect, but they would allow you to run cases on your Windows machine.
1) The most "native" option, would be to install a version of python outside of Anaconda and follow the directions for Windows at: https://rocketcea.readthedocs.io/en/latest/quickstart.html#windows-batch-file
I'm assuming you already have MinGW installed from:
https://rocketcea.readthedocs.io/en/latest/installgfortran.html#windows-10
The most trouble-free version of python is probably 3.7 from:
https://www.python.org/downloads/release/python-379/
I recommend you start with 3.7 and only after that works, try moving to more recent versions of python.
You would need to either set the path in your terminal to find the correct python EXE file, or call the full path each time (e.g. D:\python37_64\python.exe script_name.py). I like batch files for each version (e.g. py37 script_name.py)
2) Run cases on Google Colaboratory as shown in: https://rocketcea.readthedocs.io/en/latest/quickstart.html#google-colaboratory
3) Run cases on Windows 10 with WSL as shown in: https://rocketcea.readthedocs.io/en/latest/quickstart.html#windows-10-with-wsl
Please post a comment of your results (especially if you are able to solve the original FORTRAN problem in Anaconda)
I'm sorry I wasn't more help.
When I run this code:
I get
The two returns of MR=9.0 and MR=9.1 are identical. I assume this might have something to do with caching, but I haven't been able to find the error in the source code.