Mariana Vargas

Robust Estimation of Laplacian Constrained Gaussian Graphical Models with Trimmed Non-convex Regularization

The problem of discovering a structure that fits a collection of vector data is of crucial importance for a variety of applications. Such problems can be framed as Laplacian constrained Gaussian Graphical Model inference. Existing algorithms rely on the assumption that all the available observations are drawn from the same Multivariate Gaussian distribution. However, in practice it is common to find scenarios where the datasets are contaminated with a certain number of outliers. The purpose of this work is to address that problem. We propose a robust method based on Trimmed Least Squares that copes with the presence of corrupted samples. We provide statistical guarantees on the estimation error and present results on simulated data.

Introduction

Numerous modern applications rely on relational information to study a particular phenomenon or address a specific task. For example, in social network analysis, predictions about the behavior of a particular user may also depend on the behavior of its mutual friends. The need to incorporate relational information into Machine Learning algorithms calls for data processing methods that can handle highly correlated data. In this context, graphs have been widely used because of their ability to compactly encode pairwise relationships between points. While in many domains the graph structure arises naturally, in many scenarios it is often the case that no a-priori graph is readily available. The problem of estimating the graph structure that fits a collection of data can be framed as Laplacian Constrained Gaussian Graphical Model (LCGGM) inference.

In its general form, a Gaussian Graphical Model (GGM) is parameterized by a matrix called the precision matrix that encodes the conditional dependence relationships between 1 Inria Grenoble, France. Correspondence to: Mariana Vargas Vieyra <mariana.vargas-vieyra@inria.fr>.

Presented at ICML Workshop on Principles of Distribution Shift 2022,Copyright 2022 by the author(s).

variables. Most methods for estimating GGMs impose an L 1 penalization on the parameters to encourage sparisity. This is known as Graphical Lasso [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Banerjee | Sparse covariance selection via robust maximum likelihood estimation[END_REF][START_REF] Cai | Estimating sparse precision matrix: Optimal rates of convergence and adaptive estimation[END_REF]. When constrained to graph Laplacian matrices, the precision matrix of the GGM can be deemed an operator with specific properties that permits to reason about the observed data. In particular, the eigenvectors of a Laplacian matrix provide a Fourier basis of the data, and the eigenvalues can be interpreted as spectral frequencies [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF]. The study of graph Laplacians have led to the field of Graph Signal Processing [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Dong | Graph signal processing for machine learning: A review and new perspectives[END_REF].

In general, methods for LCGGM estimation encourage sparse estimators [START_REF] Ying | Nonconvex sparse graph learning under laplacian constrained graphical model[END_REF]2021;[START_REF] Lake | Discovering structure by learning sparse graphs[END_REF][START_REF] Koyakumaru | Learning sparse graph with minimax concave penalty under gaussian markov random fields[END_REF], although guarantees for non-regularized Maximum Likelihood Estimators have been explored [START_REF] Egilmez | Graph learning from data under laplacian and structural constraints[END_REF][START_REF] Ying | Minimax estimation of laplacian constrained precision matrices[END_REF][START_REF] Pavez | Laplacian constrained precision matrix estimation: Existence and high dimensional consistency[END_REF]. Other algorithms take a different route and propose estimators based on global smoothness [START_REF] Dong | Learning laplacian matrix in smooth graph signal representations[END_REF][START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF]. Finally, algorithms that allow more restrictive constraints have been studied [START_REF] Kumar | Structured graph learning via laplacian spectral constraints[END_REF]. These approaches rely on the hypothesis that all observed elements come from the same distribution. However, in practice this is often a strong assumption, as available samples might be contaminated with outliers. While robust methods have been proposed in the context of vanilla GGMs [START_REF] Yang | Robust gaussian graphical modeling with the trimmed graphical lasso[END_REF], the problem of estimating LCGGM with corrupted samples remains under explored.

In this work we propose a robust method based on Trimmed Least Squares to estimate sparse LCGGMs in scenarios where the available dataset has a relatively small proportion of outliers. We provide statistical guarantees on the estimation error and empirically validate our method on simulated data.

Background and Setting

We start by introducing the necessary notation and background on Laplacian matrices and Laplacian-constrained Graphical Models. We then describe the problem setup.

Graph Laplacian

Let us consider a graph G = (V, E, W ) where V is a set of p nodes, E ⊆ V × V is a set of edges, and W is a matrix with non-negative entries such that W ij > 0 if and only if (v i , v j ) ∈ E. The graph Laplacian is defined as L = D -W , where D is a diagonal matrix with elements D ii = p j=1 W ij . If we consider graph Laplacians associated to graphs with one connected component, then the set of valid graph Laplacians can be defined as:

S L = {Θ : Θ ij ≤ 0 if i = j, Θ1 = 0} (1)
The graph Laplacian is a crucial tool to reason about the structure of the graph and the data associated to it as they allow to generalize notions of Signal Processing to non-Euclidean, irregular domains [START_REF] Bronstein | Geometric deep learning: Going beyond euclidean data[END_REF]. In particular, graph Laplacians provide a Fourier basis of the associated data through their spectral decomposition. They also provide a quadratic form that quantifies the smoothness of the data. Global smoothness has been widely used as a criterion to assess the quality of an estimated graph [START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF][START_REF] Dong | Learning laplacian matrix in smooth graph signal representations[END_REF][START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF].

Problem setup

Let X = [X 1 , . . . , X p ] be a random vector with a Gaussian distribution, X ∼ N (0, Σ), and let X = [x T 1 , . . . , x

T n ] be a design matrix consisting of n p-dimensional observations of X. The precision matrix Θ = Σ -1 characterizes the conditional independence structure of the variables. That is, X is associated with a graph G where the nodes correspond to the elements of X and Θ ij > 0 if and only if X i and X j are connected in the graph. Let S = 1 n X T X be the sample covariance matrix obtained from the data, and let J be a constant matrix defined such that J ij = 1/p for all i, j. We define the linear operator L : R p(p-1)/2 → R p×p that maps vectors to Laplacian matrices introduced by [START_REF] Kumar | Structured graph learning via laplacian spectral constraints[END_REF] as

(Lw) ij =      -w k i > j (Lw) ji i < j -j =i (Lw) ij i = j
where k = i -j + j-1 2 (2p -j). L has an adjoint operator L * : R p×p → R p(p-1)/2 defined as

(L * W ) i = W ii -W ij -W ji + W jj .
In this work we adopt the setting for learning Gaussian Graphical Models under Laplacian constraints proposed by [START_REF] Ying | Nonconvex sparse graph learning under laplacian constrained graphical model[END_REF]. That is,

min w≥0 -log det(Lw + J) + Tr(LwS) + i r λ (w i ) (2)
where i r λ (w i ) is a nonconvex regularization function that encourages sparsity, and where we write w ≥ 0 for element-wise comparison. Note that by adding J in the first term the argument of the log determinant becomes full rank [START_REF] Egilmez | Graph learning from data under laplacian and structural constraints[END_REF].

We now assume that our data set is contaminated with samples drawn from a different distribution. This way, "good" samples come from the distribution with true parameter Θ and "bad" samples come from a different distribution. Let

[n] = {1, . . . , n}. Let G = {i : X (i) ∼ N (0, (Θ ) -1 )} be the set of indexes of "good" samples and B = [n] \ G the set of indexes of "bad" samples. We further denote g = |G| and b = |B| the number of "good" and "bad" samples respectively.

In the following we will introduce a method for robust estimation of Laplacian constrained graphical models and derive statistical guarantees on the estimation error.

Robust Estimation of Laplacian Constrained Gaussian Graphical Models

In this Section we first introduce a method inspired on Trimmed Least Squares for robust maximum likelihood estimation of Laplacian constrained graphical models with a nonconvex penalization. We then derive statistical guarantees for the estimation error.

Proposed Method

Following the same route as [START_REF] Yang | Robust gaussian graphical modeling with the trimmed graphical lasso[END_REF], we propose to solve the following objective function

min w≥0,h -log det(Lw + J) + Tr(LwS h ) + i r λ (wi) s.t. h T 1 = ĝ hi ∈ {0, 1} for all i = 1, . . . , n Lw 1 ≤ R (3) where S h = 1 ĝ i h i X (i) X (i) T
, h is a vector such that h i ∈ {0, 1} and h T 1 = ĝ, and ĝ is a hyperparameter that indicates our prior knowledge about how many "good" items we have in our dataset. Observe that the constraint Lw 1 ≤ R is rather mild in the sense that we can pick a large enough R so that w is feasible [START_REF] Loh | Regularized m-estimators with nonconvexity: statistical and algorithmic theory for local optima[END_REF].

The objective function (3) is biconvex. We can therefore alternate between w and h by fixing h and fitting w, and vice versa. Provided a fixed h, we optimize for w using the Majorization-Minimization (MM) framework [START_REF] Sun | Majorizationminimization algorithms in signal processing, communications, and machine learning[END_REF] as proposed by [START_REF] Ying | Nonconvex sparse graph learning under laplacian constrained graphical model[END_REF]. More specifically, we first define a surrogate function f k (w) = f (w|w (k) ) that approximates the objective function of Equation (3) at the point w (k) . Let z

(k) i = r λ (w (k)
i ) be the derivative of r λ at w

Input: w (0) , T , ĝ Initialize h (0) S h (0) ← 1 ĝ i h (0) i X (i) X (i) T repeat z (k) ← r λ (w (k) i ) for all i = 1, . . . , p(p -1)/2 ∇f (w (k) ) ← -L * (Lw (k) + J) -1 + L * S h (k) + z (k) Choose η (k) with line search. w (k+1) ← w (k) -η (k) ∇f (w (k) ) + if k mod T = 0 then Update h (k) S h (k) ← 1 ĝ i h (k) i X (i) X (i) T else h (k) ← h (k-1) {h (k) remains the same} S h (k) ← S h (k-1) {S h (k)
remains the same} end if until Stopping criterion satisfied earizing the regularization function i r λ (w i ) as follows:

f k (w) = -log det(Lw + J) + Tr(LwS h ) + i z (k) i w i .
We then update w (k) by

w (k+1) = w (k) -η (k) ∇f k (w (k) ) +
where [a] + = max(0, a), and where η (k) is the step size at the k th iteration that can be found with line search. Updating h when w is fixed amounts to setting h i = 1 if X (i) is among the ĝ elements with the highest log-likelihood and h i = 0 otherwise. In practice we can update h every T iterations until convergence, where T is a user input hyperparameter. The full procedure is described in Algorithm 1.

Choice of regularization function. In this work we restrict our analysis to two nonconvex type of penalties, namely MCP [START_REF] Zhang | Nearly unbiased variable selection under minimax concave penalty[END_REF] and SCAD [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]. These functions are defined as

r λ;MCP (t) = λ -t β if 0 ≤ t ≤ βλ 0 if t ≥ βλ and r λ;SCAD (t) =      λ if 0 ≤ t ≤ λ (βλ-t) β-1 if λ ≤ t ≤ βλ 0 if t ≥ βλ
However our analysis can be extended to other penalization functions as long as they fulfill a set of conditions related to their unbiasedness and sparsity promoting nature [START_REF] Ying | Nonconvex sparse graph learning under laplacian constrained graphical model[END_REF][START_REF] Loh | Regularized m-estimators with nonconvexity: statistical and algorithmic theory for local optima[END_REF]2017).

Theoretical Analysis

Optimization of Equation (3) yields local optima due to the non-convexity of the objective function. In what follows we will prove that arbitrary local optima are guaranteed to be consistent under the Frobenius norm, provided some assumptions on the data and local optima are fulfilled. We first present in detail these assumptions.

The following lemma is often refered to as restricted strong convexity and it is crucial for deriving the error bound presented in this section.

Lemma 3.1. (Ying et al. ( 2020)) Let w and ŵ be such that w, ŵ ≥ 0, Lw + J and L ŵ + J are positive definite, and Lw -L ŵ F ≤ r. Then the following holds,

-L * (Lw + J) -1 + L * (L ŵ + J) -1 , w -ŵ ≥ ( Lw * 2 + r) -2 Lw -L ŵ 2 F
Let κ = ( Lw * 2 + 1) -2 . We further assume the following for the coefficient λ. Assumption 3.2. The following holds for the regularization coefficient λ,

4 max{ L * (Lw +J) -1 -S h ∞ , τ1} ≤ λ ≤ κ- τ2 + λk √ 2 1 R .
The following lemma allows us to use Lemma 3.1 for arbitrary local optima of Equation (3).

Lemma 3.3. Let ( w, h) be any local optimum of Equation (3). Then, under Assumption 3.2 we have that

L w -Lw F ≤ 1
The proof can be found in Appendix C.

We need to make further assumptions about local optima and the dataset in order to handle scenarios where the data is contaminated with outliers. The following condition is known as mutual incoherence condition: Assumption 3.4. Let ( w, h) be any local optimum of Equation (3). Then,

Tr{(S h -S h )(L w-Lw )} ≤ τ1 w-w 1 +τ2 w-w 2 .
We finally assume a bound on the outliers and on the number of outliers with respect to the number of "good" samples.

Assumption 3.5. Let X B ∈ R b×p be the design matrix formed by the outliers. We assume X B 2 2 ≤ f (X B ) √ g log p for some function f [START_REF] Yang | Robust gaussian graphical modeling with the trimmed graphical lasso[END_REF].

Assumption 3.6. The ratio between the amount of outliers and "good" samples is such that b g ≤ 1 √ n .

Lemma 3.7. [START_REF] Nasrabadi | Robust lasso with missing and grossly corrupted observations[END_REF][START_REF] Yang | Robust gaussian graphical modeling with the trimmed graphical lasso[END_REF] If Assumption 3.5 holds then Assumption 3.4 holds with

τ 1 = f (X B ) log p g and τ 2 = f (X B ) b log p g
We can now establish the consistency of arbitrary local optima found with Algorithm 1 under the Frobenius norm.

Theorem 3.8. Let ( w, h) be any local optimum of the ob-

jective in Equation (3). Let λ = C log p n where C = 4 max 2 √ αγ + γ 2 √ 2 log p , f (X B ) , γ = 2 √ 2 L * (Lw + J) -1 ∞ , and D = p(p -1)/2. Then, under assump- tions 3.2 -3.6, if n ≥ f (X B ) + 3RC κ 2 log p it holds that L w-Lw F ≤ 1 √ 2κ 3 2 C log p n D+ b log p g f (X B )
with probability at least 1 -1/p α-2 for some α > 2.

The proof can be found in Appendix B.

Experiments

We empirically demonstrate that our method can effectively cope with the presence of outliers by conducting experiments on simulated data. Results on real world data can be found in Appendix D.

We considered two different simulated settings where the training sets were contaminated with items sampled from a distribution with a different graph structure. That is, we generated data from a mixture of Gaussians as

X ∼ p 0 N (0, (Θ ) -1 ) + (1 -p 0 )N (0, (Θ Corr ) -1 )
where Θ is the true parameter and p 0 is the probability of sampling a "good" element.

The scenarios we considered are:

• Sc1: The graph structures for Θ and Θ Corr are Barabasi-Albert graphs of degre 1 generated with different seeds. For this scenario we set p = 50, n = 250, and p 0 = 0.8. The weights of the edges were uniformely sampled from U[0.25, 2].

• Sc2: Θ is a graph with three hubs sampled from a stochastic block model with intra-cluster probability 0.05 and inter-cluster probability 0.003. Θ Corr is a star graph. We fixed p = 21, n = 105, and p 0 = 0.8. The weights of the edges were uniformely sampled from U[1, 2.5].

We compare our methods R-SCAD and R-MCP to their non-robust counterparts, SCAD and MCP. For R-SCAD and R-MCP we tunned the ĝ hyperparameter from a set of values such that ĝ represents the 80%, 85% and 90% of the data. We fixed the update interval hyperparameter as T = 50. For comparison we use the Relative Error (RE) between the estimated parameter Θ and the ground truth Θ , and the F-score defined as

RE = Θ -Θ F Θ F
and FS = tp tp + 0.5(fp + fn) respectively.

Figure 1 summarizes the mean and standard deviation of RE and FS in both scenarios for different values of λ over 100 simulations. We observe that R-SCAD and R-MCP are systematically better that their non-robust counterparts in terms of relative error and F-score.

A. Useful lemma

Lemma A.1. [START_REF] Ying | Nonconvex sparse graph learning under laplacian constrained graphical model[END_REF]) Let λ = 4α log p/nγ, where γ is defined as in Theorem 3.8. Let n ≥ 8α log p for some α > 2. Then

