P Chinda

S Chanchaona

P Brault

W Wechsatol

Mathematical Modeling of a Solid Oxide Fuel Cell with Nearly Spherical-Shaped Electrode Particles

Keywords: micro-scale model SOFC, microstructure; optimal size, overpotential, spherical-shaped particle

A micro-scale model of a Solid Oxide Fuel Cell (SOFC) involving electrochemical reactions and mass transport behaviors throughout electrodes packed with nearly spherical-shaped particles was mathematically developed. The model was validated by comparing the predicted results with already available experimental results. By minimizing the cell overpotential, the model provides the optimal size of electrode particles corresponding to the specified porosities. The optimal volumetric ratio between ionic and electronic conducting particles is also proposed. The study results substantiate the fact that SOFC overpotential could be effectively decreased by increasing the operating temperature as well as operating pressure. The present study reveals the working mechanisms of SOFC at the micro-scale level, while demonstrating the use of microstructure relations to enhance the SOFC performance.

Introduction

The solid oxide fuel cell (SOFC) is a highly efficient energy conversion system that transforms chemical energy to electrical energy and heat directly from gasified fuels by electrochemical reactions of fuels and oxidants. It can be powered with hydrocarbons such as CH 4 , C 2 H 6 and C 3 H 8 , and also with CO besides H 2 as a fuel sources. The SOFC can be applied to versatile power generation systems, such as stationary power supplies and auxiliary power sources for automotive and residential applications [START_REF] Singhal | High temperature solid oxide fuel cells: fundamentals design and applications[END_REF][START_REF] Larminie | Fuel cell systems explained[END_REF]. Researchers have been mainly focused on the SOFC fabrication with various forms, and developing new materials for commercial applications. Ramakrishna et al. [START_REF] Ramakrishna | Innovative design to improve the power density of a solid oxide fuel cell[END_REF] performed numerical analyses on a CFD-ACE package and proposed an innovative thin-walled geometry to improve the power density of SOFC. More recently, a novel 3-D SOFC was developed by Koh et al. [START_REF] Koh | Design and fabrication of three dimensional solid oxide fuel cells[END_REF] to increase the volumetric power density by increasing the volumetric surface areas of the electrodes to express its method of fabrication. Moreover, Yang et al. [START_REF] Yang | Computational analysis of thermo-fluid and electrochemical characteristics of MOLB-type SOFC stacks[END_REF] and Hwang et al [START_REF] Hwang | Computational analysis of species transport and electrochemical characteristics of a MOLB-type SOFC[END_REF][START_REF] Hwang | Detailed characteristic comparison between planar and MOLB-type SOFCs[END_REF] developed a numerical model to investigate the characteristics of a mono-block-layer built (MOLB) SOFC and they found that the MOLB-SOFC had a higher fuel/oxidant utilization than the planar type SOFC. The high utilization of the fuel/oxidant reflects the high current through the SOFC. However, the current density distributions are less uniform in MOLB-type SOFC, which is a disadvantage in this type of SOFC. The major limitation of using SOFCs is their high operating temperature. As a solution to this, and to achieve stability and economy, intermediate temperature SOFCs (IT-SOFCs) and low temperature SOFCs are subject to a growing interest. The techniques of thin electrolyte and electrode coating involving new materials have been chosen for intermediate and low temperature SOFCs. Furthermore, Yttria Stabilized Zirconia or YSZ is normally used as the electrolyte for cell operation at high temperatures, however, by reducing the electrolyte thickness it is possible to use YSZ for SOFCs at intermediate and low temperatures. New materials used in SOFCs for intermediate temperature SOFCs such as ceria based electrolyte have been proposed [START_REF] Huijsmans | Intermediate temperature SOFCa promise for the 21 st century[END_REF][START_REF] Choy | The development of intermediatetemperature solid oxide fuel cells for the next millennium[END_REF][START_REF] Ralph | Materials for lower temperature solid oxide fuel cells[END_REF][START_REF] Steele | Materials for IT -SOFC stacks 35 years R&D: the inevitability of gradualness ?[END_REF], but for overcoming the major limitations of SOFCs, a new technique of SOFC design has been investigated by many researchers: thin film deposition. This technique has been studied for a few decades with Proton Exchange Membrane Fuel Cells (PEMFC) [START_REF] Rabat | Plasma sputtering deposition of PEMFC porous carbon platinum electrodes[END_REF] and now shows immense possibility to be applied to SOFC applications [START_REF] Beckel | Thin films for micro solid oxide fuel cells[END_REF][START_REF] Morse | A novel thin film solid oxide fuel cell for micro scale energy conversion[END_REF][START_REF] Pederson | Application of vacuum deposition methods to solid oxide fuel cells[END_REF]. One of the main advantages of such a thin film fuel cell is that its ohmic resistance is low enough to maintain a high performance cell at low operation temperatures.

In a SOFC system, computer simulation based on theoretical modeling is known to be a very efficient method for predicting materials properties, reaction kinetics and unit cell performance. Different levels of model, depending on the objective, have been developed, each with different degrees of sophistication and detail. Theoretical modeling of SOFC can be divided into two levels: macroand micro-scale modeling. Detailed understanding of fluid and mass transport shows whether heat transfer can be obtained from macroscale modeling. The effects of porosity change on reactant concentration and ohmic overpotentials could also be considered in the macro-scale modeling, but the activation overpotential must be treated as a property independent of the porosity, while the electrochemical reactions must occur only at the interface between electrode and electrolyte layers. The macro-scale modeling is not a good predictor of the effect of electrode structure on the chemical reaction within SOFCs. Influence of the electrode structure on the reaction of kinetics at the three phase boundary (TPB) has recently been included in micro-scale modeling of SOFCs [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF][START_REF] Deseure | Theoretical optimization of a SOFC composite cathode[END_REF][START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF][START_REF] Ni | Micro-scale modeling of solid oxide fuel cells with micro-structurally graded electrodes[END_REF]. In literature for micro-scale modeling, the properties of porous electrodes and their effects on SOFC performance have been investigated from different aspects.

Costamagna et al. [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF] developed a theoretical microscopic model for the evaluation of the performance of an electrode formed by a mixture of electronic and ionic conductors with nearly spherical particle shape. The results of the mode1 show that the effects of morphology, i.e. the volumetric composition of the electrode and the dimensions of the particles, strongly influence the electrode resistance. The model has been applied to erbia-doped bismuth oxide/Pt or (Bi 2 O 3 ) 0.7 (Er 2 O 3 ) 0.3 /Pt or EDB/Pt and Yttria-tabilized Zirconia/Strontium-doped Lanthanum Manganite (YSZ/LSM) cathodes and YSZ/Ni anodes. Chen et al. [START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF] developed a mechanistic model for oxygen reduction at the YSZ/LSM interface; a complete micro-model for YSZ/LSM composite cathode considering all forms of overpotential was developed which established the interrelationship among the transport phenomena, electrochemical processes and the microstructure of the composite cathode. The exchange current densities of the rate-limiting steps used in the simulation and the model can be used as a tool to guide the design of high performance cathodes. Deseure et al. [START_REF] Deseure | Theoretical optimization of a SOFC composite cathode[END_REF] developed a microscopic model for a composite cathode similar to Chen et al. [START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF]. Simulation was conducted to predict the optimal design parameters, e.g. cathode thickness, particle size, particle size ratio and YSZ volume fraction, for a LSM/YSZ composite cathode. Hussain et al. [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] considered a 1-D anode-supported planar SOFC with a thin layer reaction zone in the vicinity of electrolytes. Their numerical results showed that the increase of either porosity or tortuosity of electrodes leads to worse cell performance, while better cell performance is obtained when the volume fraction of electronic conducting particles is approximately equal to that of ionic conducting particles in the reaction zone. Recently they extended their investigation to the cell design with two distinct layers of anode electrode, i.e., reaction layer and backing layer. Ni et al. [START_REF] Ni | Micro-scale modeling of solid oxide fuel cells with micro-structurally graded electrodes[END_REF] developed a mathematical model for modeling the performance of SOFCs with functionally graded electrodes at the micro-scale level. The model considered all forms of overpotentials and was able to capture the coupled electrochemical reactions and mass transfer involved in the SOFC operation. Additional modeling analyses were conducted to gain better understanding of the SOFC working mechanisms at the micro-scale level and to quantify the performance of microstructurally graded SOFCs. It was found that micro-structural grading could significantly enhance the gas transport but had negligible effects on the ohmic and activation overpotentials, especially for thick electrodes. The optimal size of particles was investigated for maximizing SOFC cell performance. The volumetric fraction between ionic and electronic conducting particles was also considered.

Unlike other previous researches, the effect of particle's geometrical shapes is included in this model. For example, porosity and tortuosity are parameters that depend on the shape of electrode particles. Normally, porous electrodes are random stacks of nearly spherical-shaped particles. Thus, this work is aimed to maximize the SOFC performance by optimizing the microstructure of the SOFC electrodes. The optimal size of particles was investigated for maximizing SOFC cell performance. The volumetric fraction between ionic and electronic conducting particles was also considered.

Experimental

Mathematical modeling

Anode side modeling

Porous composite electrodes are commonly used in fuel cells due to their extended zone for electrochemical reactions, thus improving the current output of an electrode. The performance of a porous electrode depends on, but is not limited to: (1) the activity of electro catalysts, (2) the specific surface area (m 2 m -3 ) available for electrochemical reaction, (3) the mass transport phenomena, and (4) the ionic and electronic conductivity of the electrode. Because of the complexity of mass transport, the ionic and electronic transport and electrochemical reaction, so a complicated simulation technique was employed to predict cell performance under various operating conditions, and for parametric study. The present model aims to optimize the SOFC performance by optimizing the microstructure relations of SOFC electrode by micro-scale modeling. The relation of electrode porosity and electrode particle radius size has been investigated for maximizing SOFC cell performance. Moreover, this model is considered for optimization of the volumetric fraction of ionic and electronic conducting particles of the electrodes. For the sake of simplicity, we assume that the particle radius of both ionic and electronic conducting particles is equal.

