Flow Simulation-Based Methodology for Reducing The Risk of Fuel Fire In An Aircraft’s Fuel System Enclosure

Certain elements of the fuel system are located in the equipment bay (fuel system enclosure) situated under the floor of the pilot cabin, which is marked in Fig. 1 with a red, dashed line. The complete fuel supply system also includes tanks in the upper wings and elements under the engine cowling, which were outside the scope of the present work. The location of the elements of the existing equipment bay ventilation system is shown in Fig. 1.

Fig. 1.

As shown in Fig. 1, the existing drainage and ventilation system of the equipment bay consists of an inlet A, two larger outlets with fairings, B and C, near the front bulkhead of the bay, and two smaller outlets with fairing E and D in the bottom surface, near the rear bulkhead. The larger inlet with scoop in front of inlet A exclusively serves the cabin ventilation and is not available for the ventilation of the equipment bay. The arrows F and G indicate the location of two pairs of 2 mm-diameter drainage holes in the skin of the fuselage, located on each side of the symmetry plane, in front of two bulkheads of the structure, which are intended to allow collected fluids to flow out of the bay.

The strategy adopted for minimizing fire hazard due to fuel leaks consisted of:

Identifying potential locations of leaks,

It was assumed that the most dangerous location for potential leaks was between the last two bulkheads of the equipment bay, where the electric-driven starter/reserve fuel pump is located. This space is shown in Fig. 2. The worst-case scenario could occur if the fuel vapor concentration rose above the flammability limit (1.5% volume fraction). In such a case, the vapor might be ignited by sparks from the commutator of the electric motor driving the fuel pump, or from the connections of the electric batteries. This pump, powered by an electric motor, operates primarily during engine start. After engine start, it may serve as a reserve pump, in case of failure of the main pump located close to engine, outside of equipment bay.

Fig. 2.

Based on the designers’ experience with the volume of a possible leak in given location, the maximum value of mass flow of the simulated leak was set to 2% of the maximum mass flow generated by the starter/reserve fuel pump, or 2.6 g/s.

In the event of a fuel leak, the desired scenario of for minimizing the fire hazard should facilitate removing the rapid removal of liquid fuel reaching the bottom surface of the bay, with minimal wetting of the surface. Moreover, it should also ensure the removal of fuel vapor without creating regions of vapor concentration above the flammability limit near the potential ignition sources.

The flow and phase-change simulations were conducted using the ANSYS Fluent solver, incorporating specific sub-models for different phenomena: Species, Discrete-Particle Model (DPM) and Eulerian Wall Film (EWM) model. The Species model provides the physical properties of gases (oxygen, nitrogen, fuel vapor) and of liquid fuel. The Discrete Particle Model (DPM) makes it possible to model input of the leaking fuel into the computational domain in the form of discrete particles. The input details of DPM include the coordinates of injection points, droplet diameter, direction and cone angle of the stream. After the droplets reach the bottom surface of the equipment bay, they form an Eulerian Wall Film (EWF). This is driven by gravity and the tangential stress on the upper film surface exerted by the airflow. The EWF model is coupled with the phase change model, which determines the evaporation rate and makes it possible to track the vapor flow inside the domain. In the present study, the Diffusion-Balance phase change model was used with standard settings.

The species applied in the simulations included air components (O₂ and N₂), liquid fuel, and fuel vapor. The properties of the aviation fuel AVGAS 100LL applied in the simulations (most importantly: the dependence of partial pressure of fuel vapor on temperature) were not available in the solver database, and were adopted from the literature [1], [8], [9] and used in over-writing the data for iso-octane (C₈H₁₈) from the ANSYS material database. The dependence of fuel vapor partial pressure on temperature was estimated based on the formula: 1 p=exp{ [ 0.7553−413.0T+459.6 ]S0.5log10(RVP)−[ 1.854−1042T+459.6 ]S0.5++[ 2416T+459.6−2.013 ]log10(RVP)−8742T+459.6+15.64 } \[p=\exp \left\{ \begin{array}{*{35}{l}} \left[ 0.7553-\frac{413.0}{T+459.6} \right]{{S}^{0.5}}{{\log }_{10}}\left( RVP \right)-\left[ 1.854-\frac{1042}{T+459.6} \right]{{S}^{0.5}}+ \\ +\left[ \frac{2416}{T+459.6}-2.013 \right]{{\log }_{10}}\left( RVP \right)-\frac{8742}{T+459.6}+15.64 \\ \end{array} \right\}\]

