The Portland State Aerospace Society (PSAS) is an engineering student group and citizen science project at Portland State University dedicated to building low-cost, open-source, open-hardware rockets and avionics systems. The group’s stated long term goal is to place a 1 kg cubesat into low Earth orbit with their own launch vehicle. One step needed to achieve this goal is to transition the current rocket design from a solid motor to a liquid fuel engine. The liquid propelled rocket engine project is being conducted as part of a mechanical engineering senior capstone project at Portland State University. This project is on track to being completed by June of 2016.
The complexity and cost of building a liquid fuel rocket engine typically makes such devices unobtainable for a majority of parties interested in their construction. Until recently, manufacturing processes and techniques limited the geometries available to the designer and rendered such engines cost prohibitive as options for inexpensive orbital space flight. Advances in additive manufacturing technologies provide the potential to prototype complex geometries on a lower budget and with shorter lead times which would be considered unobtainable with traditional manufacturing methods. Furthermore, bipropellant liquid fuels offer many complex engineering considerations; the full analysis of which may not be within the design ability of many amateur builds. It is therefore advantageous to develop techniques for additive manufacturing rapid prototyping to make the study of bipropellant fuels more accessible.
Explored herein is the process of designing and testing a liquid bipropellant engine on the scale of 50-1000 lbf of thrust using liquid oxygen (LOX) and ethanol as propellants. A low cost pintle injector and accompanying regeneratively cooled thrust chamber is developed using a combination of traditional manufacturing techniques and additive processes. Equations have been determined to describe the more complex geometries of the nozzle contour and the sizing of other important components. Figure 1 shows the preliminary design of the combustion chamber and nozzle with cooling channels. Heat transfer analysis is being conducted to determine the type of metal to be selected for the nozzle and cooling chamber, most likely a high-temperature steel such as inconel.
Figure 1: Combustion chamber and nozzle preliminary design
In order to achieve a simple and easy to manufacture pintle injector design, the regenerative cooling channel interface, fuel manifold and most of the injector plate will be part of the additively manufactured combustion chamber. The design will accommodate easy installation of different injector types to allow for testing of different injector and spray configurations. Cold flow tests will be conducted to study the viscous losses in 3D-printed metal cooling channels, manifold system, and pintle injector designs.
An engine test-stand has been constructed and a pressure fed fuel system with actuated valves is being developed for fuel delivery. Pressure transducers and thermocouples will be ported into the final engine design to collect pressure and temperature data in the nozzle and combustion chamber. Analytical solutions and empirical data will be compared against simulations carried out by computational fluid dynamics (CFD) simulations. CFD simulations will be conducted using Loci/CHEM and Star-CCM+ simulation software.
Parameters must be selected. The generation of force is the desire of a rocket engine. the force equation is:
$$ F = \frac{\dot{w}}{g} V_{e} + A_{e}(P_{e}-P_{a}) $$This is important when designing engines for flight, but when testing engines $P_{e}$ is optimally expanded and assumed to be atmospheric pressure (14.7 psia) and the equation becomes:
$$ F = \frac{\dot{w}}{g} V_{e} $$For engines in flight the effective exhaust velocity, $c$, can be used, and equation X can be expressed as:
$$F = c\frac{\dot{W}}{g} $$Where c is defined as
$$c=V_{e} + A_{e}(P_{e}-P_{a})\frac{g}{\dot{W}} $$which changes with altitude.
For a desired force, there must then be selected an exit velocity to find the necessary flow rate, $\dot{w}$. The theoretical exit velocity is defined as:
$$V_{e} = \sqrt{\frac{2g\gamma}{\gamma-1}RT_{i}\bigg[1-\bigg(\frac{P_{e}}{P_{i}}\bigg)\bigg]^\frac{\gamma-1}{\gamma}+v_{i}^2} $$Because the inlet velocity is very small, it is assumed to be zero, this gives the following:
$$V_{e} = \sqrt{\frac{2g\gamma}{\gamma-1}RT_{c_{ns}}\bigg[1-\bigg(\frac{P_{e}}{P_{c_{ns}}}\bigg)\bigg]^\frac{\gamma-1}{\gamma}} $$which is dependent on the propellants, chamber pressure which should be chosen for the design, $P_{i}$, and the mach number at the nozzle inlet, $M_{i}$, which must be calculated iteratively.
Pressure at the exit, $P_{e}$, of a test nozzle should be expanded to the pressure of testing conditions (assumed to be sea-level atmospheric pressure).
Flame temperature and gamma can be obtained for a given propellant combination by using NASA's CEArun tool. link to portion in document which describes how to use cearun
The propellent mixture ratio, can be used to manipulate flame temperature as well as cooling capacity. Increasing fuel; and thereby decreasing rw, will result in lower temperatures and more coolant mass flow-rate. CEArun should be used to determine the effect of a given mixture ratio.
Things to add to this section: choosing OF ratios, (maximize theoretical specific impulse, balance with temperature limits, increases in oxidizer tend to lead to higher temperatures. selecting chamber pressures should be dependent on material limits, impulse tends to increase with increases in pressure, more enthalpy is available to convert to exhaust velocity, higher pressures generally result in higher temperatures.
define constants to be used in mathmatics. discuss gas constants and their calculations.
discuss each parameter as well as nomenclature, what does inj, inlet, throat, or exit mean, show locations in diagrams.
import math
from mpmath import *
from IPython.display import display
from ipywidgets import widgets
#Define Constants
g = 32.2 # gravitational constant in ft/s^2
J = 778 # Energy conversion factor (ft-lb/Btu)
Regas = 1544 # Gas constant (ft/degR)
Rgas = 8.314 # Gas constant (J/mol/K)
Force: The desired force output of the nozzle. Increasing force will increase mass flow rate of propellants, and can be used to increase cooling capacity.
Molecular weight: CEARun Parameter
Atmospheric Pressure: Pressure of the open atmosphere at desired optimized altitude. 14.7 psi for sea level applications.
Exit Pressure: Desired expansion pressure. Equal to atmospheric pressure in an ideally expanded nozzle.
Pressure at Injector: Desired chamber pressure after drop across injector. Higher chamber pressures allow for longer super sonic flow sections and greater expansion ratios, but will result in larger cooling surface areas and heavier nozzles.
Weight mixture ratio: Optimized using CEARun for particular propellant combinations. Use of fuel rich ratios can aid in cooling by increasing fuel mass flow rate and thereby cooling capacity. Non ideal mixtures can have the effect of small changes in flame temperature.
Fuel mixture percentage: mass percent of aqueous fuel solutions. Lowering can increase cooling capacity by increasing the mass flow rate of coolant, as well as raising the specific heat of the coolant. Results in less available fuel for combustion and lower combustion efficiencies.
Ambient temperature: Local air temperature, assumed to be the holding temperature of fuel. Refridgeration of fuel can be acheived by altering ambient temperature. (this could potentially be altered by adding a refridgeration temperature differential value, or solving for needed refridgeration temperature in order to provide sufficient cooling capacity for an engine near convergence on heat transfer.
Nozzle Stagnation Temperature: Flame temperature at stagnation, assumed to be chamber flame temperature. Obtained by CEARun.
Contraction Ratio: Ratio of chamber cross sectional area to throat cross sectional area. Reasonable values may be obtained by Fig 4-9 Huang (add picture of fig 4-9 xx). Small rockets may increase avaible tooling area by increasing this value. Larger values result in larger cooling surface areas and greater cooling capacity requirements. Decreasing can aid in heat transfer convergence.
Ratio of specific heats: The ratio of specific heats of propellants. Obtained via CEARun.
#initial design considerations inputs
#consider calculated exit pressure based on desired optimized altitude and then printing results by iterating
#through various ambient pressures as a tool for designing flight engines.
#-------------------------------------------------------------------------
#also consider grouping all CEArun generated perameters together and dedicating a whole section to using the tool.
F_ = widgets.Text("500" , description="Force, F (lbf)" , width=60 )
M_ = widgets.Text("24.029" , description="Molecular Weight, (lb/lb-mol)" , width=100 )
Pa_ = widgets.Text("14.7" , description="Ambient Pressure, (lbf/in^2)" , width=60 )
Pe_ = widgets.Text("14.7" , description="Exit Pressure, Pe (lbf/in^2) " , width=60 )
Pinj_ = widgets.Text("350" , description="Pressure at injector, Pinj (lb/in^2)" , width=60 )
rw_ = widgets.Text("1.2" , description="(oxidizer/fuel) Weight mixture ratio" , width=60 )
rwfuel_ = widgets.Text("59.5" , description="Fuel mixture percentage" , width=60 )
Tamb_ = widgets.Text("551" , description="Ambient temperature, (deg R)" , width=60 )
Tcns_ = widgets.Text("5487.1" , description="Nozzle stagnation temperature, (deg R)" , width=80 )
epsilonc_ = widgets.Text("10.5" , description="Contraction ratio" , width=60 )
gam_ = widgets.Text("1.1208" , description="Specific heat ratio" , width=80 )
display(F_ )
display(M_ )
display(Pa_ )
display(Pe_ )
display(Pinj_ )
display(rw_ )
display(rwfuel_ )
display(Tamb_ )
display(Tcns_ )
display(epsilonc_ )
display(gam_ )
F = float(F_.value ) #Nozzle force
M = float(M_.value ) #Molecular Weight
Pa = float(Pa_.value ) #Atmospheric Pressure
Pe = float(Pe_.value ) #Pressure at exit
Pinj = float(Pinj_.value ) #Pressure at injector
rw = float(rw_.value ) #Weight mixture ratio
rwfuel = float(rwfuel_.value ) #Fuel mixture percentage
Tamb = float(Tamb_.value ) #Ambient temperature
Tcns = float(Tcns_.value ) #Nozzle stagnation temperature
epsilonc = float(epsilonc_.value ) #Contraction ratio
gam = float(gam_.value ) #Specific Heat Ratio
The mach number at the inlet varies with dependence on the contraction ratio. Higher contraction ratios lead to smaller mach numbers at the nozzle inlet, small engines tend to have large contraction ratios in order to be machinable, which results in small mach numbers at the inlet.
(Figure 1-10 (huang) could be displayed here xx)
#Iterative process to find mach number at the inlet.
#equation 1-23, huang
Mi_current = .3 # a reasonable guess for mach number at inlet
Mi_last = 0
while (abs(Mi_current-Mi_last)>0.0001):
Mi_last = Mi_current
Mi_current = math.sqrt(((1+(gam-1)/2*Mi_last)/((gam+1)/2))**((gam+1)/(gam-1)))/epsilonc
Mi = Mi_current
print("Mach Number at inlet: %.2f" % Mi)
Mach Number at inlet: 0.06
#Assumptions
Pcinj = Pinj #Chamber total pressure is equal to injector pressure
Minj = 0 #Mach number at injector assumed to be zero
R = Regas/M #Gas constant for the flow
Tci = Tamb #Fuel holding temperature is equal to ambient temp
Vinj = 0 #Injector velocity is zero
After determining the inlet mach number, several other parameters may also be obtained.
The nozzle stagnation pressure can be found as a function of the injector pressure, gamma, and the inlet mach number:
$$ \frac{P_{c_{inj}}}{P_{c_{ns}}} = \frac{\big( 1 + \gamma M^2_i\big)}{\bigg( 1 + \frac{\gamma - 1}{2} M^2_i \bigg)^\frac{\gamma}{\gamma-1}} $$(equation 1-14, huang)
Which may be used to find the nozzle pressure at the inlet:
$$ \frac{P_{inj}}{P_i} = 1 + \gamma M^2_i $$(equation 1-15, huang)
Further, the pressure at the throat, also known as the critical pressure, may also be determined:
$$ \frac{P_t}{P_{c_{ns}}} = \bigg( \frac{2}{\gamma + 1} \bigg)^{\frac{\gamma}{\gamma - 1}} $$(equation 1-16, huang)
It is important to note that past the throat of a de Laval nozzle, because the flow is moving faster than a pressure wave will propegate; a pressure wave due to an ambient pressure higher than that inside the nozzle will not occur. Because of this it is possible to expand the flow to pressures lower than ambient (or over-expand). Peak force production occurs when this pressure difference is negligible, and it is therefore important to consider at which altitude a nozzle should produce its peak force. It is uncommon for seperation to occur at the nozzle wall due to this overexpansion, but severe expansion angles or drastic over-expansion can lead to flow seperation.
Velocity is obtained by the conversion of propellant enthalpy to motion, exit velocity is a function of the ratio of specific heats, nozzle pressure ratio, and chamber temperature, and is represented by,
$$ V_{e} = \sqrt{\frac{2g\gamma}{\gamma-1}R T_{c_{ns}}\bigg[1-\bigg(\frac{P_{e}}{P_{c_{ns}}}\bigg)^\frac{\gamma-1}{\gamma}}\bigg] $$(1-18, huang)
The velocity at the throat is independent of pressure, and represents the speed of sound in the particular medium at the throat; it is a function of the ratio of specific heats and nozzle stagnation temperature,
$$ V_t = \sqrt{\frac{2g\gamma}{\gamma + 1}RT_{c_{ns}}} $$(1-22, huang)
Mass flow rate will depend on the desired force for a nozzle design. In order to provide a specific force requirement, and given the above calculated exit velocities; a mass flow rate is calculated,
$$ F = \frac{\dot{w}}{g}V_e $$(1-6, huang)
The propellant composition determines the specific mass flow rates of fuel and oxidizer. The fuel mass flow rate is calculated,
(mass percent formula for fuel * flow rate)
The difference between fuel and total flow rate is the remaining oxidizer flow rate,
(ox flow rate formula)
The throat temperature is a function of either the throat pressure ratio, or throat velocity ratio and ratio of specific heats and is calculated,
$$ \frac{T_i}{T_t} = \bigg( \frac{P_i}{P_t} \bigg)^{\frac{\gamma - 1}{\gamma}} = \bigg( \frac{V_t}{V_i} \bigg)^{\gamma - 1} $$(1-13, huang)
The exit temperature is calculated similarly,
$$ \frac{T_i}{T_e} = \bigg( \frac{P_i}{P_e} \bigg)^{\frac{\gamma - 1}{\gamma}} = \bigg( \frac{V_e}{V_i} \bigg)^{\gamma - 1} $$(1-13, huang)
The inlet temperature is a function of the inlet mach number and ratio of specific heats,
$$ T_i = \frac{T_{c_{ns}}}{\big( 1+\frac{\gamma - 1}{2}M^2_i \big)} $$(table 1-1, huang)
The throat area is critical to a choked flow. It is determined by first calculating the volume of gas according to the perfect gas law,
$$ 144 P_xV_x=RT_x $$(1-9, huang)
which is then applied with the mass flow rate and throat velocity.
$$ \dot{W} = \frac{A_iv_i}{144V_i} = \frac{A_xv_x}{144V_x} $$(1-11, huang)
Further, the expansion ratio is represented as a function of the exit ratio and ratio of specific heats,
$$ \epsilon = \frac{A_e}{A_t} = \frac{\bigg( \frac{2}{\gamma + 1} \bigg)^{\frac{1}{\gamma - 1}} \bigg[ \frac{P_{c_{ns}}}{P_e} \bigg]^{\frac{1}{\gamma}}}{\sqrt{ \frac{\gamma + 1}{\gamma - 1} \bigg[ 1 - \big( \frac{P_e}{P_{c_{ns}}} \big)^{\frac{\gamma - 1}{\gamma}} \bigg] }} $$(1-20, huang)
which may be applied to the area ratio to determine the requisite exit area,
$$ A_e = A_t \frac{\bigg( \frac{2}{\gamma + 1} \bigg)^{\frac{1}{\gamma - 1}} \bigg[ \frac{P_{c_{ns}}}{P_e} \bigg]^{\frac{1}{\gamma}}}{\sqrt{ \frac{\gamma + 1}{\gamma - 1} \bigg[ 1 - \big( \frac{P_e}{P_{c_{ns}}} \big)^{\frac{\gamma - 1}{\gamma}} \bigg] }} $$(1-20, huang)
Finally, the effective exhaust velocity is now calculated for engines operating at ambient pressures which are not equal to the exit pressure.
$$ c = v_e + A_e \big( P_e - P_a \big) \frac{g}{\dot{W}} $$(1-8, huang)
#Calculations
Pcns = Pcinj*(1+((gam-1)/2)*Mi**2)**(gam/(gam-1))/(1+gam*Mi**2) #Nozzle stagnation pressure (1-14, huang)
Pi = Pinj/(1+gam*Mi**2) #Pressure at inlet (1-15, huang)
Pt = Pcns*(2/(gam+1))**(gam/(gam-1)) #Pressure at throat (1-21, huang)
Ve = math.sqrt(((2*g*gam)/(gam-1))*R*Tcns*(1-(Pe/Pi)**((gam-1)/gam))) #Velocity at exit (1-18, huang)
Vt = math.sqrt(((2*g*gam)/(gam+1))*R*Tcns) #Velocity at throat (1-22, huang)
wdot = F*g/Ve #Propellant flow rate (1-6, huang)
wdotf = 1/(1+rw)*wdot #Fuel mass flow rate
wdoto = wdot-wdotf #Ox mass flow rate
Tt = Tcns*(Pt/Pcns)**((gam-1)/gam) #Temperature at throat (1-13, huang)
Ti = Tcns/(1+.5*(gam-1)*Mi**2) #Inlet temperature (Tab. 1-1, huang)
Te = Ti/(Pi/Pe)**((gam-1)/gam) #Temperature at exit (1-13, huang)
Volumet = R*Tt/144/Pt #Gas volume at throat (1-9, huang)
At = 144*wdot*Volumet/Vt #Area of throat (1-11, huang)
Wdot = At*Pcns*math.sqrt(g*gam*(2/(gam+1))**((gam+1)/(gam-1))/R/Tcns) #Theoretical weight flow (1-19, huang)
epsilon = ((2/(gam+1))**(1/(gam-1))*(Pcns/Pe)**(1/gam)/math.sqrt((gam+1)/
(gam-1)*(1-(Pe/Pcns)**((gam-1/gam))))) #Expansion ratio (1-20, huang)
Ae = At*epsilon #Area of exit (1-20, huang)
c = Ve+Ae*(Pe-Pa)*(g/wdot) #Effective exhaust velocity (1-8, huang)
#Print results
print("Pressures:")
print("Total Chamber, Pcinj (psi): %.2f" % Pcinj )
print("Stagnation, Pcns (psi): %.2f" % Pcns )
print("Inlet, Pi (psi): %.2f" % Pi )
print("Throat, Pt (psi): %.2f\n" % Pt )
print("Velocities:")
print("Throat, Vt (ft/sec): %.2f" % Vt )
print("Exit, Vt (ft/sec): %.2f" % Ve )
print("Effective, c (ft/sec): %.2f\n" % c )
print("Mass Flow Rates:")
print("Total, wdot (lb/sec): %.2f" % wdot )
print("Fuel, wdotf (lb/sec): %.2f" % wdotf )
print("Oxidizer, wdoto (lb/sec): %.2f" % wdoto )
print("Theoretical, Wdot (lb/sec): %.2f\n" % Wdot )
print("Temperatures:")
print("Fuel (holding), Tci (R): %.2f" % Tci )
print("Stagnation, Tcns (R): %.2f" % Tcns )
print("Inlet, Ti (R): %.2f\n" % Ti )
Pressures: Total Chamber, Pcinj (psi): 350.00 Stagnation, Pcns (psi): 349.33 Inlet, Pi (psi): 348.65 Throat, Pt (psi): 202.73 Velocities: Throat, Vt (ft/sec): 3464.05 Exit, Vt (ft/sec): 7804.45 Effective, c (ft/sec): 7804.45 Mass Flow Rates: Total, wdot (lb/sec): 2.06 Fuel, wdotf (lb/sec): 0.94 Oxidizer, wdoto (lb/sec): 1.13 Theoretical, Wdot (lb/sec): 2.06 Temperatures: Fuel (holding), Tci (R): 551.00 Stagnation, Tcns (R): 5487.10 Inlet, Ti (R): 5485.96
Material properties are dependent on the material and process used. Current listed properties represent materials available through i3D manufacturing processes according to their print datasheets. Values do not typically deviate to a magnitude which would result in unacceptable designs if different print materials were to be used, but it is important to check any design against material datasheet properties which represent the final process which will be used to machine the engine.
Properties are represented by the following variables: a: Thermal expansion coefficent (in/in-F) E: Modulus of elasticity (lb/in^2) k: Thermal conductivity (Btu/in^2-s-F/in) Twg: Maximum acceptable material temperature (R) v: Poisson's ratio (in/in) RA: RA surface roughness (in)
Alteration of the material index variable (MatInd) is the only change necessary to choose a material for nozzle construction.
#material properties
# Typical values for various printable materials
# Inconel 718: a = 0.000008065 , E = 25450000 , k = .000225694 , Twg = 1700 , v = 0.274 , RA = 650*10**-6 in
# Aluminum (AlSi10Mg): a = 0.000014 , E = 10200000 , k = 0.001716667 - 0.001983333, Twg = 1500, v = 0.33 , RA = 0.59-0.75 * 10**-3 in
# Cobalt Chrome: a = 0.0000069 , E = 30000000 , k = 0.000333333 , Twg = 2563 , v = 0.29 , RA = 0.39 *10**-3 in
# Titanium: a = 0.000032724 , E =16100000 , k = .000111667 , Twg = 3479 , v = 0.31, RA
# Steel
Props = [
# a E k Twg v RA
["0.000008065" , "25450000" , "0.000225694" , "1700" , "0.274" , "0.00065"], # 0 Inconel
["0.000014" , "10200000" , "0.001716667" , "1200" , "0.33" , "0.00059"], # 1 Aluminum
["0.0000069" , "30000000" , "0.000333333" , "2563" , "0.29" , "0.00039"], # 2 Cobalt Chrome
["0.000032724" , "16100000" , "0.000111667" , "3479" , "0.31" , "0.00065"] # 3 Titanium
]
MatInd = 1 # Material index
a_ = widgets.Text(Props[MatInd][0], description="Thermal expansion ratio nozzle material" , width=100 )
E_ = widgets.Text(Props[MatInd][1], description="Elastic modulus of nozzle material, E (psi)" , width=70 )
k_ = widgets.Text(Props[MatInd][2], description="Wall thermal conductivity, k (Btu-in/in^2-s-F)" , width=120 )
Twg_ = widgets.Text(Props[MatInd][3], description="Maximum Material Temperature at wall, (deg R)" , width=60 )
v_ = widgets.Text(Props[MatInd][4], description="Poisson's ratio nozzle material" , width=60 )
RA_ = widgets.Text(Props[MatInd][5], description="Surface Rougness, (inch RA)" , width=100 )
#display user input boxes
display(a_ )
display(E_ )
display(k_ )
display(Twg_ )
display(v_ )
display(RA_ )
#Convert user inputs to float
a = float(a_.value ) #Thermal expansion ratio of nozzle material
E = float(E_.value ) #Elastic modulus of nozzle material
k = float(k_.value ) #Wall thermal conductivity
Twg = float(Twg_.value ) #Maximum Wall Temperature
v = float(v_.value ) #Poisson's ratio nozzle material
RA = float(RA_.value ) #RA surface roughness of nozzle material
Certain variables are dependant on design choices, such as a choice of fuel, or a desired cooling channel area ratio. These variables may be used to manipulate resulting design dependent values.
High aspect ratios have been associated with higher cooling efficiencies (cite paper xx). Increasing the cooling channel aspect ratio will result in a smaller hydraulic diameter, and greater pressure losses through the cooling channels. Increases in aspect ratio will also lead to an increased number of cooling channels, which decreases stresses in the chamber walls due to a pressure differential on either side of a given wall. Typically, small scale rockets will not experience significant stresses due to pressure within the cooling channel, as the pressures required to cause significant stresses within such small geometries are high enough to require turbomachinery, which are generally not available to builds on the amature scale. Lowering aspect ratios can therefore be used to decrease pressure losses and make pressure fed systems more manageable.
Fuel boiling point is a material property, but may be manipulated in certain fuels, particularly those which are aqueous solutions. Dilution of fuels with water will generally change the boiling point of a fuel, though the effect is generally not linear. Boiling point is taken at standard temperature and pressure for a given fuel solution.
The specific heat of the fuel is a material property. The value used for specific heat should be that of a pure sample of the fuel.
Injector design is beyond the scope of this document, but the pressure drop across any particular design must be accounted for, and is done so here.
Heat of vaporization is a material property. In order to determine coolant capacity, the boiling point of the coolant must be calculated at the lowest pressure along the cooling channel. The heat of vaporization is used to calculate this value.
The characteristic length of the chamber is a range of emperically determined values for various propellant combinations. Accepted ranges for characterist length for various fuels can be found (consider adding table xx) in table 4-1 in Huang.
Factor of safety to protect against nucleate boiling. For pressures which will not allow for a supercritical coolant, nucleate boiling should be avoided. Designers should determine an acceptable level of risk when approaching the boiling point of the coolant. High factors of safety should also be avoided, as over cooling an engine can decrease performance significantly. Typical factors of safety range from 1.1 to 1.3 for engine applications.
The pressure of the fuel when it arrives at the fuel manifold. Higher pressures can protect against nucleate boiling by raising the temperature required for nucleate boiling.
CEArun parameter. Can be manipulated slightly by changing composition of the fuel or mixture ratio.
Determined experimentally. Adjust after measuring engine performance in testing.
When film cooling is used in conjunction with regenerative cooling, a correction factor for the adiabadic gas wall temperature must be applied. There is no reliable model for determining this, but a correction factor can be determined after measuring the temperature rise in the coolant during hot fire.
#Define User Inputs for Equation Variables, assign values to default values
ARhb_ = widgets.Text("12.5" , description="Cooling channel aspect ratio at throat" , width=60 )
Bpf_ = widgets.Text("688.67" , description="Boiling Point of Fuel, Bpf (deg R)" , width=100 )
CpFuel_ = widgets.Text("0.5876" , description="Specific heat of pure Fuel, CpFuel (BTU/lb-F)" , width=80 )
DelPi_ = widgets.Text("87.5" , description="Injector pressure drop, DelPi (lbf/in^2)" , width=60 )
Hvapf_ = widgets.Text("36800" , description="Fuel heat of vaporization, Hvapf (J/mol)" , width=80 )
Lstar_ = widgets.Text("28" , description="Characteristic Chamber Length, L* (in) " , width=60 )
nc_ = widgets.Text("1.075" , description="Nucleate boiling factor of safety, nc" , width=60 )
Pco_ = widgets.Text("550" , description="Initial fuel pressure, Pco (lbf/in^2)" , width=60 )
Tt_ = widgets.Text("5232.4" , description="Throat Temperature (deg R)" , width=80 )
eta_ = widgets.Text("0.90" , description="Combustion efficiency" , width=60 )
edafilm_ = widgets.Text("0.8" , description="Correction factor from film cooling for Taw" , width=60 )
#Display user input text boxes
display(ARhb_ )
display(Bpf_ )
display(CpFuel_ )
display(DelPi_ )
display(Hvapf_ )
display(Lstar_ )
display(nc_ )
display(Pco_ )
display(Tt_ )
display(eta_ )
display(edafilm_ )
#Convert string entries to floats
ARhb = float(ARhb_.value ) #Cooling channel aspect ratio
Bpf = float(Bpf_.value ) #Boiling Point of fuel
CpFuel = float(CpFuel_.value ) #Specific Heat of pure fuel
DelPi = float(DelPi_.value ) #Pressure drop across injector
Hvapf = float(Hvapf_.value ) #Fuel heat of vaporization
Lstar = float(Lstar_.value ) #Characteristic Chamber Length
nc = float(nc_.value ) #Nucleate boiling factor of safety
Pco = float(Pco_.value ) #Initial fuel pressure
Tt = float(Tt_.value ) #Temperature of throat flow
eta = float(eta_.value ) #Combustion efficiency
edafilm = float(edafilm_.value ) #Adiabadic wall temperature correction factor due to film cooling
#Performance Parameters
Is = F/wdot #Specific Impulse (1-28, huang)
Istc = c/g #Thrust chamber specific impulse (1-31a, huang)
cstar = (math.sqrt(g*gam*R*Tcns)/gam
/math.sqrt((2/(gam+1))**((gam+1)/(gam-1)))) #Characteristic Velocity (1-32a, huang)
Cf = (math.sqrt(2*gam**2/(gam-1)*(2/(gam+1))**((gam+1)/(gam-1))
*(1-(Pe/Pcns)**((gam-1)/gam)))+epsilon*((Pe-Pa)/Pcns)) #Thrust Coefficient (1-33a, huang)
#Print Performance Parameters
print("Specific Impulse, Is (lb s/lb): %.2f" % Is )
print("Thrust Chamber specific impulse, Istc (lb s/lb): %.2f" % Istc )
print("Characteristic Velocity, c* (ft/s): %.2f" % cstar )
print("Thrust Coefficient, Cf: %.2f" % Cf )
Specific Impulse, Is (lb s/lb): 242.37 Thrust Chamber specific impulse, Istc (lb s/lb): 242.37 Characteristic Velocity, c* (ft/s): 5325.59 Thrust Coefficient, Cf: 1.47
#Thrust Chamber Layout
Vc = Lstar*At #Chamber volume (4-4, huang)
Ac = epsilonc*At #Chamber cross sectional area (fig. 4-11, huang)
Lc = Lstar/Ac #does not account for the converging part of combustion chamber use 4-5, huang instead
#print thrust chamber layout
print("Chamber Volume, Vc (in^3): %.2f" % Vc)
print("Chamber Cross Sectional Area, Ac (in^2): %.2f" % Ac)
print("Chamber Length, Lc (in): %.2f" % Lc)
Chamber Volume, Vc (in^3): 27.35 Chamber Cross Sectional Area, Ac (in^2): 10.26 Chamber Length, Lc (in): 2.73
#Heat Transfer
CpH20 = 0.998137 #Water specific heat
Cplc = rwfuel/100*CpFuel+(1-rwfuel/100)*CpH20 #Coolant specific heat
Pr = 4*gam/(9*gam-5) #Prandtl number
mucc = (46.6*10**-10)*M**0.5*Tcns #Viscosity in the combustion chamber
mut = (46.6*10**-10)*M**0.5*Tt #Viscosity in the throat
rlam = Pr**0.5 #Laminar flow local recovery factor
rturb = Pr**0.33 #Turbulent flow local recovery factor
Reffcc = ((1+rturb*((gam-1)/2)*Mi**2)/(1+((gam-1)/2)*Mi**2))
Refft = ((1+rturb*((gam-1)/2))/(1+((gam-1)/2)))
Tawi = edafilm*Tcns*Reffcc #Adiabatic wall temperature at inlet
Tawt = edafilm*Tcns*Refft #Adiabatic wall temperature at throat
Tcc = 9/5*(math.log(Pa/(Pi+DelPi))*Rgas/Hvapf+1/(Bpf*5/9))**-1 #Critical temperature of fuel coolant
Twc = Tcc/nc #Maximum coolant wall temperature
Cpg = gam*R/(gam-1)/J #Specific heat at constant pressure
rt = math.sqrt(At/math.pi) #Radius of throat
re = math.sqrt(Ae/math.pi) #Radius of exit
rmean = rt*(1.5+.382)/2 #Mean throat curvature
sigmat = (1/((.5*Twg/Tcns*(1+(gam-1)/2)+.5)**0.68 #Correction factor for property variations across BL
*(1+(gam+1)/2)**0.12)) #specified at throat
sigmai = (1/((.5*Twg/Tcns*(1+(gam-1)/2*Mi**2)+.5)**0.68 #Correction factor for property variations across BL
*(1+(gam+1)/2*Mi**2)**0.12)) #specified at inlet
hg = ((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)
*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*sigmat) #heat transfer coefficient at throat
#film cooling estimate
Cpvc = 0.582784 #btu/lb/f
edac = .5
Gc = (hg/Cpvc/edac)/math.log((Tawt-Twc)/(Tawt-Twg))
Acool = 1 # Desired cooled surface area (in^2)
q = hg*(Tawt-Twg) #required heat flux
Tbulk = (Twc + Tci)/2 #Coolant bulk temp
t = k/q*(Twg-Twc) #Calculated Wall thickness for desired coolant wall temp
Qc = (wdotf+Gc*Acool)*Cplc*(Twc-Tci) #Coolant capacity
hc = q/(Twc-Tbulk) #Coolant side heat transfer coefficient
H = 1/(1/hg+t/k+1/hc) #Overall heat transfer coefficient
#print heat transfer
print("Coolant Specific Heat, Cplc: %.2f" % Cplc)
print("Prandtl Number, Pr: %.2f" % Pr)
print("Viscosity in combustion chamber, mucc: %.8f" % mucc)
print("Viscosity in throat, mut: %.8f" % mut)
print("Laminar Flow Local Recovery Factor, rlam: %.2f" % rlam)
print("Turbulent Flow Local Recovery Factor, rturb: %.2f" % rturb)
print("Effective Combustion Chamber Recovery Factor, Reffcc: %.2f" % Reffcc)
print("Effective Throat Recovery Factor, Refft: %.2f" % Refft)
print("Adiabatic wall temperature at inlet, Tawi: %.2f" % Tawi)
print("Adiabatic wall temperature at throat, Tawt: %.2f" % Tawt)
print("Throat radius, rt: %.3f" % rt)
print("Exit radius, re: %.3f" % re)
print("Mean throat curvature, rmean: %.3f" % rmean)
print("Correction factor for property variations across BL at throat, sigmat: %.2f" % sigmat)
print("Correction factor for property variations across BL at inlet, sigmai: %.2f" % sigmai)
print("Specific heat at constant pressure, Cpg: %.5f" % Cpg)
print("Heat transfer coefficient at throat, hg: %.5f" % hg)
print("Required heat flux, q: %.2f" % q)
print("Critical temperature of coolant, Tcc: %.2f" % Tcc)
print("Coolant capacity, Qc: %.2f" % Qc)
print("Coolant wall temperature, Twc %.2f" % Twc)
print("Coolant bulk temperature, Tbulk: %.2f" % Tbulk)
print("Wall thickness, t: %.5f" % t)
print("Coolant side heat transfer coefficient, hc: %.5f" % hc)
print("Overall heat transfer coefficient, H: %.5f" % H)
print("Film-Coolant flowrate per unit area of cooled wall, lb/in^2/s Gc: %.3f" % Gc)
Coolant Specific Heat, Cplc: 0.75 Prandtl Number, Pr: 0.88 Viscosity in combustion chamber, mucc: 0.00012534 Viscosity in throat, mut: 0.00011952 Laminar Flow Local Recovery Factor, rlam: 0.94 Turbulent Flow Local Recovery Factor, rturb: 0.96 Effective Combustion Chamber Recovery Factor, Reffcc: 1.00 Effective Throat Recovery Factor, Refft: 1.00 Adiabatic wall temperature at inlet, Tawi: 4389.64 Adiabatic wall temperature at throat, Tawt: 4379.47 Throat radius, rt: 0.558 Exit radius, re: 1.037 Mean throat curvature, rmean: 0.525 Correction factor for property variations across BL at throat, sigmat: 1.27 Correction factor for property variations across BL at inlet, sigmai: 1.40 Specific heat at constant pressure, Cpg: 0.76629 Heat transfer coefficient at throat, hg: 0.00872 Required heat flux, q: 27.72 Critical temperature of coolant, Tcc: 974.12 Coolant capacity, Qc: 341.68 Coolant wall temperature, Twc 906.16 Coolant bulk temperature, Tbulk: 728.58 Wall thickness, t: 0.01820 Coolant side heat transfer coefficient, hc: 0.15609 Overall heat transfer coefficient, H: 0.00759 Film-Coolant flowrate per unit area of cooled wall, lb/in^2/s Gc: 0.338
import sympy as sympy
from sympy.solvers import solve
from sympy import Symbol
from sympy import log
#######These should be in a reasonable input section, not burried in code.#########
Re = 3200 #Desired channel Reynold's number
muf = 2.7998705 *10**-5 #Viscosity of the fuel mixture
rhof = 0.031972653 #density of 70% ethanol (lb per cubic inch)
#solve for minimum number of cooling channels to acheive appropriate cross sectional area,
#and the appropriate channel aspect ratio
n_g = 0 #a reasonable guess for a particular geometry's number of cooling channels
n = 1
while (abs(n_g-n)>0):
n_g = n
bt = (2*(rt+t)*pi/n_g-t)
d = 2*ARhb*bt**2/(ARhb+1)/bt
Vco = Re/rhof/d*muf #Velocity for chosen Reynold's Number
A_u = (wdotf+Gc*Acool)/rhof/Vco #unobstructed cross sectional area
#find min number of channels for ideal ratio
L = Symbol('L')
n = round(solve(L*t+(A_u*L/ARhb)**0.5-(2*rt+2*t)*math.pi,L)[0])
Vco = Re/rhof/d*muf #velocity for chosen Reynold's Number
A_u = (wdotf+Gc*Acool)/rhof/Vco #unobstructed cross sectional area
#cooling channel geometry
bt = (2*rt+t)*math.pi/n-t #base width at throat
be = (2*re+t)*math.pi/n-t #base width at exit
rccht = bt/2 #cooling channel effective radius at throat
rcche = be/2 #cooling channel effective radius at nozzle exit
d = 2*ARhb*bt**2/(ARhb+1)/bt #cooling channel hydraulic diameter
A_ob = n*(bt**2*(1-math.pi/4)+ARhb*bt*t) #total obstructed area
A_c = A_u + A_ob #Total area for cooling channels to be used in parametric equations
#friction factor
if Re > 4000:
f=0.2
F=0
while (abs(f-F)>0.000000000001):
F=(f+F)/2
f=(1/(-2*math.log10(RA/3.71/d+2.51/Re/F**0.5)))**2
else:
f=64/Re
#calculate pressure drop of a cooling channel
#to do this we must:
#calculate arc length of parametric equations (or use solidworks for now...)
#calculate reynolds number of passage
#calculate friction factor for the pipe use solver of darcy FF
#use this to determine holding tank pressure
#print cooling channel geometry
print("Fuel viscosity, muf: %.8f" % muf)
print("Fuel density, rhof: %.5f" % rhof)
print("Velocity for turbulent flow, Vco: %.2f" % Vco)
print("Base width at throat, bt: %.5f" % bt)
print("Base width at exit, be: %.5f" % be)
print("Cooling channel effective radius at throat, rccht: %.5f" % rccht)
print("Cooling channel effective radius at exit, rcche: %.5f" % rcche)
print("Cooling channel hydraulic diameter, d: %.5f" % d)
print("Cooling channel unobstructed area, A_u: %.5f" % A_u)
print("Cooling channel obstructed area, A_ob: %.5f" % A_ob)
print("Cooling channel total area, A_c: %.5f" % A_c)
print("Number of cooling channels, n: %.0f" % n)
print("Friction factor of cooling channel, f: %.5f" % f)
Fuel viscosity, muf: 0.00002800 Fuel density, rhof: 0.03197 Velocity for turbulent flow, Vco: 58.38 Base width at throat, bt: 0.02522 Base width at exit, be: 0.06192 Cooling channel effective radius at throat, rccht: 0.01261 Cooling channel effective radius at exit, rcche: 0.03096 Cooling channel hydraulic diameter, d: 0.04671 Cooling channel unobstructed area, A_u: 0.68368 Cooling channel obstructed area, A_ob: 0.48169 Cooling channel total area, A_c: 1.16537 Number of cooling channels, n: 82 Friction factor of cooling channel, f: 0.02000
#wall stresses
Ste =(Pco-Pe)*rcche/t+E*a*q*t/2/(1-v)/k #combined tangential stress at nozzle exit
Stt =(Pco-Pt)*rccht/t+E*a*q*t/2/(1-v)/k
Sce =(Pco-Pe)*re/t+E*a*q*t/2/(1-v)/k #maximum compressive stress as coaxial shell design
#Print wall stresses
print("Combined tangential stress at nozzle exit, Ste: %.2f" % Ste)
print("Combined tangential stress at throat, Stt: %.2f" % Stt)
print("Maximum compressive stress as coaxial shell design, Sce: %.2f" % Sce) #does not account for support of channel walls
Combined tangential stress at nozzle exit, Ste: 32224.43 Combined tangential stress at throat, Stt: 31554.40 Maximum compressive stress as coaxial shell design, Sce: 61802.33
#heat transfer plots
%matplotlib inline
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
#this needs to be revisited to ensure adjustments are being made for film cooling, the temperature doesnt seem to be graphing quite correctly
x = np.linspace(Tbulk,Twg)
#figure 1 Max stress v. Gas side wall temp
plt.figure(1)
line, = plt.plot(x,(Pco-Pt)*rccht/(k/((((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*((1/((.5*x/Tcns*(1+(gam-1)/2)+.5)**0.68*(1+(gam+1)/2)**0.12)))))*(Tawt-x))*(x-Tcc))+E*a*q*(k/((((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*((1/((.5*x/Tcns*(1+(gam-1)/2)+.5)**0.68*(1+(gam+1)/2)**0.12)))))*(Tawt-x))*(x-Tcc))/2/(1-v)/k, '--', linewidth=2)
plt.suptitle("Figure 1")
plt.title("Max stress v. Gas-side wall temp")
plt.xlabel("Temp (R)")
plt.ylabel("Stress (psi)")
plt.ylim([0,30000])
plt.xlim([Twc,Twg])
#figure 2 Wall thickness v. gas side wall temp
plt.figure(2)
line, = plt.plot(x,k/((((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*((1/((.5*x/Tcns*(1+(gam-1)/2)+.5)**0.68*(1+(gam+1)/2)**0.12)))))*(Tawt-x))*(x-Tcc), '--', linewidth=2)
plt.suptitle("Figure 2")
plt.title("Wall thickness v. Gas-side wall temp")
plt.xlabel("Temp (R)")
plt.ylabel("Thickness (in)")
plt.ylim([0,0.020])
plt.xlim([Twc,Twg])
(906.1595326365948, 1200.0)
Nozzle geometry is a critical perameter of any engine. Given an throat area, contraction ratio, expansion ratio, as well as the rapid expansion and contraction angles surrounding the throat; a unique set of parametric equations will justify the inner nozzle contour. Furthermore, given the wall thickness calculated as a function of the heat transfer properties of the throat; there exists a unique set of parametric equations justifying the remaining cooling channel geometries. These contours are generated such that the cross sectional area of the cooling channels normal to the surface of the rocket is a constant. Scaling functions may be applied, with an inversely proportional effect on the velocity of the fluid through the cooling channels. Left unchanged, the velocity of the coolant will be highest at the throat. The cross sectional area of the channels is certainly uniform, but the hydraulic diameter is decreasing as the cooling channels increase in aspect ratio (defined as the height of the rectangular duct divided by the base).
A designer may utilize these parametric equations by importing them into computer aided drafting (CAD) software, and rotating the resulting contours about the y axis. The particular proceedure of such a rotation in any particular CAD software is beyond the scope of this document.
rcc = math.sqrt((At*epsilonc)/math.pi)
#Parametric Equations
#calculate parameters
tout = 2*t # calculate a better outer thickness or place in input parameters
theta = 45
thetaN = 45
A1 = math.tan(math.radians(90-thetaN))/(2*(rt-rt*.382*(1-math.cos(math.radians(thetaN)))))
A2 = -(.382*rt*math.sin(math.radians(thetaN)))+A1*(rt+.382*rt*(1-math.cos(math.radians(thetaN))))**2
X1 = rt+.382*rt*(1-math.cos(math.radians(thetaN)))
print("Parametric Equations:\n")
print("Internal Contour:\n")
print("Bell")
print("x(t):")
print(" t")
print("y(t):")
print(" %f * t^2 - %f" % (A1,A2) )
print("Parameters:" )
print("t1:")
print(" %f" % X1 )
print("t2:")
print(" %f\n" % (re-0.0001) )
Rbell = t
A3 = 1.382*rt
A4 = 0.382*rt
print("Throat, rapid expansion")
print("x(t):")
print(" %f * (1- (.382/1.382) * cos(t*3.14159/180))" % A3 )
print("y(t):")
print(" %f * sin(t*3.14159/180)" % A4 )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (thetaN-0.001) )
A5 = 2.5*rt
A6 = -1.5*rt
print("Throat, rapid contraction")
print("x(t):")
print(" %f * (1 - 0.6 * cos(t*3.14159/180))" % A5 )
print("y(t):")
print(" %f * sin(t*3.14159/180)" % A6 )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.001) )
A7 = -math.tan(math.pi/2-math.radians(theta))
A8 = (math.tan(math.pi/2-math.radians(theta))*2.5*rt*(1-0.6*math.cos(math.radians(theta))))-1.5*rt*math.sin(math.radians(theta))
p1 = 2.5*rt*(1-0.6*math.cos(math.radians(theta)))
p2 = rcc-1.5*rt*(1-math.cos(math.radians(theta)))
print("Linear contraction")
print("x(t):")
print(" t")
print("y(t):")
print(" %f * t + %f" % (A7,A8))
print("Parameters:" )
print("t1:")
print(" %f" % (p1+0.0001) )
print("t2:")
print(" %f\n" % (p2-0.0001) )
A9 = 1.5*rt
A10 = 1.5*rt
A11 = (math.tan(2*math.pi-(math.pi/2-math.radians(theta))))*(rcc-1.5*rt*(1-math.cos(math.radians(theta))))+(-math.tan(2*math.pi-(math.pi/2-math.radians(theta)))*(2.5*rt*(1-0.6*math.cos(math.radians(theta))))-1.5*rt*math.sin(math.radians(theta)))-1.5*rt*math.sin(math.radians(theta))
print("Inlet")
print("x(t):")
print(" %f - %f * (1 - cos(t*3.14159/180))" % (rcc,A9))
print("y(t):")
print(" %f * sin(t*3.14159/180) + %f" % (A10,A11))
print("Parameters:")
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.0001) )
Parametric Equations: Internal Contour: Bell x(t): t y(t): 1.009704 * t^2 - 0.237475 Parameters: t1: 0.619964 t2: 1.036417 Throat, rapid expansion x(t): 0.770575 * (1- (.382/1.382) * cos(t*3.14159/180)) y(t): 0.212995 * sin(t*3.14159/180) Parameters: t1: 0.001 t2: 44.999000 Throat, rapid contraction x(t): 1.393948 * (1 - 0.6 * cos(t*3.14159/180)) y(t): -0.836369 * sin(t*3.14159/180) Parameters: t1: 0.001 t2: 44.999000 Linear contraction x(t): t y(t): -1.000000 * t + 0.211144 Parameters: t1: 0.802646 t2: 1.561697 Inlet x(t): 1.806764 - 0.836369 * (1 - cos(t*3.14159/180)) y(t): 0.836369 * sin(t*3.14159/180) + -1.942055 Parameters: t1: 0.001 t2: 44.999900
print("Cooling Channel, inner:\n" )
A12 = 1/(4*A1**2)
print("Bell")
print("x(t):")
print(" t + ( %f * t / ( t^2 + %f )^0.5)" % (t,A12))
print("y(t):")
print(" %f *t^2 - %f + ( %f / ( %f *( t^2 + %f )^0.5))" % (A1,A2,-t,2*A1,A12))
print("Parameters:")
print("t1:")
print(" %f" % (X1+0.0001) )
print("t2:")
print(" %f\n" % (re-0.0001) )
print("Rapid Expansion")
print("x(t):")
print(" %f * (1- (.382/1.382) * cos(t*3.14159/180)) + ( %f * cos(t*3.14159/180))" % (A3,t) )
print("y(t):")
print(" %f * sin(t*3.14159/180) + (%f *sin(t*3.14159/180))" % (A4,-t) )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (thetaN-0.001) )
print("Rapid Contraction")
print("x(t):")
print(" %f * (1 - 0.6 * cos(t*3.14159/180)) + ( %f * cos(t*3.14159/180))" % (A5,t) )
print("y(t):")
print(" %f * sin(t*3.14159/180) + (%f *sin(t*3.14159/180))" % (A6,t) )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.001) )
A13 = math.cos(math.radians(theta))*t
A14 = math.sin(math.radians(theta))*t
print("Linear contraction")
print("x(t):")
print(" t + %f" % A13)
print("y(t):")
print(" %f * t + %f" % (A7,A8+A14))
print("Parameters:" )
print("t1:")
print(" %f" % (p1+0.0001) )
print("t2:")
print(" %f\n" % (p2-0.0001) )
print("Inlet")
print("x(t):")
print(" %f - %f * (1 - cos(t*3.14159/180)) + ( %f * cos(t*3.14159/180))" % (rcc,A9,t))
print("y(t):")
print(" %f * sin(t*3.14159/180) + %f + (%f *sin(t*3.14159/180))" % (A10,A11,t))
print("Parameters:")
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.0001) )
Cooling Channel, inner: Bell x(t): t + ( 0.018199 * t / ( t^2 + 0.245218 )^0.5) y(t): 1.009704 *t^2 - 0.237475 + ( -0.018199 / ( 2.019409 *( t^2 + 0.245218 )^0.5)) Parameters: t1: 0.620064 t2: 1.036417 Rapid Expansion x(t): 0.770575 * (1- (.382/1.382) * cos(t*3.14159/180)) + ( 0.018199 * cos(t*3.14159/180)) y(t): 0.212995 * sin(t*3.14159/180) + (-0.018199 *sin(t*3.14159/180)) Parameters: t1: 0.001 t2: 44.999000 Rapid Contraction x(t): 1.393948 * (1 - 0.6 * cos(t*3.14159/180)) + ( 0.018199 * cos(t*3.14159/180)) y(t): -0.836369 * sin(t*3.14159/180) + (0.018199 *sin(t*3.14159/180)) Parameters: t1: 0.001 t2: 44.999000 Linear contraction x(t): t + 0.012868 y(t): -1.000000 * t + 0.224012 Parameters: t1: 0.802646 t2: 1.561697 Inlet x(t): 1.806764 - 0.836369 * (1 - cos(t*3.14159/180)) + ( 0.018199 * cos(t*3.14159/180)) y(t): 0.836369 * sin(t*3.14159/180) + -1.942055 + (0.018199 *sin(t*3.14159/180)) Parameters: t1: 0.001 t2: 44.999900
print("Cooling Channel, outer:\n" )
A15 = A_c/math.pi
print("Bell")
print("x(t):")
print(" t + ( %f * t / ( t^2 + %f )^0.5) + ( t / ( t^2 + %f )^0.5) * ( -t + (t^2 + %f * sin(2*3.14159 - arctan(%f*t)))^(.5)) / (sin(2*3.14159-arctan(%f*t)))" % (t,A12,A12,A15,-2*A1,-2*A1))
print("y(t):")
print(" %f *t^2 - %f + ( %f / ( %f *( t^2 + %f )^0.5)) + ( 1 / ( %f * (t^2+ %f )^0.5)) * ( -t + ( t^2 + %f * sin(2*3.14159-arctan(%f*t)))^0.5) / ( sin(2*3.14159 - arctan(%f*t)))" % (A1,A2,-t,2*A1,A12,-2*A1,A12,A15,-2*A1,-2*A1))
print("Parameters:")
print("t1:")
print(" %f" % (X1+0.0001) )
print("t2:")
print(" %f\n" % (re-0.0001) )
print("Rapid Expansion")
print("x(t):")
print(" %f * (1- (.382/1.382) * cos(t*3.14159/180)) + ( %f * cos(t*3.14159/180)) + ( -( %f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + (( %f * (1 - (0.382)/(1.382) * cos(t*3.14159/180) ))^2 + %f * cos(t*3.14159/180))^0.5)" % (A3,t,A3,A3,A15) )
print("y(t):")
print(" %f * sin(t*3.14159/180) + (%f *sin(t*3.14159/180)) + ( -( %f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + ( (%f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) * ( -tan(t*3.14159/180) )" % (A4,-t,A3,A3,A15) )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (thetaN-0.001) )
print("Rapid Contraction")
print("x(t):")
print(" %f * (1 - 0.6 * cos(t*3.14159/180)) + ( %f * cos(t*3.14159/180)) + ( -( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + ( %f ) * cos(t*3.14159/180) )^0.5 )" % (A5,t,A5,A5,A15) )
print("y(t):")
print(" %f * sin(t*3.14159/180) + (%f *sin(t*3.14159/180)) + ( -( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) * ( tan(t*3.14159/180) )" % (A6,t,A5,A5,A15) )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.001) )
A16 = math.cos(math.radians(theta))*A15
A17 = math.tan(math.pi/2-math.radians(theta))*2.5*rt*(1-0.6*math.cos(math.radians(theta)))
A18 = -1.5*rt*math.sin(math.radians(theta))
A19 = A15 * math.cos(math.radians(theta))
A20 = math.tan(math.radians(theta))
print("Linear contraction")
print("x(t):")
print(" t + %f + ( -t + ( t^2 + %f )^0.5 )" % (A13,A16))
print("y(t):")
print(" %f *t + %f + %f + ( -t + ( t^2 + %f )^0.5 ) * ( %f ) + ( %f )" % (A7,A17,A18,A19,A20,A14))
print("Parameters:" )
print("t1:")
print(" %f" % (p1+0.0001) )
print("t2:")
print(" %f\n" % (p2-0.0001) )
print("Inlet")
print("x(t):")
print(" %f - %f * (1 - cos(t*3.14159/180)) + ( %f * cos(t*3.14159/180)) + ( -( %f - %f * ( 1 - cos(t*3.14159/180) ) ) + ( ( %f - %f * ( 1 - cos(t*3.14159/180)))^2 + %f * cos(t*3.14159/180) )^0.5 )" % (rcc,A9,t,rcc,A9,rcc,A9,A15))
print("y(t):")
print(" %f * sin(t*3.14159/180) + %f + ( %f * sin(t*3.14159/180) ) + ( -( %f - %f * ( 1 - cos(t*3.14159/180) ) ) + ( ( %f - %f * ( 1 - cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) * tan(t*3.14159/180)" % (A10,A11,t,rcc,A10,rcc,A10,A15))
print("Parameters:")
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.0001) )
Cooling Channel, outer: Bell x(t): t + ( 0.018199 * t / ( t^2 + 0.245218 )^0.5) + ( t / ( t^2 + 0.245218 )^0.5) * ( -t + (t^2 + 0.370950 * sin(2*3.14159 - arctan(-2.019409*t)))^(.5)) / (sin(2*3.14159-arctan(-2.019409*t))) y(t): 1.009704 *t^2 - 0.237475 + ( -0.018199 / ( 2.019409 *( t^2 + 0.245218 )^0.5)) + ( 1 / ( -2.019409 * (t^2+ 0.245218 )^0.5)) * ( -t + ( t^2 + 0.370950 * sin(2*3.14159-arctan(-2.019409*t)))^0.5) / ( sin(2*3.14159 - arctan(-2.019409*t))) Parameters: t1: 0.620064 t2: 1.036417 Rapid Expansion x(t): 0.770575 * (1- (.382/1.382) * cos(t*3.14159/180)) + ( 0.018199 * cos(t*3.14159/180)) + ( -( 0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + (( 0.770575 * (1 - (0.382)/(1.382) * cos(t*3.14159/180) ))^2 + 0.370950 * cos(t*3.14159/180))^0.5) y(t): 0.212995 * sin(t*3.14159/180) + (-0.018199 *sin(t*3.14159/180)) + ( -( 0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + ( (0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) * ( -tan(t*3.14159/180) ) Parameters: t1: 0.001 t2: 44.999000 Rapid Contraction x(t): 1.393948 * (1 - 0.6 * cos(t*3.14159/180)) + ( 0.018199 * cos(t*3.14159/180)) + ( -( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + ( 0.370950 ) * cos(t*3.14159/180) )^0.5 ) y(t): -0.836369 * sin(t*3.14159/180) + (0.018199 *sin(t*3.14159/180)) + ( -( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) * ( tan(t*3.14159/180) ) Parameters: t1: 0.001 t2: 44.999000 Linear contraction x(t): t + 0.012868 + ( -t + ( t^2 + 0.262301 )^0.5 ) y(t): -1.000000 *t + 0.802546 + -0.591402 + ( -t + ( t^2 + 0.262301 )^0.5 ) * ( 1.000000 ) + ( 0.012868 ) Parameters: t1: 0.802646 t2: 1.561697 Inlet x(t): 1.806764 - 0.836369 * (1 - cos(t*3.14159/180)) + ( 0.018199 * cos(t*3.14159/180)) + ( -( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) ) + ( ( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180)))^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) y(t): 0.836369 * sin(t*3.14159/180) + -1.942055 + ( 0.018199 * sin(t*3.14159/180) ) + ( -( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) ) + ( ( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) * tan(t*3.14159/180) Parameters: t1: 0.001 t2: 44.999900
print("Nozzle surface, outer:\n" )
print("Bell")
print("x(t):")
print(" t + ( %f * t / ( t^2 + %f )^0.5 ) + ( t / ( t^2 + %f )^0.5 ) * ( -t + ( t^2 + %f * sin(2*3.14159 - arctan(%f*t)) )^0.5 ) / ( sin(2*3.14159 - arctan(%f*t)) ) + ( %f * t / ( t^2 + %f )^0.5 )" % (t,A12,A12,A15,-2*A1,-2*A1,tout,A12))
print("y(t):")
print(" %f *t^2 - %f + ( %f / ( %f *( t^2 + %f )^0.5)) + ( 1 / ( %f * ( t^2+ %f )^0.5 ) ) * ( -t + ( t^2 + %f * sin(2*3.14159 - arctan(%f*t)) )^0.5 ) / ( sin(2*3.14159 - arctan(%f*t)) ) + ( %f / ( %f * ( t^2 + %f )^0.5 ) )" % (A1,A2,-t,2*A1,A12,-2*A1,A12,A15,-2*A1,-2*A1,-tout,2*A1,A12))
print("Parameters:")
print("t1:")
print(" %f" % (X1+0.0001) )
print("t2:")
print(" %f\n" % (re-0.0001) )
print("Rapid Expansion")
print("x(t):")
print(" %f * ( 1- (.382/1.382) * cos(t*3.14159/180) ) + ( %f * cos(t*3.14159/180) ) + ( -( %f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + ( ( %f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) + ( %f * cos(t*3.14159/180) )" % (A3,t,A3,A3,A15,tout) )
print("y(t):")
print(" %f * sin(t*3.14159/180) + ( %f *sin(t*3.14159/180) ) + ( -( %f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + ( ( %f * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) * ( -tan(t*3.14159/180) ) + ( %f * sin(t*3.14159/180) )" % (A4,-t,A3,A3,A15,-tout) )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (thetaN-0.001) )
print("Rapid Contraction")
print("x(t):")
print(" %f * ( 1 - 0.6 * cos(t*3.14159/180) ) + ( %f * cos(t*3.14159/180) ) + ( -( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + ( %f ) * cos(t*3.14159/180) )^0.5 ) + ( %f * cos(t*3.14159/180) )" % (A5,t,A5,A5,A15,tout) )
print("y(t):")
print(" %f * sin(t*3.14159/180) + ( %f *sin(t*3.14159/180) ) + ( -( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( %f * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) * ( tan(t*3.14159/180) ) + ( %f *sin(t*3.14159/180) )" % (A6,t,A5,A5,A15,tout) )
print("Parameters:" )
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.001) )
A21 = math.cos(math.radians(theta))*tout
A22 = math.sin(math.radians(theta))*tout
print("Linear contraction")
print("x(t):")
print(" t + %f + ( -t + ( t^2 + %f )^0.5 ) + %f" % (A13,A16,A21))
print("y(t):")
print(" %f *t + %f + %f + ( -t + ( t^2 + %f )^0.5 ) * ( %f ) + ( %f ) + ( %f )" % (A7,A17,A18,A19,A20,A14,A22))
print("Parameters:" )
print("t1:")
print(" %f" % (p1+0.0001) )
print("t2:")
print(" %f\n" % (p2-0.0001) )
print("Inlet")
print("x(t):")
print(" %f - %f * (1 - cos(t*3.14159/180) ) + ( %f * cos(t*3.14159/180) ) + ( -( %f - %f * ( 1 - cos(t*3.14159/180) ) ) + ( ( %f - %f * ( 1 - cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) + ( %f * cos(t*3.14159/180) )" % (rcc,A9,t,rcc,A9,rcc,A9,A15,tout))
print("y(t):")
print(" %f * sin(t*3.14159/180) + %f + ( %f * sin(t*3.14159/180) ) + ( -( %f - %f * ( 1 - cos(t*3.14159/180) ) ) + ( ( %f - %f * ( 1 - cos(t*3.14159/180) ) )^2 + %f * cos(t*3.14159/180) )^0.5 ) * tan(t*3.14159/180) + ( %f * sin(t*3.14159/180) )" % (A10,A11,t,rcc,A10,rcc,A10,A15,tout))
print("Parameters:")
print("t1:")
print(" 0.001" )
print("t2:")
print(" %f\n" % (theta-0.0001) )
Nozzle surface, outer: Bell x(t): t + ( 0.018199 * t / ( t^2 + 0.245218 )^0.5 ) + ( t / ( t^2 + 0.245218 )^0.5 ) * ( -t + ( t^2 + 0.370950 * sin(2*3.14159 - arctan(-2.019409*t)) )^0.5 ) / ( sin(2*3.14159 - arctan(-2.019409*t)) ) + ( 0.036397 * t / ( t^2 + 0.245218 )^0.5 ) y(t): 1.009704 *t^2 - 0.237475 + ( -0.018199 / ( 2.019409 *( t^2 + 0.245218 )^0.5)) + ( 1 / ( -2.019409 * ( t^2+ 0.245218 )^0.5 ) ) * ( -t + ( t^2 + 0.370950 * sin(2*3.14159 - arctan(-2.019409*t)) )^0.5 ) / ( sin(2*3.14159 - arctan(-2.019409*t)) ) + ( -0.036397 / ( 2.019409 * ( t^2 + 0.245218 )^0.5 ) ) Parameters: t1: 0.620064 t2: 1.036417 Rapid Expansion x(t): 0.770575 * ( 1- (.382/1.382) * cos(t*3.14159/180) ) + ( 0.018199 * cos(t*3.14159/180) ) + ( -( 0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + ( ( 0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) + ( 0.036397 * cos(t*3.14159/180) ) y(t): 0.212995 * sin(t*3.14159/180) + ( -0.018199 *sin(t*3.14159/180) ) + ( -( 0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) ) + ( ( 0.770575 * ( 1 - (0.382)/(1.382) * cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) * ( -tan(t*3.14159/180) ) + ( -0.036397 * sin(t*3.14159/180) ) Parameters: t1: 0.001 t2: 44.999000 Rapid Contraction x(t): 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) + ( 0.018199 * cos(t*3.14159/180) ) + ( -( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + ( 0.370950 ) * cos(t*3.14159/180) )^0.5 ) + ( 0.036397 * cos(t*3.14159/180) ) y(t): -0.836369 * sin(t*3.14159/180) + ( 0.018199 *sin(t*3.14159/180) ) + ( -( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) ) + ( ( 1.393948 * ( 1 - 0.6 * cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) * ( tan(t*3.14159/180) ) + ( 0.036397 *sin(t*3.14159/180) ) Parameters: t1: 0.001 t2: 44.999000 Linear contraction x(t): t + 0.012868 + ( -t + ( t^2 + 0.262301 )^0.5 ) + 0.025737 y(t): -1.000000 *t + 0.802546 + -0.591402 + ( -t + ( t^2 + 0.262301 )^0.5 ) * ( 1.000000 ) + ( 0.012868 ) + ( 0.025737 ) Parameters: t1: 0.802646 t2: 1.561697 Inlet x(t): 1.806764 - 0.836369 * (1 - cos(t*3.14159/180) ) + ( 0.018199 * cos(t*3.14159/180) ) + ( -( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) ) + ( ( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) + ( 0.036397 * cos(t*3.14159/180) ) y(t): 0.836369 * sin(t*3.14159/180) + -1.942055 + ( 0.018199 * sin(t*3.14159/180) ) + ( -( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) ) + ( ( 1.806764 - 0.836369 * ( 1 - cos(t*3.14159/180) ) )^2 + 0.370950 * cos(t*3.14159/180) )^0.5 ) * tan(t*3.14159/180) + ( 0.036397 * sin(t*3.14159/180) ) Parameters: t1: 0.001 t2: 44.999900
The flow properties at any axial location of a rocket nozzle are modeled as an ideal gas. Therefore, if the particular geometry of the nozzle is known, then the properties within a finite element of the flow may be known. By determining the flow properties of the gas, a heat transfer coefficient may also be found, as well as the heat flux into any particular channel. By assuming that the heat flux into the fluid is equal to the heat flux through the nozzle wall, it may be determined whether a particular finite volume of coolant aquires enough heat to achieve a phase change, resulting in nucleate boiling when the coolant pressures are low enough for that phase change to be from a liquid to a gas.
In some cases; where the pressure is high enough, a coolant super critical state may be achieved. In such cases, there may be advantages to allowing a phase change. It is; however, not advantageous to allow nucleate boiling, which would result in an insulative gas layer forming between the wall and the coolant, quickly leading to engine failure.
The scope of this analysis is solely to determine the amount of heat accepted by a finite volume of fluid and produce a reasonable guess as to the phase state of the coolant at the exit of the convergent section of the nozzle. This will allow for a designer to alter their original assumptions about the converging angle of the nozzle, the rapid expansion angle of the nozzle, or the length of the combustion chamber.
In general, a decrease in the surface area of the nozzle will result in less heat being transfered into the coolant. Therefore, heat transfer may be decreased by an increase of either the convergence angle or the rapid expansion angle of the nozzle. Furthermore, heat transfer may be decreased by shortening the combustion chamber. Typically, these changes will result in lower nozzle efficiencies.
Coolant capacity may be increased by either the dilution of the fuel with water (in the case that the fuel is aqueous), or by increasing the mass flow rate of fuel in the nozzle. The latter is achievable by increasing the desired force of a nozzle, or by increasing the amount of film cooled surface area within the nozzle. Coefficients modifiying the adiabadic wall temperature of the flow due to the addition of film cooling should be determined by measurements performed in testing.
import sympy as sympy
from sympy.solvers import solve
from sympy import Symbol
from sympy import log
import numpy as np
import matplotlib.pyplot as plt
#Boiling analysis
#Heat transfer along bell
#find Ax/At
x = [0]
dx = 0.01
thetaNr = math.pi/180*thetaN
thetar = math.pi/180*theta
import numpy as np
ARatioSub = []
ARatioSup = []
rx = []
Ax = []
AR = []
Px = []
Tx = []
Vx = []
ax = []
Mx = []
sigmax = []
hgx = []
Tawx = []
qx = []
res = 1000 # resolution of mach lookup matrix
MachSub = np.linspace(1/res,1,num=res)
MachSup = np.linspace(1+1/res,10,num=res*9)
PRSub = np.linspace(1/res,2,num=res*2)
PRSup = np.linspace(2+2/res,1000,num=res*998)
#build AR matricies
for Machx in MachSub:
ARatioSub.append(1/Machx*(((1+((gam-1)/2)*Machx**2)/(1+((gam-1)/2)))**((gam+1)/(gam-1)))**0.5)
for Machx in MachSup:
ARatioSup.append(1/Machx*(((1+((gam-1)/2)*Machx**2)/(1+((gam-1)/2)))**((gam+1)/(gam-1)))**0.5)
i=0
while x[i] < ((A1*re**2-A2)-(A1*X1**2-A2)):
rx.append((((A1*re**2-A2)-x[i]+A2)/(A1))**.5)
x.append(x[i]+dx)
Ax.append(math.pi*(rx[i])**2)
AR.append(Ax[i]/At)
Volumex = Ax[i]*dx
Volumet = At*dx
j=0
#match mach number to current area ratio
currentAR = 0
while currentAR < Ax[i]/At and j < len(ARatioSup)-1:
j=j+1
currentAR = ARatioSup[j]
y1 = MachSup[j]
y2 = MachSup[j-1]
x1 = ARatioSup[j]
x2 = ARatioSup[j-1]
Mx.append((Ax[i]/At-x1)*(y2-y1)/(x2-x1)+y1)
Px.append(Pcns/(1+0.5*(gam-1)*Mx[i]**2)**(gam/(gam-1)))
Tx.append(Tcns/(Pcns/Px[i])**((gam-1)/gam))
Vx.append(math.sqrt(2*g*gam/(gam-1)*R*Tcns*(1-(Px[i]/Pcns)**((gam-1)/gam))))
ax.append(math.sqrt(g*gam*R*Tx[i]))
sigmax.append((1/((.5*Twg/Tcns*(1+(gam-1)/2*Mx[i]**2)+.5)**0.68*(1+(gam+1)/2*Mx[i]**2)**0.12)))
hgx.append(((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*(At/Ax[i])**0.9*sigmax[i]))
Tawx.append(edafilm*Tcns*((1+rturb*((gam-1)/2)*Mx[i]**2)/(1+((gam-1)/2)*Mx[i]**2)))
qx.append(hgx[i]*(Tawx[i]-Twg))
i=i+1
while x[i] < ((A1*re**2-A2)-(A1*X1**2-A2))+.382*rt*sin(thetaN/180*math.pi):
yr = Symbol("yr")
rx.append(solve((x[i]-(A1*re**2-A2))**2+(yr-1.382*rt)**2-(.382*rt)**2,yr)[0])
x.append(x[i]+dx)
Ax.append(math.pi*(rx[i])**2)
AR.append(Ax[i]/At)
Volumex = Ax[i]*dx
Volumet = (At)*dx
j=0
#match mach number to current area ratio
currentAR = 0
while currentAR < Ax[i]/At and j < len(ARatioSup)-1:
j=j+1
currentAR = ARatioSup[j]
y1 = MachSup[j]
y2 = MachSup[j-1]
x1 = ARatioSup[j]
x2 = ARatioSup[j-1]
Mx.append((Ax[i]/At-x1)*(y2-y1)/(x2-x1)+y1)
Px.append(Pcns/(1+0.5*(gam-1)*Mx[i]**2)**(gam/(gam-1)))
Tx.append(Tcns/(Pcns/Px[i])**((gam-1)/gam))
Vx.append(math.sqrt(2*g*gam/(gam-1)*R*Tcns*(1-(Px[i]/Pcns)**((gam-1)/gam))))
ax.append(math.sqrt(g*gam*R*Tx[i]))
sigmax.append((1/((.5*Twg/Tcns*(1+(gam-1)/2*Mx[i]**2)+.5)**0.68*(1+(gam+1)/2*Mx[i]**2)**0.12)))
hgx.append(((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*(At/Ax[i])**0.9*sigmax[i]))
Tawx.append(edafilm*Tcns*((1+rturb*((gam-1)/2)*Mx[i]**2)/(1+((gam-1)/2)*Mx[i]**2)))
qx.append(hgx[i]*(Tawx[i]-Twg))
i=i+1
while x[i] < ((A1*re**2-A2)-(A1*X1**2-A2))+.382*rt*sin(thetaN/180*math.pi)+1.5*rt*sin(theta/180*math.pi):
yr = Symbol("yr")
rx.append(solve((x[i]-(A1*re**2-A2))**2+(yr-2.5*rt)**2-(1.5*rt)**2,yr)[0])
x.append(x[i]+dx)
Ax.append(math.pi*(rx[i])**2)
AR.append(Ax[i]/At)
Volumex = Ax[i]*dx
Volumet = (At)*dx
j=0
#match mach number to current area ratio
currentAR = ARatioSub[0]
while currentAR > Ax[i]/At and j < len(ARatioSub)-1:
j=j+1
currentAR = ARatioSub[j]
y1 = MachSub[j]
y2 = MachSub[j-1]
x1 = ARatioSub[j]
x2 = ARatioSub[j-1]
Mx.append((Ax[i]/At-x1)*(y2-y1)/(x2-x1)+y1)
Px.append(Pcns/(1+0.5*(gam-1)*Mx[i]**2)**(gam/(gam-1)))
Tx.append(Tcns/(Pcns/Px[i])**((gam-1)/gam))
Vx.append(math.sqrt(2*g*gam/(gam-1)*R*Tcns*(1-(Px[i]/Pcns)**((gam-1)/gam))))
ax.append(math.sqrt(g*gam*R*Tx[i]))
sigmax.append((1/((.5*Twg/Tcns*(1+(gam-1)/2*Mx[i]**2)+.5)**0.68*(1+(gam+1)/2*Mx[i]**2)**0.12)))
hgx.append(((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*(At/Ax[i])**0.9*sigmax[i]))
Tawx.append(edafilm*Tcns*((1+rturb*((gam-1)/2)*Mx[i]**2)/(1+((gam-1)/2)*Mx[i]**2)))
qx.append(hgx[i]*(Tawx[i]-Twg))
i=i+1
m = math.tan(thetar)
b = rt*(2.5-1.5*math.cos(thetar))-m*((A1*re**2-A2)+1.5*rt*math.sin(thetar))
while x[i] < (A1*re**2-A2)-(A7*p2+A8):
rx.append(m*x[i]+b)
x.append(x[i]+dx)
Ax.append(math.pi*(rx[i])**2)
AR.append(Ax[i]/At)
Volumex = Ax[i]*dx
Volumet = (At)*dx
j=0
#match mach number to current area ratio
currentAR = ARatioSub[0]
while currentAR > Ax[i]/At and j < len(ARatioSub)-1:
j=j+1
currentAR = ARatioSub[j]
y1 = MachSub[j]
y2 = MachSub[j-1]
x1 = ARatioSub[j]
x2 = ARatioSub[j-1]
Mx.append((Ax[i]/At-x1)*(y2-y1)/(x2-x1)+y1)
Px.append(Pcns/(1+0.5*(gam-1)*Mx[i]**2)**(gam/(gam-1)))
Tx.append(Tcns/(Pcns/Px[i])**((gam-1)/gam))
Vx.append(math.sqrt(2*g*gam/(gam-1)*R*Tcns*(1-(Px[i]/Pcns)**((gam-1)/gam))))
ax.append(math.sqrt(g*gam*R*Tx[i]))
sigmax.append((1/((.5*Twg/Tcns*(1+(gam-1)/2*Mx[i]**2)+.5)**0.68*(1+(gam+1)/2*Mx[i]**2)**0.12)))
hgx.append(((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*(At/Ax[i])**0.9*sigmax[i]))
Tawx.append(edafilm*Tcns*((1+rturb*((gam-1)/2)*Mx[i]**2)/(1+((gam-1)/2)*Mx[i]**2)))
qx.append(hgx[i]*(Tawx[i]-Twg))
i=i+1
xcenter = (A1*re**2-A2)-(A10*math.sin(0)+A11)
ycenter = rcc-1.5*rt
while x[i] < (A1*re**2-A2)-(A10*math.sin(0)+A11):
yr = Symbol("yr")
rx.append(solve((x[i]-xcenter)**2+(yr-ycenter)**2-(1.5*rt)**2,yr)[1])
x.append(x[i]+dx)
Ax.append(math.pi*(rx[i])**2)
AR.append(Ax[i]/At)
Volumex = Ax[i]*dx
Volumet = (At)*dx
j=0
#match mach number to current area ratio
currentAR = ARatioSub[0]
while currentAR > Ax[i]/At and j < len(ARatioSub)-1:
j=j+1
currentAR = ARatioSub[j]
y1 = MachSub[j]
y2 = MachSub[j-1]
x1 = ARatioSub[j]
x2 = ARatioSub[j-1]
Mx.append((Ax[i]/At-x1)*(y2-y1)/(x2-x1)+y1)
Px.append(Pcns/(1+0.5*(gam-1)*Mx[i]**2)**(gam/(gam-1)))
Tx.append(Tcns/(Pcns/Px[i])**((gam-1)/gam))
Vx.append(math.sqrt(2*g*gam/(gam-1)*R*Tcns*(1-(Px[i]/Pcns)**((gam-1)/gam))))
ax.append(math.sqrt(g*gam*R*Tx[i]))
sigmax.append((1/((.5*Twg/Tcns*(1+(gam-1)/2*Mx[i]**2)+.5)**0.68*(1+(gam+1)/2*Mx[i]**2)**0.12)))
hgx.append(((0.026/(2*rt)**0.2*(mucc**0.2*Cpg/Pr**0.6)*(Pcns*g/cstar)**0.8*(2*rt/rmean)**0.1)*(At/Ax[i])**0.9*sigmax[i]))
Tawx.append(edafilm*Tcns*((1+rturb*((gam-1)/2)*Mx[i]**2)/(1+((gam-1)/2)*Mx[i]**2)))
qx.append(hgx[i]*(Tawx[i]-Twg))
i=i+1
del x[-1]
plt.figure(0)
plt.scatter(x,rx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Radius (in)')
plt.xlim([0,math.ceil(max(x))])
plt.figure(1)
plt.scatter(x,Mx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Mach Number')
plt.xlim([0,math.ceil(max(x))])
plt.figure(2)
plt.scatter(x,Px)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Pressure (psia)')
plt.xlim([0,math.ceil(max(x))])
plt.figure(3)
plt.scatter(x,Tx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Temperature (R)')
plt.xlim([0,math.ceil(max(x))])
plt.figure(4)
plt.scatter(x,Vx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Velocity (ft/s)')
plt.xlim([0,math.ceil(max(x))])
plt.figure(5)
plt.scatter(x,hgx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Heat transfer coefficient (Btu/in^2/s/R)')
plt.xlim([0,math.ceil(max(x))])
plt.figure(6)
plt.scatter(x,Tawx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Adiabatic wall Temp of gas (R)')
plt.xlim([0,math.ceil(max(x))])
plt.figure(7)
plt.scatter(x,qx)
plt.xlabel('Nozzle Position (in)')
plt.ylabel('Heat flux (Btu/in^2/s)')
plt.xlim([0,math.ceil(max(x))])
(0, 3)
#Determine heat added to fuel
j=1
Aqx=[]
Qx =[]
lenCC = 0
while j < len(rx):
Aqx.append(math.sqrt(dx**2+(rx[j-1]-rx[j])**2)*math.pi*(rx[j-1]+rx[j]))
Qx.append(Aqx[j-1]*(qx[j-1]+qx[j])*.5)
lenCC = lenCC + math.sqrt(dx**2+(rx[j-1]-rx[j])**2)
j+=1
Qtot = sum(Qx)
print("Maximum combustion chamber length: %f in" % lenCC)
print("Current combustion chamber length: %f in" % Lc)
# print((Qc-Qtot)/(qx[len(qx)-1]*2*math.pi*(rx[len(rx)-1]))) #determine exactly what this is, looks like a ratio of heat fluxes
print("Heat transfered into coolant: %f" % Qtot)
print("Coolant Capacity: %f" % Qc)
Maximum combustion chamber length: 3.358804 in Current combustion chamber length: 2.730269 in Heat transfered into coolant: 218.211184 Coolant Capacity: 341.681397
#pressure loss in cooling passage
delP = f*rhof/2*Vco**2/d*(lenCC+Lc)
#print(f)
#print(rhof)
print("Coolant velocity: %f ft/sec" % Vco )
#print(d)
print("Cooling channel pressure drop: %f psi" % delP )
Coolant velocity: 58.380550 ft/sec Cooling channel pressure drop: 142.058534 psi
#consider saving this to a file rather than printing out.
#Properties of flow along the axial length of the nozzle
#i=0
#print("x,rx,Mx,Px,Tx,Vx,hgx,Tawx")
#while i < len(x):
# print("%f,%f,%f,%f,%f,%f,%f,%f" % (x[i],rx[i],Mx[i],Px[i],Tx[i],Vx[i],hgx[i],Tawx[i]) )
# i=i+1