P [ L * (Lw + J) -1 -L * S ∞ ≤ λ/2] ≥ 1 - 1 p α-2

B. Proof of Theorem 1

Proof. Let ( w, h) be any local optimum of Equation (3). Lemma 3.3 tells us that L w -Lw F ≤ 1. We can then apply Lemma 3.1 with r = 1. Let κ = ( Lw * 2 + 1) -2 . Then,

L w -Lw 2 F ≤ 1 κ -L * (L w + J) -1 + L * (Lw + J) -1 , w -w (4)
On the other hand we have

0 ≤ L * (L w + J) -1 -L * S h , w -w -z, w -w (5)
which follows from the fact that ( w, h) is a local optimum. Combining equations ( 4) and ( 5) we obtain

L w -Lw 2 F ≤ 1 κ L * (Lw + J) -1 -L * S h , w -w -z, w -w = 1 κ L * (Lw + J) -1 -L * S h , w -w (i) + L * S h -L * S h , w -w (ii) -z, w -w (iii)
We now proceed to bound all three terms above.

To bound (i) we apply Hölder's inequality and Assumption 3.2:

L * (Lw + J) -1 -L * S h , w -w ≤ L * (Lw + J) -1 -L * S h ∞ w -w 1 ≤ 1 4 λ w -w 1
To bound (ii), note that by Assumption 3.4 we have

L * S h -L * S h , w -w ≤ τ 1 w -w 1 + τ 2 w -w 2 ≤ 1 4 λ w -w 1 + τ 2 w -w 2
Finally, we use Cauchy-Schwarz to bound (iii) as follows:

z, w -w ≥ -z 2 w -w 2 -z, w -w ≤ z 2 w -w 2 ≤ p(p -1)/2 z ∞ w -w 2 ≤ p(p -1)/2λ w -w 2 Let D = p(p -1)/2.
Putting all together we obtain

L w -Lw 2 F ≤ 1 κ λ 2 w -w 1 + (τ 2 + Dλ) w -w 2 ≤ 1 κ 3 2 Dλ + τ 2 w -w 2 ≤ 1 κ 3 2 Dλ + τ 2 1 √ 2 L w -Lw F
where the last inequality follows from the fact that √ 2 w -w 2 ≤ L w -Lw F . It follows that

L w -Lw F ≤ 1 κ √ 2 3 2 Dλ + τ 2 C. Proof of Lemma 1
Proof. Assume L w -Lw F ≥ 1. Let w t = w + t( w -w ). Take t = 1/ L w -Lw F . Then t(L w -Lw ) F = 1 and hence we can apply Lemma 3.1:

-L * (Lw t + J) -1 + L * (Lw + J) -1 , w t -w ≥ ( Lw 2 + 1) -2 Lw t -Lw 2 F , then t -L * (L w + J) -1 + L * (Lw + J) -1 , w t -w ≥ ( Lw 2 + 1) -2 Lw t -Lw F = ( Lw 2 + 1) -2 t L w -Lw 2 F , so we have -L * (L w + J) -1 + L * (Lw + J) -1 , w t -w ≥ ( Lw 2 + 1) -2 L w -Lw 2 F . (6)
Since ( w, h) is a local optimum we have,

0 ≤ L * (L w + J) -1 -L * S h , w -w -z, w -w . (7)
Combining ( 6) and ( 7) we obtain

( Lw 2 + 1) -2 t L w -Lw 2 F ≤ L * (Lw + J) -1 -L * S h -z, w -w .
After adding and subtracting L * S h and rearranging, the right hand side becomes

L * (Lw + J) -1 -L * S h , w -w (i) + L * S h -L * S h , w -w (ii) -z, w -w (iii)
To bound (i) we apply Hölder's inequality and Assumption 3.2:

L * (Lw + J) -1 -L * S h , w -w ≤ L * (Lw + J) -1 -L * S h ∞ w -w 1 ≤ 1 4 λ w -w 1
To bound (ii), we use Assumption 3.4:

L * S h -L * S h , w -w ≤ τ 1 w -w 1 + τ 2 w -w 2 ≤ 1 4 λ w -w 1 + τ 2 w -w 2
Finally, we use Cauchy-Schwarz to bound (iii) as follows:

z, w -w ≥ -z 2 w -w 2 -z, w -w ≤ z 2 w -w 2 ≤ p(p -1)/2 z ∞ w -w 2 ≤ p(p -1)/2λ w -w 2
Putting all terms together and using our assumption on λ yields the result.

D. Results on Real-world Data

We also carried experiments on the 2019-nCoV dataset available in the R package nCov2019. The dataset collects information about 98 Chinese patients who were affected by the Covid infection. Patients are classified into those who survived the disease and those who did not.

We use the same experimental setup as in [START_REF] Ying | Nonconvex sparse graph learning under laplacian constrained graphical model[END_REF], that is, we transformed categorical features into one-hot encodings, which results in a design matrix X ∈ R 32×98 . Our goal is to learn a graph from the data where the nodes represent patients and edges encode a notion of similarity between patients. A desirable property for a graph is for it to connect nodes that belong to the same class membership, "survived" or "did not survive". Note that this behavior encodes a notion of smoothness of the graph with respect to the survival status of the patients.

To compare R-SCAD and R-MCP to their non-robust counterparts, SCAD and MCP, we estimated graphs for each method for different values of λ and collected the percentage of "bad" edges with respect to the total number of recovered edges. A "bad" edge in this context is an edge that connects two nodes belonging to different class memberships. For both robust methods we fixed the hyperparameter ĝ = 29, and T = 50. Table 3 shows that graphs estimated with the robust methods are in general better at capturing the different node clusters, in contrast to non-robust methods. R-MCP generally outperforms MCP by a considerable margin and it is on a par with MCP for λ = 0.01. Although R-SCAD lags behind SCAD for two values of λ, it outperforms SCAD by a large margin for smaller values of λ, and it is on a par with SCAD for λ = 0.1. The fact that robust methods better capture node clusters for many different values of λ demonstrates that this "grouping effect" does not depend on the sparsity of the estimated graph.

  (k) i . Then the surrogate function can be obtained by lin-Algorithm 1 Trimmed LCGGM Learning

Figure 1 .

 1 Figure 1. Mean relative error and F-score for scenarios Sc1 and Sc2, for different values of λ.

Table 1 .

 1

					Table 2.
	λ	SCAD R-SCAD	λ	MCP R-MCP
	0.0001	18.73	16.41	0.0001 18.74	17.93
	0.001	18.62	11.58	0.001	18.70	11.69
	0.01	18.03	21.35	0.01	18.15	18.67
	0.1	16.75	16.84	0.1	17.44	16.49
	0.5	10.18	14.40	0.5	11.90	7.67