The configuration of a SOFC anode can be simply modeled as shown in Figure 1a. For modeling, the coupled electrochemical reaction and mass transfer in a composite anode are illustrated by Figure 1b, and the reaction processes can be summarized as: (1) transport of reactant (H 2 ) to the reaction sites through the pores of the anode and transport of O 2-from the electrolyte to the reaction sites through the ionic conducting particles; (2) electrochemical reaction of H 2 and O 2-to form H 2 O and electrons at the active sites;

(3) transport of electrons from the active sites to the current collector through the electronic conducting particles and transport of the H 2 O product to the anode surface via the pores of the anode. The configuration of the composite cathode of a SOFC also consists of both electronic and ionic conducting particles. Therefore, the mathematical model of the cathode side can be formed in the similar manner. Charge balance in the electronic and ionic conducting particles can be written as: e,a v n,a i,a J = -A J = -J

where J e,a and J i,a is the current density (Am -2 ) in the electronic and ionic conducting particles, respectively, A v is the electrochemically active surface area per unit volume of the porous anode (m 2 m -3 ), which depends on the anode properties, and J n,a is the transfer current per unit area of reactive surface (Am -2 ). A v can be calculated with the method based on the particle coordinate number in binary random packing of spheres proposed by Costamagna et al. [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF], given as:

22 el io v el t el io el io ZZ A = sin r n n n p p Z (2)
where is the contact angle between electron and ion conducting particles, r el is the radius of the electron conducting particles, n t is the total number of particles per unit volume, n el and n io are the number fraction of electron and ion conducting particles, Z el and Z io are the coordination number of electron and ion conducting particles. Z is the total average particle coordination number in the electrode which is equal to 6 [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF] and p el and p io are the probabilities for an electron and an ion conducting particle to belong to connecting ends of the composite. All previous mentioned parameters required to calculate A v [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF]

where is a volume fraction of the electron conducing particle in an electrode, Z el -el represents the average coordination number between electronic particles and Z io-io represents the average coordination number between ionic particles [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF]

1- n = (4/3) r [n (1 n )(r / r ) ] ( 7
)
where ε is the electrode porosity and from Eq. ( 7), can be rewritten as:

3 3 el t el el io el (4/3) r n [n (1 n )(r / r ) ] + = 1 (8)
For the sake of simplicity, assume that both ionic and electronic conducting particles have the same radius (r io = r el ), thus Eq. ( 8) can be rewritten as,

3 t el 4 n r + = 1 3 (9)
Eq. ( 9) illustrates the constraint among , n t and the particle radius (r el , r io ). If the total number of particles n t is fixed, will depend solely on the particle radius. An aim of this work is to find the optimal radius of particles leading to the maximum active surface area A v . According to the Butler and Volmer equation [START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF], the transfer current density J n,a can be determined from J is the reference exchange current density of the electrode -i.e. for an anode, ref 0,a J is 400 A/m 2 and for a cathode, ref 0,c J is 1320 A/m 2 [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF]. P H2O and HO 2 0 P are respectively the partial pressures of H 2 O within the porous anode and at the anode surface.

Similarly, P H2 and H 2 0 P are the partial pressures of H 2 within the porous anode and at the anode surface. The parameter α is the charge transfer coefficient which is 0.5 [START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF], while z is the number of charges involved in a reaction. The overpotential at the anode η a [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF] is

0 0 a e,a i,a e,a i,a = ( - ) -( - ) (11)
where i,a and e,a denote the ionic and electronic potentials. The superscript 0 represents the equilibrium state obtained from Ohm's law, eff e,a e,a e,a = J and i,a i,a eff i,a = J [START_REF] Rabat | Plasma sputtering deposition of PEMFC porous carbon platinum electrodes[END_REF] where ρ eff [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] is the effective resistivity (Ωm) and J is the current density. The subscripts i and e represent ionic and electronic conductors, respectively. ρ eff can be determined by Eq. ( 13) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF].

where τ is electrode tortuosity, σ e is the conductivity of electronic conductor i.e. for an anode, e,a is 6 2 10 Sm -1 [16,21] and for a cathode, e,c is 4 1 10 Sm -1 [22]. σ i is the ionic conductor i.e. for an anode i,a and for a cathode, i,c is equal to 4

At the anode surface, J i,a = 0 and J e,a = J total . At the electrolyte interface, J i,a = J total and J e,a = 0. Therefore, two boundary conditions for the second order differential Eq. ( 16) can be written as, eff a e,a total x 0 d x = 0 , = -J dx and eff a an total i,a x t an d x = t , = J dx [START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF] The electrochemical reactions due to mass transport occur in SOFC electrodes. Mass transport through a porous medium can be determined using concepts described by Krishna et al. [START_REF] Krishna | The Maxwell -Stefan approach to mass transfer[END_REF]. For a single component, the mass transport equation can be written as,

-3 -1 i i i (y P) = -N + R (mol m s ) RT t (18)
where ε and N i represent the porosity and the rate of mass transport (or denotes as the molar flux of each species, mol m -2 s -1 ), respectively, into porous media. R i is the rate of reaction inside the porous medium. The term on the left-hand side is valid when an unsteady state is approached. The first and second terms on the right-hand side represent the diffusion rate and the rate of reaction inside the porous medium. It was assumed earlier that the diffusion process is at a steady state and that the electrochemical reactions take place at the three phase boundary (TPB) inside the electrode. Therefore, within the anode, the first and second term on the righthand side is significant, Eq. ( 18), and becomes

i i N = R (19)
The rate of mass transport, N i , generally depends on the operating conditions (reactant concentration, temperature and pressure) and the microstructure of the material (porosity, tortuosity and pore size). Three models were used to develop expressions for N i : Fick's, the Dusty-Gas and the Stefan-Maxwell models. Only, Fick's and Dusty-Gas are used to describe the transport of components within porous media but the Stefan-Maxwell Model (SMM) is a well-known mass transport model applied to nonporous media [START_REF] Suwanwarangkul | Performance comparison of Fick's, Dusty -Gas and Stefan -Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF].

Normally, to describe the transport of components within porous media the simplest method is using Fick's model that takes into account diffusion and convection transport. The convection transport term is described by the Darcy equation that relates with the permeability, viscosity of fluid, total pressure and pressure gradient [START_REF] Suwanwarangkul | Performance comparison of Fick's, Dusty -Gas and Stefan -Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF]. However, it is assumed that the total pressure change within the pore is insignificant ( dP dx 0 ). Thus, only the diffusion transport is necessary to determine the rate of mass transport by Fick's model. Another model to describe the transport of components within porous media is Dusty-Gas Model (DGM) that also includes the permeation flux and is applied to take into account the effect of total pressure gradient on mass transport. However, this study was done on the basic assumption that the total pressure is uniform over the entire depth of the porous electrode. Thus, the definite total pressure gradient term can be ignored. Consequently, only the diffusion flux is addressed for the DGM model, which means the primary mode of species transport in the porous electrode is by diffusion.

The mass transfer process inside the porous electrode is by means of diffusion and was also selected for the reason described above and by recommendations in the literature [START_REF] Suwanwarangkul | Performance comparison of Fick's, Dusty -Gas and Stefan -Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF][START_REF] Yakabe | Evaluation and modeling of performance of anode-supported solid oxide fuel cell[END_REF][START_REF] Lehnert | Modeling of gas transport phenomena in SOFC anodes[END_REF]. Suwanwarangkul et al. [START_REF] Suwanwarangkul | Performance comparison of Fick's, Dusty -Gas and Stefan -Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF] concluded that the Dusty-Gas Model (DGM) is the most appropriate model to simulate gas transport phenomena inside a SOFC electrode. The mass flux of species i can be determined by the Dusty-Gas Model (DGM) which is considered to be both Knudsen diffusion and Molecular diffusion. The DGM is the most suitable model for the H 2 -H 2 O system [START_REF] Suwanwarangkul | Performance comparison of Fick's, Dusty -Gas and Stefan -Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF]. Suwanwarangkul et al. also demonstrated that only the DGM is recommended for a multi component system. This is because it takes into account Knudsen diffusion effect as well as Graham's law of diffusion to calculate the flux ratios. Thus, the DGM was used to develop the expressions for N i in this work.

The electrochemical reactions due to mass transport occuring in SOFC electrodes is represented by the term R i. ; the rate of electrochemical reaction per unit volume of the porous electrode (mol m -3 s -1 ). The electrochemical reaction rate of H 2 O and H 2 can be expressed by Faraday's Law as:

2 H v n,a R = (-A J )/(2F) (20) 2 H O v n,a R = (A J )/(2F) (21)
The DGM includes the Stefan-Maxwell formulation and takes into account Knudsen diffusion and Molecular diffusion [START_REF] Suwanwarangkul | Performance comparison of Fick's, Dusty -Gas and Stefan -Maxwell models to predict the concentration overpotential of a SOFC anode[END_REF]. It is assumed from this model that pore walls consist of giant molecules ('dust') uniformly distributed in space. These dust molecules are considered to be a dummy, or pseudo species in the mixture. The general form of the DGM is shown in Eq. ( 22):

n j i -i j 0 i i i eff eff eff j=1,j i i,k ij i,k y N y N BP N dy 1 dP + = - P y 1 RT dx dx D D D n j i -i j 0 i i i eff eff eff j=1,j i i,k ij i,k y N y N BP N dy 1 dP + = - P y 1 RT dx dx D D D ( 22
)
where x is a vertical direction as shown in Figure 1a and y i , y j is the mole fraction of species i and j.

