Nadia Kari
email: nadia.kari@univ-paris-est.fr

Chukwuka Ozuem

Divitha Seetharamdoo
email: divitha.seetharamdoo@ifsttar.fr

Jean-Marc Laheurte

Franc ¸ois Sarrazin

Analysis of a Chiral Helix Metamaterial Using Eigenmode Expansion Method and Characteristic Mode Theory

Metamaterials are commonly associated to antennas and other microwave devices due to their unique ability to manipulate electromagnetic waves. In general, the approach used to select metamaterials for a given application is based on a classical approach by extracting the scattering parameters to evaluate the material effective properties. In this paper, we address the problem of evaluation and analysis of a specific metamaterial, a chiral helix metamaterial through the use of two modal analysis Expansion Eigenmode Method and the Characteristic Mode Theory. An interaction of the chiral helix with Circular Polarized electromagnetic plane waves is studied. The modal net stored energy of the metamaterial is also calculated. A good agreement is obtained between the net stored energy calculated by the two methods for the chiral metamaterial. The effect of polarisation on the metamaterial is also highlighted by the modal analysis. These modal approaches applied to chiral metamaterials can be of interest for the design of circularly polarised metamaterial antennas.

I. INTRODUCTION

Nowadays, the study of electrically small antennas with high performances in terms of gain, bandwidth and radiation pattern is an appealing research topic. Metamaterials (MTMs) and more particularly their potential to improve antenna performances [START_REF] Xiong | Impedance bandwidth and gain improvement for microstrip antenna using metamaterials[END_REF] are also extensively studied at microwave frequencies. Generally, the methodology used to analyse MTMs while associating them to antennas, is based on their effective permittivity and permeability. The major disadvantage of this approach is that it does not give a full physical insight on the radiation properties and nearfield coupling. Modal methods such as Eigenmode Expansion Method (EEM) [START_REF] Ramm | Theoretical and practical aspects of singularity and eigenmode expansion methods[END_REF] and Characteristic Mode Analysis (CMA) [START_REF] Vogel | Characteristic mode analysis: Putting physics back into simulation[END_REF] can however provide a better physical insight on the radiating phenomena and based on the modal currents. Indeed, metamaterials can thus be classified into magnetic or electric structures in terms of stored energy, before their association to a given antenna. The main advantage of EEM and CMA is that they are independent of excitation and they are applicable to arbitrarily-shaped structures. The net and modal stored energies of MTMs can further be calculated thus providing insight on the antenna radiation efficiency. In this paper, our objective is to analyze the properties of chiral metamaterials which could be an interesting element for circularly polarized antennas with two modal approaches: EEM and CMA. To this end, after performing a modal analysis on the metamaterial, the net stored energy of a chiral metamaterial is calculated and compared for the two modal approaches. A single-turn helix metamaterial [START_REF] Semchenko | Stored and absorbed energy of fields in lossy chiral single-component metamaterials[END_REF] is considered. The specificity of this metamaterial is that it interacts only with waves of one of the two orthogonal circular polarizations (CPs). This type of MTM is used for wide types of application from antennas and microwave to optics.

II. STORED ENERGY OF HELIX METAMATERIAL:

COMPARISON BETWEEN EEM AND CMA

In this section, our aim is to calculate the net stored energy of a chiral metamaterial helix using the two modal analysis EEM and CMA. The design of the structure consists in a single turn metallic strip helix the dimensions of R = 7.75 mm, h = 12 mm, 2r = 1.55 mm and with a pitch angle of α = 13.65 • [START_REF] Semchenko | Stored and absorbed energy of fields in lossy chiral single-component metamaterials[END_REF] as depicted in Fig. 1 the evaluation of the impedance matrix, the commercial tool FEKO is used to design and to evaluate the impedance matrix of the structure. The frequency range varies from 2.6 GHz to 2.9 GHz. The Z matrix is then extracted and the stored energy is computed based on the evaluation of the eigenvalues of both of CMA and EEM. The net stored energy W m -W e is computed based on the net stored energy formulation defined in [START_REF] Capek | Optimal composition of modal currents for minimal quality factor q[END_REF]. The first two modes for each modal analysis method are calculated and presented in Fig. 1. It is assumed that two modes are sufficient to analyze electrically small structures [START_REF] Cabedo-Fabres | The theory of characteristic modes revisited: A contribution to the design of antennas for modern applications[END_REF]. For CMA, the characteristic currents are real vectors, whereas the eigencurrents are complex vectors for EEM. Therefore we sorted the eigencurrents according to the ratio between imaginary and real parts of the complex eigenvalues of EEM. The W m -W e is plotted against frequency. When this quantity W m -W e vanishes, it indicates the resonance. The helix metamaterial resonates at 2.82 GHz for CMA. The net stored energy of mode 1 is negative and then positive after the resonance, this indicates that the energy switches from electric to magnetic at resonance. For mode 2, Fig. 1 shows that the stored electric energy is greater than the magnetic one in the whole frequency band of interest. It should be noted that the stored energy presented by EEM are complex quantities, since the eigencurrents are complex vectors [START_REF] Ramm | Theoretical and practical aspects of singularity and eigenmode expansion methods[END_REF]. Fig. 1, shows the real part of the W m -W e (the imaginary parts are not shown). The W m -W e vanishes at 2.82 GHz, which indicates the resonance. The net stored energy of mode 1 goes from electric to magnetic energy at resonance. For mode 2 the net stored energy is mostly electrical. Current distributions, at 2.82 GHz for TCM, of mode 1 switches from capacitive to inductive modes and mode 2 is a purely capacitive mode as shown in Fig. 3 (b).

(b). For the design and

III. SCATTERING PARAMETERS AND MODAL WEIGHTING COEFFICIENTS (MWC)

The scattering parameters of a periodic array of spiral helix is computed by the numerical simulator CST Microwave Studio. Periodic boundary conditions are applied to simulate the periodic structure. Three configurations are investigated. In the first case, the electric field E is polarized along the axis of the helix and the magnetic field H are along the y direction. In the second case, the H is polarized along the axis of the helix and the E is polarized along the x direction and in the third case the E is along the helix axis and H along y direction as shown in Fig. 2. The reflection and the transmission coefficients are computed with CST. The results are depicted in Fig. 2. The resonance frequency is found around 2.82 GHz for all the configurations which agree with the results of the stored energy. It is important to note that around the resonance, the components of the reflected fields have the same magnitude with a 90 • shift between the components. This is to confirm the CP of the reflected field by the helix. In order to analyze the chirality and the interaction of the helix with different excitation, the MWCs are computed for the three configurations mentioned above. As illustrated in Fig. 3. mode 1 presents a coefficient peak at the resonance frequency for all the configurations except for mode 2 which is a purely reactive mode over all the frequency band. This is can be predicted 3 that the helix has a maximum interaction for case 2. It is important to note that the optimal helix has "maximum chirality" at the resonance and highly affected when it is excited with axial magnetic field and electric field polarized along the ends of the helix, this chiral material has different response and this agrees with the results that the spiral helix emit circularly polarized wave.

IV. CONCLUSION

The properties of chiral helix metamaterial is demonstrated using two modal analysis CMA and EEM. FEKO was used to evaluate the impedance matrix and a mathematical formulation was implemented to evaluate the stored energies. The analysis of the chiral metamaterial by modal approaches shows a strong dependence on the polarization of the fields. The stored energies were also calculated for these metamaterial. It should be noted that the stored energies are independent of excitation. Further studies will be presented where the excitation can be accounted for to determine the net stored energy for a given application. The comparison between EEM and CMA are in good agreement also for chiral metamaterials, not only in terms of modal significance but also currents distributions, far fields and stored energies.

The authors acknowledge funding from I-Site FUTURE in the framework of the Western project.

Fig. 1 .

 1 Fig. 1. (a) Net Stored Energy against frequency by CMA and EEM. (b) the chiral metamaterial structure

Fig. 2 .

 2 Fig. 2. Scattering parameters of all the configurations

Fig. 3 .

 3 Fig. 3. (a) Modal Weighting Coefficients for different configurations (b) Modal currents of the helix metamaterial by CMA


Gérald Camille Leroux

Gérald Le Mestre

Christophe Jego

Patrick Adde

Michel Jezequel

A 5-Gbps FPGA prototype of a (31,29)2 Reed-Solomon Turbo Decoder

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

I. INTRODUCTION

