Ignacio Gómez-García-Toraño

Laurent Laborde
email: laurent.laborde@irsn.fr

Modelling reflooding of intact core geometries in ASTEC V2.1: improvements and validation on PERICLES experiments

The injection of water into core is a key Accident Management strategy in water-cooled nuclear reactors.

Among the flow regimes implicated during reflooding, the most significant one in terms of heat transfer removal is the transition boiling. This mechanism takes place over a short distance near the quench front, which is difficult to capture when using coarse meshes, typically in integral severe accident codes. In the ASTEC V2.1 integral severe accident code, the user can select among two concurrent two-phase thermalhydraulics modelling (five-or six-equation scheme) and for each one, two alternative reflooding models to describe transition boiling region (TBR). On one hand, the DROPLET model calculates a large heat transfer coefficient over a user-defined distance above the quench front. On the other hand, the EXPCHF model calculates the heat flux downstream by imposing the Critical Heat Flux at the quench front and prescribing an exponentially decreasing heat flux along a flow-dependent length.

This article aims at selecting the most appropriate reflooding model for the evaluation of the heat transfer associated with transition boiling in ASTEC V2.1 both for the five and the six-equation hydrodynamic model. For that purpose, experimental data from the PERICLES facility dealing with the reflooding of intact core geometries have been used.

After an exhaustive analysis of the different reflooding models and hydrodynamic models, this study has found that the DROPLET model should be selected when using the five-equation scheme, since it provides enough rod precooling above the quench front. However, the EXPCHF model should be selected when using the six-equation scheme, since it benefits from the larger amount of liquid available near the front. This last configuration gives the best results on the simulated PERICLES tests and is thus recommended for reactor applications.

Highlights

• ASTEC V2.1 reflooding model for rod-like geometries is validated against PERICLES facility.

• DROPLET reflooding alone is recommended for the five-equation scheme.

• EXPCHF reflooding model is recommended for the six-equation scheme.

• : interfacial friction between the gas and liquid phases in the six-equation scheme (Pa/m) • FNR : non-stratified interfacial friction term in the five-equation scheme (kg/m 4 )

• g: gravity (m/s 2 )

• HD: heat transfer coefficient to the DROPLET model (W/m 2 K)

• I: inertia term (Pa/m)

• L1: characteristic length of the TBR associated with the EXPCHF model (m)

• Re: Reynolds number (-)

• TBO: burnout temperature (K)

• Tliq: liquid temperature (K)

• TMFB: minimum film boiling temperature (K)

• Tsat: saturation temperature (K)

• Twall: wall temperature (K)

• v , v : gas, liquid velocity (K)

• Z: elevation (m)

• ZD: characteristic height associated to the DROPLET model (m) 2 Greek letters:

• α: void fraction (-)

• αD: threshold void fraction associated to the DROPLET model (-)

• ΔV = v -v : drift velocity (m/s)

• φ

: heat transfer along the TBR using the i reflooding model (W/m 2 ) • μ : dynamic viscosity (Pa s)

•

: density (kg/m 3 ) • σ : surface tension (N/m) Subscripts, superscripts:

• G: gas phase

Introduction

The injection of water into the core, frequently called core reflooding, is one of the prime Accident Management strategies in all current Light Water Reactors [START_REF] Hermsmeyer | Review of current Severe Accident Management ( SAM ) approaches for Nuclear Power Plants in Europe requirements[END_REF]. Reflooding starts at any time after the Reactor Control Protection System switches on the Emergency Core Cooling System. The action prevents core degradation if launched before the Core Exit Temperature exceeds a given threshold, generally 650ºC corresponding to Peak Cladding Temperatures of about 950ºC (NEA-OECD, 2010). Beyond that limit, the focus of the actions shifts towards the mitigation of the accident, ranging from the core damage minimization to prevention of vessel failure. Thus, reflooding is considered as one of the most important thermal-hydraulic phenomena in the context of nuclear safety.

As a matter of fact, this mechanism allowed the mitigation of the severe accident at Three Miles Island Unit 2, despite the core underwent significant degradation [START_REF] Broughton | A Scenario of the Three Mile Island Unit 2 Accident[END_REF].

