Coupled Feed-System and Chamber Modelling of a Gaseous-Oxygen/Kerosene Liquid Rocket Engine: Injector Stiffness, Startup Transients and Low-Frequency Stability

1. Introduction

Small bipropellant rocket engines built by university and amateur teams are almost always designed with a two-stage toolchain: an equilibrium thermochemistry code fixes the chamber conditions and nozzle geometry at a chosen mixture ratio, and a separate hydraulic network model sizes the plumbing that must deliver those propellant flows [1,2]. The two stages are usually coupled only in one direction. The thermochemistry code supplies a design chamber pressure, and the feed model treats that pressure as a fixed downstream boundary condition.

System-level, one-dimensional models of this kind are the established engineering tool for propulsion analysis. Dedicated environments exist for the purpose: ESPSS/EcosimPro is the European standard for space propulsion system simulation [3], and recent reviews survey the transient modelling approaches in use across the field [4]. General-purpose multi-domain libraries such as Simscape Fluids [5] are increasingly applied to the same task, while method-of-characteristics solvers are used where the propellant lines must be resolved as distributed systems [6], and component-level models validated against test data have been reported for the LOX/methane class [7]. Outside propulsion the same philosophy—low-order components assembled into a network and validated against transient measurements—is long established in hydraulic power systems [8], in pipeline transients [9], and in two-phase thermal devices [10].

This is a reasonable engineering simplification at the design point, and it is the standard architecture of the Simscape Fluids feed-system models that have become common in this community. It has a structural consequence that is less widely appreciated: a feed model with a prescribed chamber pressure cannot predict what happens when the feed system fails to deliver the design flow. The chamber pressure it reports is an input, not a result. Injector stiffness computed from it is the stiffness the designer assumed, not the stiffness the hardware will have. Startup transients are whatever ramp was scripted. And low-frequency feed-coupled instability—chug—is by construction impossible, because chug is precisely a feedback loop through the chamber [11,12,13,14].

Low-frequency instability has been treated with combustion time-lag models since Summerfield’s original analysis [11], and double-time-lag formulations remain in use for predicting chugging during start-up in system-level codes [15]. It remains a live concern wherever the injector pressure drop is small relative to chamber pressure, as in deep-throttled engines [16] and in small thrusters operated by pulse-width modulation [17]; analogous low-frequency modes arise in hybrid motors [18] and in cryogenic feed lines [19]. What these analyses have in common is that the chamber pressure is a state of the model, which is precisely what the open-loop feed model gives up.

The distinction matters most when the hardware does not match the design, which for a first-iteration engine is the normal case rather than the exception. The system studied here supplies two concrete instances. Its oxidiser flow turns out to be set by a quarter-inch gate valve that is choked at the operating point, so the oxidiser injector—the component sized to meter that flow—has almost no authority over it. And a cold-flow calibration of the fuel injector, reported in Section 3.2, shows the effective discharge area carried by the model to be 19 to 40 larger than the bench value. Neither discrepancy is visible in a model that holds chamber pressure fixed: the first because the consequence of the restriction is absorbed by the imposed boundary, and the second because the resulting flow error is never propagated into a pressure the rest of the system can respond to.

In this work we take an existing, working Simscape Fluids model of a gaseous-oxygen (GOX)/kerosene feed system for a nominally 100 N engine and replace its prescribed-pressure chamber with a lumped-volume chamber whose pressure is an integrated state. The contribution is threefold. First, we give the coupled formulation and its verification, together with a reproducible build procedure that derives the coupled model from the open-loop one by script rather than by hand. Second, we show what the coupling changes: the predicted operating point moves substantially, and the mechanism by which it moves is a feedback path the open-loop model cannot express. Third, we derive a closed-form low-frequency stability boundary for the coupled system and use it, together with a numerical mechanism-separation procedure, to distinguish genuine feed-coupled chug from a mixture-ratio instability that involves the feed line only incidentally.

2. Materials and Methods

2.1. Engine and Feed System

The engine is a pressure-fed GOX/kerosene thrust chamber designed for $100 \mathrm{N}$ at $p_c=10 \mathrm{bar}$. Chamber geometry follows from an equilibrium analysis at $\mathrm{O}/\mathrm{F}=1.5$: throat diameter $D_t=10.04 \mathrm{m}\mathrm{m}$ ($A_t=79.17 \mathrm{m}{\mathrm{m}}^2$), characteristic length $L^{*}=1.35 \mathrm{m}$ (chamber volume $V_c=L^{*}{A}_t=106.9 \mathrm{c}{\mathrm{m}}^3$), and nozzle area ratio $A_e/A_t=2.256$. The design mass flows are ${\dot{m}}_f=0.01898 \mathrm{k}\mathrm{g} {\mathrm{s}}^{-1}$ and ${\dot{m}}_o=0.02884 \mathrm{k}\mathrm{g} {\mathrm{s}}^{-1}$.

