Luis F Contreras-Samamé

Salvador Domínguez-Quijada

Olivier Kermorgant

Philippe Martinet

CoMapping: Multi-robot Sharing and Generation of 3D-Maps applied to rural and urban scenarios

We present an experimental study for the generation of large 3D maps using our CoMapping framework. This framework considers a collaborative approach to efficiently manage, share, and merge maps between vehicles. The main objective of this work is to perform a cooperative mapping for urban and rural environments denied of continuous-GPS service. The study is split in to 2 stages: Pre-Local and Local. In the first stage, each vehicle builds a Pre-Local map of its surroundings in real-time using laser-based measurements, then relocates the map in a global coordinate system using just the low cost GPS data from the first instant of the map construction. In the second stage, vehicles share their pre-local maps, align and merge them in a decentralized way in order to generate more consistent and larger maps, named Local maps. To evaluate performance of all the cooperative system in terms of map alignments, tests are conducted using 3 cars equipped with LiDARs and GPS receiver devices in urban outdoor scenarios of the École Centrale Nantes campus and rural environments.

I. INTRODUCTION

Urban and rural outdoor scenarios can be a real challenge to perform mapping, since environment complexity, such as terrain roughness or lack of structure or dimensions can affect negative the task performance. Other aspects as access to unexplored area, dimension of the region or communication constraints can be considered in the process. In the case of constructing maps for large areas, the idea of using a set of vehicles for building accurate maps in a reasonable amount of time is feasible [START_REF] Dinnissen | Map merging of multi-robot slam using reinforcement learning[END_REF], because a cooperative mapping extends the capability of a single robot by sharing and merging data between group members. When all the data is analysed and merged in a single computation unit, the process is called centralized multi-robot mapping [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF]. On the opposite, the goal here is to have each mobile unit build its own map and merge them upon rendezvous, which is a decentralized process [START_REF] Özkucur | Supervised feature type selection for topological mapping in indoor environments[END_REF], [START_REF] Dinnissen | Map merging of multi-robot slam using reinforcement learning[END_REF], [START_REF] Zhang | Aerial and Ground-based Collaborative Mapping: An Experimental Study[END_REF].

This approach is considered in this article in order to obtain the 3D map of an urban environment. Two vehicles, the ZOE and FLUENCE, work independently during the exploration and define a meeting point to allow direct exchange of data (such as pose, size, limits and maps) during the map reconstruction. The mapping of a rural scenario was also considered, in which was used Lidar-measurement data performed by GOLFCAR robot. The mapping was based on multi-robot approach because the GOLFCAR path was split in two, simulating two robots for the task (see Fig. 1) Luis Contreras, Salvador Domínguez and Olivier Kermorgant are from the Laboratoire des Sciences du Numérique de Nantes (LS2N), École Centrale de Nantes (ECN), 44300 Nantes, France. Emails: {name.surname}@ls2n.fr.

Philippe Martinet is with INRIA Sophia Antipolis, 06902 Sophia Antipolis, France Email: Philippe.Martinet@inria.fr This work was supported by the Erasmus Sustain-T and ALFS project. Parts of the devices used here were funded by the project ROBOTEX, reference ANR-10-EQPX-44-01. • In the first stage named "Pre-Local Mapping", each platform constructs its own map by processing range measurements from a 3D LiDAR for a six degrees of freedom (6-DOF) movement and using GPS data (GPS/GGA) only in the beginning for representing the map in robot common frame. • In the second stage named "Local Mapping", the robots define which part of their pre-local maps are shared with the other robots based on a Sharing algorithm from the CoMapping framework. Finally the registration process is executed, including an intersecting technique of maps to accelerate the process. We also propose some indicators to analyse the alignment process in both stages. Our proposal has been tested and validated in real situations. The results include maps developed with data acquired on the surroundings of the ECN ( École Centrale Nantes) campus for the urban environment and a farm for the rural experiment, corresponding to the ALFS project.

II. RELATED WORKS

An important topic of the mapping problem is the type of representation used for the environment. Generally 3 types of representations are encountered: feature, grid and topological maps [START_REF] Montemerlo | FastSLAM 2.0: An improved particle filtering algorithm for simultaneous localization and mapping that provably converges[END_REF]. Our system uses two types of map representation: 3D point cloud (feature) for the pre-processing in Pre-Local Mapping and Octree format (grid) for the exchange of maps in the Local Mapping stage later.

In a situation of collaborative mapping, the registration method from a single robot is really important. Several registration applications use Laser Range-finder sensor to build 3D maps [6] [7]. In that case, a high laser scan rate compared to its tracking can be harmful for this task, since it may lead to distortion in the map construction, in which methods based on Iterative Closest Point (ICP) [START_REF] Pomerleau | Comparing ICP variants on real-world data sets -open-source library and experimental protocol[END_REF] can be used to match laser returns for different scans. Implementations, with 2-axis and 3-axis lidar and matches geometric structures of a set of local point generated to finally get a point cloud, were presented in [START_REF] Zlot | Efficient Large-Scale 3D Mobile Mapping and Surface Reconstruction of an Underground Mine[END_REF]. In all those cases, the proposed methods used batch processing to build the maps with accuracy, hence are not applicable to realtime map construction. Regarding our Pre-Local Mapping stage, we reconstruct the map of the environment as a point cloud in real-time using a 3-axis lidar by extraction and matching of geometric features in Cartesian space based on a modified version of the registration method from [START_REF] Zhang | LOAM: Lidar odometry and mapping in realtime[END_REF]. Then our system uses GPS position data to re-localize the cloud in a global frame.