from ref. [8], where:

· p is the fuel vapor partial pressure, [psi];

The flow simulation of gas species was conducted using the ANSYS Fluent URANS solver with the k-omega SST turbulence model. As far as the Eulerian Wall Film model is concerned, the Fluent solver determines the surface distribution of the thickness and velocity of liquid film from the conservation equations for mass, momentum, and energy [10]: 2 ∂ρlh∂t+∇s⋅(ρlhV→l)=m˙s, \[\frac{\partial {{\rho }_{l}}h}{\partial t}+{{\nabla }_{s}}\cdot \left( {{\rho }_{l}}h{{{\vec{V}}}_{l}} \right)={{\dot{m}}_{s}},\] 3 ∂ρlhV→l∂t+∇s⋅(ρlhV→lV→l+DV→)=−h∇sPL+ρlhgτ→+32τ→fs−3μlhV→l+q˙→s+τ→θw, \[\frac{\partial {{\rho }_{l}}h{{{\vec{V}}}_{l}}}{\partial t}+{{\nabla }_{s}}\cdot \left( {{\rho }_{l}}h{{\overrightarrow{V}}_{l}}{{\overrightarrow{V}}_{l}}+\overrightarrow{{{D}_{V}}} \right)=-h{{\nabla }_{s}}{{P}_{L}}+{{\rho }_{l}}h\overrightarrow{{{g}_{\tau }}}+\frac{3}{2}{{\vec{\tau }}_{fs}}-\frac{3{{\mu }_{l}}}{h}{{\vec{V}}_{l}}+{{\vec{\dot{q}}}_{s}}+{{\vec{\tau }}_{\theta w}},\] 4 ∂ρlhTf∂t+∇s⋅(ρlhTfVl→+DV→)=1cp[ 2kfh(Ts+Tw−2Tm)+q˙imp+m˙vapL ]. \[\frac{\partial {{\rho }_{l}}h{{T}_{f}}}{\partial t}+{{\nabla }_{s}}\cdot \left( {{\rho }_{l}}h{{T}_{f}}\overrightarrow{{{V}_{l}}}+\overrightarrow{{{D}_{V}}} \right)=\frac{1}{{{c}_{p}}}\left[ \frac{2{{k}_{f}}}{h}\left( {{T}_{s}}+{{T}_{w}}-2{{T}_{m}} \right)+{{{\dot{q}}}_{imp}}+{{{\dot{m}}}_{vap}}L \right].\]

In equations (2)-(4), ρ_l is liquid density, [kg/m³]; h is film thickness, [m]; ∇_s is surface gradient; Vl→\[\overrightarrow{{{V}_{l}}}\] is the average fluid velocity, [m/s]; ṁ_s is the film mass source, [kg/s]. DV→\[\overrightarrow{{{D}_{V}}}\] is the advection tensor, accounting for parabolic distribution of velocity in the film, p_L is pressure exerted on the film, [Pa]; τ→fs\[{{\vec{\tau }}_{fs}}\] is viscous stress on the gas-liquid border, [Pa]; τ→θw\[{{\vec{\tau }}_{\theta w}}\] is tension due to collection or separation of droplets, [Pa]. q˙imp\[{{\dot{q}}_{imp}}\] is the heat produced by slowing down the droplets, [J/s]; L is the latent heat of phase change, [J/kg].

