Mathematical calculation of engine fuel injectors

There is a code that performs calculations of parts of fuel injectors. There are two classes in which the basic calculations take place, they are called LiquidJetInjector and GasJetInjector, respectively.

The Reynolds class stores constants for calculations.

The SimilarCalculations class is a class in which the calculations required in both LiquidJetInjector and GasJetInjector take place.

In order for the code to work correctly, first we need to pass the basic values to the SimilarCalculations class, and then to the class whose calculations we are interested in.

import math
from enum import Enum

class SimilarCalculations:
def __init__(self, diameter: float, length: float, mass_flow_rate: float, viscosity: float, density: float):
self.diameter = diameter
self.length = length
self.mass_flow_rate = mass_flow_rate
self.viscosity = viscosity
self.density = density

def injector_nozzle_area(self) -> float:
return (math.pi * self.diameter ** 2) / 4

def reynolds_number(self) -> float:
return (4 * self.mass_flow_rate) / (math.pi * self.diameter * self.viscosity)

def average_speed(self) -> float:
return self.mass_flow_rate / (self.density * self.injector_nozzle_area())

def relative_length_injector(self) -> float:
return self.length / self.diameter

class Reynolds(Enum):
LAMINAR = 2000
TURBULENT = 10000

class LiquidJetInjector:
def __init__(self, density_comb: float, sigma_fuel: float):
self.calculations = SimilarCalculations(diameter, length, mass_flow_rate, viscosity, density)

self.density_comb = density_comb
self.sigma_fuel = sigma_fuel

self.reynolds_number = self.calculations.reynolds_number()
self.injector_nozzle_area = self.calculations.injector_nozzle_area()
self.average_speed = self.calculations.average_speed()

def linear_hydraulic_resistance(self) -> float:
if self.reynolds_number < Reynolds.LAMINAR.value:
return 64 / self.reynolds_number

elif Reynolds.LAMINAR.value <= self.reynolds_number <= Reynolds.TURBULENT.value:
return 0.3164 * self.reynolds_number ** (-0.25)

else:
return 0.031

def injector_losses_inlet(self) -> float:
if self.reynolds_number < Reynolds.LAMINAR.value:
return 2.2 - 0.726 * math.exp(-74.5 * ((viscosity * length) / mass_flow_rate))
else:
return 1 + 2.65 * self.linear_hydraulic_resistance()

def injector_flow_coefficient(self) -> float:
return 1 / (math.sqrt(self.injector_losses_inlet() + self.linear_hydraulic_resistance() *
(length / diameter)))

def pressure_drop_injector(self) -> float:
return mass_flow_rate ** 2 / (2 * density * self.injector_flow_coefficient() ** 2 * self.injector_nozzle_area **
2)

def weber_criterion(self) -> float:
return (self.density_comb * self.average_speed ** 2 * diameter) / self.sigma_fuel

def media_diameter_spray_droplets(self) -> float:
return diameter * round(math.pow((27 * math.pi) / 4, 1 / 3.)) * self.weber_criterion() ** (-0.333)

class GasJetInjector:
def __init__(self, combustion_pressure: float, pressure_drop_internal_circuit: float, gas_constant_gen_gas: float,
temperature_gen_gas: float, entropy_expansion_ratio: float):
self.calculations = SimilarCalculations(diameter, length, mass_flow_rate, viscosity, density)

self.reynolds_number = self.calculations.reynolds_number()
self.injector_nozzle_area = self.calculations.injector_nozzle_area()
self.average_speed = self.calculations.average_speed()

self.combustion_pressure = combustion_pressure
self.pressure_drop_internal_circuit = pressure_drop_internal_circuit
self.gas_constant_gen_gas = gas_constant_gen_gas
self.temperature_gen_gas = temperature_gen_gas
self.entropy_expansion_ratio = entropy_expansion_ratio

def calculate_injector_pressure(self) -> float:
return self.combustion_pressure + self.pressure_drop_internal_circuit

def density_gen_gas(self) -> float:
return self.calculate_injector_pressure() / (self.gas_constant_gen_gas * self.temperature_gen_gas)

def injector_flow_coefficient(self) -> float:
return (math.sqrt(1.23 ** 2 + (232 * length) / (self.reynolds_number * diameter))) / \
((116 * length) / (self.reynolds_number * diameter))

def injector_nozzle_area_outlet(self) -> float:
return mass_flow_rate / (self.injector_flow_coefficient() * density * (
self.combustion_pressure / self.calculate_injector_pressure()) ** (1 / self.entropy_expansion_ratio) *
math.sqrt(2 * (self.entropy_expansion_ratio / (self.entropy_expansion_ratio - 1)) *
self.gas_constant_gen_gas * self.temperature_gen_gas * (1 - (
self.combustion_pressure / self.calculate_injector_pressure()) ** (
(self.entropy_expansion_ratio - 1) / self.entropy_expansion_ratio))))