Several methods have been proposed to merge maps. In [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF] and [START_REF] Contreras | Efficient decentralized collaborative mapping for outdoor environments[END_REF] techniques were proposed for 3D merging of occupancy grid maps based on octrees [START_REF] Hornung | OctoMap: An Efficient Probabilistic 3D Mapping Framework Based on Octrees[END_REF]. Their merging process refined the transformation estimate between maps by ICP registration [START_REF] Pomerleau | Comparing ICP variants on real-world data sets -open-source library and experimental protocol[END_REF]. Specifically, in [START_REF] Contreras | Efficient decentralized collaborative mapping for outdoor environments[END_REF] an ICP version was performed including an efficient technique to exchange maps between robots in order to optimize the bandwidth resources of a multi-robot network for decentralized cases. Those last cases are experimentally studied in this paper. To evaluate the alignment post-ICP, we have first to determine the transformation matrices between frames (see Figure 2). Let us say the correct transformation matrix T1 exact (provided by Ground-Truth GPS data for instance) is obtained by a matrix product, which is mathematically denoted by [START_REF] Dinnissen | Map merging of multi-robot slam using reinforcement learning[END_REF].

A. Merging Indicators

T metrics =⇒ T ICP exact = (T 2 ICP ) -1 • (T 1 2 ) -1 • T 1 exact , (1
) where T2 ICP is the ICP alignment transformation matrix and T 1 2 is the initial transformation estimate (provided by GAA low lost GPS data for instance). T ICP exact will represent the transformation matrix used as metric to evaluate the refinement. That matrix T ICP exact will be denoted T metrics . Very often, alignments are evaluated using Euler representation of 6 parameters. The matrix elements ε x , ε y and ε z could be used for the translation evaluation in the x, y, z axis; and ε roll , ε pitch and ε yaw for the orientation. However in order to reduce the complexity of alignment analyse, we propose two indicators: ε * t and ε * r . We can can have a better idea of the orientation evaluation by re-expressing T metrics in terms of an Axis-Angle representation (θ, u) as is shown in (2). In the same way, translation evaluation can be reduced to analyse only the module of the vector with components ε x , ε y and ε z (see ( 3)).

T metrics =    R 3×3 ε x ε y ε z 000 1    , R 3×3 → (θ, u) (2)
ε * t = (ε 2 x + ε 2 y + ε 2 z ) , ε * r = θ, (3)

TABLE I INDICATORS: ICP ALIGNMENT EVALUATION IN SITUATIONS OF KNOWN MAPS TRANSFORMATION

Table I shows the ICP alignment evaluation considering a known initial transformation. We also made the comparison between metrics before and after ICP. Analogously, in order to obtain indicators for correct ICP refinement, it was determined that values of ε * r 0.25 and ε * t 6.14 are considered as coherent merging results. In the real experiments where the exact correct transformation is not known, we can use the ICP matrix as a criterion of the error evaluation, in other words: T metrics = T 2 ICP . In the following sections, our metrics from (3) will be renamed as εr * and εt * , associated to experiments with unknown maps transformations. In the Pre-Local Mapping stage, inspired by the work in [START_REF] Contreras | Efficient decentralized collaborative mapping for outdoor environments[END_REF], its versatile configuration does not depend on a specific LidarSLAM method. A modified version of 1 [START_REF] Zhang | LOAM: Lidar odometry and mapping in realtime[END_REF] was used as LidarSLAM method for this article, considering its good place in the KITTI ranking 2 . Each mobile platform runs this Pre-Local Mapping node. The behaviour of this node is illustrated in Figure 3, where P represents the raw point cloud obtained by a laser scan initially. The accumulated cloud during each sweep k generates P k , which is processed by the Lidar-Odometry algorithm at a frequency around 10Hz. The algorithm receives P k and computes the lidar motion between two successive sweeps, obtaining the transformation T k . Deformation in P k is corrected using the estimated lidar motion in order to use it at a frequency of 1Hz by the Lidar-Mapping algorithm. This algorithm executes the matching and registration of the non-distorted cloud onto a map. Both algorithms (Lidar-Odometry and Lidar-Mapping) are solved with a non-linear optimization, specifically the Levenberg-Marquardt method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Finally, using the GPS information for identifying the initial vehicle pose, it is possible to coarsely project the map of each robot into a common coordinate frame. This projected cloud is denoted as the Pre-Local map. This Pre-Local mapping node can work with different kind of GPS data, either a Ground-Truth option (DGPS-RTK) or an approach most economical, for example using GPS-GGA type information (Global Positioning System Fix Data).

III. METHODOLOGY A. Pre-Local Mapping Stage

B. Local Mapping Stage

The architecture of Local Mapping Stage is described in Figure 4, where the process run on a robot referenced as "i", which shared its map with a robot "n". compression step is executed in each vehicle, using octree structure from the octomap package [START_REF] Hornung | OctoMap: An Efficient Probabilistic 3D Mapping Framework Based on Octrees[END_REF] as a format of compression, optimizing the car-communication.

An overview of this sharing technique is showed in Figure 5, in which the point clouds A and B are the Pre-Local Maps from two different cars "i" and "n" respectively. Both cars receiving simple data about the 3D limits of the maps (i.e. bounding cubic lattice of the point clouds) and the sharing node determines what part of its map will be sent to the other car. These limits are expressed as two vectors: in the first one Amin, the components represent the lowest displacement from the origin along each axis in the point cloud; and in the second one, Amax, is related to the point of highest displacement. The pseudo-code of the map sharing step is depicted in Algorithm 1. First of all, the functions as GetV alues() sort in ascending order the array of components along each axis of the vectors Amin, Amax, Bmin, Bmax and returns the 2nd and 3rd values from this sorted array, named (V 2) and (V 3) respectively. Then for each axis, the mean of (V 2) and (V 3) establishes the Cartesian coordinates (C x ,C y ,C z ) of the geometric center of the exchanged region (S). This region S is a cube whose edge length is 2L and the points from A contained in this region are extracted to generate a new point cloud A sel . At each iteration, the cube region is adjusted until the number of points from A sel is less than a defined threshold (maximum number of points desired to exchange N p max ). Once the condition is accomplished the transfer begins. All the sharing process is analogous on the other robot "n". Finally, encoding of A sel and B sel in octree format is done to reduce the usage of bandwidth resources of the network.

