An Open-source Solver to Model the Catalytic Decomposition of Monopropellants for Space Thrusters

Low thrust propulsion systems are fundamental in different operations of a small spacecraft (Han et al. 2009). For that, there are three different types of technologies: cold gas, electric, and liquid monopropellant propulsion systems. A considerable amount of literature has been published, in which comparison about these propulsion systems, stating the advantages and drawbacks of each, can be found (Tummala and Dutta 2017; Lemmer 2017; Garrigues and Coche 2011; Martinez-Sanchez and Pollard 1998; Sutton and Biblarz 2016). Fundamentally, liquid monopropellant thrusters are the most commonly used type of propulsion systems for in-space missions. In comparison to electric thrusters, the liquid monopropellant ones have the advantages to have easy design and relatively high thrust, and, consequently, quick response in maneuvers. Profitability and high specific impulse are the main advantages compared to cold gas thrusters. A typical monopropellant thruster is composed of a simple feeding and control systems (Sutton and Biblarz 2016), as shown in Fig. 1. Decomposition of the monopropellant occurs in the catalytic chamber, which, by its turn, is composed of a catalytic bed with porous catalyst particles (grains or pellets) (Jang et al. 2012; Soares Neto et al. 2003; Vieira et al. 2005). The exhaust gases are expanded in the thrust chamber, generating thrust. https://doi.org/10.5028/jatm.v12.1111 ORIGiINAL PAPER


INTRODUCTION
The liquid hydrazine is the most commonly used propellant for monopropellant thrusters. When the hydrazine is decomposed by a suitable solid catalyst, it becomes an excellent both storable fuel and source of gases (Sutton and Biblarz 2016). However, since hydrazine is toxic, mutagenic, and carcinogenic, handling requires significant safety precautions, and special materials are necessary to store it, due to the fact that it is a highly hypergolic propellant. Due to all these disadvantages, new green monopropellants, or nontoxic propellants, with similar or better characteristics, have been developed, as an alternative to hydrazine. Research efforts are performed to replace this common monopropellant (Spores 2015;Tsay et al. 2016;Schmuland et al. 2013).
In general, the most critical point in the development of monopropellant thrusters is the experimental tests. These can be very expensive and, for hydrazine, certain handling precautions need to be taken. Therefore, theoretical and numerical simulations become important tools to reduce cost, risk, and development time, as well as to optimize the design of monopropellant thrusters.
Hydrazine monopropellant thrusters have been studied for more than 60 years. Thus, a great variety of theoretical, numerical, and experimental studies can be found in the literature (Kesten 1967;Shankar et al. 1985;Schmitz and Smith 1967;Crespo 1976;Kesten 1968;Shankar et al. 1984;Hwang et al. 2012). Comparative studies are possible since the thermophysical processes that occur in a catalytic bed are similar in all monopropellant thrusters. Efforts to establish a universal model, which also includes all type of liquid monopropellant thrusters (including hydrazine and green monopropellants), can be carried out. In specific terms, the catalytic chamber problem involves reactive flow through porous media, with energy transfer and multicomponent mass transfer.
CFD allows the numerical simulation of the reactive flow in porous media. The main commercial CFD tools are ANSYS Fluent ® and COMSOL. OpenFOAM ® is an open-source CFD tool, written in C++. The code can be adapted for different problems of fluid mechanics, and the numerical approach is based on the finite-volume method (FVM) (Moukalled et al.2016;Holzmann 2016). In OpenFOAM ® , the user can effciently manage the computational meshes and the discretization process of the governing equations of fluid mechanics. The framework can solve steady-state and unsteady flows. OpenFOAM ® is not intended to be a ready-to-use software, even though it might be supposed to be used as a standard simulation package (Goisis and Osio 2011).
Specifically, the solver default for reacting flow in OpenFOAM ® is ReactingFoam. Other default solvers with similar characteristics are fireFoam and ChemFoam. Several modifications on default solvers can be performed to the adaptation of particular problems of fluid mechanics with chemical reactions. Thus, new solvers for general catalyst reactors were developed (Maestri and Cuoci 2013;Kwiatkowski et al. 2013;Horgue et al. 2015). Classes inside the code for the solution of mass diffusion transfer were developed, including the model for fluid flow configuration involving multicomponent transport phenomena (Novaresio et al. 2012;Solsvik and Jakobsen 2011;de Melo 2013). In the case of the catalytic chamber, the use of different default OpenFOAM ® solvers is necessary, including modifications and interconnections.
According to the previously reported, this work aims to present a mathematical model and an open-source toolbox in OpenFOAM ® for the CFD modeling of catalytic decomposition inside catalytic chambers of monopropellant thrusters. Detailed model hypotheses and the mathematical model are given in second and third sections, respectively. The main characteristic of the new OpenFOAM ® solver, and the implementation and description of the case studies are given in the fourth section. Different steps were followed to modify the ReactingFoam solver adapting to the proposed model. The catalytic chamber was modeled as an axisymmetric segment. Profiles of temperature as a function of time to three initial bed temperatures were simulated. Mass fraction profile of species for a 294.4 K initial bed temperature were obtained and compared. Particles with two different diameters and two porosities were simulated and compared also. The verification of the proposed model was performed comparing the results of the simulation running with the Kesten work operation conditions (1969), and the experimental and simulation results of the same work.