Considerable experimental efforts have been invested over the last decades to study the heat transfer mechanisms taking place within the core during reflooding. Rod bundle reflooding experiments at the FLECHT-SEASET [START_REF] Lee | PWR FLECHT SEASET unblocked bundle, forced and gravity reflood task data evaluation and analysis report[END_REF], FEBA [START_REF] Ihle | FEBA -Flooding Experiments with Blocked Arrays Evaluation Report[END_REF], PERICLES [START_REF] Digonnet | PERICLES Programme : boil-up, boil-off and reflooding high pressure experiments in a PWR assembly[END_REF][START_REF] Housiadas | Two-dimensional quench front progression in a multiassembly rod bundle[END_REF] and the RBHT [START_REF] Hochreiter | RBHT Reflood Heat Transfer Experiments Data and Analysis[END_REF] facilities have contributed to enhance the understanding on such mechanisms for intact core geometries using a wide range of thermal-hydraulic conditions (e.g. system pressure, injection flow rate, initial temperature, inlet subcooling, power). Likewise, rod bundle reflooding experiments at the LOFT FP-2 [START_REF] Cronenberg | Hydrogen generation behaviour in the LOFT FP-2 and other experiments: comparative assessment for migrated severe accident conditions[END_REF], CORA [START_REF] Schanz | Information on the evolution of severe LWR fuel element damage obtained in the CORA program[END_REF] and QUENCH (Sepold et al., 2001;[START_REF] Steinbrück | Synopsis and outcome of the QUENCH experimental program[END_REF] facilities have contributed to the same goal for partially degraded cores under low pressures, both conditions being representative of a reflooding scenario during the early in-vessel phase of a severe accident. As a matter of fact, a thorough analyses of the latter experiments contributed to identify criteria for the successful coolability of partially degraded cores [START_REF] Hering | Integration of New Experiments into the Reflood Map[END_REF][START_REF] Hering | Degraded core reflood: Present understanding and impact on LWRs[END_REF].

Visual studies representative of sub-channel of an intact core performed by [START_REF] Ishii | Flow visualization study of inverted annular flow of post dryout heat transfer region[END_REF][START_REF] Obot | Two-phase flow regime transition criteria in post-dryout region based on flow visualization experiments[END_REF] have shown that the governing heat transfer regimes occurring during reflooding can be split in nuclear boiling (upstream the quench front), transition boiling (just above the quench front) and dispersed film boiling (well above the quench front). Of special importance is the Transition Boiling Region (TBR), which covers a length of up to few centimetres along which there is a non-stable two-phase flow leading to large heat exchanges with the wall. This length depends on the existing flow regimes above the quench front, which depend on the void fraction at the quench front (linked to the pressure, inlet subcooling, residual power and flow rate). Low system pressures, low flow rates and/or high residual powers lead to a significant increase of the void fraction as soon as the flow at the quench front becomes saturated, the flow above the front ending up being highly dispersed [START_REF] Yadigaroglu | Modeling of reflooding[END_REF] and the TBR becoming increasingly shorter. A sketch of these experimental observations is depicted at the centre of Figure 1, the depiction of flow patterns being inspired in the figure from [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF].

The modelling of reflooding is a challenging task not only for severe accident codes, but also for bestestimate system thermalhydraulic codes like CATHARE (France), ATHLET (Germany), TRACE and RELAP5 (USA) or MARS (Korea), as evidenced in the BEMUSE benchmark [START_REF] Mendizábal | Post-BEMUSE Reflood Model Input Uncertainty Methods (PREMIUM) Benchmark: Final Report[END_REF].

As far as the modelling in integral severe accident codes is concerned, the main challenge is to model the significant heat transfer taking place along the short TBR, while keeping coarse axial meshes (20-30 cm) to limit the CPU time of a complete reactor sequence. In this sense, one of the most interesting approaches to overcome this issue was formulated was the one formulated by [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF].

Following experimental observations of Ishii and co-workers, such model states that the heat transfer along the TBR follows a decreasing exponential function (over an empirical length), the maximum value being the Critical Heat Flux (CHF) at the quench front location. This model was implemented in the two-fluid and six-equation ICARE-CATHARE mechanistic severe accident code (developed at IRSN), which, in contrast to integral codes, is dedicated to the more detailed description of circuits thermalhydraulics and in-vessel phenomena. With this tool, the model was validated using data from the PERICLES facility at CEA (France) and RBHT facility at PSU (USA) [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF]. The model is also showed to be valid for partially degraded cores, as shown by the validation of ICARE-CATHARE on the QUENCH experiments at KIT (Germany) (Chikhi et al., 2012;Chikhi and Fleurot, 2012).

However, the transposition of this modelling approach to integral severe accident codes is not straightforward. Indeed, integral severe accident codes need to be fast running tools with the highest possible accuracy to support accident management studies. Therefore, they introduce simplifications in the number of equations, fluids and involved terms in the mass, momentum and energy balance equations in comparison to best-estimate system thermalhydraulics codes, which use a two-fluid model involving the resolution of two momentum equations. The integral ASTEC V2.1 severe accident code, developed at IRSN, uses by default a two-fluid and five-equation scheme to obtain the void fraction, pressure and temperature-velocity fields for the liquid and gas phases [START_REF] Chatelard | Main modelling features of the ASTEC V2.1 major version[END_REF][START_REF] Chatelard | ASTEC V2 severe accident integral code main features, current V2.0 modelling status, perspectives[END_REF].

However, a two-fluid and six-equation approach has been recently transposed from the IRSN in-house DRACCAR code [START_REF] Glantz | DRACCAR: A multiphysics code for computational analysis of multi-rod ballooning, coolability and fuel relocation during LOCA transients. Part One: General modeling description[END_REF] as to give a more accurate prediction of the system thermalhydraulics. This is in consonance with the MELCOR V2.2 integral severe accident code developed at SNL (Humphries et al., 2017), which also uses a two-fluid and six-equation scheme for the resolution of the system thermalhydraulics.

Currently, all reactor safety analyses conducted using ASTEC V2.1 are carried out with the fiveequation scheme and considering that the heat transfer along the TBR is given by the droplet projection model (in the following referred as DROPLET model). This model sets a large, uniform and constant heat transfer coefficient above the quench front providing that there is enough water in the fluid mesh.

However, it introduces important uncertainties in the calculations, since the value of such a coefficient quite depends on the local conditions in the vicinity of the front [START_REF] Yadigaroglu | Modeling of reflooding[END_REF]. Therefore, IRSN has started to conduct analyses using the referred exponential heat flux model proposed by [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF][START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF], referred as EXPCHF model in this article (DRACCAR model in ASTEC V2.1 series).

In this context, the present article aims at presenting the improved reflooding validation results due to the new six equation description and determining the most appropriate modelling choice concerning the heat transfer along the TBR (DROPLET or EXPCHF) for the five and six-equation hydrodynamic models, keeping in mind that a five-equation description may still be requested for specific computationally expensive sequences. For such a task, code predictions have been validated using data from the PERICLES facility dealing with the reflooding of intact PWR cores [START_REF] Digonnet | PERICLES Programme : boil-up, boil-off and reflooding high pressure experiments in a PWR assembly[END_REF].

Description of selected physical models by ASTEC V2.1

The CESAR and ICARE modules

The ASTEC code (Accident Source Term Evaluation Code) aims at simulating the progression of entire severe accident sequences in nuclear water-cooled reactors, from occurrence of the initiating event up to the release of radioactive elements out of the containment. Its main field of applications are Severe Accident Management (SAM) and Level 2 Probabilistic Safety Assessment (PSA2) studies. The code is modular, each module describing a set of physical phenomena, which may take place during a severe accident. A detailed description of the scope of each ASTEC V2.1 module can be found in [START_REF] Chatelard | Main modelling features of the ASTEC V2.1 major version[END_REF]. Within this work, only the CESAR and ICARE modules are used.

The CESAR module describes the thermal-hydraulics throughout the primary and the secondary circuit, including the reactor pressure vessel. Within the core, the CESAR thermal-hydraulics modelling (Gómez-García-Toraño and Laborde, 2019) is based on a two-dimensional two-fluid and a five (by default) or six-equation approach (see 2.3 for hydrodynamics description of both models), including four differential conservation equations of mass and energy of gas and liquid. Both approaches are operative, but the latter still requires further validation, particularly for the analysis of core reflooding.

On the other hand, the ICARE module handles the simulation of all in-vessel phenomena from the heat transfers in intact geometry to the core degradation and material relocation. In this study, the following heat transfers in intact geometry are concerned: conduction and gap heat transfers inside rods, and wall to fluid heat exchanges (convection, film boiling and radiative heat transfer to droplets and steam).

Radiative heat exchanges between rods and the surrounding shroud are neglected: the temperature difference between the shroud and the rods is too low to generate a significant heat flux compared to the flux extracted by the fluid.

Modelling reflooding in ASTEC V2.1

The current reflooding model is only devoted to situations where water enters at the bottom of the core.

In terms of quench front detection and heat transfer description, the model is valid regardless the core damage state and the selected hydrodynamic model. However, in terms of hydrodynamics, the momentum equation takes a different form depending on these, which in turn affects the calculation of the phase velocities and hence, the power extracted from the core. In the following, the article presents the main features of the reflooding model putting the special focus on intact geometries and the five and six-equation schemes. The modelling features presented hereafter are operational starting from the release of ASTEC V2.1.1.6.

Quench front detection

For the calculation of the quench front position, ASTEC V2.1.1.6 follows a similar approach as the one followed in ICARE-CATHARE code [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF]. The code postulates that the wall temperature profile upstream the quench front is flat due to the strong cooling of the wall by the quench front, down to a temperature close to the saturation temperature (the value Tsat+5 was selected). Besides, it considers that the temperature profile downstream the quench front is flat with values higher than minimum film boiling temperature TMFB. Finally, the temperature profile along the TBR is approximated by a step, where the lower bound is the temperature of the first mesh below the quench front (index IQF-1) and the upper bound is the temperature of the first mesh above the quench front (index IQF+1). The quench front mesh (IQF) is defined as the first mesh (starting from the bottom of the wall) where the wall temperature is higher than Tsat+5, providing that the temperature of the mesh IQF+1 is higher than TMFB and the void fraction in the mesh IQF is lower than 0.9999. The elevation ZQF of the quench front in this mesh is computed as a function of the current temperature in the mesh IQF, ensuring a regular progression of the quench front in the mesh.

Heat transfer along the TBR

As already pointed out in the introduction, ASTEC has two different models to describe the heat transfer along the TBR, these being the DROPLET and the EXPCHF models (the latter formerly known as DRACCAR model).

The DROPLET model postulates the existence of a large heat transfer coefficient HD (by default 150 W/m 2 /K) towards the interface in all meshes over a distance ZD (by default 0.8 m) above the quench front elevation ZQF. This heat transfer is applied to evaporate the water in a given mesh, providing that its void fraction is lower than αD (by default 0.999), the heat transfer being given by Eq. 1. Therein, Tliq is the temperature of the liquid in the current mesh and Twall is the temperature of the wall in contact with the liquid. The variables HD, ZD and αD are user-dependent and have a considerable impact on model predictions. The default values were obtained following an internal study involving a sensitivity analysis on several PERICLES experiments with the CESAR five-equation scheme. The default values of the DROPLET model are the ones minimizing the relative error between the computed and the experimental evolution of the quench front progression and cladding temperature at selected elevations.

! "#$ %$&' (" ) = + % × -. / 00 -. 0 1 2 Eq. 1

As for the EXPCHF model, the heat flux is described by a decreasing exponential function (see Eq. 2), following the formulation proposed by [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF] and later on adopted by [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF]. Therein, CHF represents the critical heat flux, z the absolute elevation of the section where the heat transfer is evaluated and L1 the characteristic length over which this heat transfer is applied. Alike [START_REF] Glantz | DRACCAR: A multiphysics code for computational analysis of multi-rod ballooning, coolability and fuel relocation during LOCA transients. Part One: General modeling description[END_REF], the CHF is calculated using Zuber's correlation for pool boiling with the correction on liquid subcooling due to Ivey and Morris [START_REF] Ivey | On the relevance of the vapour-liquid exchange mechanism for subcooled boiling heat transfer at high pressure[END_REF]. This is in contrast to [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF], who calculated the CHF using the Groenenveld look-up table and additionally included a dependency with the liquid velocity and the void fraction near the quench front.

The later approach could be implemented in the future.

! "#$ (3'4567 ) = + × 8 9:74; < Eq. 2
Such modelling takes advantage of previous experiments performed by Ishii and co-workers (Obot and Ishii, 1988), who observed that the L1 associated to each flow regime (smooth, rough wavy, agitated, dispersed) depended on the capillary number according to Eq. 3. The liquid velocity is calculated at the junction between the meshes IQF and IQF-1. However, such observations did not identify an expression of L1 for the TBR, which egged on [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF] and then [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF] to propose correlations based on the capillary number with the form given by Eq. 4. In both studies, the constant took two different values depending whether Re>2000, the value in each regime being different in both studies (an order of magnitude higher in [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF]). Due to the similarity between ASTEC and ICARE-CATHARE, this work adopts the values described in [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF]. Expression given by [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF], with a C1 coefficient an order of magnitude lower, have also been tested in ASTEC but the length of the TBR zone was then very short and the reflooding (if any) was too slow in all PERICLES tests. This can be explained by two remarks.

First, the values applied in [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF] are deduced from steady-state experiments in which the TBR may be shorter than in reflooding conditions. Second, the model of [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF] accurately describes all regimes downstream the quench front (TBR but also inverted annular flow and dispersed flow in which the heat transfer can also be efficient) while only one specific regime is described in ASTEC, covering a larger region.

> = ? @ A Eq. 3 B = > C.E Eq. 4
At this point, it is worth mentioning that the liquid velocity used to calculate the capillary number contains a correction factor depending on the selected hydrodynamic model. This is done to keep L1 within the physically sounded value of few centimetres. This way, 1 -@ is used in the five-equation scheme to correct the high liquid velocities associated with high void fractions i.e. residual presence of the liquid phase. Conversely, @ is used in the six-equation scheme without correction since the latter can deal with the presence of residual phases. Therefore, the introduction of such a correction factor allows calculating a length L1 that stays within physical bounds (few centimetres) regardless of the selected hydrodynamic model.

