A Simple Expression for the Tortuosity of Gas Transport Paths in Solid Oxide Fuel Cells’ Porous Electrodes

Based on the three-dimensional (3D) cube packing model, a simple expression for the tortuosity of gas transport paths in solid oxide fuel cells’ (SOFC) porous electrodes is developed. The proposed tortuosity expression reveals the dependence of the tortuosity on porosity, which is capable of providing results that are very consistent with the experimental data in the practical porosity range of SOFC. Furthermore, for the high porosity (>0.6), the proposed tortuosity expression is also accurate. This might be helpful for understanding the physical mechanism for the tortuosity of gas transport paths in electrodes and the optimization electrode microstructure for reducing the concentration polarization.


Introduction
Emerging as an attractive device for future power plants, solid oxide fuel cells (SOFC) enjoy some advantages, such as higher energy conversion efficiency, quiet operation, fuel flexibility, and low environmental hazards.In recent years, significant strides have been made in improving SOFC performance [1][2][3].In other words, for a given operating voltage, the output current density has been increased drastically.For example, the output current density of a single SOFC can reach 3.5 A¨cm ´2 at 0.7 V [4].However, the higher the output current density is, the more serious the concentration polarization is.Thus, the accurate modeling of gas transport in electrodes is essential for the understanding and prediction of the concentration polarization and fuel cell performance.However, it is a challenging task to simulate the gas transport in electrodes accurately because the presence of the solid phase causes the transport paths of gas to deviate from straight lines.In order to take the effect of tortuous gas transport paths into account, the effective gas diffusion coefficient D eff is usually corrected by introducing the tortuosity, which is commonly expressed as [5][6][7][8][9]: where D 0 is the bulk gas diffusion coefficient, ε is the porosity of the electrode, and τ (τ = L eff /L 0 ) is the tortuosity of the gas transport paths, which is the ratio of the average actual path length (L eff ) to the straight distance (L 0 ).It should be noted that the square of tortuosity is equal to the tortuosity factor.Although Equation ( 1) is widely adopted in SOFC literature, the studies on tortuosity are limited, which includes experimental and numerical approaches [10][11][12][13].Wilson et al. [14] reconstructed the three-dimensional (3D) microstructure of an electrode employing the focused ion beam-scanning electron microscope (FIB-SEM) technique, and then the tortuosity is estimated by solving the diffusion equation in the pore phase.Based on the microstructure of the electrode obtained by FIB-SEM, Iwai et al. [15] adopted the random walk calculation and the lattice Boltzmann method (LBM), which are two different methods to calculate tortuosity factor.Although the FIB-SEM experimental method provides valuable information about tortuosity, it is not predictive and relies on specific sophisticated equipment.For the numerical approach, the 3D sphere packing model is widely adopted.The electrode microstructure is modeled by a set of the random packing spheres, and is used as the geometry structure for the transport model.The tortuosity can be quantitatively evaluated by solving the transport equation [16,17].However, the tortuosity value predicted by the 3D sphere packing model is smaller than the experimental data.In order to address this issue, the 3D cube packing model was proposed by the author [11].In this approach, the electrode representative volume element (RVE) consists of a cubic lattice (N x ˆNy ˆNz where lattice units, N x , N y and N z are respectively the number of lattice units in x, y and z directions) in which ε ˆNx ˆNy ˆNz lattice units are randomly selected as the pore-forming pore phase in the electrode, as depicted in Figure 1.Baseed on this new electrode microstructure, the tortuosity can be obtained by solving the diffusion equation.The good agreement with experimental data verifies the accuracy of this method.
Energies 2015, 8, page-page 2 electron microscope (FIB-SEM) technique, and then the tortuosity is estimated by solving the diffusion equation in the pore phase.Based on the microstructure of the electrode obtained by FIB-SEM, Iwai et al. [15] adopted the random walk calculation and the lattice Boltzmann method (LBM), which are two different methods to calculate tortuosity factor.Although the FIB-SEM experimental method provides valuable information about tortuosity, it is not predictive and relies on specific sophisticated equipment.For the numerical approach, the 3D sphere packing model is widely adopted.The electrode microstructure is modeled by a set of the random packing spheres, and is used as the geometry structure for the transport model.The tortuosity can be quantitatively evaluated by solving the transport equation [16,17].However, the tortuosity value predicted by the 3D sphere packing model is smaller than the experimental data.In order to address this issue, the 3D cube packing model was proposed by the author [11].In this approach, the electrode representative volume element (RVE) consists of a cubic lattice (  As a key factor for eff D determination, tortuosity is a critical input parameter in SOFC multiphysics models of the cell level.However, there is still a lack of an effective method to determine the tortuosity in the entire porosity scope of SOFC electrodes, which leads to a big disagreement about tortuosity.As a result, the tortuosity is usually set as an empirical value randomly or used as an adjustable parameter to fit experimental data in modeling studies [17,18], which may nullify the subsequent modeling result.
Based on the 3D cube packing model, the objective of this paper is to derive a simple and effective expression for tortuosity calculation and eliminate the uncertainty of tortuosity.The derived expression not only reveals the gas diffusion mechanism in porous electrodes, but also provides a trustworthy tortuosity evaluation criterion for SOFC modeling, which will be helpful to optimize the electrode microstructure and improve SOFC performance.

Theory
For the 3D cube packing model, in the arbitrary layer (perpendicular to z direction), there is the ε lattice unit for the pores due to the random pore distribution in electrodes [11].Likewise, the location relationship among pores in the two arbitrary adjacent layers is also similar.Consequently, the two arbitrary adjacent layers are selected as the objects, as shown in Figure 2a, which are named the n layer and the n + 1 layer, respectively, for the above layer and the below layer.Gas transporting from the n layer to the n + 1 layer mainly follows two typical paths.Figure 2b displays the first typical path, which is the simplest situation.The lattice unit in the n layer is porous, As a key factor for D eff determination, tortuosity is a critical input parameter in SOFC multiphysics models of the cell level.However, there is still a lack of an effective method to determine the tortuosity in the entire porosity scope of SOFC electrodes, which leads to a big disagreement about tortuosity.As a result, the tortuosity is usually set as an empirical value randomly or used as an adjustable parameter to fit experimental data in modeling studies [17,18], which may nullify the subsequent modeling result.
Based on the 3D cube packing model, the objective of this paper is to derive a simple and effective expression for tortuosity calculation and eliminate the uncertainty of tortuosity.The derived expression not only reveals the gas diffusion mechanism in porous electrodes, but also provides a trustworthy tortuosity evaluation criterion for SOFC modeling, which will be helpful to optimize the electrode microstructure and improve SOFC performance.

Theory
For the 3D cube packing model, in the arbitrary layer (perpendicular to z direction), there is the ε ˆNx ˆNy lattice unit for the pores due to the random pore distribution in electrodes [11].Likewise, the location relationship among pores in the two arbitrary adjacent layers is also similar.Consequently, the two arbitrary adjacent layers are selected as the objects, as shown in Figure 2a, which are named the n layer and the n + 1 layer, respectively, for the above layer and the below layer.Gas transporting from the n layer to the n + 1 layer mainly follows two typical paths.Figure 2b displays the first typical path, which is the simplest situation.The lattice unit in the n layer is porous, under which the lattice unit in the n + 1 layer is also porous.The cross-section of the first typical path can be determined from the following equation: where S 1 is the cross-section of the first path, and S 0 is the cross-section of RVE.This is because the probability of a lattice unit as a pore is porosity ε.The second typical path is displayed in Figure 2c.
Here there is a pore lattice unit connecting to the first typical path in the n layer, under which the lattice unit is solid.The cross-section of the second typical path S 2 can be written as: Equation ( 3) is complex, and it is derived by following steps.For the second typical path, it has to satisfy the following conditions: (1) There is a first typical path, the probability of which is ε 2 .
(2) In the second typical path a pore lattice unit connects to the first typical path in the n layer, under which the lattice unit is solid in the n + 1 layer.
Energies 2015, 8, page-page 3 under which the lattice unit in the n + 1 layer is also porous.The cross-section of the first typical path can be determined from the following equation: where 1 S is the cross-section of the first path, and 0 S is the cross-section of RVE.This is because the probability of a lattice unit as a pore is porosity ε .The second typical path is displayed in Figure 2c.
Here there is a pore lattice unit connecting to the first typical path in the n layer, under which the lattice unit is solid.The cross-section of the second typical path 2 S can be written as: Equation ( 3) is complex, and it is derived by following steps.For the second typical path, it has to satisfy the following conditions: (1) There is a first typical path, the probability of which is ε 2 .
(2) In the second typical path a pore lattice unit connects to the first typical path in the n layer, under which the lattice unit is solid in the n + 1 layer.A pore lattice unit in the first typical path contacts with four lattice units in the n or n + 1 layer, which are labeled as A1, A2, A3 and A4 (B1, B2, B3 and B4) in the n (n + 1) layer respectively, as displayed in Figure 2d.The probability for one of these four lattice units to be a pore and the lattice unit under this lattice unit to be a solid is In the n + 1 layer, a pore lattice unit contacts with 4ε pore lattice units on average, according to the coordination number theory.In order to avoid the appearance of the first typical path, there are only 4 − 4ε lattice units left to possibly form the second typical path, because the probability of the first typical path has been included in ε 2 .For example, the B1 lattice unit is a pore for ε = 0.25.As shown in Figure 2d, if the A1 lattice unit is also a pore, the first typical path will be formed, which is not allowed.Therefore, A2, A3 and A4 lattice units are possible pores, except for the A1 lattice unit.
Gas transported from the top layer to the bottom layer follows many tortuous paths, each of which consists of many of the first and second typical paths.As a consequence, the average length of these tortuous paths (Leff) is determined by: where 0 L is the length of RVE along the transport direction.A pore lattice unit in the first typical path contacts with four lattice units in the n or n + 1 layer, which are labeled as A 1 , A 2 , A 3 and A 4 (B 1 , B 2 , B 3 and B 4 ) in the n (n + 1) layer respectively, as displayed in Figure 2d.The probability for one of these four lattice units to be a pore and the lattice unit under this lattice unit to be a solid is rpε ´ε2 q{p1 ´ε2 qsrp1 ´εq{p1 ´ε2 qs " ε{p1 `εq 2 .
In the n + 1 layer, a pore lattice unit contacts with 4ε pore lattice units on average, according to the coordination number theory.In order to avoid the appearance of the first typical path, there are only 4 ´4ε lattice units left to possibly form the second typical path, because the probability of the first typical path has been included in ε 2 .For example, the B 1 lattice unit is a pore for ε = 0.25.As shown in Figure 2d, if the A 1 lattice unit is also a pore, the first typical path will be formed, which is not allowed.Therefore, A 2 , A 3 and A 4 lattice units are possible pores, except for the A 1 lattice unit.
Gas transported from the top layer to the bottom layer follows many tortuous paths, each of which consists of many of the first and second typical paths.As a consequence, the average length of these tortuous paths (L eff ) is determined by: where L 0 is the length of RVE along the transport direction.On average, the total cross-section area of these tortuous paths (S t ) may be evaluated by the following equation: The molar flux of gas may be written as: where N is the molar flux, D is the diffusion coefficient of gas, c up (c bot ) is the molar concentration of gas at the top (bottom) surface of RVE.

Results and Discussion
Figure 3 shows the tortuosity for different porosities established by different approaches.Clearly, the Equation ( 9) results agree very well with the experiment results, indicating that Equation ( 9) is capable of reliably predicting the tortuosity of gas transport paths in SOFC porous electrodes.However, the tortuosity predicted by the 3D sphere packing model is much smaller than the experiment results.The reasons can be found in our previous study [11].
Energies 2015, 8, page-page 4 On average, the total cross-section area of these tortuous paths (St) may be evaluated by the following equation: The molar flux of gas may be written as: where N is the molar flux, D is the diffusion coefficient of gas, cup ( cbot) is the molar concentration of gas at the top (bottom) surface of RVE.
Based on the effective media theory, the molar flux of gas can be rewritten as: (7) Comparing Equations ( 6) and (7) gives: The formula of tortuosity is obtained by comparing Equations ( 1) and ( 8), which is estimated as:

Results and Discussion
Figure 3 shows the tortuosity for different porosities established by different approaches.Clearly, the Equation ( 9) results agree very well with the experiment results, indicating that Equation ( 9) is capable of reliably predicting the tortuosity of gas transport paths in SOFC porous electrodes.However, the tortuosity predicted by the 3D sphere packing model is much smaller than the experiment results.The reasons can be found in our previous study [11].) are widely used in SOFC models.Berson et al. [16] reported that Bruggeman correction, which was only valid for ε 0.6  , introduced big errors in low  9), solid square symbols denote the results calculated by the 3D sphere packing model [17].
Energies 2015, 8, page-page porosity.For ε 0.6  , the Equation ( 9) results approximately equate to the Bruggeman correction results, which demonstrates that Equation ( 9) is also suitable for high porosity.).
However, the relationship between porosity and tortuosity is not given in Equation ( 1) (the parallel-capillary model).The tortuosity is often used as an adjustable parameter to fit experimental data.In Figure 4, the tortuosity factor is four for the parallel-capillary model, which is independent of porosity.As a result, f(ε,τ) is linear to the porosity.For ε 1  , the value of f(ε,τ) is only 0.25 according to the parallel-capillary model.This contradicts with the fact that f(ε,τ) equates to 1 for ε 1  .

