cuspaceflight / bamboo

Cooling system modelling for liquid rocket engines.
GNU Affero General Public License v3.0
16 stars 4 forks source link

Gradient discontinuity at inflection point when using Nozzle.y(x) #22

Closed dug20 closed 3 years ago

dug20 commented 3 years ago

Currently being looked into - I believe this is the reason for heat transfer peaks not being at the throat.

dug20 commented 3 years ago

It seems that (x, y) at the inflection point and (x, y) at the exit cannot both be satisfied by the parabola fitting used. There is likely a typo somewhere causing this.

dug20 commented 3 years ago

It is not possible to perfectly fit a parabola to the Rao bell nozzle data, using the current method in place. There are four constraints on the parabola (inflection point coordinates, inflection point angle, exit point coordinates, exit point angle), but a quadratic only has 3 degrees of freedom (x = ay^2 + by + c).

Previous problems arose because we were ignoring the nozzle exit coordinates when doing parabola fitting, and only used the other three constraints. This meant that the exit area given by the parabola disagreed with the exit area that the user had requested for the nozzle.

A temporary fix is in place, whereby the quadratic is only fitted to the exit coordinates, exit angle, and inflection point coordinates. This means that the inflection point angle is not used, so there will be a discontinuity in gradient at the inflection point.

One possible long term fix is to use the following three constraints for quadratic fitting:

This would completely eliminate the need for the angle at the inflection point.

The problem can be demonstrated below, where we print the inflection angle obtained from reading Rao's bell nozzle graphs, and compare it to the shape that bamboo is using internally (for the nozzle.y(x) function).

'''
Reference [1] for Rao bell nozzle graphs: http://www.aspirespace.org.uk/downloads/Thrust%20optimised%20parabolic%20nozzle.pdf
'''

import bamboo as bam
import numpy as np
import matplotlib.pyplot as plt

'''Gas properties - obtained from ProPEP 3'''
gamma = 1.264               #Ratio of specific heats cp/cv
molecular_weight = 21.627   #Molecular weight of the exhaust gas (kg/kmol) (only used to calculate R, and hence cp)

'''Chamber conditions'''
pc = 10e5           #Chamber pressure (Pa)
Tc = 2458.89        #Chamber temperature (K) - obtained from ProPEP 3
mdot = 4.757        #Mass flow rate (kg/s)
p_amb = 1.01325e5   #Ambient pressure (Pa). 1.01325e5 is sea level atmospheric.

'''Create the engine object'''
perfect_gas = bam.PerfectGas(gamma = gamma, molecular_weight = molecular_weight)
chamber = bam.ChamberConditions(pc, Tc, mdot)
nozzle = bam.Nozzle.from_engine_components(perfect_gas, chamber, p_amb, type = "rao", length_fraction = 0.8)
engine = bam.Engine(perfect_gas, chamber, nozzle)

print(nozzle)

#Notice that the gradient of the fitted parabola (which is used for nozzle.y()) does not exactly match the Rao theta_n from the graphs in Reference [1].
print(f"Rao theta_n from reading off graphs = {nozzle.theta_n}")
dydx = (nozzle.y(nozzle.Nx) - nozzle.y(nozzle.Nx - 0.001))/0.001
print(f"Rao theta_n from nozzle.y(x) = {np.arctan(dydx)}")

nozzle.plot_nozzle()
plt.show()
dug20 commented 3 years ago

Fixed by fitting parabola to (x, y) at inflection point, and gradient at inflection point. Correct exit area is then achieved by extrapolating parabola (or cone) until we reach the required exit area.