Mikael Mr

Dumortier

Mikael Dumortier
email: mikael.dumortier@iemm.univ-montp2.fr

Mikael Dumortier
email: mikael.dumortier@gmail.com

José Sanchez

Michel Keddam

Olivier Lacroix

Mikael Dumortier

Dear Mr Veziroglu

Mikael Dumortier

Mikaël Dumortier

Meng Ni
email: bsmengni@polyu.edu.hk

M M Hussain
email: mmhussain@uwaterloo.ca

Anchasa Pramuanjaroenkija
email: anchasa@gmail.com

Theoretical considerations on the modelling of transport in a three-phase electrode and application to a proton conducting solid oxide electrolysis cell

Keywords: Mathematical modelling, Electrolysis, Porous media, SOEC

This paper presents a numerical model for fuel cells and electrolysis cells that use cermets as electrodes. The continuity equations for mass and charge are established via a control volume method in order to take into account the variation of metal or ceramic fractions inside the cermet. The Dusty Gas model is used to describe diffusion inside the electrodes and the Brinkman equation is employed for the conservation of momentum inside the fluid phase. The model is tested for the protonconducting SOEC technology on a generic cell design and the calculations are performed with COMSOL Multiphysics 4.1™. A parametric analysis is carried out on a proton-conducting SOEC in galvanostatic mode in order to evaluate the influence of parameters on oxygen production across the electrode in the anodic compartment.

Nomenclature

A i

Gas reactant or product in the half-reaction at an electrode C T Total concentration (mol m -3 )

D i

Diffusion coefficient for species i (m 2 s -1 )

D ij

Diffusion coefficient of species i inside species j (m 2 s -1 )

D i,Kn

Knudsen diffusion coefficient for species i (m 2 s -1 ) Hydrogen is a promising energy carrier due to its non-polluting combustion which releases only water and a high amount of energy (120 MJ kg -1 ) [1]. However, hydrogen gas does not exist naturally on earth and therefore needs to be produced industrially. Today, the hydrogen economy is mainly dependent on hydrocarbons, as it is produced principally through steam reforming [2]. However, hydrogen production from hydrocarbons results in a concomitant production of large amounts of carbon dioxide. Several processes, such as high temperature water electrolysis (HTE), are being developed to obtain mass production of hydrogen independently from oil or to transform carbon dioxide into chemical products with high value. HTE installations can be coupled with existing power plants that release large amounts of heat, which is then used as energy input for the splitting of water [3,4]. HTE therefore requires less electric power than low temperature electrolysis.

d
One of the most promising technologies for HTE is the Solid Oxide Electrolysis Cell (SOEC) based on Solid Oxide Fuel Cell (SOFC) technology. The SOEC principle is shown schematically in Figure 1. In this study, the cell is an electrochemical membrane reactor composed of two channels separated by a membrane assembly made of a dense solid oxide layer positioned between porous electrodes. The electrodes are made of cermet, a porous composite material made of solid oxides and metals. In proton conducting cells, steam is fed at the anode and oxidized in the cermet. The reaction produces oxygen that is released into the channel, electrons that are conducted by the metal and protons that are transported into the perovskite. Having passed through the dense membrane, protons recombine with electrons at the cathode to produce pure hydrogen. Conversely, in anion-conducting cells, steam is reduced at the cathode, producing hydrogen and oxygen ions that recombine into oxygen at the anode.

As ions can only pass through dense solid membranes at relatively high temperatures, electrolysis takes place under hard operating conditions (temperature ranged between 500 °C and 1000 °C [5] and pressure from 1 to 100 bar).

As stated above, charge transport inside the solid electrolyte depends on the material that is chosen, whereas the global transport process in porous cermets is much more complex: metal allows electron transport, ions are transferred through the solid oxide phase and gases flow through open pores. The electrochemical reaction occurs at the junction of these three phases, called a triple-phase boundary.

These electrochemical membrane reactors are still at the research stage and their future development requires further efforts in numerical simulation to optimize the heat and mass transfers as well as electrolysis cell design.

Numerous mathematical models have been developed for fuel cells, especially for SOFCs [6] and PEMFCs [7,8]. Numerical modelling of electrolysis cells is a more recent topic, but the few studies that have been presented use the same approach as for fuel cells. Meng Ni et al [9] developed one of the first SOEC models in 2007; they built a 1D model for a SOEC by integrating a set of ODEs along the electrodes, with assumptions based on the results reported in the literature. This model features no differential equations and can easily be included in engineering numerical codes with user routines such as Fluent™ or Star-CD™. It is a very good compromise between precision and time calculation. In 2009 Grondin et al [10] proposed a 3D model including thermal effects, Knudsen diffusion and the generalized Butler-Volmer law. The equation set was solved using the commercial software COMSOL Multiphysics 4.1™. These models are based on regular assumptions when modelling SOFCs.

These frequently include the consideration of ideal gases [11]. However, those gases tend to deviate from the ideal gas behaviour at high pressure, especially steam. To assess the effects of this simplifying hypothesis, the compressibility factors of gases have been included in the equation of state used in this work.