It was assumed that the most critical conditions for the ventilation of the equipment bay occur during the engine start and warm-up procedure, when the airplane is standing on the ground and the ventilation of the bay is achieved solely due to external flow produced by the propeller. Due to the lack of actual geometry for the propeller, its geometry, as well as the geometry of the airplane, were re-created by three-dimensional scanning of the test-bed airplane.

The thrust obtained from the Fluent model of the running propeller on a stationary airplane was validated by comparison with the estimated thrust produced by two other four-blade-propellers from report [11], the first designated 5868-R6 with an R.A.F.6 airfoil and the second designated 5868-9 with a CLARK-Y airfoil. The thrust computed with CFD (using the Moving Reference Frame) for the scanned propeller under test conditions (stationary airplane, n=700 rotations/min) was 8% lower than the thrust estimated for the first propeller and 3% lower than that estimated for the second propeller. It was concluded that the propeller model obtained from 3D-scanning of the actual propeller could be used for simulating the airflow in the vicinity of the airplane without a significant risk of overprediction of dynamic pressure.

Simulating airflow through the equipment bay induced by the airplane propeller presents a challenge due to the presence of objects of varying sizes within the same domain, influencing the flow. These include the airplane fuselage, a large object, as well as a large number of smaller elements of various airplane systems, including elements of the fuel system in the equipment bay. A computational mesh resolving all these objects would require a large number of cells. Therefore, it was decided that the computational model would consist of three sub-models:

A “big” model encompassing a mesh for the external domain, shown in Fig. 3, including the geometry of the airplane with the running propeller, the frontal part of the engine bay, the fuselage with inlets and outlets to the equipment bay, the wings, and the tail surfaces. This mesh covered a semi-spherical domain within 50 meters from airplane, with the plane located at its centre. The purpose of this model was to determine the pressure and velocity field around the airplane, in order to use them as boundary conditions for smaller models.

Fig. 3.

Fig. 4.

Fig. 5.

Fig. 6.

Due to limitations of computational time, the degree of sensitivity of the results of fuel drainage and removal of fuel vapor to the mesh density was analyzed within a domain of approximately one-third the volume of the equipment bay: the volume between the two last bulkheads, which for this purpose was separated from the rest of the bay by a fictitious wall surface. The position of the dividing bulkhead is indicated with blue dashed line in Fig. 6. This volume had two inlets of air, one already existing, and another introduced in the project, whose position is indicated by blue ellipse in Fig. 6. This volume had two outlets of fluid and gas, the ones visible in Fig. 6 near the rear bulkhead. The bay floor at the juncture with bulkheads was modified with a step preventing sidewise spill of the fuel, visible in Fig. 7. This feature was introduced during the present study and was absent in the baseline version of the bay. The meshes used in this investigation were prepared in the Fluent Meshing tool, operating in the classic mode, which has the option of scaling the mesh size field (the set of all parameters defining the number and distribution of mesh nodes). The size field was defined by such parameters as minimum and maximum element sizes on the given surface, the element growth ratio, and definition of the boundary layer associated with the given surface. The crudest mesh was defined with maximum and minimum element sizes on the largest surface of 44.5 and 0.17mm, respectively a growth ratio of 1.2, and first boundary layer element height of 0.1mm. For increasing mesh density the size field was scaled by factors of 0.9, 0.8 and 0.75 applied to the size field, which resulted in mesh sizes of 1.5 million, 1.9 million, 2.36 million, and 2.5 million elements.

Fig. 7.

In order to simplify the modelling of the fuel film in this part of the study, the modelling of droplet flow was omitted. Instead, the fuel film originated from a predefined zone in the floor of the bay, shown in Fig. 7, where appropriate source terms for the fuel film were defined. The film was then driven by gravity and airflow towards the closest outlet. The film flux was set to 0.77g/s (approximately one-third of the value considered for the full bay), at which value the fuel film mass and its vaporization rate stabilized within several seconds. In Fig. 8 the change in time of film mass is shown for four mesh densities. Results indicate that while its mass stabilizes in time in each case, the stable values for denser meshes are lower than ones obtained for less dense meshes within the tested range. A possible explanation for this result is the increasing accuracy in determining the surface area covered by the film, achieved with smaller surface elements.

