P Chinda

W Wechsatol

S Chanchaona

P Brault

Micro-Scale Modeling of an Anode-Supported Planar Solid Oxide Fuel Cell

Keywords: Solid Oxide Fuel Cells, Micro -Scale Model, Mass Transfer, Electrochemical Reaction, Micro Scale Geometry, Spherical -Shaped Particles, Parametric Study

A microscale model of a Solid Oxide Fuel Cell (SOFC) involving the mass transfer together with the electrochemical reaction, the transportation of electrons and ions through the respective spherical shaped electron conducting and ion conducting particles inside the electrodes was mathematically developed.

Couples of useful parameters were introduced in order to represent the characteristics of the cell. The predicted cell performance was showed according to various operating and design conditions. The effects of micro scale electrode geometry on the cell performance were also taken into account. Parametric study according to the volumetric fraction of ionic and electronic conducting particles was conducted in order to examine the effects of operating conditions on the cell overpotentials. The study results substantiate the fact that SOFC overpotential could be effectively decreased by increasing the operating temperature as well as operating pressure. This present study reveals the working mechanisms of SOFC at the micro scale level, while demonstrating the use of micro scale relations to enhance the SOFC performance. The accuracy of the presented model was validated by comparing with already existing experimental results from the available literatures.

Macroscopic Approach to Alternative Microscopic Modeling

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. The supplied fuels can be hydrocarbon compounds such as CH 4 , C 2 H 6 and C 3 H 8 , as well as CO beside H 2 . The SOFC can be applied to versatile power generation systems, either stationary power generators or auxiliary power sources in automobiles, in aircrafts or even in 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]. Many researchers have been focusing on developing the SOFC fabrication techniques, and new materials to make the SOFC possible commercial applications. Ramakrishna et al. [START_REF] Ramakrishna | [END_REF] performed numerical analyses of a SOFC on a CFD -ACE package and proposed an innovative thinwalled geometry to improve the power density of the cell.

Recently, a novel 3-D SOFC was developed by Koh et al. [4] to increase the volumetric power density by increasing the volumetric surface areas of the electrodes. Yang et al. [5] and Hwang et al. [6,7] developed a numerical model to investigate the characteristics of a monoblock layer built SOFC (MOLB SOFC), and found that the MOLB SOFC provides higher fuel/oxidant utilization than the planar type SOFC. The high utilization of the fuel and oxidant reflects the high current generated by the SOFC. However, the current density distributions are less uniform in MOLB SOFC, which is a significant disadvantage. The major limitation of SOFC usage is their high operating temperature. Intermediate temperature SOFCs and low temperature SOFCs become the subject of interest to lower the operating temperature as well as the material cost of SOFCs. The techniques of thin electrolyte and electrode coating involving new materials have made the creation of intermediate and low temperature SOFCs possible. Yttria Stabilized Zirconia (YSZ), which is normally used to produce cell electrolytes operating at high temperatures, can still be used in intermediate and low temperature SOFCs. Many studies on Ceria as a possible material for fabricating electrolyte are widely conducted by many researchers [8 -11]. The thin film deposition, a technique widely used for decades in fabrication of Proton Exchange Membrane Fuel Cells (PEMFCs) [12] is now applied in SOFC fabrication [13 -15]. One of the main advantages of such a thin film fuel cell is that its ohmic resistance is low enough to maintain its high performance at low operating temperature.

In an analysis of the SOFC system, computer simulation based on theoretical modeling is known to be an efficient method to predict the cell, stack or system performance. Mathematical models that predict performance can aid in understanding and development of SOFCs. A mathematical simulation of a SOFC is helpful in examining issues such as temperatures, materials, geometries, dimensions, fuels, and fuel reformation and in determining their associated performance characteristics. When physical properties or reaction kinetics are not known reliably, they can be estimated by fitting performance data on small sizes, laboratory scale cells to a mathematical model. Thus, modeling is an important tool in design optimization, helping to answer important practical questions such as what air and fuel flow rates must be used to avoid excessive temperature or pressure drop. On the other hand, by providing answers to questions such as how much the electrical properties of the cell materials must be improved, simulations at the cell and electrode level or the macroscopic level and microscopic level modeling can guide to development of new and improved SOFCs's materials. Mathematical simulation, therefore, has the potential to guide technology development, test the significance of various design features, asses the effectiveness of developments in materials of fabrication procedures, and select optimum operating conditions from a set of feasible parameters. Various modeling approaches exist. The modeling may focus on individual thermal-mechanical, flow, chemical, and electrochemical subsystems or on coupled integrated systems. Because the subsystems are typically characterized by different length scales, modeling may also take place on different levels, ranging from the system level modeling, the cell stack level modeling and the single cell level modeling.

The viewpoint of single cell level modeling is divided into two levels those are macroscopic modeling and microscopic modeling. The macroscopic are first examined from the viewpoint of fluid and mass transport whether heat transfer phenomena that can be obtained understanding of theses transport phenomena in the cell.

The effects of porosity on reactant concentration and ohmic overpotentials could also be considered in the macroscopic modeling. The macroscopic appearances of electrodes such as tortuosity and porosity are also related to the microscopic arrangement of electrode particles. However, the activation overpotential was being treated as a property independent of the porosity, while the electrochemical reactions was assumed to occur merely at the interface between electrode and electrolyte layers. Thus macroscopic modeling can hardly predict the effects of electrode structure on the chemical reaction within SOFCs so, it is insufficient to investigate precisely the micro-structural effects. Another level model is the microscopic modeling, which is considered inside the electrodes. Influence of the electrode structures on the electrochemical reaction at the three phase boundary (TPB) has recently been included in microscopic modeling. The microscopic modeling has been simulated and applied to SOFC in recent years and has offered possibilities to enhance the SOFC performance [16][17][18][19][20]. In literatures for microscopic modeling, the properties of porous electrodes and their effects on SOFC performance have been investigated from different aspects.

Costamagna et al. [16] developed a theoretical microscopic model for the evaluation of the performance of an electrode formed by a mixture of electronic conductor and ionic conductor with the spherical particles 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 EDB/Pt and YSZ/LSM cathodes and YSZ/Ni anodes. Chen et al. [17] developed a mechanistic model for oxygen reduction at 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. [18] developed a microscopic model for a composite cathode similar to Chen et al. [17]. 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. [19] considered a 1-D anode-supported planar SOFC with thin layer reaction zone in the vicinity of electrolyte. Their numerical results shown 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 reaction zone. Recently they extended their investigation to the cell design with two distinct layers of anode electrode, e.g., reaction layer and backing layer. Ni et al. [20] developed a mathematical model for modeling the performance of solid oxide fuel cell (SOFC) 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 micro-structurally graded SOFC. 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.

Unlike other previous researches, the effect of electrodes' microscopic geometry on the cell performance of an anodesupported Solid Oxide Fuel Cell is included in this microscopic scale model. For example, porosity and tortuosity are parameters that depend on the shape and size of electrode particles. In this study, the microscopic model of an SOFC porous electrode that formed by mixture of electronic conductor and ionic conductor was developed. The model is modeled with spherical particles electrode and takes into account electronic, ionic, and gas transport together with the electrochemical reaction. The forms of overpotentials losses are considered aims to optimize the SOFC performance by optimizing the microstructure relations of SOFC electrode by microscopic scale modeling base on the random packing sphere model [16].

A microscopic scale model of a Solid Oxide Fuel Cell (SOFC) involving electrochemical reactions and mass transport behaviors throughout electrodes packed with spherical-shaped particles was mathematically developed.

The model was validated by comparing the predicted results with already available experimental results. By minimizing the cell overpotentials, the model provides the optimal value of the electrode parameters. 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. The processes to modeling can be explained below:

The Cell Geometry and the Model Assumptions

The Cell Geometry

An illustration of physical domain of a SOFC is shown in figure 1. Part (a) in figure 1 shows a structure of a unit cell with multiple flow channels, whereas part (b) in figure 1 shows the cross-sectional representation of a composite anode and also part (c) in figure 1 shows the crosssectional representation of a composite cathode.

The computational domain for anode is shown by part (b) in figure 1, which includes the land portion of the interconnect (current collector) interfacing the anode, the porous portion of the anode interfacing the flow channel, the anode and another is interface to an electrolyte. The computational domain for cathode is shown by part and (c) in figure 1, which includes the land portion of the interconnect (current collector) interfacing the cathode, the porous portion of the cathode interfacing the flow channel, the cathode and another is interface to an electrolyte. The electrolyte was made of Yttria Stabilized Zirconia (YSZ).