In telecommunications, forward error correction (FEC) is a system of error control that improves digital communication quality. Since the early 90's, channel coding families based on iterative decoding process are considered to be the most efficient coding schemes. With the invention of turbo codes [START_REF] Berrou | Near Shannon limit error correcting coding and decoding: Turbo Codes[END_REF] followed by the rediscovery of low-density parity-check (LDPC) codes [START_REF] Gallager | Low Density Parity Check Codes[END_REF], it is now possible to approach the fundamental limit of channel capacity. The turbo code family refers to two classes of codes: Convolutional Turbo Codes (CTCs) [START_REF] Berrou | Near Shannon limit error correcting coding and decoding: Turbo Codes[END_REF] and Turbo Product Codes (TPCs) [START_REF] Pyndiah | Near optimum decoding of product codes[END_REF]. The general concept of TPC is based on product codes that were constructed by the serial concatenation of two (or more) systematic linear block codes. In the last few years, many TPC decoder architectures have been designed [START_REF] Adde | Performance and complexity of block turbo decoder circuits[END_REF][START_REF] Cuevas Ordaz | New architecture for high data rate turbo decoding of product codes[END_REF]. Product codes using binary Bose-Chaudhuri-Hocquenghem (BCH) component codes have generally been chosen [START_REF] Pyndiah | Block Turbo Codes: Ten years later[END_REF]. The turbo product code concept was recently successfully extended to RS component codes [START_REF] Zhou | Reliable transmission with low complexity Reed Solomon block turbo codes[END_REF][START_REF] Le Bidan | Some results on the binary minimum distance of Reed-Solomon codes and block turbo codes[END_REF]. Actually, RS codes are a widely used subclass of non-binary BCH codes. For this reason, the code construction approach dedicated to BCH codes can be naturally applied to RS codes. It is possible to design efficient TPC architectures using RS component codes from the BCH-TPC architecture know-how [START_REF] Piriou | Efficient architecture for Reed Solomon block turbo code[END_REF]. [START_REF] Berrou | Near Shannon limit error correcting coding and decoding: Turbo Codes[END_REF] At present, the main challenge for hardware implementation of turbo decoder is to meet the demand for even higher data rates. Indeed, the increasing demand of high data rate and reliability in modern communication systems is pushing next-generation standards toward error correction schemes allowing high throughput decoding with near Shannon limit performance. Thus, data rates above 1Gbps and 10Gbps will be required for This work was supported with funds from the Brittany Region under the competitiveness pole InterAccess the next generation standards of wireless technologies (WIFI, WIMAX, DVB-H…) and data transmission over Passive Optical Networks (PON), respectively. In order to achieve high throughput and to reduce latency, parallel implementations become mandatory. In 2006, a full-parallel architecture for the turbo decoding of product codes was proposed in [START_REF] Jego | Full-parallel architecture for turbo decoding of product codes[END_REF]. The major advantage of this full-parallel architecture is that it enables the memory block between each half-iteration to be replaced by simple interconnection networks. Moreover, the latency of the turbo decoder is strongly reduced. In such a fullparallel architecture, most of the hardware complexity remains in the duplicated Soft Input Soft Output (SISO) decoders.

In this paper, the implementation onto FPGA devices of an ultra high throughput (31,29) 2 RS turbo decoder is presented. The purpose of this prototype is to demonstrate that an RS turbo decoder can effectively achieve data rates above 1Gbps. To the authors' knowledge, this is the first implementation of FEC decoder onto an FPGA target at information throughputs superior to 1Gbps. Previously, a BTC decoder included in a 12.4Gbps optical setup was described in [START_REF] Darabiha | A 3.3-gbps bitserial block-interlaced min-sum ldpc decoder in 0.13-um cmos[END_REF]. Since only a part of the transmitted data is actually coded, the information throughput of the decoder is 156Mb/s. In 2007, a LDPC code decoder ASIC implementation at an information throughput equal to 2.4Gbps was proposed in [START_REF] Mizuochi | Experimental demonstration of net coding gain of 10.1 db using 12.4 Gbps block turbo code with 3-bit soft decision[END_REF].

The remainder of the paper is organized as follows. Section II briefly recalls the basic principles of turbo decoding for RS product codes. The challenging issue of designing high data rates full parallel turbo decoding architecture is developed in Section III. Section IV describes an implementation of the resulting turbo decoder and its experimental setup onto FPGA devices. Finally, characteristics and performance in terms of BER of the prototype are detailed in Section V.

II. TURBO DECODING OF RS PRODUCT CODES

A product code is a concatenation of systematic linear block codes. The product code inherits the properties of the elementary codes that it is composed of. Let us consider two identical systematic linear block codes C having parameters (n, k, d min ), where n, k and d min stand for code length, number of information symbols and minimum Hamming distance, respectively. The parameters of the resulting product code are given by: n p = n 2 , k p = k 2 , d min p = d min 2 and R p = R 2 (code rate). Thus, it is possible to construct powerful product codes using linear block codes. When C is an RS code over Galois field GF(2 m ), we obtain an RS product code over GF (2 m ). Digital communication systems usually employ some form of binary signaling. A binary expansion of the RS product code is then required for transmission. The extension field GF(2 m ) forms a vector space of dimension m over GF [START_REF] Gallager | Low Density Parity Check Codes[END_REF]. A binary image of the A 5-Gbps FPGA prototype of a (31,29) 2 Reed-Solomon Turbo Decoder Camille Leroux, Gérald Le Mestre, Christophe Jégo, Patrick Adde and Michel Jézéquel Institut TELECOM, TELECOM Bretagne, CNRS Lab-STICC Electronics department, 29238 Brest Email: firstname.surname@telecom-bretagne.eu RS product code is thus obtained by expanding each code symbol of the product code matrix into m bits using a polynomial basis. The binary image is a product code with length mn 2 , dimension mk 2 and minimum distance at least as large as the symbol level minimum distance.

Turbo decoding process involves sequentially decoding the rows and columns of the product code matrix and by exchanging soft information between the decoders until a reliable decision can be made on the transmitted bits. In this work, SISO decoding of the RS component codes is performed at the bit-level using the Chase-Pyndiah algorithm. First introduced in [START_REF] Pyndiah | Near optimum decoding of product codes[END_REF] for binary BCH codes and latter extended to RS codes in [START_REF] Aitsab | Performance of Reed-Solomon block turbo codes[END_REF], the Chase-Pyndiah decoder consists of a softinput hard-output Chase-2 decoder [START_REF] Chase | A class of algorithms for decoding block codes with channel measurement information[END_REF] augmented by a softoutput computation unit.

III. FULL-PARALLEL RS TURBO DECODING ARCHITECTURE

With current technologies, parallel decoding architectures are the only solution to achieve information rates superior to 1Gbps. An easy architectural solution is to duplicate the elementary decoders and the corresponding interleaving memories in order to achieve the given throughput. This results in a turbo decoder with unacceptable cumulative area. Thus, smarter parallel decoding architectures have to be designed in order to better trade off performance and complexity under the constraint of a high-throughput.

A. Previous work

Many turbo decoder architectures for product codes have been previously designed. The classical approach involves decoding all the rows or all the columns of a matrix before the next halfiteration. When an application requires high speed decoders, an architectural solution is to cascade SISO elementary decoders for each half-iteration. In this case, memory blocks are necessary between each half-iteration to store channel data and extrinsic information. Each memory block is composed of four memories of mN 2 soft values. Thus, duplicating a SISO elementary decoder results in duplicating also the memory block which is very costly in terms of silicon area. In 2002, a new architecture for turbo decoding product codes was proposed in [START_REF] Cuevas Ordaz | New architecture for high data rate turbo decoding of product codes[END_REF]. The idea is to store several data at the same address and to perform semi-parallel decoding to increase the data rate. However, it is necessary to process these data by row and by column. Let us consider l adjacent rows and l adjacent columns of the initial matrix. The l 2 data constitute a word of the new matrix that has l 2 times fewer addresses. This data organization does not require any particular memory architecture. The results obtained show that the turbo decoding throughput is increased by l 2 when l elementary decoders processing l data simultaneously are used. Turbo decoding latency is divided by l. The area of the l elementary decoders is increased by l 2 /2 while the memory is kept constant.

B. Full-parallel decoding principle

All rows (or all columns) of a matrix can be decoded in parallel. If the architecture is composed of 2N elementary decoders, an appropriate treatment of the matrix allows the elimination of the reconstruction of the matrix between each decoding. Specifically, let i and j be the indices of a row and a column of the N 2 matrix. In full-parallel processing, the row decoder i begins the decoding by the soft value in the i th position. Moreover, each row decoder processes the soft values by increasing the index by one modulo N. Similarly, the column decoder j begins the decoding by the soft value in the j th position. In addition, each column decoder processes the soft values by decreasing the index by one modulo N. Actually, the full-parallel decoding of turbo product code is possible thanks to the cyclic property of RS codes. Indeed, every cyclic shift c'

= (c N-1 , c 0 , .., c N-3 , c N-2 ) of a codeword c = (c 0 , c 1 , .., c N-2 , c N-1
) is also a valid codeword in a cyclic code. Therefore, only one clock period is necessary between two successive matrix decoding operations. The full-parallel decoding of a product code matrix N 2 is detailed in Figure 1. A similar strategy was previously presented in [START_REF] Chi; Parhi | High speed VLSI architecture design for block turbo decoder[END_REF] where memory access conflicts are resolved by means of an appropriate treatment of the matrix. However, interleaving memory is still required. The elementary decoder latency depends on the structure of the decoder (i.e. number of pipeline stages) and the code length N. Here, as the reconstruction matrix is removed, the latency between row and column decoding is null.

N rows of N soft values N columns of N soft values Soft value index (i+1) = i + 1 mod N index (j+1) = j -1 mod N i j

C. Full-parallel architecture for RS product codes

The major advantage of our full-parallel architecture is that it enables the memory block of 4mN 2 soft values between each half-iteration to be removed. However, the codeword soft values exchanged between the row and column decoders have to be routed. One solution is to use an interconnection network for this task. In our case, we have chosen an Omega network. The Omega network is one of several interconnection networks used in parallel machines [START_REF] Lawrie | Access and alignment of data in an array processor[END_REF]. It is composed of log 2 N stages, each having N=2 exchange elements. In fact, the Omega network complexity in terms of number of interconnections and of 2x2 switch transfer blocks is N x log 2 N and N/2 x log 2 N, respectively. For example, the equivalent gate complexity of a 31x31 network can be estimated to be 200 logic gates per exchange bit. Figure 2 depicts a full-parallel architecture for the turbo decoding of product codes. It is composed of cascaded modules for the turbo decoder. Each module is dedicated to one iteration. However, it is possible to process several iterations by the same module. In our approach, 2N elementary decoders and two interconnection blocks are necessary for one module. An interconnection block is composed of two Omega networks exchanging the R and R iter soft values. Since the Omega network has low complexity, the full-parallel turbo decoder complexity essentially depends on the complexity of the elementary SISO decoder.

