Travelling Waves Approach in a Parabolic Coupled System for Modelling the Behaviour of Substances in a Fuel Tank

: The aim of this work was to provide a formulation of a non-linear diffusion model with forced convection in the form of a reaction–absorption system. The model was studied with analytical and numerical approaches in the frame of the parabolic operators theory. In addition, the solutions are applied to a gas interaction phenomenon with the intention of producing an inerted ullage in an Airbus A320 aircraft centre fuel tank. We made use of the travelling wave (TW) solutions approach to study the existence of solutions, stability and the precise evolution of proﬁles. The application exercise sought to answer a key question for aerospace sciences which can be formulated as the time required to ensure an aircraft fuel tank is safe (inerted) to prevent explosion due to the presence of oxygen in the tank ullage.


Introduction
In the 1930s, Fisher [1], proposed a reaction-diffusion model to describe the interactive dynamic of genes. In parallel, Kolmogorov, Petrovskii and Piskunov [2] proposed the same equation in combustion theory. In both cases, the models considered a Gaussian order two-diffusion model with a non-linear reaction of the form f (u) = u(1 − u). The authors introduced the concept of travelling wave (TW) solutions to describe the propagation features of each involved species. Afterwards, the Fisher-KPP model was subjected to remarkable mathematical research to explore further applications (see [3][4][5]). More recently, some analyses [6] have shown new patterns of formation in chemistry and biology (compared to those existing in the current literature, see the remarkable references [7,8]) for travelling waves solutions and making use of numerical algorithms.
The aim of this paper was to define a model and obtain results to predict the behaviour of interacting gases in a global medium with forced convection (advection) and diffusion. The gas interacting process, aimed to model, consists of the removal of oxygen from a fuel tank to prevent any hydrocarbon combustion in the case a generated spark or hot spot occurs [9]. One of the technical solutions raised for this purpose is known as an inerting system, which shall be understood as a fire preventive [10]. As a short description, the inerting system consists of a filter or membrane that separates the nitrogen and oxygen in the air. The nitrogen is introduced into the fuel tank while the oxygen is expelled outwards. The inerted atmosphere prevents the fire propagation as there is not sufficient oxygen concentration to produce and sustain a combustion process. One possibility to model the nitrogen-oxygen interaction is to consider the nitrogen replacement by oxygen in the fuel tank air space. In [11], the authors developed an algebraic model based on a mass balance to determine the oxygen and nitrogen concentrations in fuel tank bays. Additionally, in [12], the algebraic model was compared with a transport differential equation based on a mass balance in a differential tank element. In both cases, the modelling exercise does not introduce diffusion into the interacting species, although the forced convection is introduced by the transport (in the sense of particle movement) phenomenon modelled.
The diffusion in the fuel tank was introduced by Ghadirian et al. [13] for simulating the interactive process of fuel vapours in a fuel tank. Such diffusion was proposed to be governed by the classical parabolic homogeneous equation: where C f is the fuel concentration in the air, and z is the vertical coordinate. This problem considers only the diffusion, disregarding the natural convection terms and reaction or absorption phenomenon. The intention of this paper was to consider the forced convection together with a proposed interaction between the substances in the fuel tank.

