**1. Introduction**

While new applications and constraints have arisen in the field of space transportation in the last decade, hybrid propulsion has been actively investigated worldwide for its benefits in comparison with other chemical engines. The technology is mature enough that the space industry is starting to invest and develop applications based on it. These engines offer significant advantages [1]: they are green propellant-based and have good propulsive performances with a relative simplicity and low costs. The safety level is higher than with other types of chemical engines since the propellants are non-explosive and stored separately. Moreover, hybrid engines can be throttled, reused, and ignited or extinguished several times.

Despite these benefits, hybrid engines suffer from some disadvantages [1] that can be addressed with more or less difficulty. The fundamental difference with other types of chemical propulsion comes from the nature of the combustion which occurs through a diffusion flame located in the turbulent boundary layer close to the fuel grain surface [2]. As a consequence, the fuel regression rate in hybrid motors depends on the heat transfers at the fuel grain surface, which are not directly controllable, and is a function of the mass flux in the combustion chamber. The fuel regression rate in hybrid engines is generally low, which means it may be difficult to obtain high thrust levels. The diffusion-limited combustion is the source of relatively low combustion efficiency due to a lower degree of propellants mixing. Moreover, the increase of the grain port diameter during the combustion causes an O/F

(oxidizer to fuel ratio) shift and some performance modifications. Finally, like all chemical rocket engines, hybrid motors can suffer from pressure oscillations and instabilities that can be provoked by different phenomena such as combustion, acoustics, hydrodynamics or hybrid intrinsic behavior.

These fundamental issues have been investigated and solutions do exist depending on the mission and engine configuration. Although a low regression rate may be interesting for long duration and low thrust missions, a higher regression rate can be obtained by using liquefying fuels [3,4]. By increasing the residence time of propellants and improving their mixing, swirl oxidizer injection enhances the combustion efficiency making the hybrid engine competitive regarding that aspect [5,6]. Moreover, the vortex motion of the flow tends to stabilize the combustion and reduces pressure oscillations. In order to overcome the O/F shift due to the fuel regression, several options are possible such as maintaining a constant fuel-burning surface during combustion [7], adding an oxidizer aft-injection [8,9] or using an A-SOFT (Altering-intensity Swirling-Oxidizer-Flow-Type) engine [10] described in Figure 1.

**Figure 1.** A-SOFT (Altering-intensity Swirling-Oxidizer-Flow-Type) concept with dual front oxidizer injections [10].

The oxidizer mass flow rate influences the regression rate but also plays a major role in the thrust of the engine. As a consequence, to properly manage both the thrust and the O/F ratio in hybrid rockets, two control variables are required. In A-SOFT engines, these variables are the total oxidizer mass flow rate and the effective swirl number, which can be viewed as a ratio between axial and radial oxidizer momentums or more simply as the swirling intensity of the flow. In practice, these variables are manipulated through the use of axial and swirl injectors. The swirl oxidizer injection increases the regression rate compared to the non-swirled case [6,11] and an A-SOFT engine uses this effect to control the fuel regression rate and hence the O/F ratio. An ideal A-SOFT motor should be capable of managing the thrust and O/F ratio independently and instantaneously in a feedback control loop to adapt to any required thrust profile while maintaining an optimal operating O/F ratio. The concept has been studied within the Hybrid Rocket Research Working Group in Japan and theoretical [10,12–14], numerical [15,16], and experimental [5,11,17,18] investigations have been conducted to prove the feasibility and interest of such motors. To move forward in the development of A-SOFT engines, two major issues have to be addressed. The first one is related to sensing and measurements and the second one concerns the design of the control law.

To operate properly, an A-SOFT motor would require in-flight, instantaneous, precise, and robust measurements of the thrust and O/F ratio that are the variables to be regulated by a feedback loop control. Although the thrust cannot be measured directly in flight conditions, it can be estimated by using external parameters such as the rocket's acceleration and trajectory and/or by using the motor's data such as the combustion chamber pressure and the propellants mass flow rates. While the oxidizer mass flow rate is not difficult to measure, it is more complicated for the fuel. The main difficulty comes from the necessity of getting instantaneous measurements that are required for the feedback control in flight conditions. One of the solutions that is currently considered is to use RBS (Resistor-Based

Sensors) that are embedded in the fuel grain and deliver an electric voltage proportional to their length. The sensors burn at the same rate as the fuel and allow for measuring the fuel regression rate and mass flow rate while the engine is operating. The SPLab (Space Propulsion Laboratory) from Politecnico di Milano developed its own RBS technology [19,20]. JAXA and SPLab have initiated a collaboration in order to study these types of measurements for A-SOFT applications [21].

A closed-loop throttling of a hybrid rocket has already been demonstrated based on thrust or pressure feedback using proportional/integral control algorithms [22]. Although the results are interesting, the O/F ratio evolution was not considered at all, which is a major difference with the A-SOFT engine concept. Regulating both the thrust and the O/F ratio instantaneously and independently significantly increases the complexity of the feedback control. Due to the complex physical phenomena and their respective interactions in a combustion chamber, an efficient feedback control law must be developed for A-SOFT hybrid motors.

In that context, the goal of this work is twofold: firstly, we investigate the effect of measurement errors on the engine control by conducting an error propagation analysis regarding the thrust and the O/F ratio. Secondly, we develop a regulation law suitable for A-SOFT engine applications. In addition, we perform numerical simulations of a simplified hybrid motor using a feedback control loop to confirm that the regulation law is valid, in a case where RBS are used for the fuel regression rate estimation. From the error propagation analysis, we demonstrate that it is possible to determine the requirements of measurement precision to fit a desired precision on the thrust and O/F ratio control. Moreover, we propose two methods to help design the engine and the RBS based on propagation analysis. Finally, we perform numerical simulations with simple thrust profiles and confirm the efficiency of the proposed regulation law for these configurations.

#### **2. Materials and Methods**

#### *2.1. Fundamental Equations of A-SOFT Engines*

The regression rate law proposed by Marxman [2] is largely used to analyze hybrid motors. However, in the case of a swirl oxidizer injection, an extended version of this law can be used [11] and is given by Equation (1):

$$
\dot{r}[t] = f[\mathcal{S}\_\mathfrak{e}[t]] \mathcal{g}[\mathcal{G}\_\mathfrak{o}[t]].\tag{1}
$$

While *g* is the classical regression rate law depending on the oxidizer mass flux (Equations (2)–(4)), *f* is an extended function depending on the effective swirl number (Equations (5) and (6)):

*g*[*Go*[*t*]] = *aG<sup>n</sup> o* [*t*], (2)

$$\mathcal{G}\_o[t] = \frac{4}{\pi} \frac{\dot{m}\_o[t]}{\phi^2[t]},\tag{3}$$

$$
\dot{m}\_o[t] = \dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]\_\prime \tag{4}
$$

$$f[S\_{\mathfrak{f}}[t]] = \left(1 + S\_{\mathfrak{e}}^2[t]\right)^m,\tag{5}$$

$$S\_{\varepsilon}[t] = S\_{\mathcal{S}} \frac{\dot{m}\_{oT}^{2}[t]}{\left(\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]\right)^{2}}.\tag{6}$$

The geometric swirl number *Sg* depends on the geometry of the injectors and remains constant during the engine operation. Finally, the fuel regression rate can be expressed by Equation (7):

$$\dot{\tau}[t] = 4^n \pi^{-n} a \left( 1 + S\_{\mathcal{S}}^2 \frac{\dot{m}\_{oT}^4[t]}{\left(\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]\right)^4} \right)^m \left( \frac{\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]}{\phi^2[t]} \right)^n. \tag{7}$$

The fuel mass flow rate can be expressed by Equation (8), with *C*1 = <sup>4</sup>*<sup>n</sup>π*1−*naLρf* .

$$\dot{m}\_f[t] = \mathbb{C}\_1 \phi[t] \left( 1 + S\_g^2 \frac{\dot{m}\_{oT}^4[t]}{\left(\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]\right)^4} \right)^m \left( \frac{\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]}{\phi^2[t]} \right)^n. \tag{8}$$

The mixture ratio is then estimated by Equation (9):

$$\begin{split} \xi[t] &= \frac{\dot{m}\_o[t]}{\dot{m}\_f[t]} \\ &= \frac{\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]}{C\_1 \phi[t] \left(1 + S\_g^2 \frac{\dot{m}\_{oT}^4[t]}{\left(\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]\right)^4}\right)^m \left(\frac{\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]}{\phi^2[t]}\right)^n} . \end{split} \tag{9}$$

The thrust is evaluated by Equation (10):

$$\begin{split} \mathcal{F}[t] &= I\_{\Phi}[\mathcal{S}]g\_{0} \left\{ \boldsymbol{m}\_{oA}[t] + \boldsymbol{m}\_{oT}[t] + \boldsymbol{m}\_{f}[t] \right\} \\ &= I\_{\Phi}[\mathcal{S}]g\_{0} \left\{ \boldsymbol{m}\_{oA}[t] + \boldsymbol{m}\_{oT}[t] + \mathcal{C}\_{1}\boldsymbol{\Phi}[t] \left( 1 + \mathcal{S}\_{3}^{2} \frac{\boldsymbol{m}\_{oT}^{4}[t]}{\left( \dot{m}\_{oA}[t] + \dot{m}\_{oT}[t] \right)^{4}} \right)^{m} \left( \frac{\dot{m}\_{oA}[t] + \dot{m}\_{oT}[t]}{\boldsymbol{\Phi}^{2}[t]} \right)^{n} \right\}. \end{split} \tag{10}$$

The dependence of the specific impulse on the chamber pressure is generally of second order compared to the mixture ratio influence and has been neglected in this analysis. Consequently, the specific impulse is considered as a function of the mixture ratio only and is given by Equation (11). In that study, we consider constant nozzle and combustion efficiencies and no nozzle throat erosion:

$$I\_{sp}[\xi] = \eta\_c \eta\_n I\_{sp,th}[\xi]. \tag{11}$$

As a consequence from the previous development, it is possible to express the thrust and the mixture ratio as two functions of three independent variables: *m*˙ *oA*, *m*˙ *oT*, and *φ* (Equation (12)):