From the part (b) in figure 1, the light gray spherical particles represent the ionic conducting particles that in this case are YSZ (Yttria -Stabilized Zirconia) and the dark gray spherical particles represent the electronic conducting particles that in this case is Nickel (Ni) contain inside the anode, respectively. The parameter N i represents the rate of mass transport (or denotes as the molar flux of each species, mol m -2 s -1 ) both of reactant (H 2 ) flow from the channel to the reaction sites of the anode and product (H 2 O) flow remove from the reaction sites to the channel. In this study the transport of N i vary only in the vertical direction or in the direction x which is shown in the figure 1 (b) not vary in the horizontal direction or the direction y in figure 1 (b). The parameter x a represents the anode distance which x a = 0 means the position at the interface between the anode and flow channel and at x a = t an means the position at the interface between anode and electrolyte. In the similar manner, from the part (c) in figure 1, the light gray spherical particles represent the ionic conducting particles that in this case are YSZ (Yttria -Stabilized Zirconia) and the dark gray spherical particles represent the electronic conducting particles that in this case is LSM (Strontium -doped Lanthanum Manganite) contain inside the cathode, respectively. The parameter N i represents the molar flux of oxidant or oxygen species (mol m -2 s -1 ) flow from the channel to the reaction sites of the cathode. In this study the transport of N i vary only in the vertical direction or in the direction x which is shown in the figure 1(c) not vary in the horizontal direction or the direction y in figure 1 (c) same as the anode side. The parameter x c represents the cathode distance which x c = 0 means the position at the interface between the cathode and flow channel and at x c = t cat means the position at the interface between cathode and electrolyte.

The Model Assumptions

This model is considered only one cell. The packed spherical shaped particles inside the electrodes create the void space, the tortuous flow paths and the interfacial surface between ionic and electronic conducting particles.

The cell is assumed to operate under steady state, and to assume that the parameters vary only in the vertical direction. Temperature and total pressure is assumed to be uniform throughout the entire electrode. The total pressure gradient inside the electrode is assumed to be negligible. The convective flux is negligible in the porous anode when compared to the diffusive flux of gaseous species. This implies that the primary mode of specie transport within the porous anode is by diffusion. The mass transfer process inside the porous electrode is by means of diffusion was also selected by the recommendation in the literatures [21 -23]. Suwanwarangkul et al [21] were 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 Dusty Gas Model (DGM) that is considered both Knudsen diffusion and Molecular diffusion. The fuel or the anodic reactant is the humidified with H 2 0.95 mole fraction and H 2 O 0.05 mole fraction. The oxidant or the cathodic reactant is air with O 2 and N 2 the mole fraction of 0.21 and 0.79, respectfully. All gases involved in electrochemical reaction behave like an ideal gas, the Soret, Dofour and gravity effects are negligible. The model is negligible the momentum effect in the electrode. The electrode is a composite porous structure, it is modeled based on the assumption that the composite electrode is represented by the particles of electron and ion conducting particles packed together in random and disregards the actual geometric details of the individual particles, so and according to Costamagna et al. [16] for modeling the electrode with random stacks of nearly spherical shaped particles. The each phase of electrode is considered to be homogeneous.

Cell Performance and Cell Geometry Relation

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 electron conductivity. Because of the complex nature, a detailed simulation technique was employed to predict cell performance under various operating conditions. The present model aims to study the effect of electrodes' microscopic geometry on the cell performance of an anode -supported Solid Oxide Fuel Cell and the parametric study according to the volumetric fraction of ionic and electronic conducting particles was conducted in order to examine the effects of operating conditions on the cell performance. For simulation with the microscopic scale modeling in this work the model considered on mass transport, electrode electrochemical reactions and charge transport. The details on the modeling are described below as follow by starting on considered with the anode side, cathode side and electrolyte, respectively.

Anode Side

The configuration of a SOFC anode can be simply modeled as shown by the part (b) in figure 1. The coupled electrochemical reaction and mass transfer in a composite anode are illustrated by figure 2. 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, [START_REF] Larminie | Fuel Cell Systems Explained[END_REF] 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. Therefore, the anode is usually formed by the mixture of ionic and electronic conducting materials with proper porosity present in the structure so that gases can reach the reaction sites. In the same way the configuration of the composite cathode of a SOFC consists of both electronic and ionic conducting particles, the mathematical model can be formed in a similar manner.

Charge balance in the electronic and ionic conducting particles can be written as e,a a

J = R or e,a v n,a J = A J (1) and i,a a J = -R or i,a v n,a J = -A J (2)
where J e,a and J i,a are the current density (Am -2 ) in the electronic and ionic conducting particles. R a is the volumetric current density produced in the anode due to hydrogen oxidation reaction (A m -3 ). A v is the electrochemically active surface area per unit volume of the porous anode (m 2 m -3 or m -1 ) and J n,a is the transferred current per unit area of reactive surface (A m -2 ). The term 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 [16], given as

22 el io v el t e l io el io ZZ A = sin r n n n p p Z ( 3
)
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 [16] 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 [16,17] are related to each other as shown here,

el 3 io el n = [ +((1-)/(r /r ) )]
and

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

2 io el io 2 el el io el (Z-3)(r /r ) Z = 3+ [ n ( 1 n )( r / r ) ] (5) 0.4 2.5 el el -el p = 1-4-Z and 0.4 2.5 io io -io p = 1-4-Z (6)
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 [16,17],

2 el el el el nZ Z = Z and 2 io io io io nZ Z = Z (7)
The total number of particles per unit volume (n t ), depending on porosity and particle radius sizes (radii r el and r io for electronic and ionic conducting particles, respectively) can be determined according to Eq. ( 8) [16],

t 33 el el el io el 1- n = (4/3) r [ n ( 1 n )( r / r ) ] (8)
According to Butler -Volmer equation [18], the volumetric current density produced in the anode, R a and the transfer current density J n,a can be determined from : 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 [24], z is the number of charges involved in a reaction.

H2 2 H H2 v 0,ref H2,ref zF ( 1-)zF c aa R = A J exp -exp - a c
R is universal gas constant equal to 8.314 Jmol -1 K -1 and T is cell operating temperature (K).The overpotential at the anode η a [16,24] is

00 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.

The electronic and ionic potentials can be obtained from applying Ohm's law as Eq. (12). where ρ eff [19] 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) [19].

and eff i,a i,a = ( 1 )(1-) ( 13
)
where σ e is the conductivity of electronic conductor -i.e. for an anode, e,a is 6 2 10 Sm -1 [16,24] and for a cathode, e,c is 4 1 10 Sm -1 [25]. σ i is the ionic conductivity -i.e. for an anode i,a is 4 3.44 10 exp(-10300 T ) , for a cathode, i,c is also τ is found by applying a geometrical model for tortuosity of streamline in porous media when the particles have a spherical shape [26]. The relevant equations for the tortuosity calculation are shown in an appendix A.

From Eq. ( 11), Eq. ( 12) and Eq. ( 13), the first and second derivatives of η a are, The term A v in Eq. ( 16) relates to the electrode porosity ( ), the particle radius (r io , r el ) and a volume fraction of the electron conducing particle in an electrode ( ). 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) to determine the overpotential at the anode η a , can be written as shown in Eq. (17),

eff a a e,a total x 0 d x = x = 0 , = - J dx and eff a a an total i,a x t an d x = x = t , = J dx ( 17
)
The positions of x = x a = 0 and x = x a = t an are shown by part (b) in figure 1.

While the electrochemical reactions due to mass transport occur in SOFC electrode. Mass transport through porous medium can be determined using concepts described in [27]. For a single component, the mass transport equation can be written as