A sketch of the heat flux profiles associated to the EXPCHF and DROPLET models is depicted in Figure 1. Therein, the real flow patterns corresponding to a situation of a low-quality CHF (e.g. high flow rates, high pressures) and a high-quality CHF (e.g. low flow rates, high electrical powers) have been also

included. The solid red line is the calculated one during the simulation whereas the dashed one is the hypothetical one (if the restrictions of αD and L1 did not apply). As it is observed, the TBR length calculated by the two models is not necessarily the same. Indeed, in the DROPLET model this length is limited by ZD or, more likely, by the unavailability of water in the mesh αD, whereas in the EXPCHF the length is limited by the characteristic length L1. In the EXPCHF model, αD can also be limiting.

Heat transfer correlations in other flow regimes

As far as the rest of boiling heat transfer mechanisms are concerned, nucleate boiling below the quench front is evaluated with the Thom correlation [START_REF] Thom | Boiling in Sub-Cooled Water during Flow up Heated Tubes or Annuli[END_REF], whereas film boiling above the TBR is evaluated with the Berenson correlation [START_REF] Berenson | Film boiling heat transfer from an horizontal surface[END_REF].

Hydrodynamic aspects

The summary of the different terms considered in the momentum conservation equation for the five and six-equation hydrodynamic models without the presence of porous media were already introduced in (Gómez-García-Toraño and Laborde, 2019). For the sake of understanding, this section gives a short reminder on both numerical schemes implemented in CESAR.

For the five-equation scheme, one momentum balance equation on the mean phase velocity is solved for all junctions. Assuming that the flow along the test section is one-dimensional and that the cross section is uniform, the average momentum balance equation for the junction J along the flow direction can be written according to Eq. 5. Therein, the first term represents the axial pressure gradient, G and G the regular and singular friction terms, G the advection term, G the hydrostatic pressure gradient and H G the inertial term computed with the averaged fluid velocity.

IJ I) K G K G K G K G K H G = 0 Eq. 5
Additionally, the five-equation scheme introduces an algebraic equation for the calculation of the drift velocity between the gas and liquid phases (MN). For vertical junction, the drift velocity can be written according to Eq. 6. Therein, is a smoothing function depending on the mean void fraction, whereas FNR is the non-stratification friction term which depends on the geometry, flow configuration and other thermodynamic parameters [START_REF] Wallis | One-Dimensional Two Phase Flow[END_REF]. Such expression is evaluated if the presence of both phases is not residual (0.05 < α < 0.95), and it is set to zero (equity of gas and liquid velocities) in the rest of cases.

MN = N -N = O P - Q$
Eq. 6

The five-equation scheme is suitable for quasi-steady developed flows in pipes without strong phase changes and with 0 < < 0.7. However, for accelerating flows (where the advective term is significant), a six-equation scheme involving two sets of momentum conservation equations in each phase have better capabilities [START_REF] Bestion | System Code Models and Capabilities[END_REF]. Considering a one-dimensional flow, a uniform test section and neglecting the capillary forces (pressure of the liquid and the gas are similar and equal to P), the two momentum balance equations for the liquid and gas phases projected along the flow direction at the junction J can be written as follows:

IJ I) K K K K - K H = 0 Eq. 7 IJ I) K K K K K K H = 0 Eq. 8
It is observed that the interfacial friction between the liquid and gas phases is now explicitly considered. For a vertical non-stratified flow, the expression is given by Eq. 9. Therein, represents the interfacial friction coefficient, integrating the effects of geometry and flow configuration.

It can be obtained using a number of correlations, as explained in [START_REF] Wallis | One-Dimensional Two Phase Flow[END_REF].

= • • • MN • |MN| Eq. 9
In terms of numerical stability, the two schemes have similar behaviour while the five-equation scheme converges slightly faster in general. Numerical issues are more related to the reflooding models than to the hydrodynamic model: in the DROPLET model, the quench front progression from mesh to mesh is associated with large heat flux peaks because of different fluid conditions from one mesh to another.

This issue is less pronounced in the EXPCHF model in which the accurate quench front location is used to distribute heat fluxes on a continuous way upstream and downstream the quench front.

The PERICLES Experimental facility

Experimental setup

The PERICLES experimental program at CEA Grenoble (France) mainly aimed at improving the understanding on core thermal-hydraulics during the reflooding phase of a PWR [START_REF] Digonnet | PERICLES Programme : boil-up, boil-off and reflooding high pressure experiments in a PWR assembly[END_REF][START_REF] Housiadas | Two-dimensional quench front progression in a multiassembly rod bundle[END_REF]. The test section consists of an insulated stainless-steel shroud, which entails 368 electrically heated fuel rod simulators (FRSs) and 25 stainless steel guide tubes in a 17x17

geometry. Each FRS is electrically heated over 3656 mm and consists of three helical nichrome V wires, which are embedded into a boron nitride matrix. Stainless steel claddings surround the FRSs. The wire density is varied along the FRSs in order to obtain a PWR-representative cosine shape power distribution around the centre with an axial peaking factor of 1.6. A sketch of the test section is shown in Figure 2. In all reflooding tests, the electrical power is supplied to the FRSs, the injection being launched when maximum temperatures reach a set point and stopped when the core is refilled. Influence of pressure, heating power, injection mass flow rate and initial temperature were studied. Several tests with a wide range of pressures, mass flow rates, electrical powers and initial temperatures have been used to validate the code. The list of tests can be found in Table 1. The initial temperature at the hottest part of the bundle is 600°C for all tests. Mass flow rates are ranging between 1-5 g/cm 2 /s, which corresponds to nominal injection velocities of 10-50 mm/s (35-180 m/h). It is worth pointing out that the High-and Low-Pressure Injection Systems (HPIS, LPIS) of a French PWR 900 MWe generally consist of three trains, each of them giving nominal flow rates of 50 kg/s (HPIS) and 100 kg/s (LPIS). These correspond to injection velocities of about 9-18 mm/s (equivalent to 30-60 m/h), which are equivalent to those used in PERICLES experiments.

The assembly temperature field is mapped through thermocouples (THCs) located at different elevations along some FRSs. Pressure along the longitudinal axis can be measured through differential pressure transducers. The liquid mass flow rate injected and the steam mass flow rate exiting the bundle can be measured by two different flow meters. The quench front is assumed to reach a given axial location when the temperature of the nearest thermocouple falls below Tsat+5. No significant uncertainties can be attributed to the measurements, as pointed out in [START_REF] Digonnet | Essais de renoyage d'un assemblage complet de réacteur à eau pressurisée -Essais bouillotte et de dénoyage Essais transitoires mixtes[END_REF]. The ASTEC V2.1 model of the experimental bundle consists of one weighted rod, representing the 368 electrically heated FRSs, and one weighted hollow tube, representing the 25 guide tubes. Each FRS consists of a stainless-steel cladding and a fuel pellet, which encloses an inner and outer cylinder of boron nitride and a heated nichrome cylinder in between. The density of the nichrome cylinder varies along the active bundle length in order to reproduce the experimental power deposition. Both the FRSs and the guide tubes are surrounded by a shroud.

The bundle is axially discretized in 18 meshes giving an axial mesh height of about 200 mm, which corresponds to the current recommendation for reactor calculations with ASTEC. The axial discretization was the result of an iterative the study originally considered 30 axial meshes.

Then, axial refinement was progressively decreased down to an optimal value below which the prediction on the key figures of merit (quenching time and peak cladding temperature) started to be impaired. Concerning radial discretization, only one radial mesh has been considered when modelling the PERICLES test section. This radial discretization is justified because the tests were generally characterised by a one-dimensional quench front progression. However, some tests involving high mass flow rates and/or high heating powers and low pressures (e.g., p7) showed top flooding at the upper part of the bundle, which implies the existence of multi-dimensional effects. A finer radial discretization will be used in future validation studies on the 2-D PERICLES, where two-dimensional effects were more visible.

The temperature, pressure and mass flow rate vary according to the experimental values along the calculation. A pressure boundary condition is set at the bundle outlet, the values being varied according to the experimental data along the calculation.

Results using CESAR five-equations

Results using the DROPLET reflooding model

A set of predictions using the current reflooding model in CESAR five-equations is carried out for tests named p5 and p7, which differ only by the reflooding mass flow rate (see Table 1) and compared with experimental results. The quench front progression and the cladding temperature evolution at z=2550 mm (i.e. hottest bundle region) are shown in Figure 3.

The experimental temperature evolution hints the existence of four differentiated phases regardless the mass flow rate: heat-up phase, precooling, quenching and saturation. It is observed that the heat-up and pre-cooling phases are longer as the injected mass flow rate decreases. Indeed, a decrease of the mass flow rate reduces the quench front velocity and enhances the dispersion of the flow above the front, the last factor contributing to a shortening of the TBR.