2) Registration step with ICP: The intersecting volumes of A sel and B sel are calculated and renamed as A int and B int , which are down-sampled in order to reduce the time execution of the registration. In the feature descriptors estimation step, the surface normals and curvature of these input clouds are computed in order to improve the feature points matching, which is the most costly step of the ICP [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF] [START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF]. Normal-pointclouds A intN and B intN generated after descriptors estimation are aligned with the ICP algorithm. The ICP refines an coarse alignment between clouds, estimating the best transformation to align a source cloud B intN to a target cloud A intN by iterative minimization of a cost function. Corresponding pairs (b', a') from A intN and B intN are determined iteratively. Least squares registration is executed and the mean squared distance E is minimized with respect to translation t and rotation R:

E(R, t) = 1 N p b' N p b' i=1 a' i -(R b' i + t) 2 , (4)
where A. Urban scenarios: ECN campus case 1) For one robot: First of all, we analysed the impact of GPS quality on the 3D map generation of each vehicle.

N
For that, we took as reference some tests executed with the FLUENCE car in Table I and we performed randomly a path with the robot around the surroundings of the ECN campus. The generated map was twice re-projected on a global frame using two types of GPS data: a high-accuracy DGPS-RTK and standard low cost GPS-GGA. The map with the 1st type of GPS was assumed as referential map, using this DGPS-RTK data as correct transformation matrix T 1 exact . On the other hand, the map re-located with GPS-GGA information was considered as a map with a coarsely placement, used to fill the initial transformation T 1 2 . Results of ICP alignment between those 2 maps were addressed and substituted in [START_REF] Dinnissen | Map merging of multi-robot slam using reinforcement learning[END_REF], to obtain our metrics proposed in (3)

TABLE II ALIGNEMNT EVALUATION: OUR INDICATORS VS BENCHMARK[2] Robot ε * t ε * r ε * t [2] ε * r [2] FLUENCE
1.781047 0.003453 3.005938 0.013086

The first two values in Table II corresponds to translation and orientation evaluation. Regarding the analysis performed in Table I, we have coherent alignment results using GGA vs Ground-Truth data because the values comply with the conditions ε * t 6.14 m and ε * r 0.25 rad. For the moment, we only show our indicator results with indicators proposed by [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF], corresponding to the 3rd and 4th values in Table II. According to [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF], coherent alignments must be ε * t 0.9 and ε * r 0.21, so their metrics showed only correct alignments in rotation.

2) For multi-robot team: In that case, our proposed system was validated considering two vehicles: FLUENCE and ZOE car. For both robots, ENU (East-North-Up) coordinate system was considered as external reference of the world frame {W }, where y-axis and x-axis correspond to the North and East respectively, but coinciding its origin with the GPS coordinates [Longitude: -1.547963; Latitude: 47.250229].

Experiments were realized in urban outdoor environment for an area of approximately 1000m x 700m. For the validation, the robots build clouds from different paths (see Figure 7), running pre-local mapping process in real-time. Results of the Pre-Local Mapping of this experiment are shown in Figure 8, which also illustrates the "sharing region" determined during the map exchange process in each robot. Pre-Local maps from the FLUENCE and ZOE cars were projected on the global frame using Ground-Truth and GGA GPS data respectively. In order to simulate constraints on the network bandwidth and memory usage of the robot-team, the maximum number of points exchanged between cars was set to 410000. In the decentralized case, a meeting point for the robots was defined so they can transfer their maps and perform individually the relative registration process assuming its Pre-Local map as target cloud for alignment reference. Each mobile unit runs an intersecting algorithm, then an ICP refinement obtains an improved transform between maps. Figure 9 depicts the intersection between the shared point clouds during the ICP-alignment process for each robot. We can also see either in Table III or Table IV

B. Rural scenarios: farm case

For that case, our proposed system was validated using the GOLFCAR robot, but in order to do a collaborative mapping approach, the robot path was split in two for simulating two mobile units (see Figure 11(a)). All data come from the rural outdoor environment in an area of approximately 300m x 150m. The external reference of that world frame was also parallel to ENU coordinate system, with a GPS origin given [Longitude: -1.355357; Latitude: 46.809106].

Initial Pre-Local maps and "sharing region" for each simulated robot are exposed in Figure 11(b). Those maps were projected on the global frame using GGA GPS data in both cases. For this scenario, the maximum number of points was set to 50000. Similar to the urban case, each platform mobile executes the intersecting algorithm and then an ICP alignment (see Figure 12) representation. Table V exposes also ICP alignments results with the respective indicators. Again indicators from [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF] are not applicable on this case; on the other hand our proposed indicators confirm coherent merging results since values accomplish the conditions for translation and rotation (ε * t 6.14 , ε * r 0.25). Finally, Figure 13 shows the final 3D map reconstruction for this scenario. Tests demonstrated the influence of working with intersecting techniques, because it allows accelerating the registration by reducing the number of points for processing. In this context, experiments also showed the relevance of the map sharing algorithm by optimizing the performance of the robot-team by reducing the amount of data transferred on the network. Finally, the CoMapping framework remains an appropriate candidate to exchange and generate efficiently large maps with an approach multi-robot suitable for different kind of environments.

V. CONCLUSION AND FUTURE WORK

A cooperative framework for large scale maps applied to urban and rural environments was presented. In this framework, each robot generates its Pre-Local map using 3D-Lidar range measurements and executes individually the merging process considering initially a coarse map alignment with an optimum data exchange. The experimental results highlight the efficiency and versatility of the framework for cooperative mapping with three vehicles in different scenarios. Indicators were also proposed to demonstrate the success of the map merging process for the group of robots. Future work will be addressed to analyse the Sharing step and propose metrics to evaluate the map compression.

Fig. 1 .

 1 Fig. 1. Scheme of the CoMapping System for different cases

Fig. 2 .

 2 Fig. 2. Transformation Matrices Relations

Fig. 3 .

 3 Fig. 3. Architecture of Pre-Local Mapping Stage

