Patrice Pajusco
email: patrice.pajusco@telecom-bretagne.eu

Nadine Malhouroux-Gaffet

Relationship between physical and subspace propagation channels Simulation and experimental study

The channel capacity of a MIMO telecommunication system depends on different parameters such as the signal-tonoise ratio, the space-time characteristics of the propagation channel, the antenna array geometry and the antenna pattern. However, the propagation channel remains the strongest constraint imposed on designers of new radio air interfaces. Using a MIMO channel model, computer simulations are performed to highlight the relation between the radiation pattern of subspace channels and the physical space-time characteristics of the propagation channel. This provides a better understanding of the subspace channel origins and achievability of degrees of freedom. In addition to the simulation approach, UWB MIMO channel measurements are used. Both simulation and experimental results show that there is not a bijective relation between physical paths and subspace channels. The number of subspace channels is the same but the amplitudes and directions can be different, leading to inaccuracies in the space-time characteristics of the propagation channel. Nevertheless, the results are promising for fast-jointed DoA-DoD channel estimation.

I.

INTRODUCTION

The channel capacity of MIMO communication systems depends on different parameters, such as the signal-to-noise ratio, the space-time characteristics of the propagation channel, the antenna array geometry and the antenna patterns. Among all of these elements, only the propagation channel cannot be optimized by air interface designers. As Bonek highlights in [START_REF] Bonek | MIMO -Study Propagation First[END_REF], it is essential to thoroughly understand the subjacent phenomena in order to provide reliable models and conclusions. The aim of this article is to investigate the relationship between the physical properties of propagation channels and the subspace propagation channels obtained by the singular decomposition of the MIMO transmission radio channel. This study will allow an improvement in the definition and the understanding of the origin of the maximum capacity of MIMO communication systems.

Simulations were performed to highlight the relationship between the radiation pattern associated with each eigenmode and the space-time characteristics of the propagation channel. This study is important for engineering and simulation studies. Indeed, in a MIMO communication, the power transmitted over subspace channels generates a radiation pattern which is no longer omnidirectional. The prediction of radio coverage and interference is no longer related only to the antenna pattern, but also to the spatial characteristics of the propagation channel. This study will also provide a better physical understanding of the subspace channels' origin, their associated MIMO capacity and the corresponding degrees of freedom of the propagation channel. To complete the simulation study, an experimental MIMO UWB campaign was carried out. It is based on a measurement technique that is able to investigate the 3D spacetime characteristics of a real channel [START_REF] Malhouroux-Gaffet | Effect of propagation phenomena on MIMO capacity of wireless systems between 3 and 10 GHz " EUCAP 2010[END_REF]. Results obtained in an indoor environment will be presented and analyzed.

II. DEFINITIONS

A. Transmission and Propagation Channel

First, it is important to define the transmission channel and the propagation channel. A SISO transmission channel can be seen as a two-port network formed by the two antenna extremities of the radio link. The transfer function can be easily obtained by using, for example, a vector network analyzer. The concept of the propagation channel is more abstract because its characteristics cannot be directly measured. It can be simply defined as the transmission channel without the antenna effects.

B. MIMO Channel

A MIMO system is defined as a transmission system with T n transmitter antennas and R n receiver antennas. In a narrowband context, the propagation channel can be represented by a complex matrix H of dimension R T n n × . Each element , i j h of the matrix H represents the complex transmission gain between the receiving antenna i and the transmitting antenna j. X is the vector corresponding to the transmitted signals applied to the emitter antenna; Y is the vector corresponding to the received signals. In a noise-free configuration, the received signal is expressed as:

Y = HX

(1)

By using singular value decomposition (SVD) the matrix H can be expressed as:

{ } 1 with and =rank( )

b H i i diag b λ = = H = UΣV Σ H (2) 3 ,
,

T n v v v v ⎡ ⎤ = ⎣ ⎦ V : Unitary matrix of T T n n × dimension
containing the vectors columns v k of the input basis.

C. Channel Capacity

The channel capacity concept was introduced by Shannon in 1948 [START_REF] Shannon | A Mathematical Theory of Communication[END_REF]. In the case of a link between two antennas, the capacity in bit/s/Hz is given by equation 3.

With ρ is the signal to noise ratio. This definition highlights two fundamental properties. First, the capacity grows with signal to noise ratio. Secondly, when the signal to noise ratio is large, the increase of capacity will be very slow.

The capacity increase of a MIMO system is based on the data rate of the previous subspace channels. The number of channels is directly related to the number of significant eigenvalues of the matrix H HH . So the larger the number of significant eigenvalues is, the bigger the potential capacity. The total capacity of MIMO system is defined as the sum of the capacity of each subspace channel:

min( , ) 2 1 log 1 T R n n MIMO i i T C n ρ λ = ⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠ ∑ (4)
The capacity will depend on the power allocation strategy on each subspace channel. The optimal solution consists in using the Waterfilling algorithm. However, this solution requires the channel state information at the transmitter which is not always available.

D. Subspace Channel Pattern Radiation

The real and subspace channels in MIMO were recalled in figure 1 in function of X,Y and X',Y'. Considering the K th subspace channel which corresponds to the K th nonzero singular value of the matrix channel H , we have :

k k k k k k x v y u ' ' , Σ = (5)
That means that the K th subspace channel is obtained by weighting the transmitter array by the vector k u' and the receiver array by the vector k v' . So, the resulted subspace channel will have a gain of k k , Σ and will be orthogonal to the other subspace channels. Using a 3D complex antenna pattern of the array, it is easy to find the physical radiation pattern of this subspace channel. Considering that the antenna pattern of the i th element of the transmitter array is A θ ϕ corresponding to the K th subspace channel is given by:

