Application Center - Maplesoft
  Submit your work
Application Center Applications Theoretical Performance of a LOX-LH2 Rocket Preview

App Preview:

Theoretical Performance of a LOX-LH2 Rocket

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



Theoretical Performance of a LOX-LH2 Rocket

Introduction
 

Liquid hydrogen (at 20.27 K) and liquid oxygen (at 90.17 K) are burned in the combustion chamber (with an infinite area ratio) of a rocket.

 

 

 

The combustion products contain H2, O2, H, O, OH, HO2 and H2O2

 

Chemical

In Feed

In Products

H2

1

n1

O2

0.34974

n2

H2O

 

n3

O

 

n4

H

 

n5

OH

 

n6

HO2

 

n7

H2O2

 

n8
 

 

This application will calculate

 

• 

the adiabatic flame temperature and composition of of combustion products
• 

the pressures and temperatures in the throat and exit
• 

and the theoretical rocket performance, including the ideal specific impulse, characteristic velocity, sonic velocity and more

 

Assumptions

• 

The combustion chamber is large compared to the throat, hence the assumption of an infinite area ratio
• 

The flow composition does not change through the nozzle expansion (i.e. reaction rate is slow compared to flowrate). This is also known as "frozen" flow

 

Physical Properties
 
> 

restart:
> 

with(ThermophysicalData:-Chemicals):

 

Enthalpies

> 

h_O2   := Property("Hmolar", "O2",   "temperature" = T):
h_H2   := Property("Hmolar", "H2",   "temperature" = T):
h_H2O  := Property("Hmolar", "H2O",  "temperature" = T):
h_OH   := Property("Hmolar", "OH",   "temperature" = T):
h_H    := Property("Hmolar", "H",    "temperature" = T):
h_O    := Property("Hmolar", "O",    "temperature" = T):
h_HO2  := Property("Hmolar", "HO2",  "temperature" = T):
h_H2O2 := Property("Hmolar", "H2O2", "temperature" = T):

Entropies

> 

s_O2   := Property("Smolar", "O2",   "temperature" = T)  - R * ln(P_c/P_s):
s_H2   := Property("Smolar", "H2",   "temperature" = T)  - R * ln(P_c/P_s):
s_H2O  := Property("Smolar", "H2O",  "temperature" = T)  - R * ln(P_c/P_s):
s_OH   := Property("Smolar", "OH",   "temperature" = T)  - R * ln(P_c/P_s):
s_H    := Property("Smolar", "H",    "temperature" = T)  - R * ln(P_c/P_s):
s_O    := Property("Smolar", "O",    "temperature" = T)  - R * ln(P_c/P_s):
s_HO2  := Property("Smolar", "HO2",  "temperature" = T)  - R * ln(P_c/P_s):
s_H2O2 := Property("Smolar", "H2O2", "temperature" = T)  - R * ln(P_c/P_s):

 

Enthalpy of formation

> 

h_f_O2L  := Property("HeatOfFormation", "O2(L)");
h_f_H2L  := Property("HeatOfFormation", "H2(L)");
h_f_O2   := Property("HeatOfFormation", "O2");
h_f_H2   := Property("HeatOfFormation", "H2");
h_f_H2O  := Property("HeatOfFormation", "H2O");
h_f_OH   := Property("HeatOfFormation", "OH");
h_f_H    := Property("HeatOfFormation", "H");
h_f_O    := Property("HeatOfFormation", "O");
h_f_HO2  := Property("HeatOfFormation", "HO2");
h_f_H2O2 := Property("HeatOfFormation", "H2O2");

-12979.000

  

-9012.000

  

0.

  

0.

  

-241826.000

  

37278.206

  

217998.828

  

249175.003

  

12020.000

  

-135880.000

 (2.1)

Reference enthalpies

> 

h_r_O2   := eval(h_O2,   T = 298.15);
h_r_H2   := eval(h_H2,   T = 298.15);
h_r_H2O  := eval(h_H2O,  T = 298.15);
h_r_OH   := eval(h_OH,   T = 298.15);
h_r_H    := eval(h_H,    T = 298.15);
h_r_O    := eval(h_O,    T = 298.15);
h_r_HO2  := eval(h_HO2,  T = 298.15);
h_r_H2O2 := eval(h_H2O2, T = 298.15)