Figure 1 shows the feed system. It is drawn from the active blocks of the executed model: several components present in the Simulink diagram—a poppet main oxidiser valve, a ball valve, a relief valve and several line segments—are commented out in the model as delivered and are therefore absent here, so that the schematic and the simulated system agree.

The oxidizer circuit runs from a 10 $\mathrm{L}$ bottle at 150 $\mathrm{bar}$ through a dome-loaded regulator set to $19.11$ $\mathrm{bar}$, approximately 1 $\mathrm{m}$ of $12.5$ $\mathrm{m}$$\mathrm{m}$ line, a quarter-inch gate valve acting as the main oxidiser valve, a further metre of line, and the oxidiser injector orifice. The fuel circuit is a positive-displacement pump represented as a constant pressure rise of $10.99$ $\mathrm{bar}$, a short 20 $\mathrm{m}$$\mathrm{m}$ line, a solenoid main fuel valve, and the fuel injector orifice. Full component parameters are given in Appendix A.

2.2. Baseline Open-Loop Model

The baseline model is implemented in Simscape Fluids [5] (MATLAB R2023b) using the isothermal-liquid domain for kerosene and the gas domain for oxygen. The combustion chamber is represented by two controlled reservoirs, one in each domain, both driven by a prescribed pressure

$$p_c\left(t\right)=\left(p_{c,\mathrm{des}}-p_a\right)\left[1-exp\left(-{\displaystyle \frac{t-t_{\mathrm{ig}}}{\tau }}\right)\right]+p_a,\tau =0.2 \mathrm{s},$$

(1)

with $p_{c,\mathrm{des}}=10 \mathrm{bar}$. Thrust is reported as $F=\dot{m}{I}_{sp}{g}_0$ with a fixed $I_{sp}$.

2.3. Coupled Chamber Formulation

We replace Equation (1) with a mass balance on a constant-volume chamber holding a perfect gas at the equilibrium flame temperature and discharging through the throat. With $m=p_c{V}_c/\left(R_g{T}_c\right)$ and $\mathrm{d}m/\mathrm{d}t={\dot{m}}_{\mathrm{in}}-{\dot{m}}_{\mathrm{thr}}$,

$${\displaystyle \frac{\mathrm{d}{p}_c}{\mathrm{d}t}}={\displaystyle \frac{R_g{T}_c}{V_c}}\left({\dot{m}}_{\mathrm{in}}-{\dot{m}}_{\mathrm{thr}}\right).$$

(2)

The throat mass flow uses the choked and subsonic branches of the isentropic nozzle relation, with $\mathrm{\Pi }=p_a/p_c$ and ${\mathrm{\Pi }}_{\mathrm{crit}}={[2/(\gamma +1)]}^{\gamma /(\gamma -1)}$:

$${\dot{m}}_{\mathrm{thr}}=\left\{\begin{array}{cc}{\displaystyle \frac{p_c{A}_t}{c^{*}}},\hfill & \mathrm{\Pi }\le {\mathrm{\Pi }}_{\mathrm{crit}},\hfill \\ A_t{p}_c\sqrt{{\displaystyle \frac{2\gamma }{(\gamma -1)R_g{T}_c}}\left({\mathrm{\Pi }}^{2/\gamma }-{\mathrm{\Pi }}^{(\gamma +1)/\gamma }\right)},\hfill & {\mathrm{\Pi }}_{\mathrm{crit}}<\mathrm{\Pi }<1,\hfill \\ 0,\hfill & \mathrm{\Pi }\ge 1.\hfill \end{array}\right.$$

(3)

At steady state Equation (2) recovers the textbook result $p_c=\dot{m}{c}^{*}/A_t$, which we use as the primary verification case in Section 3.1. The associated filling time constant is

$${\tau }_c={\displaystyle \frac{V_c}{c^{*}{\mathrm{\Gamma }}^2{A}_t}}={\displaystyle \frac{L^{*}}{c^{*}{\mathrm{\Gamma }}^2}},\mathrm{\Gamma }=\sqrt{\gamma }{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{2(\gamma -1)}}},$$

(4)

which evaluates to $1.88$ $\mathrm{m}$$\mathrm{s}$ at the design condition—two orders of magnitude faster than the 200 $\mathrm{m}$$\mathrm{s}$ asserted by Equation (1).

Thrust is obtained from the vacuum thrust coefficient with an explicit ambient correction,