def checking_diameter_injector(self) -> float:
return math.sqrt((4 * self.injector_nozzle_area_outlet()) / math.pi)


To test the values at the end of the code, I used this construction (it is temporary and will not be in the final version, but I will attach it to understand the work)

#
diameter = 0.0005
length = 0.002
mass_flow_rate = 0.00712
viscosity = 0.00116
density = 1315
#
density_comb = 0.1
sigma_fuel = 0.5
#
combustion_pressure = 14.58 * 10**6
pressure_drop_internal_circuit = 10**6
gas_constant_gen_gas = 260.5
temperature_gen_gas = 777.4
entropy_expansion_ratio = 1.23
#
calc = SimilarCalculations(diameter, length, mass_flow_rate, viscosity, density)
calc_2 = LiquidJetInjector(density_comb, sigma_fuel)
calc_3 = GasJetInjector(combustion_pressure, pressure_drop_internal_circuit, gas_constant_gen_gas, temperature_gen_gas, entropy_expansion_ratio)
#
print(calc_3.injector_nozzle_area_outlet())
print(calc_2.injector_flow_coefficient())
print(calc.injector_nozzle_area())


For a more accurate understanding of the formulas and parameters written in the code, I attach a link to the source of the manual according to which the code was written. The required pages are 38-44.

The manual for the calculation of liquid propellant injectors

• Have you thought about having the SimilarCalculations class as a base class and inheriting it in the specific injector classes? This seems like a good case for inheritance.
– TomG
Commented Feb 15 at 19:02

SimilarCalculations is not a very good name. If we interpret it to be a superclass of your other two classes, it can simply be Injector. A superclass would be a simple and remove the need for a calculations member.

In your equations there's a lot of redundant parens, and not enough linebreaks. I'll demonstrate an alternative format that I find a little easier to understand.

Nearly all of your methods are simple and non-mutating. That's good! It means that they can be made @property.

Because your code is already non-mutating, you can tighten the class definition to a frozen dataclass.

You use math.pow in a place where the ** operator is simpler. Also, you write ** (-0.333) in the same line as raising to the power of 1 / 3.; the latter will be more accurate.

In your demo, you describe it as temporary - but it should be permanent! Convert it to unit tests.

Suggested

import math
from dataclasses import dataclass
from enum import Enum

@dataclass(frozen=True)
class Injector:
density: float
diameter: float
length: float
mass_flow_rate: float
viscosity: float

@property
def injector_nozzle_area(self) -> float:
return math.pi * self.diameter**2 / 4

@property
def reynolds_number(self) -> float:
return 4 * self.mass_flow_rate / (math.pi * self.diameter * self.viscosity)

@property
def average_speed(self) -> float:
return self.mass_flow_rate / (self.density * self.injector_nozzle_area)

@property
def relative_length_injector(self) -> float:
return self.length / self.diameter

class Reynolds(Enum):
LAMINAR = 2000
TURBULENT = 10000

@dataclass(frozen=True, slots=True)
class LiquidJetInjector(Injector):
density_comb: float
sigma_fuel: float

@property
def linear_hydraulic_resistance(self) -> float:
if self.reynolds_number < Reynolds.LAMINAR.value:
return 64 / self.reynolds_number
elif Reynolds.LAMINAR.value <= self.reynolds_number <= Reynolds.TURBULENT.value:
return 0.3164 * self.reynolds_number**-0.25
return 0.031

@property
def injector_losses_inlet(self) -> float:
if self.reynolds_number < Reynolds.LAMINAR.value:
return 2.2 - 0.726*math.exp(-74.5 * self.viscosity * self.length / self.mass_flow_rate)
return 1 + 2.65 * self.linear_hydraulic_resistance

@property
def injector_flow_coefficient(self) -> float:
return 1 / math.sqrt(
self.injector_losses_inlet
+ self.linear_hydraulic_resistance * self.length / self.diameter
)

@property
def pressure_drop_injector(self) -> float:
return self.mass_flow_rate**2 / (
2 * self.density * self.injector_flow_coefficient**2 *
self.injector_nozzle_area**2
)

@property
def weber_criterion(self) -> float:
return self.density_comb * self.average_speed**2 * self.diameter / self.sigma_fuel

@property
def media_diameter_spray_droplets(self) -> float:
return self.diameter * round(
(27 * math.pi / 4)**(1/3)
) * self.weber_criterion**(-1/3)

@dataclass(frozen=True, slots=True)
class GasJetInjector(Injector):
combustion_pressure: float
pressure_drop_internal_circuit: float
gas_constant_gen_gas: float
temperature_gen_gas: float
entropy_expansion_ratio: float

@property
def injector_pressure(self) -> float:
return self.combustion_pressure + self.pressure_drop_internal_circuit

