I haven't posted recently, as you'll no doubt have noticed. A great deal has happened to me since the last time I made an entry here, almost none of it related to rocketry. Those who need to know, know. Suffice it to say, I'm in a much better frame of mind now than I was and I am ready to take up my project where I left off, and begin sharing ideas with all of you.
What I want to talk about today is cooling. Apparently, the engineer who was in charge of the tube bundle development for the Rocketdyne F1 was told to "make sure it doesn't melt". The goal of rocket motor cooling can't really be stated any clearer than that, so here are my thoughts on the subject.
All of my earlier attempts to estimate cooling requirements were based on the approach given in Krzycki, and I tried to compare my engine design to similar sized ones in the literature to get some idea of the heat flux that might be encountered. I decided to adopt a slightly more stringent approach. To that end, I studied various sources on the subject, the chief ones being Huzel & Huang and Humble, Henry and Larson. I also dipped into a lot of the NTRS papers, as well as The RAE's Beta Project report.
The best amateur analysis of the regenerative cooling problem I have found online is that by G N Sortino at Fubarlabs, "Thermal Analysis of Steady State Engines" which can be found here:-
wiki.fubarlabs.org/fubarwiki/Thermal-Analysis-of-Steady-State-Engines.ashxI quickly found, as did the above author, that there are a lot of tricky variables in these equations, and it is not at all easy to know where to begin. Luckily, Huzel and Huang's approach to cooling design is eminently practical. By combining their approach with an engine concept that I felt could be fabricated, I started to make some small headway.
I began by thinking about an engine composed of a tube bundle and with a thrust of 1000lbf. I arbitrarily decided on a tube ID of 6mm, purely on a rule of thumb basis. I then ran a simulation on RPA of this hypothetical engine, the propellants being a 90% ethanol/water mix with Nitrous Oxide. Here are the leading particulars for this thrust chamber as calculated by RPA:-
#---------------------------------------------------------------------------------------------------
#
# Table 1. Thermodynamic properties
#---------------------------------------------------------------------------------------------------
# Parameter Injector Nozzle inl Nozzle thr Nozzle exi Unit
#---------------------------------------------------------------------------------------------------
Pressure 2.0684 2.0661 1.1392 0.1015 MPa
Temperature 1867.3657 1867.1488 1646.1334 995.3005 K
Enthalpy -1347.9581 -1348.3832 -1778.8090 -3060.6928 kJ/kg
Entropy 11.1984 11.1986 11.1986 11.1986 kJ/(kg·K)
Internal energy -2117.3892 -2117.7249 -2457.0646 -3469.2680 kJ/kg
Specific heat (p=const) 1.9622 1.9622 1.9364 2.5923 kJ/(kg·K)
Specific heat (V=const) 1.5496 1.5496 1.5241 2.0981 kJ/(kg·K)
Gamma 1.2663 1.2663 1.2705 1.2356
Isentropic exponent 1.2662 1.2662 1.2705 1.2270
Gas constant 0.4120 0.4120 0.4120 0.4105 kJ/(kg·K)
Molecular weight (gas) 20.1788 20.1788 20.1793 20.2543
Molar mass (gas) 0.0202 0.0202 0.0202 0.0203 kg/mol
Molar mass (total) 0.0202 0.0202 0.0202 0.0203 kg/mol
Density 2.6883 2.6856 1.6796 0.2484 kg/m³
Sonic velocity 987.0341 986.9786 928.2780 708.0312 m/s
Velocity 0.0000 29.1568 928.2780 1850.8024 m/s
Mach number 0.0000 0.0295 1.0000 2.6140
Area ratio 19.9000 19.9000 1.0000 3.3900
Mass flux 78.3036 78.3036 1559.0966 459.7134 kg/(m²·s)
Viscosity 0.0001 0.0001 0.0001 0.0000 kg/(m·s)
Conductivity, frozen 0.2122 0.2122 0.1908 0.1306 W/(m·K)
Specific heat (p=const), frozen 1.9240 1.9240 1.8840 1.7340 kJ/(kg·K)
Prandtl number, frozen 0.5869 0.5869 0.5868 0.5600
Conductivity, effective 0.2177 0.2168 0.1919 nan W/(m·K)
Specific heat (p=const), effective 1.9550 1.9610 1.8920 nan kJ/(kg·K)
Prandtl number, effective 0.5814 0.5855 0.5860 nan
#---------------------------------------------------------------------------------------------------
The first step in the heat transfer calculation was to find the gas side heat transfer coefficient, hg. Before doing this, figures for adiabatic wall temperature and the actual gas side cooled wall temperature were derived.
The adiabatic wall temperature may be found at any point in the chamber by multiplying the stagnation temperature by a stagnation recovery factor, typically 0.9. I decided to model the cooling on the throat case as this would be the most stringent. At the throat, the stagnation recovery factor is 1, so it can be seen from this and the figures above that the adiabatic wall temperature at the throat is about 1867.15 K.
Taw = 1867.15 K
The next figure to determine was the gas side cooled wall temperature, or Twg. Of course there is no way to know what this is, so the approach is to state a figure that it is within the wall material's ability to withstand. I envisaged a tube wall bundle made of 316L tubes welded together. I wanted to keep this below 400 deg C so I decided on a figure for Twg of 300 deg C or 573.15 K.
Twg = 573.15 K
Next, to calculate hg:-
hg = 0.026k (pv/u)^0.8 (1/Dt)^0.2 (Cpu/k)^0.4
Where:-
k = gas thermal conductivity
p = gas density
v = gas velocity
u = gas dynamic viscosity
Dt = throat diameter
Cp = gas specific heat at constant pressure
All of the above values can be got from the RPA table, apart from Dt. For the engine under consideration, standard design equations put the throat diameter at 45mm i.e. 45 x 10^-3m.
With these values in the hg equation, the following values resulted:-
hg = 0.00494 x 861439.4 x 1.859 x 0.0516
hg = 408.2 W/m^2K
Using hg, Taw and Twg the calculation can now be moved on further to gain the throat heat flux, q:-
q = (Taw - Twg) x hg
q = (1867.15 - 573.15) x 408.2
q = 528.2 x 10^3 W/m^2
It is now possible to calculate the temperature of the cooled side of the wall Twc, as follows:-
Twc = Twg - qt/k
Where:-
t = wall thickness
k = wall material thermal conductivity
I had decided to use metric hydraulic tubing for the bundle material. Tube of 10mm OD has a 6mm ID hence the wall thickness t is 2mm. The thermal conductivity of 316L at 300 deg C is about 16 W/m^2.
Hence:-
Twc = 573.15 - [(528.2 x 10^3 x 2 x 10^-3)/16]
Twc = 573.15 - 66
Twc = 507.15 K
This figure of 507.15 K for Twc is rather apposite, as it corresponds to 234 deg C. We know that the oxidiser, nitrous oxide, sits at about 760 psi at STP. We will want to pressurise it above and beyond this value for reasons of both safety and flow homogeneity. If we assume a figure of 800 psi as our pressurant level, then we can use this for the fuel as well. The boiling point of pure ethanol at 800 psi is about 240 deg C. The wall temperature on the cooled side is below this, and this with ethanol/water. By adjusting the coolant pressure we can tune the response to induce some nucleate boiling.
It is now possible to calculate the required coolant heat transfer coefficient to permit the calculated gas side heat transfer and temperature differential. Huzel and Huang give a practical method of doing this:-
hc = q/(Tcw - Tco)
where Tco = coolant temperature. What figure must be used for this? In Huzel and Huang it is assumed that the coolant flowing through the throat tube is in a two pass circuit and has already been through the throat once. A value of 600 deg R is assigned for Tco. This is 60 deg C. I was troubled by this figure as I felt it would be neccesary to account for the temperature increase the coolant would see as it traversed the chamber, let alone the throat. As the flow of the coolant through the tubes is turbulent and is constantly being renewed, it would seem that the instantaneous coolant temperature would be unlikely to approach the values of Tcw in the various parts of the chamber. Therefore I decided to use a "film temperature" for the Tco figure. Film temperature of the coolant can be defined as:-
Tfilm = 0.5(Tcw - Tcs)
Where Tcs = coolant static temperature, i.e. the temperature on entering the coolant circuit.
So for the above case, Tfilm = 0.5(507.15 K - 293.15 K)
Tfilm = Tco = 380.15 K i.e. ~ 100 deg C
So the hc required to achieve steady state in the cooling circuit becomes:-
hc = q/(Tcw - Tco)
hc = 528.2 x 10^3/(507.15 - 380.15)
hc = 4159 W/m^2K
Recall now that I was talking about using 316L tubes in this design. The tubes have a 10mm OD. For a 1000lbf engine, the fuel flowrate will be about 0.84 Kg/sec. The chamber diameter will be about 200mm. This means a total of 63 tubes to form the chamber. Now, the formula to calculate hc for a given flowrate and tube diameter is:-
hc = 0.023 (Cp m/A) (u/Dpv)^0.2 (k/uCp)^0.66
Where:-
Cp = coolant specific heat
m = coolant flowrate
A = coolant flow area
D = coolant flow diameter
u = coolant dynamic viscosity
k = coolant thermal conductivity
p = coolant density
The coolant transport properties (density, specific heat and so on) were all calculated at the Tco value of 380.15 K. They are detailed below:-
Density = 755.6Kg/m^3
Thermal conductivity = 0.185 w/mK
Specific heat = 3.39 x 10^3 J/KgK
Dynamic viscosity = 0.3 x 10^-3 N-s/m^2
These values were calculated for 90% ethanol/water at Tco of 380.15 K. The thermal conductivity was calculated using the Fillipov equation. Dynamic viscosity was found by using the Refutas equation to find kinematic viscosity and then multiplying this by the density.
Going back to the practical aspect of 316L tubes, I envisaged the coolant circuit as consisting of a series of tubes joined by manifolds to allow the coolant to pass up and down the chamber, thereby splitting the 63 tubes up into groups of coolant passages. It can be seen that the values to be used in the hc calculation are essentially constants. Only changing the number of tubes in a coolant passage group will alter the coolant flowrate and thus velocity. So the problem then became one of iteration to come as close to an hc value of 4159 W/m^2K as was practicable.
It was found that by splitting the 63 tubes up into 21 coolant passage groups of 3 tubes each, the following result was gained:-
21 groups of 3 tubes, each group of 3 tubes forming one coolant passage
0.84 kg/sec / 21 = 0.04 kg/sec
0.04 kg/sec through a 6mm ID tube gives a velocity of 1.87 m/sec
then:-
hc = 0.023 ((3.39x10^3 x 0.04)/2.83x10^-5) (0.3x10^-3/(6x10^-3 x 755.6 x 1.87))^0.2 (0.185/0.3x10^-3 x 3.39x10^3)^0.66
hc = 0.023 x 4791519.4 x 0.129 x 0.325
hc = 4620.3 W/m^2K
Recall that the figure required was 4159 W/m^2K. Do these calculations mean that this hypothetical system would be steady state? I'd say yes, tentatively. They are as far as I've got up to now and they feel right. I'm open to all comments and that is why I have posted them here. I'd be very grateful for any thoughts on all of this.
Carl.