Interconnection block

A module for one iteration Fig. 2. Full-parallel architecture for decoding of RS product codes.

D. Elementary SISO decoder architecture

The block diagram of an elementary SISO decoder is shown in Figure 3, where iter stands for the current half-iteration number. R iter is the soft-input matrix computed from the previous half-iteration whereas R denotes the initial matrix delivered by the receiver front-end (R iter =R for the 1 st halfiteration). W iter is the extrinsic information matrix. α iter is a scaling factor that depends on the current half-iteration and is used to weight the influence of the extrinsic information during the first iterations. The decoder architecture is structured in three pipelined stages identified as reception, processing and transmission units [START_REF] Cuevas Ordaz | New architecture for high data rate turbo decoding of product codes[END_REF]. During each stage, the N soft values of the received word R iter are processed sequentially in N clock periods. The reception stage computes the initial syndromes S i and finds the L r least reliable bits in the received word. The main function of the processing stage is to build and then to correct the N ep error patterns obtained from the initial syndrome and to combine the least reliable bits. Moreover, the processing stage also has to produce a metric (Euclidean distance between error pattern and received word) for each error pattern. Finally, a selection function identifies the maximum likelihood codeword d and the competing codewords c (if any). The transmission stage performs different functions: computing the reliability for each binary soft value, computing the extrinsic information and updating the received soft values. The N soft values of the codeword are thus corrected sequentially. The decoding process needs to access the R and R iter soft values during the three decoding phases. For this reason, these words are implemented in six Random Access Memories (RAM) of size qmN controlled by a finite state machine.

IV. IMPLEMENTATION OF AN RS TURBO DECODER IN AN ULTRA HIGH RATE COMMUNICATIONS SETUP

The purpose of this prototype is to demonstrate that an RS turbo decoder can effectively achieve information rates until 5Gbps.

A. 5Gbps experimental setup

The experimental setup is composed of a Dinigroup board [17] that includes six Xilinx Virtex-5 LX330 FPGAs. A Xilinx Virtex-5 LX330 FPGA contains 51,840 slices that can emulate up to 12 million gates of logic. It should be noted that Virtex-5 slices are organized differently from previous generations. Each Virtex-5 slice contains four Look Up Tables (LUTs) and four flip-flops instead of two LUTs and two flip-flops in previous generation devices. The board is hosted on a 64-bit, 66MHz PCI bus that enables communication at full PCI bandwidth with a computer.

Figure 4 shows the different components of the digital communication setup implemented onto six FPGAs. One FPGA is dedicated to the component implementation of the transmission part. The five other FPGAs are dedicated to the implementation of the RS turbo decoder. One decoding iteration was implemented onto each FPGA resulting in a 5 full-iteration turbo decoder. Each decoding module corresponds to a fullparallel architecture dedicated to the decoding of a matrix of 31 x 31 coded soft values. We recall here that a coded soft values over GF(32) is mapped onto 5 LLR values, each LLR being quantized on 5 bits. The decoding process needs to access the 31 coded soft values from each of the matrices R and R iter during the three decoding phases of a half-iteration as explained in section III. For these reasons, 31 x 5 x 5 x 2 = 1,550 bits have to be exchanged between the decoding modules during each clock period f 0 =37.5MHz. The Dinigroup board offers 200 chip to chip LVDS for each FPGA to FPGA interconnect. Unfortunately, this number of LVDS is insufficient to enable the transmission of all the bits between the decoding modules. To solve this implementation constraint, we have chosen to add SERializer/DESerializer (SERDES) modules for the parallel-toserial conversions and for the serial-to-parallel conversions in each FPGA. Indeed, SERDES is a pair of functional blocks commonly used in high speed communications to convert data between parallel data and serial interfaces in each direction. SERDES modules are clocked with f 1 = 8f 0 = 300MHz and operate at 8:1 serialization or 1:8 deserialization. In this way, all data can be exchanged between the different decoding modules. Thanks to 200 chip to chip LVDS and the SERDES modules, data are exchanged between the different FPGAs at a throughput of 58,125Gbps while the working frequency and the information rates are only 37.5MHz and 5Gbps, respectively.

A. The transmission part of the prototype

A Pseudo Random Generator (PRG) sends out 145 (29*5) pseudo random data streams at each clock period (f 0 ). This module is composed of flip-flops and XOR gates that generate 145 parallel outputs. A (31,29) 2 RS encoder processes the 145 data streams in parallel. The original proposed encoding scheme uses both sequential and parallel encoders. k sequential encoders process k rows and a parallel encoder processes n columns. Sequential encoders are conventional implemented by a set of linear feedback shift registers (LFSR). Alternatively, a linear block code can be encoded by multiplying the information vector k by the generator matrix. This can be easily implemented by an XOR gate arborescence which is a parallel combinatorial encoder. The classical sequential approach requires a kn memory block between row and column encoding. The original proposed encoding scheme avoids the use of memory between row and column processing. Moreover, the depicted architecture decreases the latency from kn to only one symbol. 29 conventional sequential encoders are cascaded with one parallel encoder. 4805 encoded data (equivalent to a matrix 31 x 31 of 5-bit symbols) are generated in 31 clock periods (f 0 ). The noise generator models 155 uncorrelated White Gaussian Noise (WGN) samples and adds them to the previously encoded data. We have chosen to integrate an AWGN channel from a white Gaussian noise generator adapted to hardware implementation. High accuracy, speed and low-cost hardware are reached by combining the Box-Muller and Central limit methods [START_REF] Danger | Efficient FPGA implementation of gaussian noise generator for communication channel emulation[END_REF]. Each output sample is a 5 bit vector resulting in 775 (155*5) bits to be sent in one clock period (f 0 ). The Signal to Noise Ratio (SNR) is controllable via a computer that calculates the equivalent variance square root σ and sends it to the transmission part thanks to a PCI bus.

B. The transmission part of the prototype

One decoding iteration was implemented onto one FPGA resulting in a prototype that contains 5 full-iteration turbo decoders as previously explained. Full-parallel decoder architecture has been detailed in section III. It was observed that most of the complexity remains in the duplicated SISO decoders. A complexity analysis detailed in [START_REF] Leroux | Towards Gbps turbo decoding of product code onto an FPGA device[END_REF] allowed us to select a configuration to be implemented on the SISO decoder. All the soft information within the decoder is quantized and processed with q = 5 bits (1 sign bit and 4 reliability bits). The SISO decoder architecture is structured in three pipelined stages identified as reception, processing and transmission units as described in subsection III.D. Each stage processes N = 31 symbols in N = 31 clock periods (f 0 ). In our architecture, the number of least reliable bits in a received word R iter is equal to three (L r = 3). The processing unit computes the syndrome of the N ep =8 error patterns and their metric values. It also sequentially selects the decided word d and one competing codeword c. Finally, the transmission unit calculates extrinsic information W iter and soft output R iter+1 . The multi-stage dynamically reconfigurable interconnection network omega is based on a circular permutation principle. The connecting scheme consists of circularly shifted input data. One decoding iteration need 2N = 62 SISO decoders and 2 interconnection omega networks. Between two halfiterations, two omega networks are required to transmit R and R iter values. However, the communication resources represent less than 1 % of the total complexity of one decoding iteration. A PRG module is also implemented onto each FPGA of the reception part. It generates the exact same data as in the transmitter in order to compare data before and after each decoding iteration. A BER block is used to measure the error rate comparing data from the PRG and the hard fullparallel decoder output D iter . The number of errors and the total number of decoded data are then transmitted to the computer for a SNR value thanks to the PCI bus.

V. CHARACTERISTICS AND PERFORMANCE

OF THE (31,29) 2 RS TURBO DECODER

A. Characteristics of the 5Gbps experimental setup

The total occupation rate of the FPGA that contains the transmission part is 79%. This corresponds to 164,953 Virtex-5 slice LUT-Flip Flop pairs. In addition, DSP resources for the transmission part take up 155 out of 192 (80%). No BlockRAM of 18Kbs are assigned. The PRG module needs only 1,900 slice Flip-Flops and 145 slice LUTs. The (31,29) 2 RS encoder module occupies 10,728 slice Flip-Flops and 4,471 slice LUTs. It represents only 5 % and 2% in terms of slice Flip-Flops and slice LUTs of the total design complexity, respectively. It means that the most important module of the transmission part in terms of complexity is the AWGN channel hardware emulator. This is due to the fact that the AWGN channel is composed of 155 uncorrelated White Gaussian Noise (WGN) elements. Thus, computational resources of the channel module take up 57,741 slice Flip-Flops and 148,560 slice LUTs. The total occupation rate are about 28% and 72% of a Xilinx Virtex-5 LX330 FPGA for slice Flip-Flops and slice LUTs, respectively. Moreover, 155 DSP resources are necessary for the 155 multiplication operations between each WGN sample and the variance square root value σ given by the computer. The rest of the FPGA is occupied by one SERDES module.