$$F=p_c{A}_t{C}_{F,\mathrm{vac}}(\mathrm{O}/\mathrm{F})-p_a{A}_e,$$

(5)

which at the design point returns $99.5$ $\mathrm{N}$ against the $99.33$ $\mathrm{N}$ optimum-expansion value from the design code.

2.3.1. Ignition and Extinction

A step change of the chamber boundary from ambient to flame temperature is both unphysical and numerically intractable: it causes the transient initialisation of the Simscape gas network to fail to converge. We therefore introduce a combustion progress variable $\sigma \in [0,1]$,

$$\sigma \left(t\right)=\left[1-exp\left(-{\displaystyle \frac{t-t_{\mathrm{ig}}}{{\tau }_{\mathrm{ig}}}}\right)\right]g\left({\dot{m}}_f\right)g\left({\dot{m}}_o\right),g\left(\dot{m}\right)=\frac{1}{2}\left[1+tanh3\left(\frac{\dot{m}}{{\dot{m}}_{\min }}-1\right)\right],$$

(6)

with flame-establishment constant ${\tau }_{\mathrm{ig}}=50 \mathrm{m}\mathrm{s}$. Chamber properties are blended between unreacted oxygen and equilibrium combustion by $\sigma$; the smooth gate g makes ignition conditional on both propellants being present, so loss of either extinguishes the chamber and shutdown is a physical blowdown rather than a scripted ramp. Before ignition only the gaseous oxidiser pressurises the chamber, the liquid fuel being blended into the gas-phase source term by $\sigma$.

2.3.2. Thermochemistry

$c^{*}$, $T_c$, $\gamma$ and $C_{F,\mathrm{vac}}$ are interpolated on the instantaneous mixture ratio from an equilibrium table. The gas constant is taken from the tabulated triple as $R_g={\left(c^{*}\mathrm{\Gamma }\right)}^2/T_c$, so that throat flow and chamber filling remain mutually consistent with the tabulated $c^{*}$ rather than being independently parameterised.

The table is generated with NASA CEA [20] for RP-1/${\mathrm{O}}_2$ at $p_c=10\mathrm{bar}$, finite-area combustor with $A_c/A_t=10$, equilibrium composition, at 19 mixture ratios from $\mathrm{O}/\mathrm{F}=0.10$ to $3.20$. The spacing is deliberately dense below $\mathrm{O}/\mathrm{F}=1$ because, as Section 3.3 shows, the as-built feed system drives the engine into that region, and it is the local slope $\mathrm{d}{c}^{*}/\mathrm{d}(\mathrm{O}/\mathrm{F})$ there that governs the behaviour examined in Section 3.5. All 19 cases converged.

Two consistency checks are applied to the parsed table. The molar mass implied by $R_g={\left(c^{*}\mathrm{\Gamma }\right)}^2/T_c$ is compared with the molar mass CEA reports directly, and agrees to within $1.02$ across the full range, which confirms that the chamber station has been identified correctly in every block; the finite-area listing prints both an injector and a combustion-end station, and these differ by only $0.04$ in temperature, so a misidentification would otherwise be silent. Separately, the $\mathrm{O}/\mathrm{F}=1.50$ row reproduces the project’s pre-existing independent CEA run to all printed digits ($c^{*}=1677.9 \mathrm{m} {\mathrm{s}}^{-1}$, $T_c=2575.15 \mathrm{K}$, $\gamma =1.2345$, $M=17.628 \mathrm{k}\mathrm{g} {\mathrm{kmol}}^{-1}$).

2.4. Feed-Line Dynamics

Fluid inertia is enabled in the liquid line, giving the momentum state $I\mathrm{d}\dot{m}/\mathrm{d}t=\Delta p$ with inertance $I=\ell /A$. For the fuel line, $\ell =1.05 \mathrm{m}$ and $A=3.1416\times {10}^{-4} {\mathrm{m}}^2$, giving $I=3342 {\mathrm{m}}^{-1}$. The Simscape gas pipe has no inertance term; this asymmetry is acceptable here because the oxidiser line is choked at the gate valve, which acoustically isolates everything upstream of it from the chamber, and because gas inertance is smaller by orders of magnitude.

This lumped inertance–resistance–capacitance representation is the low-order limit of the distributed water-hammer equations that govern liquid line transients [9,21]. It deliberately does not resolve the unsteady wall friction and geometry corrections that become significant when a line is modelled in detail [22,23], nor the column separation and cavitation that accompany severe transients [24]. The justification is one of scale. With the modelled kerosene properties ($\rho =800\mathrm{k}\mathrm{g} {\mathrm{m}}^{-3}$, $\beta =2.179\mathrm{G}\mathrm{Pa}$) the acoustic speed is 1650 $\mathrm{m}$/$\mathrm{s}$, so the first quarter-wave mode of the $1.05$ $\mathrm{m}$ fuel line lies near 390 $\mathrm{Hz}$—a factor of three above the 130 $\mathrm{Hz}$ phenomena examined in Section 3.5. The lumped inertance therefore places the low-frequency response correctly even though it cannot represent the wave behaviour above it.