Fig. 4 .Fig. 5 .

 45 Fig. 4. Architecture of Local Mapping Stage for one robot "i", receiving map data from another robot "n".1) Map SharingStep: Here lies the key of the efficiency of the collaborative framework[START_REF] Contreras | Efficient decentralized collaborative mapping for outdoor environments[END_REF]. Once the Pre-Local Maps are generated, the robots have to exchange their pointclouds to start the map merging task. However in many cases, the sharing and processing of big map data can penalize the team performance with respect to memory usage and execution time. The sharing technique overcomes this problem, since each mobile platform just transfers a certain volume of its Pre-Local map to the other mobile units. Before the transfer of the map Previous the transference, a map

Algorithm 1 :

 1 Selection of point cloud to share with another robot.

  p b is the number of points b'. The resultant rotation matrix R and translation vector t are applied to source cloud B intN . Again the matching between points from A intN and B intN is re-computed, until the variation of the mean square error between two consecutive iterations is less than a threshold. The final ICP refinement for n iterations is obtained by multiplying the individual transformations: T ICP = n j=1 T j . T ICP is applied to B sel to align and merge with the original cloud A, generating finally the Local Map A L . Similarly, this merging process is executed individually in the other mobile unit. IV. RESULTS Our experiments considered 3 vehicles: a Renault ZOE, a Renault FLUENCE and a GOLFCAR equipped with a Velodyne VLP-16 3D LiDAR of 360 • horizontal and 30 • vertical field of view.

Fig. 6 .

 6 Fig. 6. ZOE (left), FLUENCE (center) and GOLFCAR(right).

Fig. 7 .

 7 Fig. 7. Paths followed by the FLUENCE (green one) and the ZOE (red one) during experiments. Image source: Google Earth.

Fig. 8 .

 8 Fig. 8. Top view of unaligned Pre-Local Maps generated by FLUENCE (green one), and ZOE (red one) projected on a common coordinates system.

Fig. 9 .

 9 Fig. 9. Urban case: Alignment of the intersecting regions with ICP refinement performed in FLUENCE vehicle, when it received the ZOE map (a) Green and red maps represent the target and source clouds pre ICP, top view (b) Green and blue maps represent the target and aligned source clouds post ICP, top view.

  that our proposed indicators confirm coherent merging results because the values are in the range ε * r 0.25 and ε * t 6.14 established from Table I. On the other hand, indicators from [2] do not work for this application because their results are not even in the ranges of translation (ε * t 0.9) or rotation (ε * r 0.21) established by themselves.

Fig. 10 .

 10 Fig. 10. Top view of final 3D Local Map of ZOE car. Finally, Figure 10 depicts one of the final merging results, specifically the final 3D Local map from the ZOE projected on a 2D map in order to make qualitative comparisons.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. (a) Paths followed by GOLFCAR1 (green one), GOLFCAR2 (red one) during experiments. Image source: Google Earth. (b) Top view of unaligned Pre-Local Maps generated by GOLFCAR1 (green one), and GOLFCAR2 robot (red one) projected on a common coordinate system.

Fig. 13 .

 13 Fig. 13. Top view of final 3D Local Map for the rural environment