As far as code predictions are concerned, CESAR five-equations provides a good representation of the experimental data when mass flow rates are higher than 3.6 g/cm 2 s e.g. p7 and p12 tests. However, the quenching times become increasingly underestimated as the mass flow rate decreases below that value, as shown in the prediction of the test The behaviour is congruent with the earlier quenching of the cladding. The behaviour of the cladding temperature at z=2500 mm during the first instants of the transient is of special interest. Within that time window, the quench front remains below z=1700 mm, this meaning that the DROPLET model (which is applied over a distance ZD = 800 mm) should not have a strong impact on the cladding temperature. Nevertheless, cladding temperatures are underestimated during that time gap of the heat-up and the precooling phases. This shows that the code overestimates the water availability far above the quench front for low reflooding mass flow rates. This point will be further discussed in section 5 when comparing the five and six-equation schemes.

It is noticed that ASTEC results predict bottom flooding at all bundle elevations whatever the mass flow rate. This is in contrast to the experimental results, where bottom and top flooding may take place, especially for high mass flow rates (p7). Indeed, a higher mass flow rate enhances 2-D effects in the bundle, leading to a faster progression of the water near the wall. Eventually, part of this water falls down back to the upper core regions. This top flooding mechanism, which is relevant as the mass flow rate increases, cannot be captured by the geometrical ASTEC model. This could be already anticipated, since only one radial mesh has been used to describe the core fluid region in the current ASTEC model (section 3.2). The modelling of top flooding is out of scope for this article, but future improvements in this area are foreseen.

Comparison of DROPLET vs. EXPCHF reflooding models

The previous section has shown that the DROPLET model predicts faster refilling times than the experimental ones as the injected flow rate decreases because of the stronger precooling associated with the choice of the reflooding model and the hydrodynamic model. Despite this deviation has a lot to do with the 5-equation model, it is still interesting to study the individual influence of the DROPLET and EXPCHF models. The comparison of the experimental and computational temporal evolution of the quench front and cladding temperature, steam generation rate and water flow rate at z=1800 mm is depicted in Figure 4 for the test p5. That elevation is taken to better illustrate the differences between the two reflooding models.

The choice of the model significantly affects the quench front progression and the cladding temperature evolution. Indeed, the DROPLET model manages to evacuate the heat stored in the fuel rod, which is the reason why the quench front reaches the top of the bundle. In contrast, the EXPCHF model cannot evacuate enough heat in the vicinity of the quench front, this leading to temperature escalation and increasing difficulty for quench front progression. The figure also shows that code predictions slightly underestimate the cladding temperature during the first instants of the transient, this becoming more apparent for higher cladding elevations.

The contrast of the behaviour can be explained in terms of the length over which each model is applied.

To explain this, the attention is focused on the instant t A . At that time, the quench front is located at z=1200 mm for both simulations, but the behaviour at z=1800 mm is considerably different. Indeed, the that section. This happens because the mesh located at z=1800 mm has enough water and because it lies within the ZD length imposed by the DROPLET model. Hence, the HD heat transfer coefficient is applied. Conversely, EXPCHF calculates a weaker evaporation, since a film boiling regime is detected by the code. This is consistent with the higher liquid flow rates at z=1800 mm.

Comparison CESAR five vs. six equations

The previous section has shown that the five-equation scheme only provides reasonable results when the DROPLET reflooding model is used. The aim of this section is to determine if the six-equation scheme and the DROPLET or EXPCHF reflooding models provide a better representation of the experimental results. To discuss this appropriateness, the experimental and predicted quench front progression, cladding temperature at z=2500 mm and integral steam mass at the outlet are represented in Figure 5 for the test p5. Black scattered dots and dashed lines represent the experimental data; solid lines are associated to predictions with the five-equation scheme, whereas dotted lines with the sixequation scheme; the blue colour is associated with the DROPLET model whereas the orange one with EXPCHF.

It is observed that the EXPCHF reflooding model together with the six-equation scheme provides a good agreement with the experimental results, unlike the same model being used with the five-equation scheme. In particular, the quench front can now rise to the top of the bundle because of the stronger precooling predicted by the six-equation scheme near the quench front, which is reflected in the cladding temperature evolution. For identical reasons, predictions using the DROPLET model and the sixequation scheme are worse than those obtained when using the five-equation scheme. Moreover, for a given numerical scheme, DROPLET predicts a faster refilling than EXPCHF because of the longer extension over which this model is calculated together with its high heat transfer coefficient HD.

Cladding temperatures at z=2500 mm are underestimated at the beginning of the transient for either scheme or reflooding model, even if the quench front is well below that location. This indicates that, for the test p5 involving a low mass flow rate, either scheme predicts a greater water availability well above the quench front in comparison to the experiment. The modelling of the small-droplet disperse flow regime will require further improvements for either numerical scheme. In turn, this stronger precooling leads to an acceleration of the quench front at the second half of the bundle for almost all cases.

An analysis of the integrated steam mass along the experiments shows that almost three out of four simulations can predict the experimental steam mass at the end of the experiment. For these scenarios, the bundle can be cooled down, the initial energy being extracted by water injection. It is noticed that the use of the EXPCHF model together with CESAR five-equations is significantly lower than the experiment and the rest of simulations, which indicates difficulties in terms of core cooling. The reader can observe that the slopes in the curves of the steam mass relate to the quench front velocity. In particular, the use of the DROPLET model grows quickly at the beginning of the transient, then becomes increasingly deviated from the experimental because of the underestimation of PCTs. To better understand the differences between the five and six-equation schemes, selected axial profiles corresponding to predictions with the EXPCHF reflooding model are depicted in Figure 6 for the test p5 at t=t*. This instant is chosen because it is the last one before the quench front progression starts to differ in both simulations. Indeed, the quench front elevation is represented by a horizontal dashed blue line in all subplots of the figure, which indicates that its location is the same at t* for both simulations.

On one hand, the five-equation scheme leads to the prediction of significant liquid velocities above the quench front, since the drift velocity is set to zero for void fractions higher than 0.95 (see section 2.3).

Hence, the available liquid is instantly made available in all meshes above the quench front. Part of this water is transformed into steam (especially along the TBR), but nearly 40 % of the mass flow rate exits the bundle without contributing to the overall cooldown, as shown in axial water flow rate. On the other hand, the six-equation scheme leads to differentiated phase velocities regardless the void fraction because two equations on phase velocities are solved. The prediction of small liquid velocities, which is congruent with the lower void fraction, ensures that the liquid mass is kept near the quench front.

Therefore, the steam generation rate along the TBR is generally more intense in the six-equation scheme than that predicted by the five-equation one.

Figure 6: Comparison of CESAR five vs. six equations using the EXPCHF reflooding model on selected axial profiles at time=t* for the test p5. The quench front elevation at t* is represented by a single horizontal blue line for both simulations.

Calculations on further PERICLES reflooding experiments have been launched using a combination of the presented hydrodynamic model and heat transfer reflooding models. For each calculation, the time needed by the quench front to reach 2500 mm and the Peak Cladding Temperature is retrieved. Results can be found in Figure 7 and Figure 8. The elevation of 2500 mm is convenient because of the inexistence of top flooding for all experiments, in contrast to higher elevations. The thermal-hydraulic features of each experiment can be found in Table 1.

At first glance, neither the hydrodynamic model nor the reflooding model affect the quenching time at high pressures (p74, p75, p77, p80) or high mass flow rates (p8 and p13), predictions being close to the experimental data. Indeed, experiments involving high pressures are characterized by a relatively long inverted annular flow region above the quench front, followed by an agitated region of slugs and by a small-droplet dispersed flow region [START_REF] Ishii | Flow visualization study of inverted annular flow of post dryout heat transfer region[END_REF][START_REF] Obot | Two-phase flow regime transition criteria in post-dryout region based on flow visualization experiments[END_REF], which in turn lead to a strong heat removal from the wall. Therefore, an agreement between experimental data and code predictions exists because a significant precooling is predicted whatever the combination of the reflooding model or the hydrodynamic model. However, the agreement on the quenching time is not completely reflected in the Peak Cladding Temperature. This is related to the overestimation of the heat transfer in the dispersed flow region, as already shown in Figure 5. Future improvements on the description of the small-droplet dispersed flow region are foreseen.

The rest of tests are characterized by combinations of electrical power, system pressure and mass flow rates in such a way that the flow ends up being highly dispersed. In those cases, a precise estimation of the precooling is necessary for a correct prediction of the quench front progression. It is observed that DROPLET predicts a faster refilling than EXPCHF regardless of the hydrodynamic model because of the higher integrated heat flux along the TBR, which spans over a considerable distance ZD. As far as the hydrodynamic model is concerned, the six-equation scheme predicts a faster refilling than the fiveequation scheme because of the higher availability of water near the quench front, as shown previously in this chapter. For those tests, the agreement with the PCT is satisfactory. This is especially true for the six-equation scheme together with the EXPCHF model, which predicts a slighter precooling than calculations the DROPLET model with either numerical scheme. Yet, predicted precooling is stronger than the experimental one, as reflected in the slight underestimation of the PCTs. Finally, it is observed that the EXPCHF together with five equations clearly overestimates PCTs, which is a direct consequence of the longer quenching times.

Discussion

To support the decision on the most suitable model for either scheme, the relative error of the fuel rod quenching time at z=2500 mm (with respect to the experimental data) for each PERICLES test has been depicted in Figure 9. The red bands represent the relative error of T20%, which is the maximum deviation admissible for severe accident codes. This variable has been selected because of its importance for accident management. Indeed, it is essential to know whether a given injection can cool down the core and the time needed for such aim. Relative error on PCT is not shown, since all calculations lie within the limits of T20% regardless of the combination of heat transfer model and hydrodynamic scheme.