Materials and Methods
The methods followed along this research were based on a hybrid assessment, analytical and numerical, together with a validation exercise with real testing data. The modelling set of equations were obtained based on the mass conservation principles.
The nitrogen gas was supplied by a manifold with several nozzles located along the fuel tank (see [14] for a complete description). The gases were dynamic together with the nitrogen nozzle location and the particular tank geometries led to an heterogeneous gas distribution along the tank. This heterogeneous condition has been modelled considering the stratified two-phase flows in [15], making use of the volume of fluid (VOF) model (see [16,17] for a detailed description) by solving the phase continuity equation. In the modelling exercise developed in this paper, it was considered that the heterogeneous nitrogen discharge in complex tank geometries provoked gas dispersion. As a consequence, a diffusion process emerges in response to the differences in the gas concentration along the tank zones. The selection of an appropriate diffusion principle is of relevance and leads to a whole significant discussion (see [18] and the references therein). For our purposes, the considered set of equations is non-linear in general and aims to characterize the propagation of the interacting gases. Such interactions are of significance in tank zones not influenced by the advection, so that the nitrogen propagates by diffusion. This progressive propagation in a diffusive front has been used for modelling purposes in other fields to simulate porosity effects [19,20]. In our case, the propagation dynamic was conceived as an inherent feature introduced by the travelling waves solutions.
In addition to the diffusion, the rate of nitrogen discharge along the domain induces a forced convection within the tank airspace.
To obtain the modelling equations, we consider the following dynamics between the nitrogen concentration (N) and the oxygen concentration (Θ): initially, the fuel tank ullage is not saturated with nitrogen, so it is considered that the nitrogen time rate (N t ) is high. As long as the nitrogen concentration increases, the oxygen is expelled and the nitrogen time rate decreases. Then: where n ∈ (0, 1) shall be obtained for each particular case. In the same manner, the oxygen concentration decreases for increasing nitrogen quantities: The expressions (2) and (3) will be referred to as the reaction (Θ n ) and absorption (−N m ) in the modelling set of equations Note that both the nitrogen and oxygen are functions of space and time N(x, t) and Θ(x, t) for 0 < t < ∞, x ∈ R 3 (d = 3, in this case).
Consider a virtual domain within the fuel tank Σ ∈ R 3 . The mass conservation principle establishes that the time rate of nitrogen and oxygen concentrations (N t , Θ t ) shall be equal to the concentration flux through the boundaries of Σ (∂Σ) plus the absorption or reaction that act as independent terms.
We admit that the vectors Ψ N ∈ R 3 and Ψ Θ ∈ R 3 are the concentration fluxes along the boundary (∂Σ) per unit of area and per unit of time through the borders of the subdomain Σ. Consider the vector π ∈ R 3 as the normal local at any point of ∂Σ. We admit that the absorption and reaction terms per unit of time and volume are given by the generic functions G N (N, Θ, x, t) for the nitrogen and G Θ (N, Θ, x, t) for the oxygen.
Each Ψ N , Ψ Θ , G N , G Θ is continuous with continuous derivatives. Therefore, for sufficiently continuous initial data, i.e., in the domain, the evolving solutions N(x, t) and Θ(x, t) are continuous and smooth (Chapter 7 [21]) in virtue of the regular parabolic spatial operator that adopts the Laplacian form as it will be shown. Then: The Fick law [21] enables us to account for an expression to the vector fluxes: The diffusion between the two interacting species requires that ρ = σ. Note that the vector a ∈ R 3 introduces the advection term to account for the forced convection within the fuel tank. Hence: Note that: So that: The initial concentrations of nitrogen and oxygen are mathematically given by constant Considering the above information, the following problem P is formulated: The operator's parabolic regularity ensures the monotone evolution of each of the involved solutions (referring to Chapter 7 [21]). The decreasing and increasing condition of the solutions are, then, given by the reaction (increasing) and absorption terms (decreasing). As a consequence and preliminary, the nitrogen concentration N increases while the oxygen concentration Θ decreases. Then, and upon evolution, the solutions will satisfy N(x, t), Θ(x, t) ∈ L 1 (R 3 ), which means that the total concentration of each substance does not blow up in the domain due to the non-linearity proposed in the reaction and absorption terms.
The problem P was introduced as a baseline model to account for the different terms involved in the gases dynamics. Now, a new problem is defined accordingly to account for the asymptotic behaviour of each of the substances. As it will be shown in Section 3.3, the fuel tank inerting process tends towards stationary values for the oxygen and nitrogen concentrations. This fact leads to consider travelling waves (TW) solutions connecting two stationary concentrations, both given by the initial and the asymptotic values. Mathematically, such solutions can be modelled by a Fisher-KPP reaction type. Consequently, the new problem in the TW domain P T is: where r ≥ max{N 0 (x)} is the asymptotic solution to N. If the initial data are constant, the same condition is expressed as r ≥ N 0 . The kinematic boundary condition related to the non-slip suggests that the velocity variations are located closed to the tank walls, while far from such borders, the described dynamic behaves free with minor influence of the boundary conditions. In addition, a dedicated description of the involved gases kinematic is out of the analysis scope; nonetheless, the influence and behaviour of the convective effects (through the vector a) are considered, and particularly, so is the effect of such advection into the gases concentrations (u, v). Thus, the intention is to study the dynamic associated with the solutions at the core of the domain (sufficiently large as to consider R 3 ) not affected by the domain borders.