Table 3 .

 3 Percentage of "bad" edges in SCAD and R-SCAD (left), and MCP and R-MCP (right) for different values of λ.


Frederique Battin
email: leclerc@ensic.inpl-nancy.fr

Keywords: Premixed laminar flame, methane, n-butylbenzene, modelling, diesel fuel

To better understand the chemistry involved during the combustion of components of diesel fuel, the structure of a laminar lean premixed methane flame doped with n-butylbenzene has been investigated. The inlet gases contained 7.1% (molar) of methane, 36.8% of oxygen and 0.96% of n-butylbenzene corresponding to an equivalence ratio of 0.74 and a ratio C 10 H 14 / CH 4 of 13.5%. The flame has been stabilized on a burner at a pressure of 6.7 kPa using argon as diluent, with a gas velocity at the burner of 49.2 cm/s at 333 K. Quantified species included the usual methane C 0 -C 2 combustion products, but also 16 C 3 -C 5 hydrocarbons, 7 C 1 -C 3 oxygenated compounds, as well as 20 aromatic products, namely benzene, toluene, phenylacetylene, styrene, ethylbenzene, xylenes, allylbenzene, propylbenzene, cumene, methylstyrenes, butenylbenzenes, indene, indane, naphthalene, phenol, benzaldehyde, anisole, benzylalcohol, benzofuran, and isomers of C 10 H 10 (1-methylindene, dihydronaphtalene, butadienylbenzene). A new mechanism for the oxidation of n-butylbenzene is proposed whose predictions are in satisfactory agreement with measured species profiles in flames and flow reactor experiments. The main reaction pathways of consumption of n-butylbenzene have been derived from flow rate analyses.

INTRODUCTION

If many detailed kinetic models are available for the oxidation of mixtures representative of gasolines, they are much less numerous in the case of diesel fuels because of their more complex composition. The constituents of diesel fuel contain from 10 to 20 carbon atoms and include about 30% (by mass) of alkanes, the remaining part being mainly alkylcyclohexanes (24%), alkyldecalines (15%), alkylbenzenes (10%) and polycyclic naphtenoaromatic compounds [START_REF] Battin | Guibet Fuels and engines[END_REF]. If the oxidation of alkanes has been extensively studied, the abundance of models diminishes considerably when other families of components are considered and very few models exist for substituted cycloalkanes and aromatics compounds, except for toluene [START_REF] Battin-Leclerc | [END_REF]. The oxidation of alkylbenzenes with alkyl side-chains from C 2 to C 4 has been studied in a flow reactor in Princeton [3][4][5][6] and that of n-propylbenzene in a jet-stirred reactor at Orleans [7]. The autoignition of n-butylbenzene has also been investigated in a rapid compression machine in Lille [8][9][START_REF] Crochet | Proceedings of the 3 rd European Combustion Meeting[END_REF]. While premixed flames of ethylbenzene [START_REF] Ergut | [END_REF] and non-premixed methane flames doped with ethylbenzene and isomers of propylbenzene and butylbenzene [12] have been studied, no measurement in a premixed flame containing butylbenzene has yet been reported.

The first purpose of the present paper is to experimentally investigate the structure of a premixed laminar methane flame doped with n-butylbenzene. The use of a methane flame will allow us to be more representative of the combustion mixtures containing large hydrocarbons, such as those present in a diesel fuel, than hydrogen or unsaturated C 2 flames. Large alkanes decompose to methyl radicals and the use of a methane flame will capture the involved chemistry. This study will be performed using a lean flame for the chemistry to be better representative of that occurring in engines controlled via Homogenous Charge Combustion Ignition (HCCI), which are under development. The second objective is to use these results in order to develop a new mechanism for the oxidation of n-butylbenzene based on our experience in modeling the reactions of both alkanes [13] and light aromatic compounds (benzene [14] and toluene [15]).

EXPERIMENTAL PROCEDURE

The experiments were performed using an apparatus developed in our laboratory to study temperature and stable species profiles in a laminar premixed flat flame at low pressure and recently used in the case of rich methane flames doped by light unsaturated soot precursors [16][17][18]. A scheme of the apparatus is presented in figure 1.

FIGURE 1

The body of the flat flame matrix burner, provided by McKenna Products, was made of stainless steel, with an outer diameter of 120 mm and a height of 60 mm (without gas/water connectors). This burner was built with a bronze disk (95% copper, 5% tin). The porous plate (60 mm diameter) used to assist the flame stabilization was water cooled (water temperature: 333 K) with a cooling coil sintered into the plate. The burner could be operated with an annular co-flow of argon to favor the stabilization of the flame. This horizontal burner was housed in a water-cooled vacuum chamber evacuated by two primary pumps and maintained at 6.7 kPa by a regulation valve. This chamber was equipped of four quartz windows for optical access, a pressure transducer (MKS 0-100 Torr), a microprobe for samples taking and a thermocouple for temperature measurements. The burner could be vertically translated, while the housing and its equipments were kept fixed. A sighting telescope measured the position of the burner relative to the probe or the thermocouple with an accuracy of 0.01 mm. The flame was lit using an electrical discharge.

Gas flow rates were regulated by RDM 280 Alphagaz and Bronkhorst (El-Flow) mass flow regulators. Methane (99.95 % pure) was supplied by Alphagaz -Air Liquide. Oxygen (99.5% pure) and argon (99.995% pure) were supplied by Messer. Liquid butylbenzene, supplied by Alfa Aesar (purity 99%), was contained in a glass vessel pressurized with argon. After each load of the vessel, argon bubbling and vacuum pumping were performed in order to remove oxygen traces dissolved in the liquid hydrocarbon fuel. The liquid reactant flow rate was controlled by using a liquid mass flow controller, mixed to the carrier gas and then evaporated by passing through a single pass heat exchanger, the temperature of which was set above the boiling point of the mixture. Carrier gas flow rate was controlled by a gas mass flow controller located before the mixing chamber.

Temperature profiles were obtained using a PtRh (6%)-PtRh (30%) type B thermocouple (diameter 200 µm). The thermocouple wire was supported by an arm and crossed the flame horizontally to avoid conduction heat losses. The junction was located at the centre of the burner.

The thermocouple was coated with an inert layer of BeO-Y 2 O 3 to prevent catalytic effects [19].

The ceramic layer was obtained by dipping the thermocouple in a hot solution of Y 2 (CO 3 ) 3 (93% mass.) and BeO (7% mass.) followed by drying in a Meker burner flame. This process was repeated about ten times until the whole metal was covered. Radiative heat losses are corrected using the electric compensation method [20].

The sampling probe was constructed of quartz with a hole of about 50 µm diameter (d i ). The probe was finished by a small cone with an angle to the vertical of about 20°. For temperature measurements in the flames perturbed by the probe, the distance between the junction of the thermocouple and the end of the probe was taken equal to two times d i , i.e. to about 100 µm. The sampling quartz probe was directly connected to a heated transfer line made with a passivated stainless steel tube and heated at 423 K. This line was itself connected to a heated pressure transducer (MKS 0-100 Torr) and through heated valves to a turbomolecular pump, a pyrex line and a heated 10 ml stainless steel loop located inside a gas chromatograph. The pressure in the transfer line was always below 1.3 kPa (10 Torr) so that the pressure drop between the flame and the inlet of the probe ensured rapid quenching of reactions. The on-line connection to a chromatograph was made via a heated transfer line in order to analyse compounds above C 6 .

The analyses could be then performed according to two methods: Gas samples of compounds with a sufficient vapour pressure were directed through the pyrex line towards a volume which was previously evacuated by the turbo molecular pump down to 10 -7 kPa and which was then filled up to a pressure of 1.3 kPa and collected in a pyrex loop after compression by a factor 5 by using a column with rising mercury. Pressures in the pyrex line were measured by a MKS 0-10 Torr pressure transducer before compression and by a MKS 0-1000 Torr pressure transducer after compression. A chromatograph with a Carbosphere packed column and helium or argon as carrier gas was used to analyse O 2 , H 2 , CO and CO 2 by thermal conductivity detection and CH 4 , C 2 H 2 , C 2 H 4 , C 2 H 6 by flame ionisation detection (FID). Water was detected by TCD but not quantitatively analyzed. Calibrations were performed by analysing a range of samples containing known pressures of each pure compound to quantify. Mole fractions were derived from the known total pressure of the gas included in the pyrex loop. This method of sampling could allow us to analyze heavier hydrocarbons from allene to toluene using a chromatograph with a Haysep packed column with FID and nitrogen as gas carrier gas as in our previous study [18], but it was found more accurate and faster to use the second method described hereafter. For the compounds for which both methods could be used, a very good agreement was obtained between both measurements. Gas samples were directed towards the loop located inside the on-line gas chromatograph and previously evacuated by the turbo molecular pump down to 10 -7 kPa. The loop was then filled up to a pressure of 1.3 kPa and contents were immediately analysed using FID and helium as gas carrier gas by using a capillary column (a HP-Plot Q or a HP-1 column). This method of sampling was used to analyze methane, all the hydrocarbons from C 3 and the oxygenated products other than carbon oxides.