MODEL DESCRIPTION
In this section, the physical problem and the hypotheses to the model are synthesized. In addition, hydrazine monopropellant thruster information is provided, because it was used as a reference for verification of the proposed model. Figure 2 shows the regions commonly considered for the study of the catalytic chamber. To simplify the model, in the figure, it is possible to observe that the gas-solid region covers most of the catalytic bed, and liquid-gas and liquid hydrazine regions are negligible compared to that region. Therefore, in this model, it was only considered the gas-solid region. What is more, at the limit between the liquid-gas region and the gas-solid region, the hydrazine in the gaseous state was considered. The distribution of pellets along the catalytic bed was considered uniform (isotropic).  In addition to general regions (liquid hydrazine, liquid-gas, and gas-solid regions), the catalytic bed is composed of three parts, shown in Fig. 2: the interstitial gas, that is the fluid moving through the bulk flow column and which is in the outside of the catalyst particles (a); the fluid static on the surface of the catalyst particles (b); and the catalyst, that corresponds to the solid phase and is stationary (c). In this work, two sets of balance equations were formulated for (a) and (band c together), respectively. In the last part, only chemical species (reagents and products) on the surface of the catalyst pellets were considered. The reason is that the chemical reaction occurs quickly on such a surface.
The monopropellant decomposition occurs through a catalyst material. The Shell-405 is the most commercially successful catalyst for hydrazine decomposition (Pakdehi et al., 2014). Two reactions steps are commonly considered. When hydrazine comes in contact with the catalyst surface, it is exothermically decomposed, producing nitrogen and ammonia (Eq. 1). Subsequently, the endothermic dissociation of ammonia takes place and produces nitrogen and hydrogen (Eq. 2) (Pakdehi et al. 2014;Kesten 1968): (1) The expression for Arrhenius rate is (3) where A 1 , T a e T are the reaction constant, the activation temperature, and the temperature of the region, respectively. Arrhenius parameters for hydrazine decomposition with the Shell-405 catalyst are given in Table 1 for both homogeneous and heterogeneous reactions. Finally, in heterogeneous catalytic processes, different steps are considered for the reactions to take place, which are detailedly described in the literature, where it is concluded that the slowest step controls the global reaction rate (Fogler et al. 2002). In this work, the mass diffusion was considered to affect the rate of reaction and also the dominating step (Kesten 1967). The mass transport between interstitial and static gas on the catalyst particle surface is fundamental for the catalytic bed modeling. The stationary thin film surrounding the catalyst particles was considered for the mass and heat transfers. Therefore, thin film resistance to heat and mass transfer to (or from) the interstitial gas is fundamental.

MATHEMATICAL MODEL
The numerical modeling of the catalytic chamber requires the solution of conservation equations for reactive gas flow in porous media. The fluid was assumed to be a continuum, fluid flow, Newtonian, unsteady, multicomponent, thermally-perfect mixtures of gases, and laminar. A set of equations for the interstitial gas and catalyst particle is given in Interstitial Gas Equations and Particle Equations, respectively.