Analysis in the Travelling Waves (TW) Domain
Solutions to problem P T (15) tend to equilibrium solutions upon evolution. Such equilibrium conditions are given by the technological limits and other environmental variables (see Section 3.3). In addition, note that the increasing monotone behaviour in the nitrogen leads to N t > 0 while Θ t < 0 for the oxygen that is expelled out of the tank. Then, the equilibrium condition for the nitrogen is expressed as The oxygen concentration decreases so as to reach an asymptotic value close to the null equilibrium: To understand the involved dynamic and permit the analysis of TW evolutions, the initial conditions are given by a step function of Heaviside type: The main objective is to determine the nitrogen and oxygen concentrations starting from the given initial data and evolving to the equilibrium solutions in the TW domain. Note that the step-like initial data permit us to analyse the behaviour for positive and null initial masses, as both are defined in the x intervals of a Heaviside function, i.e., H(−x) = 0 for x ∈ (0, ∞) and H(−x) = 1 for x ∈ (−∞, 0). Along this section, solutions are obtained making use of hybrid analytical and numerical evidence.

TW Profiles
where d is a unitary vector in the TW direction, λ corresponds to the TW-speed and ϕ 1 ∈ L ∞ (R 3 ) ∩ L 1 loc (R 3 ) is the TW profile. Note that TW profiles are equivalent under translation ξ → ξ + ξ 1 and symmetry ξ → −ξ. We admit the vector in the TW propagation direction is expressed as d = (1, 0, 0) ∈ R 3 and the TW moves from −∞ to ∞. Hence: After replacement and standard operations, the problem (15) reads: Note that the TW profile ϕ 2 asymptotically tends towards zero whilst keeping the positivity condition. In addition, ϕ 1 increases to r > H(−x), so that (r − ϕ 1 ) > 0. These particular behaviours will be shown afterwards, but shall be taken into account during the coming assessments.
A linearisation of (19) close to the equilibrium condition (ϕ 1 = r, ϕ 2 = 0) is difficult in general due to the non-Lipschitz character of ϕ n 2 . Note that close to the critical point ϕ 2 = 0, ϕ n 2 > ϕ 2 , hence, the following alternative problem is formulated: (19): TW profiles resulting from (20) are lower solutions to (19). This statement will be shown afterwards, but it is relevant to ensure the positivity of any TW profile. Indeed, any positive solution to (20) will provide evidence to ensure the positivity of any upper solution in (19).
Considering r − ϕ 1 =φ 1 , then: With the boundary conditions expressing the asymptotic behaviour of solutions: The fundamental TW problem consist of finding a positive and monotone profile for a certain interval of TW speeds. Particularly, solutions are expressed in the proximity of the equilibrium points ϕ 1 = 0 and ϕ 2 = 0. When solutions approach the boundary conditions in (22), a linearisation approach can be followed, so that solutions can be exponentially expressed asφ From the second identity, it is possible to find an expression to relate the given speeds (TW speed and advection terms) with the exponential factor: Similarly, for γ 1 : Note that for ξ >> 1, the term 4Ee −γ 2 ξ → 0 + , then: Observe that the results for γ 1 and γ 2 establish that the exponential decaying terms are bigger than the advection and TW-speed terms (λ + a). Hence, any TW profiles evolves towards the asymptotic equilibrium and this evolution is convergent, independently of the advection and TW speed terms.
Note that the existence, uniqueness and comparison of solutions to problems (19) and (20) follow a standard procedure, as proposed in [21] (Chapter 7). Along the following analysis, a numerical exercise is considered to provide particular forms of TW solutions to problem P T (19). Furthermore, solutions to (20), close to the critical points (ϕ 1 = r, ϕ 2 = 0), are effectively shown to be lower solutions to (19). As solutions to problem (20) are lower compared with solutions to (19), a monotone (increasing or decreasing) TW profile and tip exist. The numerical exercise has been pursued with the following features: • The methodology is based on the Matlab function bvp4c which provides an implicit Runge-Kutta approach with an interpolant extension [22]. The collocation method requires the specification of boundary conditions, in this case, given by the stationary solutions at −∞ and +∞; • The integration interval is (0, 100), large enough to study the TW evolution in their domain, and not impacted by the boundary conditions that act as tractors; • The considered error for computation is 10 −6 ; • The integration domain has been divided into 100,000 nodes; • To make the numerical assessment tractable and without loss of generality, particular values were taken for the involved parameters: r = 2, (a + λ) = 4-and different values of n and m.
Results are represented in Figures 1-5 with the following remarks: • Case with n = 0.9 and m = 0.1-according to the results in Figures 1-3, TW profiles to (20) are lower solutions compared to those for (19). Note that ϕ 2m (as solution to (20)) is positive along the domain. Then, any TW moving with the exponential decaying term γ 2 as per (25) is positive in the whole domain. In addition, we note that the solutions to (19) are very close to solutions to the linearised (20), which permits the validation of the goodness of the linearised exercise; • Similarly, the same exercise is repeated for n = 0.2 and m = 0.9. TW profiles, as solutions to (20), are lower compared with solutions to (19) (refer to Figures 4 and 5). In the same way, note that solutions to (19) evolve close to solutions for the linearised (20). The numerical exercise was performed for an integration interval (0, 100); nonetheless, the interval was cut for representation purposes. Note that f 1 and f 2 represent ϕ 1 and ϕ 2 solutions to (19) while f 1m and f 2m are solutions to (20) .
. Figure 2. TW profiles for a + λ = 4 with n = 0.9, m = 0.1. For the sake of simplicity, r = 2. Note the lower solution ϕ 1m as the solution for the problem (20). The numerical exercise was performed for an integration interval (0, 100); nonetheless, the interval was cut for representation purposes. Note that f 1 and f 2 represent ϕ 1 and ϕ 2 solutions to (19) while f 1m and f 2m are solutions to (20).  (20). The numerical exercise was performed for an integration interval (0, 100); nonetheless, the interval was cut for representation purposes. Note that f 1 and f 2 represent ϕ 1 and ϕ 2 solutions to (19) while f 1m and f 2m are solutions to (20).  To summarize, the TW-speed expressions in (25) and (26) provide positive lower solutions to the problem (19). This result permits us to ensure that any TW profile to the original problem (15) is positive upon evolution to the equilibrium conditions N = r and Θ = 0.