One decoding iteration was implemented on each of the five other FPGAs. These FPGAs are occupied by the same design. The total occupation rate of each FPGA that contains one decoding iteration module, two SERDES modules, one PRG module and one error estimation module is 46%. This corresponds to 97,383 Virtex-5 slice LUT-Flip Flop pairs. Note that the decoding module represents only 50 % in terms of slices of the total design complexity. Actually, a decoding module for one iteration is composed of 31 x 2 = 62 elementary decoders and 2 interconnection blocks. Each elementary decoder uses information quantized on 5 bits with N ep = 8 error patterns and only 1 competing codeword. These reduced parameter values allow a decrease in the required area for a performance degradation which remains around 0.5dB as we will show in the next subsection. Thus a (31, 29) RS elementary decoder occupies 729 slice LUTs, 472 slice Flip-Flops and 3 BlockRAM of 18Kbs. An interconnection block occupies only 2,325 slice LUTs. Computation resources of one decoding iteration module take up 29,295 slice Flip-Flops and 49,848 slice LUTs. It means that the occupation rates are about 14% and 24% of a Xilinx Virtex-5 LX330 FPGA for slice registers and slice LUTs, respectively. In addition, memory resources for the decoding module take up 186 BlockRAM of 18kbits. It represents 32% of the total BlockRAM available in the Xilinx Virtex-5 LX330 FPGA. Note that one BlockRAM of 18kbits is allocated by the Xilinx tool ISE to memorize only 31x5x5 = 775 bits in our design. The occupation rate of each BlockRAM of 18kbits is then only about 4%. As input data are clocked with f 0 =37.5 MHz resulting in an input data rate of T in =5.715Gbps. However, taking into account the code rate R=0.875, the output information rate becomes T out =5Gbps. In conclusion, the implementation results showed that a turbo decoder dedicated to the (31, 29) 2 RS product code can effectively be integrated to our 5Gbps experimental setup.

The prototyping board used for our experimental setup is the major limitation in terms of information rate. Indeed, the board offers only 200 chip to chip LVDS for each FPGA to FPGA interconnect. But, it is insufficient to exchange 1,550 bits between the FPGAs in one clock period. To solve this constraint, we have implemented SERDES modules that operate at 8:1 serialization or 1:8 deserialization. This configuration requires to clock the SERDES with a frequency f 1 =8f 0 . It means that we have to take the ratio between f 0 and f 1 into account to find the frequencies. Values equal to 37.5MHz and 300MHz have been finally selected for f 0 and f 1 , respectively. However, one parallel decoding iteration decoding can be worked until a frequency f 0 equal to 85MHz onto a Xilinx Virtex-5 LX330 FPGA. As input data are clocked with f 0 = 85MHz, an information rate of 10Gbps is obtained. In particular, it means that a turbo decoder dedicated to the (31,29) 2 RS product code can effectively be integrated to the physical layer of a 10Gbps optical access network.

B. (31,29) 2 RS turbo product code performance

In [START_REF] Leroux | Towards Gbps turbo decoding of product code onto an FPGA device[END_REF], a complexity analysis of the BCH SISO elementary decoder leads to a low complexity decoder architecture for an acceptable performance degradation. As the architecture of an RS SISO elementary decoder is similar, an equivalent complexity analysis was done. Table 1 gives the characteristics of the prototype dedicated to the turbo decoding of a (31,29) 2 RS product code. The characteristics of an optimal Chase-Pyndiah algorithm in terms of BER performance are also given for comparison. By analysis the algorithm, five parameters appear to directly affect the complexity. The objective is to reach a better complexity/performance trade off in an ultra high-throughput context. The value of the scaling factor α iter depends on the current iteration in optimal Chase-Pyndiah algorithm. Keeping α iter to a power of 2 for each iteration enables to remove a multiplication operation that becomes a simple bit shifting. Therefore, the elementary decoder area is decreased by 8% for a loss of BER performance inferior to 0.1dB. Similarly, eight iterations are classically done for the turbo process. However, significant gains are obtained during the first 5 iterations. The quantification level q has also a notable impact on decoding process and has to be chosen considering the potential loss in terms of performance. The same remark can be done for N ep , L r and c. The complexity analysis led to a low complexity RS SISO elementary decoder (-30%) to be duplicated in the fullparallel turbo decoder.

Simulated performance of optimal Chase-Pyndiah algorithm and measured performance of full-parallel turbo decoder for the (31, 29) 2 RS product code are presented in Figure 5. Performance is evaluated by Monte-Carlo simulation and hardware emulation using our 5Gbps experimental setup, respectively. A cluster of computer working in parallel has given a reliable estimation for the optimal Chase-Pyndiah algorithm until a BER of 10 -8 . As the hardware emulation speeds up the simulations by a few orders of magnitude, prototype performance has been measured until a BER of 10 - 13 . The impact of decoding parameters on the performance is clearly shown. Indeed, we observe a degradation of 0.5dB in terms of BER performance for the prototype at BER= 10 -7 . Moreover, an error floor phenomenon occurs at lower SNR. It is not due to the asymptotic theoretical bound because we have applied the simple criterion introduced in [START_REF] Le Bidan | Some results on the binary minimum distance of Reed-Solomon codes and block turbo codes[END_REF]. This criterion improves by construction the binary minimum distance of the product code and thus the asymptotic performance. The parameter set that has been selected for our full-parallel turbo decoding architecture can explain this error floor. In particular, the number of N ep error patterns affects the BER performance for lower SNRs. The use of RS product codes for forward-error correction in ultra high-speed context has been investigated in this article. We have designed a full-parallel turbo decoding architecture dedicated to the (31, 29) 2 RS product code and a 5Gbps experimental setup. The purpose of this prototype is to demonstrate that an RS turbo decoder can effectively achieve information rates above 1Gbps. The turbo decoder architecture enables decoding of TPCs at information rates until 10Gbps onto FPGA devices. A first application can be envisioned in the fiber optic communication systems: 10Gbps data transmission over passive optical networks. Moreover, enhancing the study to a larger code ((63, 61) 2 RS product code) can increase throughput until 40Gbps and also be envisioned for 40Gbps line rate transmission over optical transport networks.

Fig. 1 .

 1 Fig. 1. Full-parallel decoding of an RS product code matrix.

Fig. 3 .

 3 Fig. 3. Block diagram of the turbo-decoder at the iter-th half-iteration

Fig. 5 .

 5 Fig. 5. (31,29) 2 RS turbo product code performance VI. CONCLUSION

Table 1 .

 1 characteristic comparison between prototype and optimal algorithm

		Optimal algorithm	Prototype
		features	features
	arithmetic representation	floating-point 5bits fixed-point
	number of iterations	iter = 8	iter = 5
	number of error patterns	N ep = 16	N ep = 8
	number of least reliable bits	L r = 4	L r = 3
	number of competing codewords	c = 16	c = 1
	value of the scaling factor α iter	0.2< α iter < 0.8 α iter ={0.25;0.5;1}

ACKNOWLEDGMENT

The authors thank E. Boutillon for providing the WGN IP. They are grateful to R. Le Bidan for providing the optimal bit error performance of iterative decoding of a (31, 29) 2 RS product code.


Hazem Moussa

Amer Baghdadi

Michel Jézéquel

Binary de Bruijn Interconnection Network for a Flexible LDPC/Turbo Decoder

This paper proposes a novel on-chip interconnection network adapted to a flexible multiprocessor LDPC/turbo decoder and based on the de Bruijn network. The main characteristics of this network -including its logarithmic diameter, scalable aggregate bandwidth, and optimized routing technique-allow it to efficiently support the communicationintensive nature of the two decoding techniques. We present a detailed hardware implementation of the routers and the network interfaces as well as the packet format and the routing algorithm. In order to evaluate the performance of the proposed network, a generic RTL VHDL description has been developed and synthesized with ST CMOS 0.18 µm technology. The flexibility and the scalability of this on-chip communication network enable it to be used in the emerging multi-code applications and standards. In addition, the results obtained for a 16-processor network demonstrate a major aggregate bandwidth of 296 Gbps with a relative small area of 3.56 mm².

INTRODUCTION

To reduce the error rate in digital communications, with a lower signal-to-noise ratio (closer to the Shannon limit), iterative decoding processing algorithms have recently emerged [START_REF] Berrou | Near Shannon Limit Error-Correcting Coding and Decoding: Turbo-Codes[END_REF] [START_REF] Mackay | Good error-correcting codes based on very sparse matrices[END_REF] in order to improve the error correction of the transmitted data. Indeed, thanks to appropriate channel coding (parallel recursive systematic convolutional codes for turbo codes and linear block codes with a sparse parity check matrix for LDPC codes), the noisy received information can be recovered by the decoding process.

Although these two kinds of error correcting codes have similar performance, they differ in certain cases where sometimes turbo codes are better than LDPC codes or vice versa [START_REF] Lestable | Block-LDPC codes vs duo-binary turbo-codes for european next generation wireless systems[END_REF]. This is why in many telecommunication applications, it may be desirable to have a flexible system able to support a variety of error correction standards or a variety of codes [START_REF] Gunnam | Decoding of Quasi-cyclic LDPC Codes Using an On-the-Fly Computation[END_REF] [START_REF] Lin | A dual mode channel decoder for 3GPP2 mobile wireless communications[END_REF], and in particular both LDPC and turbo codes.