The oxidiser main valve is represented by a flow-coefficient ($K_v$) model with a critical pressure ratio, which is the standard treatment for a gas restriction driven into choking; experimental characterisation of choked discharge through a gas restrictor [25] and of the onset of choking in calibrated liquid orifices [26] supports the use of a single coefficient with a switched critical branch, and it is this element that Section 3.3 identifies as controlling the whole engine.

With the momentum state present, the line initial pressure must be consistent with its source: the fuel line is therefore initialised at the pump outlet pressure of $1.2$ $\mathrm{M}$$\mathrm{Pa}$ rather than at ambient. An inconsistent initialisation is harmless in a purely resistive–capacitive line, which relaxes it in a single solver step, but in a dead-ended line with inertia it excites a lightly damped resonance that persists for the whole run, since without through-flow there is almost no dissipation to remove it.

2.5. Cold-Flow Bench Data

The fuel injector was characterised on a cold-flow bench before the model was built. The element is a swirl injector (Mejarra INJ 001) with a documented fuel orifice diameter of $1.33$ $\mathrm{m}$$\mathrm{m}$, giving a geometric area of $1.389$ $\mathrm{m}$${\mathrm{m}}^2$. The working fluid on the bench was water, not kerosene.

Flow was measured with an inline turbine flowmeter and pressure at four taps, logged asynchronously at approximately 1 $\mathrm{Hz}$ each. Two taps ($p_1$, $p_2$) follow the water flow closely (correlation with flow rate $+0.96$ and $+0.91$ over the flowing portion of the record); the remaining two are on a separate circuit that is pressurised only in the second half of the run and are excluded from the analysis.

Because the bench fluid differs from the propellant, the quantity transferred to the model is the effective discharge area $C_{d}A$ rather than $C_d$ or A separately. From

$$\dot{m}=C_{d}A\sqrt{2\rho \Delta p_{\mathrm{inj}}},$$

(7)

$C_{d}A$ is independent of which incompressible fluid was used, whereas $C_d$ and A are not separately identifiable from flow data alone. This is not merely a bookkeeping convenience: the discharge coefficient of a small orifice is not a constant of the geometry but depends on Reynolds number, on inlet chamfer and edge condition, and on whether the orifice cavitates, effects that have been resolved experimentally for circular and chamfered orifices [27] and for injector passages under fuel-injection conditions [28]. Reporting the lumped product avoids attributing the discrepancy to either factor on evidence that cannot distinguish them. Comparable direct measurement of injector mass flow, rather than inference from an assumed coefficient, is reported for self-pressurising propellants in [29].

The flowmeter and the transducers are logged on independent cycles, so during the pump ramp a flow sample can be paired with a pressure sample up to 2 $\mathrm{s}$ stale. We therefore fit only quasi-steady samples, defined as those for which both flow and upstream pressure differ from their immediate neighbours by less than 3 . This retains 66 of 117 aligned samples; Figure 3 shows the record and which samples survived.

2.6. Combustion Time Lag

Propellant crossing the injector face does not release its heat there; it must atomise, vaporise and mix. The chamber therefore burns what was injected ${\tau }_{\mathrm{comb}}$ earlier, implemented as a pure transport delay on the gas-phase source term. This is the $\tau$ of the Crocco–Cheng n–$\tau$ formulation [11,12]. A pure delay is used rather than a first-order lag because the destabilising mechanism depends on the unbounded phase roll-off of a true delay. We treat ${\tau }_{\mathrm{comb}}$ as a swept parameter rather than predicting it; constant- and variable-lag closures for the same quantity, and their effect on predicted chug boundaries, are compared in [15].

2.7. Numerical Implementation

The coupled model is generated from the open-loop model by a build script, which replaces the chamber subsystem, writes all parameters into the model workspace, and configures the solver. The model is therefore self-contained and the derivation reproducible. Simulations use ode23t with relative and absolute tolerances of ${10}^{-6}$ and a maximum step of 1 $\mathrm{m}$$\mathrm{s}$, chosen so that the $1.88$ $\mathrm{m}$$\mathrm{s}$ filling constant of Equation (4) is resolved.

3. Results

3.1. Verification