( ) ( ) ( ) 1 , , ' T X n TX T k i k i A a u i θ ϕ θ ϕ = = × ∑ (6) ( ) ( ) ( ) 1 ,
, '

R X n RX R k i k i A a v i θ ϕ θ ϕ = = × ∑ (7)

III. BEHAVIOR OF PHYSICAL AND SUBSPACE CHANNELS

A. Introduction

The matrix H of the transmission channel depends on the propagation phenomena between the transmitter and the receiver, on the characteristics of the antennas and the geometry of the array. There is, a priori, a relation between the physical properties of the propagation channel ie, that is the space-time characteristics of the propagation and the singular value decomposition. Some authors such as Loyka [START_REF] Loyka | Multiantenna capacities of waveguide and cavity channels[END_REF] had studied this relation in the case of waves guides. It was, for example, possible to compare the various propagation modes of the wave guide to the modes generated by a MIMO array set in the guide. These results are interesting because each mode of the wave guide can be considered as a plane wave undergoing some reflections. But this study does not represent a typical radio-mobile environment. From MIMO simulations, deeper investigations will be carried out by the evaluation of the directions of arrival and departure.

B. Channel Simulation Principle

The developed channel simulator uses the double directional approach introduced by Steinbauer in [START_REF] Steinbauer | The double-directional radio channel[END_REF]. This allows simulation of the antennas pattern independently of the propagation channel. Thus, the transmission channel can be known whatever the antenna array geometry. Most of the MIMO channel models are based on this formalism (802.11n, UMTS, Wimax...). The principle consists in modelling the radio channel by a set of multipaths described by their physical characteristics: direction of departure (DoD) direction of arrival (DoA), complex attenuation for each polarization, delay, Doppler. The transmission MIMO channel is then obtained from the characteristics of the multipaths, of the array geometry and the antenna pattern. An efficient and reliable implementation of this model requires optimizations detailed in [START_REF] Conrat | A Versatile Propagation Channel Simulator for MIMO Link Level Simulation[END_REF].

C. Simulation Configuration

The choice of the array geometry is very important. Indeed, most of the MIMO array geometries provide ambiguitie results. For example, the linear array has a conical ambiguity. Thus, two rays coming from different directions are perceived in an identical path by the array and consequently reduce the number of degrees of freedom of the MIMO channel. Circular arrays have been chosen for simulations because they show no ambiguity in azimuth. Simulations are carried out at 5 GHz, corresponding to the band used by WIFI systems. The two sensor arrays consist of 30 sensors vertically polarized and distributed on a circle of 10 cm radius. This configuration corresponds to a sensor spacing of / 3 λ . Initially, in order to simplify the result analysis, the multipaths modelling is restricted to the horizontal plane. Moreover, the channel depolarization is not taken into account.

Even if these simulations are not difficult to perform, it is important to report them in the paper to show some unexpected behaviour in very simple configurations. For each simulation, the geometry of the propagation channel and the radiation pattern corresponding to each singular value at both sides of the link will be represented. The amplitude of the rays is proportional to the length of the ray. For example, figure 2 represents the result of a simulation with only one path. Only one singular value is obtained and the corresponding subspace channel corresponds to a beamforming in the direction of the single path.

D. Effect of the Number of Multipaths

The channel, represented in figures 3a and 3b, consists of 2 paths with different directions and respective amplitudes of 1 and 0.5. Figures 3c andd respectively show the radiation pattern corresponding to each singular value with the transmitter array and the receiver array. The colours correspond to the index of each singular value. There is no, a priori, relation between the colours in the figures a-b and the colours in the figures c-d. The decomposition in singular values allows the identification of 2 nonzero values. When the paths are sufficiently spaced, the number of singular values and the number of propagation paths are identical. The decomposition maintains the relation between direction of departure and direction of arrival of each path. Each singular value carries a particular path and the amplitude of the singular value is approximately proportional to the power of the path. This bijection between each singular value and a physical path is not always checked. For a same geometrical configuration, the variation of the amplitude or phase of the path can modify the decomposition significantly. A singular value can be linked to several paths. For example, figure 4 gives results considering the same paths model with identical path amplitudes. Each singular value exploits simultaneously the two propagation paths. The subspace radiation patterns beam in 2 directions but not exactly in the real physical direction. The directions of pointing of each singular value present opposite errors. Consequently, the averaged radiation pattern of the singular values is in agreement with the physical paths. This remark is also applicable to the amplitude.

• 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • c) TX 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • d) RX
• 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • c) TX 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • d) RX
• 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • c) TX 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • d) RX
All these results remain valid with an unspecified number of paths but only if the paths are independent. This condition requires having paths sufficiently angularly spaced in comparison to the angular resolution of the array. This condition applies on both sides of the link and to the whole of the paths. If we consider a model with 3 paths of identical amplitudes and with two of them having close directions of arrival. The result illustrated in figure 5, shows that only two eigenvalues are corresponding to the major part of energy.

• 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • c) TX 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • d) RX

E. 3D Generalization with Polarisation