Fig. 8.

Figure 9 shows the change in vaporization rate for different mesh densities. As the mesh density increases, the maximum vaporization rate first increases, then decreases to a value close to that obtained for the smallest mesh. The changes are within approximately 22% of the maximum value, whereas the total mesh size changes by 40% of the maximum size. Figure 10 presents the mass of vapor collecting inside the bay as result of difference between the vaporization rate and outflow flux through the ventilation outlets. Except for the result for the roughest mesh, the maxima for the denser meshes are very close to each other and they occur at a similar time. For meshes bigger than 1.5 mln, after reaching maximum value the vapor mass stabilizes or tends to decrease.

Fig. 9.

Fig. 10.

The results for gas vapor collection, along with the vaporization rate and film mass development over time, indicate that for mesh densities higher than very low density used for the coarsest mesh, the investigated phenomena become qualitatively insensitive to mesh density. Most importantly, the mass of fuel vapor trapped inside the bay is relatively insensitive to mesh density. Based on results shown in Fig. 9 and Fig. 10 it was decided that mesh density corresponding to 2.36 million elements was sufficient for further simulations of the fuel vapor flow in the full volume of the equipment bay.

Regarding the sensitivity of the vapor parameters inside the equipment bay to the mesh density in the external domain, the most important potential effects concern the influence on the pressure difference between the inlets and outlets of the bay. A decreasing pressure difference may worsen the ventilation of the bay and may lead more gas vapor to be trapped inside it. Therefore, it was decided to skip simulations with varying external grid densities and to directly investigate the effect of a decreasing pressure difference driving the internal flow on the removal of fuel vapor from the bay. The reduced pressure difference was achieved by lowering pressure values at the boundary condition surfaces, which were reduced by 30% and 50% from their nominal values (164 Pa of total pressure at the inlets and -60 Pa of static pressure at the outlets). These changes in pressure difference were significant, and should exceed the effects of local pressure dependence on the external mesh density. The simulation results, shown in Fig. 11, indicate a rather moderate influence of the inlet-outlet pressure difference on the mass of the trapped gas. Closer inspection of this phenomenon revealed that the primary cause of this reduction was a decreasing vaporization rate, which accompanies the airflow at decreasing pressure difference between the inlets and outlets, as shown in Fig. 12. In Fig. 13, in turn, also demonstrates that the fuel vapor outflow rate is proportional to the vaporization rate.

Fig. 11.

Fig. 12.

Fig. 13.

The first stage of flow simulation for the full baseline bay involved achieving quasi-steady airflow conditions in the equipment bay after flow initiation, before initializing the DPM model and simulations of fuel film flow and vaporization. Flow simulations conducted for the case with a stationary airplane in the large domain (mesh types 1 and 2) produced the following values of pressure on inlet and outlet surfaces and air mass flow through the equipment bay:

The pressure values presented in Table 1 were determined on an “intermediate” model as surface-averaged values on the respective inlet and outlet surfaces and applied as boundary conditions for simulations of fuel leak and phase change. The leak source was set at a point located under a connection of elements chosen in consultation with the designers. It was modelled as an injection point of a stream of discrete particles using the Discrete Particle Model of ANSYS-Fluent. The droplets had a pre-defined diameter of 1 mm, and were subject to 1g acceleration, influenced by the airflow in the bay.

After quasi-steady airflow conditions were reached, the Discrete Particle Model was activated and droplets began to be injected into the domain. Upon reaching the floor surface, the droplets turned into an Eulerian Wall Film and moved along the floor of the bay. The transition of the droplets into the Eulerian Wall Film at 2s after the start of leak simulation is shown in Fig. 14. The fuel mass flow of the simulated leak was 2.6 g/s (approximately 2% of the electric pump’s capability).

Fig. 14.