Figure 2 presents a typical chromatogram obtained with the HP-Plot Q column. The identification of these compounds was performed using a gas chromatograph with mass spectrometry detection (GC-MS) fitted with the same column or by comparison of retention times when injecting the product alone in gas phase. Four C 3 species (propyne, allene, propene and propane), nine C 4 species (diacetylene, vinylacetylene, 1-butyne, 2-butyne, 1,3-butadiene, 1,3-butadiene, 1-butene, iso-butene and butane) and four C 5 species (cyclopentene, cyclopentadiene, 1,3-pentadiene and isoprene (2-methyl-1,3-butadiene) were observed. The butenes, mainly 1-butene and iso-butene, were not well separated. In this lean flame, it was also possible to quantify seven oxygenated species (methanol, ketene, acetaldehyde, ethanol, acrolein, propanal and acetone). For methane, C 3 species, 1,3-butadiene and butenes, calibrations were performed by analyzing each pure compound for a range of known pressures in the transfer line and mole fractions were derived from the known total pressure. For other species, we have used the fact that the response of an FID detector is linked to the Effective Carbon Number (ECN) in the molecule [21]. The detector response coefficient (DRC GP-i ) of each compound i can be calculated from the response of a molecule of close structure (DRC GP-mod ), the ECN of this molecule (ECN mod ) and that of the compound i (ECN i ):

DRC GP-i = DRC GP-mod × (ECN i / ECN mod )
ECN is obtained as the sum of the contributions of all the atoms of carbon and oxygen included in the molecules. These contributions have been taken to 1 for aliphatic and aromatic atoms of carbon, 0.95 for olefinic ones, 1.3 for acetylenic ones and 0 for carbonyl ones (C=O) and to -1, -0.6 and -0.25 for atoms of oxygen included in ethers, in primary alcohols and in tertiary alcohols, respectively [21].

FIGURE 2

Figure 3 presents a typical chromatogram obtained with the HP-1 column. The identification of these compounds was performed using a gas chromatograph with GC-MS detection fitted with the same column. Two C 6 cyclic species (methycyclopentene and methylcyclopentadiene), twelve monocyclic aromatic compounds (benzene, toluene, xylenes, phenylacetylene, ethylbenzene, styrene, allylbenzene, n-propylbenzene, cumene (iso-propyl-benzene), methylstyrenes, butenylbenzenes (compounds with C 10 H 12 as chemical formula) and butadienylbenzene) and 5 bicyclic aromatic compounds (indene, indane, methylindene, dihydronaphthalene and naphthalene) were identified as the most probable products according to their mass spectrum. Several peaks have been observed for xylenes, methylstyrenes and compounds with C 10 H 10 as chemical formula, which have been identified as butadienylbenzenes, methylindene and dihydronaphthalene. In this last case, the separate quantification of the different isomers was not possible. Five oxygenated aromatic compounds (phenol, benzaldehyde, benzylalcohol, anisole (methylphenylether) and benzofuran) were also quantified. For those compounds having a low vapour pressure, calibrations were performed in two steps. First, a calibration was made for gas-phase cyclohexane, a cyclic compound with a much larger vapour pressure, by analysing this species for a range of known pressures introduced in the transfer line.

This calibration allowed the detector response coefficient of gas-phase cyclohexane (DRC GP-CHX ) to be obtained: DRC GP-CHX = A GP-CHX / P GP-CHX with:

A GP-CHX : the area of the cyclohexane peak in the case of the gas-phase injection, P GP-CHX : the cyclohexane partial pressure in the transfer line.

The detector response coefficient (DRC GP-P ) of each gas-phase aromatic compound of interest, P, was derived from its relative detector response coefficient (DRC L-P ) compared to that of cyclohexane (DRC L-CHX ) obtained by injecting an acetone solution containing known mass of both cyclohexane (m CHX ) and the P (m P ) in liquid phase: DRC GP-P = (DRC L-P / DRC L-CHX ) × DRC GP-CHX × (M P / M CHX ) with: DRC L-CHX = A L-CHX / m CHX , A L-CHX is the area of the cyclohexane peak in the case of the liquid injection, DRC L-P = A L-P / m P , A L-P is the area of the P peak in the case of the liquid injection, M CHX the molar mass of cyclohexane, Mp the molar mass of P. This method of calibration has been used for all major aromatic compounds, except from methylstyrenes, butenylbenzenes, compounds with C 10 H 10 as chemical formula, phenol, benzylalcohol and benzofuran for which the ECN method has been used.

FIGURE 3

Calculated uncertainties on the species quantifications were about ± 5% for major hydrocarbon compounds and ± 10% for oxygen, hydrogen, carbon oxides and the minor hydrocarbon products. The detection limit of the FID is about 2 ppm.

EXPERIMENTAL RESULTS

A laminar premixed flat flame has been stabilized on the burner at 6.7 kPa (50 Torr) with a gas flow rate of 5.44 l/min corresponding to a gas velocity at the burner of 49.2 cm/s at 333 K and with mixtures containing 7.1% (molar) of methane, 36.8% of oxygen and 0.96% of n-butylbenzene corresponding to an equivalence ratio of 0.74. showing that the presence of the probe induces a thermal perturbation involving a lower temperature. Without the probe, the lowest temperatures measured the closest to the burner (0.47 mm above) were around 1095 K. Due to the thinness of this lean flame and the size of the thermocouple, it was not possible to measure lower temperatures. The highest temperatures were reached from 5 mm above the burner and were around 1960 K.

FIGURE 4

Figure 5 presents the profiles of both hydrocarbon reactants and shows that n-butylbenzene (fig 5b) is consumed close to the burner, at 1 mm height, while some methane (Fig. 5a) remains up to 2 mm. 6b) and of the main C 0 -C 2 species involved in the combustion of methane versus the height above the burner. The mole fraction of water has been obtained from a material balance on the other major species. In this lean flame, the profiles of carbon monoxide (fig. 6c) and hydrogen display a marked maximum at 2 mm height and the major final products are for a large extent carbon dioxide (fig. 6d) and water. Ethylene (fig. 7c) is the most abundant C 2 species and is produced first. It reaches its maximum concentration close to the burner around 1 mm. The profiles of ethane (fig. 7d) and of acetylene (fig. 7b) peak around 1.1 mm and 1.2 mm, respectively.

FIGURES 6 AND 7

Figure 8 presents the profiles of the observed C 3 products, with propene (fig. 8c) and propane (fig. 8d) peaking first around 1 mm above the burner, while the maxima for allene (fig. 8b) and propyne (fig. 8a) are around 1.25 mm. The peak mole fractions observed for propene and propane are more than ten times larger than that of allene and propyne. 10c) and isoprene (fig. 10d) reach their maxima first around 1 mm above the burner. The profiles of methylcyclopentadiene (fig. 10e) peak further around 1.1 and the maxima of that of cyclopentadiene (fig. 10a), the most abundant species of this series with a peak concentration of 50 ppm, and methylclopentene (fig. 10f) are the last ones around 1.5 mm.

FIGURE 10

Figure 11 displays the profiles of light oxygenated species. Methanol (fig. 11a), ketene (fig. 11b), acetaldehyde (fig. 11c) and ethanol (fig. 11d), which are intermediate products of the combustion of methane, as well as propanal (fig. 11f), are produced very early and reach their maximum concentration close to the burner, around 0.9 mm. The profiles of acroleïn (fig. 11e) and acetone (fig. 11g) peak around 1.2 mm. The most abundant of these species are methanol and acetaldehyde with peak concentrations above 150 ppm. 13a), n-propylbenzene (fig. 13b) and butenylbenzenes (fig. 13e) occur first around 1 mm above the burner, followed by that of benzene (fig. 12a), cumene (fig. 13c) and methylstyrenes (fig. 13d) around 1.3 mm and finally that of phenylacetylene (fig. 12d) around 1.5 mm. The most abundant of these species are benzene with a peak concentration of 800 ppm, toluene with a peak concentration of 800 ppm and styrene with a peak concentration of 1000 ppm)and to a lesser extent allylbenzene and butenylbenzenes with peaks concentration about 100 ppm. Xylenes, phenylacetylene and n-propylbenzene are present in low amounts with peak concentrations below 15 ppm.

FIGURES 12 AND 13

Figure 14 displays the profiles of bicyclic aromatic compounds. The maxima of the profiles of indene (fig. 14a) and the isomers of C 10 H 10 (fig. 14d) are observed first around 1 mm above the burner, followed by that of indane (fig. 14b) around 1.3 mm and finally that of naphthalene (fig. 14c) around 1.8mm. The most abundant of these species is indene with a peak concentration of 50 ppm. Indane and naphthalene are present in very low amounts with peak concentrations well below 10 ppm.

FIGURE 14

Finally figure 15 presents the profiles of oxygenated aromatic species. Benzaldehyde (fig. 15b), the most abundant of these species with a peak concentration of 150 ppm, and benzylalcohol (fig. 15c) reach their maxima first around 1 mm above the burner. The profiles of phenol (fig. 15a) and anisole (fig. 15d) peak further around 1.4 mm and the maximum of that of benzofuran (fig. 15e) is the last one around 1.5 mm. This last compound, which can be a precursor of dioxins, corresponds to a very low peak concentration of 2 ppm.

FIGURE 15

DESCRIPTION OF THE PROPOSED MECHANISM