$$
\begin{pmatrix} F \\ \zeta^{\mathbb{I}} \end{pmatrix} = \begin{pmatrix} \mathbb{F} \left[ \dot{m}\_{oA}, \dot{m}\_{oT}, \phi \right] \\ \zeta^{\mathbb{I}} \left[ \dot{m}\_{oA}, \dot{m}\_{oT}, \phi \right] \end{pmatrix} . \tag{12}
$$

#### *2.2. Error Propagation of A-SOFT Engines*

In this study, errors refer to RMS (Root Mean Square) relative errors. The RMS relative error of a function *h* depending on three independent variables [*x*, *y*, *z*] is given by Equation (13) (by using the log function properties: *∂x* = *x∂* ln *x*, and: *∂h*/*h* = *∂* ln *h*):

$$c\_{h} = \frac{\delta h}{h} = \frac{\sqrt{\left(\frac{\partial h}{\partial x} \delta x\right)^2 + \left(\frac{\partial h}{\partial y} \delta y\right)^2 + \left(\frac{\partial h}{\partial z} \delta z\right)^2}}{h}$$

$$= \sqrt{\left(\frac{\partial h}{\partial x} \frac{\delta x}{h}\right)^2 + \left(\frac{\partial h}{\partial y} \frac{\delta y}{h}\right)^2 + \left(\frac{\partial h}{\partial z} \frac{\delta z}{h}\right)^2} \tag{13}$$

$$= \sqrt{\left(\frac{\partial \ln h}{\partial \ln x} c\_x\right)^2 + \left(\frac{\partial \ln h}{\partial \ln y} c\_y\right)^2 + \left(\frac{\partial \ln h}{\partial \ln z} c\_z\right)^2}.$$

*Aerospace* **2019**, *6*, 65

By using this definition, the relative errors of thrust and O/F ratio can be written by Equations (14) and (15):

$$
\varepsilon\_F = \sqrt{\left(\frac{\partial \ln F}{\partial \ln \dot{m}\_{\phi A}} \varepsilon\_{\dot{m}\_{\phi A}}\right)^2 + \left(\frac{\partial \ln F}{\partial \ln \dot{m}\_{\phi T}} \varepsilon\_{\dot{m}\_{\phi T}}\right)^2 + \left(\frac{\partial \ln F}{\partial \ln \phi} \varepsilon\_{\phi}\right)^2} \tag{14}
$$

$$\varepsilon\_{\xi} = \sqrt{\left(\frac{\partial \ln \zeta}{\partial \ln \dot{m}\_{oA}} \varepsilon\_{\dot{m}\_{oA}}\right)^2 + \left(\frac{\partial \ln \zeta}{\partial \ln \dot{m}\_{oT}} \varepsilon\_{\dot{m}\_{oT}}\right)^2 + \left(\frac{\partial \ln \zeta}{\partial \ln \phi} \varepsilon\_{\Phi}\right)^2}. \tag{15}$$

To simplify the relations, we will consider that the oxidizer Mass Flow Rate (MFR) errors are equal: *em*˙ *oT* = *em*˙ *oA* = *eMFR*. Moreover, in order to perform a comparative analysis of the results, we define the thrust, O/F, diameter, and total normalized errors (Equations (16)–(19)) as the following:

$$
\overline{\mathcal{e}\_F} = \frac{\mathcal{e}\_F}{\mathcal{e}\_{MFR}} \, \tag{16}
$$

$$
\overline{c\_{\xi}} = \frac{c\_{\xi}}{\mathcal{e}\_{MFR}},
\tag{17}
$$

$$\overline{\mathcal{c}\_{\Phi}} = \frac{\mathcal{c}\_{\Phi}}{\mathcal{c}\_{MFR}} \, \tag{18}$$

$$\overline{e\_{total}} = \frac{\sqrt{e\_F^2 + e\_{\overline{\xi}}^2}}{e\_{MFR}}.\tag{19}$$

Consequently, the thrust, O/F and total normalized errors are written in Equations (20)–(22):

$$
\overline{\varepsilon\_F} = \sqrt{\left(\frac{\partial \ln F}{\partial \ln \dot{m}\_{oA}}\right)^2 + \left(\frac{\partial \ln F}{\partial \ln \dot{m}\_{oT}}\right)^2 + \left(\frac{\partial \ln F}{\partial \ln \phi}\right)^2 \overline{\varepsilon\_\phi}^2},\tag{20}
$$

$$\overline{\mathcal{E}\_{\xi}^{2}} = \sqrt{\left(\frac{\partial \ln \zeta\_{\flat}^{x}}{\partial \ln \dot{m}\_{\flat A}}\right)^{2} + \left(\frac{\partial \ln \zeta\_{\flat}^{x}}{\partial \ln \dot{m}\_{\flat T}}\right)^{2} + \left(\frac{\partial \ln \zeta\_{\flat}^{x}}{\partial \ln \phi}\right)^{2} \overline{\mathcal{E}\_{\Phi}^{2}}},\tag{21}$$

$$
\overline{\overline{\varepsilon\_{total}}} = \sqrt{\overline{\varepsilon\_{\overline{\Gamma}}^2} + \overline{\varepsilon\_{\overline{\zeta}}^2}}.\tag{22}
$$

As it can be seen, these normalized errors are functions of the ratio *eφ* = *<sup>e</sup>φ*/*eMFR*. The oxidizer mass flow rate measurements precision determines *eMFR* while *eφ* depends on the fuel regression rate sensors' precision. We introduce the parameter *αT* defined as the ratio of the oxidizer injections and varying in [0, 1] (0 if: *m*˙ *oT* = 0 and 1 if: *m*˙ *oA* = 0) and the function Φ in Equations (23) and (24):

$$
\hbar a\_T = \frac{\dot{m}\_{oT}}{\dot{m}\_{oA} + \dot{m}\_{oT}} \,\tag{23}
$$

$$
\Phi[\alpha\_T] = 1 + S\_{\mathcal{K}}^2 \alpha\_T^4. \tag{24}
$$

The sensitivity coefficients involved in the expressions of the normalized errors can be calculated based on Equations (9) and (10), and are given in the following relations (Equations (25)–(30)):

$$\frac{\partial \ln F}{\partial \ln \dot{m}\_{oA}} = (1 - a\_T) \left\{ \frac{n + \frac{\pi}{\xi}}{1 + \frac{\pi}{\xi}} - \frac{4m S\_g^2 a\_T^4}{(1 + \frac{\pi}{\xi}) \Phi[a\_T]} + \left( 1 - n + \frac{4m S\_g^2 a\_T^4}{\Phi[a\_T]} \right) \frac{d \ln I\_{sp}}{d \ln \xi} [\xi] \right\},\tag{25}$$

$$\frac{\partial \ln F}{\partial \ln \dot{m}\_{\text{oT}}} = a\_T \left\{ \frac{n + \frac{\pi}{5}}{1 + \frac{\pi}{5}} - \frac{4mS\_{\text{g}}^2 (a\_T - 1)}{(1 + \xi)} \frac{a\_T^3}{\Phi[a\_T]} + \left( 1 - n + \frac{4mS\_{\text{g}}^2 (a\_T - 1) a\_T^3}{\Phi[a\_T]} \right) \frac{d \ln I\_{\text{sp}}}{d \ln \frac{\pi}{5}}[\xi] \right\},\tag{26}$$

$$\frac{\partial \ln F}{\partial \ln \phi} = (2n - 1) \left\{ -\frac{1}{1 + \frac{\mathcal{I}}{\mathcal{J}}} + \frac{d \ln I\_{sp}}{d \ln \mathcal{J}}[\mathcal{\xi}] \right\},\tag{27}$$

$$\frac{\partial \ln \xi}{\partial \ln \dot{m}\_{oA}} = (1 - a\_T) \left\{ 1 - n + \frac{4m S\_\mathcal{S}^2 a\_T^4}{\Phi[a\_T]} \right\},\tag{28}$$

$$\frac{\partial \ln \tilde{\xi}}{\partial \ln \dot{m}\_{\text{oT}}} = \alpha\_T \left\{ 1 - n + \frac{4m S\_{\tilde{\xi}}^2 (\alpha\_T - 1) a\_T^3}{\Phi[\alpha\_T]} \right\},\tag{29}$$

$$\frac{\partial \ln \zeta}{\partial \ln \phi} = 2n - 1. \tag{30}$$

Finally, the specific impulse coefficient in Equations (25)–(27) is written in Equation (31):

$$\frac{d\ln I\_{sp}}{d\ln \mathfrak{g}^{\mathfrak{x}}}[\mathfrak{f}] = \frac{\mathfrak{g}}{I\_{sp}[\mathfrak{f}]} \frac{dI\_{sp}}{d\mathfrak{g}^{\mathfrak{x}}}[\mathfrak{f}].\tag{31}$$

The specific impulse reaches its maximum at the optimal mixture ratio *ξopt*. As a consequence, *dIsp*/*dξ*[*ξop<sup>t</sup>*] = 0, and *d* ln *Isp*/*d* ln *ξ*[*ξop<sup>t</sup>*] = 0, which simplifies the previous equations if the engine operates at the optimal mixture ratio.

#### *2.3. Regulation Law of A-SOFT Engines*

While measuring the engine's performances at a given time during its operation, the role of the control law is to determine the commands to be sent to the servo valves to reach the desired thrust and O/F ratio. As seen previously, the thrust and O/F ratio depend on three independent variables which are the two oxidizer mass flow rates ( *m*˙ *oA* and *m*˙ *oT*) and the fuel grain diameter. It is however possible to rewrite Equations (9) and (10) to highlight the role of the total oxidizer mass flow rate and effective swirl number ( *m*˙ *o* and *Se*) which are independent variables as well (Equations (32)–(34)):

$$F[t] = I\_{sp}[\xi] \mathcal{g}\_0 \left\{ \dot{m}\_o[t] + \mathcal{C}\_1 \phi[t] f[\mathcal{S}\_c[t]] \left( \frac{\dot{m}\_o[t]}{\phi^2[t]} \right)^n \right\},\tag{32}$$

$$\xi[t] = \frac{\dot{m}\_o[t]}{\mathcal{C}\_1 \phi[t] f\left[\mathcal{S}\_t[t]\right] \left(\frac{\dot{m}\_o[t]}{\phi^2[t]}\right)^n} \tag{33}$$

$$
\begin{pmatrix} F \\ \int \end{pmatrix} = \begin{pmatrix} F\left[\dot{m}\_o, S\_{\varepsilon\prime}\Phi\right] \\ \int \left[\dot{m}\_{o\prime}, S\_{\varepsilon\prime}\Phi\right] \end{pmatrix}. \tag{34}
$$

From these relations, we can calculate the total derivatives of the thrust and the O/F ratio (Equations (35) and (36)), which can be written in vector form in the Relation (37):

$$\frac{dF}{dt} = \left(\frac{\partial F}{\partial \dot{m}\_o}\right) \frac{d\dot{m}\_o}{dt} + \left(\frac{\partial F}{\partial S\_c}\right) \frac{dS\_c}{dt} + \left(\frac{\partial F}{\partial \phi}\right) \frac{d\phi}{dt} \tag{35}$$

$$\frac{d\xi}{dt} = \left(\frac{\partial\xi}{\partial\dot{m}\_o}\right)\frac{d\dot{m}\_o}{dt} + \left(\frac{\partial\xi}{\partial S\_c}\right)\frac{dS\_\varepsilon}{dt} + \left(\frac{\partial\xi}{\partial\phi}\right)\frac{d\phi}{dt},\tag{36}$$

$$
\begin{pmatrix}
\frac{dF}{dt} \\
\frac{d\xi}{dt}
\end{pmatrix} = \begin{bmatrix}
\frac{\partial F}{\partial \dot{m}\_o} & \frac{\partial F}{\partial S\_c} \\
\frac{\partial \xi}{\partial \dot{m}\_o} & \frac{\partial \xi}{\partial S\_c}
\end{bmatrix} \begin{pmatrix}
\frac{d\dot{m}\_o}{dt} \\
\frac{dS\_c}{dt}
\end{pmatrix} + \begin{pmatrix}
\frac{\partial F}{\partial \phi} \\
\frac{dS\_c}{dt}
\end{pmatrix} \frac{d\phi}{dt}.\tag{37}
$$

In the last relation, the variation related to the port diameter is an environmental change resulting from combustion, whereas the oxidizer mass flow rate and effective swirl number are the variables to be controlled. As a consequence, it is possible to control *dF*/*dt* and *dξ*/*dt* by measuring *dφ*/*dt*, which is the fuel regression rate, and by manipulating *dm*˙ *o*/*dt* and *dSe*/*dt*.

*Aerospace* **2019**, *6*, 65

As for the error propagation analysis, it is possible to calculate the partial derivative coefficients which are given in Equations (38)–(43). It should be noted that these partial derivative coefficients are related with the sensitivity coefficients by the relation: *∂* ln *h ∂*ln*x*= *∂h∂xxh*:

$$\frac{\partial F}{\partial \dot{m}\_o} = g\_0 \frac{(n+\xi)}{\xi} I\_{sp}[\xi] - g\_0 \left(n-1\right) \left(1+\xi\right) \frac{dI\_{sp}}{d\xi}[\xi],\tag{38}$$

$$\frac{\partial F}{\partial S\_{\varepsilon}} = g\_0 \dot{m}\_o \frac{df}{dS\_{\varepsilon}}[S\_{\varepsilon}] \left\{ I\_{\mathbb{S}^p}[\xi] - \xi \left( 1 + \xi \right) \frac{dI\_{\mathbb{S}^p}}{d\xi}[\xi] \right\} \left( f[S\_{\varepsilon}] \xi \right)^{-1} \tag{39}$$

$$\frac{\partial F}{\partial \phi} = g\_0 \left( 2n - 1 \right) \dot{m}\_o \left\{ -I\_{sp} \left[ \xi \right] + \xi \left( 1 + \xi \right) \frac{dI\_{sp}}{d\xi} \left[ \xi \right] \right\} \left( \phi \xi \right)^{-1} \,, \tag{40}$$

$$\frac{\partial \tilde{\xi}}{\partial \dot{m}\_o} = (1 - n) \frac{\tilde{\xi}}{\dot{m}\_o} \, \tag{41}$$

$$\frac{\partial \xi^{\varepsilon}}{\partial S\_{\varepsilon}} = -\frac{\xi^{\varepsilon}}{f[S\_{\varepsilon}]} \frac{df}{dS\_{\varepsilon}}[S\_{\varepsilon}],\tag{42}$$

$$\frac{\partial \hat{\xi}}{\partial \phi} = \left(2n - 1\right) \frac{\hat{\xi}}{\phi}.\tag{43}$$

Equation (37) can be solved so that the derivatives of the controlled variables are given in Equation (44). This system defines the regulation laws that could be used in an A-SOFT engine feedback loop control and which will be analyzed by numerical simulations:

$$\begin{split} \mathcal{I}\left(\frac{d\dot{m}\_{a}}{dt}\right) &= \left(\frac{\frac{\mathcal{I}}{\mathcal{S}\rho\_{1}l\_{sp}[\boldsymbol{\xi}]}\left(1+\boldsymbol{\xi}\right)}{\frac{(1-n)}{\mathcal{S}\rho\_{1}l\_{sp}[\boldsymbol{\xi}]}\dot{m}\_{a}\left(1+\dot{\boldsymbol{\xi}}\right)\frac{d\boldsymbol{\xi}}{d\boldsymbol{\xi}\boldsymbol{\xi}}\left|\mathcal{S}\_{c}\right|}\right)\frac{dF}{dt} \\ &+ \left(\frac{\dot{m}\_{o}\left\{\frac{1}{\mathcal{S}+\boldsymbol{\xi}^{2}}-\frac{1}{I\_{sp}[\boldsymbol{\xi}]}\frac{dI\_{sp}}{d\boldsymbol{\xi}}\left[\boldsymbol{\xi}\right]\right\}}{\frac{dI\_{s}}{d\boldsymbol{\xi}}[\boldsymbol{\xi}]}\left\{-\frac{n+\mathcal{I}}{\mathcal{S}+\boldsymbol{\xi}^{2}}+\frac{(n-1)}{I\_{sp}[\boldsymbol{\xi}]}\frac{dI\_{sp}}{d\boldsymbol{\xi}}\left[\boldsymbol{\xi}\right]\right\}}\right)\frac{d\boldsymbol{\xi}}{dt}+\left(\frac{0}{\frac{(2n-1)}{\mathcal{S}\rho\frac{d\boldsymbol{\xi}}{d\boldsymbol{\xi}}[S\_{c}]}}\right)\frac{d\boldsymbol{\rho}}{dt}.\end{split}\tag{44}$$

#### *2.4. Simplified Modeling of A-SOFT Engines for Simulations*

The main objective of performing numerical simulations is to analyze how an A-SOFT engine would operate autonomously with a feedback loop control system including resistor-based sensors for the fuel regression rate measurement. In particular, the control law and the RBS integration are the major points to be discussed to evaluate the feasibility of such a system.

In order to perform numerical simulations, the A-SOFT engine has been simplified and separated into six distinguishable subsystems as listed below and presented in Figure 2:


**Figure 2.** A-SOFT engine simplification for the simulations.

The oxidizer feed system (1) is not modeled in this study. More precisely, we suppose that the tank is maintained pressurized and that there is no issue regarding the liquid level either.

Ideal servo valves (2) would deliver an oxidizer mass flow rate proportional to the percentage of maximum valve travel (or valve opening). However, some types of valves may have a highly nonlinear relationship between the effective valve flow area (or mass flow rate) and the valve opening [22]. In order to take into account this effect that would be observed in an experimental setup, Equations (45) and (46) provide the relationships between the oxidizer mass flow rate and the corresponding valve opening if a circular port ball valve is used [22]. These relations are both valid for the axial or tangential injectors (A–T) if the same valves are used:

$$\dot{m}\_{o,A-T}[t] = c\_1 \left(1 - e^{-\frac{\upsilon \phi\_{A-T}[t]}{C\_3}}\right),\tag{45}$$

$$v o\_{A-T}[t] = -c\_3 \ln\left(1 - \frac{\dot{m}\_{o,A-T}[t]}{c\_1}\right) + c\_2. \tag{46}$$

*c*1, *c*2, and *c*3 are the valves' parameters and are given in Table 1. This set of parameters has been chosen theoretically so that it induces a maximum oxidizer mass flow rate of 100 <sup>g</sup>·s<sup>−</sup><sup>1</sup> as illustrated in Figure 3. In this example, the valve remains essentially closed until a valve opening *vo* around 20%, then increases rapidly up to *vo* close to 50%. Finally, the mass flow rate is only slightly changed as the valve opening approaches 100%. These parameters should be determined experimentally for a real application.

**Table 1.** Valves' parameters used in this study.


**Figure 3.** Mass flow rate profile in function of the valve opening.

The injection subsystem (3) deals with the axial and tangential oxidizer injectors. Both the total oxidizer mass flow rate and the effective swirl number have been defined previously by Equations (4) and (6) and these definitions are used for the simulations. The total oxidizer mass flow rate is the sum of the axial and tangential oxidizer mass flow rates, and the effective swirl number depends on the geometric swirl number.

The engine (4) simulation is based on three specific calculations which are the fuel regression rate, the combustion chamber and nozzle ejection pressures, and the generated thrust. To consider the most general case possible, the fuel regression rate law is assumed to be slightly different than the one used for the previous sections and is given in Equation (47):

$$\dot{r}[\mathbf{x},t] = a\left(1 + S\_c^2[t]\right)^m \left(\frac{4\dot{m}\_{tot}[\mathbf{x},t]}{\pi\phi[\mathbf{x},t]^2}\right)^n \mathbf{x}^b. \tag{47}$$

The interest of this relation is that it includes the influence of the local mass flux, of the swirl oxidizer injection, and of the axial location along the fuel grain. The coefficients *m* and *n* were chosen identical to the one obtained experimentally and used in previous studies [15]. The *b* coefficient has been chosen arbitrarily with a small value to introduce an explicit dependency on the axial position along the fuel grain. Finally, the *a* coefficient has been calculated so that this regression rate law is consistent with experimental data obtained with an A-SOFT engine [11]. For a real application, these coefficients should be determined carefully and based on experimental results. The total and fuel mass flow rates are given in Equations (48) and (49):

$$
\dot{m}\_{\text{tot}}[\mathbf{x}, t] = \dot{m}\_0[t] + \dot{m}\_f[\mathbf{x}, t], \tag{48}
$$

$$
\dot{m}\_f[\mathbf{x}, t] = \int\_{\mathbf{x'}=0}^{\mathbf{x}} \rho\_f \dot{r}[\mathbf{x'}, t] \pi \phi[\mathbf{x'}, t] d\mathbf{x'}.\tag{49}
$$

The mixture ratio is given by Equation (50):

$$
\zeta[t] = \frac{\dot{m}\_o[t]}{\dot{m}\_f[L, t]}.\tag{50}
$$

The combustion chamber pressure is estimated by using a thermochemical code, in our case a free version of RPA software [23] (Rocket Propulsion Analysis). The pressure is given by Equation (51), and calculated by an iterative resolution at every time step. A pressure error *perror* is then chosen as a convergence parameter. The theoretical isentropic coefficient *γ* is also determined during that process:

$$p\_c[t] = \eta\_c c\_{th}^\*[t] \dot{m}\_{tot} [L\_\prime t] A\_t^{-1}. \tag{51}$$

The nozzle exit pressure is then determined by solving Equation (52), in our case using a simple dichotomy algorithm:

$$\frac{A\_{\varepsilon}}{A\_{\varepsilon}} = \left(\frac{p\_{\varepsilon}[t]}{p\_{\varepsilon}[t]}\right)^{-1/\gamma[t]} \sqrt{\left(\frac{2}{\gamma[t]+1}\right)^{\frac{\gamma[t]+1}{\gamma[t]-1}}} \left(\sqrt{\frac{2}{\gamma[t]-1} \left\{1 - \left(\frac{p\_{\varepsilon}[t]}{p\_{\varepsilon}[t]}\right)^{\frac{\gamma[t]-1}{\gamma[t]}}\right\}}\right)^{-1} . \tag{52}$$

Finally, the thrust is estimated through Equation (53):

$$F[t] = A\_t p\_c[t] \sqrt{\frac{2\gamma[t]^2}{\gamma - 1} \left(\frac{2}{\gamma[t] + 1}\right)^{\frac{\gamma[t] + 1}{\gamma[t] - 1}}} \left\{ 1 - \left(\frac{p\_c[t]}{p\_c[t]}\right)^{\frac{\gamma[t] - 1}{\gamma[t]}} \right\} + \left(p\_c[t] - p\_a\right) A\_t. \tag{53}$$

The feedback loop control is based on getting instantaneous data from the measurement subsystem (5) on the engine's parameters which need to be regulated. In the A-SOFT case, these parameters are the thrust and the mixture ratio. The objective of the simulations is to study the RBS compatibility with an actual control system. To simplify the modeling, we will consider that the thrust and the oxidizer mass flow rate measurements are ideal (Equations (54) and (55)), which means the measured values are the real ones:

$$F\_m[t] = F[t]\_\prime \tag{54}$$

$$
\dot{m}\_{o,m}[t] = \dot{m}\_o[t].\tag{55}
$$

The fuel mass flow rate measurement is realized by embedding RBS in the fuel grain. Figure 4 presents the sensors' principle and the modeling which is adopted in the simulations for a sensor with four resistors.

**Figure 4.** (**a**) resistor-based sensor principle and (**b**) modeling.

While the fuel grain burns and its surface moves, parallel resistors are successively disconnected which modifies the electric voltage at the terminals of the sensor. As a consequence, the measured fuel grain diameter is a constant piece-wise function over time and so are the fuel regression and mass flow rates that are calculated based on it.

As a simplification, we consider that the sensors are ideal: a resistor is disconnected when the fuel surface reaches it. Moreover, we suppose that the resistors are located at the same axial position along the fuel grain length: *dsx* = 0. We consider that if a resistor is connected it delivers a voltage of 1 V and 0 V if not connected. This simplification does not remove any physical sense but avoids the calculation depending on real electric resistance values. The total electric voltage delivered by a sensor is then the sum of resistors' voltages. Applying this rule, the resistors' and sensor's electric voltage are calculated based on Equations (56) and (57). Examples of results describing the sensors' behavior are given in Section 3.2.3.

$$\begin{aligned} \forall \mathbf{s} \in [0, N\_{\mathbf{s}} - 1] \quad &\quad \forall r \in [0, N\_{r} - 1], \\ \quad &if : \phi[\mathbf{x}\_{\mathbf{s}}, t] < \phi\_{0} + 2\mathbf{s}\_{0} + 2r \text{ds} \\ \quad &\quad \text{else} \quad \quad \quad \quad \quad \quad \quad sr[\mathbf{s}, r, t] = 1, \end{aligned} \tag{56}$$

$$\forall \mathbf{s} \in [0, N\_{\mathbf{s}} - 1], \quad \text{ss}[\mathbf{s}, t] = \sum\_{r=0}^{N\_{\mathbf{r}} - 1} sr[\mathbf{s}, r, t]. \tag{57}$$

*Ns* and *Nr* are respectively the number of sensors in the grain and the number of resistors in a single sensor, *xs* is the axial location of a sensor *s*, *ss* and *sr* are respectively the electric signals of a sensor and a resistor.

The measured fuel grain diameter *φ<sup>m</sup>*,*<sup>s</sup>* for a given sensor is then calculated based on the electric signal in Equation (58):

$$\begin{aligned} \forall s \in [0, N\_s - 1], \\ \left[if: ss[s, t] > N\_s - 1 \\ \quad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \phi\_{m, s}[t] = \phi\_0, \\ \quad \qquad \qquad \qquad \qquad \qquad \qquad \phi\_{m, s}[t] = \phi\_0 + 2s\_0 + 2ds \left(N\_r - ss[s, t] - 1\right). \end{aligned} \tag{58}$$

Based on these measurements, the local regression rate can only be evaluated when the sensor's signal changes. By knowing the distance between two resistors and the delay for the fuel surface to move from one to another, Equation (59) provides the calculation method for the measured local regression rate:

$$\begin{aligned} \forall s \in [0, N\_s - 1], \\ \forall i: \phi\_{m,s}[t] = \phi\_{m,s}[t'] \quad & \qquad \dot{r}\_{m,s}[t] = \dot{r}\_{m,s}[t'], \\ \left[ \begin{array}{c} if: \phi\_{m,s}[t] > \phi\_{m,s}[t'] \end{array} \right] \quad & \qquad \dot{r}\_{m,s}[t] = \frac{\phi\_{m,s}[t] - \phi\_{m,s}[t']}{2\left(t - t'\right)}. \end{aligned} \tag{59}$$

In this relation, *t* and *t* are two successive times when the sensor's signal has changed and *t* > *t* . The fuel regression rate is constant when the measured diameter is not changed. As a consequence, the fuel regression rate can only be evaluated once the first resistor has been reached by the fuel surface.

The local and total fuel mass flow rates are then evaluated by Equations (60) and (61):

$$\forall \mathbf{s} \in [0, N\_{\mathbf{s}} - 1], \ \dot{m}\_{f, m, s}[t] = \pi \rho\_f \dot{r}\_{m, s}[t] \phi\_{m, s}[t] \frac{L}{N\_{\mathbf{s}}},\tag{60}$$

$$
\dot{m}\_{f,m}[t] = \sum\_{s=0}^{N\_s - 1} \dot{m}\_{f,m,s}[t]. \tag{61}
$$

Finally, the measured mixture ratio based on RBS is written in Equation (62):

$$
\zeta\_{m}^{\chi}[t] = \frac{\dot{m}\_{o,m}[t]}{\dot{m}\_{f,m}[t]}.\tag{62}
$$

The last subsystem to be modeled is the feedback control system (6). It can be divided into two parts: the first part concerns the use of a PID (Proportional Integral Derivative) controller and the second part is the integration of the control law. As demonstrated in the previous section, Equation (44) allows for calculating the derivatives of the oxidizer mass flow rate and effective swirl number in function of the thrust and mixture ratio derivatives and the regression rate. In the case of the simulation, the feedback control system must follow the thrust and mixture ratio profiles which are given. The control law is used to determine the targeted total oxidizer mass flow rate and effective swirl number to satisfy the control command. As a consequence, the required terms appearing in the control law are given by Equations (63)–(65):

$$\frac{dF}{dt}[t] = \frac{F\_{\text{common}}[t+dt] - F\_{\text{fl}}[t]}{dt} \, \tag{63}$$

$$\frac{d\underline{\xi}}{dt}[t] = \frac{\underline{\xi}\_{command}[t+dt] - \underline{\xi}\_{m}[t]}{dt},\tag{64}$$

$$\frac{d\Phi}{dt}[t] = 2\dot{r}\_m[t].\tag{65}$$

In the last relation, *<sup>r</sup>*˙*m* represents the measured regression rate and is averaged based on the local values provided by the sensors. Based on these calculations, we can evaluate the targeted oxidizer mass flow rate and swirl number to satisfy the commands through Equations (66) and (67):

*m* ˙ *<sup>o</sup>*,*target*[*<sup>t</sup>* + *dt*] = *dm*˙ *o dt* [*t*]*dt* + *m*˙ *<sup>o</sup>*,*<sup>m</sup>*[*t*], (66)

$$S\_{\varepsilon, \text{target}}[t + dt] = \frac{d S\_{\varepsilon}}{dt}[t] dt + S\_{\varepsilon, m}[t]. \tag{67}$$

The targeted axial and tangential oxidizer mass flow rates are then determined based on Equations (4) and (6), and the targeted valves opening by Equation (46). The PID error is then calculated (Equation (68)):

$$
\sigma\_{\text{grid}}[t] = \upsilon \alpha\_{A-T, \text{tar}\chi et}[t+dt] - \upsilon \alpha\_{A-T}[t]. \tag{68}
$$

The PID correction to be applied to the controlled valve is then given by Equation (69):

$$
\delta \nu o\_{A-T}[t+dt] = k\_P e\_{p\dot{u}}[t] + k\_I \int\_0^t e\_{p\dot{u}}[t']dt' + k\_D \frac{d e\_{p\dot{u}}}{dt}[t]. \tag{69}
$$

Finally, the commands to be sent to the servo valves are given as follows (Equation (70)):

$$
\omega \circ\_{A-T} [t + dt] = \upsilon \circ\_{A-T} [t] + \delta \upsilon \circ\_{A-T} [t + dt]. \tag{70}
$$

Regarding the numerical parameters, time derivatives are calculated based on an Euler explicit scheme and integrals are estimated through the trapezoids' method. In order to consider the latency that occurs in real experimental setup regarding the valves response times, pipeline delays, controller calculation time, etc., an artificial regulation frequency *freg* has been added for the numerical simulations. This frequency has been fixed to 2 Hz for all the simulations, meaning that the feedback control system sends new orders to the valves every 0.5 s. This value has been chosen to be realistic but is not based on experimental results or calculations. It is an additional parameter and its influence could be investigated in a future work. The numerical parameters that have been used for the simulations are summarized in Table 2.

**Table 2.** Numerical parameters used for the simulations.


## **3. Results**

#### *3.1. Error Propagation Analysis*

The A-SOFT concept is based on an innovative idea which is to use the swirl influence on the fuel regression rate to control the thrust and the O/F ratio at the same time but independently. However, other concepts exist to achieve this objective [7–9]. Among them, the AOA (Aft-end Oxidizer Addition) hybrid engine [8,9] is somehow similar to the A-SOFT engine. Unlike A-SOFT that uses dual front injections, the AOA engine combines a front injection with an aft-end oxidizer addition in order to regulate the thrust and the mixture ratio (Figure 5).

**Figure 5.** AOA (Aft-end Oxidizer Addition) concept with front and aft-end oxidizer injections.

Despite the concepts being different, both use two servo valves to regulate the thrust and the O/F ratio and it seems logical to compare their respective interests. To do so, the relative error analysis has also been conducted for the AOA engines and will be compared to the one conducted for the A-SOFT engines. Because the mathematical development is very similar to the one presented in the previous section and does not bring any additional information, materials and methods related to the AOA engines are given in Appendix A. An important point to note is the difference regarding the *αT* parameter: in the case of A-SOFT, it represents the ratio between the tangential and the total oxidizer mass flow rates (Equation (23)), while it represents the ratio between the front and the total oxidizer mass flow rates in the case of AOA (Equation (A15)).

#### 3.1.1. Optimal Mixture Ratio

To evaluate the error propagation, the study case is chosen based on experimental results [18]. We will consider an optimal mixture ratio configuration with corresponding parameters given below:


Figure 6 presents the normalized thrust error for three given *<sup>e</sup>φ*/*eMFR* ratios. As we can first notice, the *<sup>e</sup>φ*/*eMFR* ratio plays a major role in the error propagation. The thrust normalized error will increase with the ratio which means that, for a given oxidizer mass flow rate measurement error, the increase of the diameter measurement error will increase the thrust error as well. It is also clear that the geometric swirl number *Sg* has an influence on the A-SOFT errors. For *Sg* = 0, the minimal thrust errors are obtained for: *αT* ≈ 0.35 both for A-SOFT and AOA. If *Sg* = 0, the error profiles are symmetric around *αT* = 0.5 for A-SOFT since there is no swirl effect. We can notice that A-SOFT and AOA errors are very similar for *αT* ≥ 0.5 and become equal if *αT* = 1. In that case, the A-SOFT engine uses only the tangential injection and the AOA engine uses only the front injection (which is also tangential), the engines are consequently strictly identical. Finally, the AOA engine generates less thrust propagation errors than the A-SOFT engine for low values of *αT*. However, this result must be contrasted since operating in these conditions would imply that most of the oxidizer mass flow rate would come from the aft-end injection and may result in combustion troubles that are not considered at all in this analysis.

Figure 7 presents the normalized O/F ratio error for the three same *<sup>e</sup>φ*/*eMFR* ratios. Some of the previous conclusions are naturally shared with the previous ones: the A-SOFT symmetric profile when *Sg* = 0, the equal errors for A-SOFT and AOA when *αT* = 1 and the *<sup>e</sup>φ*/*eMFR* influence are also true regarding the mixture ratio. We can notice that the increase of *<sup>e</sup>φ*/*eMFR* ratio from 0 to 5 can multiply the relative errors on O/F ratio by a factor 6 while this factor was limited to a value lower than 1.5 in the case of thrust error propagation. Naturally, the mixture ratio errors are much more sensitive to the precision of the fuel diameter measurements since there is a direct relation between O/F ratio and fuel grain diameter. In contrast to the thrust error, the mixture ratio error for A-SOFT is significantly lower than for AOA (especially when *αT* ≤ 0.5). The main reason to explain this effect is that A-SOFT engines are more flexible and efficient to regulate the mixture ratio. Indeed, both oxidizer injections influence the fuel regression rate in A-SOFT while only the front injection can do it in AOA. At this point, the overall conclusions comparing A-SOFT and AOA are not evident, but Figure 8 which describes the total errors will help to conclude.

**Figure 6.** Normalized thrust errors.

**Figure 7.** Normalized mixture ratio errors.

Once again, some of the previous conclusions regarding the thrust and the mixture ratio are shared with the total error propagation (symmetry, case when *αT* is high, influence of *<sup>e</sup>φ*/*eMFR* ratio). From the total error point of view, it can be said that the AOA engine's error is generally larger than the A-SOFT one, especially for low *αT* (mainly due to the mixture ratio error) even though these differences tend to decrease as the diameter error increases. Based on these results, the A-SOFT concept seems more interesting regarding the total error propagation and the following discussions will concern A-SOFT engines only.

As visible in Figure 8, the A-SOFT total error increases as *Sg* increases, regardless of the *<sup>e</sup>φ*/*eMFR* ratio. Figure 9 focuses on this effect and presents the results for *Sg* ranging from 0 to 50 (*Sg* = 0, 2, 4, 6, 8, 10, 12, 15, 20, 30, 40, 50) at a given *<sup>e</sup>φ*/*eMFR* ratio. The case where *Sg* = 0 leads to the lowest errors but is purely theoretical and presents no interest since there would be no swirl control which is the basic principle of A-SOFT engines. For *αT* lower than 0.1 and higher than 0.9, the geometric swirl number has almost no influence on the results; however, these *αT* values generate the highest errors for low *Sg* and are the extreme operating conditions of the engine. The first local minimum at *αT* = 0.5 for *Sg* = 0 is constantly shifted to lower *αT* values and the error continues to increase with *Sg*. At the same time, a second local minimum appears around *αT* = 0.68 if *Sg* ≥ 6. Unlike the first one, this second local minimum error tends to an asymptotic limit and stops increasing even for high geometric swirl numbers. Moreover, as *Sg* increases, a local maximum appears for which the induced error can become larger than the extremity errors if *Sg* exceeds 30. These results tend to indicate that, for a given geometric swirl number, the total error can be minimized by operating close to optimal *αT* values, but also that some other values of *αT* may lead to significant errors and should be avoided. As *αT* varies during the engine operation, additional analysis is required to design and choose the optimal geometric swirl number.

**Figure 8.** Normalized total errors.

**Figure 9.** Influence of the geometric swirl number on the normalized total errors in A-SOFT engines for *Sg* values in [0,2,4,6,8,10,12,15,20,30,40,50].

#### 3.1.2. Optimal Design of the Geometric Swirl Number

As we have seen, the geometric swirl number has a significant influence on the error propagation in an A-SOFT engine. For a given *Sg*, there are optimal operating *αT* to minimize the errors but other values of *αT* could drastically increase these errors. In practice, the fuel regression rate is influenced by the effective swirl number provided by the oxidizer's injections and these parameters are related by Equation (71):

$$S\_{\mathfrak{e}}[t] = S\_{\mathfrak{K}} a\_T^2[t]. \tag{71}$$

For a given mission, the effective swirl number will be modified to fit as much as possible the thrust and the O/F commands while the grain is burning and the fuel port diameter is modified. The range of effective swirl number will obviously depend on the mission profiles, on the fuel and oxidizer combination, on the regression rate coefficients, and on other parameters like the feedback loop control system. As a consequence, it is not possible to chose a geometric swirl number which would be optimal unless the mission is known.

That being said, we could propose a method for designing the optimal geometric swirl number for a given mission, based on the error propagation results. As a study case, let us consider a lab-scale engine operating during one minute with the regression rate coefficients given earlier in the paper, and with a constant target thrust of 250 N. This should lead to an operating effective swirl number ranging between 10 and 50. Moreover, we will assume that *αT* increases linearly with time, which is the simplest evolution possible. Finally, we consider the case of an optimal mixture ratio:


Once the effective swirl number range is known, there is an infinite number of combinations of possible geometric swirl number and operating *αT*. By definition, *Sg* ≥ *Se*,*max* since *αT* never exceeds 1 (Equation (71)). Table 3 summarizes a few examples:


**Table 3.** Examples of possible geometric swirl number and *αT* combinations for a given effective swirl number range.

In order to find a potential optimal geometric swirl number regarding the total error propagation, two natural criteria could be defined: the first one could be to ensure that the error will never exceed a certain level at a given time (C1), and the second one could be to minimize the error over the entire operation duration (C2). If an optimal geometric swirl number *Sg*,*op<sup>t</sup>* does exist, these criteria could be expressed by:

• C1: ∀*<sup>α</sup>T*[*t*] : *etotal*[*<sup>α</sup>T*[*t*], *Sg*,*op<sup>t</sup>*] ≤ *etotal*,*max*,

$$\bullet \quad \mathsf{C2} \colon \forall \mathcal{S}\_{\mathcal{S}} : \int\_{\mathcal{I}} \overline{\varepsilon\_{total}} [\mathcal{a}\_T[t], \mathcal{S}\_{\mathcal{S}, \mathcal{opt}}] \leq \int\_{\mathcal{I}} \overline{\varepsilon\_{total}} [\mathcal{a}\_T[t], \mathcal{S}\_{\mathcal{S}}] \,.$$

As an example, Figure 10 shows that using a geometric swirl number of 40 could lead to a total error exceeding the upper limit if the engine has to operate with *αT* close to 0.3. Figure 11 illustrates the calculation of the integrated errors. For a given range of effective swirl number (here from 10 to 50), modifying the geometric swirl number will change the integrated error. The blue and red areas respectively represent the integrated errors for *Sg* = 50 and 80. The integrated error between two values of effective swirl numbers, denoted by *Se*,<sup>1</sup> and *Se*,2, is written in Equation (72):

$$
\overline{E\_{\text{total}}}[\mathcal{S}\_{t,1}, \mathcal{S}\_{t,2}, \mathcal{S}\_{\mathcal{S}'}a\_T[t]] = \int\_{t=0}^{t=t\_f} \overline{\mathcal{E}\_{\text{total}}}[a\_T[t], \mathcal{S}\_{\mathcal{S}}]dt.\tag{72}
$$

**Figure 10.** Illustration of the C1 criterion.

**Figure 11.** Illustration of the C2 criterion.

If the effective swirl number range is fixed, the *αT* evolution depends on *Sg* and can be written by Equation (73) (where *tf* is the firing duration).

$$a\_T[t] = \sqrt{\frac{S\_{\varepsilon,1}}{S\_{\mathcal{S}}}} + \frac{\sqrt{\frac{S\_{\varepsilon,2}}{S\_{\mathcal{S}}}} - \sqrt{\frac{S\_{\varepsilon,1}}{S\_{\mathcal{S}}}}}{t\_f} t. \tag{73}$$

Figure 12 provides the evolution of the integrated error by varying the geometric swirl number in the configuration described earlier. As it can be seen, an optimal geometric swirl number does exist in that case and is around 62.5.

**Figure 12.** Integrated error evolution depending on the geometric swirl number.

*Aerospace* **2019**, *6*, 65

This optimal value minimizes the integral error over the firing test duration, but, additionally, its maximum instantaneous error is lower than the ones for *Sg* = 50 or 100 (Figure 13). Consequently, it seems that choosing *Sg* = 62.5 for this firing configuration would be an interesting choice regarding the error's criteria.

As a partial conclusion, we proposed a method based on the error propagation analysis to choose an optimal geometric swirl number choice regarding the defined criteria. For a given engine configuration, an optimal *Sg* was found by minimizing the integrated total error over the firing duration.

**Figure 13.** Instantaneous total error evolution.

#### 3.1.3. Optimal Design of the Resistor-Based Regression Rate Sensors

In the previous sections, the ratio *<sup>e</sup>φ*/*eMFR* was arbitrarily fixed to study different cases. In practice, this parameter depends on the precision of the fuel regression and oxidizer mass flow rates measurement sensors. In this study, we consider using RBS to control an A-SOFT through a feedback loop system (Figures 2 and 4). Such sensors have various possible designs and one of the most important parameter is the radial distance between two resistors. Each time the fuel grain surface reaches a resistor, this resistor is burnt and the electric voltage provided by the sensor is modified. As a consequence, the fuel grain diameter measured by the sensor is a piece-wise constant function. In this discussion, we will consider that the sensors are ideal: the resistors are burnt as soon as the fuel surface reaches them. Moreover, we assume the resistors are axially aligned, and the initial distance grain surface–first resistor is equal to the distance between two resistors:


In this section, the error propagation analysis will be used to determine the design of RBS to satisfy given conditions. As a starting point, the error due to the oxidizer mass flow measurement must be assumed. Depending on the type of measurement system and sensor, this error generally ranges from 0.2% to 5%. As an intermediate possibility, we will consider an error of 2%. The next errors to determine are the maximum acceptable errors regarding the thrust and the mixture ratio. In this example, we will consider that these errors are equal to 2% as well:


These values lead to the calculation of the maximum tolerable normalized errors: *eF*,*max* = 1, *<sup>e</sup>ξ*,*max* = 1 and: *etotal*,*max* = √2. Based on the previous relations (Equations (20)–(22)), we can define a maximum normalized error on the diameter by Equation (74):

$$\overline{v\_{\phi,\max}} = \sqrt{\frac{\overline{v\_{\text{total},\max}}^2 - \left(\frac{\partial \ln F}{\partial \ln \overline{m}\_{\text{oL}}}\right)^2 - \left(\frac{\partial \ln F}{\partial \ln \overline{m}\_{\text{v}T}}\right)^2 - \left(\frac{\partial \ln \overline{\xi}}{\partial \ln \overline{m}\_{\text{s}A}}\right)^2 - \left(\frac{\partial \ln \overline{\xi}}{\partial \ln \overline{m}\_{\text{s}T}}\right)^2}{\left(\frac{\partial \ln F}{\partial \ln \overline{\xi}}\right)^2 + \left(\frac{\partial \ln \overline{\xi}}{\partial \ln \overline{\phi}}\right)^2}}. \tag{74}$$

According to this relation, there may be no solution if the maximum total error *etotal*,*max* is too small since it would conduct to a negative term in the square root. Figure 14 presents the results for several *Sg* values and for the optimal mixture ratio case.

**Figure 14.** Maximum normalized diameter error for *Sg* values in [20,30,40,50,60,70,80,100].

As it can be seen, there is a local minimum if *αT* varies around 0.2 or 0.3 depending on the geometric swirl number. The value of *<sup>e</sup>φ*,*max* tends to be reduced while the geometric swirl number increases. In other words, the higher *Sg* is, the more the constraint on the fuel regression rate sensors is severe and the more the relative error needs to be reduced. As an example, let us consider a geometric swirl number of 60: the maximum normalized error on the diameter would be around 3 if the engine operates at *αT* close to 0.21 according to the previous result. In the case of ideal RBS, the maximum error on the diameter could be evaluated by Equation (75):

$$
\varepsilon\_{\Phi} = \frac{\delta\phi}{\phi} = \frac{2ds}{\phi}.\tag{75}
$$

Since the fuel grain diameter will naturally increase during the combustion, the diameter relative error will decrease with time. Consequently, to estimate the required *ds* parameter, the initial fuel grain diameter should be chosen to calculate the maximum diameter relative error possible. In the current case, if we choose an initial fuel port diameter *φini* = 40 mm, this would lead to the following results:


In this configuration and by using the propagation error analysis, we could conclude that using a RBS radial distance lower than 1.2 mm between resistors would limit the thrust and mixture ratio relative errors to less than 2%. Naturally, increasing the tolerance regarding these errors would lead to increase the *ds* parameter, which may be required by practical manufacturing issues or other external factors.

#### 3.1.4. Non-Optimal Mixture Ratio

A-SOFT hybrid engines ideally operate close to the optimal mixture ratio to reach the best performance possible. However, while operating, the O/F ratio can become higher or smaller than the optimal value due to the feedback loop control efficiency/latency or due to other limiting factors like a long burning duration resulting in a too large fuel grain diameter. Additional analyses have been carried out for two other operating conditions, one dealing with an oxidizer-rich combustion and one dealing with a fuel-rich configuration. The results regarding the total error and the geometric swirl number influence are very similar to the optimal mixture ratio case and are reported in Appendix B. These results reveal that operating at non-optimal mixture ratio does not change the tendencies and conclusions which were given in the optimal case. Consequently, the O/F ratio variations from optimal to non-optimal values and vice versa under feedback control will have no significant impact on the errors' propagation.

#### *3.2. Numerical Simulations*

#### 3.2.1. Study Case

Several numerical simulations have been performed and the results are given in this section. Table 4 provides the parameters that have been used and are common to all the simulations. Numerical parameters can be seen in Table 2. Thermochemical data were obtained from RPA for an engine operating with gaseous oxygen as oxidizer and high density polyethylene as fuel.


**Table 4.** Numerical simulations' parameters.

In all the simulations, the servo valves are open at 30% at the initial instant. Figures 15 and 16 respectively show the natural evolution of the thrust and the mixture ratio of the engine without any type of feedback control and illustrate the performances' shift that was mentioned earlier in the paper. In this configuration, the engine was not designed to operate around optimal conditions regarding the thrust and the mixture ratio.

**Figure 15.** Thrust profile without feedback control.

**Figure 16.** Mixture ratio profile without feedback control.

#### 3.2.2. Feedback Control Law Efficiency

In the following discussion, three thrust profiles have been considered: constant, piece-wise constant, and linearly decreasing thrust. In order to operate at the optimal mixture ratio, the O/F command has been fixed constant for all cases and equal to the optimal value. In order to evaluate the efficiency of the regulation law that was developed in the previous section (complete regulation law), a second regulation law based on a very simple proportional rule has been tested (downgraded regulation law). The following list sums up the performed simulations:


The second regulation law that has been tested simply considers that the thrust is proportional to the oxidizer mass flow rate (downgraded regulation law). Rather than using Equations (66) and (67), the target oxidizer mass flow rate is calculated based on Equation (76) and the effective swirl number is then deduced to fit to the targeted mixture ratio:

$$
\dot{m}\_{o,tar\_{\mathcal{X}}\text{ct}}[t+dt] = \dot{m}\_o[t] \frac{F\_{common}[t+dt]}{F\_m[t]}.\tag{76}
$$

The resistor-based sensors parameters are given in Table 5. The sensors are distributed uniformly along the fuel grain.

**Table 5.** RBS (Resistor-Based Sensors) parameters for the simulations.


Regarding the PID controller, a proportional controller has been evaluated in this study and was sufficient to obtain interesting results. The PID parameters were obtained experimentally based on the simulations' results and are given in Table 6. This set of parameters simplifies the PID control law (Equation (69)), which is consequently written by Equation (77).

**Table 6.** PID (Proportional Integral Derivative) parameters for the simulations.


Figures 17–20 respectively present the thrust, the mixture ratio, the effective swirl number, and the total oxidizer mass flow rate results of the simulation case 1 with the complete regulation law. As it can be seen, the initial thrust is higher than the thrust command. The initial value of the valve opening was voluntarily fixed arbitrarily to evaluate the capacity of the feedback law to rapidly match the command. The instantaneous values of the valves opening are calculated autonomously by the algorithm after the first iteration and according to the regulation law. The thrust reaches the command after 5 s of combustion and the feedback control is able to maintain the thrust level at 260 N ± 1 N until 55 s. After 55 s, it can be seen that the thrust starts to decrease, which will be explained by other results given later in the paper.

The estimation of the fuel mass flow and regression rates can only be done once the first resistor of the RBS has been reached by the fuel surface. As a consequence, the feedback loop control starts

when the RBS provide the required measurements, which can take several seconds depending on the sensors' parameters and the fuel regression rate. With this set of RBS parameters and these test conditions, this delay is around 1.2 s. The initial mixture ratio is equal to 2.9 and rapidly changes to 2.4 after 3.5 s of combustion. After 4.5 s, the O/F ratio reaches 2.35 and is maintained between this value and the target 2.3 until 50 s. Similarly to the thrust, the engine cannot reach the target at the end of the test. Once the RBS are providing fuel regression rate measurements, the difference between the real mixture ratio and the measured one is kept lower than 0.04 which corresponds to 1.74% of the command value until 50 s. At this time, a significant part of the fuel grain has burnt and the sensors no longer provide measurements, which leads to an increasing error on the mixture ratio evaluation. Finally, we can notice that the measured O/F ratio is always slightly higher than the command and varies from 2.3 to 2.34, which leads to an oxidizer-rich combustion by the feedback control.

After a rapid regulation at the beginning of the test, the effective swirl number increases linearly until reaching 50. In this configuration, the geometric swirl number was fixed at 50, which means that *Se* cannot be higher than this value. Around 55 s, the swirl number is limited to 50 and the feedback law is no longer able to follow the thrust and the O/F ratio commands due to the operating conditions and the fuel grain combustion after more than 50 s. It explains why the target values are not reached at the end of the test. The oxidizer mass flow rate is rapidly reduced after the ignition because the initial valve opening value was providing a higher thrust than the target. After this phase, the oxidizer mass flow rate is maintained approximately constant until 55 s, when the conditions no longer allow to follow the commands.

**Figure 17.** Thrust profile, feedback with complete regulation law, case 1.

**Figure 18.** Mixture ratio profile, feedback with complete regulation law, case 1.

**Figure 19.** Effective swirl number profile, feedback with complete regulation law, case 1.

**Figure 20.** Oxidizer mass flow rate profile, feedback with complete regulation law, case 1.

Figures 21–24 present the results of the simulation case 2 with the complete regulation law. The main conclusions are similar to the previous case and the engine is able to follow the thrust and the O/F ratio commands. The simulations configuration are identical between case 1 and case 2 until 30 s and the thrust and the O/F ratio evolution are the same. At 30 s, the thrust command suddenly increases from 260 to 280 N and the engine's thrust reaches the new target in approximately 2 s. From 32 s to 48.5 s, the thrust level varies around 280 N ± 1 N. After 48.5 s, the thrust starts to decrease as it is the case in the previous simulation but earlier than it (it was after 55 s).

**Figure 21.** Thrust profile, feedback with complete regulation law, case 2.

**Figure 22.** Mixture ratio profile, feedback with complete regulation law, case 2.

We can notice small oscillations of the O/F ratio after the change of the thrust level at 30 s, the regulation requiring some time to fit the new conditions. The mixture ratio slightly oscillates around 2.3 ± 0.1 until 50 s when it starts to increase significantly. Because the second level of thrust target is higher than the first target, the total oxidizer mass flow rate needs to be increased and the fuel grain is burnt quicker than in case 1. As a consequence, the thrust and the O/F ratio start to move away from the targeted values earlier than in the previous case.

**Figure 23.** Effective swirl number profile, feedback with complete regulation law, case 2.

**Figure 24.** Oxidizer mass flow rate profile, feedback with complete regulation law, case 2.

Figures 25–28 present the results of the simulation case 2 with the downgraded regulation law. The results are similar than the previous ones, even if the regulation law is strongly downgraded in terms of precision. However, we can observe more oscillations regarding the thrust and the mixture ratio. The thrust reaches the command after 5 s of combustion and the feedback control maintains the thrust level at 260 N ± 4 N until 30 s with small oscillations. At 30 s, the engine's thrust reaches the new target in approximately 2 s. From 32 s to 50 s, the thrust level varies around 280 N ± 5 N. The thrust starts to decrease after 50 s.

**Figure 25.** Thrust profile, feedback with downgraded regulation law, case 2.

**Figure 26.** Mixture ratio profile, feedback with downgraded regulation law, case 2.

The initial mixture ratio is equal to 2.9 and rapidly changes to 2.4 after 3.5 s of combustion. After that time, the O/F ratio oscillates between 2.27 and 2.48 until 30 s, the amplitude of the oscillations decreasing with time. After 30 s, the oscillations' amplitude increases again due to the change of the thrust command, and the mixture ratio varies between 2.18 and 2.44 until 50 s. Results regarding cases 1 and 3 for the downgraded regulation law are given in Appendix C.

**Figure 27.** Effective swirl number profile, feedback with downgraded regulation law, case 2.

**Figure 28.** Oxidizer mass flow rate profile, feedback with downgraded regulation law, case 2.

Figures 29–32 present the results of the simulation case 3 with the complete regulation law. The thrust reaches the command after 5 s and is maintained close to the target with a difference lower than 1 N during all the test. The mixture ratio varies between 2.3 and 2.36 from 5 s to 50 s.

**Figure 29.** Thrust profile, feedback with complete regulation law, case 3.

**Figure 30.** Mixture ratio profile, feedback with complete regulation law, case 3.

**Figure 31.** Effective swirl number profile, feedback with complete regulation law, case 3.

**Figure 32.** Oxidizer mass flow rate profile, feedback with complete regulation law, case 3.

#### 3.2.3. Resistor-Based Sensors' Influence

Resistor-based sensors play a major role in the engine feedback control, Figure 33 provides the sensors' signals in case 2 with the complete regulation law. The signals are piece-wise functions and the regression and fuel mass flow rates are re-evaluated every time the signals are modified.

**Figure 33.** RBS signals, feedback with complete regulation law, case 2.

The previous simulations were performed with sensors equipped with a large number of resistors (*Nr* = 30), which may be difficult to manufacture. In order to study another configuration, this number was changed to 10 and the distance between resistors was increased from 1 mm to 3 mm (Table 7); the number of sensors was not changed.

**Table 7.** New RBS parameters for the simulations.


Figures 34–36 present the results of the simulation case 2 with the complete regulation law and the new RBS parameters.

The thrust reaches the command after 7.5 s and then slightly oscillates around 260 N ± 2 N until 30 s. At that time, the target level is changed and 2 s are needed to fit to the new command value. From 32 s to 48 s, the thrust oscillates around 280 N ± 3 N and then starts to decrease at the end of the simulation. The initial delay required to reach the first resistor of the sensors and to start the feedback control was changed from 1.2 s to 3.8 s with the new set of RBS parameters. This delay is multiplied by a factor close to 3 which approximately corresponds to the modification from 1 mm to 3 mm of the initial distance between the first resistor and the fuel grain surface. The mixture ratio is never completely stabilized during the test and varies within 2.1 and 2.6 after 5 s of combustion. The maximum error between the real O/F ratio and the measured value is close to 0.1 corresponding to 4.3% error before 30 s, but can reach up to 0.3 after 30 s, which corresponds to 13% of the target value. The error is multiplied by a factor 3 compared to the previous simulation which corresponds to the change of the distance between resistors. Although the thrust remains relatively close to the target, it is clear that the new set of sensors' parameters generated more oscillations and that the engine is more difficult to control. It is particularly true regarding the mixture ratio which directly depends on the RBS measurements.

**Figure 34.** Thrust profile, feedback with complete regulation law, case 2, *ds* = 3 mm.

**Figure 35.** Mixture ratio profile, feedback with complete regulation law, case 2, *ds* = 3 mm.

**Figure 36.** RBS signal, feedback with complete regulation law, case 2, *ds* = 3 mm.

Finally, Figures 37 and 38 present the results of the simulation case 2 with the downgraded regulation law and the new RBS parameters.

**Figure 37.** Thrust profile, feedback with downgraded regulation law, case 2, *ds* = 3 mm.

**Figure 38.** Mixture ratio profile, feedback with downgraded regulation law, case 2, *ds* = 3 mm.

The thrust reaches the command after 7.5 s and then oscillates significantly between 255 and 275 N until 30 s. The change to the second thrust level seems to stabilize the engine which operates around 280 N ± 2 N until 55 s. The mixture ratio is also oscillating strongly and varies between 2.1 and 2.7 before 30 s. It becomes more stable in the second phase of the test. These simulations clearly demonstrate that the complete regulation law provides significantly better results than the

downgraded one, especially if the RBS precision is not very high, and also illustrated the importance of the sensors design for a proper feedback control.

## **4. Discussion**

#### *4.1. Error Propagation Analysis*

The error propagation analysis has been performed in order to study the influence and importance of the measurements errors on the thrust and mixture ratio control. Two similar concepts (called A-SOFT and AOA) based on a dual oxidizer injection have been investigated and compared in terms of thrust, mixture ratio, and total relative errors. This analysis allows for calculating the errors propagation produced by measurements errors and could be used to help design future A-SOFT hybrid engines for real missions. Even though the AOA engine seems more interesting for limiting the thrust error propagation, we have seen that this theoretical result was not guaranteed. Indeed, it would suppose that the AOA engine has to operate at low values of *αT* which would induce other operating troubles that are not taken into account in this study. Despite the fact that the benefit of using an A-SOFT engine rather than an AOA engine in terms of thrust error propagation could be open to debate, there is no ambiguity regarding the mixture ratio and the total relative errors. Because both of the two oxidizer injections have an influence on the fuel regression rate in the A-SOFT engines (only one in the AOA engines), this type of motor is more flexible and efficient to control the mixture ratio compared to the AOA engines. Consequently, the A-SOFT motors significantly reduce the total error propagation and are more interesting than the AOA motors in this regard. The principal assumptions to conduct this study were to neglect the combustion chamber pressure influence on the specific impulse and to suppose a constant combustion efficiency. These hypotheses lead to a specific impulse that is a function of the mixture ratio only. Since the A-SOFT and AOA engines use a swirled oxidizer injection, the combustion efficiency should keep a high value during the operating duration. As a consequence, assuming a constant combustion efficiency is not a very strong hypothesis and should not affect the described results. The error propagation analysis has been conducted and applied to propose two design methods, described as follows:

(1) The first design method deals with the choice of an optimal geometric swirl number for the A-SOFT engine. As we have seen thanks to the analysis, an augmentation of the geometric swirl number directly increases the total relative error. As a consequence, we could conclude that the lowest *Sg* possible may be the most interesting choice. However, the error propagation increase due to the geometric swirl number augmentation is more or less significant depending on the engine operation condition (*<sup>α</sup>T*). Consequently, the design of an optimal geometric swirl number is not trivial and depends on other factors. Even though different choices could be made, we considered two optimal criteria in this study: the first one was to never exceed a maximal error limit and the second one was to minimize the total error for the entire firing test duration. We applied these criteria on a specific case where the effective swirl number varies between 10 and 50, and the analysis conducted to an optimal geometric swirl number of 62.5, different from 50 that would have been the minimal *Sg* possible. We demonstrated that the error propagation analysis could be used to define the geometric swirl number of an A-SOFT engine if the mission profile is known.

(2) The objective of the second method arising from the error analysis was to help define the resistor-based sensors. In particular, the radial distance between resistors is the most crucial parameter for the mixture ratio evaluation and must be chosen carefully. If the maximum total error which is tolerable and the oxidizer mass flow rate measurement error are known, the propagation analysis can be used to define the requirements for the radial distance between resistors, as demonstrated in this paper. Further analysis may be required to design other parameters such as the axial distance between resistors, the number of sensors, and their location along the grain. However, in order to study all the parameters, the theoretical analysis should be coupled with dedicated numerical simulations and experimental firing tests. Although the specific case of resistor-based sensors has been considered in this study, the method is interesting since it could be used to analyze other type of fuel regression rate measurements.

In this paper, we proposed and conducted independently two methodologies regarding the choice of the geometric swirl number and the design of resistor-based sensors. However, these parameters are not independent and should be designed in parallel. In practice, a few iterations of the described methods may be required to design the ideal swirl injector geometry and the corresponding fuel regression sensors in order to minimize the error propagation in an A-SOFT engine.

#### *4.2. Numerical Simulations*

The main objective of the numerical simulations was to evaluate the regulation law derived from fundamental A-SOFT equations and to apply it with the use of resistor-based sensors for the fuel regression rate measurement. The results described in the article showed that the A-SOFT engine feedback control could be performed in order to manage the thrust and mixture ratio independently and at the same time, and that using RBS was a promising solution for the mixture ratio evaluation. Since the sensors' principle is simple, it is possible to model their behavior and to study their implementation using simplified simulations. Although the simulations revealed interesting features of the feedback control, several assumptions have been made and deserve some discussion:

(1) As in the error propagation analysis, a constant combustion efficiency has been assumed to perform the simulations and particularly to calculate the engine performances. Because of the swirled oxidizer injection, this hypothesis is not very strong and the combustion efficiency should keep a high value in a real configuration. However, it would be possible to remove this assumption by using, for example, an empirical law giving the efficiency in function of the swirl number and the combustion chamber geometry. Detailed experimental analysis would then be required to obtain such an empirical relation.

(2) In our model, the regulation by the PID controller starts a few seconds after ignition since it needs the data from the resistor-based sensors and an initial configuration for the servo-valves is consequently required. In this study, we used an arbitrary 30% valve opening to check if the regulation law and controller were able to rapidly match the commands, even if the initial thrust value was far from the target one. In a real system, the initial thrust command is known and the initial valve opening would be chosen in order to reach satisfying conditions from the very beginning of the engine's operation. Depending on the RBS parameters and the type of the regulation law, 5 to 7.5 s were needed to reach the thrust target in the simulations. Indeed, if the distance between the first resistor and the fuel grain surface is increased, the required time to start measuring the regression rate is increased and the beginning of the regulation is delayed.

(3) In order to integrate a valve response time in the simulations, a regulation frequency of 2 Hz has been used. This value seems realistic but could be changed if the valve and other sources of delays are known.

(4) The mass flow rate and the thrust measurements have been considered as ideal in the simulations. As we have seen in the error propagation analysis, the mass flow rate precision is an important parameter for the efficiency of the regulation and to limit the errors' propagation. Since the thrust is one of the quantities to be regulated, its estimation is naturally essential for the feedback loop control. Further analysis regarding these measurements should be integrated in future work. For example, artificial noise could be added and numerical filtering used to reproduce real in-flight analogical filtering. The noisy-filtered signals should then be used as measured thrust and oxidizer mass flow rate into the feedback regulator.

(5) Several aspects should be highlighted and discussed regarding the resistor-based sensors. In our study, we considered perfect sensors meaning that the resistors were disconnected as soon as the fuel surface reaches them. In practice, it would imply that the RBS do regress at the exact same rate as the fuel, and that the connectors for the resistors are destroyed properly at the same time. This assumption will be studied more in detail in future work within an experimental firing test campaign including RBS. Even if the sensors are not ideal in reality, it would be possible to consider some realistic conditions in the simulations. For example, to consider a potential difference between the fuel and the sensors' regression rates, a specific evaluation of the sensors' burning rate could be added. The simulations have been performed to analyze the influence of the sensors' radial distance between the resistors. However, other parameters could also be investigated such as the number of sensors and their location. The axial distance between resistors has been considered as null for the simulations but a more realistic value could be used. This work has not been performed yet, but the numerical code developed for this study could be used in its current form in future work to analyze such parameters of the RBS. As we have seen, the initial distance between the first resistor of the sensors and the fuel grain surface plays a major role in the responsiveness of the feedback control at the beginning of the combustion. For real applications, this distance should be minimized as much as possible.

(6) In the current version of the simulation code and feedback law, the fuel regression rate is supposed to be constant between two re-evaluations. Similarly, the regression rate is supposed to be piece-wise constant along the fuel grain. These hypotheses have been used as a first step since they allowed a simple evaluation of the fuel regression and mass flow rates, but they could be modified easily and improved by using interpolations.

(7) The PID controller has been limited to a proportional one and was sufficient to provide interesting results in the case of the simple configurations that have been analyzed in this study. The PID coefficients depend on the characteristics of the system and experimental tests would be required to determine them.

Throughout the discussion of this study, the combustion chamber pressure was not considered for the feedback control. In practice, the pressure should be measured and used for several purposes. Firstly, the regulation law should include safety procedures in which the chamber pressure would play a major role. Secondly, the combustion chamber pressure could be used to evaluate the combustion efficiency of the operating engine by combining the other measurements. The knowledge of the combustion efficiency could be taken into account in the regulation law. Finally, the thrust measurements provided by an accelerometer could be completed by an estimation based on the chamber pressure and the propellants' mass flow rates.

Despite the few limitations and possible improvements that have just been discussed, the current version of the code allowed for obtaining interesting results regarding the A-SOFT feedback control based on RBS measurements. Three thrust command profiles have been investigated (constant, piece-wise constant, and linear evolution) with a constant optimal mixture ratio for the O/F ratio command. For these configurations, the regulation law derived from the A-SOFT fundamental equations has been compared to a downgraded law based on a simple proportionality relation. Although the downgraded law exhibited lower performances than the full regulation law in terms of precision and oscillations around the target value, the behavior of the engine could be considered as acceptable with this simple law if the sensors' precision is high. Additionally, two sets of parameters for the RBS have been tested. The importance of the sensors design regarding the feedback control have been demonstrated as well as the advantage of using the complete regulation law. In particular, if the precision of the RBS is reduced, the feedback control of the A-SOFT engine remains efficient with limited oscillations if the full regulation law is used while the engine becomes difficult to control in the case where a downgraded law is used.

## **5. Conclusions**

The error propagation analysis has been performed in order to study the influence of the measurements errors on the thrust and mixture ratio control in hybrid engines. The A-SOFT and AOA concepts are based on a dual oxidizer injection and have been investigated and compared in terms of thrust, mixture ratio, and total relative errors. It has been demonstrated that the A-SOFT engines were more interesting to limit the errors propagation compared to the AOA motors. The influence of the geometric swirl number has been highlighted, and a possible method for its design was

proposed based on the error propagation results. The importance of the fuel and the oxidizer mass flow rates measurements has also been discussed and the use of resistor-based sensors has been investigated. Moreover, a design methodology has been proposed to make these sensors suitable for A-SOFT applications based on the error propagation results. A control law has been derived from the fundamental equations of hybrid engines and the A-SOFT concept and has been validated by numerical simulations. We demonstrated the interest of using this law to obtain precise feedback control of the engine compared to a simpler proportional regulation law. Besides the regulation law, the numerical simulations have been used to study the integration of the resistor-based sensors in the control of the engine. This study successfully demonstrated the feasibility of A-SOFT feedback loop control by using resistor-based sensors for the fuel regression rate measurement and represents an additional step towards the development of such hybrid rocket engines.

**Author Contributions:** Conceptualization, J.M. and T.S.; methodology, J.M. and T.S.; software, J.M.; validation, J.M.; formal analysis, J.M. and T.S.; investigation, J.M. and T.S.; resources, T.S.; data curation, J.M. and T.S.; writing—original draft preparation, J.M.; writing—review and editing, J.M. and T.S.; visualization, J.M. and T.S.; supervision, J.M. and T.S.; project administration, T.S.; funding acquisition, J.M. and T.S.

**Funding:** This research was funded by KAKENHI Grant No. JP16H04594 and supported by the External Funding Acquisition Incentive System provided by ISAS/JAXA.

**Acknowledgments:** The authors are grateful to their colleagues from SPLab (Politecnico di Milano), and particularly to Christian Paravan, for their collaboration on this research.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