Interstitial gas equations
The continuity equation for the interstitial gas in a porous media with constant porosity is given by (4) In the aforementioned equation, e i is the bed porosity, ρ i is the density in the interstitial gas, u is the velocity vector, and t is the time. Temporal and convective changes are represented by the first two terms on the left-hand side of the governing equations (these are the ∂( )/∂t and the ∇( ) terms, respectively).
The individual chemical-species balance in the interstitial gas is given by (5) The variables W a,i , W a,s , r' a,i , a s , P m and k c,a correspond to mass fraction of the species in the interstitial gas, mass fraction of the species in the surface gas, rate of the homogeneous reaction of the species, specific area, mean density, and mass transfer coeffcient of the species, respectively. The last term on the left of the Eq. 5 is due to mass diffusion in the interstitial gas. The first term on the right is the homogeneous reaction, which takes place in the interstitial gas mixture. Finally, the last term on the right is due to the mass transfer of interstitial gas from/to surface gas mixtures.
Particularly, for chemical-species balance equations, it was used the saturation fraction. Saturation is a relation between the total fluid and the total porous volume. For the interstitial gas, it is given by: (6) and for surface gas, by (7) where (8) Thereby, ε is + ε sg = 1.
The momentum-conservation equation in the interstitial gas is described by (9) where p, D a , μ i , F O denotes the pressure, the Darcy coeffcient, the viscosity in the interstitial gas, and the Forchheimer coeffcient, respectively.
The third term on the left of Eq. 9 is the Brinkman correction. The first term on the right represents a loss due to the gradient of pressure. The Darcy-Forchheimer law represents viscous and inertial losses in porous media. These are the second and third terms on the right of Eq. 9, respectively. Equation 9 is known as the non-Darcy model for flow in porous media (Das et al. 2018). The Darcy coeffcient and the Forchheimmer coeffcient are (Das et al. 2018;Hicks 1970) (10) where the permeability in an isotropic and homogeneous porous media is given by (Bear 2013;Ingham and Pop 1998) Finally, the energy balance equation in the interstitial gas is represented by (13) The third term on the left is due to flux-heat loss. The first term on the right corresponds to the heat supplied to the fluid by the homogeneous reaction. The second term is due to the pressure loss, and the third is the transport of heat between (a and c) regions. The last term is due to the transport of heat between (b and c) regions. Subscripts i, c, s, ∧ e ff refer to interstitial gas, solid matrix, surface gas, and effective, respectively.

Particle equations
The convective terms (∇( )) are not significant in the balance equations for the particle. The individual chemical-species balance for the surface gas is given by (14) The first and the second terms on the right of the Eq. 14 are due to the heterogeneous reaction and mass transfer between interstitial and surface gases, respectively.
Enthalpy balance in the surface gas (region (b)) is represented by (15) The first term on the right is due to the heterogeneous reaction on the surface gas (region (b)); the second is due to the transport of the heat between regions (b and c), and the last is the transport of heat between regions (a and c).
The temperature balance in the solid matrix (region (c)), is given by (16) The second term on the left of Eq. 16 is the flux-heat transport in the solid matrix. The first term on the right is due to the heterogeneous reaction in the solid matrix (region (c)), the second is due to the transport of the heat between regions (b and c), and the last is the transport of heat between regions (a and c).