In order to generalize the previous results, it is essential to consider multipaths in 3D with arbitrary polarization. This study was achieved by replacing each element of the UCA by a ULA of 5 elements. Each element is built of two isotropic antennas in cross polarization. This provides a cylindrical array of 300 antennas with a geometrical shape represented in figure 6. Now, the simulation is taken into account the polarization and the path loss is now expressed by a complex matrix of dimension 2 2 × instead of a simple complex value. So the depolarization of a path can be defined by its matrix. In free space configuration with only one path, the MIMO channel has two singular values with the same amplitude. The number of singular values is doubled. This is due to the independence of two polarizations. The radiation pattern of the singular mode corresponding to each singular value is presented in figure 7. They follow the physical direction of the path. The blue segments, indicating the trajectory of the field E , show the orthogonality of the polarization between the two subspace radiation patterns. In the case of depolarization, the number of eigenvalues remains equal two. However, the amplitudes can be very different according to the matrix conditioning of the path. For example, if we consider a channel depolarization of approximately 10%, the module of the power gain is given by 0.9 0.1 0.1 0.9

⎡ ⎤ ⎢ ⎥ ⎣ ⎦
, and the phase is supposed to be random. The simulation results are represented in figure 8. The singular mode radiation patterns correspond quite well to the physical directions of the paths. However, the energy is no longer the same between the two subspace channels, and the polarization is arbitrary. The first eigenvalue contains most of the energy. Other simulations with different multipaths were also carried out. Even without depolarization, multipaths induce an imbalance of the eigenvalue amplitudes. Apart from the previous remarks, the conclusions obtained in 2D remain valid in 3D.

F. Conclusion

From the development of a propagation channel simulation tool including antenna arrays and the propagation channel, the relation between subspace and physical channels was investigated. All of these simulation results were obtained with circular and cylindrical MIMO arrays. This study shows that the number of eigenvalues corresponds to the number of independent paths. Independence can be defined by an angular spacing higher than the resolution of the antenna array. Apart from some particular cases, there is no "one to one" relationship between the radiation pattern of a singular value and the direction of single path. In fact, a singular value can beam simultaneously in different directions, not corresponding exactly to the physical direction. The number of eigenvalues is doubled when polarization is taken into account. The eigenvalue amplitudes are identical in the case of line-of-sight (LOS) configuration; otherwise they are usually different because of the depolarization induced by the antennas and the channel.

IV. EXPERIMENTAL SUBSPACE CHANNEL IN REALISTIC ENVIRONMENT.

A. Description of the Experimentation

Many experimental studies have been carried out to model the MIMO propagation channel. The most advanced experiments use double directional measurements. However, the complexity is so high that these experiments are simplified and provide only partial information of the channel (3D or polarization is missing for instance). This complexity is further increased in UWB configurations. A good compromise between accuracy and complexity was obtained with an original method presented in [START_REF] Pajusco | Novel array structure for space-time characterization of the UWB channel[END_REF]. This device is able to study the degree of freedom of the channel thanks to a large number of sensors and the non-ambiguity of the DoA and DoD.

B. Results in anechoic chamber

Experimentation was carried out in an anechoic chamber (figure 9). The distance between the two antenna arrays was 4.5 m. The far-field condition was met from 3 GHz to 10 GHz. Only one path is obtained with this experimental configuration. The degree of freedom is defined as the number of highest singular values whose sum is equivalent to 90% of the total power. Only two degrees of freedom were obtained by analyzing the results, whatever the frequency. This result is in agreement with the simulation results. Figure 10 shows that the directions of the subspace channel are in agreement with the position of the antennas. However, the amplitudes are not perfectly identical, and polarizations are not perfectly orthogonal. This is probably due to antenna depolarization.

Tx-SV1

Rx-SV1

Tx-SV2 Rx-SV2

Tx-SV2 Rx-SV2

C. Evaluation in a actual environment

Experimentation was performed in a meeting room in order to investigate the number of degrees of freedom in actual configuration. A 360° panoramic picture, shot from the emitter location and presented in figure 11, shows the environment. Contrary to the anechoic chamber, this LOS configuration shows multipaths. The number of degree of freedom is larger. Moreover, the number of degrees of freedom increases with the frequency contrary to their behavior in the anechoic chamber. This is due to the improvement of the MIMO array resolution when the frequency increases. This characteristic allows the identification of additional paths and thus increases degrees of freedom. An analysis shows an increase in the number of eigenvalues with the delays. This is explained by the increase of propagation phenomena when the delay increases. At the opposite, the direct path has only 2 eigenvalues.

V. CONCLUSION MIMO channel simulations using circular and cylindrical arrays highlighted the relations between physical and subspace channels. The number of singular values of the transmission channel corresponds to the number of independent paths. The independence of the paths depends on the angular resolution of the MIMO array and on the DoA and the DoD of the paths. There is no bijective relation between the subspace channels and the multipaths. This investigation was also carried out by experimentation. An original MIMO UWB device was used in an anechoic chamber and in a real environment. An analysis concerning the evolution of the eigenvalues according to the frequency and the environment was completed..Future study will investigate other representative environments, such as residential environment.

  columns u k of the output basis. Σ is the diagonal matrix containing the singular values sorted by descending order. The singular values are positive or null. The nonzero values correspond to the eigenvalues of the matrix H HH , identical to the nonzero eigenvalues of H H H . Considering ' the representation of Y in the output basis, the simplified relation ' Y' = ΣX is obtained. Σ is diagonal, so the MIMO transmission is equivalent to the transmission of parallel links, whose amplitude of each subspace channel K is given by k k , Σ .

Figure 1

 1 Figure 1 shows the two representations respectively Y = HX and Y' = ΣX'