The mechanism proposed here to model the oxidation of n-butylbenzene includes the previous mechanisms that were built by our team to model the oxidation of C 3 -C 5 unsaturated hydrocarbons [22,23,[16][17][18] benzene [14] and toluene [15]. Thermochemical data were estimated by the software THERGAS developed in our laboratory [24], which is based on the group additivity methods proposed by Benson [START_REF] Benson | Thermochemical kinetics[END_REF], apart from the heat of formation of biaromatic species which were taken from Burcat and Ruscic [START_REF] Burcat | Third Millennium ideal gas and consensed phase thermochemical database for combustion with updates from active thermochemical tables[END_REF]. The complete mechanism is available as supplementary material.

Reaction base for the oxidation of C 3 -C 5 unsaturated hydrocarbons [22,23,[16][17][18] This C 3 -C 5 reactions database was built from a review of the literature and is an extension of our previous C 0 -C 2 reactions database [START_REF] Barbé | [END_REF]. This C 0 -C 2 reactions database includes all the unimolecular or bimolecular reactions involving radicals or molecules including carbon, hydrogen and oxygen atoms and containing less than three carbon atoms. The kinetic data used in this base were taken from the literature and are mainly those proposed by Baulch et al. [28] and Tsang et al. [29].

The C 3 -C 5 reactions database [22][23]

includes reactions involving C 3 H 2 , C 3 H 3 , C 3 H 4 (allene and propyne), C 3 H 5 , C 3 H 6 , C 4 H 2 , C 4 H 3 , C 4 H 4 , C 4 H 5 , C 4 H 6 (1,3-butadiene, 1,2-butadiene,
1-butyne and 2-butyne), C 4 H 7 (6 isomers), some linear and branched C 5 compounds and well as cyclopentene and derived species, and the formation of benzene and toluene.

In this reactions base, pressure-dependent rate constants follow the formalism proposed by Troe [30] and efficiency coefficients have been included. It has been validated by modeling experimental results obtained in a jet-stirred reactor for methane and ethane, profiles in laminar flames of methane, acetylene and 1,3-butadiene and shock tube autoignition delay times for acetylene, propyne, allene, 1,3-butadiene, 1-butyne and 2-butyne. The agreement in modeling concentration profiles in laminar flames of pure methane is satisfactory at an equivalence ratio of 1.55 [16], but deteriorates for the formation of ethylene and acetylene at an equivalence ratio of 0.4. An improved version has recently been used to model the structure of a laminar premixed flame of methane doped with allene, propyne, 1,3-butadiene and cyclopentene [16][17][18].

Compared to this last version, reactions of formation and consumption of acetone have been added.

Mechanisms for the oxidation of benzene [14] and toluene [15] Our mechanism for the oxidation of benzene contains 135 reactions and includes the reactions of benzene and of cyclohexadienyl, phenyl, phenylperoxy, phenoxy, hydroxyphenoxy, cyclopentadienyl, cyclopentadienoxy and hydroxycyclopentadienyl free radicals, as well as the reactions of ortho-benzoquinone, phenol, cyclopentadiene, cyclopentadienone and vinylketene, which are the primary products yielded. Validations have been made using for comparison experimental results obtained in a jet-stirred and a plug flow reactors, profiles in a laminar flame of benzene and shock tube autoignition delay times.

The mechanism for the oxidation of toluene contains 193 reactions and includes the reactions of toluene and of benzyl, tolyl, peroxybenzyl (methylphenyl), alcoxybenzyl and cresoxy free radicals, as well as the reactions of benzaldehyde, benzyl hydroperoxyde, cresol, benzylalcohol, ethylbenzene, styrene and bibenzyl. Validations have been made using for comparison experimental results obtained in a jet-stirred and a plug flow reactors, and shock tube autoignition delay times.

Compared to these published mechanisms, reactions of formation and consumption of some products have been added. Cumene is obtained by combination of methyl and phenylethyl radicals or of phenyl and iso-propyl radicals. Methylstyrene is produced by addition of methyl radicals to styrene or by combination of phenylvinyl and methyl radicals. Anisole and benzofuran both derive from phenoxy radicals.

Mechanism proposed for the oxidation of n-butylbenzene

Table I presents the reactions of n-butylbenzene and derived species. The complete mechanism is then composed of this submechanism and of the submechanisms described above.

It involves 210 species in 1478 reactions and can be used to run simulations using CHEMKIN softwares [START_REF] Kee | Sandia Laboratories Report[END_REF]. Table II presents the names, the formulae and the heats of formation of the aromatic species included in this mechanism and containing at least 9 atoms of carbon.

TABLES I AND II

The primary mechanism contains the reactions of n-butylbenzene and of the radicals directly deriving from it. We have considered the unimolecular reactions of n-butylbenzene, the ipso-additions and the H-abstractions by oxygen molecules and radicals present with important concentrations. The rate constants of the unimolecular decompositions involving the breaking of a C-C bond (reactions 1-4) have been calculated from the modified collision theory and thermokinetic relationships [START_REF] Benson | Thermochemical kinetics[END_REF]. The rate constants of the unimolecular decompositions to give H atoms and phenylbutyl radicals (reactions 5-8) have been deduced from that of the reverse reaction, k = 1.0 × 10 14 s -1 according to Allara et al. [34]. The kinetic parameters of the bimolecular initiations with oxygen molecules (reactions 9-12) have been obtained using the correlation proposed by Ingham et al. [35]. The rate constants of the ipso-additions of hydrogen and oxygen atoms and of hydroxyl and methyl radicals (reactions 13-17) have been deduced from value proposed for the similar reactions in the cases of benzene and toluene. The rate constants for the abstractions of alkylic H-atoms (reactions 19-21, 23-25, 27-29, 31-33, 35-37) were deduced from the correlations proposed by Buda et al. for alkanes [13] and those for the abstractions of allylic H-atoms (18,22,[START_REF] Burcat | Third Millennium ideal gas and consensed phase thermochemical database for combustion with updates from active thermochemical tables[END_REF]30,34,(38)(39)(40)(41) were deduced from the correlations proposed by Heyberger et al. [39] and by analogy with similar reactions of ethylbenzene [15].

The reactions of phenybutyl and phenylpropyl radicals were mainly derived from the reactions of alkyl and alkenyl radicals generated by EXGAS software [13,33,[41][42][43][44]. In the case of phenylbutyl radicals, reactions involved isomerizations (reactions [42][43][44][45][46], decompositions by breaking of a C-C bond to form styrene, ethylene, propene and 1-butene as stable molecules, the formation of phenylbutenes by breaking a C-H bond (reactions 48, 51-52, 54-55 and 57) or by oxidation with oxygen molecules (reactions 58 to 63). Termination steps (reactions 64 to 67) were only written for the resonance stabilized 4-phenylbut-1-yl radicals: combinations with HO 2 radicals led to phenylbutoxy radicals, combinations with CH 3 radicals to 2-phenyl-n-pentane and disproportionnations with benzyl and allyl radicals to 1-phenyl-1-butene, toluene and propene, respectively. Phenylpropyl radicals were considered to react by isomerisation followed by In the case of the isomers of phenylbutene and allylbenzene, we have considered the bimolecular initiations with oxygen molecules, the ipso-additions of H-atoms and methyl radicals to give benzene and toluene, the additions to the double bond followed by the decompositions of the obtained adducts and the H-abstractions to give phenylalkenyl radicals, which can react by decompositions by β-scission, oxidations, cyclizations and combinations. The rate constants of the additions to the double bond and of the formation and the consumption of alkenyl radicals were derived from our previous work on the oxidation of alkenes [39,41,44].

The cyclizations of phenylbutenyl and phenylpropenyl radicals lead to bicyclic isomers of C 10 H 10 (dihydronaphthalene and methylindene) and indene, respectively. Their rate constant has been assumed equal to that proposed by Gierzak [46] for the cyclization of pentenyl radicals to give cyclopentyl radicals. For propylbenzene, unimolecular and bimolecular initiations, ipso-additions of H-atoms and methyl radicals and H-abstractions have been written. In the case of 2-phenyl-pentane, which has not been experimentally observed, only four globalized H-abstractions were taken into account. Indane has been considered as yielding indene by bimolecular initiations with oxygen molecules and H-abstractions by small radicals followed by β-scission decompositions of the obtained indanyl radicals.

As only very few studies concerning their oxidation have been published, the reactions of the other products are still very uncertain. For butadienylbenzene, globalized reactions have been written starting by bimolecular initiations with oxygen molecules or H-abstractions by O-and H-atoms and OH radicals, both followed by the cyclization of the obtained radicals yielding naphthalene. The reactions of the bicyclic isomers of C 10 H 10 have also been assumed as leading to naphthalene by molecular dehydrogenation or by H-abstractions. In the case of indene, the addition of OH radicals to the double bond included in the five members cycle has been written, as well as bimolecular initiations with oxygen molecules and H-abstractions by small radicals to give resonance stabilized indenyl radicals, the combinations of which has been assumed to be similar to those written by Da Costa et al. [14] for cyclopentadienyl radicals. Finally, as previously assumed in the work of Bounaceur et al. [51], the reactions of naphthalene and naphthyl radicals can be derived from those of benzene and phenyl radicals.

COMPARISON BETWEEN EXPERIMENTAL AND SIMULATED RESULTS

Simulations were performed using PREMIX from CHEMKIN II [START_REF] Kee | Sandia Laboratories Report[END_REF] using the experimental temperature profile as an input. To compensate the perturbations induced by the quartz probe and the thermocouple, the temperature profile used in calculations is an average between the experimental profiles measured with and without the quartz probe, shifted 0.4 mm away from the burner surface, as shown in figure 4. To decouple the effect due to the increase of equivalence ratio (φ) and that induced by the presence of n-butylbenzene, figures 5 to 7 and 11 display also the results of a simulation performed for a flame containing 7.1% methane and 15.64 % oxygen (with no additive) for φ=0.74, i.e. equal to that of the doped flame. As the temperature rise is mainly influenced by φ, we have used the same temperature profile as that used to model the doped flame. The profile of methane is very similar for the doped and undoped flames at the same φ. However the content in atoms of carbon is 2.35 lower in the pure methane flame compared to the doped one (the C/O ratio is equal to 0.096 in the undoped flame and to 0.225 in the seeded one) which is well reflected by the profiles of carbon oxides and ethane. While the profile of ethane is not much affected, the formation of ethylene and acetylene is increased by the addition of n-butylbenzene by a factor much larger than the increase in the number of atoms of carbon showing clearly to what extent these products derive from the decomposition of this additive. The presence of the additive does not influence the formation of methanol, but increase also considerably that of ketene, acetaldehyde and ethanol. the most abundant species, the formation of benzene, toluene and allylbenzene is very well reproduced, that of styrene is underestimated by a factor 1.5, that of benzaldehyde is overestimated by the same factor and that of butenylbenzenes is overestimated by a factor 2.5.

The poor separation of the C 10 H 12 isomers on the GC column that we used is a possible explanation for this discrepancy. While the simulated profiles of bicyclic compounds notably differ from the experimental ones, the maximum mole fractions obtained for minor products are all modeled with a factor better than 4, except from benzofuran, the formation of which is underestimated by a factor 10 showing that that a production pathway is missing in the model.

The production of phenol is overestimated by a still larger factor.

It is worth noting that the largest disagreements are obtained for several oxygenated compounds. For ethanol, ketene and phenol which are strongly overestimated by the model, it is probable that, despite the heated transfer line, some of these oxygenated compounds may be absorbed on the walls and result in lower measured values than those predicted by our model.

While the reactions of ketene are still very uncertain, those of ethanol and phenol are rather well determined with pressure effects taken into account.

In order to extend the temperature ranges of these simulations, the results of Litzinger et al. [6], which were obtained in a flow reactor at 1060 K, at atmospheric pressure, with nitrogen as bath gas, for an initial concentration of n-butylbenzene of 620 ppm and for an equivalence ratio of 0.98, have also been modeled. Figure 16 displays comparisons between the experimental and computed mole fraction of reactants and main products. This figure shows that a globally correct agreement can be observed. The consumptions of hydrocarbon and oxygen (fig. 16a) are correctly reproduced, as well as the formation of the major products, methane (fig. 16b), ethylene (fig. 16b), benzene (fig. 16c), toluene (fig. 16c), styrene (fig. 16c), allylbenzene (fig. 16c). The production of carbon monoxide (fig. 16b) is overestimated by a factor 1.5 at the longest residence times.

FIGURE 16

DISCUSSION

Figure 17 displays the main flows of consumption of the n-butylbenzene at a temperature about 1140 K corresponding to 87% reactant conversion. A large enough conversion has been chosen, so that the major ways of consumption of the primary products can be observed. Under these conditions, butylbenzene is mainly consumed by H-abstractions (40 % if its consumption) to give the resonance stabilized 4-phenylbut-4-yl radical which is the main source of styrene through a β-scission decomposition. At this temperature, this radical is consumed faster than it is produced which is only possible because of its accumulation at lower temperature.

FIGURE 17

The other important channels of consumption of n-butylbenzene are the H-abstractions leading to the other 4-phenylbutyl radicals (globally 41 % of its consumption). A part of 4-phenylbut-1yl and 4-phenylbut-2-yl radicals isomerizes to give resonance stabilized 4-phenylbut-4-yl radicals. However 4-phenylbut-1-yl radicals mainly decompose to give ethylene and 2phenyleth-1-yl radicals yielding phenyl radicals and ethylene or, for a smaller part, styrene and H-atoms. The decomposition of 4-phenylbut-1-yl radicals explains most part of the formation of ethylene. 4-phenylbut-2-yl radicals are consumed to produced allylbenzene and methyl radicals, butenylbenzene and H-atoms and 1-butene and phenyl radicals. This last channel is the main way to form 1-butene. Phenyl radicals are an important source of toluene and phenoxy radicals yielding phenol, and to a much lesser extent of benzene. Anisol and benzofuran are produced from minor reactions of phenoxy with methyl radicals and acetylene, respectively. As benzofuran is strongly underestimated by the model, it is probable that important way of formation of this compound are missing. Finally 4-phenylbut-3-yl radicals are almost completely consumed to form propene and benzyl radicals, which are the major source of toluene, benzaldehyde, ethylbenzene, benzylalcohol and n-propylbenzene. The formation of xylenes derive from toluene through the ipso-addition of methyl radicals. The largest part of the formation of propene is due to reactions of 4-phenylbut-3-yl radicals.

There are two additional minor channels of consumption of n-butylbenzene by ipso-additions.

The ipso-addition of H-atoms (6 % of its consumption) is the main source of benzene. Butyl radicals are also produced through this way and mostly decompose to form ethylene and ethyl radicals and to a lesser extent 1-butene and H-atoms. A very small fraction of them also combines with H-atoms to give butane. Ethyl radicals are a source of acetaldehyde by reaction with O-atoms and of ethanol by combinations with OH radicals. The ipso-addition of O-atoms (3.7 % of its consumption) leads ultimately to carbon monoxide, as well as to cyclopentadienyl radicals and 1-butene or to H-atom, propene and benzene.

Styrene reacts mainly by additions of OH radicals or of H-atoms and by H-abstractions. The additions of OH radicals are a source of formaldehyde and benzyl radicals and of benzaldehyde and methyl radicals. The additions of H-atoms produce resonance stabilized 2-phenyleth-2-yl radicals, the combinations of which yield cumene and ethylbenzene. Phenylvinyl radicals which are obtained by H-abstractions from styrene react mainly with oxygen molecules yielding resonance stabilized benzoyl radicals which decompose to produce phenyl radicals and carbon monoxide. Minor channels consuming phenylvinyl radicals involve the formation of methylstyrene and of phenylacetylene.

Butenylbenzenes and allylbenzene react mainly by H-abstractions to give resonance stabilized radicals the cyclization of which leads to methylindene and indene, respectively. The cyclization of the alkylic alkenyl radicals, which are formed by another minor channel consuming butenylbenzenes, produces dihydronaphthalene. The cyclization of the phenylpropyl radicals, which are obtained by addition of H-atoms to allylbenzene, explains the formation of indane (not shown in fig. 17, as it corresponds to a very small flow rate). The flow rate of the addition to allylbenzene is much smaller than that of the H-abstractions explaining the lower production of indane compared to indene. The reactions of methylindene and dihydronaphthalene lead to naphthalene or to resonance stabilized cyclopentadienyl radicals. These cyclic C 5 radicals mostly yield cyclopenpentadiene by combination with H-atoms, but also, to a lower extent, methylcyclopentadiene by combinations with methyl radicals, and acetylene or 1,3-butadiene by opening of the ring. A small part of cyclopentadienyl radicals derives also from the CO elimination from phenoxy radicals. Cyclopentene and methylcyclopentene derive from cyclopentadiene and methylcyclopentadiene, respectively, by additions of H-atoms followed by combinations with H-atoms Let us now describe the major ways of formation of the minor species which do not derives directly from butylbenzene. Ethane, propane and methanol are formed by combination of methyl radicals with themselves, ethyl and OH radicals, which explains why their formation is not much affected by the presence of butylbenzene. Allene derives from propene via the formation of resonance stabilized allyl radicals. Propyne is formed by addition of O-atoms to vinylacetylene, which is obtained from cyclopentadienone which derives from benzoquinone, a minor product of the reaction of phenyl radicals with oxygen molecules. The additions of methyl radicals to propyne yield iso-butyl radicals, the recombination of which is the main source of iso-butene.

Diacetylene and butynes derive from vinylacetylene. 1,2-butadiene is obtained from resonance stabilized 1-buten-3-yl radicals, which are produced by H-abstractions from 1-butene.

1,3-pentadiene derives from cyclopentene after H-abstractions and opening of the ring and isoprene from the addition of vinyl radicals to allene. Ketene and acetone both derive from CH 3 CO radicals obtained by H-abstraction from acetaldehyde. Propanal is produced by addition of OH radicals to 1-butene or butenylbenzenes. Acroleïn is mainly obtained by reaction of allyl radicals with O-atoms.

CONCLUSION

This paper presents new experimental results for a lean premixed laminar flame of methane seeded with n-butylbenzene, as well as a new mechanism developed to reproduce the combustion of this substituted aromatic compound, which can be considered as a model molecule of an important class of components of diesel fuels. Profiles of temperature have been measured and mole fraction profiles have been obtained for 55 identified stable species from C 0 to C 10 , including 20 aromatic products and 12 oxygenated compounds other than the reactants.

Several of these species are considered as toxic pollutants, this is the case of oxygenated compounds (acroleïn). Several of these products are also known as soot precursors, this is the case of all the aromatic compounds.

Satisfactory agreement has been obtained between experimental results and simulations, apart from 4 oxygenated species (ethanol, ketene, phenol, benzofuran). The prediction of the profiles of bicyclic aromatic compounds (indene, indane, naphthalene, benzofuran) could also be improved when a better knowledge of their chemistry will be available.

AKNOWLEDGEMENT

This work has been supported by PSA Peugeot Citroën and TOTAL

TABLE I: REACTIONS OF BUTYLBENZENE AND OF DERIVED SPECIES

The rate constants are given (k=A T n exp(-E a /RT)) in cc, mol, s, kcal units. _____________________________________________________________________________

Reactions

A

Reactions of phenylpropyl radical (Φ-C 3 H 6 -1)

Isomerisations-decomposition Φ-C 3 H 6 -1=>styrene+CH 3 3.4×10 [START_REF] Rao | [END_REF] for toluene b : Rate constant of this recombination calculated by the modified collision theory at 1200 K using software KINGAS [33]. c : Rate constant taken equal to that of the recombination of H atoms with alkyl radicals as proposed by Allara et al [34].

Φ-C 4 H 7 -1+O=>benzyl+HCHO+C 2 H 3 V 6.4×10 2.6 -1.1 i (97)
Φ-C 4 H 7 -1+OH=>Φ-C 3 H 6 -1+HCHO 1.4×10 12 0.0 -1.0 i (99) Φ-C 4 H 7 -1+OH=>Φ-C 2 H 4 -1+CH 3 CHO 1.4×10 12 0.0 -1.0 i (100) Φ-C 4 H 7 -1+CH 3 =>Φ-C 2 H 4 -1+C 3 H 6 Y 1.7×10 11 0.0 7.4 i (101) Φ-C 4 H 7 -2+O=>Φ-C 2 H 4 -1+CH 2 CHO 6.4×10 4 2.6 -1.1 i (102) Φ-C 4 H 7 -2+O=>benzyl+C 2 H 3 CHO+H 6.4×10 4 2.6 -1.1 i (103) Φ-C 4 H 7 -2+OH=>Φ-C 2 H 4 -1+CH 3 CHO 1.4×10 12 0.0 -1.0 i (104) Φ-C 4 H 7 -2+OH=>benzyl+C 2 H 5 CHO 1.4×10 12 0.0 -1.0 i (105) Φ-C 4 H 7 -2+CH 3 =>benzyl+C 4 H 8 Y 1.7×10 11 0.0 7.4 i (106) Φ-C 4 H 7 -3+O=>C 6 H 5 #+C 2 H 3 CHO+CH 3 6.4×10 4 2.6 -1.1 i (107) Φ-C 4 H 7 -3+OH=>Φ-CHO+n-C 3 H 7 1.4×10 12 0.0 -1.0 i (108) Φ-C 4 H 7 -3+OH=>benzyl+C 2 H 5 CHO 1.4×10 12 0.0 -1.0 i (109) Φ-C 4 H 7 -3+CH 3 =>Φ-C 3 H 5 Z+C 2 H 5 1.7×10 11 0.0 7.4 i (110) Metatheses Φ-C 4 H 7 -1+O=Φ-C 4 H 6 X+OH 8.8×10 10 0.7 3.3 i (111) Φ-C 4 H 7 -1+O=Φ-C 4 H 6 X(1-3)+OH 8.8×10 10 0.7 3.3 i (112) Φ-C 4 H 7 -1+H=Φ-C 4 H 6 X+H 2 5.4×10 4 2.5 -1.9 i (113) Φ-C 4 H 7 -1+H=Φ-C 4 H 6 X(1-3)+H 2 5.4×10 4 2.5 -1.9 i (114) Φ-C 4 H 7 -1+OH=Φ-C 4 H 6 X+H 2 O 3.0×10 6 2.0 -1.5 i (115) Φ-C 4 H 7 -1+OH=Φ-C 4 H 6 X(1-3)+H 2 O 3.0×10 6 2.0 -1.5 i (116) Φ-C 4 H 7 -1+HO 2 =Φ-C 4 H 6 X+H 2 O 2 6.4×10 3 2.6 12.4 i (117) Φ-C 4 H 7 -1+HO 2 =Φ-C 4 H 6 X(1-3)+H 2 O 2 6.4×10 3 2.6 12.4 i (118) Φ-C 4 H 7 -1+CH 3 =Φ-C 4 H 6 X+CH 4 1.0×10 11 0.0 7.3 i (119) Φ-C 4 H 7 -1+CH 3 =Φ-C 4 H 6 X(1-3)+CH 4 1.0×10 11 0.0 7.3 i (120) Φ-C 4 H 7 -1+benzyl=Φ-C 4 H 6 X+toluene 1.1×10 12 0.0 13.1 l (121) Φ-C 4 H 7 -1+benzyl=Φ-C 4 H 6 X(1-3)+toluene 1.1×10 12 0.0 13.1 l (122) Φ-C 4 H 7 -2+O=Φ-C 4 H 6 Y(2-4)+OH 8.8×10 10 0.7 3.3 i (123) Φ-C 4 H 7 -2+O=Φ-C 4 H 6 X(1-3)+OH 1.7×10 11 0.7 5.9 i (124) Φ-C 4 H 7 -2+H=Φ-C 4 H 6 Y(2-4)+H 2 5.4×10 4 2.5 -1.9 i (125) Φ-C 4 H 7 -2+H=Φ-C 4 H 6 X(1-3)+H 2 1.7×10 5 2.5 2.5 i (126) 31 Φ-C 4 H 7 -2+OH=Φ-C 4 H 6 Y(2-4)+H 2 O 3.0×10 6 2.0 -1.5 i (127) Φ-C 4 H 7 -2+OH=Φ-C 4 H 6 X(1-3)+H 2 O 3.0×10 6 2.0 -0.3 i (128) Φ-C 4 H 7 -2+HO 2 =Φ-C 4 H 6 Y(2- 4 2.6 -1.1 i (98)
Φ-C 4 H 6 Y(2-4)+H=Φ-C 4 H 7 -3 1.0×10 14 0.0 0.0 c (161) Φ-C 4 H 6 X(1-3)+H=Φ-C 4 H 7 -1 1.0×10 14 0.0 0.0 c (162) Φ-C 4 H 6 X(1-3)+H=Φ-C 4 H 7 -2 1.0×10 14 0.0 0.0 c (163) Φ-C 4 H 6 -1+H=Φ-C 4 H 7 -3 1.0×10 14 0.0 0.0 c (164)

Reactions of allylbenzene (Φ-C 3 H 5 Z) and derived radicals

d : Pre-exponentiel factor is the value proposed by Ingham et al [35]. Activation energy calculated at 1200 K by using software THERGAS [24]. e : Rate constant estimated by analogy with the values proposed by Baulch et al [28] for toluene.

f : Rate constant estimated by analogy with the values proposed by Robaugh et al [36] for toluene.

g : Rate constant estimated by analogy with the values proposed by Seta et al [37] for toluene. h : Rate constant estimated by analogy with the values proposed by Tappe et al [38] for toluene. i : Rate constant estimated by using the correlations proposed by Heyberger et al. [39] in the case of alkenes. j : Rate constant estimated by using the correlations proposed by Buda et al. [13] in the case of alkanes.

k : Rate constant estimated by analogy with the values proposed by Heckmann et al. [40] for toluene. l : Rate constant estimated by analogy with the values proposed for ethylbenzene by Bounaceur et al. [15].

m : Rate constant of isomerisation calculated by using the relations described by Buda et al. [13] with the activation energy calculated as the sum of the activation energy for H-abstraction from the substrate by analogous radicals (13.5 kcal/mol for a primary alkylic H-atom, 11 kcal/mol for a secondary alkylic H-atom, 6.5 kcal/mol for a primary allylic H-atom) and the strain energy of the saturated cyclic transition state (27.7 [14] for benzene or phenyl radicals. For bimolecular initiations and H-abstractions, A was multiplied by 8/6 to take into account the number of abstractable H-atoms.

TABLE II: NAMES, FORMULAE AND HEATS OF FORMATION FOR THE AROMATIC SPECIES CONTAINING AT LEAST 9 ATOMS OF CARBON INVOLVED IN THE MECHANISM OF TABLE I.

The heats of formation have been calculated by software THERGAS [24] at 300 K in kcal.mol -1 , apart from biaromatic species for which they have been obtained from Burcat and Ruscic [START_REF] Burcat | Third Millennium ideal gas and consensed phase thermochemical database for combustion with updates from active thermochemical tables[END_REF].

Species

Figure 4

 4 Figure 4 displays the experimental temperature profiles obtained with and without the probe

FIGURE 5 Figures 6

 56 FIGURE 5

FIGURE 8 Figure 8

 88 FIGURE 8

FIGURE 9 Figure 10

 910 FIGURE 9

Figures 12 and 13

 13 Figures 12 and 13 present the profiles of monocyclic aromatic products. The maxima of the

  decomposition (reaction 68) to lead ultimately to styrene, by cyclisation (reaction 69) to give indane, by β-scission decompositions (reactions 70-71) yielding ethylene and allylbenzene as stable molecules, by oxidation (reaction 72) to produce allylbenzene and by combination with H-atoms to form propylbenzene (reaction 73). The reactions of 1-butyl radicals (reactions 74 to 78) are a reduced version of those generated by EXGAS software for their high temperature oxidation. The reactions of oxygenated aromatic radicals are still very uncertain. The globalized reactions considered for butylphenoxy radicals (reactions 79-80) are derived from those of phenoxy and started as an elimination of carbon monoxide. Propylbenzylalcoxy radicals have been assumed to react by isomerizations with a rapid decomposition of the obtained radicals (reactions 81-82) or by direct β-scission decompositions (reactions 83-84). The secondary mechanism includes the reactions of the primary products which are not considered in the mechanisms of the oxidation of benzene and toluene, namely the three isomers of phenylbutene (reactions 85 to 164), allylbenzene (reactions 165 to 182), propylbenzene (reactions 183 to 200), 2-phenyl-n-pentane (reactions 201 to 204) and indane (reactions 205 to 216), as well as those of the derived products, namely butadienylbenzenes (reactions 217 to 226), the bicyclic isomers of C 10 H 10 (reactions 227 to 231), indene (reactions 232 to 243) and naphthalene (reactions 244 to 255).

Figures 5

 5 Figures 5 to 7 and 11 show that the model reproduces satisfactorily the consumption of

Figures 8

 8 Figures 8 to 11 present the comparison between experimental and simulated data for the other

Figure 1 :

 1 Figure 1: Scheme of the apparatus of low-pressure laminar premixed flame. The thick lines

Figure 2 :

 2 Figure 2: Typical chromatogram of light compounds obtained at a distance of 1.40 mm from the

Figure 3 :

 3 Figure 3: Typical chromatogram of heavy compounds obtained at a distance of 1.31 mm from

Figure 4 :

 4 Figure 4: Temperature profiles: experimental measurements performed without and with the

Figure 5 :

 5 Figure 5: Profiles of the mole fractions of both hydrocarbon reactants including a computed

Figure 6 :

 6 Figure 6: Profiles of the mole fractions of oxygen and the main oxygenated products including a

Figure 7 :

 7 Figure 7: Profiles of the mole fractions of hydrogen and C 2 species including a computed

Figure 8 :

 8 Figure 8: Profiles of the mole fractions of C 3 species including a computed comparison between

Figure 9 :

 9 Figure 9: Profiles of the mole fractions of C 4 species. Points are experiments and lines

Figure 10 :

 10 Figure 10: Profiles of the mole fractions of C 5 -C 6 non-aromatic species. Points are experiments

Figure 11 :

 11 Figure 11: Profiles of the mole fractions of oxygenated C 1 -C 3 species. Points are experiments

Figure 12 :

 12 Figure 12: Profiles of the mole fractions of C 6 -C 8 aromatic species. Points are experiments and

Figure 13 :

 13 Figure 13: Profiles of the mole fractions of C 9 -C 10 monoaromatic species. Points are experiments

Figure 14 :

 14 Figure 14: Profiles of the mole fractions of biaromatic species. Points are experiments and lines

Figure 15 :

 15 Figure 15: Profiles of the mole fractions of oxygenated aromatic species. Points are experiments

Figure 16 :

 16 Figure 16: Computed species mole fractions versus residence time during the oxidation of n-butylbenzene in a flow reactor at 1060 K, P= 1atm and Φ = 0.98, with an initial concentration

Figure 17 :

 17 Figure 17: Flow rate analysis for the consumption of the n-butylbenzene for a distance of

  Figure 1

  Figure 2

Figure 3

 3 Figure 3

  Rate constant estimated from a rate constant proposed by Perrin et al.[42] for 1-pentenyl radicals. o : Rate constant estimated from rate constants proposed by Tsang[43] for propene. p : Rate constant estimated by analogy with the values proposed by Touchard et al.[41] in the case of the high temperature oxidation of 1-pentene. q : Rate constant estimated by using the correlations proposed by Touchard et al.[44] in the case of the low temperature oxidation of 1-pentene. Rate constant estimated by analogy with the values proposed by Hippler et al.[45] for toluene. Rate constant estimated by analogy with the values proposed by Dayma et al.[50] for cyclohexyl radical Reactions and rate constant estimated by analogy with the values proposed by Da Costa et al.[14] for cyclopentadienyl radical. Reactions and rate constants obtained by analogy with the values proposed by Da Costa et al.

	kcal/mol for a 3 members ring, 26 kcal/mol for a 4 members ring, 6.3
	kcal/mol for a 5 members ring)

n : r : s : The rapid decomposition of the 3-phenyl-3propyl radicals obtained by isomerization has been assumed. t : Rate constant estimated by analogy with the values proposed by Gierczak et al.

[46]

for pentenyl radical u : Rate constant estimated by analogy with the values proposed by Baulch et al.

[47]

for phenoxy radical v : Rate constant estimated by using the correlation proposed by Glaude et al.

[48]

in the case of ethers. w : Activation energy estimated by using an Evans-Polanyi relationship proposed by Herbinet et al. [49] in the case of β-scissions of C-C and C-H. x : y : z :

Temperature profile used for simulation (shifted 0.4 mm away from the burner)


domain: spi phys info math

text:
Mathieu Flament
email: mathieu.flament@onera.frwilliam.boitier@onera.fr

William Boitier

Julien Floquet
email: julien.floquet@onera.fr

Continuous System Integration in DEVS Formalism

Keywords: Simulation de systèmes à temps continu, systèmes aérospaciaux, P-DEVS, QSS Continuous system simulation, Aerospace systems, P-DEVS, QSS

Dans cet article, nous décrivons d'implémentation des algorithmes classiques d'intégration numérique (Euler, Heun, Runge-Kutta 4) dans le formalisme DEVS, en utilisant la bibliothèque Artis-star. Ce travail nous a permis de les comparer aux méthodes de type QSS (Quantized-State Systems), dans un environnement commun de simulation basé sur DEVS et sur les mêmes cas tests, dont un exemple réaliste d'un système aéronautique.

Introduction

Through the years, Onera has developed many simulation tools to help the design and the performance assessment of various aerospace systems. These simulations integrate flight dynamics models together with guidance, navigation and control (GNC) models and command logic, in order to reproduce the behaviour of a real flight. Most of the physical and GNC models can be written as continuous-time models described by ordinary differential equations. The mature field of numerical integration methods provide a large choice of solutions for continuous system simulations. The Runge-Kutta 4 (RK4) method (fixed step or adaptative) is a very classical algorithm to simu-late aerospace systems and provides a good balance between stability, accuracy and computational cost. For less demanding needs, Euler, Adams-Bashford (AB), Adams-Bashford-Moulton (ABM) methods can also be used. Different books explain these integration methods as in [START_REF] Jedrzejewski | Introduction aux méthodes numériques[END_REF][START_REF] Isaacson | Analysis of numerical methods, john wiley and sons[END_REF].

To ensure the exactness of the computations, the numerical integration of a system of interconnected continuous-time models should be performed as a unique continuous-time model which aggregates the state of all subsystems. In terms of software architecture, this constraint comes into conflict with modularity : it's much simpler to associate a local integrator to each continuous-time sub-model and then connect them as discrete-time models. But, especially for control-loop stability analysis, where delays are critical, this trivial solution is not acceptable. RK4, ABM, and more generally, all predictor/corrector integration methods require an exchange of updated state values between submodels after each derivative evaluation.

In this context, our interest in DEVS formalism is twofold : first, a DEVS-based architecture could ease the sequencing of classical numerical integration methods distributed over a system of interconnected continuous-time models ; and second, the Quantized-State System methods (QSS), which have been developed jointly with DEVS formalism, make the difficulties inherent to classical methods simply vanish. The first formalization of QSS was done by Kofman and Junco in [START_REF] Kofman | Quantized-state systems : A devs approach for continuous system simulation[END_REF]. The research field of Discrete Event Simulation (DEVS) has advanced significantly since its early years in 1976 when B. Zeigler sets its first formalism down. However, the theory is still unexploited by the majority of Modelling & Simulation (M & S) software. For our concern, DEVS formalism, and Parallel DEVS (PDEVS) formalism even more, have potential in the M & S field. The first one has been proven as a common denominator for multi-formalism simulation [START_REF] Prähofer | System theoretic formalisms for combined discrete-continuous system simulation[END_REF][START_REF] Bernard | Embedding dev&dess in devs[END_REF][START_REF] Vangheluwe | Devs as a common denominator for multi-formalism hybrid systems modelling[END_REF], while the second one was written to guarantee distributed simulation.

The developments of DEVS formalism have brought many new integration methods [START_REF] Nutaro | A second order accurate adams-bashforth type discrete event inte-gration scheme[END_REF][START_REF] Migoni | Linearly implicit discrete event methods for stiff odes[END_REF][START_REF] Nutaro | Constructing multi-point discrete event integration schemes[END_REF], including the well-known Quantized-State System methods (QSS) presented also in [START_REF] Kofman | Quantization-based simulation of differential algebraic equation systems[END_REF]. We have performed a benchmark of these methods compared to classical ones, both in terms of numerical accuracy and computational efficiency.

In section 2, we briefly recall the classical numerical integration methods in comparison to the Quantized-State System methods. In section 3, we describe our implementation of the benchmarked integrators using the Artis-star library [START_REF] Ramat | Artis* : Devs modern c++ library[END_REF], and present the results for two applications : a simple test-case model (a linear second-order system) and a real-life aerospace system (Boeing 747 flight dynamics and its vertical autopilot).

Numerical integration methods

Classical methods

In the family of single-step methods, the solutions are based on approximation of derivatives until an order p using s stages. For instance, the Euler's method performs an approximation of first-order derivatives using one stage, making the assumption that this derivative is constant. Additional stages allow to evaluate higher-order derivatives. It is accepted that a single step method of order p needs a number of stage s ≥ p (until fifth order where the condition is strictly s > p). For instance, the classical form of Runge-Kutta is a fourth-order derivative composed of four stages.

For the multi-step methods, the solutions result from interpolation using multiple points. Those methods have the advantage to be computed in one stage. The family is represented by the AB and ABM methods declined for each order and in some case combined to create predictorcorrector (or PECE : Predict-Evaluate-Correct-Evaluate) algorithm.

Quantized-State System integration methods

The development of DEVS formalism has brought new kinds of integration method based on discrete state rather than discrete time.

A first approach was proposed by B. Zeigler and Lee in [START_REF] Zeigler | Theory of quantized systems : formal basis for DEVS/HLA distributed simulation environment[END_REF] called quantized system, though quickly put aside when Kofman and Junco proposed the quantized-state system (QSS) methods in [START_REF] Kofman | Quantized-state systems : A devs approach for continuous system simulation[END_REF] and advanced QSS method are proposed in [START_REF] Migoni | Linearly implicit discrete event methods for stiff odes[END_REF]. The DEVS description of order 1 QSS is described in numerous papers [START_REF] Zeigler | Theory of modeling and simulation[END_REF][START_REF] Kofman | Quantization-based simulation of differential algebraic equation systems[END_REF], though for the sake of clarity we will recall the behaviour of such method.

QSS =X = (x) Y = (y) S = ((z, ż), q) δ int (S, q, σ) = (z + σ ż, ż), z + σ ż, |∆Q/ ż| δ ext (S, q, X, e) = (z + e ż, x), q, σ λ(S, q, σ) = z + σ ż ta(S, q, σ) = σ
with X the set of input events ; Y is the set of output events ; S is the set of sequential states ; δ int is the internal transition function ; δ ext is the external transition function ; λ is the output function ; ta is the time advance function ; q the quantum state ; e the elapsed time since the last event ; z and ż the value of state and the associated derivative ; σ the elapsed time from the initial time of simulation ; ∆Q is the quantum size and

σ =            q -∆Q -(z + e ż) x if x < 0 q + ∆Q -(z + e ż) x if x > 0 ∞ otherwise
The properties of convergence and stability are already proven [START_REF] Kofman | A second-order approximation for devs simulation of continuous systems[END_REF][START_REF] Kofman | An approach for discrete event simulation of physical systems[END_REF], and specifically the stability property has to be noted (stability region is the entire negative plan) to compare with previously discussed methods. Under certain conditions, the error is known to be bounded [START_REF] Kofman | A second-order approximation for devs simulation of continuous systems[END_REF]. For a continuous system written in the state space form : ẋ = Ax + Bu, the error bound at time t, named E(t), can be derived (for QSS1, 2 and 3) by the following formula :

|E(t)| = |V ||Re(Λ) -1 Λ||V -1 |∆Q (1) With : V -1 AV = Λ
V and Λ being respectively the eigenvectors and matrix of eigenvalues of A.

3 Benchmark of integration methods with DEVS Some benchmarks between QSS and classical algorithms have already been performed [START_REF] Floros | Automated simulation of modelica models with qss methods : The discontinuous case[END_REF], but the CPU utilization of the different methods was not always measured with the same simulation engine. Here, thanks to the Artis-star environment developed by E.Ramat [START_REF] Ramat | Artis* : Devs modern c++ library[END_REF], we are able to implement all the evaluated integration algorithms in a common framework.

Integrator Mapping

Our C++ implementation of integrators using the Artis-Star environment is based on the following components (see figure 1). First, the continuous-time model (ODE) is described in a class that implements an interface AbstractODEModel with two main methods :

derivative() : returns

Ẋ = AX + Bu; (2)
-computeOutput() : returns

Y = CX + Du; (3)
with X the state vector, u the input vector, Y the output vector. Note that for reusability, the implementation of continuous-time models is independent of the simulation engine (DEVSbased or not).

Second, the integration algorithm (classical ones or QSS-based ones) is coded as an atomic DEVS model, by implementing the artis::pdevs::Dynamics interface. It is responsible for the different steps of calculations specific to a given algorithm. One integrator is associated with each continuous-time model of the system.

Third, in the case of a system composed of multiple interconnected continuous-time model, these integrators are coupled through a coordinator, which ensures the necessary communication between the different integrators.

FIGURE 1 -ODE as sub-module of atomic block

Implementation of classical integrators. The Euler's algorithm implementation is quite trivial : its state is composed of the current state of the continuous-time model and a discrete variable that corresponds to the calculation phase WAIT, SEND. The lifetime of these phases are respectively the integration timestep δh and 0.

On the WAIT to SEND transition, each integrator evaluates the model state derivative and updates its own state by Euler's formula, given the outputs of the other models. On the SEND to WAIT transition, each integrator computes its outputs from the updated state and publishes a message containing these outputs.

Through the coupling, this message is received by the other ones. For the classical integrator, the coupling is done by the generic coordinator already implemented in the Artis-star library (artis::pdevs::Coordinator).

The Heun's algorithm is the simplest predictorcorrector algorithm (PEC) (see figure 3), and needs a slightly more complex implementation, as the integration step are divided in two substeps. An intermediate exchange of state values between coupled integrators is required after the prediction sub-step. To do so, the Heun's integrator phases are WAIT, SEND-P, SEND-C, with respective lifetime of h, 0 and 0. On the WAIT to SEND-P transition, each integrator performs the computation of the predicted state t i + h ; on the SEND-P to SEND-C transition, they publish a message containing the outputs corresponding to the predicted state and compute the corrected state (final evaluation) at t i + h ; on the SEND-C to WAIT transition and send a message containing the outputs corresponding to the corrected state.

Prediction :

X p k+1 = X k + h.f (t k , X k , u k ) Correction :    X p k+1 = X k + h.f (t k+1 , X p k+1 , u k+1 ) X k+1 = 1 2 X p k+1 + X p k+1 )
with f (t, X, u) the derivative at time t (f (t, X, u) = A.X + B.u for LTI systems)

FIGURE 2 -Heun's algorithm

These implementation principles for Heun's algorithm can be easily extended to higher orders of PEC integration algorithms, by adding stages of prediction. We used this property to imple-ment the RK4. Althought the RK4 algorithm is not a PEC algorithm is its traditional formulation, it could be rewritten with strict equivalence in the form of a PEC algorithm with three prediction stages [START_REF] Edward | A modification of the Runge-Kutta fourth-order method[END_REF]. So, the RK4 integrator phases are WAIT, SEND-P1, SEND-P2, SEND-P3, SEND-C, and follows the same sequencing as the Heun's integrator.

Prediction 1 :

       p = h.f (t k , X k , u k ) X p k+1 = X k + 1 2 p q = p Prediction 2 :            p = h.f (t k + h/2, X p k+1 , u k+1/2 ) X p k+1 = X p k+1 + 1 2 (p -q ) q = 1 6 q
Prediction 3 :

           p = h.f (t k + h/2, X p k+1 , u k+1/2 ) - 1 2 p X p k+1 = X p k+1 + 1 2 (p -q ) q = q -p Correction :    p = h.f (t k + h, X p k+1 , u k+1 ) -2p X k+1 = X p k+1 + q + 1 6 p with f (t, X, u) the derivative at time t (f (t, X, u) = A.X + B.u for LTI systems) FIGURE 3 -modified RK4 algorithm
Implementation of QSS 1 and 2. Althought the Artis-star library already offers a QSS1 implementation, we had to re-implement QSS1 and QSS2 to match our integrator mapping architecture described above.

For the QSS1 and QSS2 integrator, we followed the work of Kofman in [START_REF] Kofman | Quantization-based simulation of differential algebraic equation systems[END_REF] with an implementation in two atomic models : one Derivative module and one Integrator-Quantifier (I-Q) module. Although, our Derivative module is supporting multi-dimensional states. As for the I-Q block, we managed to create one simple block which handles multi-dimensional states (with their respective quanta). This allows a clear separation of the tasks : the Derivative block calls the Derivative() function of the model, while the I-Q block calls the ComputeOutput() function.

Results

Simulation with a simple test-case. For testing purpose, we started with a simple a linear second-order system (Model 1) which is a common model of many dynamic processes :

(M1) ẋ = -2 ω 0 x -ω 2 0 (y -sw); (4) (M2) ẏ = x;
With sw being the input signal. This model has the advantage of simplicity while allowing testing the "coupling property" of our integration implementation. We called the "coupling property" the equivalence, in terms of numerical results, between the two following simulation graphs :

-a unique integrator attached to the 2dimensional (M1+M2) model (Figure 4), -two integrators attached each to the 1dimensional M1 and M2 models, coupled by a coordinator (Figure 5). Our first simulation plot (Figure 6) verify the correctness of our QSS1 and QSS2 implementation, in comparison to a reference result We have also checked that the "coupling property" is verified by Euler, Heun, RK4 and QSS1 integrators. This property has a great practical interest : a complex continuous-time system can be hierarchically divided in an arbitrary manner while keeping the same simulation results.

Nevertheless, the coupling property is not verified with our QSS2 implementation. The investigation is still ongoing. The QSS 2 with the slope's approximation of Kofman [START_REF] Kofman | Quantization-based simulation of differential algebraic equation systems[END_REF] is a difficulty for our implementation, as our QSS2 integrators only exchange their state estimates, and not the derivative estimates. A first analysis shows that for the "coupling property" to be verified with QSS of order n ≥ 2, the integrator blocks also have to exchange derivative estimates up to order (n-1).

Simulation with an aircraft dynamics model. We have simulated a more complex model representing the longitudinal control of a Boeing 747 (figure 7). This model includes a linearised flight dynamics model at an altitude of 40,000 ft and cruise speed of 774 ft/s, an altitude PIDcontroller and a thrust PID-controller. The system is split in 7 sub-models and their associated integrators, with a total state dimension of 12.

The testing scenario consists of simulating the aircraft response to a flight-level change during 60s.

For the fixed-step integration algorithms, the steptime is set to 0.2s. For the QSS1 and QSS2, a quantum value has to be chosen for each of the 12 variables that composed the total state of the system. These variables represent heterogeneous physical quantities (rotation rates, angles, velocity, positions) with different orders of magnitude. A first rule of thumb is to set each quantum relatively to the desired accuracy on the associated variable. A second rule of thumb is to choose smaller quanta for inner loops (fast dynamics) than for outer loops (slower dynamics). After a trial-and-error tuning of the quanta, we have obtained the altitude response in Figure 8. Using the RK4 outputs as a reference, the QSS2 gives quite a close result, while the QSS1 inhibits some significant deviation. Also, some instabilities have been observed on pitch angle using QSS2.

All of these QSS1 and QSS2 results are nevertheless highly dependent on the quantum choice.

Efficiency benchmark. The table 1 shows the numbers of integration steps performed by the different algorithms. This number is related to the number of evaluations of the derivative functions, which could be a costly computation on complex models.

Taking into account that the RK4 performs four derivative evaluations per step against only one for QSS, RK4 and QSS1 have a comparable efficiency on the simple test case, though QSS1 still performs twice as many derivative evaluations than RK4. QSS2 needs fewer evaluations as expected. For the more complex B-747 simulation, RK4 appears to be more efficient, but once again, the results for QSS algorithms are highly dependent on the quantum choice. In order to achieve a more consistent evaluation of the accuracy-computation cost ratio, we consider a simple linear first-order model. A set of simulations have been performed, varying the integration time step for RK4 and the (unique) quantum for QSS methods. The meansquare error of the response versus computation time is presented on the figure 8. It shows that the QSS methods, and especially QSS2, can achieve a better accuracy for the same computation cost. Unfortunately, this comparison method cannot be easily extended to higher order models, that needs multidimensional quanta.

Conclusion

Our paper identifies the main methods for distributed and multi-physics simulation for future work at ONERA. On a practical point of view, we focused on Parallel DEVS using the Artisstar library and mapped classical integrator in this formalism. We presented our implementation and results of simulation using Euler, Heun, RK4 and QSS1&2 integration methods in the same framework.

Our implementation ensures that a complex continuous-time system can be hierarchically split into interconnected subsystems, with an integrator associated with each subsystem, while keeping the same simulation results. This property has been verified for Euler, Heun, RK4 and QSS1 integrators and is the key for modu-lar, scalable and distributed simulations of complex continuous-time systems. This work is still ongoing for QSS2.

Concerning the relative efficiency of the evaluated algorithm, we have shown on very simple models that QSS1 and QSS2 can compete favourably with classical fixed time step integration methods, as RK4.

Nevertheless, the application of QSS methods to more complex, real-life aerospace system model has revealed the issue of adjusting a set of quanta to achieve a desired global accuracy while maintaining an acceptable computational cost. Altought theoretical results exist on the relationship between quanta and error bounds, they are not directly usable in practice, and further work needs to be done on this issue.

FIGURE 4 -

 4 FIGURE 4 -Grouped ODE of Model 1

FIGURE 5 -FIGURE 6 -

 56 FIGURE 5 -Coupled ODE of Model 1 decomposed in M1 & M2

FIGURE 7 -FIGURE 8 -FIGURE 9 -

 789 FIGURE 7 -Longitudinal control model of a B-747

TABLE 1

 1

		-Number of integration steps
		{M1+M2)} Coupled Boeing 747
			M1 & M2
	RK-4	502	502	2 999
	QSS1	3 920	3 920	39 450
	QSS2	825	2 763	28 921