A Novel Catalytic Micro-Combustor Inspired by the Nasal Geometry of Reindeer: CFD Modeling and Simulation

: A three-dimensional CFD model of a novel conﬁguration of catalytic micro-combustor inspired by the nasal geometry of reindeer was developed using the commercial code ANSYS Fluent 19.0. The thermal behavior of this nature-inspired (NI) conﬁguration was investigated through simulations of lean propane / air combustion performed at di ﬀ erent values of residence time (i.e., inlet gas velocity) and (external convective) heat transfer coe ﬃ cient. Simulations at the same conditions were also run for a standard parallel-channel (PC) conﬁguration of equivalent dimensions. Numerical results show that the operating window of stable combustion is wider in the case of the NI conﬁguration. In particular, the blow-out behavior is substantially the same for the two conﬁgurations. Conversely, the extinction behavior, which is dominated by competition between the heat losses towards the external environment and the heat produced by combustion, di ﬀ ers. The NI conﬁguration exhibits a greater ability than the PC conﬁguration to keep the heat generated by combustion trapped inside the micro-reactor. As a consequence, extinction occurs at higher values of residence time and heat transfer coe ﬃ cient for this novel conﬁguration.


Introduction
Hydrocarbon-fueled micro-combustors have been proposed to replace batteries in portable electronic devices [1,2], their main advantage being related to the high gravimetric energy density of hydrocarbons. In addition, micro-combustors/reactors in parallel arrangement can be used as an alternative to conventional devices to achieve an inherently safer operation [3].
The main issue of micro-combustors is thermal management. When reducing the reactor scale from macro to micro, the surface area-to-volume ratio increases. This may lead to an increase of heat losses towards the external environment with respect to the heat generated by combustion, thus reducing the operating window of stable combustion. In order to enlarge the stability limits, the catalytic coating of the walls of the micro-combustor has been proposed. It has been shown that the presence of a catalyst allows for a significant improvement in thermal performance with respect to homogeneous micro-combustors [2]. However, the stability of catalytic micro-combustors still remains an issue [4][5][6][7][8][9].
In order to improve the stability of micro-combustors, several alternative configurations have been proposed. In particular, to enhance the heat recirculation from the hot zones to the cold zones, reverse-flow [10], heat-recirculating [11,12] and Swiss-roll [13,14] micro-reactors were investigated.
On the basis of the localization of heat losses, a selective coating of the channels of the micro-combustor has also been proposed, thus saving catalyst: since the main source of heat losses is localized at the external surface of the outer channels, only these channels were coated with the catalyst, whereas the inner channels were left uncoated, being protected from heat losses by the outer ones [15][16][17][18]. It has been shown that, in partially coated configurations, the catalyst can be used essentially as a "pilot" of homogeneous combustion warranting ignition and, at the same time, a relatively higher temperature and, thus, improved stability [15][16][17][18][19].
With the aim of developing more stable catalytic micro-combustors, recently, Kaisare's group studied a spiral configuration by means of both two-dimensional [20] and three-dimensional [21] CFD simulations. This configuration was obtained by curling up a straight channel into a spiral, thus providing large ratios of internal heat exchange area to external heat loss area. In the spiral micro-reactor, the flow enters the spiral geometry through the inlet at the center and exits tangentially through the outlet. Comparison with the straight channel has shown that the spiral is more stable due to two mechanisms: transverse heat recirculation taking place from hot products to the cold feed between adjacent turns of the spiral, and better protection of the central combustion zone from heat losses [20,21]. The spiral is indeed characterized by a higher temperature and a higher contribution from homogeneous chemistry [21].
In this context, nature can play a significant role in inspiring more energy-efficient equipment. Magnanelli et al. [22] have shown how reindeer can inspire a novel industrial design with high heat efficiency. Indeed, reindeer in the arctic region live under very harsh conditions and may face temperatures below −40 • C. Therefore, in order to survive, conservation of body heat is crucial. The reindeer nasal geometry/mechanism for heat and mass exchange during respiration plays a key role in this respect. The perpendicular section of the reindeer nose has a spiral geometry. Air enters perpendicularly to the spiral and is efficiently and quickly heated.
Inspired by the efficient heat management of the reindeer nasal geometry, we here propose a novel spiral catalytic micro-combustor. To investigate the thermal behavior of this nature-inspired configuration, we developed a three-dimensional CFD model. Simulations of lean propane/air combustion were performed at different values of inlet gas velocity (i.e., residence time) and (external convective) heat transfer coefficient. For the sake of comparison, simulations at the same conditions were also run for a standard parallel-channel (PC) configuration of equivalent dimensions.