0.

  

-0.4957942313e-5

  

-241826.0005

  

37278.20600

  

217998.8279

  

249175.0027

  

12019.99997

  

-135880.0000

 (2.2)

Gibbs Free Energy of Formation of the combustion products

> 

Gibbs_H2 := proc(temp)
   return 0
end proc:
> 

Gibbs_O2 := proc(temp)
   return 0
end proc:
> 

Gibbs_H2O := proc(temp)
   local DeltaH, DeltaS, DeltaG:
   DeltaH := eval(h_H2O - (h_H2 + 0.5*h_O2), T = temp):
   DeltaS := eval(s_H2O - (s_H2 + 0.5*s_O2), T = temp):
   DeltaG := DeltaH - DeltaS*temp:
end proc:
> 

Gibbs_O := proc(temp)
   local DeltaH, DeltaS, DeltaG:
   DeltaH := eval(h_O - 0.5*h_O2, T = temp):
   DeltaS := eval(s_O - 0.5*s_O2, T = temp):
   DeltaG := DeltaH - DeltaS*temp:
   return DeltaG:
end proc:
> 

Gibbs_H := proc(temp)
   local DeltaH, DeltaS, DeltaG:
   DeltaH := eval(h_H - 0.5*h_H2, T = temp):
   DeltaS := eval(s_H - 0.5*s_H2, T = temp):
   DeltaG := DeltaH - DeltaS*temp:
   return DeltaG:
end proc:
> 

Gibbs_OH := proc(temp)
   local DeltaH, DeltaS, DeltaG:
   DeltaH := eval(h_OH - (0.5*h_H2 + 0.5*h_O2), T = temp):
   DeltaS := eval(s_OH - (0.5*s_H2 + 0.5*s_O2), T = temp):
   DeltaG := DeltaH - DeltaS*temp:
   return DeltaG:
end proc:
> 

Gibbs_HO2 := proc(temp)
   local DeltaH, DeltaS, DeltaG:
   DeltaH := eval(h_HO2 - (0.5*h_H2 + h_O2), T = temp):
   DeltaS := eval(s_HO2 - (0.5*s_H2 + s_O2), T = temp):
   DeltaG := DeltaH - DeltaS*temp:
   return DeltaG:
end proc:

Gibbs_H2O2 := proc(temp)
   local DeltaH, DeltaS, DeltaG:
   DeltaH := eval(h_H2O2 - (h_H2 + h_O2), T = temp):
   DeltaS := eval(s_H2O2 - (s_H2 + s_O2), T = temp):
   DeltaG := DeltaH - DeltaS*temp:
   return DeltaG:
end proc:

 

Ideal gas constant

> 

R := 8.3144:


Chamber Pressure

> 

P_c := 53.3172 * Unit(bar):


Atmospheric Pressure

> 

P_a := 1.01325 * Unit(bar):

 

Standard pressure

> 

P_s := 1.0 * Unit(bar):

 

Ratio of nozzle exit and combustion chamber throat area

> 

epsilon := 1.0:

 

Ratio of exit area to throat area

> 

AeAt := 1.58:

 

Moles of H2 and O2 fed into combustion chamber (only the ratio is important)

> 

N_H2L := 1:
N_O2L := 0.34974:

 

Mach number at throat ( = 1 for choked flow)

> 

M_t := 1:

 

Equilibrium Composition
 

Equating the moles of chemicals in the reactants and products gives

N_H2L H2 + N_O2L O2  →  n1 H2 + n2 O2 + n3 H2O + n4 H + n5 O + n6 OH + n7 HO2+ n8 H2O2

 

Balance on H atoms

2 n1 + 2 n3 + n4 + n6 + n7 + 2 n8 = 2 * N_H2L

 

Balance on O atoms

2 n2 + n3 + n5 + n6 + 2 n7 + 2 n8 = 2 * N_O2L

 

Hence the constraints are

> 