Application to a Fuel Tank Inerting Process
The phenomenon under study is related to the interaction between substances, i.e., the oxygen and the nitrogen in an aircraft fuel tank to generate an inert gas ullage. The model in Problem P has three different parameters (n, m, a) that need to be determined, making use of real flight testing activities. Such flight test data are extracted from [14] and summarized in Figure 6. Note that for the particular oxygen sampling location, we refer the reader to the discussion in the sourced reference [14].
The diffusion coefficients between the interacting substances can be considered as constant during the flight evolution: which corresponds to a flight profile at a total environmental (static and dynamic) temperature of −7 ºC. Note that the diffusion coefficients can be determined at any temperature [23]: where ρ T 1 and ρ T 2 are the diffusion coefficients at temperatures T 1 and T 2 expressed in K.
The calibration exercise to assess the parameters (n, m, a) requires the measurement of the nitrogen (N) or the oxygen concentration (Θ). Aiming to calibrate on (n, m) and for the sake of simplicity, the air is assumed to be an homogeneous mixture formed of nitrogen and oxygen. The flight test data measure the oxygen concentration in the centre tank of an Airbus A320 [14]. Figure 6 considers the case of an empty tank and one nitrogen gas generator acting (single membrane). . Note that the increase in oxygen concentration at t = 78 min is due to air ingestion in the fuel tanks during the descent phase (source [14]).
An empty tank is defined as that with low level of fuel and high level of fuel vapours. It is usual that the inerting system is constituted by several inert gas generators or air separation modules (ASM). This is important to increase the inerting capacity and reduce risks. An empty tank has a higher level of oxygen that shall be replaced by nitrogen to provide an inert condition. In addition, a test case with a single acting ASM is considered, which means that only one ASM is available to filter the nitrogen and remove the oxygen. This is conservative and leads to the most demanding configuration for the single inerting ASM which, in turn, will lead to provide higher lead times to get the inert condition.
The flight test data compiled in Figure 6 are interpreted over homogeneous mean oxygen values to define the parameters n and m. Hence, the flat problem P is: For each time step, the time derivatives N t and Θ t are determined. For this purpose, it is considered that the substance concentrations are expressed between the interval (0, 1) and the time is expressed in minutes. Such time derivatives are obtained for different discrete times and represented versus the concentration of each substance. Afterwards, the optimal curve fitting is determined as Therefore, n = 0.586. The same process for Θ t provides: Then, m = 0.025. The parameters obtained, based on flight test data, are indeed between (0, 1), as initially appointed.
Then, the forced convective parameter a is obtained. For this purpose, the basic continuity equation is used: where G is the inerted gas flow, rich in nitrogen, and D is the air density within the fuel tank, mainly constituted of nitrogen in stabilized flight. Note that the nitrogen flow is kept stable during the flight and increases during the descend phase. For the calibration exercise, it is considered the flight cruising phase (at t = 40 min), in which the nitrogen flow value is given by Figure 7 The A320 centre fuel tank is almost of rectangular shape with a mean height of 1.5 m and a mean wide of 6 m: A T = 9 m 2 . (36) Then, a can be determined as per (33) a = 0.0125 m/min. After calibration, the problem P is summarized as Once the model parameters have been obtained, TW solutions are provided for the initial concentrations given by the level of nitrogen and oxygen in the standard atmosphere: The TW solutions are formed of a front and a tip. The front carries the information transition in the media from one state to the other. This kind of behaviour exhibited can be applied to the propagation of the nitrogen gas within the fuel tank. The nitrogen TW constitutes a propagating front that shifts the tank ullage concentration by reducing the oxygen and increasing the nitrogen.
Note that there exist stationary solutions given by the initial conditions and the asymptotic approximation. Such asymptotic solutions are provided by the levels of concentration admitted by the particular technology employed (see [10] for a description). For this pur-pose, it is considered that t = 78 min (according to Figure 6 and previously to starting the descent phase) is sufficiently large time for the stationary solutions: The problem P T , analysed in the travelling wave domain, adopts the form: where the parameters n, m, ρ, σ and a are given in (37). In this case, a corrected speed λ * = λ − a is introduced to account for the advection term. Then: To determine the TW propagation speed, we first consider that the TW front and tip move with corrected speed λ * = λ − a, so that when the variable ξ = 0, the wave has propagated along the whole domain: Once the TW has finished the propagation along the tank, the inerted level corresponding to N = 0.98 is reached. At t = 80 min, the TW has reached all tank areas to provide the lowest value, measured in flight, for the oxygen concentration. In addition, x is a typical tank dimension obtained as the geometric mean of the tank shape: After assessing the TW speed, the model (40) is solved with the function bvp4c in Matlab, as shown in Section 3.2. The results, in the TW variable ξ, are compiled in Figure 8. . Travelling wave solutions. The blue solution represents the evolution of the nitrogen within the tank evolving from 0.8 (80% in volume) up to filling the whole tank 1 (100%). The red solution represents the oxygen density that evolves from 0.2 (20% in volume) down to zero. Note that the horizontal axis represents the variable ξ, while ϕ 1 = f 1 = N and ϕ 2 = f 2 = Θ.
The time required to reach a certain level of nitrogen N that ensures an inerted condition is then obtained. As an example, a concentration of 9% of oxygen in a fuel tank is sufficient to get the inert status [24]. A value of Θ = 0.09 corresponds to the following value of ξ ( Figure 8): To determine the time to inert, we consider that the advection term acts physically in the same direction as the TW-speed λ. This assumption is realistic as the forced convection drives the nitrogen to all tank areas where the gas is free to flow. Nonetheless, once the nitrogen is in complicated tank shapes, it moves by diffusion only. Hence, the forced convection (a) and the diffusive TW (λ) are summative terms that act to reach the inert condition in all tank areas (both free paths and difficult geometries). Thus, the time to obtain the inerted configuration under similar operational conditions compared to the A320 aircraft is: Globally, at t = 52.7 min, the tank will reach a value of 9% of oxygen in all areas. This value is overly pessimistic compared to the value provided in Figure 6, in which the value of 9% of oxygen is obtained at t = 44 min. Nonetheless, the values of oxygen in Figure 6 are only provided for different locations where the oxygen sensors are placed. Therefore, the modelled value of t = 52.7 min, shall be seen as an approach in which the diffusion acts as per the travelling wave front and tip to cover the whole domain.
It is particularly relevant to observe that the scale for the tank dimension (x = 3.88 m) is much higher than the scale of the TW spatial variable ξ. This fact can be explained in view of the TW features. The front and the tip represent the diffusion acting along the wave so that ξ represents a dimension of the wave interface where the front and the tip are confined.

Conclusions
The problem P T proposed with a coupled system of reaction diffusion equations to model an aircraft fuel tank inerting process has been discussed with a mathematical approach stressing aspects related with the existence, uniqueness and behaviour of solutions in the travelling wave domain. The application exercise to a centre fuel tank on an Airbus A320 set the evidence for the use of the problem P T to model inerting processes in aircraft. The information provided permitted the mathematical representation of the oxygen and nitrogen concentrations and enabled our understanding of the dynamic of the process based on a parabolic coupled system of equations.
In addition, finite values for the model parameters n, m and a have been shown to exist, and furthermore, the combination of such values has been shown to provide global solutions in the TW frame. Eventually, the model has been shown to provide an answer to a basic question asked by the engineering and physical application, i.e., the global time required to ensure a fuel tank is safe (inerted) to prevent explosion due to the presence of oxygen in the tank ullage.
Funding: This research received no external funding.