Summary
In the present study, a simple expression for the tortuosity of gas transport paths in SOFC porous electrodes is derived, which shows the dependence of the tortuosity on porosity.The proposed tortuosity expression is suitable for not only the practical porosity range of SOFC, but also high porosity (ε > 0.6).Consequently, the proposed tortuosity expression might be helpful for the optimization electrode microstructure to reduce the concentration polarization and further improve SOFC performance.S 0 Cross-section of representative volume element S 1 Cross-section of the first typical path S 2 Cross-section of the second typical path S t Total cross section area of the tortuous paths N Molar flux c up pc bot q Molar concentration of gas at the top (bottom) surface of representative volume element the number of lattice units in x, y and z directions) in which ε pore-forming pore phase in the electrode, as depicted in Figure1.Baseed on this new electrode microstructure, the tortuosity can be obtained by solving the diffusion equation.The good agreement with experimental data verifies the accuracy of this method.

Figure 1 .
Figure 1.Schematic of solid oxide fuel cells (SOFC) electrode microstructure; blue is pore phase, gray denotes solid phase.

Figure 1 .
Figure 1.Schematic of solid oxide fuel cells (SOFC) electrode microstructure; blue is pore phase, gray denotes solid phase.

Figure 2 .
Figure 2. Sketch map for the derivation of the 1 S and 2 S : (a) schematic of the location relationship among pores in the two arbitrary adjacent layers; (b) schematic of the first typical path; (c) schematic of the second typical path; (d) schematic of the lattice units connecting to the first typical path.

Figure 2 .
Figure 2. Sketch map for the derivation of the S 1 and S 2 : (a) schematic of the location relationship among pores in the two arbitrary adjacent layers; (b) schematic of the first typical path; (c) schematic of the second typical path; (d) schematic of the lattice units connecting to the first typical path.

Figure 3 .
Figure 3.The tortuosity for different porosities.The open circle symbols denote experiment results [15,19-27], line + solid triangle symbols denote the results calculated by Equation (9), solid square symbols denote the results calculated by the 3D sphere packing model[17].

Figure 3 .
Figure 3.The tortuosity for different porosities.The open circle symbols denote experiment results [15,19-27], line + solid triangle symbols denote the results calculated by Equation (9), solid square symbols denote the results calculated by the 3D sphere packing model[17].