The steady operating point of the coupled model reproduces $p_c={\dot{m}}_{\mathrm{tot}}{c}^{*}/A_t$ to better than ${10}^{-3}$%, evaluated independently from the interpolated $c^{*}$ at the simulated mixture ratio. The thrust relation, Equation (5), returns $99.5$ $\mathrm{N}$ at the design point against $99.33$ $\mathrm{N}$ from the design code, a $0.2$ agreement that confirms the $C_{F,\mathrm{vac}}$ convention and the ambient correction.

3.2. Validation Against Cold-Flow Data

Figure 2 shows the bench data fitted to Equation (7). The discharge law holds well over the range covered: $R^2=0.962$ with a mean absolute residual of $0.68$ and all points inside $\pm 5\%$ on the parity plot.

Figure 3. Bench record of the fuel injector run, showing which samples passed the quasi-steady filter (66 of 117). Flow on the left axis, upstream pressure $p_1$ on the right.

Figure 3. Bench record of the fuel injector run, showing which samples passed the quasi-steady filter (66 of 117). Flow on the left axis, upstream pressure $p_1$ on the right.

The instrumentation does not identify unambiguously which tap pair bounds the injector, so both plausible readings are reported in Table 1. They bracket the answer, and the conclusion is the same either way: the model over-predicts the fuel injector effective discharge area by 19 to 40 .

The model reaches its $C_{d}A$ through a physically implausible pair of values: an area of $3.892$ $\mathrm{m}$${\mathrm{m}}^2$ (equivalent diameter $2.23$ $\mathrm{m}$$\mathrm{m}$, against a documented $1.33$ $\mathrm{m}$$\mathrm{m}$) with a compensating $C_d$ of 0.288. Only the product enters the flow calculation, so this does not by itself invalidate the model, but it means the individual parameters cannot be interpreted physically and the agreement with the bench is worse than the plausible-looking $C_d$ suggests.

Two limitations bound how much this calibration establishes, and both are limitations of the available measurement rather than of the fit. First, the fitted range spans a factor of about 1.9 in $\Delta p$, so the square-root law is exercised but not strongly tested; the result is closer to a calibration at one operating point with scatter than to a demonstration of the exponent. Second, the quoted intervals are bootstrap estimates of precision only. The dominant uncertainty is systematic: the turbine flowmeter K-factor was taken from the manufacturer’s nominal value and neither it nor the pressure transducer calibrations were verified against a reference, so no total uncertainty can be quoted. An independent check—a timed gravimetric catch, or a calibration certificate traceable to a standard—would be required before the interval widths in Table 1 could be read as accuracy.

We nonetheless regard the direction of the discrepancy as established, because it is far larger than any plausible instrument error: closing a 19–40 gap in $C_{d}A$ by calibration error alone would require the flowmeter to be in error by a similar proportion, an order of magnitude beyond the specification of a turbine meter in its rated range. The magnitude, by contrast, should be treated as provisional.

The oxidiser circuit is not validated at all. The oxidiser bench run logged pressures but no flow rate, so no effective discharge area can be extracted from it, and the oxidiser injector area in the model remains the as-designed value carried over from the sizing calculation. This matters less than it otherwise would, for the reason developed in Section 3.3: the oxidiser flow is fixed by choking at the gate valve upstream, so the injector area has almost no influence on it, and the predicted oxidiser flow is insensitive to precisely the parameter that is unvalidated. That argument holds only while the valve remains choked. An engine of this design operated at a higher chamber pressure, or with the gate valve replaced, would move the metering authority back to the injector and would require the oxidiser cold-flow measurement before its predictions could be trusted.

3.2.1. Effect on the Coupled Prediction

Substituting the bench-calibrated $C_{d}A$ into the coupled model, holding $C_d$ at its model value and scaling the area, gives Table 2.

Correcting the injector reduces fuel flow by up to 24 and raises the mixture ratio from 0.252 to 0.333, but the engine remains severely oxidiser-starved and far below its design chamber pressure. The headline conclusion is therefore robust to the injector calibration: it is driven by the choked oxidiser valve, not by the fuel-side parameters. The oxidiser flow is unchanged in all three cases, as expected for a circuit whose flow is set by a choked restriction upstream.

3.3. Effect of Coupling on the Predicted Operating Point

Table 3 compares the open-loop and coupled models against the design point.

Both models agree that the engine is oxidiser-starved, but they disagree on the mechanism and the magnitude. In both, the oxidiser flow is set by the quarter-inch gate valve, which is choked: the measured pressure ratio across it is 0.506 against a critical value of 0.528 for $\gamma =1.4$. The oxidiser injector, in contrast, drops only $0.175$ $\mathrm{bar}$ of the available 20 $\mathrm{bar}$. The circuit is valve-limited, not injector-limited, and the injector consequently has almost no authority over the chamber.