Mass diffusion coeffcient
The mass diffusion coeffcient is a primordial parameter to calculate both the mass-flux between gases mixtures (interstitial from/ to surface gas) and mass diffusion-flux of individual species in the mixture. The diffusion-flux of individual species in the interstitial gas is important due to convective effects, which shall be considered in the equation per species of interstitial gas, and is given by (17) There are different ways to estimate the mass diffusion-flux and the diffusion coeffcient (D) in a mixture of gases. These were implemented already in a variety of codes (Novaresio et al. 2012;Solsvik and Jakobsen 2011;de Melo 2013). In OpenFOAM ® , the default diffusion model is based on a Schmidt number (Sc). This model results in analogies of the Prandtl number in heat transfer and of the Schmidt number in mass transfer, where ρDμ/Sca. In that model, a mass diffusion flow for a multicomponent mixture can be obtained by means of a general mass diffusion coeffcient of the bulk mixture. However, a more detailed per species mass diffusion coeffcient can be obtained through a specific mass diffusion coeffcient, D am .
When the mean of transport is a porous media, the diffusion of fluid through both interstices and the pores of the catalyst particles is more complicated. In a most rigorous analysis, one or more diffusion mechanisms must be considered, such as the bulk diffusion, the Knudsen diffusion, and the surface diffusion (Jakobsen 2014). Wilke-Bosanquet model established an effective diffusion coeffcient: where the effective Knudsen D eff ak and the effective bulk diffusion D eff am are and ε e τ are the porosity and the tortuosity, respectively.
A Knudsen number (the relation between the mean free path of the particle and the approximate pore diameter, K n = λ a /d pore ) is necessary to establish what diffusion coeffcient is relevant in a given porous media (i. e., Knudsen diffusion, bulk gas diffusion, or both). Knudsen diffusion predominates when K n > 10. Knudsen and bulk diffusions predominate when 0.1 < K n < 10. Bulk diffusion predominates when K n > 0.1.
The Kinetic Theory of the Gases (KTG) defines the Knudsen diffusion coeffcient as (21) where M a and d pore are the molecular weight of α-species and the mean pore diameter, respectively.
On the other hand, there are different ways to obtain the bulk diffusion coeffcient using theoretical or empirical correlations (Novaresio et al. 2012). In this work, the Chapman-Enskog theory was used, that is well-known in the transport models literature (Welty et al. 2009). Parameters for the calculation of the collision integral are according to Neufeld et al. (1972). The binary diffusion coeffcient was calculated first; subsequently, it was calculated the mean diffusion coeffcient of α-species in the mixture (Fairbanks and Wilke 1950).

Transport coeffcients
The transport coeffcients consist of empirical relations that depend on dimensionless numbers. Table 2 summarizes transport parameters used in the model. Table 2. Relation for transport parameters used in the model.

Mass transfer
Mass transfer coeffcient in m 2 Effective mass diffusion coeffcient and are effective mass diffusion coeffcients calculated for species α in the mixtures of interstitial and stagnant gases, respectively.

Reynolds number
Heat transfer

Heat transfer coefficient
Nusselt number Based on an analogy with Sherwood number for mass transfer, Nusselt number for heat transfer is the given (Wakao and Funazkri 1978).

Prandtl number
Mean parameters of the gas mixtures (Vafai 2015) Mean density

Mean dynamic viscosity
Mean Thermal conductivity Mean specific heat

CONSTITUTIVE EQUATIONS
The JANAF thermochemical tables (McBride et al. 1993;Chase Jr and Tables 1998) were used to calculate the properties of reactive mixture. The dynamic viscosity was calculated by means of the Sutherland equation. The equation of perfect gas was used in the model (ρ = p/RT). These are default options to the mixture and transport calculation in the OpenFOAM ® framework.

OVERVIEW OF THE SOLVER
The catalyticChamberFoamV.1 solves equations described in the third section of ths paper. For that, the reactingFoam solver was modified. Fundamentally, implemented modifications includes new dictionaries for entry of two gases names; and porous and transport parameters. The interstitial and surface gases were based on reactingFoam and chemFoam solvers, respectively. Dimensionless numbers (Sh, Sc, ℜ pa , Nu, and Pr), and mass and heat transport coeffcients (kc,α and h) are calculated inside the loop of the algorithm. The solution procedure follows a similar algorithm to the reactingFoam solver, and this includes the dictionaries and equations modified according to the model proposed. The main loop contains the system of equations to be solved, and two inner corrections for solving the pressure-momentum coupling, PISO (Pressure Implicit with Split Operator) and PIMPLE (merged PISO-SIMPLE) (Moukalled et al. 2016;Holzmann 2016), influenced by energy and species distributions.