A few multi-code implementations of channel decoders have recently been proposed to provide an answer to the need for flexibility.

In [START_REF] Lin | A dual mode channel decoder for 3GPP2 mobile wireless communications[END_REF], the implementation of a dual mode decoder for 3GPP2 was described, supporting the Viterbi decoding of a 256 state convolutional code and the MAP decoding of a turbo code with many shared hardware components.

The authors in [START_REF] Tarable | Mapping interleaving laws to parallel turbo and LDPC decoders architectures[END_REF] propose a generic approach to support all possible interleaving for any type of LDPC or turbo code. This approach presents a mapping algorithm avoiding the conflicting memory accesses generated by the interleaving laws of the codes. Indeed, an appropriate memory addressing scheme, based on spatial and temporal permutations, prevents conflicts occurring.

To the best of our knowledge, none of the existing works has presented a flexible architecture supporting LDPC and turbo codes except the contribution proposed in [START_REF] Chiu | Techniques for reconfigurable decoder for a wireless system[END_REF]. In this work, the authors target a reconfigurable decoder for a wireless system made of a number of configurable Add-Compare-Select decoders that exchange messages through multiple memory banks. However, the presented architecture requires a huge memory capacity in order to avoid the conflicting accesses. Besides, a complex interconnecting architecture, not presented, is required to achieve full flexibility with limited throughput penalty.

Finally, there is a wide and growing interest in realizing channel decoders supporting more than one code and many researches are being recently conducting in this direction.

On the other hand, the new emerging standards and applications require high performance in terms of throughput. For this purpose, the parallelization of the decoding algorithms has been widely investigated where several implementations, based on multiprocessor systems, have been proposed. Furthermore, the implementation issues of these multiprocessor architectures mainly come from the communication structure between the processing elements inside the decoder. The huge quantity of data (extrinsic information) exchanged by the processing elements requires an efficient interconnection structure, which represents a major challenge for flexible decoder implementation.

In this paper, we propose and evaluate an adapted de Bruijn network (with routers and network interfaces) as an efficient interconnection network of a flexible parallel LDPC/turbo decoder. The main features of this network architecture are its scalability enabling seamless trade-off between hardware complexity and available bandwidth for LDPC/turbo decoding and its flexibility allowing any interleaving communication profile to be implemented.

The remainder of the paper is organized as follows. Section II highlights the communications involved in parallel LDPC/turbo decoding. Then section III details the implementation This work was supported in part with funds from the ANR under the AFANA project architectures of the proposed interconnection network, including topology, routing algorithm, packet format, routers, and network interfaces. Finally, section IV presents the synthesis results and provides performance analysis and comparison. The paper is concluded in section V.

II. COMMUNICATION REQUIREMENTS OF FLEXIBLE PARALLEL ITERATIVE CHANNEL DECODER

In parallel turbo decoding, extrinsic information is iteratively and concurrently exchanged between component decoders. For each iteration, these exchanges become more and more massive with decoder level parallelism [START_REF] Muller | ASIP-Based Multiprocessor SoC Design for Simple and Double Binary Turbo Decoding[END_REF].

In addition, the communication traffic profile, which depends on the turbo code interleaving rule, can be considered as random if we target flexibility in order to support multi-mode and multistandard features. In this context, hardware implementations of parallel turbo decoders require the integration of a complex interconnection structure supporting the intensive interleaved memory accesses.

The two component decoders in Fig. 1 concurrently process the iterated frame, and for both, multiple processors can be designed for the parallel processing of frame sub-blocks. Concerning LDPC codes, their decoding process is directly related to the structure of the parity check matrices. Indeed, the codes can also be described by a bipartite graph, called a Tanner graph [START_REF] Mackay | Good error-correcting codes based on very sparse matrices[END_REF], made up of two types of nodes: variable nodes (or bit nodes) (VNs) corresponding to the bits of a codeword and check nodes (CNs) corresponding to the parity constraints (an edge in the graph represents a one in the parity check matrix).

Therefore, two categories of processors are necessary: VN and CN Processors (VNP and CNP). Thus, the architecture of a decoder (Fig. 2) can be defined as being a system made up of VNPs, CNPs, and an interconnection structure allowing the exchange of messages between processors. Mapping the Tanner graph onto this architecture involves: (1) associating one or more VN (resp. CN) to one VNP (resp. CNP) and ( 2) implementing Tanner graph edges through a corresponding message addressing in the interconnection structure.

However, similarly to parallel turbo decoders, partiallyparallel code-independent LDPC decoders, which offer the best compromise in terms of area, throughput, and flexibility [START_REF] Masera | Implementation of a flexible LDPC decoder[END_REF], introduce a problem of conflicting accesses to the memories containing the exchanged messages.

In addition, the computations being relatively simple, implementation issues mainly come from the communication structure between the VNPs and CNPs. Indeed, the communications rapidly become intensive because they depend on the number of nodes, the node degrees and the number of iterations.

In this situation, a flexible on-chip interconnection network must be designed with the aim of fully exploiting the parallelism of the LDPC/turbo decoder architecture by reducing the message latency, alleviating the memory conflicts and efficiently routing any permutation from the network input ports to its output ports.

However, even if several efficient network architectures have already been proposed to solve the communication problem inside flexible LDPC decoders [START_REF] Theocharides | Implementing LDPC decoding on a network-on-chip[END_REF] or flexible parallel turbo decoders [START_REF] Moussa | Butterfly and Benes-based on-chip communication networks for multiprocessor turbo decoding[END_REF], none of them have been designed to handle the communications in a decoder supporting LDPC and turbo codes at the same time. The de Bruijn network, based on the de Bruijn graph [START_REF] De Bruijn | A combinatorial problem[END_REF], was chosen due to its appropriateness to support the communications of a multiprocessor LDPC/turbo decoder. In fact, the structure of the network allows any permutation to be routed efficiently thanks to the path diversity offered by the network. The conflicting memory accesses are alleviated because the conflicting packets can be deviated appropriately until they attain the targeted processor rather than being blocked or buffered.