Figure 1 .

 1 Figure 1. Physical and subspace representation of a MIMO transmission

  with the complex components e θ and e ϕ of the radiation directions θ and ϕ . Theses diagrams include the spatial configuration of the array. We can also similarly define the radiation pattern ( ) , RX i a θ ϕ of the antenna pattern of the i th element of the receiver array. The relation between the antenna pattern ( )

Figure 2 .

 2 Figure 2. Subspace channel with one path configuaration.(a) and (b) represent respectively the arrival ray and the corresponding departure ray. (c) and (d) represents respectiveley the patterns of the arrival eigenvalue and the corresponding departure eigenvalue.

Figure 3 .

 3 Figure 3. Subspace channel with two paths configuration with different amplitudes. (a) and (b) represent respectively the 2 arrival rays and the corresponding 2 departure rays. (c) and (d) represents respectiveley the patterns of the 2 arrival eigenvalues and the corresponding departure eigenvalues.

Figure 4 .

 4 Figure 4. Subspace channel with two paths configuration with same amplitude. (a) and (b) represent respectively the 2 arrival rays and the corresponding departure rays. (c) and (d) repesents respectiveley the patterns of the arrival eigenvalues and the corresponding departure eigenvalues.

Figure 5 .

 5 Figure 5. Subspace channel with three paths configuration with same amplitude. (a) and (b) represent respectively the 3 arrival rays and the corresponding 3 departure rays. (c) and (d) repesents respectiveley the patterns of the arrival eigenvalues and the corresponding departure eigenvalues.

Figure 6 .

 6 Figure 6. Cylindrical array. Let us consider a path having the following geometrical characteristics: 135 departure θ

Figure 7 .

 7 Figure 7. Subspace radiation patterns with one path.

Figure 8 .

 8 Figure 8. Subspace radiation patterns with one path and depolarization.

Figure 9 .

 9 Figure 9. Anechoic chamber environment

Figure 11 .

 11 Figure 11. Meeting room environment


Y Feng

M Tayyab

P Boivin

A Lattice-Boltzmann model for low-Mach reactive flows

Keywords: lattice Boltzmann, reacting flows, combustion

A new Lattice-Boltzmann model for low-Mach reactive flows is presented.

Based on standard lattices, the model is easy to implement, and is the first, to the authors' knowledge, to pass the classical freely propagating flame test case as well as the counterflow diffusion flame, with strains up to extinction.

For this presentation, simplified transport properties are considered, each species being assigned a separate Lewis number. In addition, the gas mixture is assumed to be calorically perfect. Comparisons with reference solutions show excellent agreement for mass fraction profiles, flame speed in premixed mixtures, as well as maximum temperature dependence with strain rate in counterflow diffusion flames.

Introduction

Industries from the aerospace, aeronautic and automotive sectors are increasingly relying on numerical simulation tools. From the occasional use of a research and development department, these tools progressively made it to conception and production departments, where they help to continuously improve designs. The field of low-Mach external aerodynamics and aeroacoustics, have been particularly impacted by the rapid development of Lattice-Boltzmann (LB) methods [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] in the last five to ten years. From industrial benchmarks [START_REF] Marié | Comparison between lattice boltzmann method and navier-stokes high order schemes for computational aeroacoustics[END_REF][START_REF] Casalino | Lattice-Boltzmann aeroacoustic analysis of the lagoon landing-gear configuration[END_REF][START_REF] Sengissen | Simulations of lagoon landing-gear noise using lattice boltzmann solver[END_REF], these methods quickly ramped up to full scale applications: full-scale cars [START_REF] D'hooge | Application of real-world wind conditions for assessing aerodynamic drag for on-road range prediction[END_REF][START_REF] Gleason | Comparison of computational simulation of automotive spinning wheel flow field with full width moving belt wind tunnel results[END_REF], full-scale aircraft engines [START_REF] Casalino | Turbofan broadband noise prediction using the lattice boltzmann method[END_REF] and even fullscale aircrafts [START_REF] Khorrami | Computational evaluation of airframe noise reduction concepts at full scale[END_REF][START_REF] Khorrami | Simulation-based airframe noise prediction of a full-scale, full aircraft[END_REF], oftentimes with outstanding results.

Combustion modelling in the LB framework, however, remain relatively marginal within the combustion community. A hybrid finite difference lattice Boltzmann model has been presented for the simulation of low Mach number flows with significant density changes by Filippova et al [START_REF] Filippova | A Novel Lattice BGK Approach for Low Mach Number Combustion[END_REF] almost twenty years ago. An LB model for combustion modelling was then presented by Yamamoto and coworkers [START_REF] Yamamoto | Simulation of combustion field with lattice Boltzmann method[END_REF][START_REF] Yamamoto | Combustion Simulation using the Lattice Boltzmann Method[END_REF][START_REF] Yamamoto | LB simulation on soot combustion in porous media[END_REF] for steady and unsteady reactive flows. About ten years ago, a pressure-based lattice Boltzmann model was presented by Chen et al. for low Mach number combustion simulations by [START_REF] Chen | A novel coupled lattice Boltzmann model for low Mach number combustion simulation[END_REF][START_REF] Chen | A simple lattice Boltzmann scheme for combustion simulation[END_REF], based on an incompressible LB model developed earlier by Guo et al [START_REF] Guo | Lattice BGK Model for Incompressible Navier -Stokes Equation[END_REF]. An entropic lattice Boltzmann model on two-dimensional standard lattice was also presented for compressible thermal flows and extended to combustion applications [START_REF] Prasianakis | Lattice Boltzmann method for thermal flow simulation on standard lattices[END_REF][START_REF] Rienzo | Lattice boltzmann model for reactive flow simulations[END_REF]. More recently, detailed kinetics were successfully implemented in the LB framework [START_REF] Chiavazzo | Coupling of the model reduction technique with the lattice Boltzmann method for combustion simulations[END_REF][START_REF] Chiavazzo | Efficient simulations of detailed combustion fields via the lattice Boltzmann method[END_REF].