Results and Discussion
In this section, simulation results obtained for the parallel-channel (PC) configuration and the nature-inspired (NI) configuration are presented and discussed. We start from the results of simulations run at different values of inlet gas velocity, v in (i.e., residence time, τ = length of the reactor/inlet gas velocity-L/v in ) (Section 2.1). These results show differences in terms of thermal behavior between the two configurations. The nature of these differences is supposed and, then, confirmed by the results of simulations run at different values of (external convective) heat transfer coefficient, h ext (Section 2.2).
2.1. Effect of the Inlet Gas Velocity (i.e., Residence Time) Figures 1 and 2 show the maps of (top) propane molar fraction and (bottom) temperature as calculated over different cross-sections of the micro-reactor at different values of inlet gas velocity, v in (i.e., residence time, τ), for the parallel-channel (PC) configuration and the nature-inspired (NI) configuration, respectively.
There are aspects common to both configurations. In particular, from the maps of propane molar fraction, it can be seen that, as the inlet gas velocity is increased, the reaction front is progressively shifted downstream until it reaches a condition of non-complete fuel conversion that is well evident at v in = 30 m/s. Accordingly, at a fixed axial position, temperature decreases with increasing inlet gas velocity. Indeed, as the inlet gas velocity is increased, the residence time decreases and, thus, blow-out is approached. Overall, the results of Figures 1 and 2 suggest that the blow-out behavior is There are aspects common to both configurations. In particular, from the maps of propane molar fraction, it can be seen that, as the inlet gas velocity is increased, the reaction front is progressively shifted downstream until it reaches a condition of non-complete fuel conversion that is well evident at vin = 30 m/s. Accordingly, at a fixed axial position, temperature decreases with increasing inlet gas velocity. Indeed, as the inlet gas velocity is increased, the residence time decreases and, thus, blow-out is approached. Overall, the results of Figures 1 and 2 suggest that the blow-out behavior is substantially the same for both PC and NI configurations. This is confirmed when looking at Figures 3 and 4 showing the maximum wall temperature, Twall_max, the bulk In both cases, blow-out occurs at the same value of residence time (τ = 0.0033 s, i.e., v in = 15 m/s). Conversely, as far as extinction is concerned, the two configurations exhibit rather different behaviors. In particular, extinction occurs at τ = 0.5 s (i.e., v in = 0.1 m/s) in the case of the PC configuration and at a higher value of residence time, τ = 1 s (i.e., v in = 0.05 m/s), in the case of the NI configuration. These differences in extinction behavior are likely to be related to the greater ability of the NI configuration to keep the heat generated by combustion trapped inside the micro-reactor.
Catalysts 2020, 10, x FOR PEER REVIEW 5 of 14 temperature (at the exit section), Tbulk, and the propane conversion, xC3H8, as a function of the residence time, , for the PC configuration and the NI configuration, respectively.   temperature (at the exit section), Tbulk, and the propane conversion, xC3H8, as a function of the residence time, , for the PC configuration and the NI configuration, respectively.