TEST CASE DESCRIPTIONS
A typical structure of OpenFOAM ® cases includes three main folders: "0", "constant" and "system". In "0", initial and boundary conditions of main parameters are defined. The files in "constant" contains all model properties required by the OpenFOAM ® classes (chemistry, combustion, turbulence, and thermophysical), as well as additional transport and porous parameters created specifically for this model. Finally, in "system", it is defined the mesh geometry in blockMeshDict; the numeric methods of interpolation in fvSchemes; the information of iterative methods of algebraic system formed in fvSolutions; and solver running specifications in controlDict. Figure 3c shows the 2D axisymmetric computational mesh with its boundaries. The simulation can be performed over a 2D domain given that the catalytic bed has cylindrical symmetry. Preliminary convergence test was performed, aiming to test the stability of the solution concerning the adopted spatial discretization.
Spatial discretization and mesh resolution procedures consisted of (i) initial mesh definition with 50 × 10 segments (Fig. 3a); (ii) run the case; (iii) measuring of the temperature, pressure, and velocity values in two references points, one close to the wall and another close to the axis; (iv) increase mesh resolution in 2 times; (v) run the case again; (vi) measuring and comparing the previous and current temperature, pressure, and velocity values; (vii) repeat (iv), (v) and (vi) steps until the difference of temperature, pressure, and velocity values approximate to 2%. The final geometrical domain was divided into 400 segments on the axial direction and 80 segments in the radial coordinate, for a total amount of 3,200 cells (Fig. 3b). Euler implicit method and Gauss linear scheme were employed for the time discretization and spatial coordinates, respectively.