Thermodynamic closure is one of the key issues in extending the LB capabilities from low-Mach aerodynamics and aero-acoustics to reactive flows, and obtaining satisfactory results in cases with significant thermal expansion is challenging, as indicated by the low-number of combustion models in the LB context. Most of LB models mentioned above for reactive flows were specifically applied under constant pressure, and lack full coupling between thermodynamics and the flow. In particular, because their equilibrium density distribution function only depends on the local density, hydrodynamic pressure and velocity, and not on temperature, these models may fail at setting in motion fluids at rest when heat-release is applied. This aspect is however critical in representing canonical combustion phenomena, such as ignition in a mixing layer, or the dynamics of freely propagating flame in a premixed mixture.

To address this issue, more complex, multi-speed models with aid of finite volume approach were developed for subsonic and supersonic flows [START_REF] Lin | Double-distribution-function discrete boltzmann model for combustion[END_REF], including detonations. Alternatively, coupled lattice Boltzmann models on standard lattices have been investigated in low Mach thermal compressible flows [START_REF] Li | Coupling lattice boltzmann model for simulation of thermal flows on standard lattices[END_REF][START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF].

The model presented in this work is based on the latter [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF] and takes full advantage of the simplicity and practicality of standard lattice models. It is able to deal with multi-component, calorically perfect gas mixtures, and the coupling between the thermodynamics and the flow is two-way. Simplified transport properties are considered, with non-unity Lewis number set for each component. This new LB model for low-Mach reactive flows is presented in the first Section. The second Section validates the model using Cantera [START_REF] Goodwin | Cantera: An objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] as a reference on canonical combustion test cases of premixed and nonpremixed combustion.

Lattice-Boltzmann combustion modelling

Introductory comments: athermal LB method Lattice Boltzmann methods intend to solve the continuous Boltzmann equation through discretization of space, time, and velocity [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF][START_REF] Qian | Lattice bgk models for navierstokes equation[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF]. In the so-called DnQm model, a n-dimensional physical space is filled with a regular lattice (or a Cartesian grid) and velocity space is discretized on a set of lattice tensor E n m = [c 0 , ..., c i , ..., c m-1 ] . On every lattice node x, f i (x, t) denotes the density distribution of particle with velocity c i , thus the local density ρ and momentum ρu are defined as ρ = Σ i f i and ρu = Σ i c i f i .