The open-loop model stops there. The coupled model additionally captures the consequence: with $p_c$ free to fall, the fuel injector sees a larger pressure difference and delivers more fuel, which drives the mixture ratio further down, which lowers $c^{*}$, which lowers $p_c$ again. Figure 4 shows the resulting divergence between the two models.

Two further quantities follow from the coupled formulation and have no open-loop counterpart at all, because both are defined against a chamber pressure the open-loop model does not compute (Figure 5).

Figure 6 and Figure 7 show the ignition and shutdown transients resolved by the coupled model, neither of which is available from a prescribed-pressure formulation.

3.4. Linear Stability of the Coupled System

Linearising the fuel line and chamber about a steady point, with the oxidiser treated as constant because its supply is choked, gives

$$\begin{array}{cccc}\hfill I{\displaystyle \frac{\mathrm{d}\delta {\dot{m}}_f}{\mathrm{d}t}}& =-\delta p_c-R\delta {\dot{m}}_f,\hfill & \hfill R& ={\displaystyle \frac{2\Delta p_{\mathrm{inj}}}{{\dot{m}}_f}},\hfill \end{array}$$

(8)

$$\begin{array}{cccc}\hfill C_c{\displaystyle \frac{\mathrm{d}\delta p_c}{\mathrm{d}t}}& =\delta {\dot{m}}_f(t-{\tau }_{\mathrm{comb}})-G\delta p_c,\hfill & \hfill C_c& ={\displaystyle \frac{V_c}{R_g{T}_c}},G={\displaystyle \frac{A_t}{c^{*}}},\hfill \end{array}$$

(9)

whence the characteristic equation

$$\left(C_{c}s+G\right)\left(Is+R\right)+e^{-s{\tau }_{\mathrm{comb}}}=0.$$

(10)

With ${\tau }_{\mathrm{comb}}=0$ this reduces to $s^{2}I{C}_c+s(C_{c}R+GI)+(GR+1)=0$, whose coefficients are all positive: the system without a combustion lag is unconditionally stable, and enabling fluid inertia alone cannot produce chug.

On the imaginary axis $s=\mathrm{j}\omega$, writing $a=C_{c}I$, $b=C_{c}R+GI$ and $c=GR$, Equation (10) separates into $a{\omega }^2-c=cos\omega {\tau }_{\mathrm{comb}}$ and $b\omega =sin\omega {\tau }_{\mathrm{comb}}$, so that

$$a^2{\omega }^4+\left(b^2-2ac\right){\omega }^2+\left(c^2-1\right)=0.$$

(11)

A positive root—hence a stability boundary at finite ${\tau }_{\mathrm{comb}}$—exists whenever $c=GR<1$. Substituting $G={\dot{m}}_{\mathrm{tot}}/p_c$ and $R=2\Delta p_{\mathrm{inj}}/{\dot{m}}_f$ gives the compact result

$$GR=2{\displaystyle \frac{\Delta p_{\mathrm{inj}}}{p_c}}\left(1+\mathrm{O}/\mathrm{F}\right),\mathrm{so}\mathrm{instability}\mathrm{requires}\begin{array}{|c|}\hline {\displaystyle \frac{\Delta p_{\mathrm{inj}}}{p_c}}<{\displaystyle \frac{1}{2\left(1+\mathrm{O}/\mathrm{F}\right)}}\\ \hline\end{array}$$

(12)

At $\mathrm{O}/\mathrm{F}=1.52$ this evaluates to $19.8$ . The classical design rule that the injector should drop of order twenty percent of chamber pressure [1,2] thus emerges from the lumped model rather than being imposed on it, in the same spirit as the time-lag stability analyses that established it [11,12] and consistent with the pressure-drop margins reported as limiting in throttled engines [16]. Critical lags and frequencies along the boundary are given in Table 4.

3.5. Mechanism Separation

Introducing the combustion lag produced large-amplitude oscillation for ${\tau }_{\mathrm{comb}}\ge 2.5\mathrm{m}\mathrm{s}$ at every injector stiffness tested, including 125 —far inside the stable region predicted by Equation (12)—and at a frequency that depended only on ${\tau }_{\mathrm{comb}}$ and not on the feed line. Both observations are inconsistent with feed-coupled chug.

To identify the mechanism we repeated the sweep with the thermochemistry lookups frozen at their operating-point values, so that $\mathrm{d}{c}^{*}/\mathrm{d}(\mathrm{O}/\mathrm{F})=0$, and separately with fluid inertia disabled (Table 5, Figure 8).