Effect of the (External Convective) Heat Transfer Coefficient
On the basis of the results obtained from simulations run by varying the inlet gas velocity (i.e., the residence time), we have supposed that the differences between the PC configuration and the NI configuration stem mainly from their different ability to keep the heat produced by combustion trapped inside the micro-reactor. In order to confirm this supposition, we performed simulations at different values of (external convective) heat transfer coefficient. For both configurations, T wall_max , T bulk , and x C3H8 are plotted as a function of the dimensionless parameter R in Figures 5-7, respectively. The parameter R is defined as where h ext is the heat transfer coefficient, A ext is the external surface area, T in is the inlet gas temperature, Q is the inlet flow rate of propane, and ∆H comb is its heat of combustion.
the residence time), we have supposed that the differences between the PC configuration and the NI configuration stem mainly from their different ability to keep the heat produced by combustion trapped inside the micro-reactor. In order to confirm this supposition, we performed simulations at different values of (external convective) heat transfer coefficient. For both configurations, Twall_max, Tbulk, and xC3H8 are plotted as a function of the dimensionless parameter R in Figures 5-7, respectively. The parameter R is defined as where hext is the heat transfer coefficient, Aext is the external surface area, Tin is the inlet gas temperature, Q is the inlet flow rate of propane, and Hcomb is its heat of combustion.  As R increases (i.e., as the heat losses towards the external environment increase with respect to the heat produced by combustion), extinction is approached. In the case of the PC configuration, extinction occurs at R = 0.013. Conversely, for the NI configuration, extinction is not found under the simulation conditions investigated, confirming that this configuration exhibits much better thermal stability and performance than the PC configuration. As R increases, both Twall_max and Tbulk decrease; As R increases (i.e., as the heat losses towards the external environment increase with respect to the heat produced by combustion), extinction is approached. In the case of the PC configuration, extinction occurs at R = 0.013. Conversely, for the NI configuration, extinction is not found under the simulation conditions investigated, confirming that this configuration exhibits much better thermal stability and performance than the PC configuration. As R increases, both T wall_max and T bulk decrease; however, x C3H8 remains very close to 100% over the whole range of R values explored. In addition, at fixed R, both T wall_max and T bulk are higher in the case of the NI configuration. This trend can also be observed from the temperature maps of Figure 8, computed for both PC and NI configurations at two different values of R.
Catalysts 2020, 10, x FOR PEER REVIEW 8 of 14 however, xC3H8 remains very close to 100% over the whole range of R values explored. In addition, at fixed R, both Twall_max and Tbulk are higher in the case of the NI configuration. This trend can also be observed from the temperature maps of Figure 8, computed for both PC and NI configurations at two different values of R. Global assessment of Figures 5-8 thus confirms that the NI configuration is more resistant to extinction than the PC configuration due to its greater ability to keep the heat generated by combustion trapped inside the micro-reactor. Global assessment of Figures 5-8 thus confirms that the NI configuration is more resistant to extinction than the PC configuration due to its greater ability to keep the heat generated by combustion trapped inside the micro-reactor.

Mathematical Model
The thermal behavior of the two configurations of catalytic micro-combustor shown in Figure 9, the parallel-channel (PC) configuration (top) and the nature-inspired (NI) configuration (bottom), was simulated through a three-dimensional CFD model. Figure 9 also shows an image of the cross-section of a reindeer nose (at around 10 cm into the nose), from which we took inspiration for the novel micro-reactor configuration proposed here. Notice the spiral shape of the nasal turbinates. It is worth stressing that, in the case of the reindeer nose, air enters perpendicularly to the spiral. The same occurs when the fuel/air mixture enters the micro-reactor with the NI configuration. Conversely, in the case of the spiral micro-combustor investigated by Kaisare's group, the flow enters the spiral geometry through the inlet at the center and exits tangentially through the outlet [20,21], thus following a path that is completely different from that of the spiral configuration of Figure 9. Table 1 gives the values of the most important geometrical parameters for both configurations: the distance between catalytic walls, d; the thickness of the catalytic walls, w; the length of the catalytic micro-combustor, L; the total catalytic surface area, Acat; the ratio between the total catalytic surface area and the volume of the micro-reactor, Acat/V; and the external surface area, Aext. These parameters were set to the same values, or to values as close as possible, compatibly with the constraints imposed by the domain construction for the two different configurations.