-3 -1 i ii ( 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 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 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 right-hand side is significant, Eq. ( 18), therefore, becomes

ii 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 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 use 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 medium [21].

Normally, to describe the transport of components within porous media the simplest 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 [21].

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, it is assumed in this study 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. So, which means the primary mode of species transport in the porous electrode is by diffusion is considered in this work.

The mass transfer process inside the porous electrode is by means of diffusion was also selected by the reason that describe above and by the recommendation in the literatures [21][22][23]. Suwanwarangkul et al [21] were concluded that the Dusty-

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 [21]. 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 i i 0 i eff eff eff j=1,j i ij i,k i,k y N y N N dy B P 1 dP + = - P y 1 RT dx dx D D D ( 22
)
where x is a vertical direction that shown in figure 1 (b) and (c), y i and y j is the mole fraction of species i and j.

From the DGM that is considered both Knudsen diffusion and Molecular diffusion, so diffusion process within a pore can be typically divided into two diffusion mechanisms: Molecular diffusion 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 [START_REF] Bird | Transport Phenomena[END_REF]:

1 eff i eff eff im i,k 11 D = DD ( 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 theses parameters are described in [22,[START_REF] Chan | [END_REF]. It is 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 of hydrogen and water vapor or H 2 (1) -H 2 O(2) in the anode, the effective diffusion for binary systems ( eff i,B D ) becomes :

1 1 eff i,B eff eff eff eff ij 12 i,k 1,k 1 1 1 1 D = = D D D D ( 24
)
The Knudsen diffusion coefficient ( i,k D ) and effective Knudsen diffusion coefficient ( eff i,k D ) for each gas specie [START_REF] Bird | Transport Phenomena[END_REF] can be calculated by Eqs. ( 25) and (26).

i,k ii 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 ij D T 1 1 D = 0.0018583 MM p (27)
where

ij i j = ( + )/
) ( 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 listed in the table 1 [30].

Similar to the effective Knudsen diffusion coefficient, the effective binary diffusion coefficient through the porous electrode can be calculated from Eq. [START_REF] Chan | [END_REF].

eff ij ij D = D ( 29
)
Consider again on Eq. ( 22) that shown 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 Eq. ( 22) reduces to:

n j i -i j ii eff eff j=1,j i ij i,k y N y N N dy P + = - RT dx DD (30)
For diffusion with heterogeneous chemical reaction, the flux ratios are governed by reaction stoichiometry. Summing Eq. ( 30) over the n species leads to the Graham's law of diffusion in gaseous mixtures [27].

n ii i = 1 N M = 0 (31)
where M i is the molecular weight of component i. For two binary component systems of hydrogen and water vapor or (H 2 (1)-H 2 O(2) in the anode, so Eq. ( 30) becomes : N N = -M /M (Graham's law), N 1 can be written as follow : For two binary component systems of hydrogen and water vapor or (H 2 (1)-H 2 O(2) in the anode, so Eq. ( 34 (39)

1 2 1 -
In the above mentioned equations, the gas molar fraction, gas diffusion coefficients, current density, reaction rate, and the overpotential are dependent of depth in the vertical direction (x). Due to the complicated phenomena involved, an iterative scheme is used to solve the non-linear equations. The iterative computation is successful if a convergence is obtained. Therefore, the differential equation was solved numerically in MATLAB [START_REF] Gilat | Numerical Methods for Engineers and Scientists: an Introduction with Applications Using MATLAB 1 st Edition[END_REF]. After obtaining the distributions of species concentration, current densities, and overpotentials, the overall electrode overpotential can be determined [16].

an 00 t ,a e,a i,a e,a(x= 0) i,a(x = t ) ( - ) -( - ) (40)

Cathode Side

The configuration of a SOFC cathode can be simply modeled similar to anode that shown by the part (c) in figure 1 by instead the electronic conducting particle with LSM (Strontium -doped Lanthanum Manganite) and the ionic conducting particle with YSZ (Yttria -Stabilized Zirconia). The reaction processes can be summarized as at the same time of the reaction at the anode side. The oxidant (oxygen) is introduced at the cathode side where it diffuses through the porous cathode structure and reacts with electrons from the external circuit to produce O 2-at the vicinity of active sites. Then, they are transferred to the anode side through the electrolyte.

Similar to the anode, the coupled electrochemical reactions and transport behaviors at the cathode can be determined from 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. [16] that shown in above section similar to the anode part.

Similar to the anode, the boundary conditions for Eqs. ( 42) and ( 43 Because the flux of nitrogen is not participate in the electrochemical reaction on the cathode side, that only the oxygen to participate in the electrochemical at the cathode. The concentration of O 2 can be determined by self-diffusion mechanism [START_REF] Chan | [END_REF], and can be expressed below.

About Cathode self-diffusion on the cathode side,

2 2 2 2 2 O c O O O O J = -D C + y J (45) 2 2 2 2 2 O ,k(eff) O O ,k(eff) O -N (eff) D = D + D (46)
where J O2 is oxygen mol flux (mol s -1 ), C O2 is O 2 concentration (mol m -3 ) and

The Electrolyte

The dense electrolyte is used to conduct O 2-from the cathode side to electrochemically react with reactant (hydrogen) at the vicinity of the anode active sites to produce H 2 O and electrons. These electrons flow to the cathode side via an external circuit traveling through an electrical load or resistor. It typically made from Yttria -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.

Model Validation

For validating the present model we validated our results with available experimental data [33]. By using the parameters listed in table 2. The predicted cell performance is compared with the experimental data found in the literature [33], and is shown in figure 3. The performance of the cell is predicted when the cell is supplied with 95% H 2 and 5% H 2 O as fuel, operating at a temperature and pressure of 800 C and 1 atm, respectively. Oxygen composition in the ambient air is used as oxidant. The anode, the cathode and the electrolyte were made of Nickel -doped Yttria -Stabilized Zirconia (Ni/YSZ), Strontium -doped Lanthanum Manganite (LSM or La 0.8 Sr 0.2 MnO 3 ) and Yttria -Stabilized Zirconia (YSZ), with the thickness about 1 mm., 50 m. and 10 m., respectively.

It can be seen that the present model predictions agree well with the experimental results. It is worth mentioning that all the parameters used in the model validation are obtained from Rogers et al. [33] except the parameters such as particle sizes, volumetric fraction of electronic and ionic conducting particle given with the other reference that from Hussain et al. [19] due to these parameters are not showed in Roger et al's work.

Moreover, the value of these parameters (particle sizes, volumetric fraction of electronic and ionic conducting particle etc.) are typically use in many literatures [16 -20, 24 -25, 29, 32]. The value of tortuosity is varied to obtain the best agreement between the present model predictions and the experimental results shown in figure 3, since the value reported by Rogers et al. [33] is unity. The typical tortuosity value for SOFC electrodes is in the range of 2-6 [34]. Hence, the tortuosity value of 2.75 used in the present model predictions is in the typical range for SOFC electrodes and provides the best agreement with the experimental results. The R 2 -value, an indicator between 0 and 1, reveals how closely the predicted values correspond to the experimental data, is obtained for the present model as 0.99985. This model was employed to gain understanding of the coupled mass transfer and electrochemical phenomena in the electrodes of the SOFC. The present study reveals the working mechanisms of SOFC at the micro scale level, while demonstrating the use of microstructure effect to enhance the SOFC performance and the micro characteristics of the electrodes are included into the model. The predicted cell performances are showed as the cell overpotentials at various operating and design conditions. Moreover, parametric analyses were performed to investigate the effects of operating parameters on SOFC performance. 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 as in table 3 used in the simulation and the simulation results are shown below.

The Effect of Electrode Microstructure on the SOFC Overpotentials

Because one of the performance of the SOFC was known to strongly depend on the electrode properties i.e.

the volume fraction of the electron conducing particle in an electrode (φ), the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ) [35] and the electrode porosity ( ). All referring parameters are related to the electrode performance for electrochemical reaction. So, in this study the microstructure parameters relation are established and used to optimize the electrode performance by the effect of electrode microstructure on the SOFC overpotentials have been investigated. Firstly, this work is considered the effect of the ratio of particle size between the ionic and electronic conducting particle of the electrode (r io /r el ) and the effect of the porosity ( ) on the electrode performance.

Figures 4a and4b show the effect of the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) on the anode overpotential ( a ) and the effect of the porosity ( ) on the anode overpotential ( a ), respectively. For figure 4a, in the case of to vary the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ), the operating pressure and temperature are 1 atm and 1073 K. The typical operating current density of 10,000 Am -2 , the electrode porosity ( ) or the void fraction of the electrode is equal to 0.4, the volume fraction of electron conducting particles (υ) of the solid part of the electrodes is equal to 0.6.

The supplied fuel compositions are H 2 and H 2 O by 0.95 and 0.05 mole fractions respectively and oxidize by oxygen in the air by 0.21 mole fraction and 0.79 mole fraction of N 2 . The others standard parameters used in the calculation are given in table 3. For figure 4b in the case of to very porosity, all parameters for calculation are the same as figure 4a except that the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ) is equal to the unity. Figure 4a, the x -axis represents the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) at the anode side whereas y -axis represents anode overpotentials ( a ) and figure 4b, the x -axis represents the anode porosity ( ) whereas y -axis represents anode overpotentials ( a ).

Figures 5a and5b show the effect of the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) on the cathode overpotential ( c ) and the effect of the porosity ( ) on the cathode overpotential ( c ), respectively.