con1 := 2 * n1 + 2 * n3 + n4 + n6 + n7 + 2 * n8 = 2 * N_H2L:
> 

con2 := 2 * n2 + n3 + n5 + n6 + 2 * n7 + 2 * n8 = 2 * N_O2L:


Total moles of products

> 

nt := n1 + n2 + n3 + n4 + n5 + n6 + n7 + n8:

 

For a given temperature, minimizing the Gibbs Free Energy of the combustion products will give the equilibrium composition

> 

gibbs :=   n1 * (Gibbs_H2(T)   + R * T * ln(n1/nt))
         + n2 * (Gibbs_O2(T)   + R * T * ln(n2/nt))
         + n3 * (Gibbs_H2O(T)  + R * T * ln(n3/nt))
         + n4 * (Gibbs_H(T)    + R * T * ln(n4/nt))
         + n5 * (Gibbs_O(T)    + R * T * ln(n5/nt))
         + n6 * (Gibbs_OH(T)   + R * T * ln(n6/nt))
         + n7 * (Gibbs_HO2(T)  + R * T * ln(n7/nt))
         + n8 * (Gibbs_H2O2(T) + R * T * ln(n8/nt)):
      

Hence the values of n1, n2, n3, n4, n5, n6, n7 and n8 are given by the numeric solution of these equations, where L1 and L2 are the Lagrange multipliers.

> 

eqComposition := L1 * diff(lhs(con1), n1) + L2 * diff(lhs(con2),n1) = diff(gibbs, n1),
                 L1 * diff(lhs(con1), n2) + L2 * diff(lhs(con2),n2) = diff(gibbs, n2),
                 L1 * diff(lhs(con1), n3) + L2 * diff(lhs(con2),n3) = diff(gibbs, n3),
                 L1 * diff(lhs(con1), n4) + L2 * diff(lhs(con2),n4) = diff(gibbs, n4),
                 L1 * diff(lhs(con1), n5) + L2 * diff(lhs(con2),n5) = diff(gibbs, n5),
                 L1 * diff(lhs(con1), n6) + L2 * diff(lhs(con2),n6) = diff(gibbs, n6),
                 L1 * diff(lhs(con1), n7) + L2 * diff(lhs(con2),n7) = diff(gibbs, n7),
                 L1 * diff(lhs(con1), n8) + L2 * diff(lhs(con2),n8) = diff(gibbs, n8):

 

Heat Balance
 

The flame temperature is given by equating the heat of the reactants to the heat of the products

> 

H_reactants := N_H2L * h_f_H2L + N_O2L * h_f_O2L;

-13551.27546

 (4.1)
> 

H_products := + n1 * (h_f_H2   + (h_H2   - h_r_H2))
              + n2 * (h_f_O2   + (h_O2   - h_r_O2))
              + n3 * (h_f_H2O  + (h_H2O  - h_r_H2O))
              + n4 * (h_f_H    + (h_H    - h_r_H))
              + n5 * (h_f_O    + (h_O    - h_r_O))
              + n6 * (h_f_OH   + (h_OH   - h_r_OH))
              + n7 * (h_f_HO2  + (h_HO2  - h_r_HO2))
              + n8 * (h_f_H2O2 + (h_H2O2 - h_r_H2O2)):
> 

flameTemp := H_reactants = H_products:

 

Numerical Solution of Equilibrium Composition and Flame Temperature
 
> 

res:=fsolve({eqComposition, flameTemp, con1, con2}, {L1 = -1000, L2=-1000, T = 3000, n1 = 0.1, n2 = 0.1, n2 = 0.1, n3 = 0.1, n4 = 0.1, n5 = 0.1, n6 = 0.1, n7 = 0.1, n8 = 0.1})

{L1 = -17182.81432, L2 = -89531.03079, T = 3383.844425, n1 = .3061915322, n2 = 0.1788202357e-2, n3 = .6590840027, n4 = 0.3479245713e-1, n5 = 0.2147607467e-2, n6 = 0.3462930001e-1, n7 = 0.1551202423e-4, n8 = 0.5830551718e-5}

 (5.1)

 

Hence the temperature in the rocket combustion chamber is

> 

T_c := eval(T,res) * Unit(K)