Results show that the most suitable modelling of the heat transfer along the TBR depends on the choice of the hydrodynamic scheme. In particular, if the five-equation scheme is used, the best choice is to only rely on the DROPLET model with reasonable values of heat transfer coefficient and TBR length.

However, if the six-equation scheme is used, the EXPCHF model alone can cool down the wall in the quench front vicinity. Even if predictions generally lie within the admissible error band, none of the two combinations provide a perfect match for all experiments. For example, five-equations plus DROPLET have difficulties to predict p5 and p9 tests, characterized by a high electrical power, whereas sixequations plus EXPCHF have difficulties to predict p25 and p27 tests, characterized by a combination of low mass flow rates and pressures. This is related to the difficulties experienced by integral severe accident codes to provide an accurate description of precooling when the dispersed film boiling is dominant, as it is the case for the considered thermal-hydraulic conditions.

The computational costs of both modelling options are shown in Figure 10. It can be generally stated that the use of six-equations plus the EXPCHF model does not significantly impair computational costs with respect to the five-equation scheme plus the DROPLET model. In fact, it sometimes requires a lower computational cost even if both modelling options predict an identical refilling time (see p7, p8 and p25 tests simulations in Figure 7). Hence, the computational cost should not be in this case the decisive factor to choose the most adequate modelling option. However, reactor calculations should be carried out to verify this point.

The two options seem appropriate to represent the core reflooding, but the EXPCHF model together with six-equation scheme brings more physical meaning to the system. The model also allows a reduction of the user effect and hence contributes to narrow uncertainty margins, since all its main parameters have been derived from experimental observations [START_REF] Ishii | Flow visualization study of inverted annular flow of post dryout heat transfer region[END_REF][START_REF] Obot | Two-phase flow regime transition criteria in post-dryout region based on flow visualization experiments[END_REF]. It is reminded that the DROPLET model considers fixed coefficients to describe the length and heat transfer associated with the TBR. This is in contradiction to existing experimental studies, which suggest that the magnitude of precooling depends on the local thermal-hydraulic conditions near the front [START_REF] Yadigaroglu | Modeling of reflooding[END_REF].

  with the calculation of L1 (m) • Ca: Capillary number (-) • : interfacial friction factor in the six-equation scheme (kg/m 4 ) • : smoothing function for drift calculation in the five-equation scheme (-) • : advection term (Pa/m) • : hydrostatic pressure gradient (Pa/m) • : regular friction term (Pa/m) • : singular friction term (Pa/m)

  • i: reflooding model (DROPLET or EXPCHF)• J: junction under consideration • L: liquid phase • W: wall • M: mixture between L and G • Sat: saturation

  Currently, all reactor safety studies performed use CESAR in five-equation mode and reflooding calculations are launched using the DROPLET model. Nevertheless, IRSN is pursuing an intense validation effort to possibly replace the DROPLET model by the EXPCHF model which is more physical-sounded and less user-dependent.

Figure 1 :

 1 Figure 1: Sketch of axial parietal heat fluxes calculated along the TBR by the (left) DROPLET and (right) EXPCHF models. Experimental observations for low and high-quality CHF can be found at the centre and are based on (Nelson and Unal, 1992).

Figure 2 :

 2 Figure 2: Sketch of the PERICLES test facility.

Figure 3 :

 3 Figure 3: Comparison of the predicted and experimental temporal evolutions of the quench front location and the cladding temperature at z=2500 mm. Predictions correspond to the DROPLET reflooding model using CESAR five-equations on tests different flow rate.

Figure 4 :

 4 Figure 4: Influence of DROPLET and EXPCHF heat transfer models on the temporal evolution of the quench front progression and the cladding temperature, steam generation rate and water flow rate at z=2500 mm for the test p5. Predictions obtained using CESAR five-equations.

Figure 5 :

 5 Figure 5: Comparison of the predicted and experimental evolution of the quench front progression, cladding temperature at z=2500 mm, total steam mass at the outlet for the test p5.

Figure 7 :

 7 Figure 7: Predicted times to quench the fuel rods up to 2500 mm using the DROPLET or EXPCHF reflooding models together with CESAR five or six-equations vs. experimental data.

Figure 8 :

 8 Figure 8: Predicted Peak Cladding Temperature along the transient using the DROPLET or EXPCHF reflooding models together with CESAR five or six-equations vs. experimental data.

Figure 9 :

 9 Figure 9: Relative error in the prediction of the quenching time of the fuel rods up to 2500 mm. Dotted red lines indicate a relative error of T20% with the experimental data.

Figure 10 :

 10 Figure 10: Computational times associated with the quenching of the rods up to 2500 mm for two most suitable modelling options of this study (DROPLET and CESAR five equations, EXPCHF and CESAR six-equations).

Table 1 :

 1 Selected Pericles experiments for the validation of the ASTEC V2.1 reflooding model. All tests have an initial temperature of 600 ºC, which is larger than TMFB (condition required to apply the reflooding model). Tests with non-constant electrical power have been excluded from this study.

	Name	Pressure	Electrical power	Flow rate	Subcooling
		(bar)	(W/cm 2 )	(g/cm 2 s)	(ºC)
	p1	3	3.35	3.6	60
	p5	3	3.35	2.5	60
	p7	3	3.35	5	60
	p8	3	3.35	8	60
	p9	2	3.35	3.6	60
	p11	4	3.35	3.6	60
	p12	3	3.35	3.6	60
	p13	3	3.35	13	60
	p25	3	4.20	5	60
	p27	3	4.20	3.6	60
	p712	20	2	1.7	4
	p74	20	2	3	4
	p75	20	2	3	4
	p77	40	2	3	7
	p80	10	1	3	1

Conclusions

Within this article, the reference reflooding model of ASTEC V2.1, devoted to bottom flooding, has been alternately used with the five and six-equation schemes to simulate experiments of intact core reflooding performed in the PERICLES test facility. While the five-equation scheme uses a differential equation on the mean velocity and algebraic equation to calculate the drift velocity, the later solves two differential equations on the phase velocities with implicit appearance of the interfacial friction. The reference reflooding model is based on two main assumptions:

• A quench front is calculated if the temperature of a solid heated component falls below Tsat+5 and if the temperature of its upper component neighbour exceeds TMFB, providing there is enough water in the mesh.

• The heat transfer along the vicinity of the quench front is the sum of two contributions. First, an exponentially decreasing function over an empirical length, this approach being taken from [START_REF] Chikhi | Reflooding model for quasi-intact rod configuration: Quench front tracking and heat transfer closure laws[END_REF][START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF]. Second, the existence of a large and uniform heat transfer coefficient over a longer distance, this representing the projection of droplets well above the front. They are named EXPCHF and DROPLET model, respectively.

Calculations using different combinations of hydrodynamic models and heat transfer reflooding models along the TBR have been performed. The main conclusion is that the choice of the reflooding model depends on the choice of the hydrodynamic model. If the five-equation scheme is selected, the DROPLET model is the most appropriate option. However, if the six-equation scheme is used, the EXPCHF model is best, since it benefits from the larger amount of liquid available near the front.

Computational costs are not significantly impaired by using either option, but reactor calculations should be carried out to give more validity to this statement.

In turn, the present article recommends the using the six-equation approach together with the EXPCHF model for accidental sequences involving reflooding, because of the more significant physical soundness of the system together with improved validation results on PERICLES experiments.

There are still developments that can be put in place to improve the modelling in both situations.

Concerning the modelling in six-equations, improvements in the description of the Critical Heat Flux could be developed and additional terms can be introduced in the L1 length as suggested by [START_REF] Nelson | A phenomenological model of the thermal hydraulics of convective boiling during the quenching of hot rod bundles Part I: Thermal hydraulic model[END_REF]. With regards to the five-equation scheme, it makes sense to replace the fixed TBR length (ZD) by a dynamic length that depends on the capillary number and a multiplicative factor, to avoid excessive precooling far above the front.


K M D Motchon
email: koffi.motchon@univ-reims.fr

Lala H Rajaoarisoa
email: lala.rajaoarisoa@imt-lille-douai.fr

Midzodzi Komi
email: midzodzi.pekpe@univ-lille.fr

Etienne Pekpe
email: lucien.etienne@imt-lille-douai.fr

Lucien

K M Pekpe

L Etienne

S Lecoeuche
email: stephane.lecoeuche@imt-lille-douai.fr

An Algebraic Approach for Discrete Dynamic Reconstruction for Switched Bilinear Systems

Keywords: Bilinear system, switched system, regularly persistent input, mode compatibility

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

INTRODUCTION