In addition, the logarithmic diameter of the network (log 2 (N) leads to small latencies and the number of routers is reduced because it is a direct network. Moreover, its recursive structure makes it highly scalable. Finally, in [START_REF] Samatham | The de Bruijn multiprocessor network: a versatile parallel processing and sorting network for VLSI[END_REF], it has been shown that the de Bruijn network is suitable for VLSI implementation. NASA developed an 8192-processor de Bruijn network used by a Viterbi decoder for the Galileo space mission to Jupiter [START_REF] Rowley | Fault-tolerant ring embedding in de Bruijn networks[END_REF].

A. Network Architecture and Processor Mapping

The de Bruijn graph consists of N=k m nodes, where each node is identified by a vector of size m in the k-ary number system. For our network, we chose the binary system. Thus, a node v m v m- As an example, Fig. 3 presents the schematic view of a complete de Bruijn network for N=8 processors (m=3). Each processor is connected to a network interface and each network interface connected to a router.

Concerning the mapping of the decoding elements on the network graph, we chose to simply allocate the first N/2 processors to the first component decoder (in the case of turbo codes) or to the VNPs (in the case of LDPC codes).

B. Routing Algorithm and Packet Arbitration

The shortest path routing algorithm [START_REF] Francalanci | High-performance self-routing algorithm for multiprocessor systems with shuffle interconnections[END_REF] was implemented thanks to its ability to reduce latencies. However, the algorithm being sequential and in order to bring the number of cycles back to one, we modified the algorithm of [START_REF] Francalanci | High-performance self-routing algorithm for multiprocessor systems with shuffle interconnections[END_REF] by unrolling the "for" loop, and developed another method to compute the shortest path output. In order to reduce the size of routers by avoiding packet buffering when packets are in conflict, the routing algorithm was supplemented by a mechanism deflecting the "youngest" packet. This requires storing "age" information (Time Stamp TS) within the packet updated at each hop. Moreover, in [START_REF] Moussa | Butterfly and Benes-based on-chip communication networks for multiprocessor turbo decoding[END_REF], a defined priority criterion is proposed to improve the packet arbitration inside a parallel turbo decoder. This criterion represents the contribution of the packets to the convergence of the iterative process. Consequently, we propose to apply this criterion also for iterative LDPC decoding with the effect of adding a priority field in the packet header.

It is worth to note that the de Bruijn network introduces some variable latencies which can be intrinsically tolerated in this case of iterative probabilistic decoding. Indeed, if an extrinsic information is used before being updated, the new value will be taken into account at the next iteration. Furthermore, the convergence of the decoding process is improved through the contribution priority field.

C. Packet Format

Packets are divided into two parts: the header and the payload. The first one is composed of the destination router identifier, the TS value, the priority field and a one-bit flag indicating whether intra symbol permutation, as defined in certain standards like WiMAX, is applied to this packet. The payload is made up of the destination memory write address and of the extrinsic information on 64 bits (the 8 LSB for LDPC codes and all the 64 bits for turbo codes). The latter corresponds to the reception of each possible symbol of a double binary turbo code (4 x 16 bits).

D. Router

The router architecture (Fig. 4) is divided into two units: operative and control units. The first one is composed of flipflops, registers, and multiplexers implementing the packet switching and selecting the packets coming from routers and those coming from the processor. The control unit generates the multiplexer selection signals as well as the validity signals according to the routing algorithm and conflicting packet arbitration. Each port of a de Bruijn router has an associated control signal used to validate the data passing through the port. Moreover, there are two output signals corresponding to the read request of the emission FIFOs contained in the network interface.

E. Network Interface

The network interface (Fig. 5) manages the packet flow in both directions, from the processor to the network and vice versa. The interface contains a reconfigurable table (RAM) for the construction of the packet header according to the parity check matrix of an LDPC code or according to the interleaving rule of a turbo code. Indeed, the permutation of the extrinsic information is done via the addressing of the packets. Thus, the de Bruijn network can support the LDPC/turbo codes of wireless standards like WiMAX by simply programming the content of the table accordingly. Moreover, packets are stored in two FIFOs at the input and output of the network interface. Finally, two multiplexers implement the intra symbol permutation.

IV. SYNTHESIS AND RESULTS ANALYSIS

For performance evaluation, we developed and synthesized a description in generic RTL VHDL of the complete system. The case study was the WiMAX standard for turbo decoding with a network of 16 processors. About LDPC decoding, we considered a 1024 bits code also mapped onto the same network. Thus, to estimate the area consumption and the aggregate bandwidth, we used ST CMOS 0.18 μm technology. We indicate for each network its topology, the target technology, the maximum clock frequency in MHz, the area in mm², and the estimated aggregate bandwidth in Gbps.

In Table I, we compare the results obtained with two networks recently proposed and based on the Butterfly and Beneš 2N-N topologies [START_REF] Moussa | Butterfly and Benes-based on-chip communication networks for multiprocessor turbo decoding[END_REF]. As the proposed routing of de Bruijn network is more complex (shortest path computation + deflection) than those of Butterfly and Beneš 2N-N networks, the resulting critical path is longer, and thus the maximum clock frequency and aggregate bandwidth are lower. On the other hand, the de Bruijn area cost is in between because at the same time, the network interfaces are bigger (due to the routing table) and the routers are smaller and fewer. However, the de Bruijn network resolves the limitations of the Butterfly and Beneš 2N-N solutions. Indeed, the conflicting packets deflection exempts the usage of buffers. Furthermore, there is no need to pre-compute the routing paths because the choice of the routers to cross is performed dynamically. To finish, the sufficient path diversity provided by the de Bruijn network, coupled to the deflection technique allow implementing optimizations based on application-specific priorities.

Next, we examine the results obtained with other previously proposed flexible LDPC decoder architectures (Table II) for a given codeword size of 1024 bits and based on an interconnection network. Frequency, area, and aggregate bandwidth being related to the target technology, it is not possible to compare directly and accurately the values of Table II. However, the de Bruijn network presents the best efficiency ratio between aggregate bandwidth and occupied area.

Regarding scalability, the recursive structure of the Beneš, the 2D-mesh, and the de Bruijn networks makes them highly scalable. On the other hand, the MDN network lacks scalability because the network structure is randomly built by a complex algorithm based on a constructive heuristic.

Although flexible, the Beneš network [START_REF] Masera | Implementation of a flexible LDPC decoder[END_REF] requires a precalculation related to the code to configure the switches. The MDN [START_REF] Kienle | Implementation issues of scalable LDPC-decoders[END_REF] uses a general approach to network design with random graphs as basic topologies. This network lacks flexibility because all permutation patterns must be simulated prior to their use in order to determine buffer size at design-time. Finally, the 2D-mesh network [START_REF] Theocharides | Implementing LDPC decoding on a network-on-chip[END_REF] is designed for a flexible decoder. However, the described network can only manage codes with a maximum codeword size of 1024 bits. Consequently, the de Bruijn network, allowing to support any LDPC code, constitutes the most flexible solution.

V. CONCLUSION

In this paper, we have presented the architecture of a new onchip interconnection network based on the de Bruijn graph and dedicated to a flexible multiprocessor LDPC/turbo decoder. We have detailed the topology and the network architecture, the packet format as well as the routers and the network interfaces. We have also proposed to use priorities specific to LDPC and turbo codes for the arbitration of conflicting packets. The unique trade-off between flexibility and performance (scalability, frequency, area) offered by the de Bruijn network demonstrates its efficiency in the context of a flexible parallel channel decoder supporting LDPC and turbo codes.

Figure 1 .

 1 Figure 1. Extrinsic information exchanges in a parallel turbo decoder

Figure 2 .

 2 Figure 2. Extrinsic information exchanges in a partially-parallel LDPC decoder

Figure 4 .

 4 Figure 4. Router architecture

Figure 5 .

 5 Figure 5. Network Interface architecture

TABLE I .

 I ASIC SYNTHESIS RESULTS OF BUTTERFLY, BENEŠ 2N-N AND DE BRUIJN NETWORKS FOR A 16-PROCESSOR TURBO DECODER

			[11]	Proposed
	Topology	Butterfly	Beneš 2N-N de Bruijn
	Technology (µm)	0.18	0.18	0.18
	Frequency (MHz)	302	416	266
	Network Area (mm²)	4.834	2.692	3.565
	Aggregate Bw. (Gbps)	308	426	296

TABLE II

 II

	.	COMPARISON ASIC SYNTHESIS RESULTS OF LDPC
			DECODERS
			[10]	[16]	[9]	Proposed
	Topology		2D-mesh MDN Beneš de Bruijn
	Technology (µm)	0.16	0.18	0.13	0.18
	Frequency (MHz)	500	200	245	266
	Network Area (mm²)	25	7.92	0.98	3.565
	Aggregate Bw. (Gbps)	400	19.2	62.72	296

…v 1 has its two output edges connected with nodes v m-1 v m-2 …v 1 0 and v m-1 v m-2 …v 1 1.


Pierre Cordesse
email: pierre.cordesse@polytechnique.edu

Samuel Kokh

Ruben Di Battista

Marc Massot

Derivation of a two-phase flow model with two-scale kinematics and surface tension by means of variational calculus

Keywords: two-phase flow, least action principle, two-scale system, surface tension

Derivation of a two-phase flow model with two-scale kinematics and surface tension by means of variational calculus

Introduction

Two-phase flow occur in many industrial processes and can display multiple flow topologies. For example a jet atomization involves separated phases close to the injector and dispersed phase downstream with in between a mixed region where complex elements such as ligaments and rings detach from the liquid core. These multi-scale phenomena are out-of-reach for direct numerical simulations in realistic configurations due to computational costs. This leads to the derivation of reduced-order models that are used for numerical simulations.

Two approaches are found in the literature to build reducedorder models for two-phase flows:

1. coupling two models, a first one for the dispersed flow by employing an element derived from the Kinetic Based Moment Method [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods[END_REF] which treats polydispersed droplets in size, velocity and temperature, and a second one dealing with the separated phase and the mixed region with either an hierarchy of diffuse interface models [START_REF] Drui | A hierarchy of simple hyperbolic two-fluid models for bubbly flows[END_REF][START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF], front tracking methods [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows[END_REF] or some level of LES on the interface dynamics [START_REF] Herrmann | A sub-grid surface dynamics model for sub-filter surface tension induced interface dynamics[END_REF].

2. employing a unified model that encompasses any flow topology. Works in this direction are found in [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Cordesse | Derivation of a two-phase flow model with twoscale kinematics and surface tension by means of variational calculus[END_REF] where a unified model accounting for micro-inertia and micro-viscosity associated to bubble pulsation is proposed.

In the present work, the second approach is retained. In [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Cordesse | Derivation of a two-phase flow model with twoscale kinematics and surface tension by means of variational calculus[END_REF], a model was designed to account for two-scale kinemat-ics in the case of a mono-disperse flow with spherical gas inclusions by means of the Least Action Principle. Our contribution involves mainly three matters: first we propose a definition of the interface using a Probability Density Function, second we propose to drop some restricting hypotheses on the topology of the bubbly flow used in [START_REF] Drui | A hierarchy of simple hyperbolic two-fluid models for bubbly flows[END_REF][START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF], third we aim at not only accounting for two-scale kinematics but also for two-scale surface tension effects.

Description of a surface by means of a PDF

Our analysis relies on the key idea that the two-phase interface is not a geometric locus. Instead, we propose to use a probability density function (PDF) f that will characterize the properties of the interface S . We will use this probabilistic description to distinguish small and large scales of the interface.

We suppose given a probability density function (PDF)

f : (t, x, ϕ, H) ∈ Ω → [0, 1], with Ω ⊂ [0, +∞) × R d × R × R.
The quantity f (t, x, ϕ, H) dxdϕdH represents the probability that at a position x and the instant t the interface can be described by a level set function that takes the value ϕ and the field of interface mean curvature takes the value H.

In the following we shall note ξ = (ϕ, H) and for the sake of simplicity we shall omit the time dependency of f in the notations. Thanks to zero-th moments of f we can then define fields over the whole space that characterize the interface as follows:

ϕ(x) = ϕ f (x, ξ ) dξ , (1a)
H(x) = H f (x, ξ ) dξ . ( 1b
)
The set ϕ -1 ({0}) defines a surface that we consider to be a probable geometrical locus for the interface S . Similarly, the H(x) will provide a probable value of the mean curvature of S when x lies in the vicinity of ϕ -1 ({0}). Let us note that thanks to (1a) it is possible to define a unit normal field

x ∈ R d → n ∈ R d by setting n(x) = ∇ϕ(x)/ ∇ϕ(x) .
Definition (1a) also provides an alternate definition H alt for the mean curvature with H alt (x) = div(n). In the following we shall consider that both definition are equivalent, i.e.

H alt = H.

We note 1 A the characteristic function of a set A ⊂ R d and R + = [0, +∞). Considering that the interface S separates two different fluid, we can define α the volume fraction of the fluid 1 as

α(x) = 1 R + (ϕ )f (x, ξ ) dξ . (2)
Let us underline that the approach that is described above allows to retrieve the classic definition of an interface by a surface. Indeed, given a surface defined by a level set function x → φ, defining the measure

f (x , ξ ) dx dξ = δ(ϕ -φ(x ))δ H -div ∇ φ(x ) ∇ φ(x )
allows to retrieve ϕ = φ and H = div (∇ φ/ ∇ φ ). Nevertheless, we emphasize that there is no physical reason to consider such case. Indeed, a sharp interface is just a model for finitely fine transition zone.

The PDF f is supposed to be able to encompass all the scales l of S , l ∈ ]0, +∞[. If obtained by a DNS, it would be then limited to the mesh minimum size, such that l ∈ ]l DN S , +∞[.

Two-scale description of a surface

Filtering of the PDF We now want to introduce a separation between large and small scales of the interface S . Let l c > l be the cutoff characteristic length scale, such that l ∈ ]l, l c ] (resp. l ∈ ]l c , +∞[) corresponds to the small (resp. large) scales of S .

We introduce a new PDF f c characterizing all the scales of S up to the cutoff length scale l c , l ∈ ]l c , +∞[. Let us note B c (x) = {x / xx < l c /2}. Since f c does not account for any scale below l c , we propose to operate a spatial averaging over B c (x) by considering the distribution defined by

1 |B c (x)| 1 Bc(x) (x )f (x , ξ) dx .
In order to define f c we propose further filter the values of f : since f c has undergone a spatial averaging, it seems reasonable to upper-bound the possible mean curvature H described by f c with a value H max . Consequently we choose to set

f c (x, ξ) = 1 |B c (x)| Bc(x) 1 R +(H max -H )f (x , ξ) dx . (3)
For the sake of consistency, we assume in the following that when l c → 0 and H max → +∞ we have f c → f in some sense. Furthermore, we propose to define H max by the curvature of the ball of radius l c /2, i.e. H max = 2/l c .

Using the moments of f c , we can define a filtered level set ϕ c , filtered mean curvature H c and a filtered volume fraction, α c by setting

ϕ c (x) = ϕ f c (x, ξ ) dξ , (4a)
H c (x) = H f c (x, ξ ) dξ (4b) α c (x) = 1 R + (ϕ )f c (x, ξ ) dξ . (4c)
We now notice that α

= 1 R + (ϕ )(f -f c + f c ) dξ , which
yields a natural decomposition of α as

α = α + α c , with α = 1 R + (ϕ )(f -f c ) dξ / (5)
Thus α is the superposition of a contribution coming from f c and a fluctuating part α. Let us emphasize that α ∈ 2 illustrates this remark by showing the effect of the spatial averaging on the ball B c (x). In such example, α equals either

[-1, 1] while α ∈ [0, 1] and α c ∈ [0, 1]. Figure
1 or 0, α c ∈ [0, 1] but α(x) < 0 when α(x) = 0.
We can now also define the interfacial density area inside the ball B c by integrating the local mean curvature, H over B c (x), more precisely we set In a similar manner, we define Σ c the filtered interfacial area density, Σ c , as

Σ(x) = 1 |B c (x)| ϕ =0 H f (x , ξ )1 Bc(x) (x )dx dξ . ( 6
)
Σ c (x) = 1 |B c (x)| ϕ =0 H f c (x , ξ )1 Bc(x) (x )dx dξ . (7)
Consequently, using similar lines as for the volume fraction, Σ can be split into two contributions, Σ c and a fluctuating part Σ defined by

Σ = Σ -Σ c . (8)
As opposed to the volume fraction, the interfacial area density, being defined with respect to the B c , is not affected by the spatial averaging on B c . Therefore, we have

0 ≤ Σ c ≤ Σ 0 ≤ Σ ≤ Σ. (9)
In the same lines, we define the filtered mean curvature H c and decompose the mean curvature into two contributions,

H = H c + H. ( 10
)
As a result of the filtering proposed in Equation (3) on the mean curvature maximum value, we have

H c < H max = 2/l c (11)
It is possible to characterize the difference between α and α c . Indeed, we have

f (x) -f c (x) = 1 |B c | Bc(x) (f (x) -f c (x ))dx . (12)
Since x ∈ B c (x), if we suppose that f is smooth with respect to the space variable, a Taylor expansion yields

f (x) -f c (x ) = ∂f ∂x l c r + O l 2 c (13) where x -x = l c r, r ∈ B c (0). Thus f (x) -f c (x) = O l 2 c , (14)
and finally we obtain

α(x) = α c (x) + O l 2 c (15)

Two-scale two-phase flow interface

We now turn to the derivation of a two-phase flow model using our two-scale interface description. We assume that the pressure and velocity of each component are at equilibrium. We consider a set of bulk variables: the total density ρ, the mass fraction of fluid 1 Y and the velocity v, then we consider variables that characterize the position geometry of the interface. In classic approaches like [START_REF] Kapila | Two-phase modeling of deflagration-todetonation transition in granular materials: Reduced equations[END_REF] this characterization is achieved by using the volume fraction which can only capture large scale interface features: small scales are lost. Such model cannot account accurately for complex two-phase flow topology variations across several scales.

Our goal is two-fold: first we aim at enriching the description of the interface by solving supplementary geometric variables such as the interfacial density area, Σ and the mean curvature H. Second, we propose to account for a two-scale kinematic evolution.

We suppose the two-phase flow interface to be defined by a PDF f that describes all the scales of the interface. Let us consider a cut-off length scale l c . We can now express variables that characterize the interface in terms of f and f c .

The Lagrangian energy L of the system is the difference between the kinetic energy, K and the potential energy, U , of the system. We separate bulk and interfacial contribution of the potential energy, denoting the bulk potential U b and the interfacial energy U i . The bulk kinetic energy contribution K is defined through the bulk variables only, and therefore no distinction between small and large scales are needed.

K = 1 2 ρ v 2 (16)
We suppose that the bulk potential energy takes the form

U b = ρe(ρ, Y, α).
Let us note that U b depends on bulk variables but also on α.

In order to account for surface tensions effects we also add a potential interfacial energy U i . If we have access to all scales of the interface, we can express the potential interfacial energy by

U i = 1 2 σ ∇α ( 17
)
where σ is the surface tension coefficient of the mixture and α, defined in Equation ( 2), is accounting for small and large scales. After filtering small scales, ∇α c cannot describe small scale surface tension. By (4c) and ( 5), we have

∇α = ∇α c + ∇α. ( 18
)
This suggests that the interfacial energy, accounting for scale separation can be expressed as

U i = 1 2 σ ∇α c + W, ( 19
)
where W is a subscale interfacial energy W that we intend to model using the geometric characterization of the interface.

We suppose that small scale interfacial energy variations are involved when there is a variation of the fluctuation of interfacial density area Σ. Under such variations we suppose that surface tension forces are oriented along the normal n to the interface and proportional to the local mean curvature. We also assume that for small scales the interface can only undergo displacement δh n that are oriented along the normal to the interface. The infinitesimal energy variation involved with the infinitesimal work of such forces read the interface δh n, reads

δW = γHδh, ( 20
)
where γ is a constant characterizing the surface tension properties of the mixture. In [START_REF] Cordesse | Derivation of a two-phase flow model with twoscale kinematics and surface tension by means of variational calculus[END_REF], relationships realted to geometric variables describing the interface have been derived when the interface is subjected to a small displacement δh n. These relations can be expressed as follows

δ Σ = -2H Σδh, δα = Σδh, (21)
where δ characterizes an infinitesimal variation. Let us emphasize that (21) involves α, and not α. This is due to the fact that the small scale interface is only accurately described by α as shown in Figure 2. Therefore the inifinitesimal variation of small scale interfacial energy takes the form

δW = 1 2 γ δ Σ Σ . (22)
This suggests that

W = 1 2 γ ln Σ Σ 0 (23)
where Σ 0 is defined as the interfacial density area of B c (0), Σ 0 = 6/l c , i.e. the interfacial density area associated to the maximal curvature that can be capture at any point x by α c . The limit case Σ(x) → 0 in (23) has to be discarded as in this case there is no interface in the vicinity of x. We propose to detect such situations using the values of H. This suggests to adopt an alternate definition for W by considering

W = γ(H) 2 ln Σ Σ 0 , (24)
γ(H) = γ if H ≥ H max , 0 otherwise. ( 25
)
We now turn to the definition of the small scale kinetic energy K i by setting

K i = 1 2 m(D t h) 2 . ( 26
)
The coefficient m has the dimension of a mass. Injecting the Equation (21), it yields

K i = 1 2 m (D t α) 2 Σ2 . ( 27
)
We thus propose to define the interfacial energy as

K i + U i = 1 2 σ ∇α (28) = 1 2 m (D t α) 2 Σ2 + 1 2 σ ∇α c -γ(H) ln Σ Σ 0 (29)
where Equation ( 28) gives the form of the interfacial energy in terms of geometric variables that cover the whole range of scales whereas Equation (29) shows the decomposition of the interfacial energy induced by the introduction of a cutoff length scale and therefore is defined with respect to filtered and fluctuating geometric quantities.

Relation (15) suggests that we can substitute α by α c in the bulk potential energy. Therefore when accounting for a scale separation between small and large scale, the Lagrangian of the system is now a function of ρ, Y, v, α c , Σ and can be written as

L ρ, Y, v, α c , Σ = 1 2 ρv 2 + 1 2 m (D t α c ) 2 Σ2 + 1 2 σ ∇α c -γ(H) ln Σ Σ 0 -ρe(ρ, Y, α c ). (30)
From now on, we shall assume the mean curvature x → H to be a fixed given field. Discarding this restriction will be the matter of future works. To simplify notation, we will now drop the notations for fluctuating and filtered quantities on the volume fraction and the interfacial density area and suppose m, σ and γ to be constant.

Extremization of the Action

We now follow classic lines of the Least Action Principle. Consider B(t) ⊂ R 3 the volume occupied by the fluid for t ∈ [t 0 , t 1 ]. Let X ∈ B(t 0 ) be the Lagrangian coordinates associated with the reference frame at instant t = t 0 , then we note (t, X) → ϕ L the position of the fluid particle whose position is X at t = t 0 . If (t, x) → b is any Eulerian field it can be associated with the Lagrangian field (t, X) → b L by setting b(ϕ L (X, t), t) = b L (X, t). The flow can be fully characterized by (t, x) → (ρ, v, Y, α, Σ) or equivalently by x → (Y, α, Σ) and (t, X) → ϕ L if ϕ L complies with the mass conservation equation.

For a given transformation of the medium x → (Y, α, Σ) and

(t, X) → ϕ L , let (t, x, λ) → (Y λ , α λ , Σ λ ) and (t, X, λ) → ϕ L be a family of medium transformations parametrized by λ ∈ [0, 1]. We note Ω(λ) = (t, ϕ L (t, X, λ))|X ∈ B(t 0 ), t ∈ [t 0 , t 1 ]
and we require these fields to satisfy constraints pertaining to mass conservation

∂ρ λ ∂t + ∇ • (ρ λ v λ ) = 0, ∂ρ λ Y λ ∂t + ∇ • (ρ λ Y λ v λ ) = 0 (31)
10 th International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019 supplemented by constraints regarding the topology evolution

D t Σ λ + 2H Σ λ D t h λ = 0, D t α λ -Σ λ D t h λ = 0, (32)
and finally classic boundary constraints

(Y λ , α λ , Σ λ )(t, x, λ = 0, 1) =(Y, α, Σ)(t, x), (33a) ϕ L (X, t, λ = 0, 1) =ϕ L (X, t), ( 33b
) (Y λ , α λ , Σ λ )(t, x, λ) |(t,x)∈∂ Ω(λ) =g(Y, α, Σ)(t, x), (33c)
ϕ L (X, t, λ) |(t,X)∈∂([t0,t1]×B(t0)) =ϕ L (X, t). (33d)
Following standard lines, this family of transformation yields a family of infinitesimal transformations defined as follows

δ λ ϕ(t, ϕ L (t, X)) = ∂ ϕ L ∂λ t,X
(t, X, λ = 0), (34a)

δ λ b(t, x) = ∂ b ∂λ t,x (t, x, λ = 0), (34b)
for b ∈ {ρ, Y, v, α, Σ}. Let us now define the Hamiltonian action A associated with Ω for the family of transformations

(t, x, λ) → (Y λ , α λ , Σ λ ) and (t, X, λ) → ϕ L A (λ) = Ω(λ) L(ρ λ , Y λ , v λ , D t h λ , α λ , Σ λ , ∇α λ ) dxdt. (35)
The Least Action Principle states that a physical transformation of the system verifies

dA dλ (0) = 0. (36)
Relation (36) will provide the motion equations of the flow.

In order to obtain a set of partial differential equations, we need to express dA /dλ . Using definition (34b) we can write

dA dλ (0) = Ω(0) ∂L ∂ρ δ λ ρ + ∂L ∂Y δ λ Y + ∂L ∂v δ λ v + ∂L ∂(D t h) δ λ (D t h) + ∂L ∂α δ λ α + ∂L ∂Σ δ λ Σ + ∂L ∂(∇α) δ λ (∇α) dxdt. ( 37
)
Applying (34b) with the constraints (31) allows to express following relations between the infinitesimal variations

δ λ ρ = -∇ • (ρδ λ ϕ) , (38a)
δ λ Y = -∇Y T δ λ ϕ, (38b)
δ λ v = D t (δ λ ϕ) -δ λ ϕ T ∇ v, (38c)
δ λ (D t h) = 1 Σ δ λ (D t α) - D t h Σ δ λ Σ. (38d)
Recasting relations (38) into (37) provides

Ω(0) [A T δϕ + B δα + C δΣ] dxdt = 0, (39)
A T = ∂ t ∂L ∂v + ∇ • [( ∂L ∂v T v T + (∇v) T ∂L ∂v T + ∂ t 1 Σ ∂L ∂D t h ∇α + ∇ • 1 Σ ∂L ∂D t h (∇α) T v + 1 Σ ∂L ∂D t h (∇v) T ∇α -ρ ∇ ∂L ∂ρ T + ∂L ∂Y ∇Y, (40a)
B = ∂L ∂α -∇ • ∂L ∂∇α -∂ t 1 Σ ∂L ∂D t h -∇ • 1 Σ ∂L ∂D t h v , (40b)
C= ∂L ∂Σ - 1 Σ ∂L ∂D t h D t h. (40c)
We can conclude that the Least Action Principles applied to the Lagrangian energy defined by (30) yields the following equations of motion

A = 0, B = 0, C = 0. (41)
Let us further express the equations of motions into a more familiar form. With the definition (30) of L one then has

∂ ρ L = |v| 2 2 -e -ρ∂ ρ e, ∂ Y L = -ρ∂ Y e, (42a)
∂ v L = ρv, ∂ (Dth) L = mD t h, (42b)
∂ Σ L = - γ Σ , ∂ α L = -ρ∂ α e, (42c)
∂ (∇α) L = σ ∇α ∇α . (42d)
We obtain

A= ∂ t (ρv) + ∇ • ρvv T + ρ∇v • v + ∂ t m Σ D t h∇α +∇ • m Σ D t h (∇α) T v + m Σ D t h(∇v) T ∇α -ρ∇[1/2|v| 2 -e -ρ∂ ρ e] -ρ∂ Y e∇Y, (43a)
B = ∂ t m Σ D t h + ∇ • m Σ D t h v + ρ∂ α e + ∇ • σ ∇α ∇α , (43b)
C = - γ Σ - m Σ (D t h) 2 . ( 43c
)

Final form of the system

We define the pressure p of the two-phase medium and the partial pressures p k of each phase by where ρ is the mixture density defined as ρ = α 1 ρ 1 + α 1 ρ 2 , α k is the volume fraction of phase k = 1, 2, ρ k the partial density. Then by injecting relations (42a) into (41) one obtains the system

p = ρ 2 ∂e ∂ρ , p k = ρ 2 k ∂e ∂ρ k , (44)
∂ρ ∂t + ∇ • [ρv] = 0, (45a)
∂ρY ∂t + ∇ • [ρY v] = 0, (45b)
∂ρv ∂t + ∇ • (ρvv T ) + p + γ 2 I d (45c) + ∇ • σ ∇α∇α T ∇α - 1 2 σ ∇α I d + γ Σ I d = 0, D t α -- γ m Σ = 0, (45d)
D t 1 Σ + 1 √ -mγ p 2 -p 1 + ∇ • σ ∇α ∇α = 0. (45e)
System ( 45) is a generalization of the system found in [START_REF] Drui | A hierarchy of simple hyperbolic two-fluid models for bubbly flows[END_REF][START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF] and degenerates towards it when considering the interfacial area density as a function of the volume fraction only.

In the momentum equation (45d), the terms function of the volume fraction gradient are common terms found in the literature [START_REF] Blanchard | A large scale multi-fluid/dispersed phase approach for spray generation in aeronautical fuel injectors, proceeeding[END_REF]. Equation (45e) is the transport equation of the fluctuating interfacial density area. In a steady state regime, Equation (45e) yields the classic Poisson equation

p 2 -p 1 + ∇ • σ ∇α ∇α = 0 (46)
It is important to notice that System (45) is valid for any flow topology as opposed to the system found in [START_REF] Drui | A hierarchy of simple hyperbolic two-fluid models for bubbly flows[END_REF][START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF] only valid for dispersed flow.

Conclusions

In the present work, we a definition of a two-phase interface by means of a PDF that departs from the classic geometrical definition. Then we derived two-phase flow model accounting for small and large scale separation of the interface description by means of supplementary convected geometric variables.

The Least Action Principle yields a model that accounts for two-scale kinematics and two-scale surface tension through the introduction of density area flucutations. This is an extension of previous work that was able to account for a monodisperse spherical bubbly flow at the small scale [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Cordesse | Derivation of a two-phase flow model with twoscale kinematics and surface tension by means of variational calculus[END_REF].

In future works, we will try to encompass new effects of the interface dynamics such as stretching through the evolution of the mean curvature that was assumed constant in the present work.

Remark 1 Figure 1 :

 11 Figure 1: Large and small scale interface dynamics

Figure 2 :

 2 Figure 2: Filtering of the volume fraction using the spatial averaging on the ball B c , α = 0 α = 1.

  10 th International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019

th International Conference on Multiphase Flow, ICMF 2019, Rio de Janeiro, Brazil, May 19 -24, 2019