Hacène Fouchal
email: hacene.fouchal@univ-reims.fr

Alain Ninet
email: alain.ninet@univ-reims.fr

Crestic

Partial Signature for Cooperative Intelligent Transport Systems

Keywords: C-ITS, VANETs, security, authen-ticity

On C-ITS (Cooperative Intelligent Transport Systems) vehicles send and receive sensitive messages informing about events on roads (accidents, traffic jams, etc, ..). The authentication of these messages is highly recommended in order to increase the users confidence about this system. This authentication ensures that only messages coming from trusted vehicles are accepted by receivers. An adapted PKI (Public Key Infrastructure) for C-ITS provides certificates for each vehicle. The certificate will be used to sign messages. This principle is used within deployed C-ITS solutions over the world. This solution is easy to implement but has one major flaw: each message needs to be sent with its signature and its certificate. The size of the message to send becomes high. In the meantime, for many C-ITS use cases, each message is sent many times for robustness reasons. The communication channel could be overloaded. In this paper, we propose to split the signature into some equal parts. When a message has to be sent, it will be sent with one of these parts. A receiver will save the received message with its actual part. For each reception, it will collect the remaining signature parts until all the signature parts are received. Our solution is implemented in a C-ITS architecture working through Bluetooth protocol using the advertising model. The solution is applicable for vehicle speeds reaching 130 km/h. We have proved, through a set of real experimentations, that our solution is possible.

I. Introduction

The deployment of C-ITS is a hot challenge in the era of smart cities. Urban roads and highways are often crowded. The dissemination of relevant information about traffic jams, accidents and road works help users to reduce their travel and to travel in safer conditions. A dedicated WIFI has been designed for connected vehicles: IEEE 802.11p (denoted also ETSI ITS-G5). The deployment of ITS-G5 hotspots (denoted RSU-Road Side Units) is on-going in many cities. V2V (vehicle-to-vehicle) and V2I (vehicle to infrastructure) communications become possible. The penetration rate of such technologies is still very low which makes the coverage of such technology very limited. In a previous work [START_REF] Thomas | Intelligent Transport System Based on Bluetooth[END_REF] we have used the BLE protocol (Bluetooth Low Energy) on smartphones used as receivers on vehicles. In this paper, we use the same environment for our experimentations.

We need to consider the authentication of senders in order to avoid untrusted drivers to notify fake events. It is recommended to handle PKI in charge of distributing certificates to all trusted vehicles in the eco-system. All messages are sent with their signature computed using their certificates. The certificate is also sent within the payload message.

Over ITS-G5, the message size could reach 1500 bytes. Most of usual notification messages (DENM: Decentralized Event Notification Messages) have a payload size between 200 and 300 2 bytes where the size of the signature and the certificates is about 300 bytes.

Many solutions have been proposed to reduce this size; some of them have been deployed in ETSI standards as sending certificates only ondemand, sending the hash of the certificates.

The main contribution of this paper is to propose an additional mechanism: splitting the signature into many equal parts (partial signature). In each message only a partial signature is sent. When a message has to be sent, it will be sent with one of these parts. A receiver will save the received message with its actual part. For each reception, it will collect the remaining signature parts until all signature parts are received. Since messages are sent many times, we only modify the signature part in each sent message.

We have proven in [START_REF] Thomas | Intelligent Transport System Based on Bluetooth[END_REF] that BLE could be used with efficiency with vehicles, our present contribution is to experiment our splitting mechanism within the BLE protocol. Indeed, in the version 4.2 of this protocol (and later versions), the packet size could be 32 bytes for the advertising option. We have already proposed an architecture based on this protocol where the signature size is 16 bytes. But the whole signature size is 64 bytes (using the ECDSA algorithm: Elleptic Curve Digital Signature Algorithm). For our experimentations we split the signature into 4 parts. We have shown that for a speed of 130 km//h the vehicle is able to receive more than 20 messages during the receiving slot time (nearly 6 secondes). It is able to rebuilt the original message with all received messages and to verify the message signature even if some of these parts are dropped. This paper is composed as follows: Section II describes some works dedicated to vehicular networks and security of exchanged messages. In section III we present the unsual security mechanism handled in C-ITS. Section IV presents the experimental environment of the study. V is dedicated to detailed processes used in this solution. Section VI give some conclusions and ideas about future enhancements of the study. [START_REF] Bhover | V2X communication protocol in VANET for co-operative intelligent transportation system[END_REF] presents another alternative to WAVE/DSRC. It mixes three technologies: Wi-Fi Direct, ZigBee and Cellular Network. Wi-Fi Direct is used as a direct link between nodes. ZigBee is used to connect roadside sensors and Cellular Network for long distance communications. [START_REF] Jeong | A Hybrid V2X System for Safety-Critical Applications in VANET[END_REF] uses the same principle, however this paper does not take into account the hybridation between ITS-G5 and Cellular Network.

II. Related works

In [START_REF] Santa | Deployment of vehicular networks in highways using 802.11p and IPv6 technologies[END_REF], the authors propose a new network architecture deployed in Spain on vehicles switching between 802.11p and 3G, depending on the availability of RSUs.

[5] presents a detailed study on performance evaluation of IEEE 80211.p networks versus LTE (Long Term Evolution) vehicular networks. The authors have measured some performance indicators as the end-to-end delay for both networks in different scenarios (high density, urban environments, etc.). Many important issues have been measured as network availability and reliability. The authors have proved through simulations that LTE solution meets most of the application requirements in terms of reliability, scalability, and mobility. However, IEEE 802.11p provides acceptable performance for sparse network topologies with limited mobility support.

[6] studies throughput of VANETs with unidirectional traffic for different conditions and transmission ranges of wireless equipments. All studied vehicles are randomly connected. The paper gives few results of simulation studies achieved on NS-2 toolbox. They have measured performances indicators in case of congestion.

In [START_REF] Lu | Evaluation of vehicular communication networks in a car sharing system[END_REF], an evaluation of vehicular communication networks through car sharing scenarios is detailed. They adopted a specific mobility model which has been imported to a simulator. They have worked on a grid Manhattan network and they have observed some performance parameters such as delay, packet loss, etc. They have shown that vehicular communication are realistic under some conditions. Some recent studies have proposed various solutions either for privacy [START_REF] Kim | UDiPP: A framework for differential privacy preserving movements of unmanned aerial vehicles in smart cities[END_REF], [12], [START_REF] Fouchal | Taladrez A security scheme for wireless sensor networks[END_REF], [START_REF] Sedjelmaci | A survey on pseudonym changing strategies for Vehicular Ad-Hoc Networks[END_REF] 3 and authenticity [START_REF] Zeadally | Vehicular ad hoc networks (vanets): status, results, and challenges[END_REF], [START_REF] Fouchal | Taladrez A security scheme for wireless sensor networks[END_REF].

III. Security architecture A. Preliminaries

In the area of C-ITS (Cooperative Intelligent Transportation System), a protocol stack has been defined and standardized by the ETSI standardisation institute in Europe [START_REF]Conference on Computing, Networking and Communications (ICNC): Wireless Networks[END_REF]. Over the Transport-Networking layer (defined as geonetworking layer), the Facilities layer has been designed in order to be an efficient interface between the application layer (close to the driver and the vehicle sensors) and the Transport-Networking layer. Many types of messages are provided by this layer. In this paper will focus only on one main message, DENM (Decentralized Environmental Notification Message). A vehicle sends DENM messages in order to inform of any event (i.e. accident, traffic jam, etc.). The event could be triggered automatically thanks to the connexion to the vehicle CAN (Controller Area Network) Bus. In fact, some smart rules are run on the OBU (On-Board Unit) which generates automatically appropriate messages. It could also generates manual messages for sensitive cases as animal on the road.

B. Privacy and authenticity

C-ITS security considers two main aspects:

• Authenticity: this aspect allows to consider only messages coming from trusted drivers. • Privacy: this aspect allows to protect drivers data and avoid driver tracking. Each involved vehicle in the eco-system has to subscribe to the PKI server. The subscription allows to provide a long term certificate (LTC) to the vehicle (signed by a root PKI).