The following 5a and 5b, all the conditions and parameters for calculation same as figures 4a and 4b. Also in figure 5a, x -axis represents the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) at the cathode side whereas y -axis represents cathode overpotentials ( c ) and figure 5b, the xaxis represents the cathode porosity ( ) whereas y -axis represents cathode overpotentials ( c ). For both anode and cathode that show in figures 4a and 5a, the minimum overpotentials were found at the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) of about 1 to 1.25. The results agree well with the results from Costamagna et al.'s work [16]. As shown in Costamagna et al's work [16], the ratio between the particle radii is near to unity, the high performance electrode occurs at the ratio of particle size between the electronic and ionic conducting particle (r el /r io ) between 0.8 to 1 or it means that ratio of particle size between the ionic and electronic conducting particle (r io /r el ) between 1 to 1.25 occur at the volumetric fraction of the electron conducting particle (υ) between 0.2 to 0.65. It should be underlined that, in general, the particles do not have an exactly spherical shape, nor is the radius exactly defined, on the contrary there is a distribution around a mean value. This can influence the percolation thresholds in a real electrode.

Similarly, the effect of electrode porosity that show in figures 4b and 5b, the minimum overpotential is found at the porosity between 0.4 to 0.5. Reducing porosity that means has more solid particles are available to facilitate the transport of electronic and ionic charges which can increase the cell performance. In addition with decrease porosity can increase in active surface area for electrochemical reaction. However with decrease porosity is less space available for gas transport. The combined effects for minimization of the overall result in an optimal porosity for the electrode. Because the anode in an anode supported SOFC is thicker than the cathode so its overpotential must be higher than the cathode and the anode side has more effect on cell performance than the cathode.

The effect of volumetric fraction of the electron conducting particle (υ) of the solid part of the electrodes on the overpotentials ( a ) at various ratio of particle size between the ionic and electronic conducting particle (r io /r el ) at the anode side is shown in figure 6a. The effect of υ was studied at the typical operating current density of 10,000 Am -2 .The operating pressure and temperature are 1 atm and 1073 K, the electrode porosity ( ) or the void fraction of the electrode is equal to 0.4. The supplied fuel compositions are H 2 and H 2 O by 0.95 and 0.05 mole fractions respectively and oxidize by oxygen in the air by 0.21 mole fraction and 0.79 mole fraction of N 2 . The others standard parameters used in the calculation are given in table 3. The effect of volumetric fraction of the electron conducting particle (υ) of the solid part of the electrodes on the overpotentials ( a ) at various porosity ( ) on the anode side is shown in figure 6b. For figure 6b in the case of to very porosity, all parameters for calculation are the same as figure 6a except that the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ) is equal to the unity. Figure 6a, the xaxis represents the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) at the anode side, whereas the yaxis represents the anode overpotentials ( a ). Also in figure 6b, the xaxis represents anode porosity ( ) whereas yaxis represents anode overpotentials ( a ). Both results indicate that the minimum anode overpotentials occurred at υ which are equal to 0.5, are found at the ratio of particle size between the ionic and electronic conducting particle (r io /r el ) of about 1 to 1.5 and the at the porosity ( ) between 0.35 to 0.45.

The changes in a with υ demonstrate 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 [16]. When equals to 0.5, the volume fraction of ion conducting particles (1 -) of the solid part of the electrodes becomes equal to the volume fraction of electron conducting particles (υ) of the solid part of the electrodes 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 layers are functions of . Increasing beyond 0.5 increases the effective electronic conductivity but decrease the effective ionic conductivity, which results in increasing the ohmic overpotential.

The Effect of Operating Parameters on the SOFC Overpotentials

Figures 7 and8 show the effect of operating temperature and the effect of operating pressure on the overpotentials ( a ) at various ratio of particle size between the ionic and electronic conducting particle (r io /r el ) and at various porosity ( ), respectively. Both results are considered at the anode side. For figure 7a, in the case of to vary the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ), the typical operating current density of 10,000 Am -2 , the electrode porosity ( ) or the void fraction of the electrode is equal to 0.4, the volume fraction of electron conducting particles (υ) of the solid part of the electrodes is equal to 0.5.

The supplied fuel compositions are H 2 and H 2 O by 0.95 and 0.05 mole fractions respectively and oxidize by oxygen in the air by 0.21 mole fraction and 0.79 mole fraction of N 2 . The operating pressure is at 1 atm and the operating temperature is varied. The others standard parameters used in the calculation are given in table 3. For figure 8a, all of the calculation parameters are the same as in the case of figure 7a except the operating temperature is at 1073 K and the operating pressure is varied. For figure 7b in the case of to very porosity, all parameters for calculation are the same as figure 7a except that the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ) is equal to the unity. In the similar cased also for figure 8b in the case of to very porosity, all parameters for calculation are the same as figure 8a except that the particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ) is equal to the unity.

Figures 7a and 8a, the xaxis represents the ratio of particle size between the ionic and electronic conducting particle (r io /r el ), whereas the yaxis represents the overpotential ( a ). Also, figures 7b and 8b, the x axis represents the porosity ( ), whereas the yaxis represents the overpotential ( a ). It can be seen from figure 7 that increasing the temperature can result in increasing in the reference exchange current density and the ion conductivity of electrolyte, which in turn increase the electrochemical reaction rate. In other words, cell performance increases with the increases of temperature. Increasing the temperature not only enhances the rate of electrochemical reaction at the reaction sites but also increases the electronic and ionic conductivities of the electron and ion conducting particles, respectively. Therefore, the overpotentials of electrodes and ohmic overpotentials of electrolyte decrease with increasing temperature which in turn minimizes the cell potential loss or the overpotentials and thereby enhances the cell performance. These phenomena indicate that 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. So, it is important to develop new materials with high ionic conductivity and high catalytic activity at intermediate temperatures and for higher performance can use the optimum relationships of electrode microstructure parameters to enhance cell performance so cells can perform at higher temperatures, or developed new materials for utilizing cells at lower temperatures.

Figure 8 shows that at the operating pressure equal to 1 atm, 3 atm and 5 atm. It can be seen that at high pressure, the electrode overpotential is lower than that 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 creates other difficulties, such as limitation on material selection, gas sealing, and mechanical strength of the cell components [START_REF] Minh | Science and Technology of Ceramic Fuel Cells[END_REF].

Gas Composition Distribution in Electrodes

The distribution of gas composition along the electrode depth was studied at the typical operating current density of 10,000 Am -2 . The electrode porosity ( ) or the void fraction of the electrode is equal to 0.4, the volume fraction of electron conducting particles (υ) of the solid part of the electrodes is equal to 0.5. The particle size ratio of the ionic to electronic conducting particle of the electrode (r io /r el ) is equal to the unity. Figure 9 shows the important characteristic on gas distribution in an anode -supported SOFC.

The model analysis shows the cell performance strongly depend on the electrode microstructure, i.e. the volumetric composition of the electrode, the dimensions of the particle both of ionic and electronic conducting particles and the void volume of the electrode. The combined microstructure effects in a monotonic can decrease in electrode overpotential and hence the better cell performance. 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 microscopic model SOFC is feasible and has great potential to improve the SOFC performance by simply modifying the electrode microstructures.

Conclusion

A micro scale model of a Solid Oxide Fuel Cell (SOFC) involving the mass transfer in the electrode together with electrochemical reaction, the transport of electrons and ions in the respective electronconducting and ion conducting particles in the electrodes packed with spherical shaped particles was mathematically developed.

A theoretical model has been presented for the evaluation of the performance of an electrode formed by a mixture of electronic conducting and ionic conducting particles. The effects of micro scale electrode geometry on the cell performance were also taken into account. Parametric study according to the volumetric fraction of ionic and electronic conducting particles was conducted in order to examine the effects of operating conditions on the cell overpotentials.

The model analysis shows the cell performance strongly depend on the electrode microstructure, i.e. the volumetric composition of the electrode, the dimensions of the particle both of ionic and electronic conducting particles and the void volume of the electrode. The combined microstructure effects in a monotonic can decrease in electrode overpotential and hence the better cell performance. 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. 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. This paper reveals that micro-structural electrode is very promising to enhance the SOFC performance and worthy of further investigation. Aside from their electric performance, other properties associated with micro-structural should be further investigated. The model predicted results are validated with the experimental data found in the literature. A good agreement is obtained between the predicted values and the measured data.