Parameter
Value Distance between catalytic walls, d (mm) 0.783 d  w Figure 9. Computational domains of (top) parallel-channel (PC) configuration and (bottom) nature-inspired (NI) configuration; the magnified image with red edges shows some geometrical details; an image of the cross-section of a reindeer nose (at around 10 cm into the nose) is also shown (from Magnanelli et al. [22]).
It is worth stressing that, in the case of the reindeer nose, air enters perpendicularly to the spiral. The same occurs when the fuel/air mixture enters the micro-reactor with the NI configuration. Conversely, in the case of the spiral micro-combustor investigated by Kaisare's group, the flow enters the spiral geometry through the inlet at the center and exits tangentially through the outlet [20,21], thus following a path that is completely different from that of the spiral configuration of Figure 9. Table 1 gives the values of the most important geometrical parameters for both configurations: the distance between catalytic walls, d; the thickness of the catalytic walls, δ w ; the length of the catalytic micro-combustor, L; the total catalytic surface area, A cat ; the ratio between the total catalytic surface area and the volume of the micro-reactor, A cat /V; and the external surface area, A ext . These parameters were set to the same values, or to values as close as possible, compatibly with the constraints imposed by the domain construction for the two different configurations. The three-dimensional CFD model adopted in this work simulates the coupling of the fluid flow and the chemical processes at the gas-solid interface and in the gas phase for lean propane/air combustion. It solves the mass, momentum, chemical species, and energy conservation equations in the fluid (coupled with the ideal gas equation), along with the energy conservation equation in the solid walls. Laminar flow was assumed. (The validity of this assumption is discussed at the end of this section.) Steady-state simulations were run. The model equations are listed below (vector notation was used).

Fluid Phase
Mass where ρ is the density and where p is the static pressure and τ is the stress tensor given by where µ is the molecular viscosity and I is the unit tensor. (The second term on the right-hand side is the effect of volume dilation.) Chemical species-for the ith species (an equation was solved for the N-1 species, where N is the total number of chemical species present in the system): where Y i , → J i , and R hom,i are the mass fraction, the diffusion flux, and the net rate of production by homogeneous (i.e., volumetric) chemical reaction of the ith species, respectively. → J i is given by where D is the diffusion coefficient. Energy with where h is the sensible enthalpy defined as where C p,i is the specific heat of the ith species and T ref is the reference temperature (298.15 K). In Equation (7), k is the fluid thermal conductivity and S h is the source of energy due to volumetric chemical reaction.