TABLE III RELATIVE

 III ICP TRANSFORMATION FOR FLUENCE CAR Quantitative refinement results are shown in Tables III and IV. ICP transformations were converted in an Euler parametrization (x, y, z, roll, pitch, yaw) in meters and radians. TableIIIcorresponds to the refinement process for the FLUENCE car, when it received the map from ZOE, and that map is aligned to the own pre-local map from FLUENCE.

	Cars	x	y	z	roll	pitch	yaw
	F-Z	0.62208 εr * = 0.009814 and εt * = 5.186073 (Our indicators) -1.07885 5.03433 -0.00907 0.00294 εr * = 3.028573 and εt * = 6.735261 (Jessup's indicators [2])	-0.00232

TABLE IV RELATIVE

 IV ICP TRANSFORMATION FOR ZOE CAR = 0.007873 and εt * = 6.079665 (Our indicators) εr * = 3.023563 and εt * = 6.814636 (Jessup's indicators [2])

	Car	x	y	z	roll	pitch	yaw
	Z-F	0.27346 εr	0.48724	-6.05394 0.00346	0.00691	-0.00145

*

  = 0.046647 and εt * = 0.640166 (Our indicators) εr * = 3.119865 and εt * = 1.014078 (Jessup's indicators [2])Regarding the ICP transformation, quantitative alignment results relatives to GOLFCAR1(G1) are shown in TableV. All the ICP transformations are also expressed in Euler

				TABLE V
		RELATIVE ICP TRANSFORMATION FOR GOLFCAR(1)
	Cars	x	y	z	roll	pitch	yaw
	G1-G2	0.1461 εr	0.3570	0.5108 -0.01526	-0.04404 0.00138

*

LOAM: https://github.com/laboshinl/loam_velodyne

KITTI ranking: http://www.cvlibs.net/datasets/kitti/ eval_odometry.php

ACKNOWLEDGMENT

The authors gratefully acknowledge support of this work by Lionel Génevé and all the members of the ECN and LS2N institutions.


Rim Christophe Proust
email: christophe.proust@ineris.fr

Ben Moussa

Mohamed Guessasma

Saleh Khashayar

Jérôme Fortin

Proust Christophe

Ben Moussa Rim

Guessasma Mohamed

Fortin Jérôme

Potential accelerating effect of thermal radiation in dust flame propagation : some experimental evidence

Introduction

The role of thermal radiation in premixed flame propagation has been a matter of debate for decades. And it is not only a challenging scientific point, it has significant practical implications. For instance, a route to explain the Buncefield explosion ( [4]) was the implication of tiny particles raised by the blast and promoting flame acceleration through enhanced heat exchanges by thermal radiation in the flame front. In dust explosion protection, the flame is implicitly supposed to propagate like a in a gaseous mixtures but if thermal radiation is dominant for some dusts, many aspects concerning the way to mitigate the explosions for those particular dusts would need to be revised ( [START_REF] Proust | Amplification des effets des explosions sous l'effet du rayonnement thermique . XIVeme congrès of SFGP[END_REF]).

A review of the various aspect linked to this question was proposed recently ( [START_REF] Moussa | Physical mechanisms involved into the flame propagation process through aluminium dust clouds : a review[END_REF]) suggesting in particular that aluminium dust flames could be subject to heat radiation dominated flame propagation regimes. A very similar point was made for the specific case of premixed gases seeded with particles ( [START_REF] Liberman | Effects of thermal radiation heat transfer on flame acceleration and transition to detonation in particle-cloud hydrogen flames[END_REF]). But to date, experimental evidence is severely lacking.

It is intended to provide in this communication, preliminary experimental evidence of the potential promoting role of thermal radiation is dust flame propagation.

Experimental setup

Infrastructure

Experiments were done using the flame propagation tube (figure 1) already used in a number of studies (for instance [START_REF] Proust | Flame propagation and combustion in some dust-air mixtures[END_REF] The experimental chamber is a vertical tube (length 1.5 m and diameter 10 cm), filled with the dust cloud by the bottom thanks to a fluidized bed device. The suspension is rising slowly upwards until complete filling (duration 30 seconds). This way, the suspension may be very homogeneous. The cloud is ignited (spark) at the bottom after having stopped the flow, shut the upper part of the tube with a gate valve and removed the suspension generator. The flame is then propagating upwards freely without being thrusted by the expansion of the burnt products which are vented out through the open end. The dust concentration is controlled by weighting the dust suspension generator before and after the tests, and measuring the volume of the gas flow.

Measurements

High speed video is used to obtain an accurate measurement of the flame speed. Besides a set of transducers is used and installed at about 1 m from the bottom end of the tube (figure 1 right) :

• A fast response fluxmeter. This device (Captec technology) is a thin polymer sheet on each side of which a collection of thermocouples are inserted. On side of the sheet is glued on a cold wall acting as a reference. The voltage between both sides of the sheet is proportional to the heat received thermal flux. In this particular version of the device, optically shielded and bare strips alternate is such a way that only the radiative flux is measured. The response time is lower than 0.1 s. The wavelength absorption range is 0.1 to 12 µm i.e. over nearly all the absorption spectra of CO 2 and H 2 O. The absorptivity is 0.5. Note the fluxmeter is "looking" in the direction of the flame while the photodiode is transversal to the direction of the propagation.

Fluxmeter Ionisation gage

Collimating tube of the photodiode

Pitot system

• A pitot sensor pointing towards the flame front and able to measure the flow velocity ahead of the flame. This technique coupled with the flame velocity measurement (video) was used to obtain the burning velocity in situations where the "tube method" was not applicable as explained hereafter.

Burning velocity

As shown by Andrews and Bradley ( [START_REF] Andrews | Determination of laminar burning velocities : a critical review[END_REF]), a number of methods exist to try and measure the laminar burning velocity of premixed flames. In tubes open at the ignition and closed at the other end, the mixture ahead of the flame is supposed to be at rest so that the burning velocity can be deduced applying the mass conservation on a control volume limited by a cross section upstream, the flame area downstream and the joining portion of the tube (Figure 2). This technique is in practice accurate only if the flame area is easy to determine which is not often the case is rather large tube. In this case, an alternative consists in measuring at the same point of the flame, preferably on the top, the flow velocity and the flame velocity.

The difference between both would give the burning velocity provided the angle of the flame relative to the direction of the propagation is known. A comparison between both methods is shown on figure 3 in the case of premixed methane-air mixtures. A good agreement is found.

2.4

Combustible mixtures A large number of experiments were performed using methane-air mixtures with and without particles in suspension (Alumina and SiC). Some tests were done with aluminium dust air mixtures. The characteristics of the dusts are given in table 1.

Sf Su At Sf Af Su ⋅ = ⋅ Sf-U Su α ( ) ( ) α cos ⋅ - = U Sf Su

Results

Heat fluxes measurements were presented earlier ( [START_REF] Proust | Thermal radiation in dust flames[END_REF]). The measured values are recalled here in view of comparison with the flame dynamics. The reference situation is that of methane air mixtures and the results are presented graphically on figure 3. Clearly, the amount of energy lost by radiation is not large but may play some role in the near limit situation since it is believed that flame may extinguish as soon as the amount of energy lost by the flame amounts about 10%. When alumina particles are added, a continuous reduction of the burning (and flame) velocity is observed. But at the same time, the radiated heat flux diminishes which might simply imply that the flame temperature diminishes. The reduction of the temperature should be linked with the absorption of heat by the particles (conductive exchanges). Extinction is observed when the radiated heat flux amounts that of the limit CH 4 -air flame. Surprisingly, alumina particles do not seem to increase the radiated heat flux. A closer look to the optical properties of this material, reveals that, in the spectral range of methane-air flame (a few microns), the absorption/emission coefficient is as low as 0.1 to 0.2 which may explain this unexpected result. SiC particles have a much higher absorption/emission coefficient (0.7) and a different behaviour is observed (figure 5). And a different result is obtained. It seems that for small particle concentrations (below 150 g/m3) the laminar burning velocity increases so as the radiated heat flux showing a link between both parameters. For higher particle concentrations, the amount of heat absorbed (by conduction) by the solid increases which leads to a reduction of the flame temperature, burning velocity and flux. There seems to be a possibility for a promoting effect of the radiated heat flux even with rather low additional flux, since a 10 to 20% increase of the laminar burning velocity seems to be associated to a 10% increase of the radiated heat flux. It is however possible that the radiated heat flux could be underestimated because of the deposits of powder on the sensor. Note that the temperature of the methane-air flame is rather low and much more effect would be expected with a large burning temperature like in aluminium dust air flames. Alumina and aluminium are poor emitters as compared with SiC but it can be expected that the much larger combustion temperature in the flame (3300°C as compared to 1900°C) would more than compensate for the lower emissivity. Tests were performed with increasing dust concentration (figure 6). Remember that the tube is open at the ignition end (bottom). Clearly when the concentration of dust is rather low, perhaps below 500 g/m 3 , the flame speed reaches a constant value, 0.4 m/s at 80 g/m 3 , 0.7 m/s at 350 g/m 3 , suggesting a burning velocity on the order of 0.3 m/s (in line with available data). But above a certain concentration threshold, a tremendous flame acceleration is observed : the flame speed increases from a few m/s in the ignition zone up to 50 m/s after one meter of propagation. Although the present authors already witnessed abnormally violent explosions in this same configuration, with aluminium dust, it is the first time they managed to record the phenomenon.

Conclusion and perspectives

In this paper, some evidence of the promoting role of thermal radiation on the flame propagation process in dust cloud is provided. It seems to confirm some theoretical predictions ( [START_REF] Liberman | Effects of thermal radiation heat transfer on flame acceleration and transition to detonation in particle-cloud hydrogen flames[END_REF]). One potential consequence of the interaction of thermal radiation in the propagation process would be fast flame self acceleration ( [START_REF] Moussa | Physical mechanisms involved into the flame propagation process through aluminium dust clouds : a review[END_REF]). The present data seem also to confirm this view. Furthermore, the manifestation of the self acceleration seems to depend on the amount of dust in the cloud. As far as aluminium particles are concerned, not less than 500 g/m 3 which is well above the stiochiometric conditions. This trend was suggested earlier in theoretical studies ( [START_REF] Deshaies | Radiative transfer as a propagation mechanism for rich flames of reactive suspensions[END_REF]). However, this work needs to be pursued first to confirm this self acceleration process and second to investigate further the condition of appearance (size of the device, dust concentration, particle size effect,…). In parallel a theoretical study is in progress to try and explain these findings.

Figure 1 .

 1 Figure 1. Experimental setup

Figure 2 :

 2 Figure 2 : laminar burning velocity methods

Figure 3 .

 3 Figure 3. radiated heat flux as function the total heat flux released (methane air flame) 1

1Figure 4 .

 4 Figure 4. methane air flame (stoichiometry) seed with alumina particles

Figure 5 .

 5 Figure 5. methane air flame (stoichiometry) seed with SiC particles

Figure 6 :

 6 Figure 6 : self acceleration of flames propagating in a aluminium dust-air cloud

Table 1 :

 1 main characteristics of the dusts.

	Nature	Specific mass (kg/m 3 ) Mass median diameter (µm) Sauter mean diameter (µm)
	Aluminium oxide 3950	80	40
	Silicium carbide 3210	30	15
	Aluminium	2700	20	10


Yann Grégoire
email: yann.gregoire@ineris.fr

Christophe Proust

Emmanuel Leprette

Didier Jamois

Monitoring of a dust explosion in a 10 m 3 vessel

Introduction

To mitigate the effects of an industrial dust explosion, various explosion mitigation systems may be used (vent panels, explosion suppressing systems, …) [START_REF] Eckhoff | Dust Explosions in the Process Industries[END_REF]. Those need to be certified as dust explosion protection devices according to the ATEX directive and the certification is based on the experimental demonstration of their efficiency. However, such demonstration has to be performed in specific condition that must be representative of their intended use. Today, no standardized procedure exists for generating an explosive dust could in explosion volumes larger than 1 m 3 . Furthermore, large variations due to a scale effect have already been identified between the standardized test in the 20 l Siwek sphere or the ISO6184 1 m 3 vessel [START_REF] Proust | Measuring the violence of dust explosions with the ''20 l sphere'' and with the standard ''ISO 1m 3 vessel'' Systematic comparison and analysis of the discrepancies[END_REF]. However industrial installations are often composed of several enclosures of much wider capacities and the generation of controlled large scale dust explosions is a complicated problem. In this paper we aim to present some of the latest works performed at INERIS in view of generating controlled dust explosions in a 10 m 3 volume and getting a better understanding on the dust flame behavior in such conditions.

Experimental setup

To ease the understanding of the mechanisms at stake, a parallelogram shaped vessel, with a transparent front face, with a cross section of about 1,5 x 1,5 m and a length of 4 m was built (Fig. 1). 2). Details concerning the dust injector are not intended in the present study, but one can summarize its functioning: it is made of a transparent 70 l cylinder pressurized up to 10 bar during the dust injection (after a particle-air mixing phase), which is performed within about 10 seconds through a DN45 pipe. The dust is injected at center of the right wall of the confinement. Most of the tests consisted in monitoring the dust injection without ignition of the cloud to measure the dust cloud repartition, concentration and turbulence. Controlling those three parameters is a fundamental for the generation of reproducible explosions. To achieve this goal, measurements are performed at two scales: a local scale, through the "classical" use of sensors placed at various locations of interest, and a global scale where innovative image processing algorithms are used on the complete video sequences of the dust injection. The local sensors consist in three turbulence gauges and three opacimeters placed at about mid-range in height between the central axis of the vessel and its roof, at the both ends (right and left, 40 cm from the walls) and the middle of the structure. Also a pressure gauge is placed in the dust injector. Alternatively, the global scale measurement consists in optical flow computation and analysis on the cinematographic records of the events (Photron high speed camera at 500 images/s in a 1024 x 464 pixels² window). Explosion tests were also performed (see Fig. 1). In such case, the wall at the left side of the confinement (Fig. 2) was replaced by an inertialess plastic foil (opening when the explosion overpressure exceeds 50 mbar). The ignition of the reactive mixture is performed with 2 5 kJ chemical igniters (similarly as the standard 1 m 3 vessel case).

Phenomenology of the generation of the dust cloud

A complete description of the phenomenology is out of the scope of the present abstract. The present discussion is limited to specific features of the cloud development that are helpful for the understanding of the following parts of this study. The dust dispersion can be summarized with the following schematic: The cloud formation in the confinement can be summarized as follows: first the combustible mixture is injected at a high velocity (≈ 240 m/s). A conical turbulent jet is formed, and impacts the left wall. The mixture moves then towards the right side at a lower entrainment velocity close to 10 m/s. Assuming the highly turbulent jet spread with a half angle of 11.8° [START_REF] Turner | Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows[END_REF], one can evaluate the delay to cover the whole volume of the vessel: 1,25 s. As observed in the cinematic records, the whole 10 m 3 volume is covered within 1,3 s. Such eye-measurement is inaccurate but the result remains consistent with the jet theory. In the meantime, turbulent structures of intensity u' are induced by the reverse movement of the jet and its interaction with the walls, they are expected to homogenize the local turbulence and dust concentration fields. To evaluate the turbulence intensity u' and scale Lt, it is possible to use the free jets models from Hinze [START_REF] Hinze | Turbulence[END_REF]. In the present configuration, the maximal u' velocity should be close to 7 m/s (similarly the integral scale of the turbulence Lt is close to 130 mm, see [START_REF] Proust | Differentiated routes for the simulation of the consequences of explosions[END_REF] for further details on the calculation). It is expected to fall down to 1 m/s within 7 s. This model usually provides results in good accordance with experiments, it has also been compared with other experiments performed at INERIS and published in the scientific literature [START_REF] Roux | Adaptation, extension et promotion de la norme NFU54-540 de décembre 1986, rapport final pour le compte du Ministère de l'Industrie-SQUALPI[END_REF]. At last, with the current dimensions of our system, and the injection pressure, it is estimated that the diffusion time over the whole volume, by the turbulence only would be on the order of 10 s which is comparable to the total time of the complete dust injection. Our objective is to obtain a homogeneous cloud in terms of dust concentration but also with a minimal local turbulence intensity to ensure most of the particles remain in suspension (the sedimentation velocity is expected to be close to 1 m/s in the present case). Note that knowing the turbulence parameters u' and Lt, it is possible to estimate the turbulent combustion velocity St of the dust with models available in the literature, such as Gülder's model [START_REF] Gülder | Turbulent premixed flame propagation models for different combustion regimes[END_REF] (details on the application of the model to large dust cloud are available in [START_REF] Proust | Turbulent Flame Propagation in large Dust Clouds, 11th ISHPMIE[END_REF]).

Local turbulence and concentration measurements

Measuring the turbulence in a dust air mixture is subject to serious difficulties as the application of standard laboratory equipment's such as LDA or hot wire anemometry may be challenging. INERIS has been working for years on the development of alternative techniques [START_REF] Ch | Formation, inflammation, combustion des atmospheres explosives (ATEX) et effets associés[END_REF]. A Pitot tube technique based on a very refined concept of Mc Caffrey gauges [START_REF] Mc Caffrey | A robust bidirectional low-velocity probe for flame and fire application[END_REF] was implemented. Extensive testing in a reference jet proved that such transducers are able to detect eddies as small as 2 cm with a peripheral velocity of 0,2 m/s. Note in particular that the devices can do measurements even in very dusty mixtures (up to 500 g/m 3 at least). Tests were performed in several conditions of dust concentration from 0 to 500 g/m 3 . Typical records are compared with the optical flow method in paragraph 5.

To support flame propagation, a sufficient fuel concentration is needed. In the present case the fuel consists in dispersed particles of micrometric dimensions (of grain size distribution comprised typically between 10 and 100 µm, roughly spherical) of organic material. INERIS worked extensively on concentration measurements in dust clouds over the past 10 years [START_REF] Ch | Formation, inflammation, combustion des atmospheres explosives (ATEX) et effets associés[END_REF]. Our sensors are typically made of light source and a receptor separated by a few tenth of mm. Several devices were tested for the light source (laser, LED) as well as for the receptor (photodiodes, solar panels, …) but also the material of the diaphragms of the systems (aluminum, Teflon), chosen to prevent dust accumulation. A major issue is that the high dust concentration levels (>200 g/m 3 ) in the highly turbulent cloud quickly cloaks the system, making the measurement impossible. To overcome the dust accumulation problem, a cell-phone buzzer was added in the device. Results are promising as even in a very dense dust cloud (4 kg/ m 3 ), the receptor blank level (no dust between the light source and the receptor) change by less than 10 % between before and after the experiment. Another issue that we were confronted to is the limitations of the original Beer-Lambert absorbance law which supposes an extinction coefficient that only depends on the wavelength. However, it is very likely that the micrometric particles will diffuse the light differently as a function of their diameter, thus changing the light intensity on the receptor. A theoretical model was proposed by Proust [START_REF] Ch | Formation, inflammation, combustion des atmospheres explosives (ATEX) et effets associés[END_REF] to take this phenomenon into account:

𝜏 = Ln ( 𝐼 0 𝐼 (𝑥) ) = 3 2 × 𝑄 𝑒𝑥𝑡 ×𝜎 𝛿×𝜌 × 𝑥 (1)
where I0 is the reference light intensity, I the light intensity on the receptor when the light ray passes through a dust cloud of thickness x, Qext an extinction coefficient (close to 2 [START_REF] Ch | Formation, inflammation, combustion des atmospheres explosives (ATEX) et effets associés[END_REF]), σ the dust concentration, ρ the thermodynamic density of the particles and δ their diameter. The system was calibrated using laser granulometry, right next to the opacimeter. Note that a quick dimensional analysis on equation ( 1) shows the physics behind the model: the absorbance remains proportional to the effective absorption area. At last some measurements are performed in the 10 m 3 vessel. An example is provided in Fig. 4: Such measurement is still relatively inaccurate but it gives a satisfying order of magnitude of the local dust concentrations. These sensors are currently being miniaturized to limit their effect on the dust flow.

5-Global flow measurements through optical flow computation

The concept of Optical Flow was introduced by a psychologist, Gibson in 1950 [START_REF] Gibson | The Perception of the Visual World[END_REF] as a description of the movement perception. Robust algorithms to measure it, emerged in the beginning 1980's. These algorithms, based on differential methods, proposed by Horn-Schunck [START_REF] Horn | Determinig optical flow[END_REF] and Kanade-Lucas [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] are now at the basis of most of the more recent algorithms. Optical flow is widely used for motion estimation (mostly in robotics and computer vision) and in video compression. To our knowledge it is progressively being transferred to other domains such as fluid mechanics: however, most research articles related to this topic are entirely focused on the optical flow computation algorithm or its performances (robustness to noise, accuracy) rather than the actual analysis of the physics behind the results obtained (see for example [START_REF] Heitz | Spatio-temporal correlation-variational approach for robust optical flow estimation[END_REF]). Several optical flow algorithms exist, the one we selected for the present is an improved Horn-Schunck formulation based on the works of Sun [START_REF] Sun | A quantitative analysis of current practices in optical flow estimation and the principles behind them[END_REF]. The measurement has some limitations such as it is limited in the front plane of the vessel, or it shows poor performance in the presence of noise or occlusions. However, the code computes the velocity field between 2 images over the whole video, and provides, non-intrusively, for each pixel of each of the images, the velocity vector coordinates. Such calculation is quite demanding as we work on more than 3000 pictures of 1024 x 464 pixels and for each image the optimized yet complex algorithm has to solve several linear systems of up to 950000 equations. We believe it is one of the reasons why such algorithms were barely used until now in fields other than computer science.

Direct comparison between the macroscopic velocities (i.e. front velocities of the jet) found by the optical flow method and a rough measurement on the video gave similar results. Because this method was used for the first time, our first attempt was to ensure its results were in good agreement with the course of the experiment. Here again the agreement is satisfying as the same order of magnitude is found for the sensors and the optical flow method (note that the flow measurement is performed in the plane surface located at 75 cm before the turbulence gauge, so higher values of the velocity were expected). At last an interesting feature of the optical flow measurements was exploited. It is the dense nature of the results, which provide a large number of velocities over each image. This was used to compute statistically the standard deviation of the velocity field through multiple convolution products of the velocity matrix with square windows between 7 and 101 pixels² (see Fig. 5, right). The standard deviation of the local velocity field represents the turbulence intensity [START_REF] Hinze | Turbulence[END_REF].

Here the same order of magnitude as estimated in paragraph 3 with the phenomenological analysis is found (7 m/s). Further evolutions of the overall turbulence in the vessel can be tracked with an acceptable accuracy, showing a drastic fall of the turbulent velocities after 5 s. Consequently, this tool seems very promising for the statistical turbulence analysis. Next steps of the study will include tentative of measuring the integral scale of the turbulence through optical flow calculation.

Further developments

Measuring the turbulence intensity and scale, as well as the dust concentration is helpful to determine when and where the cloud can be ignited. In the present example (with 500 g/m 3 of wheat flour), ignition is performed in the jet close to the injection location, 5 s after the beginning of the dust injection. In those conditions the dust concentration and turbulence cannot fully diffuse and homogenize through the whole volume. Also there is a pressure gradient due to the jet, but this is the best compromise so far to prevent the sedimentation of the dust. A possible solution to this problem would thus be to a have two injectors rather than this single one or an additional flow mixing device in the vessel. The transparent face is also of a great help to analyze the flame behavior: a conical flame development is observed, it grows until it reaches the left wall, then goes backwards despite the vent opening, following the dominant flow velocities of the unburnt material. Optical flow analysis will also be performed on these tests, in view of characterizing numerically the flame growth and velocity. Implementing this flame surface growth in INERIS phenomenological tool EFFEX for confined explosion [START_REF] Proust | EFFEX: a tool to model the consequences of dust explosions[END_REF], make it possible to model the explosion with an increased accuracy. At later time explosions performed in this vessel with a transparent face with different kinds of particles will also aim to select the best test parameters (ignition delay, injection pressure) before testing of protective equipment in imposed conditions of dust nature or concentration, KST, ...

Figure 1 .

 1 Figure 1. Explosion in a 10 m 3 vessel with a transparent face

Figure 2 .

 2 Figure 2. 10 m 3 vessel with a transparent face Organic dust (wheat flour, cornstarch, …) are often selected as fuel in INERIS applications. Those are injected into the vessel with a specifically designed injector based on fluidized particle bed technique (visible on the right side in Fig.2). Details concerning the dust injector are not intended in the present study, but one can summarize its functioning: it is made of a transparent 70 l cylinder pressurized up to 10 bar during the dust injection (after a particle-air mixing phase), which is performed within about 10 seconds through a DN45 pipe. The dust is injected at center of the right wall of the confinement. Most of the tests consisted in monitoring the dust injection without ignition of the cloud to measure the dust cloud repartition, concentration and turbulence. Controlling those three parameters is a fundamental for the generation of reproducible explosions. To achieve this goal, measurements are performed at two scales: a local scale, through the "classical" use of sensors placed at various locations of interest, and a global scale where innovative image processing algorithms are used on the complete video sequences of the dust injection. The local sensors consist in three turbulence gauges and three opacimeters placed at about mid-range in height between the central axis of the vessel and its roof, at the both ends (right and left, 40 cm from the walls) and the middle of the structure. Also a pressure gauge is placed in the dust injector. Alternatively, the global scale measurement consists in optical flow computation and analysis on the cinematographic records of the events (Photron high speed camera at 500 images/s in a 1024 x 464 pixels² window). Explosion tests were also performed (see Fig.1). In such case, the wall at the left side of the confinement (Fig.2) was replaced by an inertialess plastic foil (opening when the explosion overpressure exceeds 50 mbar). The ignition of the reactive mixture is performed with 2 5 kJ chemical igniters (similarly as the standard 1 m 3 vessel case).

Figure 3 .

 3 Figure 3. Schematic of the dust cloud formation in the 10 m 3 vessel

Figure 4 .

 4 Figure 4. Concentration measurements at the back and entrance of the 10 m 3 vessel

  Secondly it has been decided to compare the velocities measured by our Pitot turbulence sensors (see paragraph 4) with the velocity field in the close surrounding of the sensors (Fig 5, left):

Figure 5 .

 5 Figure 5. Left: Optical flow calculation of the turbulence intensity in the vicinity of a turbulence gauge compared with the gauge measurements. Right: Turbulence intensity estimated in the whole domain with the optical flow method and variable convolution windows