Appendices

Appendix A: Tortuosity Model

The geometrical model for tortuosity of streamline in three dimensional porous media when the particles have the shape as spherical which used in this SOFC modeling obtained by Meijuan Yun et al. [26]. The using equations are shown below. where the parameter ε is the porosity for the spherical particle packed electrode and d is the gap size between spherical particles. Charge transfer coefficient, [24] 0.5

List of Symbols

Conductivity of electronic conductor of anode, e,a (Sm -1 ) [16] 2e6

Conductivity of electronic conductor of cathode, e,c (Sm -1 ) [16] 1e4

Conductivity of ionic conductor of anode, i,a (Sm -1 ) [16] 3.44e4*(exp(-10300/T))

Conductivity of ionic conductor of cathode, i,c (Sm -1 ) [16] 3.44e4*(exp(-10300/T))

Faraday constant, F (C mol -1 ) [16] 96485 Fig. 3 The comparison between model predictions and experimental results of Rogers et al [33].

Fig. 4 (a) The ratio of particle size between the ionic and electronic conducting particle effect on anode overpotentials. (b) Porosity effect on anode overpotentials.

Fig. 5 (a) The ratio of particle size between the ionic and electronic conducting particle effect on cathode overpotentials. (b) Porosity effect on cathode overpotentials.

Fig. 6 (a) The effect of volumetric fraction of the electron conducting particle on the overpotentials at various ratio of particle size between the ionic and electronic conducting particle on the anode side. (b) The effect of volumetric fraction of the electron conducting particle on the overpotentials at various porosity on the anode side.

4 3 .

 3 44 10 exp(-10300 T ) and τ is the electrode tortuosity.

2 M

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

  The positions of x = x c = 0 and x = x c = t cat are shown by part (c) in figure 1.

  The supplied fuel compositions are H 2 and H 2 O by 0.95 and 0.05 mole fractions respectively and oxidize by oxygen in the air by 0.21 mole fraction and 0.79 mole fraction of N 2 . The operating pressure is at 1 atm and the operating temperature is 1073 K. The others standard parameters used in the are given in table 3. The part (a) of figure 9 shows distribution of H 2 and H 2 O in the anode region and also part (b) shows distribution of O 2 molar fraction distribution in cathode region. Figure 9, the xaxis represents the anode distance whereas yaxis represents H 2 and H 2 O molar fraction for part (a) of figure 9 and O 2 molar fraction for part (b) of figure 9. The part (a) of figure 9, the concentrations of H 2 decreases considerably from anode surface towards the electrolyte interface and, the concentration of H 2 O decreases from the electrolyte interface towards the anode surface. The part (b) of figure 9 the concentrations of O 2 small decreases from cathode surface towards the electrolyte interface. The part (a) of figure 9, they are indicated the importance of gas transport in the thick anode. In addition the thin cathode shows a very low mass transfer resistance as can be seen in part (b).

AvJJti

  Active surface area per unit volume [m 2 m -3 or m -1 ] Transfer current per unit area of reactive surface [A m -2 ] Reference exchange current density for H 2 oxidation [A m -2 ] Reference exchange current density for O 2 reduction [A m -2 ] M Molecular weight [kg kmol -1 ] n el Number fraction of electronconducting particle in the electrode n io Number fraction of ionconducting particle in the electrode n Total number of particles in the electrode per unit volume [m -3 ] N Molar flux of species i [mol m -2 s -1 ] p el Probability of percolation of electronconducting particles in the electrode p io Probability of percolation of ionconducting particles in the electrode P Pressure [atm]

ai 2 Reaction order for hydrogen oxidation reaction O 2 Reactionc

 22 conducting particles in the electrode [m] r io Radius of ionconducting particles in the electrode [m] R Universal gas constant, 8.314 [J mol -1 K -1 ] R Volumetric current density produced due to hydrogen oxidation [A m -3 ] R Reaction rate [mol m -3 s -1 ] r Mean value of the pore radii [m] Molar fraction of species i [mol mol -1 ] Z Total average particle coordination number in electrode which equal to 6 Z el Coordination number of electron -conducting particles in the electrode Z io Coordination number of ion -conducting particles in the electrode G Net standard Gibbs free energy of electrochemical reaction at 1 atm and 25 ºC Charge transfer coefficient H Cathode activation overpotential [volt] Contact angle between electron and ionconducting particles in the electrode Electronic conductivity [S m -1 ] Tortuosity e Electronic potential [volt] i Ionic potential [volt] Volume fraction of electron conducting particles in the electrode Tables