@property
def density_gen_gas(self) -> float:
return self.injector_pressure / (self.gas_constant_gen_gas * self.temperature_gen_gas)

@property
def injector_flow_coefficient(self) -> float:
return math.sqrt(
1.23**2 + 232*self.length / (self.reynolds_number * self.diameter)
) / (116 * self.length) * self.reynolds_number * self.diameter

@property
def injector_nozzle_area_outlet(self) -> float:
return self.mass_flow_rate / (
self.injector_flow_coefficient * self.density * (
self.combustion_pressure / self.injector_pressure
) ** (1 / self.entropy_expansion_ratio) *
math.sqrt(
2 * self.entropy_expansion_ratio / (self.entropy_expansion_ratio - 1) *
self.gas_constant_gen_gas * self.temperature_gen_gas * (
1 - (
self.combustion_pressure / self.injector_pressure
) ** (
(self.entropy_expansion_ratio - 1) / self.entropy_expansion_ratio
)
)
)
)

@property
def checking_diameter_injector(self) -> float:
return math.sqrt(4 * self.injector_nozzle_area_outlet / math.pi)

def test() -> None:
common = {
'diameter': 0.0005,
'length': 0.002,
'mass_flow_rate': 0.00712,
'viscosity': 0.00116,
'density': 1315,
}

liquid = LiquidJetInjector(
**common, density_comb=0.1, sigma_fuel=0.5,
)
gas = GasJetInjector(
**common,
combustion_pressure=14.58e6,
pressure_drop_internal_circuit=1e6,
gas_constant_gen_gas=260.5,
temperature_gen_gas=777.4,
entropy_expansion_ratio=1.23,
)

assert math.isclose(gas.injector_nozzle_area_outlet, 8.279311695661635e-10)
assert math.isclose(gas.injector_nozzle_area, 1.9634954084936206e-07)
assert math.isclose(liquid.injector_flow_coefficient, 0.9105406506118756)
assert math.isclose(liquid.injector_nozzle_area, 1.9634954084936206e-07)

if __name__ == '__main__':
test()


naming

Carefully choose identifiers. These are terrific names:

• GasJetInjector
• LiquidJetInjector

So is Reynolds. But a "similar" prefix is not a good choice. It describes a relationship between others, rather than what it is. Minimally prefer a "common" prefix. And making each class name a noun is good, but this is not exactly a "calculations" object.

From just reading your names, a JetInjector class seems most appropriate. Or perhaps Injector.

two things in one

The __init__ ctor accepts five parameters. Given that there's no """docstring""" I assume we have MKS SI units. The first two are perfect: they describe an Injector. The next three describe dynamic fluid conditions, and maybe belong in a Fluid or Flow @dataclass. Or perhaps they are the nominal operating point this injector was designed for? Except that you didn't # comment on that, nor cite a reference where I could tell what is usual in this situation. They're certainly not "do not exceed" parameters, since they don't start with a max_ prefix.

Or maybe this truly is a single thing: an Experiment object.

property descriptors

injector_nozzle_area() is lovely -- area is a property of the injector, fully described by the diameter. Consider making it a dynamically computed attribute of the object:

    @property
def nozzle_area(self) -> float:
return (math.pi * self.diameter ** 2) / 4
...
print(my_injector.nozzle_area)


Similarly for relative_length.

I found that the "injector_" prefix was perhaps a bit redundant. Similarly for other "_injector" suffixes and prefixes.

In contrast, average speed and the dimensionless Reynolds number come from the combination of a physical injector and a dynamic flow. Maybe that combination is an OperatingPoint class which holds both? Or maybe just leave them as functions which accept both as args.

piecewise continuous

The linear_hydraulic_resistance() function is a nice piecewise model. It is not at all clear that it is continuous, and that seems worth commenting on in the source code. Or cite a reference, and let that author sweat the details. That way we'd also learn what the cryptic density_comb denotes.

EDIT: Section 2.1.11 explains it is the density of combustion products.

I'm assuming ATM that sigma_fuel is surface tension in N/m, but it would be nice to either explicitly spell that out or cite a reference that does.

lint

Layout of this source code is idiosyncratic, changing from function to function. Recommend you run \$ black *.py against it before sharing it with colleagues.

automated tests

There aren't any. The print()'s are nice enough, but they are not self evaluating. I could not run them, eyeball them, and proclaim, "Oh, yes, those figures are definitely right!".

It's not hard, write a test suite.

correctness

Absent citations, or the URL of a description for the laboratory setup where you observed and recorded fluid properties, it isn't feasible for me to tell whether this code correctly computes desired figures or not. Cite your references, already.

This code achieves some of its design goals.

In its current form I would not be willing to delegate or accept maintenance tasks on it.

• Thanks for the valuable comments! I added a link to the tutorial on which the code was written, I hope this will clarify the difficult points with its understanding Commented Feb 16 at 16:08