From the DGM that considers both Knudsen diffusion and Molecular diffusion, the diffusion process within a pore can be typically divided into two diffusion mechanisms that are Molecular and Knudsen diffusion. Molecular diffusion is dominant for large pore sizes and high system pressures while Knudsen diffusion becomes significant when the mean-free path of the molecular species is much larger than the pore size. Therefore the effective diffusivity of species i, ( eff i D ) can be written by combining the effective molecular diffusion ( eff im D ) and effective Knudsen diffusion ( eff i,k D ) coefficients as follows Bosanquet formula [27]:

1 eff i eff eff im i,k 1 1 D = D D (23)
eff im D and eff i,k D depend on the microstructure of the porous electrode (porosity, particle size and tortuosity) and on the operating conditions (temperature and pressure). The calculation of these parameters are described in Ref. [START_REF] Yakabe | Evaluation and modeling of performance of anode-supported solid oxide fuel cell[END_REF][START_REF] Chan | Polarization effects in electrolyte/electrode-supported solid oxide fuel cells[END_REF]. It should be noted that eff i D given in Eq. ( 23) is derived from the assumption of equimolar counter diffusion of reactant and product species. For binary component systems H 2 (1) -H 2 O(2), the effective diffusion for binary systems ( eff i,B D ) becomes:

1 1 eff i,B eff eff eff eff ij i,k 12 1,k 1 1 1 1 D = = D D D D (24)
The K nudsen diffusion coefficient ( i,k D ) and effective Knudsen diffusion coefficient ( eff i,k D ) for each gas specie [27] can be calculated by Eqs. ( 25) and [START_REF] Lehnert | Modeling of gas transport phenomena in SOFC anodes[END_REF].

i,k i i 2 8RT T D r 97.0r 3 M M ( 25
) eff i,k i,k D D (26)
The binary diffusion coefficient is calculated using the Chapman-Enskog equation [START_REF] Bird | Transport phenomena[END_REF].

2 ij 1/ 2 3/2 ij ij D i j T 1 1 D = 0.0018583 p M M ( 27
)
where ij i j = ( + )/2 is the collision diameter (in A ), p is total pressure (in atm), M i and M j are molecular weight of species i and j, ij D is the collision integral based on the Lennard -Jones potential [START_REF] Bird | Transport phenomena[END_REF], which can be obtained from Eq. ( 28).

)
where k is the Boltzmann constant which is equal to 1.38066 x 10 23 J K -1 and ij is the characteristic Lennard-Jones energy given by

1/2 ij i j ( ) .
The parameters used for calculation of D ij are listed in the Table 1. The effective binary diffusion coefficient through the porous electrode can be calculated from Eq. [START_REF] Tod | Thermodynamics and transport properties of gases for use in solid oxide fuel cell modeling[END_REF].

eff ij ij D = D (29)
Considering again Eq. ( 22) above, the second term on the right-hand side of Eq. ( 22) is called the permeation flux and is applied to take into account the effect of total pressure gradient on mass transport. However, it is assumed in this study that the total pressure is uniform over the entire depth of the porous anode. Thus, the definite total pressure gradient term can be ignored. Consequently, only the diffusion flux is addressed such that Eq. ( 22) reduces to:

n j i -i j i i eff eff j=1,j i i,k ij y N y N N dy P + = - RT dx D D (30)
For diffusion with heterogeneous chemical reactions, the flux ratios are governed by reaction stoichiometry. Summing Eq. ( 30) over n species leads to Graham's law of diffusion in gaseous mixtures [START_REF] Krishna | The Maxwell -Stefan approach to mass transfer[END_REF].