Hybrid systems [START_REF] Goebel | Hybrid dynamical systems: modeling, stability, and robustness[END_REF] constitute a mathematical tool useful for modelling processes whose behaviour involves both continuous and discrete dynamics. Switched systems represent an important subclass of these systems. The continuous dynamic of a switched system is represented by a collection of differential (or difference) equations called modes, while its discrete dynamic is described by a piecewise constant function orchestrating the switching between the operating modes. The piecewise constant function called the switching signal indicates at each time the active mode of the system. This information is crucial for instance in stability analysis of switched systems [START_REF] Liberzon | Switching in systems and control[END_REF][START_REF] Lin | Stability and stabilizability of switched linear systems: a short survey of recent results[END_REF] and for the synthesis of asynchronous switching controllers [START_REF] Yuan | A novel lyapunov function for a non-weighted L 2 gain of asynchronously switched linear systems[END_REF][START_REF] Etienne | Stability analysis for switched uncertain nonlinear systems with dwell time and delay in the active mode detection[END_REF].

A growing amount attention has been paid in the last decades to the reconstruction of the switching signal from the input-output measurements. Studies on observability of switched systems, that is, the determination of conditions guaranteeing the uniqueness of the reconstruction of the switched signal and the initial continuous state from the input-output data are addressed in [START_REF] Babaali | Observability of switched linear systems in continuous time[END_REF]; [START_REF] Gómez-Gutiérrez | On the observability of continuous-time switched linear systems under partially unknown inputs[END_REF]; [START_REF] Lou | Necessary and sufficient conditions for distinguishability of linear control systems[END_REF]. In [START_REF] Babaali | Observability of switched linear systems in continuous time[END_REF] existence problem of a control input for the identification of the switching signal, refered as observability for some input, is analysed and algebraic conditions based on the modes Markov pa-rameters are provided for disturbance-free continuous-time LSSs (linear switched systems). In [START_REF] Baglietto | Active mode observability of switching linear systems[END_REF], the design of control inputs allowing the estimation of the switching signal is examined for discrete-time LSSs subject to bounded measurement noises. In [START_REF] Lou | Necessary and sufficient conditions for distinguishability of linear control systems[END_REF], for disturbance-free continuous-time LSSs, rank conditions are established to characterize the observability of the switching signal for all inputs. In [START_REF] Gómez-Gutiérrez | On the observability of continuous-time switched linear systems under partially unknown inputs[END_REF], for the class of continuous-time LSSs with disturbances, a geometric characterization of observability for some or all control inputs is provided in terms of maximal invariant space is provided. All these observability studies are based on the characterization of the distinguishability property of the operating modes. This property informs on the ability of the modes to generate different output signals for some or all control inputs.

Regarding the concrete problem of reconstruction of the switching signal from the input-output measurements, various algorithms have been proposed in the literature in the last few years. For instance, for discrete-time LSSs, modelbased diagnostic methods are used in [START_REF] Halimi | A tutorial on model-based modes detection and discernibility for switched affine discrete-time systems[END_REF] and [START_REF] Hakem | Switching time estimation and active mode recognition using a data projection method[END_REF]. The approach consists in designing a residual depending on the input-output data that have particular properties for the active mode and that are sensitive to the change of mode. A similar approach is also considered in [START_REF] Fliess | Realtime estimation of the switching signal for perturbed switched linear systems[END_REF] for continuous-time LSSs. To avoid derivative measurements in the continuoustime case, distribution framework is used in [START_REF] Tian | Switching time estimation for linear switched systems: an algebraic approach[END_REF] and [START_REF] Mincarelli | Active mode and switching time estimation for switched linear systems[END_REF]. In presence of distur-bances, reduced-order observers approaches are proposed in [START_REF] Van Gorp | Hybrid sliding mode observer for switched linear systems with unknown inputs[END_REF] and [START_REF] Yang | Active mode identification and continuous state estimation for switched linear systems with unknown inputs and slow switching signal[END_REF] for the reconstruction of the switching signal and the continuous state under suitable assumptions on the dimension and the detectability of the modes. Optimization-based methods have been also proposed for the identification of both the discrete and the continuous dynamics of LSSs, for more details, we suggest the readers to see the survey paper [START_REF] Garulli | A survey on switched and piecewise affine system identification[END_REF].

To the best of our knowledge, while several algorithms are proposed in the literature for estimating the switching signal of LSSs, development of switching signal reconstruction methods for switched systems with nonlinear continuous dynamics stills an open problem. In this context, the switching signal estimation problem is addressed in the present paper for continuous-time switched bilinear systems that are used for modelling numerous nonlinear systems [START_REF] Shu | Stabilizing controller design for state-based switching bilinear systems[END_REF], including thermal dynamics in buildings [START_REF] Chasparis | Supervisory system identification for bilinear systems with application to thermal dynamics in buildings[END_REF] for instance. Regularly persistent control input data [START_REF] Bornard | Regularly persistent observers for bilinear systems[END_REF] are considered. The proposed method is based on a compatibility test between the input-output measurements collected on the system and the dynamical behaviour of the modes. For this analysis, we use the distance from the output of the system to the allowable output signals of the modes associated to the control input. Indeed, we prove that when using an output signal generated by a regularly persistent control input, in disturbancefree case this distance is zero only for modes that are compatible with the measurement data. An analysis of the distinguishability property of the operating modes is also made in the paper. This property guarantees the uniqueness of the switching signal estimation from the input-output data. On this point, we mainly focus on the distinguishability of the modes w.r.t. the control input used for the estimation of the switching signal. For regularly persistent inputs, we specify a numerically verifiable condition for modes distinguishability. This condition provides a characterization of regularly persistent inputs guaranteeing the unique determination of the switching signal. It also generalizes the controlled distinguishability condition proposed recently in [START_REF] Motchon | Robust discerning controls for the operating modes of linear switched systems subject to bounded unknown inputs[END_REF] for linear operating modes.

Notice that a distinguishability condition for bilinear modes have been also established recently in [START_REF] Motchon | Necessary and sufficient condition for controlled distinguishability of continuous-time bilinear systems[END_REF]. Operating modes relative degrees are used to establish this condition. However as for controlled distinguishability studies of [START_REF] Babaali | Observability of switched linear systems in continuous time[END_REF] the condition of [START_REF] Motchon | Necessary and sufficient condition for controlled distinguishability of continuous-time bilinear systems[END_REF] only allow to verify whether there exists a control input making distinguishable the operating modes.

The paper is organized as follows. First, the switching signal estimation problem for bilinear switched systems is formulated in Section 2. The formal definition of inputoutput data compatibility with the dynamical behaviour of the modes as well as the concept of modes distinguishability w.r.t. an input are also presented in this section. Results related to compatibility tests are provided in Section 3.1 and the condition for distinguishability of modes w.r.t. persistent inputs is presented in Section 3.2. The proposed algorithm for estimating the switching signal which is based on the compatibility tests is provided in Section 3.3. The effectiveness of this algorithm is confirmed in Section 4 through a numerical example.

Notation: Throughout the rest of the paper, ⊗ denotes the Kronecker product symbol. The set

L 2 [τ 1 ; τ 2 ], R l consists of all functions ϕ : [τ 1 ; τ 2 ] → R l such that ϕ L2 := τ2 τ1 ϕ (τ ) 2 2 dτ < ∞ where • 2 is the Eu- clidean norm. The L 2 -distance from a function ψ ∈ L 2 [τ 1 ; τ 2 ], R l to a set S ⊂ L 2 [τ 1 ; τ 2 ], R l is denoted d L2 (ψ, S ) := inf ϕ∈S ψ -ϕ L2 . For a matrix M ∈ R l×l , M denotes its transpose matrix. If M 0 (i.e. M is symmetric positive definite), • M is the weighted norm defined by w M = √ w M w.

PROBLEM FORMULATION AND PRELIMINARIES

Problem formulation

Let's consider the switched bilinear system described on [t 0 ; ∞) by the following equations:

ẋ = A σ x + N σ (u ⊗ x) + B σ u, y = C σ x, (1)
where x ∈ R n is the state vector, u ∈ R is the control input vector, y ∈ R m is the output vector and σ : [t 0 ; ∞) → P ⊂ N is a piecewise constant function continuous from the right, modelling the switching signal of the system i.e. the signal indicating at each time the active mode of the system. For every p ∈ P,

A p ∈ R n×n , B p ∈ R n× , C p ∈ R m×n and N p = N (1) p • • • N ( ) p ∈ R n×n (with N (r) p ∈ R n×n , r = 1, 2, . . . ,
) are known matrices defining the continuous dynamic of the system. The switching times of the system are elements of the strictly increasing sequence {t k } k∈N ⊂ [t 0 ; ∞). There are no jumps in the state vector evolution. It is also assumed that σ admits a minimum dwell-time, that is: Assumption 1. There exists τ d > 0 such that for every

k ∈ N, t k+1 -t k ≥ τ d .
Given a final observation time T > t 0 + τ d , the problem discussed in this paper is concerned with the reconstruction of the switching signal on [t 0 ; T ] using the control input signal and the output signal collected on the system during the observation time [t 0 ; T ].

Compatibility between input-output data and modes

The operating modes of system (1) are the following bilinear systems:

∀p ∈ P, ẋp = A p x p +N p (u ⊗ x)+B p u, y p = C p x p .
(2) Given τ ≥ t 0 , the output vector generated by mode p ∈ P at a time instant t ≥ τ , when starting at τ from an arbitrary state x o p ∈ R n is denoted by y p t, τ , x o p , u and can be expressed explicitly as follows:

y p t, τ , x o p , u = C p Φ u p,τ (t) x o p + Ψ u p,τ (t) , (3)
where Φ u p,τ (t) ∈ R n×n defined by