A vehicle could sign its messages using this LTC by means of a TPM (Trusted Plateform Module) which embeds cryptographic processes and private keys with safety. Authenticity is very well ensured. But with such a mechanism, privacy is not guaranteed at all. Indeed, from a unique LTC we could extract private informations which give the opportunity to external observers to track drivers. In order to ensure privacy, it is usual to work with a set of pseudonym certificates (PCs) for small periods. A vehicle should have a pool of PCs which should be up to date in order to be able to switch to another PC.

IV. BLE environment

In this section, we first propose the experimentation environment used in this paper. BLE advertising mechanism is used in the experimentations of our solution. The principle is to send beacons (16, 32, 64 or 128 bytes). BLE has been designed in order to advertise data wide areas in order to inform users about any product (in particular in restaurants, shopping centers, museums, ...). The communication is unidirectional without any connection between the sender and receivers. The adverser broadcasts an « Advertising » BLE message composed of Beacon code (which is a value fixed by the sender about the product to advertise) and any other relevant information. The receiver should run BLE scanner able to accept any BLE beacon and is able to wake-up a dedicated application according to the Beacon code value.

The senders denoted BLETransmitter, are deployed on roads and are connected to road operator management servers which monitor them by sending all relevant informations on roads.

A. Adapted notification message

We have defined a message format (inspired from the standard DENM format) containing the minimum information required to notify any event on the road. We only keep the mandatory fields defined in the ETSI standard. Our format is defined as follows:

CC Lat Long H T OS S

Each field is defined as follows: Cause Code CC) or Event Cause Code mixed with its sub cause code (A translation table has been defined to match with the DENM standard codes) ; Latitude Lat / Longitude Long for The GPS Location ; Heading (H) -defines the angle between the vehicle's direction and the North side ; Timestamp (T) -defines the instant where the message has been produced ; Offset (OS) -defines the offset 4 of the signature part ; Signature(S) -defines the signature part of the message.

V. A proposed authentication solution

In this section, we give details of our contribution: splitting a signature into equal parts and insert a partial signature in each message to send.

A. Principle

We consider a message M to send. The signature S is provided by the signature algorithm embedded on the C-ITS station. We split S into N parts such as S = S1, S 2 , ..., S N . Each part S i could be included at anytime in the whole message to send.

B. Communication algorithms

The algorithm 1 describes the actions to be run by each sender. When a message is valid, it has to be sent in a continuous way. But at each time, the signature field changes. One part among all parts is selected randomly and is inserted in this field of message. The sending of messages ends when the message is no more relevant (when it is not valid for instance).

The algorithm 2 describes the actions to be run by each receiver. The receiver catchs the first message and keeps receiving the same message but selects the received part of the signature, it is sorted in a temporary variable. The next time the same message is received, the signature part is also selected. When all the signature parts are received, the verification of the signature is run. If the signature is valid, then the message is accepted. We have experimented the message sending with various vehicle speeds. Figure 1 shows how many times the message is received in the slot window of the sender. In fact this number depends on the protocol range. Indeed, the BLE protocol is supposed to have range of 100-500 meters. When a fixed sender broadcasts its message, it will cover a range of distance D. A vehicle driving through D will receive the message as long as it in D. When the speed increases, the time slot where the vehicle is in D is shorter, then the number of received beacons is lesser. In our experimentation, we need to split the signature into 4 parts in order to fit with the defined structure of the adapted notificationpacket. Then if the vehicle speed is less than 130 km/h, this decomposition is possible and the vehicle is able to rebuilt the entire signature.

VI. Conclusion

In this paper we have proposed a simple solution for reducing the signature size of messages in the area of C-ITS. Indeed, the main idea is to split the signature in small parts. We only insert one part in each sent message. The messages will be repeated many times but with other signature parts.

The receiver will collect all messages and rebuilt the entire signature with all different parts. When all parts are recovered, the message could be verified and could be computed.

We have experimented this mechanism with the BLE protocol. This protocol is suited for such environment: beacons are sent continuously (the frequency could be 10 hz). We have shown that for the BLE protocol the vehicle speed has an impact on the obtained results. The splitting is different. More the speed is high, less the number of signature parts is.

We intend to investigate if we could verify the validity of the message without getting all the signature parts. Indeed, in this study we assume that all different messages are properly received. As a future works, we intend to test better the scalability of our system by launching simulations. In the meantime we intend to analyse the security performances of signing and verifying processes.

Algorithm 1 :

 1 Sender Algorithm 1: Begin -M : Message -S : Signature parts= S = S1, S 2 , ..., S N -V T : Validity slot time of the event 2: while (V T is valid ) do 3: sends M + Random(S, N ) 4: end while 5: End Algorithm 2 : ReceiverAlgorithm 1: Begin -M i : Message i -RS : ReceivedSignatureSave contains the signature parts 2: while (MessageVerified( M i ) == FASE ) do

Figure 1 :

 1 Figure 1: Number of beacons vs. speed

International Conference on Computing, Networking and Communications (ICNC): Wireless Networks 978-1-7281-4905-9/20/$31.00 ©2020 IEEE

International Conference on Computing, Networking and Communications (ICNC): Wireless Networks

Acknowledgement

This work was made possible by EC Grant No. INEA/CEF/TRAN/A2014/1042281 from the INEA Agency for the SCOOP project. The statements made herein are solely the responsibility of the authors.


Monica Ashokumar

Clément Guichet

Jean-Luc Schwartz

Takayuki Ito

Relationships between the effect of orofacial somatosensory stimulation in speech perception and speech production performance

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

➢ Orofacial somatosensory inputs modify the perception of speech sounds [START_REF] Ito | Somatosensory Function in Speech Perception[END_REF][START_REF] Trudeau-Fisette | Auditory and Somatosensory Interaction in Speech Perception in Children and Adults[END_REF] . This is related to the role of somatosensory system in production, transferred to perception through sensory-motor relationships in the human brain. Thus, somatosensory effect in perception may vary based on production ability.

Methods

➢ We examined whether somatosensory effect in speech perception was correlated with production ability in corresponding vowels.

•

Participants: Nineteen French native speakers.

➢ Speech Production Test: French words, 'Dé' for /e/ and 'Deux' for /ø/ were recorded.

• Production index: Difference in the first, second and third formant frequencies (F1, F2 and F3) between /e/ and /ø/.

Results

➢

2 Aim➢

 2 Somatosensory effect on speech perception • Orofacial somatosensory stimulation significantly increased the amount of /e/ responses (F (1,18) = 7.42, p < 0.05) Relationships between the effect of orofacial somatosensory stimulation in speech perception and speech production performance Monica Ashokumar 1 , Clément Guichet 1 , Jean-Luc Schwartz 1 , Takayuki Ito 1,To investigate the relationships between the somatosensory effect in speech perception and speech production performance. Acknowledgement This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No 860755 Contact: monica.ashokumar@grenoble-inp.fr ➢ Somatosensory effect in speech perception: Vowel identification test with an 8-member /e/-/ø/ continuum was carried out. • Somatosensory stimulation associated with facial skin deformation was applied with the presentation of auditory stimulus. • Perception index: Difference in categorical boundary between the conditions with and without somatosensory stimulation (Skin-stretch and Control).