Solid Phase
Energy where k w is the solid thermal conductivity and T w is the solid temperature. At the inlet, a fixed flat velocity profile was assumed. For species and energy, Danckwerts boundary conditions were used. At the exit, the static pressure was imposed, and far-field conditions were specified for the remaining variables.
At the fluid-wall interface, a no-slip boundary condition was assigned (the fluid has zero velocity relative to the boundary) and coupled to the species balances (the mass flux of each gas species due to diffusion to/from the surface is balanced with its rate of consumption/production on the surface): and the energy balance: where ω h is the superficial heat production rate. Newton's law of convection was used at the outer surface of the walls: where h ext is the (external convective) heat transfer coefficient, T w,ext is the temperature at the external wall surface, and T a,ext is the external temperature (300 K). The reaction rate of homogeneous propane combustion was calculated according to the single-step reaction rate by Westbrook and Dryer [23]: where the activation energy is expressed in J/kmol and the species concentrations in kmol/m 3 . A single-step reaction rate was also assumed for catalytic propane combustion [24]: where the activation energy is expressed in J/kmol and the concentration of propane in kmol/m 2 . For the molecular viscosity, the temperature dependence reported by Canu [25] for nitrogen was adopted. The fluid specific heat and thermal conductivity were calculated by a mass fraction-weighted average of species properties. The specific heat of each species was calculated as a piecewise fifth-power polynomial function of temperature.
The model equations were discretized using a finite-volume formulation on the structured meshes shown in Figure 10. The mesh convergence was checked, and halving the element size resulted in differences within 10%.

 
. R (16) where the activation energy is expressed in J/kmol and the concentration of propane in kmol/m 2 . For the molecular viscosity, the temperature dependence reported by Canu [25] for nitrogen was adopted. The fluid specific heat and thermal conductivity were calculated by a mass fraction-weighted average of species properties. The specific heat of each species was calculated as a piecewise fifth-power polynomial function of temperature.
The model equations were discretized using a finite-volume formulation on the structured meshes shown in Figure 10. The mesh convergence was checked, and halving the element size resulted in differences within 10%. The spatial discretization used first-order schemes for all terms, except for the diffusion terms that were treated with a second-order central difference scheme. Parallel computations were performed by means of the segregated solver of the CFD code ANSYS Fluent 19.0 (https://www.ansys.com/products/fluids/ansys-fluent).
Simulations were run at different values of inlet gas velocity, vin (i.e., residence time,  = length of the reactor/inlet gas velocity-L/vin), and heat transfer coefficient, hext. The effect of vin was explored by setting hext = 20 W/m 2 K, whereas the effect of hext was explored by setting vin = 0.5 m/s. Further simulation conditions are detailed in Table 2. In all the simulations, the flow is laminar. A Reynolds number, Re, can be defined assuming the inlet velocity as the mean flow velocity and the distance between catalytic walls as the characteristic length: The highest value of Re is obtained at the inlet of the micro-combustor, i.e., at the lowest value of temperature, since as temperature increases, the gas density decreases and the molecular viscosity The spatial discretization used first-order schemes for all terms, except for the diffusion terms that were treated with a second-order central difference scheme. Parallel computations were performed by means of the segregated solver of the CFD code ANSYS Fluent 19.0 (https://www.ansys.com/products/ fluids/ansys-fluent).
Simulations were run at different values of inlet gas velocity, v in (i.e., residence time, τ = length of the reactor/inlet gas velocity-L/v in ), and heat transfer coefficient, h ext . The effect of v in was explored by setting h ext = 20 W/m 2 K, whereas the effect of h ext was explored by setting v in = 0.5 m/s. Further simulation conditions are detailed in Table 2. Table 2. Further simulation conditions.

Parameter Value
Inlet gas temperature, T in (K) 300 Inlet fuel equivalence ratio (-) 0.5 Solid thermal conductivity, k w (W/m K) 18 1 1 Walls were assumed as made of silicon carbide (SiC).
In all the simulations, the flow is laminar. A Reynolds number, Re, can be defined assuming the inlet velocity as the mean flow velocity and the distance between catalytic walls as the characteristic length: The highest value of Re is obtained at the inlet of the micro-combustor, i.e., at the lowest value of temperature, since as temperature increases, the gas density decreases and the molecular viscosity increases. At the inlet temperature of 300 K and at atmospheric pressure, the density and viscosity of the fluid (this latter is calculated with the formula reported by Canu [25] for nitrogen) are equal to 1.18 kg/m 3 and 1.78 × 10 −5 Pa s, respectively. Thus, in the velocity range investigated (0.05-30 m/s), the highest value of Re is equal to around 1560, which is well below the limit value of transition to turbulent regime.

Conclusions
A three-dimensional CFD model of a novel catalytic micro-combustor was developed. The geometrical configuration of the novel micro-combustor was inspired by the nose of reindeer, which efficiently heats air, thus allowing their survival in the arctic region in spite of very harsh temperature conditions (even below −40 • C). The thermal behavior of this nature-inspired (NI) configuration was investigated through simulations of lean propane/air combustion performed at different values of residence time (i.e., inlet gas velocity) and (external convective) heat transfer coefficient. For the sake of comparison, simulations at the same conditions were also run for a standard parallel-channel (PC) configuration of equivalent dimensions.
Numerical results have shown that blow-out is not significantly affected by the geometrical configuration. Conversely, extinction, which is dominated by competition between the heat losses towards the external environment and the heat produced by combustion, occurs at higher values of residence time and heat transfer coefficient in the case of the NI configuration. This is due to the fact that the NI configuration exhibits a greater ability than the PC configuration to keep the heat generated by combustion trapped inside the micro-reactor.
These results confirm that nature can play a key role in inspiring energy-efficient configurations to improve the thermal stability and performance of catalytic micro-combustors.