Freezing $c^{*}$ removes the oscillation by five orders of magnitude, to the numerical noise floor. Removing fluid inertia changes it by less than $0.8$ . The instability is therefore driven entirely by the $c^{*}(\mathrm{O}/\mathrm{F})$ gain and not at all by feed inertia: it is not chug. The loop is a mixture-ratio oscillation—rising chamber pressure reduces fuel flow by back-pressure, raising O/F, raising $c^{*}$, raising chamber pressure further—closed by the combustion lag, with the feed line an incidental participant.

Genuine chug is correctly absent, because the engine’s fuel injector is far stiffer than the boundary of Equation (12) requires.

4. Discussion

4.1. A Prescribed Chamber Pressure Is Not a Boundary Condition

The most consequential finding is structural rather than numerical. Prescribing $p_c$ does not merely fix one variable; it severs the only path by which the feed system learns that it is not delivering the design flow. The open-loop model of this engine reports a fuel injector stiffness of $20.1$ , almost exactly the design intent, while the coupled model of the same hardware reports 125 . Neither number is a modelling error: they are answers to different questions. The first is “what stiffness would the injector have if the chamber reached 10 $\mathrm{bar}$?”; the second is “what stiffness does it have?”.

This matters because injector stiffness is normally computed precisely to decide whether an engine is at risk of feed-coupled instability, and the open-loop answer is systematically optimistic in exactly the situation where the question is worth asking—when the chamber pressure falls short. A model that assumes the design chamber pressure will report the design stiffness almost regardless of the hardware.

Section 3.5 illustrates a diagnostic that is cheap and underused in lumped multi-physics models. When a coupled model oscillates, the oscillation is usually attributed to whichever loop the modeller was thinking about. Here the natural attribution—feed-coupled chug, since the feed line carries a momentum state and the chamber a pressure state—was wrong, and two observations gave it away before any further computation: the instability persisted at injector stiffnesses far inside the analytically stable region, and its frequency was insensitive to the feed line.

Freezing a constitutive relation at its operating-point value, so that its gain is removed while everything else is untouched, then settles the question directly. It requires no new formulation, and unlike disabling a component it does not perturb the steady state. We suggest it as a routine check whenever a coupled model produces an instability whose mechanism is being inferred rather than demonstrated.

4.2. Sensitivity to the Thermochemistry Table

Because the coupled model drives the engine to a mixture ratio near 0.25, far from the design value, it operates on a part of the $c^{*}(\mathrm{O}/\mathrm{F})$ curve that a table built around the design point would represent poorly. During development we ran the model against an interpolated table anchored on the single design-point equilibrium result available at the time, before replacing it with the full nineteen-point sweep of Section 3.3. The comparison is instructive, because it separates which conclusions are robust to the thermochemistry and which are not.

The interpolated table was in error by up to 24 in $c^{*}$ and by nearly a factor of two in $T_c$ at the fuel-rich end, and—more importantly—gave a local slope $\mathrm{d}{c}^{*}/\mathrm{d}(\mathrm{O}/\mathrm{F})$ at the operating point about 65 larger than the computed one. Since Section 3.5 identifies that slope as the gain closing the oscillation loop, this was the one quantity most likely to overturn a conclusion.

It did not. Replacing the table moved the operating point modestly—chamber pressure from 4.96 to $5.34$ $\mathrm{bar}$, thrust from 41.1 to $45.4$ $\mathrm{N}$, mixture ratio from 0.245 to 0.252—and moved the oscillation amplitude at ${\tau }_{\mathrm{comb}}=2.5\mathrm{m}\mathrm{s}$ from $30.7$ to $26.4$ , consistent with the reduced gain. The onset lag, the frequencies, and every qualitative conclusion were unchanged: the engine remains severely oxidiser-starved at roughly half its design chamber pressure, and the oscillation remains a mixture-ratio instability that survives removal of feed inertia and vanishes when the thermochemistry is frozen.

The lesson we would draw is not that the table did not matter, but that its influence was confined to the quantities one would expect it to set. A reader adapting this approach should establish the same separation explicitly rather than assume it, particularly if their operating point sits on a steeper part of the $c^{*}$ curve than this one does.

4.3. Limitations

Beyond the thermochemistry table, three limitations bear directly on the results reported here. The chamber is a single well-stirred volume, so it cannot represent axial or transverse acoustic modes and the stability analysis is confined to the low-frequency regime; the intermediate- and high-frequency regimes require the formulations collected in [13,14]. The liquid line is represented by a single lumped inertance and capacitance, which places its resonance correctly in order of magnitude but does not resolve distributed wave behaviour. And the oxidiser side carries no inertance at all, a limitation of the gas-domain pipe model; this is defensible here only because the circuit is choked at the gate valve, and would not be for an engine whose oxidiser injector had genuine authority.