Fig. 1 (

 1 Fig. 1 (a) A structure of a unit cell with multiple channels. (b) A cross sectional representation of a composite anode computational domain. (c) A cross sectional representation of a composite cathode computational domain.

Fig. 2

 2 Fig. 2 Schematic of a SOFC anode and electrochemical reaction at three phase boundaries (TPB).

Fig. 7

 7 Fig.7(a) The effect of operating temperature on SOFC overpotentials at various ratio of particle size between the ionic and electronic conducting particle on the anode side. (b) The effect of operating temperature on SOFC overpotentials at various porosity on the anode side.

Fig. 8

 8 Fig.8(a) The effect of operating pressure on SOFC overpotentials at various ratio of particle size between the ionic and electronic conducting particle on the anode side. (b) The effect of operating pressure on SOFC overpotentials at various porosity on the anode side.

Fig. 9 (

 9 Fig. 9 (a) H 2 and H 2 O composition distribution along anode depth (b) O 2 molar fraction distribution in cathode

  The DGM is the most suitable model for the H 2 -H 2 O system[21]. Suwanwarangkul et al is 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. Then the Dusty -Gas model (DGM) was used to developed the expressions for N i in this work. The electrochemical reactions due to mass transport occur in SOFC electrode that represents by the term R i or R i is 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 :

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 Dusty Gas Model (DGM) that is considered both Knudsen diffusion and Molecular diffusion.

  2 is the collision diameter (in °A ), p is total pressure (in atm), M is the collision integral based on the Lennard -Jones potential[START_REF] Bird | Transport Phenomena[END_REF], which can be obtained from Eq. (28).

				i and M j are molecular
	weight of species i and j,
	1.06036	0.19300	1.03587	1.76474
	=
	Dij
	0.15610
		kT	kT	kT
	kT
	exp( 0.47635	) exp( 1.52996	) exp( 3.89411
		ij	ij	ij
	ij

ij D

  Equation[START_REF] Minh | Science and Technology of Ceramic Fuel Cells[END_REF] is in the form of an ordinary differential equation, which can be solved by involving the following initial conditions. In Eqs. (37) and (38), Dalton's law of an ideal gas,

																								2 y = P H H	2	H 2 0 P and the initial conditions
	are following as:
	IC.1 :	y				y																	(37)
		H					H ,bulk
			22
				x 0
		dy								1	y
	IC.2 :	H	22 H ,bulk total -RTJ =			1						(38)
										eff							eff
		dx				2FP
										DD
										H O H
				x = 0						22				H
																		2,k
	It can be obtained the mole fraction of H 2 O as :
																								)
	becomes
										1
				1-y			dy
	N = -	22 HH P1 +										(35)
	H
	2			eff					eff
		RT								dx
				DD
				H o H		H ,k
				2		2		2
	where	H = 1 -M /M	H O	. The parameters	H M , and	M	HO	are the molecular weight of H 2 and H 2 O,
							22									2					2
	respectively. The	D	eff	and	D	eff		are the effective Knudsen diffusion coefficient of H 2 and the effective
						H ,k		H o H
							2			22
	binary diffusion coefficient of H 2 -H 2 O [29] respectively.
	By substituting Eq. (35) into Eq. (19), the governing equation of H 2 concentration becomes
	H 2 d y 2	2		eff			eff 1	+	eff 1 -y	H	2	1	dy	H	2	2	-	n,a A J RT v	eff 1	1	eff	y	H	2	= 0	(36)
																dx					2FP
	dx		D			D			D														D	D
				H O H 2	2	H ,k 2		H O H 2	2										H	2,k	H O H 2	2

  Sr 0.2 MnO 3 ) and Yttria Stabilized Zirconia (YSZ). The dimension of cell components and standard parameters used in calculations are given in table3. The electrode porosity and tortuosity are estimated by the equations in an appendix A[26].

	4 Results and Discussion
	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 and oxidize by oxygen in the air. The anode, the cathode and the electrolyte were made of Nickel
	doped Yttria	Stabilized Zirconia (Ni/YSZ), Strontium	doped Lanthanum Manganite (LSM or
	La 0.8

Table 1

 1 List of Lennard-Jones potentials used in calculation (taken from Tod et al.,[30])

	Gas species	CH 4	H 2 O	CO 2	CO	H 2	N 2	O 2
	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

Table 2

 2 Parameters used for model validation (source: Rogers et al.[33])Conductivity of electronic conductor of anode, e,a (Sm -1 )71428.57Conductivity of electronic conductor of cathode, e,c (Sm -1 ) 5376.34Conductivity of ionic conductor of anode, i,a (Sm -1 ) 0.64Conductivity of ionic conductor of cathode, i,c (Sm -1 )

	Operating temperature, T ( C)

Table 3

 3 Standard parameters used in the SOFC model calculation

	Parameters	Values
	Operating temperature, T ( C)	800
	Operating pressure, P (atm)	1
	H 2 and H 2 O in fuel mixture inlet mole fraction	0.95, 0.05
	O 2 and N 2 in air inlet mole fraction	0.21, 0.79
	Anode thickness, t a ( m.)	2000
	Anode tortuosity, a	1.38
		(Eq. in an appendix A)
	of anode, ε a	0. 4
		(Eq. in an appendix A)
	Anode pore radius (r p ) (μm.)	0.5
	Cathode thickness, t c ( m.)	60
	Cathode tortuosity, c	1.38


Fatima Benbouzid-Sitayeb

Ismaïl Ammi
email: ammi_smail@yahoo.fr

Christophe Varnier
email: .cvarnier@ens2m.fr

Noureddine Zerhouni
email: zerhouni@ens2m.fr

An Integrated ACO Approach for the Joint Production and Preventive Maintenance Scheduling Problem in the Flowshop Sequencing Problem

Keywords: ACO, Production, Preventive maintenance, Integrated approach, Joint scheduling, Flowshop

In this paper, an integrated ACO approach to solve joint production and preventive maintenance scheduling problem in permutation flowshops is considered. A newly developed antcolony algorithm is proposed and analyzed for solving this problem, based on a common representation of production and maintenance data, to obtain a joint schedule that is, subsequently, improved by a new local search procedure. The goal is to optimize a common objective function which takes into account both maintenance and production criteria. We compare the results obtained with our algorithm to those of an integrated genetic algorithm developed in previous works. The results and experiments carried out indicate that the proposed ant-colony algorithm provide very effective solutions for this problem.

I. INTRODUCTION

Production scheduling and preventive maintenance planning decisions are inter-dependent but most often made independently. Both activities conflict since preventive maintenance activities consume potential production time, but delaying preventive maintenance because of production demands may increase the probability of machine failure, this interdependency seems to be overlooked in the literature.

One finds in the literature two joint scheduling strategies which aim is to solve conflicts between production scheduling and preventive maintenance planning [START_REF] Lee | Scheduling jobs and maintenance activities on parallel machines[END_REF]. The sequential strategy consists of two steps: First the scheduling of the production jobs then the insertion of the maintenance tasks, taking the production schedule as a strong constraint. The integrated one consists of simultaneously scheduling both maintenance and production activities based on a common representation of these two activities.

As it has been mentioned, the interaction between production and preventive maintenance, particularly their joint scheduling, is relatively little studied and rather recently in the literature. Little work has been carried out in which preventive maintenance scheduling and flowshop scheduling are jointly considered [START_REF] Aggoune | Minimizing the makespan for the flowshop scheduling problem with availability constraints[END_REF][START_REF] Cassady | Minimizing job tardiness using integrated preventive maintenance planning and production scheduling[END_REF][START_REF] Benbouzid | Resolution of joint maintenance/production scheduling by sequential and integrated strategies[END_REF][START_REF] Ruiz | Considering scheduling and preventive maintenance in the flowshop sequencing problem[END_REF].

Since the last decade, attempts are being made to solve combinatorial optimization problems using Ant Colony Optimization (ACO) algorithms. ACO is a metaheuristic allowing finding the best solution of hard-to-solve optimization problems [START_REF] Dorigo | Optimization, learning and natural algorithms[END_REF][START_REF] Dorigo | The ant system: Optimization by a colony of cooperating agents[END_REF]. The aim of ACO algorithm is to find a path (a sequence) of minimal cost while respecting the constraints [START_REF] Dorigo | Ant Colony Optimisation: A new meta-heuristic[END_REF].

ACO has been applied successfully to solve scheduling problems. The problem dealt with is flowshop scheduling problem with the objective of minimizing the makespan. Stüetzle [START_REF] Stuetzle | An ant approach for the flowshop problem[END_REF] compared ACO algorithm to some basic heuristics for the Flow Shop Problem (FSP) showing that this approach is very promising for the FSP. Rajendran and Ziegler [START_REF] Rajendran | Ant-colony algorithms for permutation flowshop scheduling to minimize makespan/total flowtime of jobs[END_REF] propose two ant-colony optimization algorithms and compared them to heuristic solutions given by Taillard [START_REF] Taillard | Benchmarks for basic scheduling problems[END_REF]. The comparison shows that the two proposed ant-colony algorithms perform better. Ying and Liao [START_REF] Ying | An ant colony colony system for permutation flowshop sequencing[END_REF] applied the ACO algorithm to permutation FSP and compared its performance with other meta-heuristics such as genetic algorithm, simulated annealing, and neighbourhood search. Their computational results showed that ACO was a more effective metaheuristic for the n/m/P/C max permutation flowshop scheduling problem. Hu and al. [START_REF] Hu | flowshop rescheduling problem under ruch orders[END_REF] proposed an ACO algorithm to solve flowshop rescheduling problem. Mutation operations are employed to enhance the ACO performance avoiding the algorithm convergence to local optimum. Numerical experiments showed that the proposed algorithm has high computation efficiency.

The term "maintenance" is nonexistent in the previously cited literature. As in industry, production scheduling and preventive maintenance planning are typically treated independently in the production systems and also in operations research literature.

Purpose of this work is to present an integrated ACO approach to solve the joint production and preventive maintenance scheduling problem using ACO algorithm. The proposed integrated ant-colony algorithm, called INTACO, is based on a common representation of production and maintenance data. The aim is to optimize a common objective function which takes into account both maintenance and production criteria. The outline of this paper is as follows: section 2 presents the context of the study: production and maintenance data and the objective functions to optimize, section 3 introduces the proposed algorithm, discusses some details of the new local search procedure used and presents the most important parameter settings. Afterwards, in section 4, convergence properties of the algorithm are investigated so as to obtain a tuned set of parameters. Then we present the obtained simulation results with some benchmark problems which we compare to another metaheuristic developed in previous work. Finally, in section 5, some conclusions are made.

II. MODELS AND PROBLEM STATEMENT

In the following, we will first present the production and the maintenance data. Then the common objective function to optimize One of the most frequent production scheduling problems is the FSP. This problem can be stated as follows: each of n jobs 1,…, n has to be processed on m machines 1,…, m in that order. The processing time of job i on machine j is p ij . The processing times are fixed, nonnegative, and may be 0 if a job is not processed on some machine. Further assumptions are that each job can be processed on only one machine at a time, the operations are not preemptable, the jobs are available for processing at time zero and setup times are sequence independent. Here we consider the Permutation Flow Shop Problem (PFSP) i.e., the same job order is chosen on every machine. The objective then is to find a sequence, i.e., a permutation of the numbers 1,…,n that minimizes the completion time C max (makespan) of the last job. Following the Pinedo's notation [START_REF] Pinedo | Scheduling: theory, algorithms and systems[END_REF] this problem can be denoted by F/prmu/C max . The PFSP is an NP-Hard problem [START_REF] Kan | Machine scheduling problems: classification, complexity and computations[END_REF].

The preventive maintenance [START_REF] Biroloni | Reliability engineering, theory and practice[END_REF] can be stated as follows: each machine is maintained at fixed time intervals previously established. The preventive maintenance tasks are periodic interventions occurring every T periods and each occurrence depends on the one preceding it on the same machine. A maintenance task consists of elementary operations which processing time p' is evaluated with more or less certainty. Moreover, the periodicity T of these tasks can vary in a tolerance interval noted [T min , T max ]. This interval gives some flexibility to plan maintenance tasks while respecting the production constraint, disturbing the least possible the production schedule, and respecting the maintenance equipment periodicity.

A machine j can be subject to several different maintenance tasks which will be repeated periodically. Let be:

-M ij : the maintenance task i on the machine j.

-T ij: periodicity of the maintenance task M ij .

-Tmin ij : earliest time separating two consecutives occurrences of the maintenance task M ij ; -Tmax ij : latest time separating two consecutives occurrences of the maintenance task M ij ; -p' ij : processing time of the task M ij . It is supposed to be known and constant.

The earliness and the tardiness of the k th occurrence of the maintenance task M ij is computed as follows:

-t' ijk : execution time of the k th occurrence of the maintenance task M ij .

-E' ijk : Earliness of the k th occurrence of the maintenance task M ij . E' ijk = max (0, t' ijk +p' ij +T minij -t' ijk +1) -L' ijk : Tardiness of the k th occurrence of the maintenance task M ij . L' ijk = max (0, t' ijk+1 -t' ijk -p' ij -T maxij )

The objective of optimization must be a compromise between the target objective maintenance and production functions.

The production objective function is the completion time of the last job on the last machine (C max ). One will note f 1 this function: 1 (1) The maintenance objective function respects the maintenance periods which influence the constraints of the production system. One will note f 2 this function:

f 1 =C max =Max(c ij )
jk m j Maxj k jk L E f ' ' 1 1 2 + = ∑ ∑ = = 2 (2)
To optimize the two criteria, we taken into account the following common objective function:

f obj =α f 1 +β f 2 (3)
(α,β) are weights which will measure the respective contributions of production and maintenance in the common objective function. They can depend on the number of tasks, processing time or problem size. They are independent and non complementary.

III. DEVELOPPEMENT OF A NEW INTEGRATED ACO ALGORITHM (INTACO)

We now present the development of the new ant-colony algorithm, called INTACO, which extends the ideas of Ant System (AS) algorithm by Dorigo and al. [START_REF] Dorigo | The ant system: Optimization by a colony of cooperating agents[END_REF]. INTACO incorporates simultaneously production jobs and maintenance tasks in the construction of a joint ant sequence, uses new mutation operators based on those proposed in [START_REF] Taillard | Some efficient heuristic methods for the flowshop sequencing problem[END_REF][START_REF] Osman | Simulated annealing for permutation flow shop scheduling[END_REF][START_REF] Dannenbring | An Evaluation of Flow Shop Sequencing Heuristics[END_REF] adapted to joint production and preventive maintenance sequence allowing the generation of the joint ant sequence neighbourhood, and a new local search technique in the AS algorithm.

In this section, we will first present the proposed representation of a joint schedule. Then we introduce the new integrated ACO algorithm (INTACO)

A. Representation of a joint schedule

We propose to code a joint schedule as a structure with two fields: the first one is a sequence S that represents the production jobs execution order. The second one is a matrix M that represents the maintenance tasks insertion sites. The element M[i,j] represents the insertion of the i th maintenance task of the j th machine in the sequence S.

Example:

The execution of the tasks on the three machines according of the preceding example is the following: Machine 1: M 01 ,P 1 ,M 1 1 ,P 9 ,P 3 ,P 8 ,M 21 ,P 5 ,P 6 ,M 31 ,P 7 ,P 4 ,P 2 Machine 2 : P 1 ,M 02 ,P 9 ,M 12 ,P 3 ,P 8 ,P 5 ,M 22 ,P 6 ,P 7 ,P 4 ,P 2 Machine 3 : M 0 3 ,P 1 ,P 9 ,P 3 ,P 8 ,M 13 ,P 5 ,P 6 ,P 7 ,M 23 ,P 4 ,M 33 ,P 2 B. Construction of a joint ant-sequence and its improvement by a local search scheme We model the problem as a search for the best path in a particular graph G <N, V> where N is the set of nodes representing joint solutions, and V represents the difference between the common objective function referred to in formula 3 of two connected solutions (Si, Sj).

Thus, the integrated strategy consists of simultaneously scheduling both maintenance and production activities. The first cycle of the algorithm commences with joint sequence generated randomly and proceeds by construction, moving from a solution to another. A number of ants, which may be thought of as independent agents, are set in motion, each ant representing a joint sequence.

Each ant i starts from a joint initial solution S 0 i which is randomly generated. At each step, a set of new joint solutions S k i is generated by applying mutation to the current one, which represents its neighbourhood N(S 0 i ). During her move, each ant will leave trail of pheromone on her path.

Each ant starts the search with a random solution. For us a solution can be complete or partial. A complete solution represents a joint production and preventive maintenance schedule. In this case, a complete initial solution is obtained after the insertion of maintenance tasks, using one of the heuristics developed by Benbouzid and al. [START_REF] Benbouzid | Resolution of joint maintenance/production scheduling by sequential and integrated strategies[END_REF], on a production schedule generated randomly.

The move from a solution to another in the neighbourhood can be done by mutation in the production sequence, the maintenance matrix, or on both at the same time. We define two types of mutation allowing generating neighbour solutions from the current one. The first one is relative to new mutation reported by Benbouzid and al. [START_REF] Benbouzid | New Genetic Operator for Solving the Joint Production and Maintenance Scheduling: Application to the Flow Shop Problem[END_REF] and the second one is based on new ideas that are developed in the current work.

(a) Benbouzid and al. [START_REF] Benbouzid | New Genetic Operator for Solving the Joint Production and Maintenance Scheduling: Application to the Flow Shop Problem[END_REF] developed a couple of new mutation: Random mutation and k-points vertical mutation. Random mutation consists in shifting randomly towards the left or the right one or more maintenance tasks. The principle of the k-points vertical mutation consists on generating randomly k positions in the production sequence S, thus defining k+1 groups. Then the sites of the production jobs which are in each group i are permuted with the group k+2-i. However the maintenance tasks sites inside each group are preserved.

(b) We propose two new mutations called IIM (Integrated Insertion-Moves) and ISM (Integrated Swap-Moves). It consists in performing the classical insertion-moves and swap-moves [START_REF] Taillard | Some efficient heuristic methods for the flowshop sequencing problem[END_REF][START_REF] Osman | Simulated annealing for permutation flow shop scheduling[END_REF][START_REF] Dannenbring | An Evaluation of Flow Shop Sequencing Heuristics[END_REF] on the production sequence S while updating the maintenance matrix M under the changes operate on the sequence S.

Having thus generated a joint ant sequence of production jobs and maintenance tasks, some changes are applied in order to improve the solution. Local search for the FSP usually starts from some initial feasible sequence and then tries to improve it by small changes. In the case of the joint production and maintenance scheduling problem in permutation flowshop, it is to be noted that the 1-point vertical mutation suggested by Benbouzid and al. [START_REF] Benbouzid | New Genetic Operator for Solving the Joint Production and Maintenance Scheduling: Application to the Flow Shop Problem[END_REF], has emerged as the best improvement for all benchmark problems tested.

C. Updating of trail intensities

Starting from a joint initial sequence, INTACO makes use of trail intensities in order to determine the solution to be appended next step.

At each step, a set of new solutions is generated by applying mutations on the current one, as defined above. When the ant (k) is on one solution (S k i ), it must choose a solution S k j among the set of solutions of its neighbourhood N(S k i ). This choice is made on the basis of the heuristic information (η ij ), and the pheromone information (τ ij) , according to the following rule [START_REF] Bonabeau | Swarm Intelligence: From Natural to Artificial Systems[END_REF]:

(4)
The heuristic information is computed using SPT (Shorting Processing Time) algorithm as follows:

∑ + = r jr ij p 1 1 η 3 (5)
When an ant move to a new solution, chosen in the neighbourhood of the current one, a local update of the pheromone trail is made. When all ants finished a cycle, a global update of the pheromone trail is done on the path of the best joint solution found in the cycle. The best joint sequence so far found is retained, and a new cycle is started. The pheromone update rule is applied as follows:

τ ij new = (1-ρ) * τ ij old + ρx1/f obj.

4

(6) Once all ants have built their solutions, extra pheromone is added on the path taken by the ant which found the best 3 pjr is the processing time of job j on machine r 4 fobj is the objective function referred in formula 3, of the solution given by the best ant.

k i S ) ( ( ( ( ) ) ) ) = = = = t p k ij (τ ij (t)) α * η ij β ∑ l∈J k i ( τ il (t)) α * η ij β 0 solution.
Pheromone is added only when the common objective function f obj of a solution is improved after mutations. Since the neighbourhood of the current solution is a set of complete joint solutions generated according to the operators detailed above, the construction of the pheromone and heuristic information matrix will be done dynamically, as describe below.

Pheromone and heuristic information

IV. EXPERIMENTAL RESULTS

In this section, we present the results of a series of computational experiments conducted to test the effectiveness of the INTACO, for an integrated resolution of the joint production and preventive maintenance scheduling problem in permutation flowshops. The test problems for evaluating the proposed ant-colony algorithm are generated by following the procedure given by Taillard [START_REF] Taillard | Benchmarks for basic scheduling problems[END_REF] for permutation flowshop scheduling problems.

To our knowledge, there is no instance for the preventive maintenance scheduling problem. Therefore, we developed a processor to generate benchmarks in preventive maintenance. The used parameters are: the number of machines and maintenance task parameter (T, T min and T max ). We generated only one maintenance type per machine for each problem. Moreover, the processing time of a maintenance task is identical for all its occurrences. The objective functions taken into consideration are ( § 2): the production objective function f 1 is the minimization of makespan; the maintenance objective function f 2 is the minimization of the sum of maintenance tasks advances and delays. The contribution of production and maintenance objective functions in the common objective function f obj are equal to 1 (α =1 and β =1). The proposed common weighted global objective function will allow tackling the problem in a simplified way.

For evaluating the performance of the INTACO with respect to the common objective function defined in section 2, the results of a research study by Benbouzid and al. [4] are referred to. Benbouzid and al. developed a couple of new heuristic and metaheuristic methods (NEH heuristic, Taboo search and GA) considering the sequential and integrated strategies. They considered the benchmark problems of Taillard [START_REF] Taillard | Benchmarks for basic scheduling problems[END_REF], and reported the best-heuristic solutions for these benchmark problems with respect to the common objective function. It is to be noted that the integrated GAs has emerged to be the best for all benchmark problems.

A. Parameters tuning

Usually, the first set of tests defines for which values of {α, β} the algorithm converges to the known optimum of a selected problem. In our case, no competitive results for the joint production and preventive maintenance scheduling problems exist, with state-of-art algorithms. For this reason, we use the Taillard PFSP 20x10 [START_REF] Taillard | Benchmarks for basic scheduling problems[END_REF]. However, we consider the relative increases in the common objective function yielded by the proposed ant-colony algorithm, relative to the best makespan value reported by the Taillard PFSP 20x10 with upper bound stated at 1659, when comparing the solutions. The best one is the one which gave the least increase.

Table 1 illustrates some of the simulations that have been carried out, using ant colonies, with different sets of {α, β} parameters. Pheromone was stored on nodes. We don't include the case α =0 and β=0 because every ant simply branches randomly from one node to the other, for which the algorithm is expected not to find the optimum. 1 shows that the set of values for which an optimum value for the common objective function is reached is α= 4 and β =1.

Once established the values of α and β, the influence of the parameters that define the pheromone update can be investigated. Figure 1 shows the performance of the algorithm for different values of the evaporation constant ρ. From the case of ρ>0.1, the algorithm does not always find the optimum. This is due to the fact that when the evaporation rate is too high, only the edges belonging to the path with the early encountered local minimum receive pheromone update every cycle. The rest of the edges have their pheromone rapidly decreasing to zero with increasing number of cycles, due to the evaporation. At each new cycle, this effect is amplified and the ants will be more and more attracted to the edges that receive pheromone update. The ants will be searching in the small neighbourhood of the local minimum.

We noticed that the simulations that gave us the best results were the ones where the number of ants equals the number of jobs. Therefore, all ants start the resolution with joint schedules where the first job is different from one solution to another. All ants are scattered on the solutions space.

Figure 2 shows the performance of the algorithm for different values of the maximum number of cycles NC max . The figure 2 shows that as the number of cycle's increases, the objective function improves until the number of cycles reaches 600 cycles. Above this number, the common objective function increases. Thus, the maximum number of cycles can be state at 600.

B. Performance analysis of the INTACO

After optimization of the parameters in the previous section for the specific Taillard 20x10 flowshop problem, this section shows some simulation results of the INTACO to some more complex benchmark problems, with the values of parameters tuned above.

In ACO, as in many other metaheuristics, the performance critically depends on the appropriate combination of local search with the solution construction mechanism. Therefore, we performed a detailed analysis of the effectiveness of the INTACO without and with the improvement procedure proposed in section III.5 at two levels: on the solution given by each ant (L1) and on the best solution given by the algorithm (L2). The results of evaluation are presented in Table 2. Overall, it is observed from the results that the improvement procedure at L1 performs better than this at L2. The main reasons for such a superior performance are due to the fact that the solution of each ant is improved, at each cycle.

We defined, in section III.3, two types of mutation allowing generating neighbour solutions from the current one. The evaluation of the performance of the proposed mutation is presented in table 3.

Evaporation constant Common objective function

Overall, it is observed from the results that the 1-point vertical mutation gives the best results. This higher performance is mainly due to the effectiveness of mutation on the production sequence while respecting the maintenance tasks initial position. This factor assumes greater significance in enhancing the performance of the obtained joint schedule in the case of mutation with 1-point than in the case of two points.

Table 4 shows the comparison between the INTACO to another meta-heuristic: integrated GAs developed in [START_REF] Benbouzid | Resolution of joint maintenance/production scheduling by sequential and integrated strategies[END_REF]. GAs solutions are obtained after 100 executions of the method. The best result is saved, as well as the associated parameters. The following parameters are the same for all the executions of the genetic algorithms: crossover rate: 0.7; mutation rate: 0.01; renewal strategy: N_best, the replacement is done between the selected population and the crossed one. Population size: between 50 and 100. The crossover and mutation operators used are: k-point horizontal crossover and k-point vertical mutation.

However, since the GAs are known to prematurely converge, the only way to get the best solution was to use a strategy called "sharing" [START_REF] Bonabeau | Swarm Intelligence: From Natural to Artificial Systems[END_REF]. This strategy allows diversifying the population, avoiding the GAs to get trapped in a local optimum. We use sharing with the following parameters: α Sharing = 0.99 and β Sharing = 4. Instead the use of "sharing" the INTACO performs better than the integrated GAs for most of the studied benchmarks.

Conclusion

Compared to other metaheuristics such as genetic algorithms, simulated annealing and taboo search, relatively few attempt have been made to solve scheduling problems using ant-colony algorithms.

In this paper, we investigated the problem of the joint production and preventive maintenance scheduling in permutation flowshop using ant-colony algorithm. We proposed a new algorithm, called INTACO, which tries to solve this problem, based on a joint representation of production and maintenance data. The goal is to optimize a common objective function which takes into account both maintenance and production criteria.

For evaluating the performance of the proposed ant-colony algorithm, the results of a recent study by Benbouzid and al. [4] are referred to. The proposed algorithm is enhanced by a local search procedure to yield high quality solutions. The proposed ant-colony algorithm performs better than the integrated genetic algorithms referred to. Future research will be to investigate a bi-criterion approach to solve this problem.

Fig. 1 .

 1 Fig.1. Influence of different setting of the evaporation coefficient to reach the optimum. (NCmax=1000, α=4, β=1, m=50).

Fig. 2 .

 2 Fig. 2. Influence of different setting of the average number of cycles to reach optimum ( α=4, β=1, ρ=0.1).

TABLE I OPTIMUM

 I FOR DIFFERENT VALUES OF α AND β. RESULTS FOR EACH PAIR (α ,β) ARE GIVEN WITH MAXIMUM NUMBER OF CYCLES NCMAX=1000, ρ =0.1, m=20.

		β = 1	β = 2	β = 3	β = 4
	α = 1	1756	1764	1793	1788
	α = 2	1754	1772	1780	1791
	α = 3	1747	1744	1755	1770
	α = 4	1740	1755	1749	1762
	α = 5	1755	1742	1747	1789
	α = 6	1769	1771	1786	1792
	α = 7	1792	1769	1777	1755
	α = 8	1818	1801	1767	1771
	α = 9	1801	1783	1793	1755
	α = 10	1797	1793	1796	1774
	Table

TABLE IV COMMON

 IV OBJECTIVE FUNCTION GIVEN BY THE INTEGRATED GA AND INTACO

	Benchmark	GA	INTACO*	INTACO
	7×7	6898	7436	7375
	8×8	8720	9396	9201
	15×5	605	634	598
	20×5	1222	1261	1221
	20×10	1686	1740	1614
	20×15	2203	2385	2002
	30×15	2745	3092	2632
	50×10	3698	3777	3183
	75×20	5851	6308	5320
	100×10	571	654	528

*: without the improvement procedure