The development over time of the area wetted by the fuel film is shown in Fig. 15. The fuel film can be seen to spread towards the sides of the equipment bay along the surface of the last bulkhead.

Fig. 15.

Results of the simulated, approximately 68 second-long flow of leaked fuel show steady growth of the mass of fuel trapped inside the equipment bay. The wetted area grows toward the surface of the last bulkhead, and then sidewise. The reason for this flow direction is visible in Fig. 16, where the profile of the bottom surface of the equipment bay is shown in the form of a colour map. It can be seen that: a) the outlets of the bay are located at a distance (approximately 5cm) in front of the last bulkhead, due to presence of a stiffener belt, approximately 2-3 mm thick (indicated also in Fig. 15), b) that the central bottom surface is sloped downwards in the direction of the fuselage length, and c) the central part of the bottom surface is several mm below the level of outlets located at its sides. This causes part of the leaked fuel to be trapped inside the bay, which then vaporizes. The growth over time of the mass of trapped fuel liquid and trapped fuel vapor is presented in Fig. 17.

Fig. 16.

Fig. 17.

The round holes in the bottom surface are enlarged versions of the drainage holes shown in Fig. 15, introduced in the later phase of the work.

Fig. 18.

The steady increase in the mass of gas vapor trapped inside the equipment bay leads to a rising concentration of vapor, which locally becomes higher than the flammability level (1.4% volume ratio, 5% mass ratio). The spatial distribution of the gas concentration was obtained using the cell register functionality in the ANSYS Fluent solver, which can be used to indicate cells with a chosen variable having prescribed value or being in a prescribed range. Fig. 19 shows two cell registers, with cells of mass concentration of fuel vapor above the flammability level in two time moments. The spatial distribution of these cells is consistent with the vector visualization of airflow in the vertical plane containing the inlet to the equipment bay, presented in Fig. 20.

Fig. 19.

Fig. 20.

It can be seen that – due to the location of the single ventilation inlet close to the last bulkhead, just below the cabin floor and in front of the “cavity step” in the floor – the air entering the bay almost immediately loses its energy and is directed downwards with low velocity, creating a slow-rotating vortex in the part of the bay just below the inlet, while in the other part of the bay the air is almost at a stand-still.

In order to minimize the likelihood of an on-board fire due to such a fuel leak, an improvement in ventilation and removal of liquid fuel from the bottom surface of the bay was sought. The measures to achieve this included:

- adding inlets to equipment bay,

The position of additional inlets to the equipment bay is shown in Fig. 21. The size of the inlets was limited by the distance between stringers and framers of the fuselage. Their position was selected such that the mutual interference between them should be minimal, and they should fulfil the additional task of providing cooling air for some additional electric subsystems in the equipment bay, not mentioned in detail in this paper. Each of the additional inlets had a scoop and a stream diverter, as shown in Fig. 22. All the inlets were located on the left side of the fuselage, due to presence of large exhaust pipe, disturbing the flow, on its right side.

Fig. 21.

Fig. 22.

In order to improve removal of leaking fuel and fuel vapor, additional outlets were proposed. The computational model was created with a large number of optional outlet surfaces, which could be switched on or off by a change in the type of boundary condition on each of the surfaces. The additional outlets, visible in Fig. 23 and in Fig. 24, were 20mm in diameter. This approach was chosen in order to determine the minimum number of outlet surfaces necessary for achieving proper ventilation by means of a change in the local boundary condition in the solver, on a surface prepared for this task, rather than by means of a more time-consuming modification of the CAD model and the creation of a new computational mesh.

Fig. 23.

Fig. 24.

First, the results concerning the additional functionality of providing cooling, fresh air towards selected sub-systems are presented. These computations were performed on an “intermediate” model, with the computational domain both inside the bay and outside the airplane, as shown in Fig. 24. The velocity contours in cross-sections, where these elements are present, are shown in Fig. 25 and Fig. 26, as contours of the Z velocity (horizontal, negative for air moving to the right side). It can be seen that the additional inlets effectively fulfil their role. In Fig. 26, the existing, inefficient inlet is also visible, just below cabin floor.