A single time relaxation process is used to model the collision term of Boltzmann equation in this study, which is the so-called Bhatnagar-Gross-Krook (BGK) model [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral onecomponent systems[END_REF][START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF].

f i (x + cδt, t + δt) -f i (x, t) = - 1 τ [f i (x, t) -f eq i (x, t)] ( 1
)
where τ is the relaxation parameter, δt is the time increment and δt = δx/c. f i (x, t), f i (x + cδt, t + δt) are the distribution functions associated with the i th discrete velocity c i , and f eq i is the i th equilibrium distribution function given by [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF][START_REF] Qian | Lattice bgk models for navierstokes equation[END_REF]

f eq i = ρw i 1 + c i • u c 2 s + (c i • u) 2 2c 4 s - u 2 2c 2 s ( 2
)
where w i is the weight coefficient associated to discrete velocity c i and c s is sound speed.

Using the Chapman-Enskog multiscale technique [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF], the equivalent Navier-Stokes equations can be recovered as

∂ρ ∂t + ∇ • (ρu) = 0, (3)
∂ρu ∂t + ∇ • (ρuu) = -∇p + ∇ • ρν(∇u + (∇u) T ) , (4)
where the kinematic viscosity ν is related to the relaxation parameter through

τ = ν c 2 s δt + 1 2 . ( 5
)
And, the pressure is related to the density by the equation of the state p = ρc 2 s . The equation of state used in this base model is clearly not suitable for reactive flows, as no temperature is considered. The most simple acceptable equation of state for reactive flows is probably the ideal gas law, relating the gas pressure and temperature through

p = ρrT, (6)
where r = R/W , with R the gas constant, and

1/W = k Y k /W k ,
f eq i = f (0) i = ρw i 1 + c iα u α c 2 s + A (2) αβ Q (2) iαβ 2c 4 s + A (3) αβγ Q (3) iαβγ 6c 6 s , A (2)
αβ = u α u β + (θ -1)c 2 s δ αβ , Q (2)
iαβ = c iα c iβ -c 2 s δ αβ , A (3)
αβ = u α u β u γ + (θ -1)c 2 s [uδ] αβγ , Q (3)
iαβ = c iα c iβ c iγ -c 2 s [cδ] αβγ , (7)
where

[cδ] αβγ = c α δ βγ + c β δ αγ + c γ δ αβ , δ αβ is the Kronecker symbol and θ is the non-dimensional temperature θ = rT c 2 s = RT c 2 s k Y k W k (8)
The moment of equilibrium distribution function on the nearest neighbor type lattices, hereafter referred to as standard lattices (D2Q9, D3Q19, D3Q27, etc.), are then

i f (0) i = ρ, (9)
i f (0) i c iα = ρu α , ( 10
) i f (0) i c iα c iβ = ρrT δ αβ + ρu α u β , (11)
i f (0) i c iα c iβ c iγ = ρrT (u α δ βγ + u β δ γα + u γ δ βα ) + ρu α u β u γ + Ψ αβγ , (12)
where Ψ αβγ is a deviation term due to defect of symmetry of standard lattices for the third order moment. This deviation can be corrected by introduction of an extra force term s i = Q iαβ ∂ ∂xγ Ψ αβγ in the lattice Boltzmann equation. The LB equation with a force term in the second order accuracy scheme is given as [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice boltzmann method[END_REF]

f i (x + cδt, t + δt) = f i (x, t) -1 τ [f i (x, t) -f eq i (x, t)] (13)
+(1 -1 2τ )s i
where the external force term s i reads [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF]

s i = Q ixx ∂ ∂x [ρu x (1 -θ -u 2 x )] + Q iyy ∂ ∂y [ρu y (1 -θ -u 2 y )], (14)
in two dimensions.

The mass and momentum conservation equations derived through Chapman-Enskog expansion from this model remain unchanged compared to the classical athermal version (3) and ( 4), with the exception of the equation of state and the relation between the relaxation time and the viscosity [START_REF] D'hooge | Application of real-world wind conditions for assessing aerodynamic drag for on-road range prediction[END_REF], which now reads

τ = ν θc 2 s δt + 1 2 . ( 15
)
The approach proposed here is hybrid: coupled with this lattice Boltzmann description, temperature T and mass fractions Y k equations are solved following a classical finite difference method

∂T ∂t + u α ∂ ∂x α T = 1 ρ ∂ ∂x α (ρD T ∂T ∂x α ) + ω h ρc p (16)
∂Y k ∂t + u α ∂ ∂x α Y k = 1 ρ ∂ ∂x α (ρD k ∂Y k ∂x α ) + ω k ρ (17)
D k and D T are respectively the k th species and thermal diffusivities.

As to simplify the introduction of this model, we assume calorically gases, with a constant mixture heat capacity. Transport properties are specified through constant Prandtl number and constant Schmidt number for each species following [30]

D T = ν P r , D k = ν Sc k . ( 18
)
Physical and lattice units for length, time, mass and temperature are related through a reference length scale L 0 , the physical sound speed c sp for space and time, a reference density ρ 0 for mass, and a reference temperature T 0 . Accordingly, the physical variables can be converted from the quantities in lattice units (l subscript) as

δx = L 0 /N i , δt = δx. c s c sp , t = N t .δt, u = u l . δx δt , (19)
µ = µ l ρ ρ l δx 2 δt
where c sp is the reference physical sound speed. It could be set to the fresh gases sound speed, but in low-Mach flows, it is more convenient to keep it as a free parameter. This allows to freely accelerate the convergence rate, by varying the time step somewhere between a velocity based time-step (as in low-Mach solvers) and a sound-speed based time-step (as in compressible solvers). We recommend to adopt a high reference temperature (here T 0 = 3000K) as to restrict the variations of θ -1 (let us remind that θ = rT /c 2 s ) and improve numerical stability [START_REF] Li | Three-dimensional non-free-parameter lattice-Boltzmann model and its application to inviscid compressible flows[END_REF][START_REF] Feng | A compressible lattice Boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF].

Model implementation

A strong advantage of the present model is its simplicity. The algorithm we implemented is illustrated in Fig. 1. In essence, two solvers are simultaneously running:

• the Lattice Boltzmann (LB) solver, for mass and momentum conservation Eqs (3) and (4)

• the Finite Differences (FD) solver, for temperature T and mass fractions Y k evolution through Eqs ( 16) and ( 17)

Start Initialization ρ, u, Y k , T Compute f eq i (ρ, u, θ) (7) s i (ρ, u, θ) (14) τ (ρ, θ) (15)
Streaming and collision step [START_REF] Yamamoto | LB simulation on soot combustion in porous media[END_REF] Update (ρ, u) Compute right handside of Eqs. ( 16) and ( 17)

Runge Kutta integration

Update (Y k , T )

Output LB FD time-step At the end of each time-step, the two solvers exchange macroscopic values:

f i = f eq i (ρ, u, θ) (Y k , T ) ( ρ , u ) θ = f ( Y k , T ) ( 8 )
(ρ, u) is sent from the LB to the FD solver, and the gas composition and temperature (T, Y k ) from the FD to the LB solver.

The LB solver, being of very classical form as well as built on standard lattices, can share most of the features of standard LB methods. In particular, a regularized BGK collision model was employed to enhance the robustness [START_REF] Latt | Lattice boltzmann method with regularized precollision distribution functions[END_REF][START_REF] Jacob | A new hybrid recursive regularized bhatnagar-gross-krook collision model for lattice-boltzmannmethod based large-eddy simulation[END_REF], albeit not strictly necessary for the test cases proposed here. The gradients in the correction term ( 14) are numerically implemented by an isotropic second order central difference scheme [START_REF] Thampi | Isotropic discrete laplacian operators from lattice hydrodynamics[END_REF][START_REF] Kumar | Isotropic finite-differences[END_REF].

The FD solver's implementation is also very straight-forward. Equations ( 16) and ( 17) are integrated through a second-order Runge-Kutta scheme in time and the same isotropic finite difference scheme was adopted for the convective and diffusive terms [START_REF] Thampi | Isotropic discrete laplacian operators from lattice hydrodynamics[END_REF][START_REF] Kumar | Isotropic finite-differences[END_REF].

Boundary conditions

A non-equilibrium bounce back scheme with a regularization procedure [START_REF] Malaspinas | General regularized boundary condition for multi-speed lattice Boltzmann models[END_REF] is employed at the inlets for the LB solver, allowing to set all velocity components and the pressure (or density). (T, Y k ) are hard-coded at inlets for the FD solver.

Outlets are treated as zero-gradient boundaries on microscopic distribution functions f i as well as T and Y k variables.

Validation through canonical combustion problems

This section aims at validating the combustion model developed in a D2Q9 Lattice-Boltzmann framework against solutions from Cantera [24], a

well-validated open source tool designed to compute quickly and efficiently.

It incorporates and executes fundamental thermodynamics, chemical kinetics and species transport calculations.

Premixed combustion: freely propagating flame

Let us consider a multi-component gas consisting of

C 3 H 8 , O 2 , H 2 O, CO 2 , N 2 ,
with Schmidt numbers listed in Tab. 1, along with the parameters from the viscosity's power-law,

µ = µ 0 T T 0 β , (20)
and the Prandtl number, defining now entirely the species and thermal diffuses properties through [START_REF] Rienzo | Lattice boltzmann model for reactive flow simulations[END_REF]. In the above T 0 = 298K is the reference temperature of the viscosity power-law, not to be confused with the LB reference temperature (3000K here). The mixture's heat capacity is assumed Having introduced these kinetic, thermodynamic and transport properties, let us now consider the most popular laminar flame, the so-called freely propagating flame. A pseudo one-dimensional domain is now initialized as

µ 0 1.
T, p, u, Y C 3 H 8 , Y O 2 , Y CO 2 , Y H 2 O =        300K, 1atm, 0, 0.06, 0.21, 0, 0 , if x < L/2
3070K, 1atm, 0, 0, 0, 0.18, 0.1 , else.

,

the fresh gases corresponding to a stoichiometric mixture of propane and air in ambient conditions (1atm., 300K), and the burnt gases being set to the adiabatic temperature at ambient pressure. The right side consists of the combustion products of the fresh gases combustion, with 3070K being the adiabatic flame temperature. Let us remind here that we are considering constant C p and a reduced number of species, which both contribute to the rather high burnt gas temperature.

While Cantera [START_REF] Goodwin | Cantera: An objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] includes an straight-forward automated procedure for computing the flame speed, our LB model does not. Computations were initialized as a Heaviside function between the fresh and burnt states. The velocity U f is imposed at the inlet, with U f set arbitrarily, given the flame does not leave the computational domain before convergence.