n ii i = 1 N M = 0 (31)
where M i is the molecular weight of component i. For two binary component systems (H 2 (1)-H 2 O(2) Eq. ( 30) becomes :

1 2 1 -1 2 1 eff eff i,k 12 N y N y N dy P + = - RT dx D D (32)
Rearranging Eq. (32) gives

2 1 2 1 1 1 eff eff i,k 12 y y N N dy 1 P N + = - RT dx D D (33)
Because y 2 = 1 -y 1 and

2 1 1 2
N N = -M /M (Graham's law), N 1 can be written as follows :

v + - = 0 dx 2FP dx D D D D D 2 2 2 2 2 2 2 2 1 2 2 2 1 1 2 2 2,k 2 n,a H H H H eff eff eff eff eff H O H H ,k H O H H H O H A J RT d y dy y 1 -y 1 v + - = 0 dx 2FP dx D D D D D 2 2 2 2 2 2 2 2 1 2 2 2 1 1 (36)
Equation ( 36) is in the form of an ordinary differential equation, which can be solved by involving the following initial conditions.

In Eqs. ( 37) and (38), by recalling Dalton's law of an ideal gas,

Eq. ( 16) can be rewritten again as, The term A V in Eq. (40) relates to the electrode porosity and the particle radius as in Eq. [START_REF] Choy | The development of intermediatetemperature solid oxide fuel cells for the next millennium[END_REF]. By combining Eqs. ( 9), ( 17) and (40) and thus varying the size of the electrode particles in the anode, the optimal size of the ion and electron conducting particles can be found.

Cathode side modeling

Similarly to the anode, the coupled electrochemical reactions and transport phenomena at the cathode can be determined as below: Because the flux of nitrogen does not participate only oxygen participates in the electrochemical reaction at the cathode. The concentration of O 2 can be determined by a self-diffusion mechanism [START_REF] Chan | A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness[END_REF], and can be expressed below.

For oxygen self-diffusion on the cathode side, For one-dimensional diffusion, The boundary condition for the transporting of O 2 Eq. ( 53) is 22 0 OO P = P at x = 0.

The cathode overpotential can be expressed as Eq. ( 54) [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF].

cat 00 t c i,c e,c i,c(x= t ) e,c (x=0) , = ( - ) -( - ) (54)

The electrolyte modeling

A dense electrolyte is used to conduct O 2-and typically is made from Yttrial-Stabilized Zirconia (YSZ). The overpotential of the dense electrolyte can be determined by Ohm's law,

elec total elec elec = J R t ( 55
)
where R elec and t elec are the resistivity (Ωm) and thickness (m) of the electrolyte respectively.

Parameters used in the model

SOFC modeling needs a number of parameters for calculation. We considered the model at steady state condition, with an operating pressure and temperature of 1 atm and 800 C. The supplied fuel compositions are H 2 and H 2 O which are oxidized by oxygen in the air. The dimensions of cell components and standard parameters used in the calculations are given in Table 2.

Mathematical model validation

To confirm the accountability of our numerical scheme, we validated our results with available experimental data [START_REF] Rogers | Validation and application of a CFD-based model for solid oxide fuel cells and stacks[END_REF]. The parameters used for validation are listed in Table 3. The predicted cell performance is compared with the experimental data found in the literature [START_REF] Rogers | Validation and application of a CFD-based model for solid oxide fuel cells and stacks[END_REF], and is shown in Fig 2 . The performance of the cell is predicted when the cell is supplied with 95% H 2 and 5% H 2 O, while operating at a temperature and pressure of 800°C and 1 atm, respectively. Oxygen in the ambient air is used as oxidant. It can be seen that the present model predictions agree well with the experimental results. Most parameters used in the model validation are obtained from Rogers et al. [START_REF] Rogers | Validation and application of a CFD-based model for solid oxide fuel cells and stacks[END_REF]. However, some parameters such as particle size, volumetric fraction of electronic and ionic conducting particle recommended by Hussain et al. [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] are used in place of unavailable parameters in Roger et al's work. The value of these parameters (particle sizes, volumetric fraction of electronic and ionic conducting particle etc.) recommended by Hussain et al. are typically used in other researches [16-22, 28, 30]. The value of tortuosity is varied to obtain the best agreement between the present model predictions and the experimental results as shown in Fig. 2. The typical tortuosity value for SOFC electrodes is in the range of 2 -6 [START_REF] Iwata | Performance analysis of planar-type unit SOFC considering current and temperature distributions[END_REF]. Hence, the tortuosity value of 2.75 used in the present model predictions is still in the typical range. Charge transfer coefficient, [START_REF] O'hayre | Fuel cell fundamentals[END_REF] 0.5 Conductivity of electronic conductor of anode, e,a (Sm -1 ) [START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF] 2e6

Conductivity of electronic conductor of cathode, e,c (Sm -1 ) [START_REF] Chan | Cathode micro model of solid oxide fuel cell[END_REF] 1e4

Conductivity of ionic conductor of anode, i,a (Sm -1 ) [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF] 3.44e4*(exp(-10300/T))

Conductivity of ionic conductor of cathode, i,c (Sm -1 ) [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF] 3.44e4*(exp(-10300/T)) Faraday constant, F (C mol -1 ) [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF] 96485

Results and Discussion

This numerical study was employed to gain an understanding of the coupled mass transfer and electrochemical phenomena in the electrodes of SOFCs. The present study reveals the working mechanisms of SOFC at the micro-scale level, while demonstrating the use of microstructure relationships to enhance the SOFC performance. Moreover, parametric analyses were performed to investigate the effects of operating parameters on SOFC performance. The optimal volumetric ratio between the ionic and electronic conducting particles is also proposed in this study. At present, the most favorable anode-supported configuration with thick anode and thin cathode/electrolyte is considered. The fuel, air composition, and the model parameters from the data available in the literatures used in the simulation are listed in Table 2. The simulation results are shown below in figs 3-8.

The effect of electrode microstructure parameters on the SOFC overpotentials

In this section, the electrode porosities and particle radius sizes are analyzed. From Eq. ( 9), it is indicated that the electrode porosity ( ) and particle radius size (r el , r io ) are mutually dependent. For simplicity in our study, we then defined that r io and r el are of equal radius. Firstly, we considered the effect of electrode porosity and its corresponding particle radius on the electrode overpotential. For each case, the volume fraction of electron conducting particles (υ) of both porous electrodes is equal to 0.6. Figure 3 shows the effect of porosities and corresponding particle radiipair on anode and cathode overpotentials. For the anode at this case, the optimal corresponding particle radius size is about 0.42 m at the an pair of about 0.59. Similarly for the cathode at this case the optimum corresponding particle size is about 0.43 m at the cat pair of about 0.59.

In general, the smaller the size of particle the larger the active surface area for electrochemical reaction of the electrode (A V ) (indicated by Eq. ( 2) and Eq. ( 7)). Therefore, reducing particle size can increase the electrode reactivity and thus reduce the electrode activation overpotential. However, when the particle size is too small, the reactant concentration at the reactive sites will approach zero and thus the concentration overpotential will increase considerably. Similarly, when considering electrode porosity, reducing porosity can reduce the ohmic overpotential as more solid particles are available to facilitate the transport of electronic and ionic charges. In addition, the activation overpotential is also reduced with decreasing porosity due to the increase in reactive . 0.75 Contact angle between e -and O 2-conducting particles, ( º ) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 15 Radius of electronic conducting particle, rel (μm) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 0.1 Radius of ionic conducting particle, rio (μm) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 0.1 Volume fraction of electronic conducting particle, υ [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 0.5 Volume fraction of ionic conducting particle, 1-υ [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 0.5 Reference H2 concentration, cH2,ref (mol m -3 ) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 10.78 Reference O2 concentration, cO2,ref (mol m -3 ) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 2.38 Reference exchange current density for H2 oxidation, [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 1320

2 H 0,ref J (A m -2 )
Reference exchange current density for O2 reduction, 2 O 0,ref J (A m -2 ) [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 400

Reaction order for H2 oxidation, γH2 [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 0.5 Reaction order for O2 reduction, γO2 [START_REF] Hussain | Mathematical modeling of planar solid oxide fuel cells[END_REF] 0.5 Source: Rogers et al. [START_REF] Rogers | Validation and application of a CFD-based model for solid oxide fuel cells and stacks[END_REF] Comparison of effect of electrode corresponding particle radius size on electrode overpotential when the electrode has different volume fraction of electron (v el ) and ion conducting particle (v io ) surface area. However, the concentration overpotential increases with decreasing porosity due to diminishing of space available for gas transport. The combined effects of porosity on ohmic, activation, and concentration overpotentials for minimization of the overall overpotential result in an optimal porosity for the electrode. The effect of increasing porosity when fixed solid fraction of the electrode means that the particle radius of the electrode decreases as shown in Eq. ( 9). This results in the reduction of the active surface area available for the electrochemical reaction as indicated by Eq. ( 2), so the reaction rate decreases resulting in an increase of the activation overpotential. Thus, the electrode porosities and corresponding particle radiiare mutually dependent in their effect on the overall SOFC overpotential. Therefore, minimizing the overall overpotentials resulting in an optimum of porosity and particle radius is very necessary and useful to enhance performance. Because the anode in an Anode-Supported SOFC is thicker than the cathode, the overpotential must be higher than the cathode, and the anode side will have more effect on cell performance than the cathode side for this kind of the SOFC.

Corresponding

Figure 4 shows effect of volume fraction of electron conducting particles (υ) on overpotentials. We considered the various fractions of υ at the various pairs of electrode porosities and corresponding particle radiuses. For simplicity, we defined r io as equal to r el . Because of our work, we considered an Anode -Supported SOFC, and mainly studied the anode side. Thus from Figure 4(a) and (b), for this case the results indicate that the minimum anode overpotential occurred at υ which is equal to 0.5 which occurs at an of about 0.57 at the corresponding pair of anode particle radius r an of about 0.43 μm. The result demonstrates the fact that the largest active area is achieved when the dimensions and volume fractions of the ion conducting and electron conducting particles are equal [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF]. When υ equals 0.5, the volume fraction of ion conducting particles (1-υ) becomes equal to the volume fraction of electron conducting particles resulting in the largest active area for electrochemical reaction and hence the best performance. Moreover, the effective ionic and electronic conductivities in the reaction zone are functions of the volume fraction of electron conducting particles; increasing υ beyond 0.5 increases the effective electronic conductivity but decreases the effective ionic conductivity, which results in increasing the ohmic overpotential and thereby decreasing the cell performance.

Figure 5 shows the distribution of H 2 molar fraction and H 2 O molar fraction in the anode region when the electrode properties of both anode and cathode are fixed at an and cat of about 0.59, υ of about 0.5, r of about 0.42 μm for anode and 0.43 μm for cathode; the other parameters for calculation are listed in Table 2. The concentration of H 2 decreases considerably from the anode surface towards the anode and electrolyte interface while the concentration of H 2 O increases considerably from the anode surface towards the anode and electrolyte interface, indicating an importance of gas transport in thick anodes.

Figure 6 shows a mapping relation between the anode porosity ( an ) and the anode corresponding particle radius (r an ) pairs that effect the anode activation overpotentials when all of the calculations are fixed and the design parameters are kept similar to the previous results listed in Table 2 except the value of υ on the anode and cathode side are at 0.5. From figure 6, it can be seen that the lower anode overpotential occurs at the anode porosity ( an ) of about 0.4 to 0.7 and the corresponding particle radius pairs ranges between 0.4 -0.5 μm. Ni et al. [START_REF] Ni | Micro-scale modeling of solid oxide fuel cells with micro-structurally graded electrodes[END_REF] indicated that the SOFC with nearly spherical particle electrodes have radii ranging from 0.012 to 5.0 m in order to achieve porosity ranges between 0.3 to 0.7. For an electrode porosity of less than 0.3, the extent of overlap between spherical particles increases with decreasing porosity, causing the active surface area (A v ) to decrease considerably with decreasing porosity [START_REF] Chen | Simulation of a composite cathode in solid oxide fuel cells[END_REF][START_REF] Ni | Micro-scale modeling of solid oxide fuel cells with micro-structurally graded electrodes[END_REF]. Conversely, the increasing of the porosity then increases the void fraction and decreases the solid fraction of the porous electrode resulting in the reduction of the active surface area (A v ) available for the electrochemical reaction. Moreover, the effective ionic and electron conductivities of the porous electrode decrease with the increase of porosity, which results in the increase of ohmic overpotential. Although the concentration overpotential decreases with the increase of porosity due to the increased mass transport rates the cell performance decreases due to increased ohmic overpotential with porosity. Furthermore, for a porosity higher than 0.7, the electrode will suffer from poor particle connectivity, poor percolation and thus poor performance and causing the collapsing of electrodes [START_REF] Ni | Micro-scale modeling of solid oxide fuel cells with micro-structurally graded electrodes[END_REF][START_REF] Lee | Viable image analyzing method to characterize the microstructure and the properties of the Ni/YSZ cermet anode of SOFC[END_REF]. Kenney and Karan [START_REF] Kenney | Engineering of microstructure and design of a planar porous composite SOFC cathode: a numerical analysis[END_REF] have indicated that the desirable porosity for high SOFC performance is in the range of 0.3-0.5 and also that an increase in porosity tends to smoothen the current density distribution. These observations can be explained by relating how porosity affects the effective ionic/electronic conductivities on the active reaction sites. The observed reduction in the average and maximum current density values with an increase in porosity can be attributed to the reduction in solid volume fraction, which results in less particle-particle contacts. This translates into a smaller number of active three phase boundary (TPB) for electrochemical reaction sites and, therefore, lowering the reaction rate for a given overpotential. Additionally, the increased porosity also reduces the effective electronic and ionic conductivity due to reduced solid fraction available for electron or ion transport. Further, an increase in porosity results in a facilitation of reactants diffusion, which minimizes the reactants transport limitation under the interconnect region and helps smoothen the current density distributions. Kenney and Karan [START_REF] Kenney | Engineering of microstructure and design of a planar porous composite SOFC cathode: a numerical analysis[END_REF] also reported that the best porosity to smoothen current density distribution is 0.5. Shi and Xue [START_REF] Shi | CFD analysis of a symmetrical planar SOFC with heterogeneous electrode properties[END_REF] show the relationship between the cell current density and the electrode porosity where the other parameters were fixed. Their results indicated that the high current density is in the porosity range of 0.5-0.8 and one can see that the maximum current densities are obtained when porosity reaches the values of about 0.75 under the pressure loss in the electrode of 1500 Pa. The extreme values of porosity are 0.2 and 0.9. These phenomena can be interpreted using the corresponding mass concentration distributions. When electrode porosity is 0.2, the permeability of the porous electrode reaches its minimum value and mass transport resistance reaches the highest value. As a result, the mass transport through the porous electrodes is difficult even with increased pressure. When the porosity increases to 0.9, the permeability of the porous electrode reaches its highest value and the mass transport resistance reaches the lowest value. However, the reactive surface area decreases due to the high electrode porosity leading to a low current density in this case. Therefore, for the electrodes consisting of spherical particles, porosity between 0.3 and 0.7 is regarded as a practical range for good performance and good particle connectivity, but for a better performance of each cell the porosity must be matched with each corresponding pair of electrode particle radius size (indicated by Eq.( 9)). Previous work did not regard the porosity and corresponding particle radius pairs [START_REF] Costamagna | Micro-modeling of solid oxide fuel cell electrodes[END_REF][START_REF] Deseure | Theoretical optimization of a SOFC composite cathode[END_REF][START_REF] Ni | Micro-scale modeling of solid oxide fuel cells with micro-structurally graded electrodes[END_REF][START_REF] Chan | Anode micro model of solid oxide fuel cell[END_REF][START_REF] Chan | Cathode micro model of solid oxide fuel cell[END_REF][START_REF] Chan | A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness[END_REF][START_REF] Kenney | Engineering of microstructure and design of a planar porous composite SOFC cathode: a numerical analysis[END_REF][START_REF] Shi | CFD analysis of a symmetrical planar SOFC with heterogeneous electrode properties[END_REF]. They showed only that the performance i.e. current density and power density of the cell depended only on the particle radius sizes or on the electrode porosity. The results from our work are consistent with the practical range of both electrode porosity and the particle radius, but our results regarded the relational pairs of the electrode porosity and the corresponding pairs of particle radius sizes. As such our results can be effective and useful for design and optimization of the electrode structure for maximizing cell performance.

The effect of operating parameters on the SOFC overpotentials

Figure 7(a) and (b) show the effect of operating temperatures on the anode overpotentials at various corresponding electrode particle radius sizes and at various electrode porosities, respectively. The electrode porosity and particle radius relationships are indicated by Eq. ( 9). For all calculations, the electrode properties, both of anode and cathode, are fixed at υ of about 0.5. Again, the operating and other design parameters are kept similar to previous results listed in table 2, except the operating temperature (T), and υ. Figure 7(a) and (b), considered the anode side of the cell because we were interested in the anode supported SOFC that the anode has more effect than the cathode.

From the results Figure 7(a) and (b), it can be seen that reducing the temperature can considerably increase the cell overpotential. Decreasing temperature can result in reduction in reference exchange current density and the ion conductivity of electrolyte, which in turn increases the ohmic and activation overpotentials. Therefore, the overpotentials of electrodes and the ohmic overpotential of electrolyte increase with decreasing temperature. These phenomena indicate that a high temperature is desirable for SOFC from the viewpoint of energy efficiency. However, high temperature limits material selection and also causes problems associated with SOFC long-term stability, such as thermal mismatch and detachment of electrode layers from the electrolyte. It is important, therefore, to develop new materials with high ionic conductivity and high catalytic activity at intermediate temperatures. For higher performance we could utilize the optimum relationships of electrode microstructure parameters to enhance cell performance at higher temperatures, or develop new materials for utilizing cells at lower temperatures. Figure 8(a) and (b) show the effect of operating pressure on the anode overpotential at various corresponding electrode particle radius sizes and the effect of operating pressure on the anode overpotential at various electrode porosities, respectively. The electrode porosity and particle radius relationships follow Eq. ( 9). For all calculations, the electrode properties, both of anode and cathode, are fixed at υ of about 0.5. Again the operating and other design parameters are kept similar to the previous results that are listed in table 2, except the operating pressure (P) and υ. From both results, it can be seen that at high pressure, the electrode overpotential is lower than at low pressure. With the increase of pressure, the reactant concentration at the reaction sites increases. This in turn enhances the rate of electrochemical reaction and the rate of mass transport resulting in the minimization of electrode overpotentials, and hence better performance. However, increase in pressure results in other problems, such as limitation on material selection, gas sealing, and mechanical strength of the cell components [START_REF] O'hayre | Fuel cell fundamentals[END_REF]. The modeling analysis shows that reducing the particle's size increases the volumetric active surface area, but adversely impedes the gas transport. There exist optimum particle sizes for both electrodes to achieve minimum overpotentials. Similar to the particle size effect, reducing electrode porosity has the opposite effect on activation and concentration overpotentials, and the ohmic overpotential can be decreased with decreasing porosity as more particles are available for electronic/ionic charge conduction. The combined effect results in a monotonic decrease in electrode overpotential with decreasing electrode porosity. From the relationship indicated by Eq. ( 9), the porosity and particle radius are mutually dependent on them for better cell performance, not only dependent on the particle radius or the porosity, but on the optimum pair between them. In our study, the optimum pair is found in a wider range that provides the cell performance closely to the maximum cell capability. Moreover, the optimal values of volumetric fractions of electronic (υ) and ionic (1-υ) conducting particles in the electrode are also proposed in our work. In our study, for using the optimum value, it can enhance the cell performance and is best with the optimum pair of the electrode porosity and the corresponding electrode particle radius. In addition, the study results substantiate the fact that SOFC overpotential could be effectively decreased by increasing the operating temperature, as well as the operating pressure. Therefore, the concept of micro-model SOFC is feasible and has great potential to improve the SOFC performance by simply modifying the electrode microstructures.

Conclusions

A micro-scale mathematical model of a Solid Oxide Fuel Cell (SOFC) has been developed to simulate the coupled electrochemical reactions and mass transfer involved in the SOFC operation. This model was developed for modeling the performance of SOFC with spherical-shaped particle electrodes. The important features of this model aim to determine the optimum relationships between porosity and electrode particle radius size or as it is called, the corresponding electrode particle radius, due to these two parameters being mutually dependent. The optimum pair of corresponding electrode particle radius can minimize the cell overpotentials. Apart from the optimum relationship of electrode particles size and porosity, this model is optimized for the electrode solid fraction between the volumetric fraction of ionic and electronic conducting particles. Furthermore, optimal values of the volumetric fraction of ionic and electronic particles were found to minimize electrode total overpotentials. At the optimal volumetric fraction of ionic and electronic conducting particles of the electrodes, parametric analyses were performed to investigate the effects of operating parameters on SOFC overpotentials. The results substantiated the fact that SOFC overpotentials could be effectively decreased by increasing the operating temperature and increasing the operating pressure. The present study reveals the working mechanisms of SOFC at the micro-scale level and demonstrates the promise of the use of microstructure relationships to enhance the SOFC performance by simply modifying the electrode. The model-predicted results are validated with the experimental data found in the literature. An excellent agreement is obtained between the predicted values and the measured data.

Nomenclature

A v Electrochemically active surface area per unit volume (m 2 m -3 ) D ij Binary diffusion coefficient of the species i (m 2 s -1 ) D i,k Knudsen diffusion coefficient of the species i (m 2 s -1 ) F Faraday constant 96485 (C mol -1 ) J e Electronic current density (A m -2 ) J i Ionic current density (A m -2 ) J n Transfer current per unit area of reactive surface (A m

Figure 1 .

 1 (a) Model configuration of a SOFC anode; (b) Electrochemical reaction at three phase boundaries (TPB) inside the anode.

  anode, the boundary conditions for Eqs. (42) and (43) are:

Figure 2 .

 2 Figure 2. The comparison between the model predictions and experimental results of Rogers et al. [32].

Figure 3 .

 3 Figure 3. Corresponding particle radii and porosities effect on anode (a,b) and cathode overpotentials (c,d).surface area. However, the concentration overpotential increases with decreasing porosity due to diminishing of space available for gas transport. The combined effects of porosity on ohmic, activation, and concentration overpotentials for minimization of the overall overpotential result in an optimal porosity for the electrode.The effect of increasing porosity when fixed solid fraction of the electrode means that the particle radius of the electrode decreases as shown in Eq. (9). This results in the reduction of the active surface area available for the electrochemical reaction as indicated by Eq. (2), so the reaction rate decreases resulting in an increase of the activation overpotential. Thus, the electrode

  Comparison of effect of electrode corresponding particle radius size on electrode overpotential when the electrode has different volume fraction of electron (v el ) and ion conducting particle (v io )Corresponding particle radiuses, r an ( m.)

Figure 4 .Figure 5 .

 45 Figure 4. Anode corresponding particle radii and porosities effect on anode overpotentials at various volume fractions of electron (υ); (a) corresponding particle radii versus anode overpotentials; (b) porosities versus anode overpotentials.

Figure 6 .

 6 Figure 6. Effect of anode porosity and corresponding anode particle radius on anode overpotential when the volume fraction of the electron conducting particle in both of electrodes (υ) are equal to 0.5.

Figure 7 .

 7 Figure 7. Effect of operating temperature on SOFC overpotentials; (a) corresponding particle radius versus anode overpotentials; (b) porosities versus anode overpotentials. parameters to enhance cell performance at higher temperatures, or develop new materials for utilizing cells at lower temperatures.Figure8(a) and (b) show the effect of operating pressure on the anode overpotential at various corresponding electrode particle radius sizes and the effect of operating pressure on the anode overpotential at various electrode porosities, respectively. The electrode porosity and particle radius relationships follow Eq. (9). For all calculations, the electrode properties, both of anode and cathode, are fixed at υ of about 0.5. Again the operating and other design parameters are kept similar to the previous results that are listed in table 2, except the operating pressure (P) and υ. From both results, it can be seen that at high pressure, the electrode overpotential is lower than at low pressure. With the increase of pressure, the reactant concentration at the reaction sites increases. This in turn enhances the rate of electrochemical reaction and the rate of mass transport resulting in the minimization of electrode overpotentials, and hence better performance. However, increase in pressure results in other problems, such as limitation on material selection, gas sealing, and mechanical strength of the cell components[START_REF] O'hayre | Fuel cell fundamentals[END_REF].

Figure 8 .

 8 Comparison of effect of cell operating pressure on electrode overpotential at varius corresponding electrode particle radius sizeCorresponding particle radiuses, r an ( m.) Effect of operating pressure on SOFC overpotentials; (a) corresponding particle radius versus anode overpotentials; (b) porosities versus anode overpotentials.

J

  Reference exchange current density of the electrode (A m -2 ) H 2 ref J Reference exchange current density for H 2 oxidation (A m -2 ) O 2 ref J Reference exchange current density for O 2 reduction (A m -2 ) M Molecular weight (kg kmol -1 ) N Molar flux (mol m -2 s -1 )

  are related to each other as shown here,

	el n =	3 [ +((1-)/(r /r ) )] io el	and io n = 1 n	el	(3)
	el Z = 3+	[n	el	2 (1 n )(r / r ) ] el io el Z-3	and	io Z = 3+	[n	el	2 el io el io el (Z-3)(r /r ) (1 n )(r / r ) ] 2	(4)
	el p = 1-4-Z el -el	2.5	0.4				and	io p = 1-4-Z io -io	2.5	0.4

Table 1 .

 1 List of Lennard-Jones potentials used in calculation of D ij (taken from Tod et al.,[START_REF] Tod | Thermodynamics and transport properties of gases for use in solid oxide fuel cell modeling[END_REF]).

	Gas species CH4	H2O	CO2	CO	H2	N2	O2
	a i	3.758 2.641 3.941 3.69 2.827 3.798 3.467
	a i	148.6 809.1 195.2 91.7	59.7	71.4 106.7

  , and the parameters 1 M and 2 M are the molecular weight of species 1 and 2 respectively.

									1
	1 N = -	P RT	1 1-y eff 12 DD 1,k eff 1 +	1 dy dx		(34)
	where 2 H N = -= 1-M /M For two binary component systems (H 2 (1)-H 2 O(2) Eq. (34) becomes 12 22 2 2 2 1 HH eff eff H o H H ,k 1-y dy P1 + RT DD dx (35)
	where	22 H O = 1-M /M H	. The parameters	2 M , and H	M	2 HO	are the molecular weight of H 2 and H 2 O, respectively. The
	D	2 eff H ,k	and	D	22 eff H o H	are the effective Knudsen diffus tituting Eq. (35) into Eq. (19), the governing equation of H 2 concentration
	becomes
	d y	H						1 -y	H	2	dy	H	n,a A J RT	1	y	H
						eff		eff		eff		eff	eff
						H O H	H ,k		H O H	2	H	2,k	H O H	2

Table 2 .

 2 Standard parameters used in this study.

	Parameters	Values
	Operating temperature, T ( C)	800
	Operating pressure, P (atm)	1
	H2 and H2O in fuel mixture inlet mole fraction	0.95, 0.05
	O2 and N2 in air inlet mole fraction	0.21, 0.79
	Anode thickness, ta ( m)	2000
	Anode tortuosity, a [21]	4.5
	Anode pore radius (rp) (μm.)	0.5
	Cathode thickness, tc ( m)	60
	Cathode tortuosity, c [21]	4.5
	Cathode pore radius, rp (μm.)	0.5
	Electrolyte thickness, telec ( m)	60
	Contact angle between e -and O 2-conducting particles, (º) [16, 19]	15
	Reference anode exchange current density, ref 0,a J (A m -2 ) [19]	400
	Reference cathode exchange current density, ref 0,c J (A/m -2 ) [19]	1320
	Total average particle coordination number in electrode, Z [16,19]	6

Table 3 .

 3 Parameters used for model validation.

	Operating temperature, T ( C)	800
	Total pressure, P (atm)	1
	Fuel composition, xH2, xH2O	0.95, 0.05
	Air composition, xO2, xN2	0.21, 0.79
	Anode conductivity, a (S/m)	71428.57
	Cathode conductivity, σc (S/m)	5376.34
	Electrolyte conductivity, κelec (S/m)	0.64
	Anode electrode layer thickness, ta ( m)	1000.0
	Cathode electrode layer thickness, tc ( m)	50.0
	Electrolyte thickness, te ( m)	10.0
	Porosity of anode and cathode, ε	0.375
	Tortuosity of anode and cathode, τ	2.75
	Pore radius of anode and cathode, rp (μm.)

Acknowledgements

Authors would like to acknowledge with gratitude the financial support from the National Research University, the Thailand Higher Education Commission and the Thailand Research Fund.

n el Number fraction of electron conducting particle n io Number fraction of ion conducting particle n t Total number of particles per unit volume P Pressure (atm) p el Probabilities for an electron conducting particle to belong to connecting ends of the composite p io Probabilities for an ion conducting particle to belong to connecting ends of the composite R i Reaction rate (mol m -3 s -1 ) r Mean value of the pore radii (m) r el Radius of e -conducting particles ( m) r io Radius of O


Keywords: JEL Classification : D43, L13, L22 vertical integration, vertical equilibrium, industry growth, successive Cournot oligopoly, double mark-up effect

This paper illustrates the effect of market size on the decision of whether or not firms should vertically integrate or disintegrate. We use a model of two successive stages of production with Cournot competition in each stage. In this model, firms choose to specialize (either upstream or downstream) or to integrate the two stages, before making their production decisions. The decision of whether or not to integrate or specialize depend on the trade-off between "escaping from" the double marginalization problem or the gain from specializing on the production stage in which the firm is more efficient. We show (using simulations) that more firms choose to be vertically integrated as the valuation of the final product or the number of consumers increase, unless the number of firms increases proportionately.

Introduction

Existing economic theories have provided little analysis of vertical disintegration while focusing on integration. This is surprising as Adam Smith indirectly mentioned the issue in his analysis of the division of labor. Actually, the few contributions to vertical disintegration are built on Adam Smith's proposition that the division of labor is limited by the extent of the market [START_REF] Young | Increasing returns and economic progress[END_REF], [START_REF] Pigou | The Economics of Welfare[END_REF], [START_REF] Stigler | The division of labor is limited by the extent of the market[END_REF]). In applying Smith's theorem, [START_REF] Stigler | The division of labor is limited by the extent of the market[END_REF] argued that vertical disintegration is the typical evolution of a growing industry, whereas a declining industry must be characterized by vertical integration. Stigler offers a testable proposition, regarding industry growth and vertical integration. His vertical separation argument is that as the market demand for an industry's product increases, it is more economic for the firms in the industry to purchase some of their inputs from specialized firms, who can reap economies of scale from large scale operation.

Conversely, as the industry demand falls, firms will produce in-house inputs that were formerly outsourced. In other words, vertical disintegration is the typical evolution of an industry in growth, while vertical integration that of an industry in decline.

The following empirical research is investigating the validity of Stigler's proposition. Tucker and Wilder (1977) examined 54 American manufacturing firms, and [START_REF] Levy | Testing stigler's interpretation of 'The division of labor is limited by the extent of the market[END_REF] explored census data for 38 industries from 1963, 1967 and 1972 : they found some support to Stigler's conjecture. Subsequently, [START_REF] Wright | Vertical disintegration and the life cycle of firms and industries[END_REF] tested Stigler's hypothesis with data on 407 investment withdrawals in the UK between 1977 and 1979. They showed that vertical disintegration is positively correlated with industry growth. By contrast, [START_REF] Stuckey | Vertical Integration and Joint Ventures in the Aluminium Industry[END_REF] found opposite results from the study of the aluminum industry. Thus, the results are inconclusive.

In this paper, we construct a model of vertical equilibrium allowing theoretical examination of Stigler's hypothesis, within a framework of Cournot competition. We examine equilibria where integrated and specialized firms coexist, and we confront the vertical equilibrium to modifications of the market size.

In this context of successive Cournot oligopolies, [START_REF] Salinger | Vertical mergers and market foreclosure[END_REF] shows, in a two stages model with fixed proportions, that vertical integration can result from the double mark-up effect. But he simply imposes the coexistence of integrated and specialized firms, without studying the integration game. [START_REF] Gaudet | Upstream-downstream specialization by integrated firms in a partially integrated industry[END_REF] do examine the vertical equilibrium in a model of successive Cournot oligopolies and show that complete integration emerges in most circumstances.

Depending on the number of firms at each stage, they show that there are few complete disintegration equilibria, and that only one specific configuration leads to the realistic coexistence of integrated and specialized firms. However, they don't introduce the motivations of vertical disintegration1 : we believe that explains why vertical integration appears to be dominant in these models. So we do assume that there exist economies from specialization of the firms.

Actually, as pointed out by [START_REF] White | The division of labor is limited by the standardization in the market : a note on Stigler's vertical integration theme[END_REF], the dynamics of vertical integrationdisintegration cannot derive exclusively from the existence of economies of scale, contrary to what was suggested by [START_REF] Stigler | The division of labor is limited by the extent of the market[END_REF]. If there are economies of scale for the production of an input, these cost savings could also be realized internally by the integrated firm. So, vertical disintegration comes from economies from performing a limited set of tasks, as opposed to economies from repeatedly performing any one activity. [START_REF] Perry | Vertical equilibrium in a competitive input market[END_REF] do introduce economies from specialization, in a two stage model of vertical equilibrium with fixed proportions. These economies come from final demand fluctuations which tend to increase profits of specialized firms relatively to integrated firms. He characterized a vertical equilibirum where upstream, downstream and integrated firms coexist. However, by explaining vertical integration from economies of synchronization, this model doesn't allow to examine the relationship between vertical equilibrium and industry growth.

Finally, [START_REF] Perry | Vertical integration and industry growth[END_REF] propose a vertical equilibrium model where vertical disintegration comes from differences in the firm's competencies. In their model, each firm is endowed with a separate cost function for each stage, and the firms which are more efficient in upstream production are less efficient for the downstream production. They assume perfect competition downstream and imperfect competition at the upstream stage, which explains the benefits of vertical integration. For a given linear final demand, a vertical equilibrium can be defined, for which the firms make a choice between upstream or downstream specialization, and confronted to shifts in the demand size. They conclude that Stigler's hypothesis was verified only in specific circumstances.

Our paper is related to some extent to that of [START_REF] Perry | Vertical equilibrium in a competitive input market[END_REF] and [START_REF] Perry | Vertical integration and industry growth[END_REF]. We also introduce economies of specialization, to explain the coexistence of vertically integrated and vertically disintegrated firms 2 . We follow [START_REF] Perry | Vertical integration and industry growth[END_REF] for the main hypothesis of a continuum of agents with different skills.

However, their assumption on competition is not well founded ; because of the limited total number of firms, it seems to be more logical to assume imperfect competition on each stage like in [START_REF] Gaudet | Vertical integration, foreclosure and profits in the presence of double marginalization[END_REF]. We also pose a simple linear final demand, and consider that firms compete in quantities after having decided whether to integrate both stages or not.

The model considers the presence of an exogenous number n of firms that play a two stage game. In the first stage, firms decide whether to specialise in either the upstream or the downstream production, or to be vertically integrated. In the second stage, upstream and downstream production decisions are taken by the firms that compete in quantities at both stages of production. The rest of the paper is organized as follows. Section 2 describes the model. Assuming a linear final demand, section 4 characterizes the vertical equilibrium. This vertical equilibrium depends on horizontal oligopoly equilibrium at each stage of the process (section 3). Section 5 tests Stigler's hypothesis through a simulation while analyzing the impact of a modification in the size of demand on the proportion of vertically integrated firms. Section 6 concludes.

The Model

The model comprises two successive stages of production (an intermediate and a final one) and three types of producers : integrated, downstream and upstream firms. An integrated firm produces both the final good and the intermediate good, with the latter used as an input in the production of the final good. An upstream firm produces only the intermediate good, while a downstream firm produces the final good and buys the intermediate good on the open market. The firms produce with fixed-coefficient technology : one unit of the intermediate good is required to produce one unit of the final good. Lastly, each firm can choose to operate in one or both of the two stages.

The firms differ in their efficiency in producing intermediate and final goods.

We assume this not only because it is usually observed in reality (see [START_REF] Gaudet | Upstream-downstream specialization by integrated firms in a partially integrated industry[END_REF] for the world oil industry), but also because skills and specialization choices are closely correlated and because these differences of skills explain diseconomies of vertical scope. Following [START_REF] Perry | Vertical integration and industry growth[END_REF], we assume that firms that are more efficient in downstream production are less efficient in upstream production, and vice versa. Efficiency is defined as lower total costs for producing all outputs levels. We pose, as they do, a cost parameter θ specific for every firm, where 0 , such that firms with a 1 ≤ ≤ θ θ near one have a relative advantage in upstream technology, whereas firms with a small ϑ have a relative advantage in downstream technology. Industry is composed of an exogenous number of firms, . If is not too large, every stage of the production process is characterized by an oligopoly : we assume Cournot competition. It is not possible in our model to permit free entry or exit of firms, nor horizontal integration, because only the most efficient firms would exist in equilibrium. We consider that the firms, each characterized by a specific Every firm faces constant marginal costs and fixed costs (decreasing average costs). This assumption is not critical to our model and is compatible with the conjecture of Cournot competition. For the purposes of our research, we assume that θ affects fixed costs. Had we assumed that θ affects marginal costs the results of our story would have been the same, but the equations would have been much more complicated. So we consider that all firms have the same marginal cost for each stage, the difference in efficiency being reflected in fixed cost3 .

θ θ θ χ θ ϕ y y C u + = / ) ( (2)
We integrate neither economies nor diseconomies of vertical integration : both exist, but they are likely to compensate each other in practice. So we consider the total cost function of the integrated firms is the simple algebraic sum of the upstream and downstream costs. Thus, if is the output of the final goods of an integrated firm θ x θ , its total cost is :

θ θ θ χ θ ϕ θ x c f x C i ). ( / . ) ( + + + = (3)
On the demand side, we assume that the final consumers , with exogenous total number i I , all behave the same way and have a linear individual demand function :

p β α q i - =
With y the aggregate final output of specialized firms and x that of integrated firms, the aggregate inverse demand function is :

) I β y x ( β / p(x,y) + - = α (4)
or, with

β α / = a and ) . /( 1 I b β = , y) b(x a p(x,y) + - = (4bis)
Growth in demand must mainly be captured by a rise in the number I of consumers, associated with a fall of . A rise of (and of b a α ) corresponds with an increase in the value granted by each consumer to the final good5 . It can therefore be regarded as a form of (indirect) growth of the market size. Modifying β raises no interesting question because it impacts simultaneously a and b.

To formalize Stigler's idea, we work on a vertical equilibrium which is defined by ( θ θ , ), with 0 < θ < θ < 1. The firms with an index θ such as θ θ < ≤ 0 produce only the intermediate goods, the firms with an index θ such as

θ θ θ < ≤
are integrated, and the others are firms specialized in the downstream production. Theoretically this equilibrium can be defined if the oligopoly markups are significant. This will result of the firm's production choice (based on their relative profits) which depends on their index of efficiency. Because of the previous assumptions, the profits of the downstream firms must decrease when θ increases and the profits of the upstream firms must increase. A firm highly qualified for the downstream (upstream) production will maximize its profit while specializing according to its comparative advantage although it undergoes the upstream (downstream) oligopoly mark-up. If a firm is fairly effective with both technologies (θ around 0.5), profit maximization can lead to an integration of the two stages, giving up the benefit of specialization (weak in it's situation) so as to avoid the double marginalization. To determine this vertical equilibrium, the horizontal equilibrium is characterized at each stage in section 3.

. Horizontal Equilibrium

The Cournot horizontal equilibrium determines the profit-maximizing outputs for integrated firms, downstream firms and upstream firms. Let the number of firms and the demand size be such as : i) the final equilibrium price p* exceeds the marginal cost (which is also the variable average cost) of integrated firm ( ).

Otherwise, some firms would not be profitable.

OPTIMUM FOR INTEGRATED FIRMS

The integrated firms must choose between two possibilities:

-To withdraw intermediate market and to use all their intermediate production for their final production.

-To take part in the intermediate market, by selling intermediate goods to the downstream firms, and/or buying these goods to the upstream firms.

In our model, this choice depends on many circumstances. For an integrated firm, selling intermediate goods to the downstream firms can increase the scale of its upstream production, and thus reduce the upstream average cost (because of existence of a fixed cost). But, it introduces more competition upstream, and therefore lowers the intermediate price addressed to the downstream specialists who are direct competitors of integrated firms. Buying intermediate goods to upstream specialists presents at opposite the advantage of increasing the price charged by the upstream specialists and thus of degrading the competitiveness of downstream firms compares to integrated firms. However, that increases the provisioning price of the integrated firm (since it can produce at a marginal cost lower than the upstream market price),

and that reduces the level of their in-house upstream production which increases the average upstream cost.

Gaudet and Long [START_REF] Gaudet | Vertical integration, foreclosure and profits in the presence of double marginalization[END_REF] show that no general answer can be brought concerning these problems, even if one supposes very simple cost functions (constant marginal cost and not of fixed cost) and absence of competence differentiation. But according to [START_REF] Salinger | Vertical mergers and market foreclosure[END_REF], it is especially in the presence of increasing average cost that the integrated firms can have sometimes interest to take part in the intermediate market.

So let us pose that the integrated firms use all and only their intermediate output for their final production. If is the production of the integrated firm θ and if is the total production of the other integrated firms (so we have ),

θ x 0 x x x x = + 0 θ
then the profit function of the firm θ can be expressed as :

ϕ θ χ π θ θ θ θ - - + - + = f x c x y x x p i ). ( ).
, ( 0The output firm θ sets is :

y x b c a x - - - - = / ) ( χ θ (5)
Firms have the information on the distribution of firms on [0,1]. So they know that the number of integrated firms is simply )

.( θ θn , and that firms differ only in their fixed costs. This means that they know that each integrated firm θ will produce the same output )

.

( θ θ θ - = n x x . So it comes : b c a y n ) ( ) ) .( 1 ( 1 χ θ θ - - = +       - + x (6)
This gives the output supplied by integrated firms given the output of specialized firms and the vertical equilibrium. 6) and ( 10), simultaneously define the equilibrium outputs and given the vertical boundaries( * x * y θ θ , ). Let's normalize to unity the cost parameters c and χ , which raise no interesting questions and don't impact our conclusions.

] 1 )) 1 ( ² ( ² [ ) ( )) 1 ( 1 ( ) 2 ( * n n b n a n x + + - + - - + -f n y r p d . θ θ θ π d (17)
A full vertical equilibrium exists if 0 < θ * < * θ < 1. There can be no integrated firm (θ * = * θ ) or only integrated firms in equilibrium (θ * = 0 et * θ = 1), but these cases are not interesting for our problem which is to test Stigler's idea.

Finally, to define the vertical equilibrium, we need to solve the simultaneous equations system (4), ( 8), ( 11), ( 12), ( 16) and ( 17). Analytical comparative statics on the vertical equilibrium are intractable, so we do numerical comparative statics over a wide range of values for market size parameters.

Vertical equilibria and industry growth

The purpose of this model of vertical equilibrium is to test the Stigler's hypothesis according to which growing industries would be characterized by vertical disintegration, and, conversely, declining industries would be characterized by vertical integration. We will thus evaluate the impact on the vertical equilibrium of modifications in the value of the model's parameters that are correlated with market size. Industry growth can of course be captured by changing the demand parameters a and b. However, on the supply size, it is important to evaluate the impact of an increase in the total number of firms n that can come with demand growth.

We focus only upon the cases where there is a full vertical equilibrium, i.e

where three types of firms exist. We consider consistent (relative at the marginal costs and at the demand size) values for fixed costs 7 : we posit f = 40 and ϕ = 10.

INCREASE IN THE NUMBER OF FIRMS

The increase in the number of firms in the industry is accompanied by a fall of the relative and absolute number of integrated firms. The integrated stage contracts from above and below (see table 1).

A rise in the total number of firms generates an additional competition on each stage that lowers oligopoly markups and increases supplied quantities8 for the three types of firms (with Cournot conjecture). The subsequent increase in intermediate demand works in the opposite direction for upstream markups but, in our model, the former effect is stronger than the latter. We checked by simulation that the final and intermediate prices drop and that total output increases with a growing number of firms. The integrated firms close to old equilibrium boundary θ * specialize downstream (rise in θ ) to benefit from the fall in intermediate price. Those close to old equilibrium boundary θ * give up downstream technology (decrease in θ ) 6 We can well observe a rise of π u and a fall of π d when θ increase ; the simulation shows us that π i is higher than π u and π d for firms with θ near 0.5 for a wide range of the parameters. 7 A rise in the fixed costs for downstream technology leads to a vertical equilibrium for which the number of downstream firms decreases but especially for which the number of integrated firms very strongly falls ; this can be explained by the fact that the integrated firms have a higher fixed cost for downstream technology (because their index θ is higher). A rise in ϕ affect identically integrated firms and upstream firms (thus θ does not change) but does not affect downstream firms, which increases the incentive to specialize downstream (rise of θ ).

because of the decreasing downstream markup and of the increased intermediate demand. Integrated stage disappears if total number of firms becomes large enough. So with the additional competition on each stage when the number of firms rise, it becomes less advantageous to be vertically integrated to escape double marginalization. The firms will have a growing interest to specialize according to "their comparative advantage" in terms of skills. This is an interesting preliminary result, but an increase in the number of firms cannot be regarded as a reliable indicator of the market size. It is the rise in the number of consumers and/or in the value which they grant to the product, which leads to industry growth.

GROWTH IN THE NUMBER OF CONSUMERS

By fixing the demand parameters β and α , and the number N of firms, we examine industry growth which occurs by increasing number of consumers I (Table 2).

With the rise of the number of consumers, the slopes of the direct and derived inverse demand curves decrease, the produced quantities x and y increase in the same proportion, but the prices p and r are not modified firstly (because demand is linear).

Thus, the rise in the proportion of integrated firms does not come initially from an increase in the upstream oligopoly margin. If we look, on the other hand, at the profit functions of the three types of firms, we observe that a increase in the number of consumers induces a stronger rise of profit for integrated firms than for upstream and downstream firms. This results of the higher fixed costs in the integrated firms than in specialized firms. Actually, if the production scale increases and when marginal costs are constant, the profit increase is stronger for firms with higher fixed costs. If it is this scale effect which explains the initial increase in the proportion of integrated firms, this tendency induces modifications of the competing situation for each stage that will impact on final equilibrium. Indeed, the subsequent fall in the number of upstream firms is accompanied by a rise in the oligopoly markup and thus induces downstream firms near to the old equilibrium boundary θ * to integrate upstream technology. But on the other hand, the decrease in the number of independent downstream firms (accentuated by the former effect) tends to lower their intermediate demand. This moderates the increase in the oligopoly price. We verified by simulation, that the first effect is stronger than the second.

Proposition 1 : A growing number of consumers leads to an increasing proportion of firms choosing to be vertically integrated.

This result comes however partly from the scale effect. To cancel this effect and to focus on the strategic behaviors, we next consider that the number of consumers and the number of firms increase proportionately.

PROPORTIONAL GROWTH IN THE NUMBER OF CONSUMERS AND FIRMS

The results are reversed if the number of firms increases proportionately with the growth in the number of consumers. We do it by fixing n/I as I increases, which implies a constant average individual output for each firm 9 (Table III).

The decrease occurs both from integrated firms specializing upstream and downstream 10 . This means, in first analysis, that the impact from increasing the number of firms (which lowers markups) overrides the impact of increasing the number of consumers. This can be interpreted as follows. By assuming a proportional growth, "the scale effect" for each firm is cancelled, so the "competition effect" dominates : the rise of the industry size induces a rise of competition, therefore a decrease in the markups (despite the demand growth) and finally a growing specialization (see 5.1). Upstream F (share) 54% 47% 40% 34% 23% 14% 5%

The explanation of this correlation is partly the same as in the case of a rise of the number of consumers I (high fixed cost for integrated firms) but not only. A rise of α results here indeed in an increased final price, and also in an increased intermediate price, because oligopoly markups increase on each stage. Downstream firms benefit as well as the integrated firms from the growth in the value to consumers through the rise in the final price and quantities ; but they also face the rise in the intermediate price. Thus, they will not benefit as much as the integrated firms from the rise of α.

Those that are close to θ will be incited to integrate the upstream stage to catch the markup. This induces a decrease in the derived demand for the intermediate input which make it profitable for upstream firms near θ to integrate forward. This tendency is reinforced by the fact that upstream firms benefit only partly, like the downstream firms, of the rise of α. They can catch downstream markup while integrating forward.

Proposition 3 : A growth in the value the consumers grant to the final product results in an increasing proportion of firms choosing to be vertically integrated.

A real problem of double marginalization appears here : specialized firms on each stage choose prices and quantities without considering the impact on the profits of the other stage firms. By eliminating the double marginalization, vertical integration greatly modifies the structure of the intermediate market and thus the margins of the specialized firms.

Conclusion

This paper has examined, in a model of successive Cournot oligopolies, how a modification of the market size could impact on the degree of vertical integration in an industry. This model differs from most of vertical equilibrium models because it introduces explicitly the benefits of vertical disintegration which derive from differences between the firms in their efficiency in producing intermediate and final goods. Vertical integration is motivated by the double mark-up effect. Given a linear final demand, we can then define a vertical equilibrium where firms choose to operate on one or both of the two stages according to their skills. This equilibrium is affected by industry growth. If industry growth is defined as an increase in the value granted to the final good, then the net impact of growth is a higher proportion of vertically integrated firms, because of the subsequent higher mark-ups on each stage. The net impact is the same if industry growth is defined as an increase in the number of consumers, because integrated firms with higher fixed costs benefit more from the rising production scale. However, the result is reversed if we suppose, in order to cancel the scale effect, a proportional increase in the number of firms and consumers.

A lower proportion of firms choose to integrate both stages, because of the diminishing impact of double mark-up effect when the number of firms increases. So, only one case of demand growth conforms to Stigler's hypothesis of positive

  so that the distance between two neighboring firms is constant for a given . Entry or exit can only be addressed by increasing or decreasing . All the firms have the information on the distribution of n n θ .

  intermediate equilibrium price exceeds the marginal cost of upstream firms ( final price p* exceeds the marginal cost of downstream firms ( + > * *

  the total production of the other downstream firms, then the downstream firm θ 's profit is : the number of downstream firms is simply θ . n and that each downstream firm θ will produce the same output θ θ . n y y d = . Considering the first order condition from (7), and substituting for p(x,y), equation (8) defines the inverse demand function of the intermediate good for each upstream firm :

Table I :

 I Increasing the number of firms (with β = 0.02, α =2, I =25)

	n	10	15	20	25	30	35
	θ	0.15 0.23 0.29 0.32 0.35 0.37
	θ	0.93 0.77 0.66 0.59 0.53 0.49
	Integrated Firms (share) 78% 54% 37% 27% 18% 12%
	Downstream F (share) 15% 23% 29% 32% 35% 37%
	Upstream F (share)	7% 23% 34% 41% 47% 51%

Table II :

 II Growth in the number of consumers(with β = 0.02, α =2 , n = 15)

	I	5	10	15	25	40	60
	θ	0.30 0.26 0.24 0.21 0.18 0.15
	θ	0.49 0.61 0.70 0.81 0.91 0.97
	Integrated firms (share) 18% 34% 46% 61% 72% 82%
	Downstream F (share) 30% 27% 24% 21% 19% 15%
	Upstream F (share)	52% 39% 30% 18% 9% 3%

Table III :

 III Proportional growth in the number of consumers and firms11

	n	6	12	18	30				42			60	78
	I	10	20	30	50				70 100 130
	θ	0.09 0.18 0.24 0.31 0.36 0.4 0.43
	θ	0.98 0.85 0.78 0.70 0.66 0.6 0.56
	Integrated firms (share) 90% 68% 54% 39% 30% 20% 13%
	Downstream F (share) 8% 17% 24% 31% 36% 40% 43%
	Upstream F (share)	2% 15% 22% 30% 34% 40% 44%
	Proposition 2 :
	(with β = 0.02, α =2, and n/I = 0.6)
				q	=	α	(	I	/	n	)	-	β	(	I	/	n	).	p

9

Because the average individual output for each firm can be written :

A proportional growth in the number of consumers and firms leads to an increasing proportion of firms choosing to be vertically integrated.

	5.4 GROWTH IN THE DEMAND PRICE
	Another way of considering demand growth is to assume an increase in the value
	granted to the final good by each consumer. A rise of α without modification of β
	(equation 4) results in a fall of the consumption's price elasticity and in a rise of
	demand for a given price. Table 4 reveals that it results in a growing number of
	integrated firms.

However, we can't assert that the number of consumers and the number of producers are positively correlated, particularly in a model where average costs are continuously decreasing. So we can't clearly conclude on the validation of Stigler's thesis at this stage.

Table IV :

 IV Growth in the value to consumers (with β =0.02, I = 25, n = 15).

		α 0.8	1	1.2 1.4 1.8 2.2 2.8
	θ	0.37 0.34 0.31 0.28 0.23 0.19 0.13
	θ	0.46 0.53 0.6 0.66 0.77 0.86 0.95
	Integrated firms (share) 8% 19% 29% 38% 54% 67% 82%
	Downstream F (share) 38% 34% 31% 28% 23% 19% 13%

Because their main purpose is to study foreclosure strategy at equilibrium and not really the vertical integration degree.

So the issue here is based on a different model that of[START_REF] Stigler | The division of labor is limited by the extent of the market[END_REF] or[START_REF] White | The division of labor is limited by the standardization in the market : a note on Stigler's vertical integration theme[END_REF] who predict that all firms will become more integrated or less integrated as the market increases.

[START_REF] Perry | Vertical integration and industry growth[END_REF] assume that the average cost curves are increasing for downstream production, decreasing upstream and U-shaped for the integrated firms. It is not clear why to a given stage corresponds a given cost curve. So we choose not to consider that the form of the upstream and downtream costs curves differ. They also assume that the difference in efficiency is reflected in fixed cost for downstream and marginal cost for upstream production, which is not well motivated either.

The variable cost of buying intermediate goods will be introduced into the profit function.

This leads to a decrease in the price elasticity of consumption which is independent of I.

The integrated firms are indeed in competition with the downstream firms, therefore even their number falls, the number of their competitors increases (cf table 1).

We checked that it is also accompanied by a drop in intermediate and final prices.

. Substituting for from (7) and maximizing, we have :

The boundaries θ * and * θ are defined by conditions ( 16) and ( 17) 6 .

correlation between industrial specialization and industry growth in our two stages model with Cournot competition.

Finally this model can also contribute to research on the effect of vertical integration on the final price (and on the consumer's surplus) in the context of successive oligopolies. In our model, an increase in the proportion of integrated firms induces ceteris paribus a decrease in the final price.