Fig. 25.

Fig. 26.

For the modified ventilation system the first run of simulations of airflow and flow of leaking fuel and its vapor was conducted with all the additional outlets active. The impact point of the droplets on the bottom surface was at the same distance from the plane of the last bulkhead, as in the first simulation for the baseline case. In Fig. 27, Fig. 28 and Fig. 29 it can be seen that the outflow of the leaking fuel and fuel vapor stabilized within a few seconds. Figure 30 presents the surface distribution of film thickness, while Fig. 31 presents spatial distribution of cells with vapor mass concentration exceeding the flammability level. It is evident that the leaking fuel is being carried away from the bay through the closest outlets, and the region of mass concentration of vapor is restricted to the very close neighbourhood of the outlets

Fig. 27.

Fig. 28.

Fig. 29.

Fig. 30.

Fig. 31.

Two more stringent simulation tests were next conducted. Since perforating the bottom surface of the equipment bay would obviously decrease its structural strength, in the first of the two tests the number of additional outlets was reduced to one in each skin strip between the longitudinal stringers. In the second test, additionally, the pressure boundary conditions corresponded to the lowest rotation velocity of the propeller, 481 rpm, occurring for approximately 30 seconds during engine start and warm-up (a decrease from 700 rpm). The approximately 15-second period of this phase was simulated, during which time the flow variables concerning the liquid fuel and vapor mass outflow, as well as the mass of fuel and vapor remaining inside the bay stabilized. The results are presented in Fig. 32 and Fig. 33. Figure 34 shows the volume cells with mass concentration above the flammability level (in red), and cells with vapor mass concentration starting from the lower value of 50% of flammability level. The visualisations were obtained for propeller rpm=481. Comparison of the two concentration patterns shows the direction where the vapor cloud may expand in the event of an increasing volume of leaking fuel. This direction is along the wall, towards other outlets.

Fig. 32.

Fig. 33.

Fig. 34.

Fig. 35.

The simulations of the flow and phase change in the equipment bay indicate that the modification most crucially needed is the elimination of zones of where liquid fuel originating from leaks in connectors of the elements of the fuel system may collect and vaporize. This may be achieved by creating sloped fillings above the strengthening belts of the bulkheads and adding additional openings in the floor, shielded by external fairings against backflow of air. These modifications are relatively easy to fabricate, however, they might necessitate additional considerations regarding the structural strength of the surfaces involved. Additional modifications include air inlets which could direct the airflow inside the bay towards the proper outlets in the floor, preventing the circulation of a mixture of air and fuel vapor inside the bay. The presence of longitudinal stringers in the floor of the bay favours fuel film flow towards the outlets. Additionally, the fact that fuel vapor is heavier than air aids in its removal through the proposed outlets.

The proposed structural modifications need to have their effectiveness verified in real working conditions. For safety reasons, this may be carried out without spilling inflammable fluids. The flow pattern of leaking fluid may be visualized with dyed water, and this technique should verify the elimination of places where fluid collects before exiting the bay, while air velocity can be effectively measured by sensors, such as thermo-anemometric sensors, situated at points of interest. Much as for the liquid fluid flow pattern, the flowpaths of gas in the bay can be verified by a smoke-visualisation technique.

Regarding the simulation techniques, the combination of the Discrete Particle Model with Eulerian Wall Film and Phase Change modelling, as used herein, has been shown to offer a convenient and effective solution for modelling of leaks of arbitrary origin, which may be easily changed in space without modification of the computational mesh. The wall film model only operates on surface elements of the mesh and is therefore a very effective modelling technique, faster than Volume of Fluid, which may serve as an alternative, but requires a locally dense volume mesh to resolve the height of the layer of the liquid and its flow. However, the development of the film flow and vapor flow over time is a process that consumes computational time and memory, which restricted the number of cases analysed.