Four further simplifications bound how far the approach generalises to other hardware. The nozzle enters only through a tabulated $C_{F,\mathrm{vac}}$ and the exit area, so three-dimensional effects—separation, side loads, plume structure—lie outside the formulation [30,31], and at small scale the measured performance of a convergent–divergent nozzle can depart appreciably from its one-dimensional value [32]. The chamber is adiabatic: no heat is transferred to the wall or to a coolant circuit, which for a regeneratively cooled engine is a first-order coupling rather than a correction [33]. The propellant sources are idealised as a constant pressure rise and a bottle with prescribed blowdown, so tank ullage dynamics and sloshing are absent [34]. And the liquid is a single non-flashing phase, appropriate for kerosene at these conditions but not for a cryogenic propellant near saturation, where flashing and phase change dominate the line behaviour [35].

Finally, the experimental support is one-sided. The fuel injector is calibrated against bench data; the oxidiser injector is not, because the oxidiser bench run recorded no flow rate, and the fuel calibration itself carries no traceable uncertainty. We have argued in Section 3.2 that the conclusions survive this—the oxidiser flow is set by choking upstream of the injector, and the fuel-side discrepancy is too large to be instrumental—but the argument is one about robustness, not a substitute for the measurement. A gravimetric cross-check of the fuel bench and an oxidiser cold-flow run with mass flow instrumentation are the two experiments that would convert the present validation from partial to complete, and they are the natural next step for this hardware.

5. Conclusions

Replacing a prescribed chamber pressure with a lumped chamber whose pressure is a state converts a propellant feed model from a hydraulic bookkeeping exercise into one that can answer the questions such models are built for. For the engine studied here the two formulations disagree not marginally but qualitatively: the open-loop model reports a fuel injector stiffness of $20.1$ , essentially the design intent, while the coupled model of the same hardware reports 125 . Both are correct answers to different questions, and only the second is about the hardware.

The following results do not depend on the thermochemistry table and are considered established.

1.

The coupled formulation reproduces the analytical steady state $p_c=\dot{m}{c}^{*}/A_t$ to better than ${10}^{-3}$%, and the thrust relation agrees with the independent design code to $0.2$ . The chamber filling time constant is $L^{*}/\left(c^{*}{\mathrm{\Gamma }}^2\right)=1.88\mathrm{m}\mathrm{s}$, two orders of magnitude faster than the 200 $\mathrm{m}$$\mathrm{s}$ the open-loop formulation asserted.

2.

A linear analysis of the coupled feed line and chamber gives a closed-form low-frequency stability boundary, $\Delta p_{\mathrm{inj}}/p_c<{\left[2(1+\mathrm{O}/\mathrm{F})\right]}^{-1}$, which evaluates to $19.8$ at the design mixture ratio. The classical twenty-percent injector-stiffness rule therefore emerges from the lumped model rather than being imposed on it. Feed inertia and injector resistance alone give an unconditionally stable characteristic equation: a combustion time lag is necessary for chug, and enabling inertia by itself cannot produce it.

3.

Freezing a constitutive relation at its operating-point value is an effective and inexpensive way to identify which loop drives an instability in a coupled model. Applied here it showed that the oscillation observed on introducing a combustion lag is driven entirely by the $c^{*}(\mathrm{O}/\mathrm{F})$ gain and not at all by feed inertia, which changes it by $0.19$ : it is a mixture-ratio instability, not chug.

4.

Cold-flow calibration of the fuel injector gives an effective discharge area of 0.80 to $0.94$ $\mathrm{m}$${\mathrm{m}}^2$ against the $1.119$ $\mathrm{m}$${\mathrm{m}}^2$ carried by the model, an over-prediction of 19 to 40 . Correcting it reduces the predicted fuel flow by up to 26 without altering the qualitative conclusion, which is governed by the choked oxidiser valve rather than by fuel-side parameters. The direction of this discrepancy is secure—it is far larger than any credible instrument error—but its magnitude rests on an uncalibrated flowmeter, and the oxidiser injector is not validated at all.

The engine-specific operating point follows: the engine reaches $p_c=5.34\mathrm{bar}$ against a design value of 10 $\mathrm{bar}$, at a mixture ratio of 0.252 against 1.52, producing $45.4$ $\mathrm{N}$ against $99.3$ $\mathrm{N}$. The oxidiser circuit is the cause: its flow is fixed by a choked gate valve and is insensitive to every other parameter varied in this study, including a $\pm 40$ change in the fuel injector discharge area and the replacement of the thermochemistry table. Restoring the design point requires opening that restriction by roughly a factor of three in $K_v$; no fuel-side change achieves it.