3383.844425*Units:-Unit(K)

 (5.2)

The equilibrium composition of the combustion products are (in moles)

> 

mol_H2   := eval(n1 ,res);
mol_O2   := eval(n2 ,res);
mol_H2O  := eval(n3 ,res);
mol_H    := eval(n4 ,res);
mol_O    := eval(n5 ,res);
mol_OH   := eval(n6 ,res);
mol_HO2  := eval(n7 ,res);
mol_H2O2 := eval(n8 ,res);

.3061915322

  

0.1788202357e-2

  

.6590840027

  

0.3479245713e-1

  

0.2147607467e-2

  

0.3462930001e-1

  

0.1551202423e-4

  

0.5830551718e-5

 (5.3)
> 

mol_total := mol_H2 + mol_H2O + mol_O2 + mol_H + mol_O + mol_OH + mol_HO2 + mol_H2O2

1.038654444

 (5.4)

 Mole fractions in the combustion products

> 

molFrac_H2   := mol_H2/mol_total;
molFrac_O2   := mol_O2/mol_total;
molFrac_H2O  := mol_H2O/mol_total;
molFrac_H    := mol_H/mol_total;
molFrac_O    := mol_O/mol_total;
molFrac_OH   := mol_OH/mol_total;
molFrac_HO2  := mol_HO2/mol_total;
molFrac_H2O2 := mol_H2O2/mol_total;

.2947963435

  

0.1721652824e-2

  

.6345556085

  

0.3349762506e-1

  

0.2067682355e-2

  

0.3334053997e-1

  

0.1493473053e-4

  

0.5613562578e-5

 (5.5)

 

Ideal Performance of an Infinite Area Ratio Rocket
 
> 

with(Units[Simple]):

 

Ideal gas constant

> 

R := 8.3144 * Unit(J/mol/K):

 

Gravity

> 

grav := 9.81*Unit(m/s^2):

 

Molecular weight of the combustion products

> 

Mw_mix :=   molFrac_H2   * Property("MolarMass", "H2",   useunits)
          + molFrac_H2O  * Property("MolarMass", "H2O",  useunits)
          + molFrac_O2   * Property("MolarMass", "O2",   useunits)
          + molFrac_H    * Property("MolarMass", "H",    useunits)
          + molFrac_O    * Property("MolarMass", "O",    useunits)
          + molFrac_OH   * Property("MolarMass", "OH",   useunits)
          + molFrac_HO2  * Property("MolarMass", "HO2",  useunits)
          + molFrac_H2O2 * Property("MolarMass", "H2O2", useunits)

12.71562490*Units:-Unit(g/mol)

 (6.1)

Specific heat capacity (at constant pressure) in the combustion chamber

> 

Cp_c_mol :=   molFrac_H2   * Property("Cpmolar", "H2",  "temperature" = T_c)
            + molFrac_H2O  * Property("Cpmolar", "H2O",  "temperature" = T_c)
            + molFrac_O2   * Property("Cpmolar", "O2",   "temperature" = T_c)
            + molFrac_O    * Property("Cpmolar", "O",    "temperature" = T_c)
            + molFrac_H    * Property("Cpmolar", "H",    "temperature" = T_c)
            + molFrac_OH   * Property("Cpmolar", "OH",   "temperature" = T_c)
            + molFrac_HO2  * Property("Cpmolar", "HO2",  "temperature" = T_c)
            + molFrac_H2O2 * Property("Cpmolar", "H2O2", "temperature" = T_c)

50.01896921*Units:-Unit(J/(mol*K))

 (6.2)

Specific heat capacity (at constant volume) in the combustion chamber

> 

Cv_c_mol := Cp_c_mol - R

41.70456921*Units:-Unit(J/(mol*K))

 (6.3)

Isentropic expansion coefficient in the chamber

> 

Gamma_c := Cp_c_mol / Cv_c_mol

1.199364246

 (6.4)

Mach number at exit

> 

M_e := fsolve(AeAt = ((Gamma_c + 1)/2) ^ (-(Gamma_c + 1)/(2 * (Gamma_c - 1))) * (1 + 0.5 * (Gamma_c - 1) * M_e^2)^((Gamma_c + 1)/(2 * (Gamma_c - 1)))/M_e, M_e = 1)