After initial establishment of the flow, the flame speed can be estimated through

S L = U f - ρ f U f -ρ b U b ρ f -ρ b , (22)
where the f subscript indicate the fresh gases, and b denotes the burnt gases properties (respectively the first and last points of the simulation domain).

Mass fractions and temperature profiles in the stoichiometric premixed flame are presented in Fig. 2, showing excellent agreement between our method and the Cantera-computed reference. x(m) T (K) As an additional validation, the equivalence ratio dependence of the flame speed, estimated through Eq. ( 22) is presented in Fig. 3 and compared with the Cantera reference value. The left side of the plot indicates convergence histories, through the time evolution of the velocity evaluated as [START_REF] Li | Coupling lattice boltzmann model for simulation of thermal flows on standard lattices[END_REF]. For most computations, U f = 0, but 0.4 and 0.8m/s were also tested at ϕ = 1, as to check the consistence of the model. As indicated by the convergence histories, the obtained flame speed is independent of the inlet velocity U f and does indeed converge to Cantera's predictions. Equivalence ratio Flame speed (m/s) time (s) Figure 3: Freely propagating flame: Variation with equivalence ratio of the flame speed in propane/air mixtures at 300K and atmospheric pressure. Cantera reference (plain dark line), and present model (dashed colored line). On the secondary x-axis are plotted the convergence histories for the LBM computed flame speed S L estimated through Eq. [START_REF] Li | Coupling lattice boltzmann model for simulation of thermal flows on standard lattices[END_REF]. U f = 0 for all ϕ = 0.5, 0.6, 0.7, 0.8, 0.9 and 1 (thin lines). Dotted and dashed lines histories for ϕ = 1 correspond to U f = 0.4 and 0.8, respectively.

Mass fractions

Non-premixed combustion: counter-flow diffusion flame

As a concluding validation, counter-flow diffusion flame is now considered. This test case, consisting of two counter flows (oxidizer and fuel), can in principle be solved as a 1D problem with an adequate choice of variables [START_REF] Liñán | The asymptotic structure of counterflow diffusion flames for large activation energies[END_REF][START_REF] Williams | Combustion Theory[END_REF], as done in Cantera [START_REF] Goodwin | Cantera: An objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF]. For validation purposes, however, it is sometimes convenient to solve a 2D problem. A two-dimensional rectangular L × l is adopted as geometry domain for our LBM computation. Top and bottom boundaries are inlet, whereas the left boundary is a symmetry, and the right is a pressure outlet. Kinetic, thermodynamic and transport properties are identical to those of the premixed flame presented above.

The fuel mixture with uniform density ρ f , speed U f and temperature T f = 300K is injected from the upper boundary at atmospheric pressure.

On the opposite inlet, air is uniformly injected at ρ ox , U ox , T ox = 300K and atmospheric pressure. Fuel and oxidizer mixtures have the same N 2 mass fraction Y N 2 = 0.767, the remainder being propane and dioxygen, respectively. The set up for L/l = 2 is illustrated in Fig. 4 through the superposition of streamlines and temperature profiles, for mass flow rates