1➢F3

  GIPSA-lab, CNRS, Univ. Grenoble Alpes, Grenoble INP, Grenoble, France, 2 Haskins Laboratories, New Haven, USA ➢ Correlation analysis Conclusion Results indicated that the participants who have a large difference between /e/ and /ø/ showed large somatosensory effect in speech perception. ➢ Somatosensory effect in speech perception can be ascribed to speech production performance. Experimental setup with somatosensory stimulation /e/ or /ø/? F1-F2 vowel space for /e/ and /ø/ • A reliable correlation with F2 (p<0.05) and marginal correlation with F3 (p=0.083) but no correlation with F1 (difference between /e/-/ø/ (in Hz) r = 0.41, p = 0.083 Production Index Production Index F1 difference between /e/-/ø/ (in Hz) F2 difference between /e/-/ø/ (in Hz) 17, p = 0.49 r = 0.48, p = 0.036


Chaabouni Salma
email: salma.chaabouni@ensat.fr

Pech Latché Alain

Jean Claude

Bouzayen Mondher

Tomato Aux/IAA3 and HOOKLESS are important actors of the interplay between auxin and ethylene during apical hook formation

Keywords: apical hook, hormone signaling, auxin/ethylene cross-talk, tomato, transcription factors

Plants implement differential cell growth as an adaptation process in order to direct their development in a way that allow them to better cope with the environmental conditions. This process requires the complex integration of multiple hormone signalings, though, a lot remain to be known about the mechanisms and the molecular actors that take part in this hormonal dialogue. We have previously shown that Sl-IAA3, an Aux/IAA gene, is a molecular link between auxin and ethylene responses in tomato plants. We show here that the expression of Sl-IAA3 in etiolated seedlings is restricted to the inner side of the apical hook, opposite to that of the HOOKLESS gene whose loss-offunction mutation results in the absence of hook formation.

We propose a model on how auxin and ethylene modulate the expression of Auxin Response Factor 2 (ARF2) via IAA3 and HLS protein to regulate hypocotyl bending.

The coordination of plant developmental processes and growth adaptation rely on complex interplay between individual hormone and non-hormone signalings. The phytohormones auxin and ethylene are essential regulators of plant development and the mechanisms governing the interactions between these two hormones are becoming better understood even though, so far, only few molecular players of this cross-talk have been identified. [START_REF] Alonso | Five components of the ethylene-response pathway identified in a screen for weak ethylene insensitive mutants in Arabidopsis[END_REF][START_REF] Chilley | The POLARIS peptide of Arabidopsis regulates auxin transport and root growth via effects on ethylene signaling[END_REF][START_REF] Swarup | Ethylene upregulates auxin biosynthesis in Arabidopsis seedlings to enhance inhibition of root cell elongation[END_REF] The ethylene-mediated regulation of auxin biosynthesis occurs through the activation of WEI2/ASA1 and WEI7/ASB1, anthranilate synthase subunits that catalyse the first step in tryptophane biosynthesis. The reciprocal effect of auxin on ethylene biosynthesis through the activation of several ACC synthases has also been described. [START_REF] Stepanova | A link between ethylene and auxin uncovered by the characterization of two root-specific ethylene-insensitive mutants in Arabidopsis[END_REF] A comprehensive study combining physiological, genetic and genomic approaches uncovered a simple mechanistic model for the interaction between the two hormones in roots. Indeed, in addition of acting independently on the same target genes, ethylene and auxin can reciprocally regulate each other's biosynthesis and influence each other's response pathways. [START_REF] Stepanova | Multilevel interactions between ethylene and auxin in Arabidopsis roots[END_REF] This model provides a likely explanation for the strong ethylene response defects observed in auxin mutants. In accordance with these data we have recently shown that phenotypic responses to the downregulation of Sl-IAA3 gene in tomato include alterations to the classical auxin-regulated processes of apical dominance and hypocotyl elongation as well as to typical ethylene responses such as apical hook formation in etiolated seedlings and leaf epinasty in light-grown plants. [START_REF] Chaabouni | Sl-IAA3, a tomato Aux/ IAA at the crossroads of auxin and ethylene signalling involved in differential growth[END_REF] These results suggest that Sl-IAA3 is an integral regulator of auxin and ethylene responses in tomato plants.

The induction of apical hook formation in Arabidopsis represents one of the best described examples of auxin-ethylene cross-talk in plants and the hook is a result of differential cell elongation on opposite sides of the hypocotyls. [START_REF] Lehman | HOOKLESS1, an ethylene response gene, is required for differential cell elongation in the Arabidopsis hypocotyl[END_REF][START_REF] Raz | Regulation of differential growth in the apical hook of Arabidopsis[END_REF][START_REF] Li | Convergence of signaling pathways in the control of differential cell growth in Arabidopsis[END_REF] An Arabidpsis mutant lacking differential growth in the apical region of the hypocothyl (hook-less1) has been identified and it was proposed that the mutated gene encodes a key regulator that integrates ethylene and auxin signaling pathways during apical hook formation. [START_REF] Lehman | HOOKLESS1, an ethylene response gene, is required for differential cell elongation in the Arabidopsis hypocotyl[END_REF] Subsequently, a partial suppressor of hls1 phenotype resulting from a loss-offunction mutation in the ARF2 gene (Auxin Reponse Factor) was isolated. [START_REF] Li | Convergence of signaling pathways in the control of differential cell growth in Arabidopsis[END_REF] Interestingly, our recent data indicate that antisensemediated inhibition of Sl-AA3, a tomato Aux/IAA gene, results in exaggerated hook formation associated with a downregulation of ARF2 expression. By contrast, accumulation of Sl-HLS transcripts is not altered in the AS-IAA3 plants suggesting that the exaggerated hook formation in the transgenic lines does not involve an alteration in Sl-HLS expression.

To further investigate potential interactions between Sl-HLS and Sl-IAA3 in controlling hook formation we analyzed the spatial expression of the Sl-HLS and Sl-IAA3 in the apical hook by native promoter-reporter constructs. We generated transgenic lines expressing a 1.3 kb fragment of the Sl-HLS promoter fused to the GUS reporter gene and assessed the GUS staining in etiolated tomato seedlings treated with ethylene. Remarkably, ProHLS::GUS staining was restricted to the outer side of the hook curvature, whereas the Sl-IAA3 promoter drove GUS staining exclusively on the inner side (Fig. 1). These data suggest that Sl-IAA3 acts as a repressor of auxin/ ethylene-mediated cell elongation on the inner surface of the apical hook and/or conversely that SI-HLS1 is involved in promoting cell elongation on the outer surface. Sl-IAA3 and Sl-HLS genes provide therefore tissue-specific markers for the inner and outer sides of the apical hook, respectively, and the corresponding promoters could be useful to target the ectopic expression of transgenes to a specific side of the hook.

We postulate that Sl-IAA3 and Sl-HLS may act in parallel pathways both of them involving ARF2 as a downstream component, or Sl-HLS may act upstream of Sl-IAA3 to downregulate its expression which might explain why Sl-IAA3 is not expressed in the upper side of the hook where the expression of Sl-HLS is high. Figure 2 depicts a model mechanism on how the convergence of ethylene and auxin signaling impacts differential growth of etiolated seedlings in an IAA3-dependent manner. IAA3 serves as central integrator of ethylene and auxin. In this model, HLS-dependent and IAA3-dependent regulation of ARF2, a negative regulator of the differential auxin response, leads to enhanced differential growth and exaggerated hook curvature. In conclusion, we suggest that the process of hook formation requires an interplay between HLS, IAA3 and ARF2 proteins.

Continued efforts are now engaged to determine the posttranslational regulation of IAA3 by auxin and ethylene and to test the potential interaction between IAA3 and ARF2 proteins.

The Universal Genome Walker Kit (Clontech Laboratories, Inc., Palo Alto, CA, USA) was used to isolate 1.3 kb of the Sl-HLS gene promoter region. The Sl-HLS promoter was then fused to the b-glucuronidase (GUS) reporter gene in the plp100 binary vector [START_REF] Szabados | New plant promoter and enhancer testing vectors[END_REF] and used for stable tomato transformation [cv. MicroTom]. Growth conditions were performed as described previously. [START_REF] Chaabouni | Sl-IAA3, a tomato Aux/ IAA at the crossroads of auxin and ethylene signalling involved in differential growth[END_REF] For histochemical GUS analysis, ProHLS::GUS and ProIAA3::GUS transgenic lines were incubated at 37°C for 5 h with GUS-staining solution (100 mM sodium phosphate buffer, pH 7.2, 10 mM EDTA, 0.1% Triton and 1 mM 5-bromo-4chloro-3-indolyl-b-D-glucuronic acid). Following GUS staining, samples were washed several times to extract chlorophyll using a graded ethanol series.

Figure 1 .

 1 Figure 1. The expression of Sl-IAA3 and Sl-HLS genes takes place on opposite sides of the hook. Tomato ProHLS::GUS and ProIAA3::GUSseedlings were dark-grown for 5 days and then treated for 48 h with 10 μl l -1 of ethylene. The images are representative of at least three independent experiments with n > 30 seedlings per experiment.

Figure 2 .

 2 Figure 2. A model mechanism describing the different players of the auxin/ ethylene interplay leading to differential growth in etiolated seedlings.


A Ponsich

C Azzaro-Pantel
email: catherine.azzaropantel@ensiacet.fr

S Domenech

L Pibouleau

F Pigeonneau

A systemic approach for glass manufacturing process modeling

Keywords: glass manufacturing, reactor networks, mathematical programming optimization

Due to the difficulty in defining accurate models to represent the physical phenomena in glass furnaces, a large number of studies consider the process as a black-box system. The systemic approach proposed in this paper relies on a decomposition of the furnace into two entities, i.e., the combustion chamber and the glass bath. Next, a reactor network structure is automatically designed to reproduce the hydrodynamic features in the considered zone. These computations are based on the minimization of the quadratic difference between the network output and Residence Time Distributions (RTD) obtained through CFD simulations.

The optimization results provide a network structure showing good agreement with the RTDs computed by CFD. In addition, they highlight, in both furnace zones, some useless volumes that do not participate to the global flow. These results stand as a good basis to extend the study to the computations of kinetics, energy or environmental criteria, which in turn may be optimized.

Introduction

During the last decades, glass industry has gained a great experience in the operating mode of glass melting furnaces, thus leading to the implementation of accurate and efficient processes. The study of such furnaces has constituted a major research activity of many companies. However, the extreme temperature conditions existing within the furnaces make it very difficult to observe and understand the occurring physico-chemical phenomena. Consequently, only a few models have been developed to understand the global behavior of the furnace and to interpret the process trends.

The main purpose of this study is to present a functional analysis of some glass manufacturing processes, particularly focusing on a description of the furnaces. In a more general and long-term perspective, the underlying idea consists in proposing generic models that would fit very well with the furnace behavior, in such a way that some performance criteria, defined according to the glass industry needs, could be optimized. These criteria might be oriented, for instance, towards technico-economic, environmental or energy requirements. In this context, the present study only constitutes a feasibility analysis, in order to lay the foundations of a more complete work. The paper is organized as follows: in the second section, a synthesis of the literature devoted to glass processes is proposed. This preliminary analysis leads to the definition of the chosen methodology, i.e. a systemic approach based on reactor networks modeling in Section 3. Section 4 shows the first results of the work, while some conclusions and perspectives are drawn in the last section.

Literature review

This literature review gives a global roadmap of the various works dedicated to the representation of glass manufacturing furnaces and their behavior. It is worth mentioning that even if the number of research works dealing with the physical phenomena occurring in the glass furnace is important, the variety of studies is quite few diversified. Indeed, it appeared that the mathematical modeling of these phenomena is based on two kinds of representation modes. The former, which constitutes the great majority of them, lies on fluid dynamics computations. The second research domain does not aim at accurately describing the involved internal mechanisms, but mostly views the furnaces as 'black-box' systems and only analyses inputs and outputs, leading to either learning / training techniques, or reactor network approaches.

Computational Fluid Dynamics (CFD)

Most published works on glass manufacturing furnaces modeling deal with CFD techniques, developed in order to represent locally and globally the hydrodynamics, heat or mass transfer phenomena taking place in the furnaces. For instance, studies such as in [START_REF] Carvalho | Modelling and optimization of the no formation in an industrial glass furnace[END_REF] try to model and optimize the chemical development of NOx pollutants, emitted in the laboratory smokes, in the combustion chamber. The optimization variables are the combustible composition, the air pre-heating conditions (combustive) and the combustion reaction stoichiometry. In [START_REF] Flesselles | Phénomènes de transferts dans les fours de verrerie industriels[END_REF], convection models are adopted to represent the convection intensity in the glass bath, while the effect of the addition of elements to the furnace structure.

In the CFD simulation framework, the furnace geometry, the melted glass physical properties, the influence of spatial temperature distribution and many other factors constitute parameters. However, the use of CFD techniques to represent the involved physical and chemical phenomena with both a good accuracy and acceptable computational times, is currently a challenging issue. Thus, the aim of our study is to propose an alternative computational strategy that would avoid such CFD calculations.

In spite of the improvements both in discretization methods (finite volume/finite elements) and computational powers, the main drawback of the CFD approach still lies in the computational times, which may be, in some cases, prohibitive. In that context, it is difficult to consider some phenomena such as, for instance, the kinetics of the chemicals involved in the process. The models then have to be simplified in order to integrate them in the CFD computations.

Learning methods

In a great variety of research areas, several methods were developed to improve the understanding of systems viewed as 'black-box' systems. The learning methods are based on the availability of data measured on real systems and then use inference mechanisms to build rules extrapolating the observed trends. In the literature devoted to the study of glass manufacturing systems, the learning methods result mainly from the Artificial Intelligence area.

• Artificial Neural Networks (ANNs) consist in generating analytical expressions of the system outputs according to the input data. These analytical expressions may be linear or not, and lead to the formulation of predictive models. In [START_REF] Carnahan | Predicting glass furnaces output using satistical and neural computing methods[END_REF], neural networks are implemented to compute the weekly production of a glass furnace according to the production planning. Several tools are used (linear or non-linear regressions, ANNs systems with retro-propagation algorithms, etc.), as well as a 175 elements database, to accurately predict the furnace output for a given input. Even though the mentioned project does not exactly aim at the description of the physical phenomena inside an furnace, it shows that possible applications can be reached; some precisions will be given on this point in the next section.

• Expert Systems methodology is to reproduce the cognitive mechanisms of an expert in a particular area. It lies on a database of facts and rules, and on an inference motor, using the database to produce new facts and thus providing a decision-support tool.

• Fuzzy-logic based learning methods present some similarities with expert systems. They are based on a fuzzy partition of the definition ranges of the input parameters and on inference processes that allow the generation of extrapolated system responses, in a similarly fuzzypartitioned output domain. These methods were mainly applied to the case of control systems of the furnaces. The control parameters are then tuned through an optimization step carried out with Genetics Algorithms, considering the following criteria: energy (fuel consumption), quality (number of bubbles or of not-melted elements in the resulting glass production) (see [START_REF] Pina | An operation system for industrial processes : application to a glass furnace, 9-12 juillet[END_REF] and [START_REF] Pina | A glass furnace operation system using fuzzy modelling and genetic algorithms for performance optimisation[END_REF]). The details of the fuzzy logics approach and of the used genetic algorithm are fully explained in [START_REF] Costa | The application of fuzzy logic in automatic modelling of electromechanical systems[END_REF].

Learning method class is obviously quite attractive since system outputs can be obtained without an accurate knowledge of the physics of the involved phenomena (of course, the knowledge of system outputs can constitute the real challenge). Nevertheless, its inherent drawback is the necessary availability of databases related to the studied problem: the database size will then condition the quality of the resulting inference model. It is worth recalling that, moreover, the database will be used not only to build and train the model, but also to validate later its behavior and evaluate the accuracy of the proposed extrapolation.

Reactor network models

Reactor network modeling is a commonly used technique in Chemical Engineering, adopted usually to model various kinds of flows inside a closed container. In a way similar to the learning techniques, it does not require an accurate physical understanding of the involved phenomena. Unitary operating items from the Chemical Engineering area, such as Continuous Stirred Tank Reactors (CSTR), Plug-flow Reactors (PR), idle zones, bypass or recycling streams etc., are organized in a stream network, in such a way that the experimental output of the real system is reproduced. According to the information found in a literature overview, two types of approach are possible:

• The systemic approach relies on the knowledge of experimental input and output signals. From the formulation of mass balance equations for each unit or node in the network, an optimization problem is generally deduced involving two objectives: (i) minimizing the gap (quadratic difference) between the model and the experiment; (ii) designing a model that is as simple as possible, i.e. including the lowest number of elementary items. A key-point of the application of this method is of course the choice of an appropriate solution technique. Two kinds of variables are involved, either structural (existence of an element of the network) or parametrical (size of the reactors, flow-rates of the internal streams) characteristics of the network. There are two classes of experimental signals, either static or dynamic:

-Static case : the input stream has known and fixed flow-rates and concentrations. The reactor network thus tries to approximate the output concentration value, assuming a steady state for the whole system. This kind of study was illustrated in [START_REF] Montastruc | Modélisation et optimisation dun réacteur en lit fluidisé de déphosphatation deffluents aqueux[END_REF] with the modeling of a fluidized bed used for the cleanup of aqueous effluents and in [START_REF] Hocine | Identification de modèles de procédés par programmation mathématique déterministe[END_REF] with the representation of a decantation bed and of a natural gas distribution network. The used optimization methods were respectively a hybrid technique (genetic algorithm and sequential quadratic programming) and a Mathematical Programming technique (Branch & Bound) implemented in the SBB solution module, integrated within the GAMS environment ( [START_REF] Brooke | GAMS Users Guide[END_REF]).

-Dynamic case : the signal is subjected to a variation of the input concentration of any tracer. Dirac pulse and step injection are usually used to obtain the Residence Time Distributions, which are typical of each type of input variation. This kind of studies thus focuses on transient state analysis, which involves the formulation of differential mass balances. These models were used in [START_REF] Laquerbe | Modélisation des écoulements dans un local ventilé par une approche systémique[END_REF] (simulated annealing) and [START_REF] Hocine | Identification de modèles de procédés par programmation mathématique déterministe[END_REF] (GAMS-SBB) for the representation of flows in a ventilated room. The experimental signals are actually outputs obtained by CFD calculations, or arbitrarily designed to solve academic examples.

• The topological approach is mainly related to the gas flows due to combustions, which typically take place in the combustion chamber of the furnace, above the glass bath ("laboratory") (see [START_REF] Faravelli | A new procedure for predicting nox emissions from furnaces[END_REF] and [START_REF] Falcitelli | Modelling practical combustion systems and predicting nox emissions with an integrated cdf based approach[END_REF]). The development of this approach is based on the observation that the CFD methods are computationally very greedy. Particularly, they prevent from considering, in addition to the hydrodynamic aspects, the chemical kinetics of combustion in a way different from extremely simplified models, which do not provide satisfying results (as compared with the experimental trends). The main idea of this method then consists in decoupling both phenomena by representing, in a first step, the flows characteristics by a reactor network, and then by applying to this latter the realistic and complex chemical kinetics of contaminant production,such as NOx [START_REF] Faravelli | Kinetic modelling of the interactions between no and hydrocarbons in the oxydation of hydrocarbons at low temperatures[END_REF]. To find the reactor network corresponding to the studied system, the methodology first determines a mapping of the combustion chamber, focusing particularly on the spatial distributions of some parameters that perform a key-role in the above-mentioned chemical kinetics: the temperature (T) and the combustion air/fuel stoichiometry (σ). The hydrodynamics computed by CFD thus enables to generate T-fields and σ-fields. Then, a group of zones is determined within the combustion chamber, where the keyparameters T and σ present homogeneous values. This operation is manually initialized and is then run automatically in such a way that both mean and standard deviation values of T and σ within each zone do not exceed a maximal allowed variation (see [START_REF] Falcitelli | Modelling practical combustion systems and predicting nox emissions with an integrated cdf based approach[END_REF]). Finally, these homogeneous zones can be represented as ideal reactors: the nature of each reactor (CSTR or PR) is deduced from an analysis of velocity fields inside the corresponding zone. A uniform velocity distribution is interpreted as a continuous stirred tank reactor (CSTR) while an unidirectional velocity field typically represents a plug-flow reactor (PR). Besides, the determined zones are composed of the elementary volumetric cells that were initially used for the space discretization during the CFD calculations. An analogy can then be empirically established between the combustion chamber mapping and the reactor network, in order to complete the definition of this latter:

-Each zone volume corresponds to the associate reactor volume.

-The mass exchanges between each zone, computed as the sum of the mass exchanges between cells belonging to different zones, correspond to the flow-rates of the network streams.

-The mean values (computed through the formula proposed in [START_REF] Falcitelli | Modelling practical combustion systems and predicting nox emissions with an integrated cdf based approach[END_REF]) of temperature and stoichiometry on each zone correspond to the temperature and the composition in each ideal reactor, thus allowing to carry out the kinetics computations of NOx production.

The numerical results provided in [START_REF] Frassoldati | Determination of nox emissions from strong swirling confined flames with an integrated cfd-based procedure[END_REF] highlighted the validity of the kinetics calculated following this methodology, which proved to be in agreement with the experimental data available for real combustion chambers.

Both of the above-mentioned reactor network strategies seem relevant, but they differ on two particular issues. First, the topological approach has the advantage to associate each homogeneous zone to a physical region in the furnace, and therefore to the operating parameters (temperature and composition) that are necessary for the kinetics computations. Second, [START_REF] Frassoldati | Determination of nox emissions from strong swirling confined flames with an integrated cfd-based procedure[END_REF] showed that the topological method is more accurate to predict the combustion kinetics, when the number of reactors in the network is high: this involves high computational times to determine a sufficiently extended network and then to carry out the kinetics calculations over it. In comparison, even if the systemic approach does not focus on the operating parameters of each reactor, it considers the network simplicity as an important aspect, which renders easier the following kinetics calculations. Thus, a reliable method, associating operating parameters to the reactors, would guarantee a large efficiency to the systemic approach.

Preliminary analysis

A soda-lime-silica glass is obtained by heating raw materials composed mainly by silica, SiO 2 and carbonated elements, Na 2 CO 3 , and CaCO 3 . Thermal energy comes from combustion chamber just above the glass bath and in a least part from electric energy directly in the glass bath. The bubbling device is generally used to increase locally the heat flux (see for more details [START_REF] Trier | Glass furnaces, design construction and operation[END_REF]). This chemical process leads to a huge production of bubbles. The quality requirements need to remove this bubble population by a chemical process activated by temperature [START_REF] Pigeonneau | Coupled modelling of redox reactions and glass melt fining processes[END_REF][START_REF] Pigeonneau | Mass transfer of rising bubble in molten glass with instantaneous oxidation-reduction reaction. Part 1: Determination of Sherwood number[END_REF]. So, the size and number of bubbles at the end of the furnace and the quantity of unmelted particles as well are important to know. The environmental requirements needs also to know the quantity of pollutant elements as SO x , NO x , and CO 2 .

In this study, two accuracy levels were considered. On the one hand, the furnace is viewed as a single "black-box" system: the input and output

GLASS MANUFACTURING PROCESS

Air (m Air , T)

Gas/Fuel (m G/F , T)
Raw materials (SiO 2 , Na 2 CO 3 , CaCO 3 ) + T

Bubble injectors

Glass : T,m Glass , quality (bubbles, non-melted elements)

Burnt gas: T, composition (SO x , NO x , CO 2 )

GLASS MANUFACTURING PROCESS

Air (m Air , T)

Gas/Fuel (m G/F , T)
Raw materials (SiO 2 , Na 2 CO 3 , CaCO 3 ) + T

Electric energy

Glass : T,m Glass , quality (bubbles, non-melted elements)

Burnt gas: T, composition (SO x , NO x , CO 2 )

Figure 1: Global description of the glass manufacturing furnace.

of the system are defined such as illustrated in Fig. 1. For this type of representation, the learning methods seem to be the only ones that can be applied. The drawback, already mentioned in the previous section, is the requirement of gathering a sufficient data number to create a reliable database on which the training may be carried out. This involves a large experimental campaign, performed on many real glass manufacturing plants, including the determination of ways of measuring the variables and criteria defined in Fig. 1, which is not a simple task under the extreme conditions governing the furnaces.

On the other hand, within the perspective of applying a reactor network modeling approach, a more accurate decomposition of the whole system is proposed in Fig. 2. The glass manufacturing process is decoupled into two entities, i.e., glass bath and combustion chamber. This decomposition would then allow taking advantage of the CFD simulations. These are carried out separately on each part of the system, imposing particular boundary conditions to the interface glass bath -combustion chamber. Therefore, the objective would be to generate two reactor networks, each one dedicated to the representation of a furnace section. This methodology would require, in order to model a whole furnace, the availability of data obtained from CFD simulations carried out on both zones of a same furnace.

As a result, the systemic approach seems to be a promising strategy to solve this kind of problems. In this framework, this paper constitutes the first step of the global methodology. Such as mentioned in the previous section, it is necessary to integrate kinetics or energy aspects to the network obtained from the CFD computations, which constitutes long-term research perspectives. The preliminary step which is presented here, has the objective to validate the feasibility of this strategy, by determining reactor networks by specific Residence Time Distributions, obtained by simulations carried out in Saint-Gobain Recherche Center. The numerical computations have been performed by use of a commercial CFD software tool. The glass is assumed as an incompressible fluid where the buoyancy force due to the thermal dilatation is taken into account under Boussinesq approximation [START_REF] Gray | The validity of the Boussinesq approximation for liquids and gases[END_REF]. The heat and mass transfer is mainly produced by the free convection due to the large scale thermal gradient occurring on the top surface of the bath (see for more details [START_REF] Pigeonneau | Natural convection of a high Prandtl number fluid in a long cavity heated from above[END_REF]). The thermal dependence of the dynamical viscosity is taken into account using an exponential behavior as a function of temperature [START_REF] Scholze | Glass. Nature, Structures and Properties[END_REF]. The radiative process is also very important in the glass industry since temperature reaches 1600 • C in the glass bath. A glass is a semitransparent material where the radiative transfer is treated as a diffusion process for which the radiative conductivity is a function of the cube of temperature [START_REF] Viskanta | Heat transfer in semitransparent solids[END_REF]. The thermal boundary conditions between the glass and furnace walls are determined by the knowledge of the design of furnace. The raw material is seen as a thermal sink and a mass inlet. The computation is carried out by a specific fitting of temperature obtained numerically and in the industrial plant. The residence time distribution is achieved after the CFD calculation. The RTD is determined by looking the answer of Heaviside function. The time derivative gives directly the RTD.

Systemic approach implementation

The set of the following computations presented in this section aims at validating the systemic approach by reactor networks implemented in this study. The objective is to show that a structure of reactor network and the associated parameters can be derived from the experimental signal corresponding to the hydrodynamics phenomena in the glass manufacturing furnaces. This experimental signal was obtained by CFD simulations that were carried out, for both combustion chamber and glass bath, by the 'Glass Melting' research team of Saint-Gobain Recherche.

Details of the method

The main principles of the flow model identification were briefly recalled in the bibliography review section. All the aforementioned studies deal with dynamic case, i.e. the (CFD-based) experimental signal is a Residence Time Distribution (RTD) resulting from a classical input variation (Dirac pulse or step injection): it is thus necessary to reproduce, with an adapted model, the output signal evolution throughout time. It is worth noting that, despite the existence of various chemical compounds in the furnace (in both combustion chamber and glass bath, such as exhibited in figure 1), only one is considered in the framework of this study. This chemical is thus viewed as a tracer (involving no chemical reaction) for which concentration variations are observed in order to get the Residence Time Distributions used for the following computations.

As it was mentioned previously, the model equations represent the mass balances at each node of the network and for each ideal reactor involved in the structure. However, because of the dynamic nature of the problem, the considered mass balances imply differential formulations: of course, the time scale discretization results in a number of equations that must be verified at each instant t = t 1 , . . . , t F . The method begins defining a superstructure, which gathers all the elementary operations likely to constitute a final network that accurately models the studied flow. Binary variables y i are associated to each potential unit: they are equal to 1 if the unit exists, and to 0 otherwise. The variables representing the flow-rate and concentration in each stream (resp. Q j and C j ), as well as the ones associated to the reactor volumes (V i ) are real variables. For the sake of illustration, figure 3 provides a hypothetical superstructure.

Therefore, the optimizing tool will have to modify and schedule the initial superstructure in order to minimize a criterion that takes into account (i) on the one hand, for each time discretization slot, the distance between the CFD output concentration Cf exp and the modeled one Cf mod ; (ii) on the other hand, the structure complexity, i.e. the number of elementary units belonging to the final network. This objective can thus be formulated such as the following equation:

y 1 y 2 -V 2 y 3 -V 3 y 4 -V 4 y 5 -V 5 y 6 -V 6 y 7 Q 1 , C 0 Q 2 , C 0 Q 2 , C 1 Q 3 , C 0 Q 3 , C 2 Q 3 , C 3 Q 4 , C 0 Q 4 , C 4 Q 4 , C 5 Q 5 , C f Q in , C f Q in , C in y 1 y 2 -V 2 y 3 -V 3 y 4 -V 4 y 5 -V 5 y 6 -V 6 y 7 Q 1 , C 0 Q 2 , C 0 Q 2 , C 1 Q 3 , C 0 Q 3 , C 2 Q 3 , C 3 Q 4 , C 0 Q 4 , C 4 Q 4 , C 5 Q 5 , C f Q in , C f Q in , C in Figure 3: Superstructure example. minZ = k[ t F t=t 1 (Cf exp (t) -Cf mod (t)) 2 ] + nj j=1 y j (1)
In this equation, k is a constant scalar balancing the relative weights of each term within the global criterion. For instance, the higher k, the greater the importance given to the model-experiment distance minimization; in the same time, the optimization tool will disregard the structure complexity reduction. For the following computations, the k value was determined after a preliminary sensitivity analysis study which showed very clear variations (between a very simple network that did not optimize well the quadratic difference criterion and complex solutions that did not bring significant improvements of the quadratic difference criterion) and allowed an easy choice of the coefficient: the k value is taken equal to 5 • 10 3 .

Besides, the model is submitted to equality constraints that consist of differential mass balances on the network (a complete description of the equations is available in Appendix 1). It must be pointed out that some variables and constraints are subjected to the existence of the element (stream or reactor) they are associated with: this requires a dynamic constraint handling and a valid bounding technique for the real variables within the optimization problem.

The resulting problem is of Mixed Integer NonLinear Programming (MINLP) nature since it involves both binary and continuous variables, with nonlinear equations. This kind of problem is of the most complex ones because of the discretization of the search space due to the integer variables. At the same time, nonlinearities may cause local optima in non-convex cases. It was shown in [START_REF] Hocine | Identification de modèles de procédés par programmation mathématique déterministe[END_REF] that the most efficient solution technique is a Mathematical Programming Branch & Bound procedure (implemented in the SBB module of the GAMS environment [START_REF] Brooke | GAMS Users Guide[END_REF]). This method was therefore chosen for this study.

The structural and parametric identification of a network requires the following data:

• Global reactor volume (i.e. the sum of the ideal reactor volumes in the network is equal to the real container volume simulated by CFD);

• Flow-rate (static) and concentration (dynamic) of the input and output streams;

• Residence Time Distribution on the considered tank (combustion chamber or glass bath).

It must be pointed out that the RTD curves provided by the CFD simulations correspond to a Dirac pulse. However, the mass balance differentiation for Dirac pulses involves complex mathematical formulations. Since the Dirac pulse corresponds to the step injection derivative, the response (output signal) is obtained by a simple integration of the former signal response over time. Thus, for the following calculations, the output signal was integrated in order to get a response to a step injection, whose magnitude is equal to 1. This results in a simpler model formulation.

Necessarily , the RTD curves computed from CFD computations have a large number of points corresponding to the time discretization (about 20,000 points for each RTD). Yet, the chosen formulation leads to write the mass balance equations for each time slot in the model. Thus, to avoid a huge number of constraints that would heavily penalize the solution procedure, the interval between each time slot was increased in order to finally use only 100 points for the time discretization (sensitivity tests that are not presented here showed that this number allowed a correct solution accuracy without slowing down the computations).

The two next subsections present the computational results obtained on the glass bath and the combustion chamber respectively.

Computations on the glass bath

Figure 4 shows the RTD associated to the glass bath in its both initial (Dirac pulse) and integrated (step injection) form. The value of the glass bath total volume is 798.05 m 3 , and the input flow-rate is equal to 12.18 m 3 /h. Table 1 sums up the superstructures tested to reproduce the flow associated to the above RTD, as well as the corresponding solutions and the quadratic distance between model and CFD-based data. The curves for these models are available in Appendix 2. It appears that superstructure 3 is the most adapted to provide a modeled output similar to the simulated one, with a quadratic difference approximately equal to 10 -6 .

An immediate comment is that, whatever the initial used superstructure, the solutions propose a separate reactor (CSTR or PR), without any stream to feed it, characterized by a volume varying from 501 to 517 m 3 . This volume is commonly called 'idle' volume, or 'useless' volume. It is worth underlying that this volume is not likely to represent a particular unique zone in the glass bath, but more certainly a set of regions, spread over the whole volume, where local dynamic behaviors (for instance, very local recirculation streams) do not participate to the global flow.

Nevertheless, even though it is commonly admitted by glass manufacturing experts that there exists some zones that are no involved in the global flow, the 'useless' volume highlighted by these computations fills a large part of the total volume (between 63 and 65 percent, depending on the considered solution). Moreover, the repeatability of this result, obtained independently from the initial superstructure, seems to confirm the trend. Two additional analyses were carried out in order to justify it.

The former study aims at proving that this is not a numerical problem that constrained the optimizing procedure to keep trapped on this solution. Thus, similar computations are carried out, only by modifying the global reactor volume: 500 m 3 instead of the initial value of 798.05 m 3 . The result is identical to the previous one (structure and volumes of the ideal reactors), except for the 'useless' volume that is equal to 205.46 m 3 : the difference between the two obtained 'useless' volumes is thus equal to the difference between the two global volumes (798.05 -500 = 298.05 m 3 ). This confirms the purpose of the 'useless' volume, which can be seen as a slack variable for the global volume constraint:

nj j=1 V j = V T (2)
This first step validates the existence of a 'useless' volume. The latter study focuses on the physical meaning of this trend. Two simulations are now carried out (with any spreadsheet software) in order to observe the response of a single CSTR reactor. Its volume is chosen equal to 800 m 3 (the global volume) and 300 m 3 (the 'useful' volume) successively. Figure 5 gives an illustration of both outputs, compared to the CFD simulation: it clearly shows that, even if the curve corresponding to the 300 m 3 reactor does not exactly fit the curve obtained by CFD, it constitutes a much better approximation than the 800 m 3 reactor curve. Thus, the 'useless' volume determined by network identification is justified.

The solutions obtained from the different superstructures have various reactors arranged in series. Additional computations (not presented here) were performed to verify that the order of the CSTRs does not have any effect on the output concentration.

Combustion chamber

The RTD used in this section is not associated to the whole combustion chamber but to a cross section (perpendicular to the glass flow in the below bath) that includes one of the fuel-air injector of the considered 'float' furnace: the zone in front of this injector is the location where the flame develops, up to the smoke outlets towards the corresponding regenerator. Figure 6 presents the RTD curve associated to this combustion chamber section in its both initial (Dirac pulse) and integrated (step injection) form. The total section volume is equal to 90.63 m 3 and the flow-rate is equal to 1.942 m 3 /s. This flow-rate is higher than in the glass bath case, which will imply lower characteristic residence times.

Just in the same way as for the glass bath study, various initial superstructures were tested and their respective solutions are presented in table 2. The quadratic differences shown in table 2 highlight that the first solution is undoubtedly the best one. Furthermore, the curves showing CFD simulation vs. model curves are available in Appendix 3 and confirm the low consistency of the other solutions (2 to 4) towards the 'experimental' curve. The second 'swelling' observable in the original Dirac pulse RTD reproduces an important recirculation stream, which is identified whatever the used superstructure. However, this trend causes a sudden change in the slope of the integrated RTD curve, which is not reproduced by superstructures 2 and 3. Besides, the slight response time-shift with respect to the input injection time (minimum residence time) can be easily assimilated to a plug-flow effect, and can only be modeled by at least 3 CSTRs set up in series (this explains why superstructure 4 is unable to reproduce this effect).

Moreover, a second trend that is similar in all four determined solutions is the existence of an almost zero flow-rate stream (direct way for structures 2 and 3 and in reverse way for structures 1 and 4). This stream passes through various CSTRs that fill a global volume of between 76 and 79 m 3 . Once more, this observation indicates that a 'useless' volume can be considered, with a very important relative size towards the global volume (about 85 percent). A more accurate analysis of the first superstructure solution will prove it and lead to further network simplifications.

Firstly, a simulation is carried out to determine the response of solution 1when the weak flow-rate stream and the two associated CSTRs are removed. This leads to consider an 'idle' volume equal to 76.26 m 3 , and the 0.01 m 3 /s stream is added to the recycling stream (Q = 1.12 + 0.01 = 1.13 m 3 /s). The simulation then provides a quadratic difference value equal to 2.09 • 10 -3 (instead of the initial 6.75 • 10 -6 ). Even though the accuracy loss seems quite important, the response curve given in Figure 7 proves that the main features of the experimental output signal are respected (minimum residence time, slope variation). This confirms the 'useless' volume existence. Besides, the magnitude of this useless volume towards the global volume can be explained more easily than in the glass bath case: the studied flow is a flame development flow, which does not necessarily fill a large volume since the main energy transfer is realized through radiation. The large size of the 'useless' volume is therefore less surprising here. Moreover, a second simplification can also be proposed: since the second plug-flow reactor volume is very low with respect to the global container one (0.72 m 3 ), this reactor can also be removed, adding the 0.72 m 3 to the 'useless' volume and the corresponding stream flow-rate to the recycling stream (now equal to 1.25 m 3 /s) The quadratic difference in the resulting solution is equal to 2.27 • 10 -3 , which represents an insignificant increase (and the response curve is almost identical to that shown in figure 7). A second simplified structure of solution 1, illustrated in figure 8, is thus obtained. This former provides a very satisfactory output. This network was not found by the optimization module because of the higher quadratic difference, but it may be obtained by increasing the weight associated to the structure complexity reduction in the optimization criterion, i.e. decreasing the k factor in equation [START_REF] Carvalho | Modelling and optimization of the no formation in an industrial glass furnace[END_REF]. Finally, in order to confirm the validity of the found solution, an optimization step is carried out, using the final determined solution as the initial superstructure. The same solution structure is found, which thus confirms the optimality of the resulting network. Only the parametrical variables (volumes and flow-rates) are modified, to obtain a solution (figure 9) that slightly improves the quadratic difference (1.09 • 10 -3 ).

Q = 2.
A last simulation was run to show that the plug-flow reactor is essential. Figure 10 then proves that, by removing the PR and the associated stream, the output curve does not reproduce the minimum residence time, which is one of the main features of the experimental response. The above-mentioned plug-flow effect is thus demonstrated.

Conclusions and perspectives

This study presents a new way of modeling glass manufacturing processes, based on a systemic approach commonly used in Process Systems Engineering. The involved processes exhibit complex phenomena occurring in the glass manufacturing system. The lack of sufficient experimental data renders difficult the use of learning methods. Thus, the proposed methodology first relies on the decomposition of the whole furnace in two entities, i.e., the combustion chamber and the glass bath. Then, it uses simulations carried out by CFD tools as experimental data to optimize a reactor network to obtain a simplified representation of the hydrodynamics in each furnace zone.

The proposed framework is an alternative computational strategy to designoriented CFD computations: the system representation is a "black-box" model, where physical and chemical phenomena occurring inside the furnace (which have an obvious influence on the fluid dynamics) are not taken into account explicitly. The advantage of the proposed method is that it only considers input/output informations in order to build a model that reproduces faithfully the Residence Time Distribution. The presented study validates this first step, while future work must be devoted to the extension of the model and to the consideration of kinetics computations. More precisely, the • The preliminary computations provided very satisfactory results, since in furnace region cases, a reactor network giving an output similar to the experimental response can be identified. This proves the method ability to reliably modeling the hydrodynamics phenomena in the glass bath or over a section of the combustion chamber. Moreover, the obtained results exhibit some trends proper to the glass manufacturing furnace, and, more particularly, the existence of 'useless' volumes that do not participate to the global flow. These numerical trends were confirmed by further studies, and therefore, they should lead glass manufacturing experts to go into the details of the physical reasons of the mentioned behavior (strong recirculation streams, particular flows local to narrow zones • • •).

• Besides, the literature review developed in section 2 showed that a quite similar methodology was used in the dedicated literature : the reactor network representing the flow in the considered container was yet derived from a topological approach. This strategy is particularly useful to carry out detailed and complex kinetics computations to study the contaminant production (particularly the NOx species that appear in the combustion chamber). In this sense, the systemic approach drawback is clearly the issue of associating operating parameters to the reactors in the network, in order to perform the above-mentioned kinetics computations in a second step. The clear advantage of the systemic approach over the topological one is the complete automation of the modeling and solution process: only rough data (global volume and inlet flow-rate) are needed, as well as an adapted initial superstructure. Neither variable initialization nor preliminary manual design of the possible network are necessary and the final reactor network is automatically generated by the optimization process. Moreover, the computational time necessary to find a solution (a few seconds) with the proposed method is much lower than the times reported in the studies applying the topological methodology (several minutes). Finally, the systemic approach guarantees an optimal solution in terms of similarity between model and (CFD-based) experimental data, and furthermore includes the concept of model simplicity that might be useful and CPU-time saving for further computations (such as, for instance, the previously mentioned kinetics calculations).

The perspectives of this work should particularly focus on the following points:

• In addition to the mass balances formulated on the reactor network, energy balances should be also written in order to identify the temperatures associated to each reactor. In this sense, experimental verification of energy balances, available for some real furnaces, would lead to a more complete structural and parametrical identification (over either the glass bath or the combustion chamber).

• Another way to produce this kind of structure would be to apply a temperature variation (Dirac pulse or step injection) on the furnace input. The CFD simulations of the output temperature would then provide another signal on which to carry out new model identification by systemic approach.

• Besides, it is worth noting that energy balances formulated on the reactor network require some knowledge about chemistry compositions, and thus about the kinetics (to determine the reaction enthalpies). A first step could thus use simplified kinetics mechanisms (for instance, the Zeldovich mechanisms in the combustion case). In a second step, a more complex kinetics mechanisms could be integrated to the reactor network (see [START_REF] Faravelli | Kinetic modelling of the interactions between no and hydrocarbons in the oxydation of hydrocarbons at low temperatures[END_REF]).

Finally, even though the application zone of the mentioned perspectives seems to be mostly the combustion chamber, the ultimate objective is to build a global representation of the glass manufacturing furnaces, by associating the reactor network to both furnace entities. The critical point of this strategy would then be located on the glass bath -combustion chamber interface.

Appendix 1 MINLP formulation for reactor network modeling

The initial data necessary to generate a reactor network that fits the experimental-CFD output curve over a discretized time period {t 1 , . . . , t F } are the following ones:

• the volume of the considered container V T ,

• the inlet flowrate Q tot and concentration C in for each instant t (here, C in (t) = 1, t = t 1 , . . . , t F ),

• the oulet concentrations Cf exp (t), t = t 1 , . . . , t F (obtained by CFD simulation),

• an initial superstructure.

The above-mentioned superstructure involves the unitary elements that will be characterized by the equations of the model. In addition, general constraints such as the upper bound of the container volume, as well as the inlet and outlet mass balances, have to be defined whatever the superstructure.

Finally, a dynamic bounding for reactors volume and streams flow-rate variables must be implemented, in order to set the associated values equal to zero when the considered items do not exist (i.e. the corresponding binary variables are equal to 0). The used variables are:

• y j is a binary variable denoting the existence of an item (j = 1, . . . , nj),

• V j is a continuous variable representing the volume of item j (j = 1, . . . , nj),

• Q k is a continuous variable representing the flow-rate of stream k (k = 1, . . . , nk),

• C i (t) is a continuous variable representing the concentration at instant t (t = t 1 , . . . , t F ) and at a specific point of the superstructure named i (i = 1, . . . , ni),

• Cf mod (t) is a continuous variable representing the ouput concentration at instant t (t = t 1 , . . . , t F ), generated by the model.

The objective function accounts for (i) the minimization of the modelexperimental quadratic distance, for each instant t; and (ii) the complexity of the resulting network, calculated as the sum of items involved in the final structure:

minZ = k[ t F t=t 1 (Cf exp (t) -Cf mod (t)) 2 ] + nj j=1 y j . (3)
v The constraints are described straightforward:

• Constraint on the global container volume.

nj j=1 V j = V T (4)
• Global mass balance at the inlet node.

Q tot + k in Q k = kout Q k (5)
Where k in / k out represent the indexes for streams entering/leaving the inlet node, apart from the inlet stream.

• Partial mass balance on the inlet node.

∀t ∈ {t 1 , . . . , t F },

Q tot C in + k in Q k C i (t) = kout Q k C i (t) (6)
Where k in / k out represent the indexes for streams entering/leaving the inlet node, apart from the inlet stream. i and i represent the points in the superstructure associated with the streams entering (k in ) or leaving (k out ) the inlet node.

• Partial mass balance on the outlet node.

∀t ∈ {t 1 , . . . , t F },

k in Q k C i (t) = Q tot solution1Cf mod (t) + kout Q k C i (t) (7)
Where k in / k out represent the indexes for streams entering/leaving the outlet node, apart from the outlet stream. i and i represent the points in the superstructure associated with the streams entering (k in ) or leaving (k out ) the outlet node.

• Mass balance on a CSTR module. ∀t ∈ {t 1 , . . . , t F }, ∀j ∈ {1, . . . , nj}, C i (t) = hQ k z 1 + V j (C i (t-1)-C i (t-1))+C i (t-1)

Where k represent the index of the stream passing through the CSTR module j; i is a superstructure point before the CSTR module and i is a superstructure point after the CSTR module; h is a small time slot used for differentiation and z 1 is a small constant.

• Mass balance on a plug-flow module.

∀t ∈ {t 1 , . . . , t F }, ∀j ∈ {1, . . . , nj}, C i (t) = C i (t)sigmoid(z 2 (t-

V j z + Q k )) (9)
Where k represent the index of the stream passing through the plugflow module j; i is a superstructure point before the plug-flow module and i is a superstructure point after the plug-flow module. Function sigmoid is used for the representation of a plug-flow effect when a step injection is used as an input signal. z 2 is constant used to tune the sigmoid function.

• Dynamic bounding of the items volume. ∀j ∈ {1, . . . , nj}, V j ≤ y j V T [START_REF] Laquerbe | Modélisation des écoulements dans un local ventilé par une approche systémique[END_REF] • Dynamic bounding of the stream flow-rates.

∀k ∈ {1, . . . , nk}, Q k ≤ M

j k y j (11)
Where j k represent the indexes of the items that are crossed by stream k. The parameter M is that of the classical big-M technique, used to define an upper bound on a variable while avoiding the formulation of bilinear constraints (which strongly penalize the solution process).

Appendix 2

Computations on the Glass Bath

Figure 2 :

 2 Figure 2: Decoupling the glass manufacturing furnace.

Figure 4 :

 4 Figure 4: Glass bath (initial and integrated) RTD.

Figure 5 :

 5 Figure 5: CSTR simulations with V = 300 m 3 and V = 800 m 3 .

  pulse RTD -disc. step inj.

Figure 6 :

 6 Figure 6: Combustion chamber (initial and integrated) RTD.

Figure 7 :

 7 Figure 7: Response for solution 1 first simplification (Simpl 1 ).

3 Figure 8 :

 38 Figure 8: Final structure,after simplifying solution 1.

3 Figure 9 :

 39 Figure 9: Final solution after optimization.

Figure 10 :

 10 Figure 10: Proof of the essential presence of the PR in the final structure.

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Comparison between CFD and model results for Superstructure 1, Glass Bath.

Fig. 13 :Fig. 14 :Fig. 15 :Fig. 16 :

 13141516 Fig. 13: Comparison between CFD and model results for Superstructure 3, Glass Bath.

Fig. 17 :Fig. 18 :

 1718 Fig. 17: Comparison between CFD and model results for Superstructure 3, Combustion Chamber.

Table 1 :

 1 Glass bath results, according to the initial superstructure.

Table 2 :

 2 Combustion chamber results, according to the initial superstructure.

List of Figures