TOPOLOGY OPTIMIZATION APPLIED TO DESIGN OF SOLID OXIDE FUEL CELLS By Xiankai Song A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering--Doctor of Philosophy 2014 ABSTRACT TOPOLOGY OPTIMIZATION APPLIED TO DESIGN OF SOLID OXIDE FUEL CELLS By Xiankai Song A topology optimization method is used to identify optimal designs of the cathode microstructure and the anode support in a solid oxide fuel cell (SOFC). Two dimensional and three dimensional models are considered. A 2D topology optimization model is developed to minimize the cathode resistance. A simplified analysis model is used in computations. Results highlight the importance of the cathode geometry in the performance. Optimal geometric features are found to depend on the material properties and various geometric parameters. To improve upon the accuracy available from a purely 2D model, a 3D finite element model is established to make an accurate prediction of the cathode resistance. A detailed 3D microstructure is reconstructed from images obtained using the 3D focused ion beam-scanning electron microscopy. A 3D topology optimization formulation is set up to minimize cathode resistance. The effect of the material properties on the geometric features is investigated. Improvements up to 35% are achieved by properly organizing the cathode microstructure. The thermal stress problem of the anode support in the SOFC stacks is also of great interest. Fuel flow, heat transfer, thermo-mechanical and electrochemical processes are considered in a coupled model. A Weibull distribution evaluating the probability of failure is used as a measure of the strength of the anode. A new material model is developed aiming at the topology optimization of the anode strength. Optimal designs for two types of objective functions, including minimum thermal compliance and minimum probability of failure, are obtained. It is observed that the designs obtained using the two objective functions can improve the performance for 10%. ACKNOWLEDGMENTS I have been supported and encouraged by many colleagues and friends throughout my PhD time at Michigan State University. I would like to acknowledge those who granted their precious help to me. Foremost, I would like to express my deepest gratitude to my academic advisor, Prof. Alejandro R. Diaz. He has taught me everything I need in research, and make me well prepared for my future work. He has provided me the best research environment I have ever had. He has taken me around the world to meet the professionals in the optimization groups. He has also been an exceptional friend and a great mentor. During my three years’ PhD work, I had the honor to work with Prof. Andre Benard, who has kept me ecstatic with his numerous insightful and creative suggestions. I am very appreciative of Prof. Jason Dale Nicholas for the shared knowledge of solid oxide fuel cells and enthusiasm towards our collaborative work. Extended thanks to my fellow students, Dr. Kazuko Fuchi, Mark Daniel Gaustad, Joonho Lee, Tingli Cai, Smruti Panigrahi and Chao Cheng for their continuous encouragement and invaluable friendship. Special thanks to Prof. Ronald Averill, who agreed to serve on my Ph.D. committee. Finally, I am grateful to my mother, Prof. Naiwen Ye, for always standing by my side, understanding me and motivating me towards a meaningful life. iii TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................ vi LIST OF FIGURES .................................................................................................................... vii KEY TO SYMBOLS OR ABBREVIATIONS ........................................................................... ix Chapter 1 Introduction ................................................................................................................. 1 Chapter 2 A 2D model for shape optimization of solid oxide fuel cell cathodes ...................... 7 2.1 Simplified analysis model of infiltrated composite cathode ............................................. 7 2.2 The simple model as a periodic single column reference design ...................................... 9 2.3 Finite element model for 2D topology optimization problem .........................................11 2.4 Objective function ........................................................................................................... 15 2.5 Constraints ...................................................................................................................... 15 2.6 The optimization problem ............................................................................................... 16 2.7 Filters .............................................................................................................................. 17 2.8 Techniques to reject inadmissible designs ...................................................................... 17 2.9 Sensitivity Analysis ......................................................................................................... 20 2.10 Examples ....................................................................................................................... 21 2.10.1 Designs with different material properties ......................................................... 22 2.10.2 Design with varying footprint B and constant perimeter/footprint ratio ........... 24 2.10.3 Design with varying perimeter P and amount of material W while keeping the footprint constant .......................................................................................... 25 2.10.4 Design with varying amount of material ........................................................... 25 2.11 Conclusions ................................................................................................................... 26 Chapter 3 A 3D finite element model for cathode microstructure of solid oxide fuel cell cathodes .................................................................................................................... 28 3.1 Analysis model of infiltrated composite cathode ............................................................ 28 3.2 Finite element formulation of 3D cathode microstructure .............................................. 31 3.3 Generation of 3D microstructure mesh from experimental data .................................... 33 3.4 Examples of the simulation of the 3D cathode resistance .............................................. 38 3.5 Experimental validation and comparison with a simple estimation model .................... 40 3.6 Other examples ............................................................................................................... 43 3.7 Conclusions ..................................................................................................................... 45 Chapter 4 Topology optimization for 3D cathode microstructure of solid oxide fuel cell cathodes .................................................................................................................... 46 4.1 Problem Setup ................................................................................................................. 46 4.2 Material model ................................................................................................................ 48 4.3 Constraints ...................................................................................................................... 52 4.4 Optimization Problem ..................................................................................................... 53 4.5 Sensitivity Analysis ......................................................................................................... 53 4.6 Numerical Examples ....................................................................................................... 55 iv 4.6.1 Designs with varying the material properties by changing work temperature .... 57 4.6.2 Design with varying the arrangement of the cathode roots with constant footprint............................................................................................................... 63 4.7 Conclusions ..................................................................................................................... 66 Chapter 5 Thermal strength design for anode support structure of solid oxide fuel cell electrode ................................................................................................................... 68 5.1 Analysis model ................................................................................................................ 69 5.1.1 The fuel flow model ............................................................................................. 72 5.1.2 The electrochemical reaction boundary condition ............................................... 73 5.1.3 The hydrogen/heat transport model ..................................................................... 76 5.1.4 The thermal structural model ............................................................................... 77 5.2 Material model ................................................................................................................ 79 5.3 Finite element model....................................................................................................... 83 5.4 Objective functions ......................................................................................................... 84 5.5 Constraints ...................................................................................................................... 84 5.6 Optimization problems.................................................................................................... 85 5.7 Sensitivity analysis.......................................................................................................... 87 5.8 Examples ......................................................................................................................... 89 5.8.1 Minimum the thermal compliance ....................................................................... 90 5.8.2 Minimum the probability of failure ..................................................................... 91 5.9 Conclusions ..................................................................................................................... 93 Chapter 6 Conclusions ................................................................................................................ 95 BIBLIOGRAPHY ....................................................................................................................... 98 v LIST OF TABLES Table 2.1 Performance of optimized designs for different materials, where R can be multiplied by l/ 0 to recover units (Fig. 2.6) .............................................................. 23 Table 2.2 Performance of optimized designs for different footprints B and PB/B = 9, where R can be multiplied by l/ 0 to recover units (Fig. 2.7) ............................................... 24 Table 2.3 Performance of optimized designs for different values of the perimeter bound PR, where R can be multiplied by l/ 0 to recover units (Fig. 2.9) .................................... 25 Table 2.4 Performance of optimized designs for different values of the amount of material bound WR, where R can be multiplied by l/ 0 to recover units (Fig. 2.9) .................. 26 Table 3.1 Bulk conductivity of CGO (Gd0.1Ce0.9O2) .................................................................. 39 Table 3.2 Scaled ASR values for 20 nm MIEC LSCF (La0.6Sr0.4Co0.8Fe0.2O3) particle diameter....................................................................................................................... 39 Table 3.3 Diffusivity and concentration of oxygen gas in air at 1atm [57] .................................. 40 Table 4.1 Performance of optimized design at different temperatures (Fig. 4.7-4.8) ................... 58 Table 4.2 Tortuosity and surface areas of optimized design at different temperatures (Fig. 4.7) .............................................................................................................................. 62 Table 4.3 Performance of optimized design for different root arrangements (Fig. 4.9-4.10) ....... 66 Table 4.4 Improvement of optimized design for different root arrangements compared to experimental data (Fig. 4.9-4.10) ................................................................................ 66 Table 5.1 The performance of the anode support design for the thermal compliance design ...... 91 Table 5.2 The performance of the anode support design for the probability of failure design ..... 92 vi LIST OF FIGURES Figure 1.1 Schematic arrangement of a SOFC ............................................................................... 1 Figure 2.1 Geometry of the reference design used for comparison .............................................. 10 Figure 2.2 A periodic computational domain is used in the cathode model as shown in a but it is extended as in b to allow the exploration of the entire space through topology optimization ................................................................................................. 12 Figure 2.3 Sketch of periodic cell. ................................................................................................ 14 Figure 2.4 Detecting neck regions that may trap internal holes ................................................... 18 Figure 2.5 Cathode model may predict that (b) is worse than either (a) or (c)............................. 19 Figure 2.6 Optimized designs for different values of h0/ 0. As h0/ 0 increases, the perimeter becomes less important and the shape of the scaffold begins to dominate. Performances reported in Table 2.1 ........................................................... 23 Figure 2.7 Optimized designs for different footprints B and PB/B = 9. For a large footprint, the cathode shape becomes more convoluted and departs further from a columnar shape, not being confined by the design space. Performances reported in Table 2.2.................................................................................................................. 24 Figure 2.8 Optimized designs for different values of the perimeter bound PR. As expected, a large perimeter corresponds to a greater performance improvement. Performances reported in Table 2.3 ............................................................................ 25 Figure 2.9 Optimized designs for different values of the amount of material bound WR. Performances reported in Table 2.4 ............................................................................ 26 Figure 3.1 The sketch of a 3D cathode microstructure ................................................................. 29 Figure 3.2 The manufacture process of a SOFC composite cathode sample [14] ........................ 33 Figure 3.3 The reconstruction process using 3D Focused Ion Beam-Scanning Electron Microscopy. ................................................................................................................. 35 Figure 3.4 An example of a black-white image of a 9μm thick scaffold. ..................................... 36 Figure 3.5 The 3D mesh of the 9 micron thick cathode microstructure ....................................... 38 Figure 3.6 The comparison of RP data for 20nm LSCF particles. ................................................ 41 Figure 3.7 The comparison of RP data for 20nm LSCF particles with a small oxygen gas vii diffusivity. ................................................................................................................... 42 Figure 3.8 The comparison of RP data for 20nm LSCF particles with a very large oxygen gas diffusivity.............................................................................................................. 43 Figure 3.9 The comparison of RP data obtained from different methods for 20nm LSCF particles ....................................................................................................................... 44 Figure 4.1 The design domain of the 3D microstructure problem ................................................ 48 Figure 4.2 The arrangement of the element face number ............................................................. 50 Figure 4.3 Penalty in the surface resistance model may make design (b) performs worse than (a) and (c) ............................................................................................................ 50 Figure 4.4 The second derivative can determine the trend of the first derivative......................... 52 Figure 4.5 The arrangements of the fixed roots ............................................................................ 56 Figure 4.6 The comparison between the analysis domain in Chapter 3 and the design domain in Chapter 4 (in white rectangle) ................................................................... 57 Figure 4.7 The optimized design at different work temperatures. Performance is reported in Table 5.1 ...................................................................................................................... 59 Figure 4.8 The comparison of RP data for 20nm LSCF particles with a small oxygen gas diffusivity. ................................................................................................................... 61 Figure 4.9 The optimized design for different root arrangements at 400°C. Performance is reported in Table 4.3 ................................................................................................... 64 Figure 4.10 The optimized design for different root arrangements at 700°C. Performance is reported in Table 4.3 ................................................................................................... 65 Figure 5.1 The cross section of the fuel supply channel and the electrode................................... 70 Figure 5.2 The analysis domain (a) and design domain (b) for the anode support ....................... 71 Figure 5.3 The i- curve and the i- curve ............................................................................ 75 Figure 5.4 The material models for thermo-mechanical problem ................................................ 82 Figure 5.5 The reference flat design with a constant thickness of 600μm, the black pixels represent the anode support; the white pixels represent the fuel flow path ................ 90 Figure 5.6 The optimal cross section in example 5.8.1 for the minimum probability problem ... 90 Figure 5.7 The optimal cross section in example 5.8.2 for the minimum probability problem ... 92 viii KEY TO SYMBOLS OR ABBREVIATIONS SOFC Solid Oxide Fuel Cell FIB-SEM Focused Ion Beam-Scanning Electron Microscopy MIEC Mixed Ionic Electronic Conducting Ji 2 atomic flux (A/cm ) electrical potential (V) Faraday’s constant (C/mol) 0 ionic conductivity (S/cm) 2 RS cathode surface resistance (Ωcm ) h0 interface reaction coefficient (S/cm ) 2 ∞ driving potential (V) 2 RP cathode resistance (Ωcm ) total current flow (A) l non-dimensional element size non-dimensional perimeter non-dimensional amount of material non-dimensional unit cell footprint 2 g diffusivity of oxygen gas (m /sec) 2 interface current generation (A/cm ) ix 2 interface oxygen gas consumption (mol/(sec∙cm )) hc oxygen exchange coefficient (S∙cm/mol) z charge number of oxygen ion n0 Avogadro constant (mol ) -1 2 viscosity of hydrogen gas (m /sec) Flow velocity (m/sec) pressure (Pa) 3 Fuel density of fuel (kg/m ) cc V close circuit potential (V) oc V open circuit potential (V) act  activation potential (V) universal gas constant (J/(mol∙K)) temperatrue (K, cH 2  H2 ) 3 concentration of hydrogen gas (mol/m ) 2 diffusivity of hydrogen gas (m /sec) cHeat the specific heat (J/(kg∙K)) 2  Heat heat conductivity (m /sec) hair convection coefficient for air flow (m/sec) x Pf probability of failure E0 Young’s modulus (Pa) displacement (m) strain stress (Pa) -1  t coefficient of thermal expansion (K )  S Weibull’s characteristic strength (Pa) Weibull’s modulus CT thermal compliance (J) CS static compliance (J) total current generation (A) 3 QF total fuel flow rate (m /sec) xi Chapter 1 Introduction A solid oxide fuel cell (SOFC) is an electro-chemical device that converts chemical energy into electric energy and heat directly, without the process of combustion [ 1][ 2 ]. Like other electrochemical devices, SOFCs, illustrated schematically in Fig. 1.1, output electrical power by generating electricity as electrons flow through an external circuit and linking spatially-separated redox reactions. Compared to other types of fuel cells, a distinguishing feature of SOFCs is the use of a solid electrolyte that conducts the oxygen ions from cathode to anode for the completion of the reactions occurring at the electrodes, which permits more flexibility in fuel cell designs [2]. The efficiency of SOFC can be as high as 60%-90% [3]. These advantages make SOFC a focus in recent fuel cell research. Figure 1.1 Schematic arrangement of a SOFC For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. In a SOFC, there are voltage losses associated with ohmic resistances in the electrodes and active 1 polarizations on the air-cathode interface and the fuel-anode interface. However, the voltage losses in a cathode usually dominate the entire stack so it is more efficient to focus on improving the cathode performance [3]. Often, this is accomplished by designing better composite cathodes. In one approach, conducting electrocatalyst nano-particles of a mixed ionic electronic conducting (MIEC) material are deposited into a porous ionic-conducting scaffold to for a composite cathode, and the composite cathode has received a lot of attention [4]-[8] . This is the arrangement for the cathode modeled in chapters 2-4, and the focus is on improving performance by optimizing the ionic conducting scaffold cathode microstructure. Attempts in designing the cathode geometry from macro to micro scales have shown progresses in reducing the cathode resistance. Chan et al.[9] performed a sensitivity analysis for a generic macro-homogeneous polarization model of a solid oxide fuel cell to study the effect of the thickness of the fuel cell components. The performance of composite cathodes as a function of various electrode microstructural parameters has also been the subject of investigation, as they typically show lower overpotential losses than single-phase electrodes [10]-[13]. Recently, detailed micromodeling of the performance of SOFCs was proposed and used to minimize the cathode total resistance [14][15]. In related work [16], the cathode-gas is modeled as a homogenized porous medium, and the author’s goal is to optimize the shape of the interface between this medium and the electrolyte. Attempts were also made by Song et al. in [17]-[19] to improve the cathode performance by designing the micro and macro scales cathode structures. As an electrochemical device that works at a high temperature (500~1000 ℃ ), the thermo-mechanical behavior also plays a significant role in the performance of SOFC. The typical anode-supported electrode design is composed by a thick layer (≈500m) of porous support on the anode side the SOFC electrode and a thin active layer (≈50m ) of cathode, 2 electrolyte and anode are printed on the anode support [20]. The porous media used in SOFC electrodes exhibit a large spread of bulk failure strength due to the random nature of their inherent flaws from manufacturing [21]. The failure stress of these materials is not a well defined number, so the mechanical failure by the probability of failure [22], which was expressed by Weibull, is utilized to assess whether the material will fail at a given stress or not. The mechanical strength of the electrode is dominated by the performance of anode support, and it was investigated experimentally in [23]. Finite element analysis which considers both the electrochemical and thermo-mechanical processes in SOFC was also performed in [24]-[27]. Further works [28][29] evaluated the probability of failure using Weibull method, and concluded that the electrode has a high probability of failure when large tensile stress exists. The anode strength was also been optimized by imposing changes to production parameters, such as powder milling and sintering temperature [30]. Topology optimization is a design methodology that is able to find the optimal geometries of structures that achieve the design goals. It was first proposed by Bendsøe and Kikuchi in the context of structural design [31]. The geometries of structures are represented by an effective density of a material in each element in the finite element models, and the range of the effective density is extended from the binary values to a continuous domain. The material properties of each element are computed using a homogenization method and used in the constitutive equations. After progressing for more than two decades, the use of topology optimization has been extended to material design [32] with applications in different physics [33]. The formulation of the simulation model is of heat transfer type, and the application of topology optimization in heat transfer type problems is extremely interesting, especially when the problem involves design dependent boundaries. In such problems it is not possible to define material 3 boundaries a priori — they are the subject of the optimization — and therefore some boundary conditions cannot be set directly. Problems such as this have been discussed by various authors in the context of heat conduction. For instance, Yoon and Kim [34] proposed an “element connectivity parameterization” method to consider design-dependent effects for thermal boundary conditions. Bruns [35] investigated a 2D problem with convection both on the top and one the boundary of the structure. Iga et al.[36] introduced a smeared-out hat function to extract the side convection. Due to the strong demand in industry, thermo-mechanical problems also attract many concerns. The earliest research in the area of topology optimization with thermal loads was by Rodrigues and Fernandes, who used a homogenization method to minimize the compliance of structures with combined temperature and mechanical loading [37]. Li et al. [38] applied the evolutionary structural optimization (ESO) method with element thicknesses as design variables. Using topology optimization method, Sigmund and Torquato generated structures with extreme thermal expansion properties [39], and Sigmund also developed thermal micro-actuators [40]. Cho and Choi [41] performed an adjoint design sensitivity analysis for topology optimization of weakly coupled thermo-elastic problem. In recent years, thermo-elastic topology optimization has been performed using to discuss the effect of the interpolation scheme in density-based methods for thermal loading [42]. For the thermo-strength design, Pedersen and Pedersen proposed an objective function based on uniform energy density instead of the structural compliance, and produced superior results focusing on strength design [43][44]. This dissertation focuses on two aspects of the performance of SOFC. Chapters 2-4 investigate the 2D and 3D finite element modeling and topology optimization method in micro-scale (nm~m level) for the design of cathode microstructures. Chapters 5 establishes a structural- 4 fluid-thermal coupled finite element model for the macro-scale (mm level) geometry of SOFC fuel supply channel and anode supported electrode, and discusses the topology optimization of the thermo-mechanical performance. Short descriptions of the chapters are below. In Chapter 2, the 2D finite element model with the consideration of the ohmic resistance in the ionic conducting scaffold and the activation polarization in the MIEC layer that surrounds the scaffold is established. The focus is on designing the shape of the scaffold and the MIEC layer to minimize the cathode resistance. Although the 2D designs are not able to be used as an example for the 3D scaffold microstructures, the 2D optimal results give a schematic direction to the cathode microstructure design by showing the effect of the material properties and geometric factors on characteristic scaffold feature sizes. In Chapter 3, the 3D finite element model is set up in a similar manner to the 2D model in chapter 2, and the 3D model is validated by experimental data. A detailed 3D model of a cathode microstructure sample is reconstructed from a sequence of images recording the cross sections of the cathode microstructure. The cathode resistance is calculated by the 3D finite element model, and the results are compared with experimental measured cathode resistance. The good consistency between the simulation and experimental results at low temperature domain makes a solid basement for the following optimization work. In Chapter 4, the topology optimization problem for the design of 3D cathode microstructures is discussed. A self-locked phenomenon caused by the use of penalized design dependent load is investigated, and this problem is solved by introducing a new material model that depends on both the first and second order spatial derivatives of the design variables. The cathode resistance of the optimal geometries is compared with the cathode resistances of the experimental sample to show the improvement. The 3D design is able to give direct and detailed guide to the cathode 5 microstructure design. In Chapter 5, an integrated finite element model which simulates the fuel flow, hydrogen gas transport, conductive-convective heat transfer, electrochemical reactions and structural deformation of the macro-scale anode-supported SOFC electrode designs is set up, and a specific topology optimization model is set up for the improvement of thermo-mechanical strength. The SOFC performances such as the current generation, flow rate, static compliance, and thermo-mechanical deformation of the electrode are calculated. A failure model is set up based on the dependency of the thermo-mechanical performance on the porosity of the porous media, and the probability of mechanical failure is used to measure the structural strength. The probability of mechanical failure is minimized by altering the geometries of the fuel flow path, modifying the path of hydrogen gas transport and heat transfer in the electrode and adjusting the structure of the anode support. Although the 3D model can accommodate several design concepts, only the design for the cross section of the anode support structure is explored in this dissertation, focusing on the concepts that are easier to manufacture. The design that optimizes the probability of failure is compared with the thermo-mechanical compliance design. 6 Chapter 2 A 2D model for shape optimization of solid oxide fuel cell cathodes An infiltrated composite cathode consists of an ionic conducting scaffold and a coating of MIEC applied on the surface of the scaffold. Finite element method has been used to analyze the performance of the cathodes. The use of finite element method also makes the efficiency compute of the sensitivity of the performance of the cathode with respect to the design, which is essential in topology optimization in chapter 2 and chapter 4. A first attention at optimizing the cathode performance in this work seeks to improve the current generation/reduce the cathode resistance by designing the microstructure of the cathode scaffold using a 2D topology optimization model. A simplified finite element model and a 2D description of the cathode geometries are used to estimate the cathode performance. The goal is to identify the structure that maximizes the current generation/minimizes the cathode resistance. The model of the SOFC is reduced to a periodic, 2D conduction problem over a rectangular representative cell with design-dependent ionic transfer boundary conditions. Special techniques are introduced to avoid physically inadmissible designs that would otherwise be allowed by the 2D model. The designs with different material properties and resource restrictions are investigated. 2.1 Simplified analysis model of infiltrated composite cathode The performance of the cathode considered here can be modeled in a manner similar to the models presented in Nicholas and Barnett [45], Shah et al. [8] and Tanner et al. [13]. The 7 approach is summarized below. The modeling of the composite cathode begins by assuming that the oxygen ionic flux density across the ionic conductor is a dilute thermodynamics process that can be described by the Nernst–Planck equation [45] - (2.1) In EQ. 2.1 the index i denotes the type of particle, Ji represents the atomic flux, per particle, i the mobility, ci the concentration, Di the diffusivity, and velocity. the charge the convective is Faraday’s constant. The first term in EQ. 2.1 corresponds to the electro-migration term which represents the flow driven by the electrical potential. The second and third terms are, respectively, diffusion and convection terms. In the case of low current densities and high porosities, all three terms can be considered independently [13]. The atomic flux in fuel cell can work in two functions: one is the fuel cell function, in which the current in driven by the diffusion of oxygen ions; the other is the electrolyer function, in which the current is driven by an applied driving potential. In the measurement of cathode resistance, the fuel cell works in the electrolyer function, and in this case, the diffusion and convection terms are zero or negligible. Thus, if one considers only one type of ionic charge carrier (such as oxygen vacancies), the effect of migration can be expressed simply as where (2.2) is the ionic conductivity in the ionic conductor. In EQ. 2.2, the current flow is proportional to the potential gradient, so the resistance of the ionic conductor is ohmic. For a steady state flow, by the continuity condition of the flow, (2.3) 8 and thus the potential distribution in the ionic conductor is governed by Laplace’s equation (2.4) The discrete electrocatalyst nano-particles deposited on the ion conducting scaffold are modeled as a continuous mixed conductor (MIEC) coating. In this coating, the resistance can be considered as a uniform surface resistance of magnitude Rs[13][45], so that the current flux density across the coating, , can be written in the form of Ohm’s law - In EQ. 2.5 the coefficient is the driving potential, (2.5) is the electric potential in the mixed conductor, and is a property of the MIEC coating. Due to the continuity of current density at the interface between the scaffold and the mixed conductor coating, the effect of the mixed conductor can be written as a convection-like (or Robin) boundary condition at the scaffold boundary by combining EQ. 2.2 and EQ. 2.5 into (2.6) In summary, in the simplified cathode model, the electric potential where the ionic conductor occupies the domain applied on the boundary of the ionic conductor the bottom of the electrolyte base is the solution to in (2.7) on (2.8) , an ionic transfer boundary condition is , and a constant potential is applied at . The formulation of this problem is similar to heat transfer problem with a convection boundary condition. 2.2 The simple model as a periodic single column reference design A reference design with a series of 1D periodic rectangular column has been used to provide a 9 good estimation of the cathode resistance (Fig. 2.1). The estimated cathode resistance is used as reference data for the comparison purpose in this chapter. The analytical solution of the cathode resistance of the reference design was proposed by Tanner, Fung and Virkar [13], and the solution was modified to handle nano-composite electrodes by Nicholas and Barnett [7]. h r electrolyte B Figure 2.1 Geometry of the reference design used for comparison The cathode resistance of the rectangular column design can be calculated by RP  rRI 1  B 1  B c( )r (1  )  (1  c)( )  r 2 2 2 2 1 c  1 c  (2.9) where (2.10) - (2.11) (2.12) is the conductivity of the ionic conducting material and RI is the surface resistance. For an 10 infiltrated cathode, Nicholas and Barnett [7] showed that RI should be revised to be RS, which is (2.13) where AS is the surface area of the ionic-conducting scaffold and AI is the surface area of the infiltrated mixed conductor particles. 2.3 Finite element model for 2D topology optimization problem For an arbitrary design with non-regular geometries, the finite element method can be used to calculate the cathode resistance. The weak form associated with the governing equation is where and and - are suitable sets of admissible functions defined of (2.14) . To complete the model, the composite cathode is assumed to be periodic in the horizontal direction with a periodicity characterized by an (irreducible) unit cell, so periodic boundary conditions are applied on the left and right boundaries, denoted by P as illustrated in Fig. 2.2(a). In the model, , the region occupied by the scaffold material, consists of two distinct components: the cathode and the electrolyte, and both are assumed to have the same conductivity . However, only the scaffold is subjected to design. 11 (a) (b) Figure 2.2 A periodic computational domain is used in the cathode model as shown in a but it is extended as in b to allow the exploration of the entire space through topology optimization As the shape of the cathode is unknown, the topology optimization problem involves design dependent boundaries, at which the boundary conditions cannot be set directly. This type of problems has been discussed by various authors in the context of heat conduction. For instance, Yoon and Kim [34] proposed an “element connectivity parameterization” method to consider design-dependent effects for thermal boundary conditions. Bruns [35] investigated a 2D problem with convection both on the top and on the boundary of the structure. Iga et al. [36] introduced a smeared-out hat function to extract the side convection. For the purpose of topology optimization, the shape of the cathode is characterized by the distribution of cathode (ionic conducting scaffold) material within the (rectangular) unit cell the computational domain (Fig. 2.2(b)). The area in or electrolyte (i.e., the set that is not occupied by scaffold material ) is occupied by gas. In an analysis suitable for topology optimization, the weak form in EQ. 2.14 is extended to the whole computational domain Glowinski et al.[46] validate this extension by showing that there exist functions and defined on , such that 12 , . (2.15) and Functions and admissible spaces are, respectively, and , function satisfying -periodic extensions of and is the characteristic function of only on properties (conductivity) in (2.16) in in suitable , and (see Wells and Zhou [47] for details). Material are characterized by , which is assumed to have the form (2.17) Ω Where is a . Note that no distinction is made between the conductivity of scaffold and electrolyte material. Typically, they have the same conductivity Discretization of the domain is piece-wise constant in . is by rectangular, four node elements. The conductivity with discontinuities only at element boundaries. Denoting the conductivity in the e-th element by resistance coefficient across the boundary and letting represent the (inverse) surface of element e, and substituting EQ. 2.15 and EQ. 2.16 into EQ. 2.14, EQ. 2.14 can be rewritten as - Ω - (2.18) Ω Noting that the second term in EQ. 2.18 can be converted to an integral in the whole analysis domain Ω That is because Ω - Ω - coincides with the unit boundary normal vector has nonzero values only in a small neighborhood of 13 (2.19) of on , and , in writing EQ. 2.18 we have assumed that coincides with element boundaries. In the same manner as the single column model in section 2.2, is assumed to be a rectangular periodic cell. (Fig. 2.3) symmetric symmetric design domain gas periodic periodic cathode root Figure 2.3 Sketch of periodic cell. One can introduce a design variable like approach [31], let : to control the material distribution in be also piece-wise constant in element e, and express material properties and .We use a SIMP with values within as - (2.20) and In EQ. 2.20 and EQ. 2.21 and - - (2.21) represent, respectively, the density and the (inverse) surface resistance coefficient across the j- th side of element e ( j = 1, 2,3 or 4) and p is the penalty factor in the SIMP model [31]. Note from (14) that material properties across the j−th side of element e and that if , e.g.,  = 0.001. 14 if there is no change in is differentiable at - 2.4 Objective function The goal of the optimization problem is to find cathode shapes that minimize resistance at the base of the cell boundary . The cathode resistance is calculated by Nicholas and Barnett [45] as (2.22) where on is the driving potential, B is the footprint (length of (Fig. 2.2), d is the thickness of the electrolyte base, and conducting material. , B, ), is the outward normal is the conductivity of the and d are all design independent, so minimizing is equivalent to maximizing the denominator of the first term in EQ. 2.22, (2.23) which is the total current density flow through . This is the objective function in this chapter. 2.5 Constraints Regardless of the shape of the cathode, the model in Section 2.2 predicts a smaller resistance for larger perimeter and larger amount of scaffold material. Larger perimeter is better simply because the current comes from the ionic transfer boundary condition. More material is better because a more massive scaffold offers less resistance to current. In order to extract the effect of shape alone, isoperimetric constraints are imposed on the perimeter and amount of material. The perimeter of the cathode can be expressed as (2.24) where is the length of the j-th side of element e. 15 The total amount of scaffold material can be measured using the element densities as (2.25) where is the element area. The perimeter and material constraints are used to control the size of the features allowed in the solution, a limitation associated with the manufacturing process and our modeling assumptions. The smallest feature is limited by the size of the particles of ionic conducting material used to build the scaffold. In addition, the assumption that the MIEC is uniformly spread over the scaffold also imposes a constraint on the smallest feature size. We use the upper bounds on the perimeter and amount of material (respectively, PR and WR) to parameterize solutions with different feature sizes: smaller bounds PR (larger bounds WR) lead — on average — to larger features. Using these bounds to control feature size is an arbitrary choice that works well in this case. A filter has a similar effect, but it typically leads to solutions with more “gray” areas. 2.6 The optimization problem Based on the objective function in Section 2.4 and the constraints in Section 2.5, the optimization problem is: Find that Maximize (2.26) Subject to W  WR P  PR and (discretized) state equations, discussed in Section 2.9 below. In this problem WR and PR are 16 prescribed upper bounds on the total amount of material and perimeter, respectively. The set is used simply to define regions 0 and 1 where the conducting material is prescribed. Specifically, (2.27) In EQ. 2.27 0 and 1 are subsets of  where is prescribed to have value 0 (gas) or 1 (scaffold or electrolyte), respectively. To break the symmetry of the periodic cell and to assign the spatial frequency of the cathode column, an electrolyte root is prescribed in the design domain (Fig 2.3). Also, to prevent the columns to connect together and form physically inadmissible design, a gas region is assigned around the design domain. 2.7 Filters A density filter with a linear weight [48] is used to prevent the appearance of checkerboard patterns and to control feature size. A smoothed Heaviside projection function [49] is used to facilitate convergence to binary solutions. The variable used in the computation of material properties in EQ. 2.24-25 are then 0. where is the filtered density and - 0. 0. 0. (2.28) is a scaling factor that determines the sharpness of the projection. 2.8 Techniques to reject inadmissible designs Because of simplifications introduced by the cathode model, the optimization problem may result in physically inadmissible SOFC designs, as explained below. As is the case with 17 checkerboard suppression schemes, special steps are required to penalize the formation of such designs. The first type of inadmissible design is a hole inside the body of the cathode without access to gas, and it is simply because it cannot be supplied with ambient gas needed for the surface electrochemical reaction. However, the cathode model may identify an internal hole as a desirable feature since the ionic transfer boundary condition in the inner hole would produce current, even when the hole is not in contact with ambient gas. Large internal holes cannot be removed by standard filters and thus a special technique is needed to prevent their formation. There are two mechanisms of forming a hole: One is by closing a “neck” region, trapping gas inside an area occupied by a scaffold, a situation illustrated in Fig. 2.4(a); the other is by reducing density in areas surrounded by the cathode material. Both of them can be identified by morphological techniques (Fig. 2.4(b)), and the internal holes can be eliminated by preventing the cathode branches from trapping a gas region and adding a negative current source to the reducing density areas surrounded by cathode material. Figure 2.4 Detecting neck regions that may trap internal holes The strength of the negative heat source is - 18 (2.29) where is a prescribed constant, (typically set to q = 50h0∞). Once an area reaches an effective density , elements in its interior are fixed to and removed from the design space by including e in 1. The artificial negative current source used in EQ. 2.29 penalizes the presence of gas (white) that is fully trapped by scaffold material (black). A negative source absorbs current and thus it reduces the flow of ions to the electrolyte base. Since less ion flow at the base means that performance is degraded, designs with negative current are less desirable. By EQ. 2.29, the negative current can only be removed by turning trapped gas into scaffold material. The second type of inadmissible design is the disjoint cathode material. To conduct current the cathode material must be singly connected and connected to gas and electrolyte. However, the cathode model may accept local optima which are not connected. To illustrate, consider the three designs in Fig. 2.5. Figure 2.5 Cathode model may predict that (b) is worse than either (a) or (c) Design (c) forms a better conduction path to the electrolyte than design (a) and should be the preferred solution. A change from (a) to (c) requires that the conductivity in the gap be increased, e.g., by increasing slightly, as shown in design (b). However, if the conductivity within the gap is increased only slightly, design (b) may actually generate less current than either (a) or 19 (c). This is because the increased current caused by higher conductivity may be canceled by a simultaneous reduction in ionic transfer through boundaries near the gap. If this were to happen, the gradient of the objective function with respect to will be positive and the gap will not close. To avoid this, a small amount of artificial surface ionic transfer with coefficient - of the order (2.30) is allowed in all elements e in the design domain. l is the element length. This term will increase ionic transfer through the gap and make design (b) better than (a). Since surface ionic transfer may influence the cathode shape, it is gradually removed as iterations progress, i.e., the magnitude of is decreased to zero gradually with iterations. This method is not the best way to resolve the single connection problem, and an improved method will be discussed in Chapter 4. 2.9 Sensitivity Analysis Upon assembly of element matrices, EQ. 2.18 becomes ( where (2.31) is the vector of values of the potential at the nodes. Matrices , , and vector are assembled from element-level quantities (2.32) (2.33) (2.34) where represents the vector of shape functions of element e. Quantities and , added to the equilibrium equations, are needed when the special techniques introduced in Section 2.8 20 are applied to the model. They represent an artificial source of heat and are assembled from element-level quantities (2.35) Material properties and (2.36) associated with element e are related to the design variables through EQ. 2.24-25, respectively. Properties the artificial heat source are related to and associated with through EQ. 2.29 and EQ. 2.30, respectively. Upon discretization, the objective function EQ. 2.26 becomes (2.37) where is assembled from - (2.38) The sensitivity of f is obtained from where (2.39) is the solution of the adjoint problem (2.40) 2.10 Examples In the following examples we investigate the influence of changes in material properties, allowable perimeter, and the amount of cathode material used. Since infiltration typically produces MIEC nano-particles 50–70 nm in diameter, the smallest feature size that can modeled before the assumption that the MIEC is uniformly spread over the scaffold is violated may be of the order of 300 nm, so the element size is selected in the same order of the smallest feature size. In the examples all length dimensions are expressed as multiples of element length. Except as noted, square elements with l=1 are used. Material properties and dimensions can be scaled by 21 the ratio . Except as noted, we set h0/0 = 2.510 . This results in material properties -4 consistent with values reported in Nicholas and Barnett [7] and Shah et al. [8]. Except as noted, in all examples the height of the electrolyte base is d = 6 and the fixed gas layer above the electrolyte and on the cell sides (see Fig. 2.3) is 2. The cathode root (i.e., the area of fixed density that connects the scaffold material to the electrolyte) is × . Note that in scaling the problem one needs to be mindful of the limitations of the analysis model. Since infiltration typically produces MIEC nano-particles 50–70 nm in diameter, the smallest feature size that can modeled before the assumption that the MIEC is uniformly spread over the scaffold is violated may be of the order of 300 nm. Finally, all designs are forced to be left-right symmetric. All optimization problems are solved using the method of moving asymptotes of Svanberg [ 50 ]. Results reported are after post-processing, applying a threshold at to remove a very small number of intermediate density elements left in the optimal solutions. 2.10.1 Designs with different material properties In this example the unit cell is 60 × 40(l). The constraint bounds are PR  360 and WR  0.25W where W is the area of . Figure 2.6 and Table 2.1 illustrate the effect of varying the ratio on the optimal solution. Larger ratios are associated with higher performance electrocatalyst material. In the solutions shown in Fig. 2.6a, b and c, is smaller, surface resistance dominates scaffold resistance and therefore the optimized designs tend to exhibit larger perimeter and more and longer branches. The ratio (b) used in the solution shown in Fig. 2.6d is larger, the electrocatalyst material is more effective, and focus shifts from the perimeter to the body of the scaffold. 22 Performance is improved in this range by emphasizing the path of the ions through the scaffold rather than by increasing the perimeter. The increased flux through the boundary that results from the improved coating material calls for a well-designed path to the electrolyte at the bottom to improve performance. h0/k0=1.25×10-4-4 h0/σ0=1.25×10 (a) h /k =2.5×10-4-4 h00/σ00=2.5×10 (b) h0/k0=1.25×10-3-3 h0/σ0=1.25×10 (c) h0/k0=2.5×10-3-3 h0/σ0=2.5×10 (d) Figure 2.6 Optimized designs for different values of h0/ 0. As h0/ 0 increases, the perimeter becomes less important and the shape of the scaffold begins to dominate. Performances reported in Table 2.1 Table 2.1 Performance of optimized designs for different materials, where R can be multiplied by l/ 0 to recover units (Fig. 2.6) 4 −3 −3 Design h0/σ0×10 (a) 1.25/l 1.33 1.80 26 (b) 2.5/l 0.72 1.07 33 (c) 12.5/l 0.17 0.35 51 (d) 25/l 0.097 0.21 54 RP×10 23 RREF×10 Improvement (%) 2.10.2 Design with varying footprint B and constant perimeter/footprint ratio In this example the unit cell has a height of 60(l) and varying footprint B. The bound on the amount of material is WR  0.25W . Fig. 2.7 and Table 2.2 illustrate the effect of varying B while keeping the ratio PB /B = 9 on the optimal solution. B=20 B=20 B=40 B=40 B=60 B=60 (a) (b) (c) Figure 2.7 Optimized designs for different footprints B and PB/B = 9. For a large footprint, the cathode shape becomes more convoluted and departs further from a columnar shape, not being confined by the design space. Performances reported in Table 2.2 Table 2.2 Performance of optimized designs for different footprints B and PB/B = 9, where R can be multiplied by l/ 0 to recover units (Fig. 2.7) −3 −3 RP×10 RREF×10 Design B/l Improvement (%) (a) 20 0.46 0.56 18 (b) 40 0.49 0.80 39 (c) 60 0.54 1.05 49 As shown in Table 2.2, RP remains nearly constant in all optimized designs obtained while varying B while keeping the ratio PB/B constant. However, the optimized design performs better than a reference design of the same perimeter and footprint. This improvement relative to the reference design is more pronounced as the footprint increases, indicating that cathode shape is 24 more important for larger footprints. 2.10.3 Design with varying perimeter P and amount of material W while keeping the footprint constant In this example the unit cell is 60×40(l). Figure 2.8 and Table 2.3 illustrate the effect of varying PR ,with WR  cW where c=0.167, 0.25, 0.33, for PR/l equal to 260, 360, and 460, respectively. (a) PR=260 (b) PR=360 (c) PR=460 Figure 2.8 Optimized designs for different values of the perimeter bound PR. As expected, a large perimeter corresponds to a greater performance improvement. Performances reported in Table 2.3 Table 2.3 Performance of optimized designs for different values of the perimeter bound PR, where R can be multiplied by l/ 0 to recover units (Fig. 2.9) −3 −3 RP×10 RREF×10 Design PR/l Improvement (%) (a) 260 0.94 1.19 21 (b) 360 0.72 1.07 33 (c) 460 0.57 1.04 45 As shown in Table 2.3, RP decreases with increasing perimeter. The percentage of improvement compared with the reference design also increases with perimeter. 2.10.4 Design with varying amount of material In this example the unit cell is 40 ×60(l). The cathode root has a width of 14. The upper bound 25 on the perimeter is PR = 360. Fig. 2.9 and Table 2.4 illustrate the effect of varying the bound on -4 the amount of material, WR. Results show that for the material used, h0/0 = 2.510 , this is a perimeter-dominated problem. Since the perimeter is the same in all solutions, the resistance is also the same. Increasing the amount of material does not affect the answer. W =0.167W R (a) Ω W =0.25W R Ω W =0.333W R (b) Ω (c) Figure 2.9 Optimized designs for different values of the amount of material bound WR. Performances reported in Table 2.4 Table 2.4 Performance of optimized designs for different values of the amount of material bound WR, where R can be multiplied by l/ 0 to recover units (Fig. 2.9) Design WR/W −3 RP×10 −3 RREF×10 Improvement (%) (a) 0.167 0.49 0.91 46 (b) 0.25 0.49 0.80 39 (c) 0.333 0.49 0.73 33 2.11 Conclusions The designs found in this chapter were very intricate and obviously difficult to realize. However, the results are useful in that they provide insights on the characteristics of optimal microstructures and guidance on what features to emphasize in cathode design. Improvements of nearly 50 % were found, compared to a simple column design. The perimeter was found to be of the utmost importance in optimizing the microstructure. However, as the ionic transfer efficacy 26 of the MIEC particle coating improves, the shape of the scaffold also becomes important. The studies point to the importance of the microstructure in the performance of SOFCs and suggest that significant improvements are possible by proper organization of the microstructure. In addition, the 2D optimization model was set up based on the methodology that uses 2D geometries to evaluate the performance of 3D microstructures. However, the complicated connections between different parts of the 3D scaffold microstructure cannot be simulated by the 2D model, and the 2D designs are also not able to improve these connections. A 3D simulation model and a 3D topology optimization model are proposed in Chapters 3-4 to resolve these problems. 27 Chapter 3 A 3D finite element model for cathode microstructure of solid oxide fuel cell cathodes In this chapter, a 3D finite element formulation is developed to overcome the limitations in the 2D simple model introduced in section 2.2. The evaluation of the cathode resistance of a porous cathode microstructure is discussed, with an emphasis on the effect of the surface resistance of the MIEC coating, the bulk resistance of the ionic conducting scaffold and the gas diffusion in the pores. To validate the 3D finite element model, a detailed mesh is reconstructed from the experimental data from 3D Focused Ion Beam-Scanning Electron Microscopy, and the simulation results with different sets of material properties are compared with experimentally measured results. 3.1 Analysis model of infiltrated composite cathode In the cathode microstructure, the oxygen gas is supplied on the top of the cathode, and it flows through the holes inside the porous cathode to the mixed conductor at the surface of the cathode scaffold. Electrochemical reactions, which consume oxygen gas and generate current, happen in the mixed conductor. The current flows through the cathode scaffold to the electrolyte at the bottom of the cathode (Fig. 3.1). 28 Γtop Air flow path Ω Γω Cathode Scaffold (ω) z x y Γ0 n0 Figure 3.1 The sketch of a 3D cathode microstructure An analysis model is performed below to simulate the current generation in the MIEC layer and the current flow through the ionic conducting scaffold to the bottom. This simulation model is to be used as a preparation of the 3D optimization model in the next chapter. The analysis domain  is selected to be the region occupied by the solid cathode microstructure and the gas pores. The solid cathode region The current flow in is a sub-domain of , and corresponds to the fluid region. is modeled in a similar manner with the 2D case, in which the current is described by the electrical potential distribution, and it is governed by Laplace’s equation, as in chapter 2 in In. EQ. 3.1 is the ionic potential and (3.1) is the conductivity of the ionic conducting material. 29 The oxygen gas flow is diffusion dominated, so the concentration distribution also satisfies Laplace’s equation in where is the oxygen gas concentration and At the gas-cathode interface (3.2) is the diffusivity of the oxygen gas. , the electrochemical reaction absorbs oxygen gas from the air and releases oxygen ion to the solid cathode, so the current generated by the electrochemical reaction is modeled using a convection-like oxygen exchange boundary condition. The dependency of the ionic potential distribution and the current generation density obeys Ohm’s law at low current density [13], - where is the current generation per unit area of from the Nernst voltage, and on , (3.3) is the driving potential which comes is the interface resistance. The previous 2D model assumes that the oxygen gas supply is abundant, so the interface model is treated independent of the oxygen gas concentration. With the introduction of oxygen gas consumption, the interface resistance also depends on the concentration of oxygen gas by the Butler-Volmer equation [51] RS     i0 0. where is the cathodic pre-exponential coefficients, c ref is the reference concentration, (3.4) - (3.5) is the concentration of oxygen gas, is the cathodic activation energy, R is the universal gas constant and T is the temperature. All of these quantities are assumed to be constants that are 30 independent of the designs. EQ. 3.5 can be linearized for the convenience of calculations, and EQ. 3.3-3.5 can be simplified as - on (3.6) - where is an approximation of the oxygen exchange coefficient. The corresponding oxygen gas loss is - on (3.7) where z is the particle charge of oxygen ion, qe is the electron charge and n0 is the Avogadro constant. 3.2 Finite element formulation of 3D cathode microstructure For a specified cathode microstructure, a finite element method can be used to calculate the cathode resistance. Weak forms associated with the governing equations are - - (3.8) and - - (3.9) where and and and and are suitable sets of admissible functions defined of are suitable sets of admissible functions defined of 31 , and . A constant potential is applied at the bottom concentration of oxygen gas at the top of the analysis domain , and the is assigned to be the same as that in the standard atmosphere of the analysis domain , as illustrated in Fig. 3.1. In preparation for the optimization work in chapter 4, both of the sub domains for current flow and gas diffusion are extended to the whole analysis domain , in the same way as in chapter 2. Also, the analysis domain is discretized uniformly by identically shaped rectangular elements, and only the element boundaries can be considered as the boundary of the cathode - - (3.10) (3.11) , Functions , , suitable admissible spaces are, respectively, extensions of , , and , , and in . Noting that (3.12) (3.13) where and . The cathode resistance is calculated by measuring the current flow through the bottom boundary of the cathode (3.14) where is the driving potential, A is the footprint (area of (Fig. 3.1). 32 ), is the outward normal on 3.3 Generation of 3D microstructure mesh from experimental data The micro-scale finite element model of the infiltrated composite cathode has been set up in section 3.2. Before applying the finite element model to the optimization work in chapter 4, the simulation result of the 3D finite element model should be validated by the experimental results. The most convenient way for the validation is to construct a 3D mesh of the 3D microstructure used in the experiment and compare the simulation results with the experimental results, and these require the reconstruction of the 3D microstructure of the cathode. First, the manufacturing process of the composite cathode is reviewed [14]. The cathode scaffold is made by (a) compressing and firing the particles of the ionic conducting material (CGO, etc.) into a pellet electrolyte base and (b) screen printing the ionic conducting particles onto the pellet and firing to make the particles merge into a porous scaffold. After that, (c) the electro-catalytic nano-particles (LSCF) are infiltrated into the scaffold, and they are fired to form a coating on the surface of the scaffold. The final step (d) is to print a current collector layer on the top, and oxygen gas can go through this layer to the cathode. In the experimental environment, one of the most common methods of reducing the surface resistance of the MIEC layer is to use smaller electro-catalytic nano-particles. The figures used below are from the work of Nicholas et al.[14] . (a) Produce 98%+ Dense Ce0.9Gd0.1O1.95 (CGO) Pellets by Firing at 1450° C for 6 hrs Figure 3.2 The manufacture process of a SOFC composite cathode sample [14] 33 Figure 3.2(cont’d) (b) Screen Print ~15mm Thick CGO Scaffold Layers, Fire at 1100 ° C for 1 hr (c) Infiltrate Nitrate Solutions into CGO Cathode, Gel Solutions at 80 ° C, Fire at 800° C for 1 hr (d) Screen Print The Current Collector Layer, Fire at 800° C for 1hr. The randomness of the particle arrangement in (b) makes it difficult to give a good and detailed prediction to the resulting structure. Fortunately, the 3D microstructure can be reconstructed based on the experimental photos of the micro-scale scaffold cross sections. 2-D images of the electrode microstructures have been obtained from Dr. Scott Barnett’s research group at Northwestern University [52][53]. The experimental photos can be made by (a) taking an image of the present surface of a pellet by 3D Focused Ion Beam-Scanning Electron Microscopy; (b) milling a slice of the pellet off from the side, and taking the image of the newly created surface; (c) stacking the surfaces together and interpolating the solid region between the 34 surfaces. The whole process is shown schematically in Fig. 3.3. The figure and data used in the Γtop Ω reconstruction courtesy of Prof. Nicholas’s group in Michigan State University [54]. Taking image of the new cross section Taking image of the original cross section Remove one layer Porous cathode microstructure (a) (b) Stacking the cross sections to reconstruct the structure (c) Figure 3.3 The reconstruction process using 3D Focused Ion Beam-Scanning Electron Microscopy. 35 Γω Γtop Ω 9μm z y Figure 3.4 An example of a black-white image of a 9μm thick scaffold. Image courtesy of Prof. Nicholas’s group in Michigan State University [54] Fig. 3.4 shows an example of a black-white image of a 9 microns thick scaffold. The thickness of each slice is 24.8nm (x direction). Each pixel in the photo has a size of 15.63nm×15.63nm (y and z directions). The black pixels represent the ionic conducting particles, and the white pixels represent the void domain. Because the sizes of the particles are around 200-300nm, which are much larger than the resolution level and the slice thickness, the images provide an accurate description of the microstructure. The image pixels can be mapped into a 3D finite element mesh to reconstruct the microstructure. For the ease of the mapping, and also to make preparations for the following optimization in chapter 4, a uniform rectangular meshing is selected. Two types of elements are used: one is for the ionic conducting material; the other is for the void domain, which corresponding to black and white pixels in the photos as shown in Fig. 3.4, respectively. Possibly easiest way of the mapping the black-white pixels into the two types of elements is a one to one mapping, in which each element corresponds to be one pixel in the photos, and the size of the element is determined by the size of the pixel and the thickness of each slice. 36 However, a one to one mapping will result in a very fine mesh, and the cost of the simulation can be prohibitive. To reduce the cost of the simulation, noting that the size of each pixel is much smaller than the size of the ionic conducting particle, pixels can be grouped, and each group is mapped into one element. In the following simulation, the sizes of the elements are selected to be 49.6nm×46.9nm×46.9nm, which corresponds to a pixel group of 2×3×3 (in x, y, z directions, respectively). In each pixel group, the black pixels have a grey scale 0, and the white pixels have a grey scale 255. If the average grey scale in one group is smaller than half of 255, the corresponding element is considered to be black, and otherwise it is assigned to be white. The resulting 3D mesh of the 9 microns’ thick cathode microstructure is shown in Fig. 3.5. The ‘black’ elements show the geometry of the scaffold, and the interfaces between the ‘black’ and the ‘white’ elements are the surfaces of the scaffold. 37 μm 9 8 7 6 5 n0 4 3 2 1 z μm 0 0 0.5 μm 0 1 0.5 x y Figure 3.5 The 3D mesh of the 9 micron thick cathode microstructure 3.4 Examples of the simulation of the 3D cathode resistance In the following example, the cathode resistance of the scaffold is estimated by using the finite element method introduced in section 3.2 and the 3D reconstructed microstructure shown in Fig. 3.5. The ionic conducting material is selected to be CGO, and the material for MIEC is selected to be LSCF. The use of smaller LSCF particles can reduce the cathode resistance. Two particle sizes are 38 selected for the simulation. The surface resistance of the mixed conductor (LSCF) and the conductivity of the ionic conducting material (CGO) are shown in Tables 3.1-3.2 [54]-[56]. Table 3.1 Bulk conductivity of CGO (Gd0.1Ce0.9O2) Temperature (°C) Temperature (1000/T(K)) Conductivity (S/cm) 400 1.485552997 1.76×10 450 1.382838968 3.69×10 500 1.293410076 6.90×10 550 1.214845411 1.14×10 600 1.145278589 1.77×10 650 1.083247576 2.67×10 700 1.027590813 3.86×10 -3 -3 -3 -2 -2 -2 -2 Table 3.2 Scaled ASR values for 20 nm MIEC LSCF (La0.6Sr0.4Co0.8Fe0.2O3) particle diameter 2 Temperature (°C) Temperature (1000/T(K)) ASR (Ω×cm ) 400 1.485552997 229.651 450 1.382838968 70.044 500 1.293410076 25.262 550 1.214845411 8.727 600 1.145278589 4.134 650 1.083247576 1.837 700 1.027590813 0.919 39 Table 3.3 Diffusivity and concentration of oxygen gas in air at 1atm [57] Temperature (°C) Temperature (1000/T(K)) σg (m /sec) 400 1.485552997 7.65×10 450 1.382838968 8.52×10 500 1.293410076 9.42×10 550 1.214845411 1.03×10 600 1.145278589 1.13×10 650 1.083247576 1.23×10 700 1.027590813 1.33×10 2 3 cmax (mol/m ) -5 3.58 -5 3.33 -5 3.11 -4 2.92 -4 2.75 -4 2.60 -4 2.47 The 9μm cathode sample reconstructed in section 3.3 is used. The 3D mesh has 28×20× 193=108080 elements. The porosity denoted by vol% of the sample is 49.15%. The oxygen gas concentration is prescribed to be cmax at the top of the analysis domain Γtop. Because the maximum values of σg and cmax shown in Table 3.3 are at most twice of the minimum values, and the change is much smaller than those of the conductivity and ASR shown in Tables 3.1-3.2 which can be as much as 20~200 times, the change of σg and cmax versus temperature is neglected. The values of σg and cmax at 600°C are selected for all the temperatures in the following section except as noted. 3.5 Experimental validation and comparison with a simple estimation model From a plot of the cathode resistance versus operating temperature, one can observe a significant increase of the cathode resistance as the temperature decreases. This reflects the main difficulty in the development of low temperature SOFC. The results are compared to the experimental results measured by the group of Prof. Nicholas [54], and the comparison is shown in Fig. 3.6. At 40 the low temperatures, a difference of less than 5% between the experimental data and the results of 3D model is observed. The calculated cathode resistance becomes significantly lower than the experimental result near the high temperature end, where the difference can be up to 50%. 10 RP (Ωcm2) 1 3D model -4 2 =1.13×10 m /sec 0.1 Experiment 0.01 1 1.1 1.2 1.3 1.4 1.5 1000/T (1/K) Figure 3.6 The comparison of RP data for 20nm LSCF particles. Experimental data courtesy of Prof. Nicholas’s group in Michigan State University [54] The difference in the high temperature region shows that there can be some physical effect ignored in the model that cannot be neglected at high temperature. As the large difference only happens at high temperature (high oxygen gas consumption), one possible reason is that the supply of oxygen gas is not abundant for the electrochemical reactions to take place, and the diffusion of oxygen gas in experiment is weaker than it is assumed by the 3D model. To investigate the effect of the oxygen gas diffusion, we calculate the cathode resistance using a low diffusivity -6 2 =5×10 m /sec, which is about 5% of the data 41 -4 2 =1.13×10 m /sec used in Fig. 3.6. As shown in Fig. 3.7, using the reduced diffusivity value, the maximum difference between the simulated cathode resistance and the experimental data is less than 12%, which is much smaller compared to the 50% maximum in Fig. 3.6. This suggests that the diffusivity used in the 3D model should be smaller than the values shown in Table 3.3. That is because the gas diffusion in nano-scale pores works differently from the diffusion in a large room. The diffusivity shown in Table 3.3 only works in the case that the size of the pores is much larger than the mean free path of the motion of the gas molecules. However, in the cathode microstructure, the size of the pores is comparable to the mean free path, which normally has a dimension of 100nm, so the pores block the gas diffusion and make the effective diffusivity smaller. 10 RP (Ωcm2) 1 3D model -6 2 =5×10 m /sec 0.1 Experiment 0.01 1 1.1 1.2 1.3 1000/T (1/K) 1.4 1.5 Figure 3.7 The comparison of RP data for 20nm LSCF particles with a small oxygen gas diffusivity. Experimental data courtesy of Prof. Nicholas’s group in Michigan State University [54] 42 3.6 Other examples 2 Here a simulation is performed using a very high oxygen gas diffusivity 10000 times larger than the data =1m /sec, about -4 2 =1.13×10 m /sec used in Fig. 3.6. By using a very high diffusivity, the gas supply becomes abundant, and the effect of oxygen gas diffusion is removed. Results are shown in Fig. 3.8. The cathode resistances with high diffusivity ( smaller than the cathode resistances for 2 =1m /sec) are -4 2 =1.13×10 m /sec. However, the small gap between the two series in Fig. 3.8 shows the effect of oxygen gas diffusion is not strong for diffusivity values in the range shown in Table 3.3. 10 RP (Ωcm2) 1 3D model -4 2 =1.13×10 m /sec 0.1 3D model 2 =1m /sec 0.01 1 1.1 1.2 1.3 1000/T (1/K) 1.4 1.5 Figure 3.8 The comparison of RP data for 20nm LSCF particles with a very large oxygen gas diffusivity The last example in this section compares the simulation results of the 3D model to the 2D 43 simple reference model in chapter 2. The comparison is shown in Fig. 3.9. The cathode -4 2 resistances are calculated using the diffusivity =1.13×10 m /sec. The two series of cathode resistance are close to each other, and at most temperatures the resistance of the 3D model is larger than the resistance predicted by the 2D reference model. The small differences between the two series of results happen because the scaffold material we use has a large conductivity and the effect of scaffold resistance is not very remarkable. For other scaffold materials with lower conductivity, the difference between the cathode resistance calculated by the 3D model and the cathode resistance from the 2D model will be larger, and one can expect that the cathode resistance calculated by the 3D model will be much higher than the cathode resistance from the 2D reference model. 10 RP (Ωcm2) 1 3D model 2D reference model 0.1 0.01 1 1.1 1.2 1.3 1000/T (1/K) 1.4 1.5 Figure 3.9 The comparison of RP data obtained from different methods for 20nm LSCF particles 44 3.7 Conclusions The cathode resistance calculated by the 3D finite element model decreases rapidly as the temperature increases, but higher temperatures are difficult to accommodate in industrial applications. The comparison between the simulation results and the experimentally measured results shows a small difference (<5%) at low temperature region, but the difference becomes large (50%) at high temperatures. The large difference comes from the reduction of oxygen gas diffusivity flowing through nano-scale pores. The difference can be reduced using smaller oxygen gas diffusivity in the 3D finite element model. The accuracy of the 3D finite element model at low temperature is higher than the simple model used in chapter 2, so the 3D finite element model provides a better foundation for the topology optimization model of the microstructure in the next chapter. 45 Chapter 4 Topology optimization for 3D cathode microstructure of solid oxide fuel cell cathodes The optimal designs of 2D cathode microstructures obtained in chapter 2 qualitatively show the effect of the material properties and the geometric factors on the optimal feature and the performance. However, the 3D microstructure can have more complicated 3D connections than between different cathode scaffold parts than the 2D model and this can change the bulk resistance of the scaffold, and also the effect of oxygen gas consumption is not considered. In this chapter, a 3D topology optimization model is set up to design the 3D microstructure of the SOFC cathode. The simulation is performed based on the finite element model introduced in chapter 3. The optimization problem is setup in a similar manner as the design problem of 2D cathode microstructure in chapter 2. The ionic conducting material is distributed in a 3D design domain to form a 3D cathode scaffold, and the MIEC is treated as a design-dependent boundary condition of the scaffold. The topology of the scaffold is optimized to minimize the cathode resistance/maximize the current generation. 4.1 Problem Setup The 3D finite element model of the SOFC cathode has been established and validated in chapter 3. In preparation for topology optimization, the analysis sub domains for both the current in cathode and the gas diffusion in the pores were extended to the whole analysis domain 46 , and was discretized by uniform rectangular elements. The boundary of the scaffold was assigned to be along the element boundaries. By assigning a design variable to each element within the analysis domain of ionic conducting material can be represented by modifying , the distribution , and the design can be changed by . The 3D design domain is selected to be a sub domain of the analysis domain 4.1. The analysis domain , as shown in Fig. is a representative cell in a SOFC cathode microstructure. The height of the representative cell is the same as the cathode thickness ( z- direction), and the cathode is composed by a series of periodic cells in x- and y- directions. As in the 2D case in section 2.6, a fixed region is assigned at the bottom of the analysis domain, and in this region the geometry cathode and gas sub-domains are pre-assigned, i.e. the geometry is fixed artificially to show the effect of the cathode root patterns from which the cathode microstructure grows. In contrast with the 2D case, the cathode scaffold is permitted to touch the side boundaries of the analysis domain, because the 3D structure and the gas pores can have off-plane connections. The optimization problem is setup as follows. The objective of the problem is to find a 3D distribution of the ionic conducting material described by the design variable of the elements with the design domain shown in Fig. 4.1, to minimize the cathode resistance (maximize current generation). The current flow is measured at the bottom boundary of the design domain. The detailed finite element formulation is described in chapter 3. 47 Γtop Ω Design Domain 9.05μm z y x 0.0188μm Θ (Fixed) n0 Γ0 0.694μm 0.469μm Figure 4.1 The design domain of the 3D microstructure problem 4.2 Material model As introduced in chapter 3, the discretized forms of the governing equations are  (4.1)  (4.2) where is assigned to be along the element boundaries. 48 In the same manner as in chapter 2, a SIMP like approach is used to extend the binary material properties to a continuous function of the design variable, let constant in with values : be also piece-wise within element e, and express the conductivity as - (4.3) - (4.4) where p is the penalty factor in the SIMP model. In addition, a design dependent (inverse) surface resistance coefficient replace is introduced to to show the change in design.  (4.5)  (4.6) In the method introduced in chapter 2, a penalized model can be used to describe the (inverse) surface resistance coefficient - where and - - (4.7) represent, respectively, the density and the (inverse) surface resistance coefficient across the j-th face of element e ( j = 1, 2, 3, 4, 5 , 6), as illustrated in Fig. 4.2. if there is no change in material properties across the j-th side of element e. 49 6 3 2 4 1 5 z y x Figure 4.2 The arrangement of the element face number Note from EQ. 4.7, the penalty is used to resolve the non-differentiability problem of - at - . By using the penalty, if , e.g., grow = 0 e = 1 (a) becomes differentiable and at . grow = 0.2 grow = 1 e = 1 e = 1 (b) (c) Figure 4.3 Penalty in the surface resistance model may make design (b) performs worse than (a) and (c) Further investigation shows that the use of EQ. 4.7 may accept local optima and prevent the design from improving, which is similar with the local minimum problem discussed in section 2.8. To illustrate, consider the three designs in Fig. 4.3. Design (c) has a larger surface area than design (a) and should be the preferred solution. A change from (a) to (c) requires increasing 50 slightly, as shown in design (b). However, if within the gap is increased only slightly, design (b) may actually generate less current than either (a) or (c). This is because the increase in the surface area of the scaffold may be canceled by a simultaneous reduction the (inverse) surface resistance coefficient. If this were to happen, the gradient of the objective function with respect to will be positive and the design will never grow. To avoid this, a modification of the material model base on both of the first and the second order derivative of - is employed. The idea is to use the value of the second order derivative to determine the trend of the change in - the first order derivative. The second derivative of - - - Even if - is defined as - - - (4.8) - - - (4.9) - - - (4.10) , when the corresponding second derivatives are not zero, as shown in Fig. 4.4, if the second derivative has a non-zero value, the trend of change of - can be determined. In this case, the penalty can be removed and design (b) can perform better than (a). To avoid the influence of small perturbation of the second derivative, the penalty in EQ. 4.7 is only removed when the absolute value of the corresponding second derivative is greater that 0.1. 51 Figure 4.4 The second derivative can determine the trend of the first derivative 4.3 Constraints In the SOFC experimental work, porosity is one of the most common parameter that helps to describe the properties of the porous cathode microstructure. A smaller porosity means more ionic conducting material in the microstructure, and that usually leads to a decrease in scaffold resistance. To make a fair comparison that shows that the improvement only comes from the design of the geometry pattern, a minimum porosity is assigned to the design. The porosity vol% can be computed by the density distribution. (4.11) where (4.12) (4.13) The minimum porosity is selected to be 49.15%, which is the same as the porosity of the 9 microns thick cathode example. 52 4.4 Optimization Problem Using the cathode resistance calculated by EQ. 3.14 as the objective function and the constraints in EQ. 4.11, the optimization problem is: Find that Minimize (4.14) Subject to (4.15) and (discretized) state equations, discussed in section 4.5 below. In this problem prescribed lower bounds on the porosity. The set is is used simply to define regions 0 and 1 where the conducting material is prescribed. Specifically, (4.16) In EQ. 4.16 0 and 1 are subsets of where is prescribed to have value 0 (gas) or 1 (scaffold or electrolyte), respectively. In the manner of section 2.7, a density filter with a linear weight is used to prevent the appearance of checkerboard patterns and to control feature size. A smoothed Heaviside projection function is used to facilitate convergence to binary solutions. 4.5 Sensitivity Analysis Upon assembly of element matrices, EQ. 4.5 and EQ. 4.6 becomes 53 ( (4.17) (4.18) where is the vector of values of the potential oxygen gas concentration at the nodes. Matrices , and , is the vector of values of the , , and vector , are assembled from element-level quantities (4.19) (4.20) (4.21) - (4.22) (4.23) - (4.24) where Ne represents the vector of shape functions of element e. Upon discretization, the objective function EQ. 4.14 becomes (4.25) where is assembled from - (4.26) The sensitivity of RP is obtained from (4.27) 54 where is the solution of the adjoint problem (4.28) 4.6 Numerical Examples In the following examples we investigate the influence of changes in material properties and the arrangement of cathode roots. In the examples, the element size in the finite element mesh in the same as the element size extracted from experimental data in section 3.3, and the height of the whole mesh is also the same as the example corresponding to 20nm diameter LSCF particles in section 3.3. To reduce the computational cost, the length and width of the analysis domain are both half of the analysis domain used in chapter 3, as shown in Fig. 4.1. Since infiltration typically produces MIEC nanoparticles 20–50 nm in diameter, the smallest feature size that can modeled before the assumption that the MIEC is uniformly spread over the scaffold is violated may be of the order of 200 nm. As introduced in section 4.1, the design domain is set up upon some fixed root. Three arrangements of the roots are selected as shown in Fig. 4.5. The arrangements have five, four and one roots, respectively, but the total footprints (area) are the same. The foot prints of the arrangements all are 50% of the bottom boundary reconstructed microstructure in chapter 3. 55 , and the ratio is consistent with the 3D (a) Five roots (b) Four roots (c) Single root Figure 4.5 The arrangements of the fixed roots All optimization problems are solved using the method of moving asymptotes of Svanberg [50]. Different from the 2D designs in Chapter 2, the 3D design is not forced to be symmetric to permit more complicated designs to be generated. Results reported are after post-processing, applying a threshold at 0.5 to remove the intermediate density elements left in the optimal solutions. To save time and resources, the width and length of the design domain are selected to be half of the size of the analysis domain used in Chapter 3. To make sure the design domain is large enough to be used as a periodic cell in the cathode microstructure, the cathode resistance at are calculated using the whole reconstructed microstructure in Sec. 3.3 and a quarter of the microstructure with half the width and length (shown in Fig. 4.6). The cathode resistances are 2 2 3.17Ω×cm and 3.23Ω×cm , which are very close. The good consistency validates the selection of the size of the design domain. 56 μm 9 8 7 6 5 n0 4 3 2 1 z μm 0 0 0.5 μm 0 1 0.5 x y Figure 4.6 The comparison between the analysis domain in Chapter 3 and the design domain in Chapter 4 (in white rectangle) 4.6.1 Designs with varying the material properties by changing work temperature In this example the porosity constraint bounds is set to 49.15%. The root arrangement in Fig. 4.5(a) is used in this example. Fig. 4.7 and Table 4.1 illustrate the effect of varying work temperature on the optimal solution. The material properties 57 and RS are determined by the work temperature, and higher temperature is associated with higher performance ionic conducting and electrocatalyst materials. The oxygen gas diffusivity -4 2 =1.13×10 m /sec at is used for all the design examples. Table 4.1 Performance of optimized design at different temperatures (Fig. 4.7-4.8) Bulk Temp Conductivity (°C) σe (S/cm) -3 400 1.76×10 450 3.69×10 500 6.90×10 550 1.14×10 600 1.77×10 650 2.67×10 700 3.86×10 -3 -3 -2 -2 -2 -2 2 RS(Ω×cm ) 1/ (RS×σe) Optimized RP Reference RP (Ω×cm ) (Ω×cm ) 2 2 Improvement (%) 229.651 0.99 2.51 3.16 21 70.044 1.55 0.836 1.10 24 25.262 2.30 0.336 0.457 27 8.727 4.02 0.14 0.199 30 4.134 5.47 0.0733 0.1076 32 1.837 8.16 0.0385 0.0583 34 0.919 11.3 0.0222 0.034 35 58 μm 9 μm 9 μm 9 8 8 8 7 7 6 6 5 5 4 4 4 3 3 3 2 2 2 1 1 1 0 0 0 0 0.4 0.5 0 μm (a) 400°C 7 6 5 0 0.4 (b) 450°C 0.5 0 μm 0 0.4 0.5 0 μm (c) 500°C Figure 4.7 The optimized design at different work temperatures. Performance is reported in Table 5.1 59 Figure 4.7 (cont’d) μm μm 9 μm μm μm 9 9 9 8 8 8 8 7 7 7 7 6 6 6 6 5 5 5 5 4 4 4 4 3 3 3 3 2 2 2 2 1 1 1 1 0 0 0.5 0.4 0 (d) 550°C μm 0 0 0.5 0.4 0 μm 0 0 (e) 600°C 0.5 0.4 0 μm 0 (f) 650°C 60 0 0.4 0.5 0 μm (g) 700°C 10 Optimal Design 1 RP (Ωcm2) Experiment Critical cathode resistance for industrial application 0.1 0.01 1 1.1 1.2 1.3 1000/T (1/K) 1.4 1.5 Figure 4.8 The comparison of RP data for 20nm LSCF particles with a small oxygen gas diffusivity. Experimental data courtesy of Prof. Nicholas’s group in Michigan State University [54] As shown in Fig. 4.8 and Table 4.1, the optimal results have remarkable improvement compared with experimental data. Because the large difference between the simulation model and experimental data (shown in Chapter 3) at high temperature region, the optimal designs in high temperature range is not reliable, but the designs make also improvement in low temperature 2 range. The lower optimal cathode resistances are closer to the critical value 0.1Ωcm for industrial application. The ratio 1/(RSe) [55] shown in Table 4.1 has a similar effect with the ration used in chapter 2, and the increase of 1/(RSe) with the increase of the work temperature shows the 61 improvement of the MIEC performance outweighs the improvement of the bulk conductivity. In the solutions shown in Fig. 4.7, (a)-(c), the ratio 1/(RSe) is smaller, surface resistance dominates scaffold resistance and therefore the optimized designs tend to have thinner current paths through the scaffold (shown in red circles in Fig. 4.7 (a)-(c)) and span the whole design domain to form smaller features and larger surface areas (shown in blue circles in Fig. 4.6 (a)-(c)). As the ratio 1/(RSe) increasing from 0.99 to 11.3, the electrocatalyst material is more effective, and the focus of the solutions shown in Fig. 4.7 (d)-(g) gradually shifts from the surface area to the body of the scaffold. The path of the current through the scaffold rather becomes more significant than the surface reactions, and the scaffold features become thicker to form a better current path (shown in red circles in Fig. 4.7 (d)-(g)). The increased flux through the boundary that results from the better material performances at high temperature calls for a well-designed path to the electrolyte at the bottom, and to achieve this the heights of the column in the design shown in Fig. 4.7 (d)-(g) become shorter than the height of the design domain to make thick columns from the same amount of material (shown in blue circles in Fig. 4.7 (d)-(g)). Table 4.2 Tortuosity and surface areas of optimized design at different temperatures (Fig. 4.7) Temp (°C) 2 Surface Area (cm ) -7 400 3.68×10 450 3.61×10 500 3.51×10 550 3.30×10 600 3.24×10 650 700 -7 -7 -7 -7 -7 3.06×10 -7 2.97×10 62 Tortuosity 1.48 1.43 1.41 1.38 1.35 1.32 1.28 Another parameter that can show the geometric feature of the cathode microstructure is the tortuosity, which measures twist of the cathode surface. In this chapter tortuosity is evaluated by (4.29) where is tortuosity, element e, and is the surface area of the cathode microstructure at the boundary of is the projection of the elemental surface area to the vertical z-direction (not including the part that parallel to x-y plane). The numerator of EQ. 4.29 measures the real paths from the top of the cathode to the bottom, and the denominator measures the ideal shortest paths (directly along z-direction) from the top to the bottom. The reference tortuosity value of the reconstructed microstructure is 1.34, and the values of the tortuosity of the optimal microstructures are shown in Table 4.2. As the temperature decreases, the surface area of the tortuosity of the optimal design increases to make good use of the surface area and reduce the surface resistance of the cathode. However, at high temperatures, the tortuosity becomes smaller than the reference value, and the design focuses more on generating a thicker path to reduce the scaffold resistance. 4.6.2 Design with varying the arrangement of the cathode roots with constant footprint In this example the designs with different root arrangements are discussed. All the three arrangements shown in Fig. 4.5 are used. Fig. 4.9-4.10 and Table 4.2 illustrate the effect of varying the arrangement while keeping the footprint as 50% As shown in Table 4.3, on the optimal solution. remains nearly constant for all the three arrangements. At a high work temperature (higher 1/(RSe)), the optimized design for the five roots’ arrangement begins 63 to perform slightly better than the single roots’ arrangement. However, the arrangements of the roots show no significant influence to the cathode performance, and the improvements made by designing the cathode geometry remains the same (Table 4.4). μm 9 9 μm 9 8 8 8 7 7 7 6 6 6 5 5 5 4 4 4 3 3 3 2 2 2 1 1 1 μm 0 0 0.4 0.5 0 μm (a) Five roots 0 0 0.4 0.5 0 μm (b) Four roots 0 0 0.4 0.5 0 μm (c) Single root Figure 4.9 The optimized design for different root arrangements at 400°C. Performance is reported in Table 4.3 64 μm 9 μm 9 8 8 7 7 6 6 5 5 4 4 3 3 3 2 2 2 1 1 1 μm 9 8 7 6 5 4 0 0 0.4 0 (a) Five roots 0.5 μm 0 0 0.4 (b) Four roots 0.5 μm 0 0 0 0.4 0.5 μm 0 (c) Single root Figure 4.10 The optimized design for different root arrangements at 700°C. Performance is reported in Table 4.3 65 Table 4.3 Performance of optimized design for different root arrangements (Fig. 4.9-4.10) Temp (°C) Rp for Five Roots Rp for Four Roots 2 2 Rp for Single Root 2 (Ω×cm ) (Ω×cm ) (Ω×cm ) 400 2.50 2.51 2.49 700 0.0222 0.0223 0.0228 Table 4.4 Improvement of optimized design for different root arrangements compared to experimental data (Fig. 4.9-4.10) Temp (°C) Improvement for Five Roots Improvement for Four Roots Improvement for Single Root 400 20% 20% 20% 700 35% 35% 33% 4.7 Conclusions While the 2D designs can only give schematic guidance on the characteristics of 3D optimal microstructures, the 3D design emphasizes in direct design of the intricate cathode microstructure, and from the 3D design examples one can extract the optimal feature for different work temperature. Improvements of 21-35 % were found, compared to the experimentally measured data of the current cathode sample. Because the scaffold components in the 3D designs are permitted to grow in all the directions, the branches of the 3D designs are able to form complicated connections between them and form porous microstructures different from the 2D tree-like design. For the designs at low operating temperatures, the surface resistance of LSCF coating outweighs the bulk resistance of CGO scaffold, so the optimal designs have larger surface areas and tortuosities to make the geometric features thinner and span the whole height of the design domain to generate a larger surface within a minimum porosity. As the temperature increases, the bulk resistance of CGO scaffold begins to perform a more significant role than the surface resistance of LSCF coating. Therefore, the optimal geometries has smaller surface areas 66 and tortuosities and prefers to distribute at a lower height to form a thicker ionic conducting path to the bottom, and the resulting cathode height can be smaller than the height of the design domain. The cathode root arrangement only has very little influence on the cathode performance. 67 Chapter 5 Thermal strength design for anode support structure of solid oxide fuel cell electrode In this chapter, a finite element model for an anode-supported SOFC electrode is established in consideration of a coupled process of fuel supply, electrochemical reactions, heat transfer and structural deformation. The fuel flow in the channel with design-dependent geometry is modeled using a mixed formulation of steady-state Navier-Stokes equation and Darcy equation. The resulting velocity field is used as a known parameter in the calculation of the fuel transport. The fuel transport is described by the convection-diffusion equation. The current generation of the electrochemical reactions is calculated by the Butler-Volmer equation, and it is coupled with the fuel transport by a fuel consumption boundary condition. The heat is generated from the electrochemical reactions, and the heat exchange between the electrode and the flows results in a temperature gradient in the electrode. The temperature change results in a thermo-mechanical deformation of the electrode. The porous media used in SOFC electrodes are brittle ceramic materials, and the porous microstructure is composed by randomly arranged ceramic particles. The random nature in the microstructures results in the randomness of their inherent flaws from manufacturing, so the ceramics exhibit a large spread of bulk failure strength. Therefore, the failure stress of these materials is not a well defined number, and a probabilistic metric, which is called the probability 68 of failure, is utilized to assess whether the material will fail at a given stress or not. As the electrode is usually sintered at a higher temperature, tensile thermal stresses will be generated inside the electrode at the relatively lower work temperature. Experimental tests [57][58] show that thermal stress has a strong effect on the mechanical failure of the fuel cells. In the cases discussed in this chapter, the thermal stress is much higher than the stress caused by pressure load, and it is the dominate cause of mechanical failure. The failure probability of a structure is the probability of a fast fracture of structure in the instant after a load is applied. The weakest link theory can be used to calculate the probability of failure of a structure. When the weakest link fails, the whole structure will fail, so the structural strength depends on the weakest link. Weibull [22] gave an expression of the probability of failure, and this expression has been validated in the prediction of the structural failure of the SOFC electrode [28][29]. In this chapter, this expression is adopted as the measure of the structural strength. A topology optimization model is setup to design the macro-scale shape of the electrode and the fuel supply channel to optimize the thermal compliance and minimize the probability of failure. Although the 3D model can accommodate several design concepts, only the cross section design is considered, focusing on the concepts that are easier to manufacture. In addition, the thermal compliance is also a commonly used objective function for thermo-mechanical design. Further, as the probability failure increases with the amount of material, the minimum amount of material problem is also solved. The optimal design under these three objective functions are compared and discussed. 5.1 Analysis model This chapter focuses on the electrode design on the anode side, and the details of the oxygen gas 69 flow on the cathode side are not considered. A typical anode supported electrode design contains an active layer where the electrochemical reactions take place, and an anode support. The active layer is a thin layer (~50nm) printed on the top of the anode support, and it is composed by three layers: a cathode, an electrolyte and an anode layer. The active layer is printed upon the relatively thick (~600nm) anode support, as shown in Fig. 5.1. Interconnect Air channel Active Electrode Anode Support Fuel channel z Interconnect y Figure 5.1 The cross section of the fuel supply channel and the electrode As the active layer in the electrode is much thinner compared to the anode support, the effect of the active layer is simplified to be a boundary condition on the anode support layer, and the analysis only concerns the region occupied by the anode support and the fuel flow. Since it is difficult to make a curved active layer, the active layer is assumed to be a plane, and only the shape of the anode support is designed. The analysis model is composed by four parts: fuel flow, electrochemical reactions, fuel/heat transport and structural deformation. The analysis domain is selected to be a portion of the rectangular fuel supply channel (shown in Fig 5.2(a)). The porous electrode region unknown sub-domain of , and is an corresponds to the fuel flow region. The active layer of the electrode is assigned to be a plane on the top of the analysis domain, and the effect of the electrochemical reactions can be simulated as a boundary condition of the electrode. 70 Active electrode layer Analysis Domain (Ω) Active Layer Γ0 ξ=0 ω Γw Γω Γw Γin Γout Interconnect Ω\ω ξ=1 (a) Analysis Domain (Ω) Active Layer Interconnect Design ξ=? Domain (b) Figure 5.2 The analysis domain (a) and design domain (b) for the anode support Since the change in the concentration of oxygen gas is not significant, the air flow channel that locates above the active layer is not included in the analysis domain, and the cooling effect of the air flow is simulated as a convection boundary condition on the top boundary of the analysis domain. Aiming at the establishment of a topology optimization model, the analysis sub domain for both of the fuel flow and the anode support deformation are extended to the whole analysis domain. 71 The geometry of the electrode is described by a design variable field which has a value 0 in the electrode domain and 1 in the fluid domain. 5.1.1 The fuel flow model The fuel pumped into the fuel channel is a mixture of hydrogen and water steam. The flow is driven by a prescribed pressure drop . The flow is normally with a medium Reynolds number (1