Φ u p,τ (t) = exp t τ {A p + N p [u (s) ⊗ I n ]} ds , (4)
is the transition matrix and Ψ u p,τ (t) ∈ R m is given by:

Ψ u p,τ (t) = C p t τ Φ u p,s (t) B p u (s) ds. (5) Remark 1. Notice that the transition matrix Φ u p,τ (t) veri- fies Φu p,τ (t) = (A p + N p [u (t) ⊗ I n ]) Φ u p,τ (t) and Φ u p,τ (τ ) = I n .
The allowable output signals that can be generated by mode p ∈ P during an observation time [τ ; τ ] ⊆ [t 0 ; T ] and when the input u is applied to the switched bilinear system is therefore [START_REF] Motchon | Robust discerning controls for the operating modes of linear switched systems subject to bounded unknown inputs[END_REF]. The existence of right-pseudo inverse of modes observability Gramian is a key tool to establish this result. In the bilinear case, the observability Gramian of each mode p on [τ ; τ ], given by

Y u p ([τ ; τ ]) = y p •, τ , x o p , u | [τ ;τ ] : x o p ∈ R n , (6) where y p •, τ , x o p , u | [τ ;τ ] is the restriction of y p •, τ , x o p , u on [τ ; τ ]. Remark 2. In the linear case (N p = 0 n×n ), output spaces Y u p ([τ ; τ ]) are closed in L 2 ([τ ; τ ], R m ) , • L2 , see for in- stance Appendix A of
G u p (τ , τ ) = τ τ Φ u p,τ (s) C p C p Φ u p,τ (s) ds (7)
is nonsingular when the control input is universal [START_REF] Bornard | Regularly persistent observers for bilinear systems[END_REF] for mode p ∈ P on [τ ; τ ]. So using a reasoning similar to that in Appendix A of [START_REF] Motchon | Robust discerning controls for the operating modes of linear switched systems subject to bounded unknown inputs[END_REF], one can see that when the control input u is universal for the modes, the output spaces Y u p ([τ ; τ ]) are also closed in the bilinear case.

To identify the active mode, it will be important to verify the compatibility between the system's input-output data and the dynamical behaviour of the operating modes (or simply, the M-IO Compatibility). The formal definition of the concept of M-IO Compatibility is expressed as follows [START_REF] Cocquempot | Fault detection and isolation for hybrid systems using structured parity residuals[END_REF]: Definition 1. (M-IO Compatibility). The input-output data u and y collected on the switched system (1) are said to be compatible with mode p ∈ P during the observation time [τ ; τ ], and we write for simplicity p ⊕(u, y) on

[τ ; τ ], if y |[τ ;τ ] ∈ Y u p ([τ ; τ ]).
Compatibility tests allow to identify modes that can produce the output data used to infer the discrete dynamic of the switched system.

From distinguishability to active mode identification

A worst case where it is difficult to detect the active mode of the switched system using the input-output data is the situation where the data are compatible with at least two operating modes, that is:

∃ (p, q) ∈ P × P, y |[τ ;τ ] ∈ Y u p ([τ ; τ ]) ∩ Y u q ([τ ; τ ]).
To better understand this problem, an analysis of controlled distinguishability of modes [START_REF] Babaali | Observability of switched linear systems in continuous time[END_REF][START_REF] Motchon | Robust discerning controls for the operating modes of linear switched systems subject to bounded unknown inputs[END_REF], recalled in Definition 2 is required. Definition 2. (Controlled distinguishability). Two modes p ∈ P and q ∈ P of system (1) are said to be con-

trolled distinguishable w.r.t. u on [τ ; τ ] if Y u p ([τ ; τ ]) ∩ Y u q ([τ ; τ ]) = ∅. The input u is called in this case a discerning input of the modes on [τ ; τ ].
Characterization of controlled distinguishability is important to know a priori whether the active mode is uniquely identifiable or not from the input-output data. Augmented systems of modes [START_REF] Gómez-Gutiérrez | On the observability of continuous-time switched linear systems under partially unknown inputs[END_REF][START_REF] Motchon | Robust discerning controls for the operating modes of linear switched systems subject to bounded unknown inputs[END_REF] are main tools for addressing operating modes distinguishability problem1 . In the bilinear case which is the setting of the present paper, the augmented system associated to any pair (p, q) ∈ P × P of modes is the bilinear system defined by the following equations [START_REF] Motchon | Necessary and sufficient condition for controlled distinguishability of continuous-time bilinear systems[END_REF]: ẋp q = A p q x p q +N p q (u ⊗ x p q )+B p q u, y p q = C p q x p q (8) where x p q (t) = x p (t) x q (t)

∈ R 2n and y p q (t) ∈ R m are the state vector and the output vector,

A p q = diag (A p , A q ), B p q = B 1 B 2 , C p q = [C p -C q ] and N p q = N (1) p q • • • N ( ) p q with N (r) p q = diag N (r)
p , N (r) q , r = 1, 2, . . . , ; we recall that N (r) p and N (r) q are the rth block matrix of N p and N q , respectively. The output vector at the time instant t ≥ τ generated by the augmented system when the input u is applied to the augmented system, starting from the state x o p q at the time τ is denoted by y p q t, τ , x o p q , u . An explicit expression of this vector can be obtained similarly to (3)-( 5). The set of all output signals generated by the augmented system associated to modes p and q is then given by

Y u p q ([τ ; τ ]) = y p q •, τ , x o p q , u | [τ ;τ ] : x o p q ∈ R 2n . ( 9
) Remark 3. For every distinct modes p ∈ P and q ∈ P, Y u p q ([τ ; τ ]) is also closed in L 2 ([τ ; τ ], R m ) , • L2 when the observability Gramian G u p q (τ , τ ) of their augmented system on [τ ; τ ] is nonsingular. This matrix is explicitly given by G u p q (τ , τ ) = τ τ Φ u p q,τ (s) C p q C p q Φ u p q,τ (s) ds with Φ u p q,τ (s) = exp s τ {A p q + N p q [u (θ) ⊗ I 2n ]} dθ the transition matrix of the augmented system.

Considering the augmented system allows to recast the distinguishability analysis into output-zeroing problems [START_REF] Gómez-Gutiérrez | On the observability of continuous-time switched linear systems under partially unknown inputs[END_REF]. For controlled distinguishability, a formulation of the output-zeroing problem based on the admissible output space of the augmented system is provided by the following result (Motchon and Pekpe, 2019): Proposition 1. Two modes p ∈ P and q ∈ P of the switched system (1) are controlled distinguishable w.r.t.

u on [τ ; τ ] iff 0 [m] / ∈ Y u p q ([τ ; τ ]
), where 0 [m] denotes the zero output signal.

MAIN RESULTS

A criteria for testing M-IO Compatibility

For testing M-IO Compatibility, we propose the following geometric characterization: For simplicity in the presentation, we define ∆ u,y p (τ , τ

) := d L2 y, Y u p ([τ ; τ ]) . ( 11
)
Proposition 3 gives a numerical method for the computation of ∆ u,y p (τ , τ ), needed for applying the M-IO Compatibility result established in Proposition 2. Proposition 3. Assume that u is universal for mode

p ∈ P on [τ ; τ ]. Then ∆ u,y p (τ , τ ) = Ψ u p,τ -y 2 L2 -Υ u,y p,τ (τ ) 2 H u p (τ ,τ ) , (12) where H u p (τ , τ ) = G u p (τ , τ ) -1 , Ψ u p,
τ is given by ( 5) and Υ u,y p,τ (τ ) ∈ R n is defined by