ρ f U f = ρ ox U ox = 0.2kg m 2 s -1 . The inter-nozzle distance is L = 20mm,
corresponds to a strain of 16.5s -1 . It is interesting to note that no specific treatment was found necessary in the inlet/outlet corners, whereas with more conventional methods authors typically extend the domain as to have an outlet far from the inlet [START_REF] Shum-Kivan | Simulation des grandes echelles de flammes de spray et modélisation de la combustion non-prémélangée[END_REF].

Mass fractions, axial velocity and temperature profiles are reported in Fig. 5, showing an outstanding agreement. The error in maximum temperature is less than 1K. Somewhat larger errors can be seen on the axial velocity profile, to be attributed to the fact that the two-dimensional setup is not strictly equivalent to the 1D problem solved in Cantera.

It becomes then of interest to check whether the dependence of the maximum temperature with the strain rate is correctly reproduced with the LB model. The 2D configuration leading to an unstable stagnation line when the domain is large, simulations were ran on a thinner domain, with L/l = 100.

Results, albeit slightly altered in accuracy remain close to the reference, with an error of 15K in maximum temperature in the conditions of Fig. 5. The dependence of the maximum temperature with strain rate in this thin 2D domain are reported in Fig. 6, showing a very good agreement.

The × symbol at the lowest strain rate corresponds to the result on the large domain simulation (L/l = 2), showing the errors to come from the simulation domain rather than the model. The strain rate at extinction is 305s -1 with our LB model, whereas Cantera sets it at 310s -1 . These results conclude the validation of our flow model.

Concluding remarks

A new LB combustion model has been presented. It is, to the knowledge of the authors, the first to pass the standard canonical combustion test cases of freely propagating premixed flame and counter-flow diffusion flame with strains up to extinction. In practice, the number of independent variables is the number of lattice discrete velocities (9 in our D2Q9 implementation), plus the temperature, and a N -1 or N mass fractions (N being the number of chemical species), as is typical in combustion modelling [START_REF] Poinsot | Theoretical and numerical combustion[END_REF].

The main assumptions are:

• low Mach number

  linking the mean molecular weight W to the species molecular weights W k and the species mass fractions Y k .A hybrid thermal LB modelConsidering the third order Hermite expansion of Maxwell-Boltzmann distribution, one can obtain the equilibrium distribution function in discrete Gauss-Hermite space[START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the navier-stokes equation[END_REF]

Figure 1 :

 1 Figure 1: The algorithm proposed consists of a Lattice-Boltzmann (LB) solver coupled with a classical Finite Differences (FD) solver. Data exchanges between the two solvers are clearly identified.

Figure 2 :

 2 Figure 2: Freely propagating flame: Temperature profile (thick line), C 3 H 8 (△), O 2 (▽), CO 2 (•) and H 2 O ( ) mass fractions. The fresh gases consist of a stoichiometric mixture of propane and air, at 300K and atmospheric pressure. Cantera reference (plain line), and present model (dashed line).

Figure 4 :

 4 Figure 4: Temperature profiles in a counterflow diffusion flame (dotted lines), superimposed with streamlines (dashed lines). Fuel mixture (Y C3H8 = 0.21, Y N2 = 0.79) is injected from the top nozzle, and air from the bottom. Conditions:ρ f U f = ρ ox U ox = 0.2kg m 2 s -1 , T f = T ox = 300K, atmospheric pressure.

Figure 5 :

 5 Figure 5: Counter flow diffusion flame.ρ f U f = ρ ox U ox = 0.2kg m 2 s -1 (a) C 3 H 8 (△), O 2 (▽), CO 2 (•)and H 2 O ( ) mass fractions profiles (b). Temperature (thick lines) and axial velocity profiles, . Atmospheric pressure, temperature of fresh gases 300K. Cantera reference (plain dark line), and present model (dashed line).

Figure 6 :

 6 Figure 6: Maximum temperature in a counter flow diffusion flame as a function of the strain rate. Reference computed with Cantera (line), and our LB model with L/l = 100 (+), or L/l = 2 as in Figs. 4 and 5 (×).

Table 1 :

 1 Power-law's viscosity coefficients (SI units), Prandtl number and Schmidt numbers for each species.Last but not least, a one-step Arrhenius kinetic model is introduced,following global reaction C 3 H 8 + 5 O 2 -→ 3 CO 2 + 4 H 2 O, associated with the kinetic rate ω = k.C C 3 H 8 .C O 2 .e -Ea/RT , where k = 9.9×10 13 cm 3 mol -1 s -1 , E a = 30kcal/mol, and C i are the molar concentrations for species i.

		782 × 10 -5
	β	0.690
	P r	0.682
	S C 3 H 8 1.241
	S O 2	0.728
	S CO 2 0.941
	S H 2 O 0.537
	S N 2	0.690

to be constant C p = 1008J/K/kg.

Acknowledgements B. Cuenot, P. Sagaut and L. Vervisch are gratefully acknowledged for fruitful discussions. We also acknowledge -for supporting Y. Feng at M2P2 -the project CLIMB, with the financial support of BPIFrance (Project No. P3543-24000) in the framework of the program "Investissement d'Avenir: Calcul Intensif et Simulation Numérique". This work was granted access to the HPC resources of Aix-Marseille Université funded by the project Equip@Meso (ANR-10-

• simplified transport properties, specified by species individual Schmidt numbers • constant mixture heat capacity.

The last hypothesis, adopted for the sake of readability, is the most limiting in practice and is to be relaxed in our next work in preparation.