In terms of electrokinetics, the Butler-Volmer equation is often used to describe the charge transfer at the junction of the three phases. Two forms of this equation are regularly employed in models: the classical Butler-Volmer equation [11], that accounts for the activation effects, and the generalized Butler-Volmer equation [10,[START_REF] Coeuret | Elements de génie électrochimique, Technique et documentation[END_REF], that accounts for both activation and concentration effects. The first equation is obtained from the second by considering that the concentrations of chemical species do not vary across the electrode.

Inside volumetric electrodes such as cermets, reaction sites are distributed within the entire volume and the electrochemical reaction can take place anywhere in the electrode. However, due to the different transport mechanisms inside the cermet, the reaction rate is higher in some regions of the electrode. Several studies have shown that the reaction rate is higher near the interface between the electrode and the electrolyte. Thus, the reaction domain is often treated as a boundary condition [9,10,[START_REF] Hussain | [END_REF]. This assumption is discussed in this paper regarding the choice of the Butler-Volmer equation.

Finally, some studies have considered electrodes in which the size of the particles and the volumetric fraction of each phase can change in depth [14,15]. This affects the properties of the medium. To ensure that the effects of such a composition within the model are rigorously taken into account, this paper demonstrates mathematically the local equations inside a three phase electrode, based on a control volume averaging approach [16].

The model aims to be sufficiently general to be applied for both SOECs and SOFCs.

Indeed, the fundamental equations describing the phenomena inside SOFCs or SOECs are the same, as one technology is the reverse of the other. From a numerical calculation point of view, only the boundary conditions differ between the two problems. This model is tested throughout a parametric analysis on oxygen production along the anode of a proton-conducting cell, involving every type of transport phenomenon. For this we carried out calculations coupling the electrochemical production of oxygen with the mass transport inside the anodic compartment taking in consideration the coupled contribution of diffusion and convection through the Peclet number (Pe).

2 System of equations

General assumptions

The following general assumptions were considered when designing the model:

• Steady state • Isothermal • Single-step electrochemical transfer
• No deformation of the structure due to pressure effects

• Spherical particles

• Compressibility coefficients dependent on pressure only

• Negligible effects of gravity (weight, natural convection)

• No interfacial ohmic resistance between grains

Electrode media description

The electrodes studied in this paper are cermets, which are considered as one of the most interesting candidates as electrode material for high-temperature steam electrolysis [17] . The cermet is a mixed media made of three phases: an electron-conducting phase, which is usually a metal, an ion-conducting phase, which is usually a ceramic, and a gas phase that allows the transport of reactants and reaction products. The metallic and ceramic phases are motionless and transport charge by migration. Conversely, the gas phases do not carry any charge and transport chemical species by diffusion and by convection. These transport phenomena inside the electrode are modelled mathematically by continuity equations obtained by applying the volume averaging method [16] to a representative volume of electrode. This volume comprises both a solid and a fluid phase and defines a mesoscopic volume V; this is large enough to smooth the morphological complexities and small enough to capture the global transport properties (i.e. typically a few pore lengths). The scale of the control volume must be smaller than the averaging length of the physical quantities in order to consider that they vary slowly inside the volume.

Normally the balance cannot be written directly on the whole control volume because it contains more than one phase. Therefore, the aim of the demonstration in this paper is to build local equations for the representative volume by writing a balance on the volume of each phase (see Figure 2). Each of these control volumes is bounded by a control surface S i , in which normal vector n r is outer-pointing. The volumetric fraction of each phase can vary inside the electrode. They are written as:

= ; = ; = (1)
with:

+ + = 1 (2)
Here V io , V e and V g are the volume of each phase in the representative volume V (see Figure 2)

ε g is the total porosity. These fractions can depend on the position inside the electrode.

Chemical reaction modelling

The half-reaction in the cermet can be written as:

+ + ∑ = 0 (3)
It involves an ion Io of charge k, gas products or gas reactants A i , and electrons, where ν e and ν i are the stoichiometric coefficients for electrons and the species i, respectively. For water oxidation for example:

→ 1 2 + 2 + 2 2 - = e ν , 1 + = k , 1 2 + = O H ν , 1 2 + = O ν
This reaction implies that there are charge and mass transfers between the three phases, which are assumed to take place at the triple-phase boundary. The geometrical nature of this zone of transfer is not well known. It can therefore be modelled as a point, a line, a surface or a volume. In calculations where the electrodes are considered as planes, this zone is logically expressed as a surface. In our case, however, it is easier to visualize the charge transfer within a reactive volume than in a reactive surface. Therefore the reactive zone V TP is expressed in cubic metres. It is assumed to be negligible compared with the volume of other phases and can be considered as a fourth phase inside the cermet. The volumetric fraction of V TP is called the specific reactive volume and is noted ε TP .

! = "# (4)

Transfer rate variables

The transfer is initiated by an external electric current or voltage. The rate of the local transfer depends on the kinetics of the charge conduction transport and gas transport competing inside the volumetric electrode. The transfer rate of the electrochemical reaction is represented locally by two variables:

• i r: the volumetric transfer current (in A m -3 ). This variable expresses the charge quantity that is transferred from one phase to another per second and per unit of volume. It is very useful for establishing a relationship between electric current and molar production.

• r r : the volumetric molar transfer rate (in mol s -1 m -3 ). This variable expresses the number of moles of electrons that are implied in the half reaction.

The two quantities are related by Faraday's constant F:

i r = F r r ( 5
$ % = $ % & ' ( ) * ( ) * + exp / 0 1 2 34 - ( 6 ( 6 + exp
The electric potential difference between the metallic phase and the ceramic phase is called the overpotential and is given by [18]:

3 = 7; -; & 9 -7; -; & 9 (7)
with:

; & = < = > + 1 ; ; & = 0 (8)
As a common assumption for fuel cells [19], the charge transfer coefficient α is considered equal to 0.5.

Some studies have considered other forms of the Butler-Volmer equation that do not feature pre-exponential terms containing out-of-equilibrium concentrations [6,20,21]. These forms are normally used when the activation overpotential and the concentration overpotential are calculated separately. In this study, the activation and concentration effects are both contained inside Equation (6). Other authors [START_REF] Hussain | [END_REF] have considered another expression of the overpotential that does not involve the equilibrium potential 0 e Φ . One argument for preferring the generalized form and the expression of the overpotential described by Equation ( 7) is to account for the electrochemical reaction having to slow down in the case of steam exhaustion, i.e. when x red ≈ 0. Moreover, when no electric current is applied, the generalized form is reduced to a Nernst equation form, which is not the case for the other forms:

; -; = ; & -; & + 2 1 ?@ A ( 6 ( ) * + ( 6 + ( ) * B (9)

Charge, mass and conservation

The method that is presented in this paragraph is used for the conservation of charge and mass. The generic form of a conservation equation for a scalar quantity is:

Stock + Flux + Source = 0 (10)
Stock is the temporal variation term, flux is the sum of fluxes of the quantity that enter or leave the control volume and source corresponds to the source of quantity inside the control volume. The method is first applied to the conservation of charge inside the ceramic phase and then extended to mass and charge inside the metallic phase. The demonstration consists of expressing the mass and charge transfer at a global scale with integral equalities; we then use the volume averaging method to approximate the integrals and obtain the continuity equations. As an example, we will detail the method for the conservation of charge inside the ceramic phase and then apply it to the conservation of charge and mass inside the other phases.

Integral form

The conservation of charge inside the ceramic phase is written as:

C CD ∭ F GH + ∬ J K GL @ MK = 0 N ( 11
)
where io ψ is the ionic charge density (C m -3 ) and J K the ionic charge flux density (A m -2 ).

This integral equation is expressed on the control volume of the ionic conductor only. Several geometrical operations are therefore needed to express a local equation on the control volume that includes all the physical phases. To express the flux term, the control surface of the phase is split into two surfaces, as shown in Figure 3:

S io = S io (int) + S io (ext) ( 12
)
S io (ext) = Sio ∩ S ( 13
)
S io (int) = Sio -(Sio ∩ S) ( 14
)
Neglecting the contact area between the conductors gives:

S g (int) = S e (int) + S io (int) (15)

Phase source terms

During the electrochemical reaction, charge is transferred between the ceramic phase and the metallic phase. On the whole, however, no charge is created inside the control volume, hence:

O = -O = -P % H ! (16)
χ e and χ io are the charge production terms (C.s -1 ) inside the metal and the ceramic respectively.

This production takes place inside the reactive volume V TP .

By neglecting the mass of electrons, the change in mass inside the gas is equal to the mass of the ions transferred from the ceramic phase to the gas phase. Hence:

Q = - R S T 1 = -U V $ % H ! ( 17
)
where φ g is the mass source (kg s -1 ) in the gas and M I the molar mass of the ion.

Flux term

Following the additivity of integration on intervals, the flux term can be written on each sub-surface:

∬ J K GL @ MK = ∬ J K GL @ MK N 7 (D9 + ∬ J K GL @ MK N 7 WD9 N (18)
The last term, integrated on S io (int), corresponds to the charge exiting the electronic conductor and entering the reactive volume because of the reaction. It can be treated as a source term:

∬ J K GL @ MK N 7 WD9 = -O (19)
The normal vector of the volume is outer-pointed, hence the minus sign. The term integrated on S io (ext) must be expressed on the control surface S using integration by substitution. Before that, another decomposition of the geometrical surfaces is needed.

The sum of the external surfaces is equal to the surface of the representative volume:

S = S g (ext) + S e (ext) + S io (ext) (20)
For each phase, the external surface fraction is defined as:

X = N 7 (D9 N (21)
Assuming that this fraction does not vary with a change in the control surface, the differential form of this equation gives:

GX = YN 7 (D9 YN (22)
Using integration by substitution, the expression of the integral on the whole external surface is given by:

∬ J K GL @ MK N 7 (D9 = ∯ J K X GL @ MK N (23)
This integral can be written on the representative volume by using Ostrogradsky-Gauss theorem [22]:

∯ J K X GL @ MK N = ∭ [. 7J K X 9 GH (24)
Ultimately, considering the definition of the representative volume (see §2. in the Cartesian base.

Internal variation term

Using integration by substitution:

C CD ∭ F GH = C CD ∭ F GH (26)
And again considering the definition of the representative volume we obtain:

C CD ∭ F GH ≈ C7a b 9 CD H (27)

Equations for mass and charge

Ultimately, the continuity equation of the ionic charge inside the representative volume is obtained by coupling equations ( 9), ( 19), ( 25) and ( 27) and dividing by V:

C7a b 9 CD + [. 7J K X 9 = P % ! (28)
The conservation of mass and the conservation of electronic charge inside the representative volume are built by following the same analysis. The continuity equation for electronic charge conservation is therefore written as:

C7a b 9 CD + [. 7J K X 9 = -P % ! ( 29
)
where ψ e is the electronic charge density (C m -3 ) and J K is the electronic charge flux density

(A m -2 ).
The electronic potentials are related to charge fluxes by Ohm's law:

For the ceramic phase:

J K = -c [; (30)
For the metallic phase:

J K = -c [; (31)
In Equations ( 30) and ( 31), the conductivity terms σ io and σ e reflect at the same time the intrinsic contribution of the materials and the non-explicit topography of the grain-to-grain contacts. No charge or mass is created inside the control volume, so the source terms are equal to zero. The next paragraph will show that source terms appear when dealing with the equation for phases.

The continuity equation for mass can then be written as:

C7d b 9 CD + [. ef g MK X h = -U V $ % ! ( 32
)
where ρ g is the density (kg.m -3 ) and g MK the velocity vector inside the fluid phase (m s -1 ).

Mass conservation for chemical species

The mathematical treatment for species is more complex. The analysis starts with the conservation of mass in the gas phase, which is given by:

C CD ∭ f GH + ∬ f g MK GL N @ MK = -U V $ % H ! (33)
The mass flux across the control surface is equal to the sum of the mass flux of the species:

∑ ∬ U i j MMMK GL @ MK N = ∬ f g MK GL N @ MK (34)
Using the same analysis as above for equations [START_REF] Darwent | National Standard Reference Data Series[END_REF] and [START_REF] Saul | [END_REF], the local equation for species is:

C7d b 9 CD + ∑ U [. eX i j MMMK h = -U V $ % ! (35)
The idea now is to express a continuity equation of mass for each species. The sum of individual volumetric molar productions is linked with the mass source term:

∑ U $ %7 9 = -U V $ % ! (36)
Here r r(i) is the volumetric molar production for species i inside the control volume (mol m -3 s - 1 ). Hence:

C7d b 9 CD + ∑ U k[. eX i j MMMK h -$ %7 9 l = 0 (37)
The stoichiometry of the reaction gives:

% )7 9 = % )7m9 m (38)
and the conservation of mass during the reaction gives:

∑ U n n + U V = 0 (39)
In this equation the mass of the electrons is ignored. Ultimately, using Equations ( 38) and (39), the molar production of the species i simplifies to:

$ %7 9 = $ % ! (40)
Assuming steady state and coupling, Equation (40) and Equation ( 37) give:

[. eX i j MMMK h = $ % ! (41)
Note: It is possible to obtain the equation in non-steady state by using the relationship between the density of the gas mixture and the concentration of the species presented in §2.7.

Diffusion model

The objective is to relate the concentration of species to their fluxes. As the model is applied to a porous media, the Knudsen diffusion regime must be envisaged if the pore radius is small enough. Moreover, the model must take into account the multi-component and nonideal nature of the mixture. The Dusty Gas Model (DGM) has been shown to be a relevant model for describing the diffusion of species inside a volumetric electrode [23]. The matricial form of the DGM has been reported by Krishna and Wesselingh [24]. Solving this system often requires the use of numerical methods; however, for a mixture composed of only two species, it is possible to reduce the DGM to a simpler form, close to Fick's equation:

i j MMMK = -o [p q + p . q . g MK (42)
where i j MMMK is the molar flux of the species i (in mol m -2 s -1 ), D i the diffusion coefficient of the species i (in m 2 s -1 ), x i the volumetric fraction of the species i and C T the total concentration of the mixture (mol m -3 ), g MK is the velocity vector in the mixture.

The non-ideality matrix used in the Dusty Gas model involves compressibility factors for gases. In this study, it is assumed that the compressibility factor for the species i has little dependence on its molar fraction. Hence the non-ideality matrix reduces to the identity matrix:

[Γ] = I d ( 43
)
By inverting the matrix system for two species, the expression of the diffusion coefficient for the species i is given by: o = o n vv is the effective binary diffusion coefficient between species i and j (m 2 s -1 ). It is calculated with the formula adapted from [START_REF] Bird | Transport Phenomena[END_REF] for non-ideal gases with SI units and compressibility factors:

o n vv = 8.41 × 10 z b { | m }/• S m €/• ! • m • ' m (46)
with

U n = 8 S 8 S m ⁄ ⁄ (47)
σ ij and ω ij are evaluated with the data and formulae presented in [START_REF] Bird | Transport Phenomena[END_REF] and [START_REF] Poling | The Properties of gases and liquids[END_REF]. Note that the mixture composed of steam and oxygen at the anode is supposed to be a polar gas, hence the calculations of its kinetic properties differ from those of hydrogen, which is non polar. o ,"W vv is the effective Knudsen diffusion coefficient for the species i (m 2 s -1 ). It is calculated with the formula presented in [24]:

o ,"W vv = Y … z b { † ‡2 ˆS ( 48
)
The tortuosity is evaluated with the RUC model [START_REF] Du Plessis | [END_REF]:

‰ = b e8 b he8 Š 8 b h b /8 e8 b h •/} 4 (49)

Momentum

The conservation of momentum for flows throughout porous media is very well described by the Darcy-Brinkman equation, expressed in our case [16].

Porous media

The different variables required to describe the geometric aspects of the porous media are taken from the work presented by Chen et al. in [31]. The surface fraction and the pore radius are given by:

X = b % ⁄ b % ⁄ b % ⁄ b % ⁄ ; $ = z / 8 b % ⁄ b % ⁄ 4 (51)
Following the analysis of Costamagna et al [32] on the application of percolation theory for SOFCs, an expression is built to estimate the value of the specific reactive volume.

ε TP = ε reaction P TP (52) P TP is the probability that the three-phase boundary is electrochemically reactive and reaction ε is the geometric extent of the reactive volume. The works in the literature usually define the reactive area as a line and calculate the triple boundary layer as a proportion of the contact perimeter between ceramic and metallic grains. Methods to calculate those quantities are given in [31].

In this study ε reaction is considered as a specific volume, extending from the junction line of the three phases to a characteristic "reaction length", l 0 , which represents the distance between grains within which the steam molecules can react:

% ••D W = z ' W7% ,% 9 % } ? & ( 53
)
This length is not known theoretically and is evaluated as l 0 = 10 -8 m, which is approximately a hundred times the length of molecular bonds [START_REF] Darwent | National Standard Reference Data Series[END_REF]. The value of this length or the choice of the nature of the triple boundary zone has no impact on local production of mass or charge on the control volume.

Properties of the gas mixture

Compressibility factors are obtained via a linear regression from the data taken from [START_REF] Saul | [END_REF][35][36][37]:

Z H2O = -4.90×10 -9 P + 1 (54) Z H2O = 2.23×10 -9 P + 1 (55) Z O2 = 2.82×10 -9 P + 1 (56)
The compressibility factor and the molar volume are related to the properties of the chemical species using Kay's rule [START_REF] Poling | The Properties of gases and liquids[END_REF]:

H '' = ∑ p H ' (57)
' ' = ∑ p ' (58)
The conservation of mass then gives:

U ' = ∑ p U (59)
Coupling the ideal gas adapted for real gases with the previous equations gives:

f ' = " " S " | " (60)
Given the lack of efficient models to calculate the viscosity of polar gas at high pressure [START_REF] Poling | The Properties of gases and liquids[END_REF], the viscosity of the mixture is assumed to be the viscosity of the predominant fluid species in the electrode.

Range of validity of the model

This model has a wide range of applications: it can be used for any type of SOFC or SOEC that respects the assumptions presented in §2.1. The model is tested on a protonconducting SOEC whose principle is shown in Figure 1.

Numerical model

The planar SOEC is built in 2D. Indeed, the hydraulic regime considered is laminar;

hence the 3D effects due to turbulence are ignored in the z direction.

Geometry

The electrolyser geometry considered here is presented in Figure 4. The channels are 0.4 mm wide, the electrodes are 0.15 mm wide and the electrolyte is 0.1 mm thick. The channels are longer than the electrode to observe the mixing along the channel. This geometry is not intended to represent a real SOEC; it is designed in a simple manner to focus on the physical phenomena. The model has been solved using COMSOL Multiphysics 4.1™. The species balance equation and the diffusion equation is solved using the "Species transport in porous media" node. Mass and momentum conservation are solved with the "Free and porous media flow" node. Charge conservation is solved using the Poisson's equation in the "PDE node".

As COMSOL Multiphysics 4.1™ does not allow the introduction of boundary conditions at the electrode/electrolyte interface for the ionic and the electronic potential, the electrodes and the electrolyte are modelled as a single domain for the charge equation. In this domain, the volume fractions vary as a function of space. Therefore, the electrolyte is defined where ε io = 1, the channel is defined where ε g = 1, and the electrodes are defined where ε e ≠ 0.

Instead of declaring multiple domains with multiple sets of equations, the charge balance is solved with one set of equations all over the domain. The counterpart of such a strategy is that the mesh must be very refined where the volume fractions vary dramatically, i.e. at electrode interfaces. In the electrolyte, TP ε cannot be calculated and must be specified. Any value can fit, as the exchange current density and thus the mass production term is equal to zero.

Numerically speaking, however, the value of TP ε must be kept at a very low value to avoid oscillations of the electric potential at the boundary between the electrodes and the electrolyte.

The reactor dimensions are presented in Table 1.

Boundary conditions

The boundary conditions are presented graphically in Figure 4 and specified in Table 2.

For the calculations, inlet velocity and pressure are considered as the same at the cathode channel and at the anode channel.

The calculation is performed in galvanostatic mode and the global reaction rate is imposed by the value of the electric current; this is equal to the product of the current density input by the surface to which it is applied. In addition, the inlet velocity is calculated in order to avoid steam depletion, which means that there will always be enough steam available for the electrochemical reaction. The minimum mass flux needed to avoid steam depletion at electrodes is given by:

•--• ˜= ™ S š • › ( š • › + V 1 š • › ™ (61)

Mesh

Various meshes have been used for the calculation. All of them are refined near the boundaries where the concentration and hydraulic gradients are high. The mesh must also be more refined inside the electrodes, as the distribution of reaction kinetics along them is unknown. One of the meshes used contains 20696 elements with an average element of 0.9 (0.96 for the triangular elements only).

Reference calculation

Physical constant and other parameters presented in Table 3 are kept constant in the calculations. The values of the physical quantities whose variation are studied are compiled in Table 4.

Results and discussion

The mass production of a chemical species is imposed by the inlet electric current:

•-= U V 1 (62)
This equation shows that in galvanostatic mode, the electric current is the only parameter influencing mass production. However other parameters such as particle radius, metal fraction and inlet oxygen molar fraction will have an influence on the distribution of this production across the electrode. The relevant variable to study the influence of parameters on this distribution is the product of the specific reactive volume ε TP by the transfer current density i r .

Indeed, Equations (5) (40) and (63) show that this product is related both to the local mass production and to the inlet electrical current I. The integral of this product across the cermet is fixed by the electrical current but its local value depends on the parameter values.

= oe

! P % GH • %' D • )" ž (63)
This product has therefore been chosen as the most suitable variable for a parametric analysis to evaluate the importance of parameters on the distribution of oxygen production across the anode. This analysis has been carried out along the anode, where every type of transport phenomenon is involved, unlike the cathode, where there is no multicomponent gas transport.

It can be used in electrode design as a guide to choose the properties of cermets and optimize the yield of the electrolyser or the durability of electrodes. The analysis started with the electrochemical parameters (form of Butler-Volmer equation, electrical conductivity); we then studied the influence of structural parameters (volume fractions and particle radius), ending with the process operating parameters.

Form of the Butler-Volmer equation

When using the Butler-Volmer law with no concentration factor, the oxygen production is close to the electrolyte/electrode interface. In this case, the energy loss of the system is minimal if the protons are conveyed over the shortest distance possible, as the ionic conductivity of ceramics is lower than the electrical conductivity of metal. When using the generalised Butler-Volmer equation, the concentration effects compete with the ohmic drop effects. The oxygen production rate along the anode becomes oxygen concentration dependent, which makes the problem far more complex. The fact that the distribution of oxygen production along the anode moves to the electrode/channel interface (see Figure 5) is explained by the fact that the difference between oxygen and steam concentration decreases as steam penetrates deeper into the electrode. As the two profiles are different, the assumptions concerning pre-exponential terms containing out-of-equilibrium concentrations in the Butler-Volmer equation should be justified by comparing concentration with electrical potential effects.

Effect of conductivity

Figure 6 shows that the distribution of oxygen production inside the electrode is less dependent on the value of the conductivities than on the value of their ratio. Indeed, the system tends to distribute oxygen production over the anode in order to minimize energy losses caused by Joule's effect.

• When σ e > σ io , Joule's effect is mainly caused by the conduction of protons inside the ceramic. The energy loss in the system is minimal at the shortest distance covered by protons,i.e. when they are produced near the electrode/electrolyte interface.

• When σ e < σ io , Joule's effect is mainly caused by the conduction of electrons inside the metal. The energy loss in the system is minimal at the shortest distance covered by electrons, i.e. when they are produced near the electrode/channel interface.

• When σ e = σ io , the distribution of oxygen production along the electrode is not symmetrical because of the concentration effects implied by the generalized Butler-Volmer law.

Effect of metal volume fraction

Metal and ceramic volume fractions, along with particle radius are involved in the definition of the specific reactive volume ε TP and the specific metal surface fraction s e . As a result, a change in the fractions has an impact on percolation and on the reaction rate. In this calculation, the porosity ε g remains constant, equal to 1/3. Also, the metal volume fraction ε e varies across the electrode according to different trends (see Figure 7). In each case, the mean value of the metal volume fraction ε e is equal to 1/3.

Figure 8 shows that the specific reactive volume ε TP does not vary monotonously at the same time as the metal fraction ε e , which is consistent with the models taken with the percolation theory for P TP in §2.7. When a critical volume fraction of metal is reached, the chain of metal grains or ceramic grains that allow the conduction of charges is interrupted.

The electric or the ionic current can no longer be conveyed and the specific reactive volume ε TP becomes equal to zero. For example, oxygen production for the metal fraction profile called decreasing + shows three regions (Figure 9). Near the channel, with a high metal fraction, electrons are conducted in the metal phase, but there is no reaction because the ceramic fraction is too low to reach the ceramic percolation threshold. In the middle of cermet, both metal and ceramic percolate, electrons and protons can be conducted in their respective phases and the electrochemical reaction occurs. Near the electrolyte, with a low metal fraction, electrons are not conducted through metal as it no longer percolates, but the protons are conducted to the electrolyte as the ceramic percolates. Hence, along with the knowledge of percolation threshold, the metal fraction or the ceramic fraction can be adjusted in order to control the reaction rate inside the electrode. This limits significant and localized heat release or electrode assembly delaminating due to peaks of oxygen production near the interfaces.

Effect of grain radius

Grain radius has a significant impact on the properties of the porous medium, such as the specific reactive volume ε TP , the surface fraction, grain coordination and pore diameter. In this study, the radius of ceramic grain is fixed and the radius of metallic grains is a multiple of the latter. Figure 10 shows that when this ratio increases, the distribution of oxygen production along the anode flattens.

Indeed, when the radius increases the specific surface of metal s e decreases. Thus, the area of metal available for conduction of electrons gets lower. Because the electric current is fixed, both the electronic current density and the ohmic losses increase. As the effective conductivity of metal gets lower, the system transfers the charge quicker from the metal to the ceramic. This explains why the distribution of oxygen production flattens when the radius of metal grain increases.

Effect of inlet oxygen volume fraction

It is important to study the impact of inlet oxygen volume fraction. Indeed, producing steam by vaporization of liquid water is costly in terms of energy. Instead of feeding the electrolyser with pure steam, it could be possible to recycle the steam that has not reacted and which is enriched in oxygen (see Figure 11).

Figure 12 shows that increasing the inlet oxygen volume ratio lowers oxygen production in the centre of the electrode and increases oxygen production at the sides of the electrodes. There are two reasons for this phenomenon: firstly, when the oxygen fraction increases, the difference of concentration that drives the reaction decreases and concentration effects are therefore less significant than electrical effects. That is why oxygen production becomes higher on the electrolyte side. Secondly, under the operating conditions used for this calculation, the diffusion coefficient of oxygen is lower than steam. Hence the transport of gases is more difficult inside the electrode and consequently the system loses less energy if the reaction occurs without crossing the electrode. That is why oxygen production becomes higher on the channel side. This calculation was performed with an inlet velocity of 0.001 m s -1 in order to avoid steam depletion.

Effect of operating pressure

Operating pressures from 10 bar to 30 bar have been tested in the calculations.

However, operating pressure does not have a major influence on the distribution of oxygen production across the anode.

Effect of inlet velocity

The parametric calculation showed that the value of velocity does not have an important influence on the distribution of oxygen production across the anode for the values that have been tested (0.0001 m s -1 to 0.1 m s -1 ).

However, depending on the value of the Peclet number inside the anodic channel, the oxygen that is released out from the porous anode can flow back to the inlet or be carried out by the water vapour flow (see Figure 13). The expression of the Peclet number is given by:

OE = Ÿ + r › • (64)
If Pe < 0.02, then diffusion transport overcomes convection transport and oxygen flows back to the inlet. On the contrary, if Pe > 0.02, then convection transport overcomes diffusion transport and oxygen is carried by the flow to the outlet. As a result, the value of velocity must be high enough not only to avoid depletion, but also to avoid this counter-diffusion effect. The condition at the inlet velocity is therefore given by:

¡ & ≥ •£p A™ S š • › ( š • › + d š • › V 1 š • › ™ , &.& r › • B (65)

Conclusion

This paper has presented a series of studies that helps to predict a wide range of transport and transfer phenomena inside an anion-conducting SOEC and a proton-conducting SOEC.

• A set of continuity equations was rigorously demonstrated for the conservation of quantities and for the DGM applied for two species. This set of equations also takes into account the non-idealities of gases under strong operating conditions and the concentration coefficient inside the Butler-Volmer equation.

• A numerical model considering the electrodes and the electrolyte as a single reacting medium in which the metal and the ceramic fractions vary was successfully implemented in the commercial software COMSOL Multiphysics 4.1™.

• A parametric analysis showed that the electrochemical reaction can be distributed all over the anode. Hence the thin layer assumption is not always valid and must be checked numerically. A parametric analysis showed that the parameters of the process have an influence on the distribution of oxygen production along the anode and that electrode structural parameters such as metal fraction or grain radius have a major impact on the oxygen production profile across the anode.

The calculations also showed that the necessity to include the concentration coefficients inside the Butler-Volmer equation and that inlet velocity is an important parameter to be monitored. It must be high enough to avoid steam depletion and to avoid the diffusion of produced oxygen back to the inlet. Three points must be mentioned: firstly, the interaction between the ionic conductivity of ceramics and pressure is not yet taken into account in the model. Secondly, a parameter without influence on the distribution of oxygen production can be influential on the overpotential or on other phenomena like velocity on the diffusion of oxygen in the channel. It is thus important to determine the significance of the quantity that is monitored in the problem. Thirdly, no thermal effects are considered for now. Their consequences on transport phenomena will be assessed in a future work.

pFrr 0 ∆

 0 Abbreviations

∇

  is the Del operator. It is written [

For x j ≈ 1

 1 (j being the solvent), this expression reduces to Bosanquet's formula, which is often used when the Knudsen diffusion is taken into account[23,25,26]:

  $ % ! 4 g MK (50)This equation is already implemented in COMSOL Multiphysics 4.1™. The permeability κ is evaluated via the Kozeny-Carman equation[30].

Figure 1 .

 1 Figure 1. Principle of proton conducting and anion conducting SOECs.

Figure 2 .

 2 Figure 2. Control volume as a sum of single phase volumes.

Figure 3 .

 3 Figure 3. Internal and external surface of the ionic conductor.

Figure 4 .

 4 Figure 4. Electrolyser regions and boundaries.

Figure 5 .

 5 Figure 5. Effect of concentration coefficients in Butler-Volmer equation on oxygen

Figure 6 .

 6 Figure 6. Effect of conductivity ratios on oxygen production across the anodic cermet. X-

Figure 7 .

 7 Figure 7. Presentation of the different metal volume fractions along the anodic cermet tested.

Figure 8 .

 8 Figure 8. Profile of specific reactive volume along the anodic cermet respect to the metal

Figure 9 .

 9 Figure 9. Effect of metal fraction on oxygen production along the anodic cermet. X-axis:

Figure 10 .

 10 Figure10. Effect of radius ratio on oxygen production along the anodic cermet. X-axis:

Figure 11 .FigureFigure 12 .

 1112 Figure 11. Configuration for steam recycling

Figure 13 .

 13 Figure 13. Oxygen flux vectors and Peclet number along the anode showing counter diffusion of oxygen back to the inlet which is located at the ordinate Y = -9×10 -4 m.

Acknowledgements

This work was supported by the Agence Nationale de la Recherche.


Mouna Baklouti
email: mouna.baklouti@enis.rnu.tn

Philippe Marquet
email: philippe.marquet@lifl.fr

Jean-Luc L Dekeyser
email: jean-luc.dekeyser@univ-lille1.fr

Mohamed Abid
email: mohamed.abid@enis.rnu.tn

FPGA-based many-core System-on-Chip design

Keywords: Field Programmable Gate Arrays, Intellectual Property, Single Instruction Multiple Data, System-on-Chip, intensive signal processing

published or not. The documents may come

Introduction

Modern embedded systems tend to be more and more sophisticated with the integration of multiple functionalities in the same system, often implemented on a single chip, called System-on-Chip (SoC). Adding to that, the wide spread of data-intensive applications, such as multimedia applications, medical imaging, numerical filtering, radar or sonar signal processing, etc., requires powerful architectures with higher execution performances. With the huge number of transistors available in today chips, and the stagnation of clock frequencies due to power dissipation issues, chips with multiple processor cores are becoming more commonplace, in particular Single Instruction Multiple Data (SIMD) on-chip architectures with intensive parallel computations are possible. SIMD architectures are widely recognized as being well-suited for media-centric applications like image and video processing applications [START_REF] Gentile | Impact of Pixel per Processor Ratio on Embedded SIMD Architectures[END_REF] because they can efficiently exploit massive data parallelism available with minimal energy [1]. As most multimedia algorithms operate on relatively independent groups of data, such as sound frame samples, rows of pixels, or video frames, it is often easier to use a data-oriented perspective to parallelize these algorithms. SIMD on-chip systems offer many advantages. The SIMD processing presents a computational efficiency originating from the reduced overhead of control operations. The cost of instruction and address decoding is amortized over the many processing elements [START_REF] Raab | Low Power Design of the X-GOLD SDR 20 Baseband Processor, Proceedings of the Design[END_REF].

Adding to that, the SIMD parallel architecture enables higher performance for intensive data-parallel applications at lower power consumption [START_REF] Abbo | A programmable smart camera architecture[END_REF][START_REF] Dulong | Recherche d'une exploitation nergtique optimale des ressources de calcul dans un systme de vision sur puce[END_REF].

The design of these high-performance embedded systems for signal processing applications is facing the challenges of not only increased computational demands but also increased demands for adaptability to different requirements for these applications.

Various SIMD on-chip implementations have been proposed. Most of them are dedicated solutions for a specific application (many processors with short memory/small amount of processor with large memory/a given interconnection network). This normally results in good performance for the targeted application; however the performance of other applications may not be so good due to the diversity of parallel application requirements. Their design approaches also necessitate a long development time if some optimizations or modifications need to be made.

From the above observations, there is clearly a need of a generic SIMD prototype that can satisfy different application requirements. The proposed design methodology is based on IP assembly to easily and rapidly imple-ment different SIMD configurations. The provided parallel architecture is programmable, parametric, extensible and configurable; presenting high performances for data-parallel algorithms. FPGA (Field Programmable Gate Arrays) are targeted devices to implement the proposed design. Compared to ASIC (Application Specific Integrated Circuit), FPGA requires much less implementation costs and offers more flexibility. FPGA can also be re-wired and remotely reconfigured at any time. The inflexible and costly ASICs give FPGA-based solutions an upper hand in terms of implementation flexibility and cost effectiveness. This paper presents an SIMD massively parallel processing System-on-Chip. This system is an FPGA IP-based programmable system. Its architecture is flexible and can be customized to satisfy the requirements of data-intensive parallel applications. This opens a rich design space to the user and allows greater area/performance tradeoffs. We propose a hierarchical design implemented at Register Transfer Level (RTL) using the VHDL (VHSIC Hardware Description Language) language.

The remainder of this paper is organized as follows. The next Section presents some significant work related to the design issues of SIMD parallel on-chip architectures. It also deals with some proposed IP-based design frameworks. Section 3 introduces the SIMD massively parallel platform model. Section 4 describes the IP assembly approach. Section 5 highlights the data-parallel programming. In Section 6, performance evaluation is given through prototype implementation and video processing computations. Section 7 discusses the efficiency of the proposed SIMD on-chip system and gives comparisons with other implementations. Finally, Section 8 concludes this paper with a brief outlook on future work.

Related work

SIMD architectures have been proposed, implemented, and thoroughly studied for almost 50 years [START_REF] Lokhmotov | Programming and compiling for embedded SIMD architectures[END_REF]. Early SIMD machines were designed to provide supercomputing levels of performance [START_REF] Bouknight | The illiac IV system[END_REF][START_REF] Dorband | The Massively Parallel Processor: Architecture and Application[END_REF][START_REF] Fountain | The CLIP7A image processor[END_REF][START_REF] Kim | MasPar MP-2 PE chip: A totally cool hot chip[END_REF][START_REF] Tucker | Architecture and Applications of the Connection Machine[END_REF]. However, the end of the 90s saw the decline of SIMD machines due to their high fabrication cost, the impossible instruction broadcast at a high clock frequency and their scalability limits in terms of clock distribution and synchronous communications. Recently, the increase in the integration density and the use of novel design approaches as IP reuse [START_REF] Rowson | Interface-Based Design[END_REF] make possible the implementation of complex systems with reduced cost. Advances in technology have led to a wide range of SIMD systems: from SIMD arrays integrated into memory chips [START_REF] Yamashita | A 3.84 GIPS integrated memory array processor with 64 processing elements and a 2-MB SRAM[END_REF] to SIMD arrays exploiting sub-word parallelism in multimedia instruction set enhancements [START_REF] Nguyen | Exploiting SIMD parallelism in DSP and multimedia algorithms using the AltiVec technology[END_REF][START_REF] Peleg | Intel MMX for multimedia PCs[END_REF] (Intel SSE, AMD 3dNOW, Motorola AltiVec).

In this work, we focus on FPGA-based massively parallel architectures. Some proposed examples are presented in the next paragraphs.

A scalable streaming-array (SSA) architecture is presented in [START_REF] Sano | Scalable Streaming-Array of Simple Soft-Processors for Stencil Computations with Constant Memory-Bandwidth[END_REF]. SSA is an unidirectional linear array of pipelined-stage modules (PSMs), implemented on multiple FPGAs. Each PSM, implemented as an independent SIMD processing unit, consists of programmable simple PEs connected with a bi-directional 1D torus network. The PE's implementation is tailored to perform stencil computations. This work proposes a dedicated architecture for high-performance stencil computations on FPGAs. Moreover, a multichip solution would mean higher developmental effort and a larger power expenditure by the system compared to single-chip. In [START_REF] Ang | High Speed Video Processing Using Fine-Grained Processing on FPGA Platform[END_REF], an FPGA-based pixel array processor is presented. The hardware architecture comprises of primitive pixel processors that use bit-serial arithmetic to compute. Each processor processes a single pixel and is connected in a 2-dimensional mesh topology to form the overall array processor. The carried implementation is dedicated to perform Laplacian filtering on a 40 by 40 pixel gray scale video.

Tanabe et al. [START_REF] Tanabe | A study of an FPGA based flexible SIMD processor[END_REF] propose an FPGA-based SIMD processor. An opencore processor that belongs to SuperH processor series was modified to include more number of computation units and then performs SIMD processing. A JPEG case-study was presented however no idea about programming was given. The achieved clock frequency of the SIMD processor is too low (55 MHz) while many soft CPUs run faster these days. An SIMD architecture with 31 processing elements (PEs) is implemented on Xilinx FPGA [START_REF] Sayed | Towards an H.264/AVC HW/SW Integrated Solution: An Efficient VBSME Architecture[END_REF]. It is dedicated to run a motion estimation application. Another FPGA-based SIMD architecture with 30 PEs is proposed [START_REF] Rosas | SIMD Architecture for Image Segmentation using Sobel Operators Implemented in FPGA Technology[END_REF]. It is specific for edge detection performing the convolution of a mask over an image. An SIMD programmable processor dedicated to image processing sensors is implemented on FPGA, based on a self designed soft IP processor cell PicoBlaze [START_REF] Schurz | A Programmable Parallel Processor Architecture in FPGA for Image Processing Sensors[END_REF]. The 1-bit PEs can only execute simple computations on binary images. These work remain appropriate solutions. While these application specific systems deliver good performance and efficiency, they lack the programmability and versatility required to support the changing standards of multimedia. Their design is also time consuming, and costly process, increasing the overall production costs. Although some implementations are based on IP design, the designers do not always find the suitable IPs that match the application. Some IPs have excess performance and some do not have enough performance.

Many researchers propose soft vector processor architectures to accelerate data parallel applications. Yiannacouras et. al. [START_REF] Yiannacouras | Fine-Grain Performance Scaling of Soft Vector Processors[END_REF] described an implementation of the VESPA FPGA-based vector processor. The proposed processor can be customized in terms of vector chaining and vector lanes. It consists of a scalar MIPS-based processor coupled with a parameterized vector coprocessor based on the VIRAM vector instruction set. In [START_REF] Christopher | VEGAS: Soft Vector Processor with Scratchpad Memory[END_REF] a new soft vector architecture, named VEGAS, is presented. VEGAS consists of a standard NiosII/f processor, a vector core, a DDR2 controller plus external memory, and a DMA engine. Vector and DMA instructions are implemented as custom instructions in the regular Nios II instruction stream. For programming, the programmer needs to manually manage vector address registers, which is error-prone. In [START_REF] Severance | VENICE: A Compact Vector Processor for FPGA Applications[END_REF], a soft vector processor for FPGA applications, named VENICE, is presented. VENICE can be used to accelerate tasks that fit the SIMD programming model. It is implemented as an accelerator to the Nios II fast processor. In this case, the Nios executes all control flow instructions and issues instructions to the VENICE vector core as one or two tandem Nios custom instructions. The drawback of the VENICE implementation is its dependency on the Nios processor and hence can only be exploited in Altera FPGAs. The design isn't so flexible to fit various processor architectures. Generally, the disadvantages of vector processing is its inefficiency when dealing with irregular parallelism and the memory can easily become a bottleneck especially if data is not mapped appropriately to memory banks.

Our work considers SIMD on-chip implementation more generally and thoroughly explores a large design space. The proposed massively parallel SoC extends these work since it provides a configurable and parametric architecture in order to be suited to a wide range of data-parallel applications.

The IP-based methodology followed to generate one configuration provides a short development time and enough flexibility. An IP library is provided with the RTL design to help the designer selecting the suitable IP that delivers the highest performance for his targeted application. The proposed implementation is parametric and hierarchical making any modification and regeneration of the hardware parallel configuration easy. This design methodology provides lots of flexibility and enables to explore SIMD architectures. In the next section, the basic SIMD SoC platform model is described.

SIMD Parallel SoC Model

In this section, the basic SIMD parallel SoC model is exposed. The work deals with the very typical and simple SIMD model inspired from traditional SIMD systems mainly the famous MasPar [START_REF] Blank | The MasPar MP-1 architecture[END_REF]. The parallel system is designed for FPGA prototyping. Based on FPGA programmability and by using replication of IPs effectively, massively parallel architectures can achieve high performance at low cost. The proposed system supports the following features:

• a main processor, the Array Controller Unit (ACU), connected to instruction ACUIns and sequential data memory ACUData. It synchronously controls the whole system.

• a parametric set of elementary processing elements (PEs), each one connected to a local data memory PEM (each PE has its own private memory). All PEs perform parallel computations.

• ACU/PE interface: this interface is devoted to send control/execution orders to all PEs. It is modeled as a bus connecting the ACU with all PEs.

• a configurable neighboring interconnection network to connect PEs with their neighbors.

• a Network-on-Chip (NoC), working as a global router to perform pointto-point communications through different types of connections.

• PE activity mechanism: each PE executes the parallel instruction if its activity bit is set to '1'. This bit is controlled by the ACU via control instructions. This is especially useful in the case of conditional instructions.

• OR Tree: it computes logic global OR of all PE activity bits to test the status of PEs. It allows the ACU to know if at least one PE is active. The following section details the parallel SoC implementation.

IP-based massively parallel SIMD design

The proposed massively parallel RTL implementation is hierarchical and parametric, allowing the designer to build different SIMD SoC hardware architectures, in order to select the best configuration with the highest performances for a given data-parallel application. This implementation is based on IP assembly. For this purpose, an IP library is provided with the design, containing different IPs that can be used to generate different SIMD configurations.

IP-Library

The following paragraphs detail the provided IPs. While some IP construction, such as memories, does not require much attention, some others, such as processors, are much more complex. The following paragraphs deal with particularities of each IP block and its integration into the SIMD SoC. To build ACU and PE processors, two methodologies are proposed. The first one, called the processor reduction methodology, is based on refining the main processor in order to obtain a small reduced one. In the second methodology, called the processor replication methodology, the PE is chosen to be the same processor as the ACU to reduce the design time and facilitate the architecture building.

Processor reduction methodology

The PE is a reduced processor derived from the same ACU processor and only responsible of executing instructions. It receives decoded instruction from the ACU, thus its instruction memory is eliminated. The reduction simplifies the PE by eliminating its instruction logic and instruction memory, and thus saves millions of gates and hundreds of megabytes of memory in the overall system. As example, Table 1 compares between the processor miniMIPS and its reduced version in terms of consumed logics.

It is clearly shown that the reduced miniMIPS occupies smaller area (approximately the half surface) than the miniMIPS. Thus, integrating reduced processors allows to gain surface in order to implement a massively parallel system. The connection between ACU and PEs is a bus of micro-instructions delivered from the ACU decode stage to the input of the PE execute stage [START_REF] Baklouti | A design and an implementation of a parallel based SIMD architecture for SoC on FPGA[END_REF]. The reduction methodology was tested on the miniMIPS processor (Fig 3) [START_REF] Opencores | minimips overview[END_REF] and on the OpenRisc processor [START_REF] Opencores | OpenRISC processor[END_REF][START_REF] Dammak | Soft-core reduction methodology for SIMD architecture: OPENRISC case study[END_REF]. It enables to construct smaller and less complex PE based on IP reuse, and so integrate a large number of PE on a single chip.

Processor replication methodology

In this methodology, ACU and PE are similar (the same processor IP) in order to facilitate and accelerate their design. Compared to the reduction methodology, each PE receives parallel instructions from the ACU (Fig 4).

The replication methodology offers a large gain in the development time.

However, we are unable in this case to integrate a large number of PEs in a single chip. In this case, we have to choose a smaller processor that can be fitted in large quantities into the FPGA.

Network IPs

The massively parallel system integrates two communication networks:

a regular one to connect each PE to its neighbors, and an irregular one to assure point-to-point communications [START_REF] Baklouti | Reconfigurable Communication Networks in a Parametric SIMD Parallel System on Chip[END_REF]. The designer can use none, one or both routers to build the needed parallel configuration depending on the application requirements. The communications are managed through communication instructions that will be described later.

Global network

The global network is designed to perform any communication pattern be-

tween

Data-Parallel programming

The SIMD SoC is programmed by a single instruction stream partitioned into sequential and parallel instructions. It is programmed in a manner similar to inline assembly in C. However, C macros are used to simplify programming and make instructions look like C functions without any run time overhead.

A sequential instruction is carried out by the ACU as in the usual sequential architecture. A parallel instruction is executed in a synchronous manner by all active PEs, each PE taking its operands from its local memory and storing the result in this same memory (or may be in its own local registers). The massively parallel system instruction set is derived from the processor IP's instruction set used in the design. We distinguish different instructions classified into five main groups: arithmetic (add/sub/mul/div), memory access, communications, jump and branch, and system control. The first three types are executed both by ACU and PE, whereas the two others are only executed by the ACU. The arithmetic and memory instructions are duplicated for the PE. So that, in the instruction set we distinguish for example between a sequential addition and a parallel addition. This is accomplished through implementing the parallel instruction in the processor instruction table while modifying its decoding. Communication instructions are encoded from the processor instructions mainly load LW and store SW instructions. Other system control instructions are also added and rely on LW and SW instructions too. These instructions consist on storing/loading data to/from special registers or memory addresses to accomplish their functions. Examples of these instructions control the NoC's switching mode, the PE's activity state, etc. Table 2 details these instructions' coding.

Experiments: video processing-based algorithms

In this section, different massively parallel on-chip configurations are designed and tested to run video processing-based algorithms. Performance p addi r1,r1,adr l.addi r1,r1,adr p LW reg,0(r1) l.lwz reg,0x0(r1) P GET STATUS (reg,ident) read activity bit of the PE identified p lui r1,0x9 l.movhi r1,0x9 by "ident".

p ori r1,r1,0 l.addi r1,r1,ident p addi r1,r1,ident l.lwz reg,0x0(r1) p LW reg,0(r1) P SET STATUS (val,ident) modify activity bit of the PE p lui r1,0x9 l.movhi r1,0x9 (identified by "ident") p ori r1,r1,0 l.addi r1,r1,ident by "val" value.

p addi r1,r1,ident l.addi r2,r0,val p addi r2,r0,val l.sw 0x0(r1),r2 p SW r2,0(r1) P GET IDENT (reg) get identity p lui r1,0x2 l.movhi r1,0x2 p ori r1,r1,0 l.lwz reg,0x0(r1) p LW reg,0(r1) GET OR TREE (reg) read the value of the OR Tree addi r1,r0,0x9005 l.addi r1,r0,0x9005 LW reg,0(r1) l.lwz reg,0x0(r1) Mbyte SSRAM. Briefly, the design is composed of a camera TRDB D5M equipped with 5 Mega-Pixel CMOS sensor and working with 15 frames per second, an SDRAM to store pixel data coming from the camera, the parallel architecture to process the pixels, an SRAM to store computed data and the 20 800×RGB×480 LCD TRDB LTM displayer working at 33.2 MHz (Figure 6).

Capturing and displaying video

In this design, the video processing system carries out the capture, buffering and video display. The SIMD parallel architecture is configured and explored with different number of PEs and with the global NoC to manage the needed I/O communications since I/O transfers present a performance challenge for this kind of high performance computing systems. To validate the configurability of the processor IP, two processors are chosen: miniMIPS [START_REF] Opencores | minimips overview[END_REF] and OpenRisc [START_REF] Opencores | OpenRISC processor[END_REF]. The two processor design methodologies (reduction and replication) are also tested. The NoC is implemented based on a crossbar interconnection network in order to assure a rapid data-parallel transfer to all PEs. It is chosen because it is a non blocking fully connected network in comparison with a multi-stage network, which is blocking and to the bus that has a limited bandwidth. Based on the parametric VHDL implementation, it was easy to generate and synthesize these different parallel VHDL configurations. The FPGA based synthesis results are presented in terms of consumed LEs as well as total memory bits in the case of the two methodologies. The consumed power and the maximum obtained frequency are also illustrated (Table 3 and Table 4).

Tables 3 and4 show, as expected, that the processor reduction methodology allows implementing a large number of PEs on a single chip. With the processor replication methodology we can not put 16 miniMIPS PEs on the CycloneII FPGA, whereas with the reduction methodology up to 32 miniMIPS reduced PEs can be integrated on the same FPGA device.

Comparing between the two used processors, we notice that the miniMIPS occupies smaller area than the OpenRisc and slightly has a reduced power consumption. Another important issue is the consumed memory. From Table 4, we notice that when doubling the number of integrated PEs, the total memory blocks will be 1.6% higher. The major SoC constraint is the limited amount of on-chip memory. Results also show that the power consumption slightly increases when doubling the number of integrated PEs because the consumed power depends on the number of components in the architecture as well as the interconnection node capacitance, circuit voltage and switching frequency [START_REF] Curd | Power Consumption in 65 nm FPGAs[END_REF].

In this design, each PE reads one pixel at a time from the SDRAM through the global network, then sends it in the same manner to the SRAM. According to previous results, the reduced miniMIPS is chosen to build hardware SIMD configurations for the following experiments.

Figure 8 shows the processing time results obtained on reduced PEs (8 and 16) with different image sizes. Figure 8 highlights that the processing time scales with the image size. For example, the time needed to process one image of size 640×480 is 2.6 higher than the time spent to process an image of size 320×240. This demonstrates the performance and the scalability of the system since it depends on the number of processors and the size of the data being worked on.

RGB to YIQ color conversion

RGB to YIQ color conversion is used in many encoders where the RGB inputs from the camera are converted to a luminance (Y) and two chrominance information (I,Q). This algorithm explores multiply/accumulate capability, which is especially dedicated to SIMD and VLIW architectures. This conversion's type can be accomplished using the following equations: An image filter convolution is tested in the following subsection.

Y I Q = 0.

Image filter convolution

An image's convolution with a mask is tested. As a case study, a sharpen filter (3x3 kernel) is selected:

0 -1 0 -1 5 -1 0 -1 0 .
The convolution of the mask can be done in parallel following the SIMD functioning so that multiple pixels can be processed at a time. The image convolution algorithm is quite similar to the color conversion algorithm since they perform image processing computations. Therefore, according to previous results we consider that we need more than 16 PEs to assure real-time processing.

The two tested parallel SIMD configurations are composed of an ACU, 32 PEs (each PE is connected to its local data memory) and a crossbar based NoC. The second configuration also integrates a 2D torus neighboring interconnection network compared to the first one. Table 5 shows the needed FPGA resources for the two designs obtained via Quartus tool. It is clear that the second parallel configuration couldn't be fitted on the Cyclone II FPGA.

To compute N convolved pixels at a time, each PE reads five pixels at a time (corresponding to the 5 non-zero filter coefficients), does the same computations and stores the resultant pixel. Each PE needs to have at least a local memory with 150 words (600 bytes). The LCD displayer begins reading the resultant data from the SRAM after the initialisation phase and the processing of the 32 first pixels in order to assure real-time processing.

The processing of one pixel doesn't exceed 30.12 ns since it is the elapsed time to read one pixel by the LCD. In this algorithm, to process one pixel each PE does 5 memory loads, 1 multiplication, 4 negations, 4 additions and 1 memory store. Through prototyping, a pixel-processing time equal to 28.75 ns < 30.12 ns is reached.

This implemented algorithm doesn't consider the neighboring pixels between PEs. So, a second algorithm version is implemented to take into account these pixels. In this case, the design can not be fit in the FPGA (table 5). The simulation results have shown a pixel-processing time equal to 47.5 ns > 30.12. Thus, an SIMD SoC integrating more than 32 PEs and an FPGA with more resources than the Cyclone EP2C70 are needed to respond to real-time constraints.

All these experiments show the scalability of the proposed parallel system, demonstrates that its peak performance scales linearly with the number of PEs and proves the efficiency of the proposed design methodology to explore different parallel configurations in order to choose the best one.

Performance

For performance comparisons, we adopted the image filter convolution benchmark. The tested parallel SoC is based on replicated NiosII processors, running at 200 MHz. Our system is firstly compared to a special purpose FPGA-based 2D systolic architecture dedicated to implement window-based image processing algorithms [START_REF] Torres-Huitzil | Real-time image processing with a compact FPGA-based systolic architecture[END_REF]. Table 6 summarizes the obtained results.

Results show the efficiency of the proposed system since it achieves better processing time results compared to a 2D systolic architecture even with smaller number of PEs (32 compared to 49).

The SIMD system is then compared to hardware acceleration based on the C-to-HW (C2H) acceleration compiler of the Nios processor [START_REF] Altera | Nios II C2H Compiler User Guide[END_REF]. We ran the application on the NiosII/f processor at 200 MHz to establish the baseline performance. The two executions are compared in Table 7. Compared to a dedicated HW acceleration, the proposed parallel and programmable SIMD system achieves a better speed up. This demonstrates the efficiency of the SIMD SoC with the advantage of flexibility compared to a HW based implementation. We also compare our system to the NIOS implementation based on customized SIMD instruction extension [START_REF] Etiemble | Introducing image processing and SIMD computations with FPGA soft-cores and customized instructions[END_REF]. Table 8 shows experimental results.

The results are presented with the Cycle per Pixel metric (CPP), which is the total number of clock cycles divided by the number of pixels, running the convolution on 256x256 gray image.

The results clearly demonstrate that our SIMD system achieves a speed up equal to 4.8 compared to the sequential execution performed by the ACU. This speed up is more than 2 times greater than when using customized SIMD instructions with the NiosII processor.

Table 9 compares between our system composed of miniMIPS reduced

PEs and the FPGA-based SIMD processor described in [START_REF] Tanabe | A study of an FPGA based flexible SIMD processor[END_REF] composed of one control unit and a number of computation units (CU) in terms of power consumption and frequency. As illustrated, our system achieves good results since it presents a higher frequency and a lower power consumption compared to the SIMD processor. The proposal processor reduction methodology is easier to design the SIMD architecture with a parametric number of PEs compared to the used dedicated implementation [START_REF] Tanabe | A study of an FPGA based flexible SIMD processor[END_REF] that makes adding more computational units to the SIMD processor a heavy task.

The preceding performance comparison confirms the promising advantages of the proposed FPGA-based SIMD architecture. Its performance has been evaluated for image convolution benchmark with good results that validate the proposed high performance architectural model. In this paper, we have shown that an FPGA-based SIMD design built with soft-core processors is powerful to compute image processing. One of the main advantages of the proposed architecture is its configurability. The architecture is scalable and flexible enough to support several applications.

Considering multimedia applications that comprise not only SIMD operations but also operations suitable for exploiting parallelism too, we can benefit from the proposed system. In fact, in vision algorithms for example, some tasks are performed on the objects found to analyse their quality or properties in order to make decisions on the image contents. It appears that SIMD type of architectures is not very efficient for this treatment. A DSP is often more appropriate. So, we can replace the ACU by a DSP. Our proposed system can also assure task parallelism by parallelizing the program between the ACU and one PE. The proposed many-core SIMD system can also work as an accelerator for another system (which can be a unique processor or an MPSoC) or as a coprocessor of the main CPU reducing its workload and accelerating the data-parallel expensive computation. In this case, it executes the data-parallel parts of the program, and the ACU can process any task of the program (while the PEs are working). This manner may introduce some control and communication overhead, which can be ameliorated via tailored components. As example, to reduce data access latency, we can integrate a

Conclusion

This work presents an SoC platform to facilitate generating SIMD massively parallel on-chip architectures. It proposes an IP-based implementation to build these architectures leading to the design of large and complex architectures with lower costs and higher performance. Different SIMD configurations can be derived from the generic SIMD SoC model. The proposed massively parallel system is characterized by its flexibility allowing matching the design with the application as well as improving its performances and satisfying its requirements. This flexibility allows a designer to choose the area/performance of a parallel system without laborious hardware design, and can better meet application requirements. The proposed RTL implementation makes any hardware modification easy in order to adapt the system to a wide range of data dependent algorithms. All these features strongly contribute to the increase of the designer's productivity. Through prototyping results, the user is able to choose the appropriate SIMD SoC configuration satisfying his needs and meeting performance requirements.

This work opens an interesting topic for future research and development on parallel applications. Future work will be to provide a high-level design framework. An automatic exploration level can help the designer to generate the most appropriate and efficient massively parallel configuration for a given application. Another future work is to study the dynamic runtime use of reduction and replication proposed methodologies.

[1] A.

Figure 1

 1 Figure 1 highlights the proposed parallel platform with a focus on its IPbased design. The used IPs are mainly memories (RAM, ROM), processors and routers. The SIMD model is flexible and can be customized to target diverse applications. The whole system is synchronously controlled by the ACU, which is responsible for fetching and decoding instructions, computing addresses and scalar data values, issuing control signals to the PE array, and monitoring the status of the PE array. A data-parallel program is composed of sequential

Figure 1 :Figure 2 :

 12 Figure 1: IP-based SIMD SoC platform

Figure 3 :

 3 Figure 3: Processor reduction methodology applied on miniMIPS processor

Figure 4 :

 4 Figure 4: Processor replication methodology applied on miniMIPS processor

  hardware components (PE-PE, ACU-PE, PE-I/O device). It provides global communication between all PEs and forms the basis for the parallel I/O system. It mainly consists of a communication mode manager and an interconnection network IP. The mode manager is responsible of establishing the needed communication mode.The internal interconnection network transfers data from sources to destinations and can have different types (shared bus, full crossbar, Delta Multistage Interconnection Network). The NoC IP is parametric (in terms of size and internal network) to offer powerful integration in various systems[START_REF] Baklouti | Scalable mpNoC for massively parallel systems Design and implementation on FPGA[END_REF].Neighborhood regular networkThis network performs neighbor communications between PEs and can have different topologies (linear array, ring, mesh, Torus, and Xnet (a twodimensional mesh with extra diagonal links)). It may connect the data in a specific direction depending on the destination as well as the distance, which defines the number of paths between every couple of PE sender and PE receiver. It is composed of a controller and several routers (equal to the number of PEs in the system). The architecture of the router IP is depicted in Fig 5.The data communicated through the router are not stored in FIFO registers since the neighborhood communications are synchronous and all performed in a single direction at a given time. An SIMD interconnection function is a bijective function, thus the data transfer occurs without conflicts. In every communication, each router only activates the required link to the needed direction. The other links will be disabled, reducing the power consumption.

Figure 5 :

 5 Figure 5: Regular router IP architecture

Figure 6 :

 6 Figure 6: FPGA-based Video Processing Chain

This

  NoC performs the PE-I/O communication mode: PE-SDRAM(I/O) and PE-SRAM(I/O). The displayer is directly connected to the SRAM to 22

Figure 7 :

 7 Figure 7: Experimental results to process one image of size 800×480 pixels

Figure 8 :

 8 Figure 8: Experimental results to process different sized images on reduced miniMIPS PEs

  The execution time results as well as the speed up values are presented in Figure9depending on the number of PEs. The number 1 in the xaxis means that the SIMD configuration only contains the ACU to measure the sequential execution time. As shown in the previous experiment, the time reduces when increasing the number of PEs. Compared to the first application, only computing operations are added while maintaining the same communication transfers through the NoC. Thus, we clearly notice that the speed up is almost the same. To reach real-time, a minimum of 8 PEs is needed in this application. Each pixel processing should not exceed 30.12 ns.

Figure 9 :

 9 Figure 9: Color conversion experimental results

Table 1 :

 1 Consumed logic of the miniMIPS processor and the reduced miniMIPS processor

	Processor	ALUTs	Registers
	miniMIPS	3339	1936
	reduced miniMIPS	1265	1318
	4.1.1. Processor IP
	The proposed design is based on Open hardware concept. So, open-
	source processor IPs are mainly chosen to implement SIMD configurations.

Table 2 :

 2 New instructions

	Macro ASM	Description		Coding
			miniMIPS	OpenRisc
	SET MODE NOC (reg)	NoC mode instruction: select comm.	SW reg,0x9003(r0) l.addi r1,r0,0x9003
		mode (reg value).		l.sw 0x0(r1),reg
	Modes 0,1 and 4:	NoC SEND instruction: transfer	p addi r1,r0,dest	l.addi r1,r0,dest
	P NOC SEND	data reg to the needed	p addi r1,r1,adr	l.addi r1,r1,adr
	(reg,dest,adr)	destination specified by dest.	p SW reg,0(r1)	l.sw 0x0(r1),reg
	Modes 2 and 3:		p addi r1,r0,adr	l.addi r1,r0,adr
	P NOC SEND (reg,adr)		p SW reg,0(r1)	l.sw 0x0(r1),reg
	Modes 0,1 and 4:	NoC RECEIVE instruction: read	p addi r1,r0,src	l.addi r1,r0,src
	P NOC REC (reg,src,adr)	data from appropriate sender.	p addi r1,r1,adr	l.addi r1,r1,adr
			p LW reg,0(r1)	l.lwz reg,0x0(r1)
	Modes 2 and 3:		p addi r1,r0,adr	l.addi r1,r0,adr
	P NOC REC (reg,adr)		p LW reg,0(r1)	l.lwz reg,0x0(r1)
	P REG SEND	Neighboring SEND instruction: transfer	p addi r1,r0,dir	l.addi r1,r0,dir
	(reg,dir,dis,adr)	data (of reg) from source to	p addi r1,r1,dis	l.addi r1,r1,dis
		destination.	p addi r1,r1,adr	l.addi r1,r1,adr
			p SW reg,0(r1)	l.sw 0x0(r1),reg
	P REG REC (reg,dir,dis,adr)	Neighboring RECEIVE instruction:	p addi r1,r0,dir	l.addi r1,r0,dir
		read data from appropriate	p addi r1,r1,dis	l.addi r1,r1,dis
		sender.

Table 3 :

 3 Synthesis results -Processor reduction methodology

	Nb	Processor	Logic Utilization		Total memory		Power	Fmax
	PEs	IP	Comb.	registers	% ACU	PE	%	cons.
			functions		LE (bytes)	(bytes)	Mem	(mW)	(MHz)
	4	miniMIPS	14176	4762	23	4096	4096	24	852	89.21
	8	miniMIPS	21425	6006	34	4096	4096	33	940	87
	16	miniMIPS	35673	8366	55	4096	4096	52	1296	85.82
	32	miniMIPS	61130	13006	93	4096	2048	66	1969	83.04
	4	OpenRisc	19674	7678	30	4096	4096	15	911	188.61
	8	OpenRisc	34667	13534	53	4096	4096	22	1205	186.06
	16	OpenRisc	64270	25160	98	4096	4096	36	1999	185.22

Table 4 :

 4 Synthesis results -Processor replication methodology

	Nb	Processor	Logic Utilization		Total memory		Power	Fmax
	PEs	IP	Comb.	registers	% ACU	PE	%	cons.
			functions		LE (bytes)	(bytes)	Mem	(mW)	(MHz)
	4	miniMIPS	21971	9243	37	4096	1024	11	1207	93.12
	8	miniMIPS	42059	16863	71	4096	1024	18	1799	90.91
	4	OpenRisc	18491	10229	41	4096	1024	13	1298	192.52
	8	OpenRisc	40185	23966	91	4096	1024	22	1921	189.97

Table 5 :

 5 Synthesis results for the convolution SIMD designs

	Neighb.	Logic Utilization		Total memory
	network	Combin.	reg	% ACU	PE	%
		functions			(bytes)	(bytes)
	none	61130	13006	93	2048	600	22
	configured 2D Torus	81796	22478 119	2048	600	25

Table 6 :

 6 Comparison between SIMD SoC and 2D systolic architectures

	System	Image size	Number of rows	Processing Time
			processed in parallel	(Ms)
	SIMD SoC	512x512	15	3
	(32 PEs)	(gray-level)
	2D systolic	512x512	15	5
	(49 processors) (gray-level)

Table 7 :

 7 Comparison between SIMD SoC and C2H HW acceleration

	System	Speed up Frequency (MHz)
	SIMD SoC	17x	184
	C2H	13.3x	95

Table 8 :

 8 Comparison between SIMD SoC and SIMD instruction extension

	SIMD instructions	CPP	SIMD system	CPP
	extension
	C version	34.7	ACU version	62.66
	C version with SIMD instructions	17.4	SIMD version 12.85
	Speed up	2	Speed up	4.8

Table 9 :

 9 Comparison between SIMD SoC and[START_REF] Tanabe | A study of an FPGA based flexible SIMD processor[END_REF]

	Proposed SIMD SoC Fmax Power (mW)	[24]	Fmax Power (mW)
	4 PEs	89.21	852	4 CUs	56.2	1149
	8 PEs	87	940	8 CUs	55.3	1233
	16 PEs	85.82	1296	16 CUs	52.7	1385

  Gentile, and D. Scott Wills, Portable Video Supercomputing, IEEE Trans. Comput. Volume 53, Issue 8, Aug. 2004.

DMA controller between the memories of the PEs and the main memory. Some work have been proposed to exploit SIMD architectures for computer vision algorithms like the SIFT (Scale-invariant feature transform) algorithm. In [START_REF] Su | Novel Parallel Approach for SIFT Algorithm Implementation[END_REF], authors propose a parallel approach for SIFT algorithm implementation using a block filtering technique in a Gaussian convolution process on the SIMD Pixel Processor. They demonstrate that their system can perform real-time processing with good performance. In [START_REF] Nieto | SIMD/MIMD dynamicallyreconfigurable architecture for high-performance embedded vision systems[END_REF], authors propose a hybrid SIMD/MIMD architecture for image processing.

The architecture is composed of a signal-processing CPU and an image coprocessor. This latter supports two working modes SIMD and MIMD. In [START_REF] Maresca | Parallel Computer Vision on Polymorphic Torus Architecture[END_REF] authors demonstrate that the interconnection network used in the SIMD architecture influences the performance of computer vision applications. Considering a set of basic vision tasks namely convolution, histogramming, hough transform, extreme point identification, etc. authors show the performance and the advantages of the polymorphic torus communication when executing these tasks. In [START_REF] Wu | Speeding up SIFT algorithm by multi-core processor supporting SIMD instruction sets[END_REF], a parallelization of SIFT using multi-core architecture with per-core SIMD support is shown. Authors present good results compared to the State-of-the-Art parallel SIFT algorithms. So, we conclude from the aforementioned implementations that we can use a pure SIMD system while optimizing/modifying the algorithm to take benefit from the HW architecture, or tailor the architecture to the requirements of the algorithm.