.4108552984

 (6.5)

 

Throat temperature

> 

T_t := T_c * (1 + (Gamma_c - 1)/2 * M_t^2)^(-1)

3077.111425*Units:-Unit(K)

 (6.6)

Exit temperature

> 

T_e := T_c * (1 + (Gamma_c - 1)/2 * M_e^2)^(-1)

3327.848224*Units:-Unit(K)

 (6.7)

 

Specific heat capacity (at constant pressure) at throat

> 

Cp_t_mol :=   molFrac_H2   * Property("Cpmolar", "H2",   "temperature" = T_t)
            + molFrac_H2O  * Property("Cpmolar", "H2O",  "temperature" = T_t)
            + molFrac_O2   * Property("Cpmolar", "O2",   "temperature" = T_t)
            + molFrac_H    * Property("Cpmolar", "H",    "temperature" = T_t)
            + molFrac_O    * Property("Cpmolar", "O",    "temperature" = T_t)
            + molFrac_OH   * Property("Cpmolar", "OH",   "temperature" = T_t)
            + molFrac_HO2  * Property("Cpmolar", "HO2",  "temperature" = T_t)
            + molFrac_H2O2 * Property("Cpmolar", "H2O2", "temperature" = T_t)

49.24896198*Units:-Unit(J/(mol*K))

 (6.8)

Isentropic expansion coefficient at throat

> 

Gamma_t := Cp_t_mol / (Cp_t_mol - R )

1.203114424

 (6.9)

Throat  pressure

> 

P_t := P_c * (1 + (Gamma_c - 1)/2 * M_t^2)^(Gamma_c / (1 - Gamma_c))

30.10276292*Units:-Unit(bar)

 (6.10)

Exit pressure

> 

P_e := P_c * (1 + (Gamma_t - 1)/2 * M_e^2)^(Gamma_t / (1 - Gamma_t))

48.21043268*Units:-Unit(bar)

 (6.11)

Chamber gas density

> 

rho_c := P_c * Mw_mix/(R * T_c)

2.409704230*Units:-Unit(kg/m^3)

 (6.12)

Throat gas density

> 

rho_t :=  P_t * Mw_mix/(R * T_t )

1.496132060*Units:-Unit(kg/m^3)

 (6.13)

Sonic velocity in chamber and throat

> 

sonicVelocity_c := sqrt((Gamma_c * P_c/rho_c))

1629.023487*Units:-Unit(m/s)

 (6.14)
> 

sonicVelocity_t := sqrt((Gamma_t * P_t/rho_t))

1555.864176*Units:-Unit(m/s)

 (6.15)

Throat velocity for an isentropic nozzle

> 

V_t := sqrt(2*T_c*R*Gamma_c*(1-(P_t/P_c)^((Gamma_c-1)/Gamma_c))/(Mw_mix*(Gamma_c-1)))

1553.437424*Units:-Unit(m/s)

 (6.16)

Ideal specific impulse

> 

Isp_ideal := V_t / grav

158.3524388*Units:-Unit(s)

 (6.17)

Ideal specific impulse as defied by NASA CEA

> 

Isp_ideal_NASA := V_t

1553.437424*Units:-Unit(m/s)

 (6.18)

Ideal specific impulse in a vacuum.

> 

Isp_vac := V_t + P_t / (rho_t * V_t)

2848.654810*Units:-Unit(m/s)

 (6.19)

Characteristic velocity (C-Star)

> 

Cstar := sqrt(1 / Gamma_c * ((Gamma_c + 1)/2) ^ ((Gamma_c + 1)/(Gamma_c - 1)) * R / Mw_mix * T_c)

2294.054024*Units:-Unit(m/s)

 (6.20)
> 

 

Legal Notice: © Maplesoft, a division of Waterloo Maple Inc. 2018. Maplesoft and Maple are trademarks of Waterloo Maple Inc. This application may contain errors and Maplesoft is not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact Maplesoft for permission if you wish to use this application in for-profit activities.