Υ u,y p,τ (τ ) = τ τ Φ u p,τ (s) C p Ψ u p,τ (
η p (τ ) ν p (τ ) = 0 2n , (14b) where
Λ p = A p 0 n×n C p C p -A p , Π p = diag B p , -C p , (15)
and

Γ p = Γ (1) p • • • Γ ( ) p 0 2n×2nm , (16)
with

Γ k p = diag N (r) p , -N (r) p , r = 1, 2, . . . , . (17)
Proof. It follows from the definitions ( 6) and ( 11) that ∆ u,y p (τ , τ ) 2 = inf

x o p ∈R n F u,y p x o p where F u,y p x o p = y p •, τ , x o p , u -y 2 L2 . It is a simple exercise to verify that ∆ u,y p (τ , τ ) = F u,y p x o p with x o p = -H u p (τ , τ ) Υ u,y p,τ (
τ ), which is equivalent to (12). Now, we will prove that Ψ u p and Υ u,y p verify (13). For every t ≥ τ , Ψ u p,τ (t) defined by ( 5) can be rewritten as follows:

Ψ u p,τ (t) = C p Ψ u p (t) , (18)
where Ψ u p (t) = t τ Φ p,s (t) B p u (s) ds ∈ R n . By observing the explicit expression (3)-(5) of the output vector of the modes, one can see that Ψ u p is the solution of the differential equation Ψu

p (t) = A p Ψ u p (t) + N p u (t) ⊗ Ψ u p (t) + B p u (t) , (19a) Ψ u p (τ ) = 0 n . (19b) Furthermore, it is straightforward to verify that Υ u,y p,τ (t) = Φ u p,τ (t) Υ u,y p (t) , ∀t ≥ τ , (20)
where Υ u,y p (t) =

t τ Φ u p,s (t) C p C p Ψ u p (s) -y (s) ds with Φ u p,s (t) = exp - t s A p + [u (θ) ⊗ I n ] N p dθ .
Observing the expression of Υ u,y p (t) and Φ u p,s (t) as well as the explicit expression (3)-( 5) of the output vector of the modes, it is also a simple exercise to check that Υ u,y p is the solution of the differential equation Υu,y

p (t) = -A p Υ u,y p (t) + M p u (t) ⊗ Υ u,y p (t) + (21a) C p C p Ψ u p (t) -C p y (t) , Υ u,y p (τ ) = 0 n , (21b)
where

M p = -N (1) p N (2) p • • • N ( ) p
. Now, stacking up the differential equations ( 19a) and (21a) with initial conditions (19b) and (21b), and using the relation

N p u (t) ⊗ Ψ u p (t) M p u (t) ⊗ Υ u,y p (t) = Γ p u (t) y (t) ⊗ Ψ u p (t) Υ u,y p (t) , one gets that Ψ u p (t) Υ u,y p (t) is solution of (14a)- (14b)
and by uniqueness of the solution of this differential equation, we have Ψ u p (t) = η p (t) and Υ u,y p (t) = ν p (t). Therefore we deduce from ( 18) and ( 20) that (13) holds and this concludes the proof.

Characterization of modes distinguishability

To characterize controlled-distinguishability of modes, the following necessary and sufficient condition is proposed: Proposition 4. (Controlled distinguishability). Let u be an universal input of the augmented system associated to a pair (p, q) ∈ P × P of distinct modes i.e. an input such that G u p q (τ , τ ) is nonsingular. Then modes p and q are controlled distinguishable w.r.

t. u on [τ ; τ ] iff ∆ u,0 [m] p q (τ , τ ) := d L2 0 [m] , Y u p q ([τ ; τ ]) > 0. ( 22
)
Proof. Since u is universal on [τ ; τ ] for the augmented system associated to modes p and q, then 22) holds. Consequently, we deduce from Proposition 1 that condition ( 22) is necessary and sufficient for the controlled distinguishability of modes p and q w.r.t. the input u on [τ ; τ ].

Y u p q ([τ ; τ ]) is also closed in L 2 ([τ ; τ ], R m ) , • L2 ; see Remark 3. This implies that 0 [m] / ∈ Y u p q ([τ ; τ ]) iff condition (
Remark 5. An explicit expression of ∆ u,0 [m] p q (τ , τ ) can be obtained similarly to (12). With this explicit expression, Proposition 4 then generalizes to bilinear operating modes, Lemma 6 of [START_REF] Motchon | Robust discerning controls for the operating modes of linear switched systems subject to bounded unknown inputs[END_REF] that characterizes controlled distinguishability for linear operating modes. Furthermore, the value of ∆ u,0 [m] p q (τ , τ ) can be obtained using a numerical method as in Proposition 3.

Reconstruction of the switching signal

The notation t k , k ∈ N will be used for the estimate of t k . By convention t 0 = t 0 . The estimated switching times { t k } k∈N will be determined successively. For the tests of controlled distinguishability and M-IO Compatibility, we suppose that: Assumption 2. ∀p ∈ P, G u p t k , t k + τ d 0 and ∀ (p, q) ∈ P × P, p = q, G u p q t k , t k + τ d 0.

This hypothesis expresses a regularly persistence condition of the excitation u for the operating modes and the augmented systems [START_REF] Bornard | Regularly persistent observers for bilinear systems[END_REF]. Furthermore, for simplicity, we assume that u is a discerning control for the modes: Assumption 3. ∀ (p, q) ∈ P × P, p = q, modes p and q are controlled distingushable w.r.t u on [

t k ; t k + τ d ].
A procedure for designing regularly persistent discerning inputs may be proposed based on the result of Proposition 4.

Offline algorithm for the estimation of σ Inputs: t 0 , T , dwell-time τ d , data u |[t0;T ] and y |[t0;T ] , and small scalars ε 1 > 0 and ε 2 > 0.

Step 1. Initialisation: Fix k = 0 and t k = t 0 .

Step 2. Detection of the active modes p k : ∀p ∈ P, compute ∆ u,y p t k , t k + τ d using ( 12)-( 17). Identify the mode p k satisfying the relation ∆ u,y p k t k , t k + τ d ≤ ε 1 and go to Step 3.

Step 3. Estimation of the switching times t k+1 :

Fix

a k = t k and b k = T . While |a k -b k | > ε 2 , compute c k = (a k + b k ) /2. If ∆ u,y p k t k , c k > ε 1 then b k ←-c k Else, a k ←-c k . End if. τ k+1 ←-b k End while
Step 4. Reinitialisation: until all active modes are detected i.e. until ∆ u,y

p k t k , T ≤ ε 1 , k ←-k + 1
and t k ←-t k+1 , and go to Step 2.

Notice that since for all p ∈ P, ∆ u,y p t k , t k + τ d are obtained with numerical errors using ( 12)-( 17), in

Step 2, for the detection of the active mode, inequality ∆ u,y

p k t k , t k + τ d ≤ ε 1 is used instead of the relation ∆ u,y p k t k , t k + τ d = 0. In Step 3, t k+1 is computed as t k+1 = inf{t : t > t k , ∆ u,y p k t k , t > ε 1 } using a procedure similar to the bisection method. Condition ∆ u,y p k t k , c k > ε 1 means that p k is not active on the whole interval [ t k ; c k ].
The system has switched on the time interval [ t k ; c k ] and the switching time t k+1 is then located in this interval.

AN ILLUSTRATIVE EXAMPLE

Let's consider the switched bilinear system (1) with the three operating modes (P = {1, 2, 3}) described by the following matrices: A 1 = -0.19 0.4 -0.1763 -1.19

, A 2 = -0.25 0.5 -0.75 -1.5 , A 3 = -0.146 0.5 -0.17 -1.25 We choose τ d = 2s and ε 1 = ε 2 = 10 -3 . We have verified successively on intervals [ t k ; t k + τ d ) that Assumption 2 of regularly persistence of the input as well as Assumption 3 of distinguishability are satisfied. For page limitation, these numerical results are not presented. Estimated switching times and actives modes obtained with the proposed method are given in Table 1 (see columns 3 and 6). These results prove the efficiency of the proposed method.

, B 1 = B 2 = B 3 = 1.7 -1.9 1.8 -1 , C 1 = [1 0], C 2 = [1 1], C 3 = [0 1], N 1 = -0.
All the active modes are well estimated with a minimum delay of 0.00028 second and a maximum delay of 0.056 second (see column 4 of Table 1). The different compatibility tests that provide the active mode detection and switching times estimation results of Table 1 are presented in Table 2.

CONCLUSION

A method for detecting the active modes and estimating the switching times for continuous-time switched bilinear systems is proposed in disturbance-free situation. It is based on a compatibility test between the input-output data and the dynamical behaviour of the operating modes.

The method applies only for data generated by persistent inputs that distinguish the input-output behaviour of the modes. A characterization of these inputs is also provided in terms of a necessary and sufficient condition for modes controlled distinguishability. Based on the established distinguishability condition, future works will focus on the experiment design problem for the proposed switching signal estimation algorithm i.e. the synthesis of persistent discerning control inputs of the modes. An extension of the results to bilinear systems with disturbances or measurement noises is also envisaged.

Proposition 2 .

 2 (M-IO Compatibility). Let u be an universal input of mode p ∈ P on [τ ; τ ]. Then p ⊕ (u, y) on [τ ; τ ] ⇐⇒ d L2 y, Y u p ([τ ; τ ]) = 0. (10) Proof. The proof follows from the fact that a point is an element of a close set iff the distance from the point to the set is zero. The closure property of Y u p ([τ ; τ ]) is discussed in Remark 2. Remark 4. Notice that when Y u p ([τ ; τ ]) is not closed, d L2 y, Y u p ([τ ; τ ]) = 0 does not necessary implies that data u and y are compatible with mode p.

Fig. 1 .

 1 Fig. 1. Output signal y generated by u and σ

Table 1 .

 1 Numerical results for the estimation of switching times and active modes

	Iteration k	t k	t k	t k -t k	σ (t k )	σ t k
	0	0	0	0	2	2
	1	3.52	3.5213	0.0012708	1	1
	2	8.45	8.4528	0.0028076	3	3
	3	12.25	12.25	0.0002756	2	2
	4	15.2	15.256	0.055753	3	3
	5	17.5	17.504	0.0040873	1	1

Table 2 .

 2 Compatibility tests for the modes

	Iteration k	p = 1	∆ u,y p	t k , t k + τ d p = 2	p = 3
	0	137.51	6.0657 × 10 -8	80.949
	1	1.6514 × 10 -6		1.3401	3.5772
	3	1.9893		2.4465	1.8625 × 10 -7
	3	0.70282	0.00027665	4.821
	4	25.764		6.3962	3.8743 × 10 -11
	5	2.2373 × 10 -8		1.0863	11.585

The augmented system associated to two modes is a dynamical system whose state vector is the concatenation of the modes state vectors and its output is the difference of their output signals.