Case
Kesten's work includes the comparison between available experimental information about a typical flow in a catalytic bed with computational simulation. A series of calculation of a 222.2 N (50lbf) nominal thrust hydrazine monopropellant thruster were provided by the author (Kesten 1969). The catalytic bed packing has a 60.96 mm radius and 91.44 mm length. The 2D axisymmetric model geometrical dimensions were taken based on this catalytic bed (Fig. 3). The steady-state pressure was taken as 1.3785 × 10 6 Pa (200 Psia) the state-stable mass flow rate as 31 kg/m 2 s, and the initial pressure as 1.01352 × 10 5 Pa.
The catalyst bed packing consisted of 3.627 mm diameter spherical particles (25 to 30 mesh) in the first 6.35 mm of the catalytic bed, and 3.175 mm × 3.175 mm (1/8 × 1/8") cylindrical particles for the remainder of the bed. Kesten compared profiles of temperatures into the time, at a point inside the catalytic bed, and set up for 294.4 K, 527.7 K, and 788.8 K initial bed temperatures, respectively. He has obtained the profiles of mole fraction on the axial distance and for 3 seconds of the simulation. Kesten's work conditions were simulated in the present effort. In this work, we simulated spherical particles in all the catalyst bed packing, one of the cases having 3.627 mm diameter and the other 3.175 mm diameter. Details about of porosity inter-(e i ) and intraparticle (e s ) are not available in Kesten's work. However, the Rocket Research Corporation (Schmitz and Smith 1967) provides approximate values around 0.336 -0.45, respectively. Porosity interparticle, that is, the bed void fraction, can be obtained using an analytic relation valid to d pa /d bed ≤ 0.5 (Dixon 1988;Rusten et al. 2007): The code includes two options for the porosity interparticle (e i ) calculated with Eq. 22, or based on an experimental value. Since the porosity data is not available in Kesten's work, we used the method of porosity value calculated, as shown in Table 3. Table 4 provides all the other physical input parameters for all cases. In the present work, the initial and boundary conditions, which can be set up in folder "0" of the case, were chosen so that the numerical simulation was in complete agreement with the conditions reported by Kesten (1969). Air, as an inert species, with mass fraction equal to 1, occupying the internal field (which is a default option to begin the running case with a given combustion model in OpenFOAM®), 1.0135 × 10 5 Pa of pressure, and with initial bed temperatures (299.4 K, 527.7 K, or 788.8 K), provided the initial conditions to run the simulations. For the boundary conditions were established so that, at the inlet, the fraction of hydrazine mass is equal to 1, hydrazine enters with a velocity of 3.6 m/s, and 400 K of temperature. The wall was considered adiabatic, and the velocity was established with the noSlip option. At the outlet, the set temperature was the same as the ambient temperature, and the pressure was adjusted (1.25 × 10 6 Pa) so the final value was equal to the stable operation condition of Kesten's work (1.3785 × 10 6 Pa).

RESULTS AND DISCUSSION
Temperature plotted in function of time at a fixed axial position of 36.57 mm inside the catalytic chamber is shown in Figs. 4, 5, and 6 for the initial bed temperatures, set to 294.4 K, 527.7 K, and 788.8 K, 36.57 K corresponding to a point in the symmetry axis. These graphs also show a comparison between theoretical and experimental results obtained by Kesten, and the results obtained in the present work (Cases 1 and 2). The Kesten's results were used only to verify the implemented code concerning to a random axial point inside the mesh that represents the catalytic bed. This point, however, is a measuring reference in Kesten's work. Since Kesten's work provides experimental and theoretical results, this work pretends to attain similar results.
Temperature profiles for the present work were similar to the ones found in Kesten's work. However we observed that temperature stability were reached seconds before what was predicted by Kesten's model. Cases 1 and 2 were used to verify the effect of the particle diameter size in the temperature. In the first 0.5 seconds of the simulation, our results were closer to those in Kesten's work in case 2. After the first second, case 1 resembled more closely to Kesten's experimental and simulated results. It was possible to observe that the larger the catalyst particle diameter, the higher the temperature along the catalytic bed. In terms of application, the size particles comparison is key to determine the more adequate bed configurations and detect damages caused by different loads in the catalytic bed.
Temperature profiles describe the two reaction steps involved in the catalytic process (Kesten 1969;Shankar et al. 1984;Makled and Belal 2009). First, hydrazine decomposition is larger than the corresponding ammonia decomposition, and the net reaction is exothermic. Therefore, temperature increases until the endothermic decomposition of ammonia takes over.    Figure 7 shows the mass fractions profile at the axial position at around 3 seconds on the simulation with 294.4 K as initial bed temperature (which is associated to Fig. 4). Theoretical results obtained by Kesten and in the present work are shown. Mass fraction results in this work are very different than the ones in Kesten's work, where hydrazine mass fraction was shown to be decaying exponentially along to the axial position, until it is completely consumed. In this work, the results have shown that this not occurs. Some hydrazine rest remains inside the catalytic bed. Since mass fraction profile obtains from Kesten's work a part that is not experimental data, it is very diffcult to determinate whether or not the work is correct. Even today, we need experimental results to know details about mass fraction information.
In conformity with the previously stated, the closest match of mass fraction for ammonia, nitrogen, and hydrogen was found after the 0.04 m of axial distance. The differences between the two works can be attributed to how the transport coeffcients were calculated. In Kesten's work, mass diffusion coeffcients were based on some inlet values per chemical species, whilst, in this work, mass diffusion coeffcients of chemical species were calculated during the time simulation, as shown in the section about the mass diffusion coefficient. Another difference is that the transport coeffcients were calculated using the most updated equations for the study of fixed catalytic beds. These were different than the ones used in Kesten's work.
Axial distance (m)  Figure 7. Axial profiles of species mass fraction at 3 seconds of simulation, being 294.4 K the initial bed temperature and case 1 porous media parameters.

CONCLUSION
A CFD modeling for catalytic chambers of monopropellant thrusters were proposed. The model were implemented in OpenFOAM ® framework and named catalyticChamberFoamV.1. Physical and thermochemical properties of the flow in a porous media can be included in the code.
Said code allows the prediction of mass fraction of chemical species involved in catalytic processes, as well as of fields of temperature along the axis of the catalytic bed. It also works as a function of time for a monopropellant thruster in ground testing conditions, which are all fundamental parameters to evaluate the performance and aid the design of these thrusters.
A relevant literature reference was used for verifying the code and similar transient temperature profiles, and less close axial profiles of mass fractions were found. Although this work uses hydrazine as a monopropellant thruster to test the code, it can be used with other monopropellants, and catalysts with different porous media characteristics.