Solid-Oxide Fuel Cell Electrode Microstructures: Making Sense of the Internal Framework Affecting Gas Transport Citation Hanna, Jeffrey (2010) Solid-Oxide Fuel Cell Electrode Microstructures: Making Sense of the Internal Framework Affecting Gas Transport. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/NRY2-2K63. https://resolver.caltech.edu/CaltechTHESIS:06072010-090420975 Abstract Optimal electrodes for solid-oxide fuel cells will combine high porosity for gas diffusion, high phase connectivity for ion and electron conduction, and high surface area for chemical and electrochemical reactions. Tracer-diffusion simulations are used to gain a better understanding of the interplay between microstructure and transport in porous materials. Results indicate that the coefficient of diffusion through a porous medium is a function of the details of the internal geometry (microscopic) and porosity (macroscopic). I report that current solid-oxide fuel cell electrodes produced from high-temperature sintering of ceramic powders severely hinder gas transport because the resulting structures are highly tortuous, complex three-dimensional networks. In addition, poor phase connectivities will assuredly limit ion and electron transport, as well as the density of active sites for power-producing reactions. With new access to a wide range of technologies, micro- and nano-fabrication capabilities, and high-performance materials, there is a new ability to engineer the fuel cell electrode architecture, optimizing the physical processes within, increasing performance, and greatly reducing cost per kilowatt. Even simple packed-sphere and inverse-opal architectures will increase gas diffusion by an order of magnitude, and provide a higher level of connectivity than traditional powder-based structures. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Tracer diffusion; SOFC; Triple-phase boundary; Tortuosity; Anode microstructure; Inverse opal Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Mechanical Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Shepherd, Joseph E. Group: Explosion Dynamics Laboratory Thesis Committee: Shepherd, Joseph E. (chair) Hunt, Melany L. Haile, Sossina M. Blanquart, Guillaume Defense Date: 27 May 2010 Record Number: CaltechTHESIS:06072010-090420975 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:06072010-090420975 DOI: 10.7907/NRY2-2K63 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 5923 Collection: CaltechTHESIS Deposited By: Jeffrey Hanna Deposited On: 08 Jun 2010 15:50 Last Modified: 16 Jan 2021 00:11 Thesis Files Preview PDF - Final Version See Usage Policy. 19MB Repository Staff Only: item control page

Solid-Oxide Fuel Cell Electrode Microstructures: Making Sense of the Internal Framework Affecting Gas Transport Thesis by Jeffrey Hanna In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 2010 (Defended May 27, 2010) ii c 2010 Jeffrey Hanna All Rights Reserved iii To everyone who has had to fight cancer and to those who continue their fight. iv Acknowledgments If only one paragraph of this gets read, I have to thank all those who have shaped who I am today and those who continue to do so. If you think about it, that goes back a long way for each of us: family, old friends and new friends, teachers I had as a kid (one of whom I must note by name, my high school physics teacher, Mr. Khan), experiences. Without all of that (and more), I would not even be writing anything that follows. I have been blessed throughout my life to meet some of the most amazing people and to be able to form close relationships with many of them. With complete respect and affection, thank you all. Let me first acknowledge my parents for always supporting me and encouraging me to pursue whatever I choose. For my entire life, they have offered nothing by love and patience, and continue to be the paramount example of how to treat others, how to live life, how to work hard, to persevere, to find strength in your weakest moments, and to make the best out of any situation. They understand that we all make mistakes, that things do not always go as planned, and that our true character shows in how we pick ourselves up, recover, and move on. That has been an invaluable lesson I have had to rely on throughout my time at Caltech. They taught me that life is not fair, and we have certainly had our share of experiences as a family and with others close to us. I will never forget January 7, 2007 for the loss of our dog Bailey, who was as much a part of the family as any one of us. That is surely one of the most (if not the most) painful experiences I had at Caltech, and life in general. We all miss you Bailey, we miss you so much. My older brother, Nick, and younger sister, Christine, also deserve thanks for so much, even the teasing. There is too much to say here about family and only limited space, so I guess beyond all the thanks I owe I want to say “We did it.” I finally finished school, and it was truly a group effort indeed. There are many people I have met at Caltech deserving of acknowledgment. My old roommate, v Won Tae Joe, definitely made the first year more enjoyable despite all the work. Won Tae continued to live with me for several more years, and now I am officially addicted to Korean food. Andrew Tchieu soon joined Won Tae and me, and has been a constant source of sarcastic humor, and a good friend ever since. The three of us have definitely battled tough times here at Caltech, and somehow we all managed to stick it out. Emily McDowell has been my connection to all the more interesting gossip at Caltech. I seriously think I ended up with some of the best people to live with at Caltech, and I am thankful to have these friendships as I move on in life. Dom Rizzo is the newest roommate, and probably one of the better cooks you will find (aside from my mom, of course). I also finally met someone here that is as much of a beer snob as I am. Aside from roommates, my recent officemates Pablo Abad-Manterola and Dave Pekarek have made the office seem like a second home. I guess it kind of takes on that feel when you have a lounge chair and a couch in your office. Many thanks also to the 4:30 pm Friday Beer Hour crew and all the people I have met through that. Again, so much to say and so little space, so thanks to everyone that got me away from the office and reminded me that life truly is not about science. While science has certainly improved life (and most certainly has been used to destroy life, unfortunately), too many people at Caltech believe science is life. Life is about the interactions, relationships, and experiences we are exposed to on the way, and science is simply the common thread that brought us all here in the first place. I must also thank many of the faculty and staff here at Caltech and beyond who have helped me with my research. My previous advisor, David Goodwin, is one of the most intelligent people I have ever met. An unfortunate diagnosis and battle with brain cancer prevented him from continuing to serve as my academic advisor, but I take with me an indelible mark of his passion for science and desire to create a world where we are more responsible with the resources we have on Earth and the way in which we use them to create the energy we rely on every day—please, continue fighting. A deserving thanks to all those who stepped in during this unfortunate set of circumstances, including my new advisor Joe Shepherd who has been just incredible with guiding me and talking to me about many issues at Caltech. Thanks to my group members Francesco Ciucci and Vaughan Thomas who also faced the same problems (I am hesitant to call our issues “problems” because they are trivial in the whole scheme of things)—the three of us managed to support each other and give help or suggestions vi when we could. And thanks to my Master’s advisor from Colorado School of Mines, Bob Kee, who also stepped in during this time and has offered a great deal of support and guidance. Tim Colonius was perhaps the first to make sure we knew support was there and was instrumental in making sure we were able to continue with our work. Finally, many thanks to my committee members, Professors Guillaume Blanquart, Melany Hunt, and Sossina Haile for being generous with their time and knowledge. Other Caltech faculty and staff who deserve my gratitude include Rob Phillips, professor in APh/BMB/ME; John Doyle, professor in CDS/EE/BE; and Coach Clint Dodd (you can find him in the little office next to the pool, or somewhere on the pool deck). Strangely, I know nothing of Rob Phillips outside of a handful of encounters either at, or walking to or from the Caltech pool/gym. From our first conversation, it is obvious that he is genuinely concerned with student life here at Caltech, and the well-being of the student community. I met John Doyle through triathlon, and he has become a great friend. He is also one of the most intelligent people I know, and he and Rob Phillips both understand that work needs to be balanced with activity outside of the office or lab. John is very enthusiastic about his work, and equally dedicated to his fitness—just do not follow him on any trail runs. Coach Dodd is perhaps the most likable, good-humored, and most entertaining person on the Caltech campus. I will not even begin to ramble about how much fun it is to be around him. Anyone reading this should seek out those mentioned above and get to know them—they will certainly improve your quality of life here at Caltech. I certainly must give my thanks and appreciation to the ME 19 ab Fluid Mechanics class, which I TAed for the fall and winter terms of the 2009–2010 academic year. I never imagined that I would have so much fun and would be able to connect with so many students. To all of those who regularly attended office hours, or stopped by my office, or even made an effort to say hello outside of the classroom, you made me feel like all the effort and all the time I committed was worth every moment. I wish you all the very best as you continue through Caltech and beyond, and nothing but success in the rest of your lives. I encourage you to find something you love doing and run with it—many of you will go on to do great things. I must also give thanks to the class professor, Chris Brennen, for allowing me the wonderful experience to serve as TA, to offer my input, and just what a pleasure the opportunity has been. I hope to one day stand in his position at the head of a class with all the enthusiasm and passion vii he devotes to teaching and the concern he demonstrates about whether students are learning. Though I have long had the suspicion that I would love teaching, the opportunity to serve as TA of ME 19, and the students and their personalities, have more than confirmed my belief. There are three people from Caltech that I cannot possibly sum up in a few brief sentences. To call them best friends would still understate their relationship to me and the impact they have had. Ricky Lim, I will use your own words here: “I’ve found a great friend and a brother. We will be there/anywhere at a moment’s notice for one another. I don’t know anyone I trust more with any matter insignificant or life altering.” (June 16, 2005). Scott Hersey, you are good people. It is literally impossible for me to sum up our experiences or conversations. I invested my friendship units wisely— you are a brother to me. Perhaps in another lifetime we will find ourselves in Cormac McCarthy’s, All the Pretty Horses. And lastly, deservingly so, Kakani Young has been with me from almost the beginning of my time here. We have endured a lot together. For all our experiences good and bad, and without taking away from anyone else, I must say that no matter what happens in the future, you alone made the entire last five years or so worth the effort. My apologies for the long-winded and lengthy messages. I truly believe there are many deserving of my thanks, and probably many I failed to mention here. I understand that words are not the proper way to show appreciation—that can only be accomplished through actions. However, this space provides me the opportunity to express my gratitude and document (however briefly) some of the experiences I have had at Caltech. That said, I again say thank you to all those that have shared this experience with me. Finally, this work was supported by the DoD Multidisciplinary University Research Initiative (MURI) program administered by the Office of Naval Research under Grant N00014-05-1-0339. viii Abstract Optimal electrodes for solid-oxide fuel cells will combine high porosity for gas diffusion, high phase connectivity for ion and electron conduction, and high surface area for chemical and electrochemical reactions. Tracer-diffusion simulations are used to gain a better understanding of the interplay between microstructure and transport in porous materials. Results indicate that the coefficient of diffusion through a porous medium is a function of the details of the internal geometry (microscopic) and porosity (macroscopic). I report that current solid-oxide fuel cell electrodes produced from high-temperature sintering of ceramic powders severely hinder gas transport because the resulting structures are highly tortuous, complex three-dimensional networks. In addition, poor phase connectivities will assuredly limit ion and electron transport, as well as the density of active sites for power-producing reactions. With new access to a wide range of technologies, micro- and nano-fabrication capabilities, and highperformance materials, there is a new ability to engineer the fuel cell electrode architecture, optimizing the physical processes within, increasing performance, and greatly reducing cost per kilowatt. Even simple packed-sphere and inverse-opal architectures will increase gas diffusion by an order of magnitude, and provide a higher level of connectivity than traditional powder-based structures. ix Contents 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Organization of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Kinetic Theory and Gas Diffusion 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Statistical Methods and Distribution Functions . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Maxwell-Boltzmann and the Velocity Distribution . . . . . . . . . . . . . . . . . 10 2.2.1.1 Sampling Molecular Speeds from the Equilibrium Distribution . . . . . 18 Scattering of Molecules from a Surface . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.2.1 Sampling from the Angular Probability Distributions . . . . . . . . . . 22 2.2.2.2 Sampling from the Speed Distribution at a Surface . . . . . . . . . . . . 23 Distribution of Free Paths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.3.1 Sampling from the Free-Path Distribution . . . . . . . . . . . . . . . . . 26 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.1 Normal Diffusion in Capillaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.2 Diffusion as a Random Walk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.3.3 Diffusion in Porous Media and Free-Molecule Diffusion in Capillaries . . . . . . . 33 2.2.2 2.2.3 2.3 3 Mass Transport in Porous Structures 37 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Simulation Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 x 3.3.1 Tracer Diffusion in Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.2 Tracer Diffusion Between Aggregates of Sphere Particles . . . . . . . . . . . . . . 47 3.3.2.1 Generation of Porous Solids Using Packed Spheres . . . . . . . . . . . . 48 3.3.2.2 Monte-Carlo Simulation of Tracer Molecules in Particle-Based Geometries 49 3.3.3 Tracer Diffusion in Three-Dimensional Voxel-Based Geometries . . . . . . . . . . 3.3.3.1 3.5 Creating Three-Dimensional Cermet Data from Two-Dimensional FIBSEM Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.3.2 Surface Normals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.3.3.3 Monte-Carlo Simulation of Tracer Molecules in Voxel-Based Geometries 57 General Issues With Tracer Simulations in Systems With Complex Geometry . . 59 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.4.1 Diffusivity Results for Cylindrical Pores and Other Capillary Models . . . . . . . 59 3.4.2 Diffusivity Results for Packed-Sphere Models . . . . . . . . . . . . . . . . . . . . 60 3.4.3 Diffusivity Results for Particle- and Voxel-Based SOFC Models . . . . . . . . . . 62 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.3.4 3.4 52 4 Turning Porous Media Inside Out: The Pore-Size Distribution 67 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4.1 Short and Long Cylinders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4.2 Inside a Single Closed Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4.3 Monodisperse Simple Cubic Packing of Spheres . . . . . . . . . . . . . . . . . . . 75 4.4.4 Monodisperse Body-Centered Cubic Packing of Spheres . . . . . . . . . . . . . . 79 4.4.5 Monodisperse Face-Centered Cubic Packing of Spheres . . . . . . . . . . . . . . . 79 4.4.6 Random Monodisperse Sphere Particles . . . . . . . . . . . . . . . . . . . . . . . 81 4.4.7 Real Structures: Particle- and Voxel-Based Reconstruction of a Solid-Oxide Fuel Cell Anode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 xi 4.4.8 4.5 Moments of the Chord Distribution . . . . . . . . . . . . . . . . . . . . . . . . . 86 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5 Piecing Together the Puzzle: SOFC Anode Microstructure 91 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.3 The Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6 An Alternative Vision: Architecturally Engineered Fuel Cell Structures 121 6.1 Summary and Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.2 A First Glance Inside a New Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.3 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Bibliography 128 xii List of Figures 1.1 Two cells of a planar anode-supported SOFC stack. . . . . . . . . . . . . . . . . . . . . . 3 1.2 One cell of a tubular anode-supported SOFC. . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Details of the membrane-electrode assembly microstructure affect transport and chargetransfer processes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 SEM images of typical fuel cell structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 A view of a three-dimensional reconstructed SOFC anode. . . . . . . . . . . . . . . . . . . 6 2.1 Representative control-volume elements in physical and velocity space. . . . . . . . . . . . 11 2.2 Plot of the probability distribution function for molecular speeds, f (c0 ). . . . . . . . . . . 18 2.3 Spherical polar coordinates in physical space. . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 Molecules crossing a surface dA in direction dΩ with speed c0 must lie in the volume element c0 cos θdAdt. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Three-dimensional polar plot of the angular distribution functions. . . . . . . . . . . . . . 22 2.6 Plot of the probability distribution function for the polar angle, h(θ) = 2 cos θ sin θ, for molecules incident on or leaving a surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 23 Comparison of the probability distribution functions for a Maxwell-Boltzmann speed distribution, f (c0 ), and the distribution of speeds across a surface, f s (c0 ). . . . . . . . . . . . 24 2.8 Plot of the probability distribution function for free paths, f (`) = e−`/λ . . . . . . . . . . . 26 3.1 Knowledge and understanding at certain scales informs design and optimization at other scales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 40 Transmission probability of molecules through a diffuse-reflecting cylindrical tube of length L and radius a for Knudsen flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 xiii 3.3 Transmission probability of molecules through a 90◦ cylindrical elbow calculated using a Monte-Carlo algorithm for Knudsen flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 46 Transmission probability of molecules through a cylindrical annulus calculated using a Monte-Carlo algorithm for Knudsen flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.5 Classic examples of ordered, isotropic, monodisperse packed spheres. . . . . . . . . . . . . 48 3.6 Randomized array of uniform-size sphere particles at 39.6% porosity, and a packing of various-sized spheres at 17.5% porosity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.7 Pictorial illustration of gas behavior under different flow conditions. . . . . . . . . . . . . 50 3.8 FIB-SEM data as a “.tif” image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.9 Image comparison with original image and reduced image. . . . . . . . . . . . . . . . . . . 53 3.10 Image file and resulting particle array. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.11 Voxel- and particle-based reconstructions of a SOFC anode. . . . . . . . . . . . . . . . . . 55 3.12 Illustration of how the surface normal is determined for a two-dimensional pixel system. . 57 3.13 Diffusibility versus porosity for various computational porous media composed of uniform sphere particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 60 Diffusibility versus porosity for various experimental and computational porous media composed of sphere particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.15 The data from figure 3.14 with included results from the SOFC models. . . . . . . . . . . 65 4.1 Probability distributions of chord lengths for molecules “bouncing” around inside cylinders of various aspect ratios L/R. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 72 Cumulative probability distributions of chord lengths for molecules “bouncing” around inside cylinders of various aspect ratios L/R. . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3 Effect of cylinder aspect ratio, L/R, on mean chord length. . . . . . . . . . . . . . . . . . 74 4.4 Probability distribution of chord lengths for molecules “bouncing” around inside a sphere. 75 4.5 Cumulative probability distribution of chord lengths for molecules “bouncing” around inside a sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 75 Probability distributions for chord lengths in the void space between simple cubic arrangements of spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 xiv 4.7 Cumulative probability distributions for chord lengths in the void space between simple cubic arrangements of spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Probability distributions for chord lengths in the void space between body-centered cubic arrangements of spheres. 4.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Probability distributions for chord lengths in the void space between a random arrangement of spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 80 Cumulative probability distributions for chord lengths in the void space between facecentered cubic arrangements of spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 78 Probability distributions for chord lengths in the void space between face-centered cubic arrangements of spheres. 4.11 78 Cumulative probability distributions for chord lengths in the void space between bodycentered cubic arrangements of spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 77 81 Cumulative probability distributions for chord lengths in the void space between a random arrangement of spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.14 Particle-based model of the reconstructed SOFC anode. . . . . . . . . . . . . . . . . . . . 83 4.15 Probability distribution for chord lengths in the pore phase of a particle-based threedimensional SOFC reconstructed anode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.16 Voxel-based model of the reconstructed SOFC anode. . . . . . . . . . . . . . . . . . . . . 85 4.17 Probability distribution for chord lengths in the pore phase of a voxel-based three-dimensional SOFC reconstructed anode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Cumulative probability distribution for chord lengths in the pore phase of a particle-based three-dimensional SOFC reconstructed anode. . . . . . . . . . . . . . . . . . . . . . . . . 4.19 85 86 Cumulative probability distribution for chord lengths in the pore phase of a voxel-based three-dimensional SOFC reconstructed anode. . . . . . . . . . . . . . . . . . . . . . . . . 86 4.20 Mean chord lengths of the pore space in sphere-based structures. . . . . . . . . . . . . . . 87 4.21 Mean-squared chord lengths of the pore space in sphere-based structures. . . . . . . . . . 88 4.22 A measure of the deviation of a chord-length distribution from an exponential distribution, based on the first and second moments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 xv 5.1 Expanded view of the triple-phase boundary region in a membrane-electrode assembly. . 92 5.2 Three-dimensional visualization of the networks constituting the Ni phase. . . . . . . . . 93 5.3 Triple-phase boundary curves of the voxel anode are made up of line segments. . . . . . . 96 5.4 Triple-phase boundary curves of the particle anode are made up of discretized circular arcs. 97 5.5 Distribution of network phase volumes for Ni in the three-dimensional SOFC reconstructions.108 5.6 Distribution of network phase volumes for YSZ in the three-dimensional SOFC reconstructions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.7 Distribution of network phase volumes for pores in the three-dimensional SOFC voxel reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.8 Highlights of the “isolated” and “dead-end” networks for the Ni phase in the voxel anode. 111 5.9 Highlights of the “isolated” and “dead-end” networks for the Ni phase in the particle anode.112 5.10 Highlights of the “isolated” and “dead-end” networks for the YSZ phase in the voxel anode.113 5.11 Highlights of the “isolated” and “dead-end” networks for the YSZ phase in the particle anode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.12 Highlights of the “isolated” and “dead-end” networks for the pore phase in the voxel anode.115 5.13 Highlights of the “isolated” networks for the pore phase in the particle anode. . . . . . . 116 5.14 Distribution of triple-phase boundary length for the three-dimensional SOFC reconstructions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.15 Highlights of the “isolated”/inactive, “dead-end”/unknown, and “across”/active triplephase boundary curves for the three-dimensional SOFC voxel anode. . . . . . . . . . . . . 118 5.16 Highlights of the “isolated”/inactive, “dead-end”/unknown, and “across”/active triplephase boundary curves for the three-dimensional SOFC particle anode. . . . . . . . . . . 118 5.17 Three-dimensional visualization showing the dependence of escape probability on number of simulated molecules and number of maximum allowed collisions. . . . . . . . . . . . . . 119 5.18 Sensitivity analysis to determine a sufficient simulation time for molecules to escape the particle-based three-dimensional SOFC anode. . . . . . . . . . . . . . . . . . . . . . . . . 120 6.1 Conceptual design of an engineered, high-performance membrane-electrode assembly. . . 123 6.2 Example of inverse opals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 xvi 6.3 The oxide framework. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.4 Diffusibility results for a uniform periodic simple cubic array of inverse opals. . . . . . . . 125 xvii List of Tables 3.1 Breakdown of the size and number of particles used in the particle reconstruction of a SOFC anode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Summary of Knudsen transport data for the three-dimensional SOFC particle-reconstructed anode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 64 Volume, surface area, and connectivity information for the networks of each phase in the three-dimensional SOFC voxel -reconstructed anode 5.2 63 Summary of Knudsen transport data for the three-dimensional SOFC voxel -reconstructed anode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 56 . . . . . . . . . . . . . . . . . . . . . 100 Volume, surface area, and connectivity information for the networks of each phase in the three-dimensional SOFC particle-reconstructed anode . . . . . . . . . . . . . . . . . . . . 100 5.3 Summary of the triple-phase boundary density, activity, and curve lengths in the threedimensional SOFC voxel -reconstructed anode . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.4 Summary of the triple-phase boundary density, activity, and curve lengths in the threedimensional SOFC particle-reconstructed anode . . . . . . . . . . . . . . . . . . . . . . . 103 6.1 Comparison of Knudsen transport data for inverse opals and the three-dimensional SOFC reconstructed anodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 1 Chapter 1 Introduction “[The fuel cell] is almost magical in its elegance and simplicity, and it is astonishing that this process has not yet been commercialised to supplant the inefficient and polluting combustion heat engines which currently dominate our civilization.” [54] As of the writing of this thesis, the need for cheaper, cleaner, renewable, and more reliable energy sources and conversion methods is arguably at the forefront of our society. Rising energy costs, limited natural resources, a reliance on those resources in unstable geographic regions, rising CO2 emissions, and the threat of global climate change have moved money and minds to deliver more efficient ways of extracting energy from already limited supplies, and to find more environmentally conscience energy sources or conversion devices. A continuum of work spans the realm of energy-related topics, each with its own attractions and drawbacks, as well as present-, short-, and long-term goals. There is no consensus on a right or wrong solution, but the scientific community and those far removed are of a common understanding that something must change. Much effort is focused on improving current combustion systems, using bioinspired fuels, harnessing renewable sources like water, wind, and solar, and/or combining one or more of these areas with new electrochemical conversion devices. Solid-oxide fuel cells (SOFCs) are the most efficient such devices known to convert chemical energy in a fuel directly into electrical power [54], without burning the fuel. Originally, SOFCs were designed to compete with large power-generation units, like central power stations. In these large-scale systems with cogeneration of power and useful heat, efficiencies as high as 75% are projected. These fuel cells operate 2 at relatively high temperatures, 600 to 1,000 ◦ C, making start-up and shutdown slower and generally more difficult. Because of this, SOFCs are not particularly suited for applications requiring rapid response or power cycling (e.g., automotive power), but during the last 10–15 years, the realization has steadily dawned that SOFCs can work well in small, portable, residential and auxiliary power systems [54]. Unlike proton exchange membrane (PEM or hydrogen) fuel cells, SOFCs can operate on many different fuels; they work very well on hydrogen, but also operate nearly as well on methanol, ethanol, methane, diesel reformate, and more. As a result of their high efficiency and fuel flexibility, SOFCs hold a promising position for near-term power generation for the following key reasons: • high efficiencies directly correlate to a significant reduction in greenhouse gas emissions compared to more traditional coal-burning power systems; • the fuel flexibility of SOFCs allows operation on emerging fuels such as biofuels, coal-derived syngas, and more—this is particularly promising for near-term adoption as well as being a significant solution to the long-term, post-fossil-fuel, hydrogen economy; • further, SOFCs can serve as a CO2 capture technology through the use of synthetic fuels produced by reaction of hydrogen with atmospheric CO2 [5]; • SOFCs can integrate directly with existing hydrocarbon infrastructure; • SOFCs can provide power from milliwatt to large-scale megawatt power-generation systems, providing an unprecedented range of applicability. Despite the clear promise of SOFCs, none of this will come to pass unless we develop SOFCs that are substantially cheaper, per kilowatt of generating capacity, than current-generation SOFCs. To bring cost per kilowatt down (the Department of Energy target of $400/kW requires an order-of-magnitude reduction!), it is necessary to either dramatically lower manufacturing costs per cell, or dramatically improve the power output per cell (or both). While some performance improvements are possible by optimizing the overall, or big-picture, macroscopic-level parameters (stack design, interconnect locations, etc.), these improvements are incremental in nature and are actively being pursued by industrial SOFC developers. Whatever the macroscopic architecture, the heart of any SOFC is the membrane-electrode 3 Fig. 1: Two cells of a planar anode-supported solid-oxide fuel-cell Fig. 2: One cell of a tubular anode-supported solid stack. The flow channels, which have characteristic cross sections Figure 1.1. Two cells of a planar anode-supported SOFC stack. The flow channels, which have Fuel flows through the characinside of the tube, which is on order the orderofofaa square millimeter, are formed the interconnect. porous anode. The dense is electrolyte and poro teristic cross sections on the square millimeter, areinformed in the interconnect. ThethinMEA The membrane-electrode assembly is fabricated with a dense elecfabricated on the outside of the tube. The outside o fabricated with a dense electrolyte sandwiched between porous electrodes. From Kee et al. [31]. trolyte sandwiched between porous electrodes. ode) is exposed to air. Wires are wrapped around the tube to form the cathode interconnect and the is formed by a metal insert that is bonded to the p chemistry and physical processes withinlayers fuel cells. Al- by assembly (MEA), consisting of porous anode and cathode separated dense ceramic eleconeaend of the tube. though aspects of the fluid flow, gas-phase chemistry, and thermal catalytic chemistry will be familiar to readers who trolyte. The opportunity for game-changing relies aspects on ourofunderstanding of thedimensions MEA, on the order small. With channel are expert in combustionbreakthroughs science and technology, ter and low flow rates (especially on the fuel the electrochemistry will be new. Fuel cells rely on the couteristic Reynolds numbers are on the order of pled interaction between thermal and electrochemical pro- and specifically, the interplay between electrode microstructure, gas transport, electrochemical reacFigure 2 illustrates an anode-supported tu cesses. Thus, the optimal design and implementation of sign. Here the fuel flows on the inside of the new systems will benefit from an improved understanding tions responsible for powerofproduction. Attempts to understand SOFC electrode performance havetubes are packed air on the outside. Multiple these interactions. ber through which air flows. Air and exces zone ataffects the tube exit. Tube been limited by a lack of 2. data and knowledge which directly Solid-Oxide Fuel Cells describing the microstructure,a combustion an advantage insofar as dissimilar-material s less troublesome than in planar systems. Ho Figure 1 illustrates theelectrodes physical layout a planar anodegas transport and electrochemistry. Optimal willof combine high porosity gas collection diffusion,can be more tro losses infor current supported SOFC. The fuel-cell stack is composed of layers cause there are longer current paths through of cells. Each cell is composed of a membrane-electrode ashigh-resistance electrodes. In any case, the un high connectivity for ion/electron conduction, and high surface area for chemical and electrochemical sembly (MEA), which consists of a dense electrolyte sandtrochemical and reacting flow processes are s wiched between an anode and a cathode. The anode is the nar and tubular systems. electrode on the fuel side and the cathode is on the air side. reactions. The electrolyte is a dense polycrystalline ceramic that is 2.1. Electrolyte impervious to any gas transport, but conducts oxygen-ions 2− Although much SOFC O research focuses on developing materials, it is also . The architecture illustrated in Fig.new 1 has high-performance defined flow Electrolytes should be good ion conductors channels that are formed in an interconnnect. The anode tle or no electronic conductivity. The electro and cathode are connected bysignificantly an external circuit through known that electrode structure and microstructure impact SOFC be performance. Thetransport. work In an SOFC t impervious to gas which electrons may flow. The voltage of an individual cell is an oxide with mobile oxygen ions O2− . T ranges from about one volt at open circuit to around onemongas electrolyte is zirconia (ZrO2 ) with 8% y in this thesis is aimed to half expose the relationship between transport. This volt at maximum power density. Themicrostructure system voltage and as a dopant (yttrium-stabilized zirconia, or can be increased to higher levels by stacking a number of atoms substitute for Zr in the lattice, but hav cells in series, as stacking batteries flashlight.framework is accomplished in two separate broadmuch tasks: (1) study thein ainternal various ratherof than +4. To modeled compensate the charge, on Figure 1 shows a counter-flow situation, where the fuel must be removed for every two Y atoms adde flows in one direction and the air flows in the other. Both YSZ has atovacancy concentration of 4% in th microstructures spanning from seemingly simple, uniform-size sphere more complex co-flow and cross-flow geometries are also ordered used. Other lay- packings, lattice. Since the addition of Y dopants also outs do not use defined channels, leaving the space between relatively open fluorite crystal structure (whi electrode and the interconnect relativelymeshes open with random and distributed-sizetheparticle-based models, to voxelated of actual SOFC electrodes; and temperatures), doping is only stable at higher screen-like structures forming the electrical connection becan move relatively easily through the crystal tween the porous electrode and the interconnect structure. Theproperties ionic conductivity σ is known for m (2) simulate gas diffusion within microstructure to ofpredict effective transport and lead Figure 1the illustrates only the interior the stack. There must electrolyte materials and is usually expressed also be some sort of manifold system that distributes flow ! " channels and between layers of cells molecular [4]. The interσ0 of −Eel to better understanding ofinto thethecomplex relationship between mechanisms and influence exp . σ= connect should have very low electrical resistance, leading T RT to nearly uniform electric potential. The overall system size the physical structure of the pore on network. For YSZ, σ0 ≈ 3.6 × 105 S·K/cm and Ee depends the net power to be delivered. For systems in J/mol [5]. It is evident that the conductivity v the kilowatt range, planar cell dimensions are on the order with temperature. of tens of square centimeters and flow channels are usually 2 4 o cells of a planar anode-supported solid-oxide fuel-cell Fig. 2: One cell of a tubular anode-supported solid-oxide fuel cell. flow channels, which have characteristic crossOne sections Figure 1.2. cell of Fuel a tubular anode-supported SOFC. Fuel isflows through flows through the inside of the tube, which formed from a the inside of the tube, er of a square millimeter, are formed in the interconnect. porous anode. thinthin dense electrolyte and porous which is formed from a porous anode.The The dense electrolyte andcathode porousare cathode are fabricated on brane-electrode assembly is fabricated with a dense elec-tube. fabricated on the of outside the tube. outside of tube (caththe outside of the The outside the of tube, i.e.,The cathode, is the exposed to air. Wires are wrapped dwiched between porous electrodes. ode) is exposed to the air. Wires are interconnect wrapped aroundand the outside of terminal is formed by around the outside of the tube to form cathode the anode the tube to form the cathode interconnect and the anode terminal a metal insert that is bonded to the porous anode at one end of the tube. From Kee at al. [31]. is formed by a metal insert that is bonded to the porous anode at y and physical processes within fuel cells. Alone end of the tube. spects of the fluid flow, gas-phase chemistry, and 1.1 Background atalytic chemistry will be familiar to readers who small. With channel dimensions on the order of a millimet in combustion science and technology, aspects of ter and low flow rates (especially on the fuel side), characochemistry will be new. Fuel cells rely on the cou1.1 illustrates physical layout of a planar Reynolds numbers are on anode-supported the order of ten. SOFC, while figure 1.2 shows raction between thermal and Figure electrochemical pro- theteristic Figure 2 illustrates an anode-supported tube-based deThus, the optimal design and implementation of sign. Here the fuel flows on the of the tube and the an anode-supported design. The anode is ainside ceramic-metallic composite (cermet), generally ems will benefit from an improved understanding tube-based air on the outside. Multiple tubes are packed into a chamnteractions. through whichzirconia air flows.(YSZ), Air and excess mix is in composed of strontiummade of nickel (Ni) and ber yttria-stabilized while the fuel cathode a combustion zone at the tube exit. Tube systems have Oxide Fuel Cells an advantage as dissimilar-material sealing canlayers be is a dense electrolyte doped lanthanum manganate (LSM).insofar Sandwiched between the electrode less troublesome than in planar systems. However, ohmic 1 illustrates the physical layout of a planar anodelosses in current collection can be more troublesome bed SOFC. The fuel-cell stack is composed of layers (YSZ) that has the ability to conduct oxygen ions (O2− ) at through high temperatures (SOFCs generally operate cause there are longer current paths the relatively Each cell is composed of a membrane-electrode ashigh-resistance electrodes. In any case, the underlying elecMEA), which consists of a dense electrolyte sand◦ in the range 600–1,000 C). Collectively, these three make the MEA. trochemical and reacting flow pieces processes are up similar in pla- Mass transfer of the fuel etween an anode and a cathode. The anode is the nar and tubular systems. on the fuel side and the cathode is on the air side. (e.g., ceramic H2 ) occurs trolyte is a dense polycrystalline thatfrom is the fuel channel, through the porous anode, to electrochemical reactive sites at or 2.1. Electrolyte us to any gas transport, but conducts oxygen-ions near the defined anode/electrolyte interface. These reactive sites are known as triple-phase boundaries (TPBs), e architecture illustrated in Fig. 1 has flow Electrolytes should be good ion conductors, but have litthat are formed in an interconnnect. The anode tle or no electronic conductivity. The electrolyte must also formed at thethrough intersection of gas, oxygen-ion conductor (electrolyte, YSZ), and electron conductor (Ni). ode are connected by an external circuit be impervious to gas transport. In an SOFC the electrolyte ectrons may flow. The voltage of an individual cell is an oxide with mobile oxygen ions O2− . The most comom about one volt at open circuit to around Products of theonecharge-transfer reactions (i.e., gas-phase products like H2added O, as well as electrons) must mon electrolyte is zirconia (ZrO2 ) with 8% yttrium at maximum power density. The system voltage as a dopant (yttrium-stabilized zirconia, or YSZ). The Y creased to higher levels by stacking a number of from be transported away the substitute TPB. Gas-phase and atoms for Zr in products the lattice,diffuse but have +3convect valence,back through the porous eries, much as stacking batteries in a flashlight. rather than +4. To compensate the charge, one oxygen ion 1 shows a counter-flow situation, where the fuel for every two atoms added. Therefore, anode where they mix must with be theremoved fuel. Electrons areYcarried through the anode metal phase to the one direction and the air flows in the other. Both YSZ has a vacancy concentration of 4% in the oxygen subnd cross-flow geometries are also used. Other laylattice. Since thegraphically addition of Y also stabilizes the interconnect. The process is illustrated in dopants figure 1.3. ot use defined channels, leaving the space between relatively open fluorite crystal structure (which without Y rode and the interconnect relatively open with doping is of only stable at labs, higheressentially temperatures), the vacancies Apart from work handful university all SOFC development work is a variation ke structures forming the electrical connection be- in acan move relatively easily through the crystal via hopping. e porous electrode and the interconnect structure. The ionic conductivity σ is known for most common illustrates only the interior of on the basic stack. themes There must more than 40 years old. Forand example, the Ni/YSZascermet anodes used in the majority electrolyte materials is usually expressed ome sort of manifold system that distributes flow ! " hannels and between layers of [4]. The inter-do not differ qualitatively from −Edeveloped σ0 those el ofcells modern SOFCs by Spacil(1) in 1964 [56]. The modern exp . σ= should have very low electrical resistance, leading T RT uniform electric potential. The overall system sizeanode, YSZ electrolyte, and LSM cathode—is the result of engineering compromises SOFC—a Ni/YSZ For YSZ, σ0 ≈ 3.6 × 105 S·K/cm and Eel ≈ 8 × 104 on the net power to be delivered. For systems in J/mol [5]. It is evident that the conductivity varies strongly att range, planar cell dimensions are on the order with temperature. f square centimeters and flow channels are usually 2 5 Figure 1.3. Details of the membrane-electrode assembly microstructure affect transport and chargetransfer processes. Image courtesy of Robert Kee, Colorado School of Mines. that took into account manufacturability, material compatibility (e.g., matching of thermal expansion coefficients), performance, and cost. These trade-offs were established for the previous era of manufacturing technology, when high-temperature sintering of ceramic powders was at the forefront of ceramics fabrication methods. Now, we have access to a wide range of technologies, micro- and nanofabrication capabilities, and high-performance materials largely unknown in previous decades. Current SOFC electrodes are typically produced from high-temperature sintering of powders embedded with pore formers. The resulting structures are highly tortuous, complex three-dimensional (3D) networks, optimal neither for gas transport nor electrochemistry. Not until the recent work of Wilson et al. [65] have we been able to visualize and analyze the microstructure of a SOFC cermet anode. The results of a 3D reconstruction of a cermet anode using focused ion-beam milling and scanning electron microscopy (FIB-SEM) are shown in figure 1.5. The microstructure is highly convoluted and only partially connected—only 63% of the TPB regions are well connected to the metal, oxide, and pore networks. A large proportion (19%) of the TPB is formed by short segments disconnected from one or more networks, and therefore inactive for electrochemistry. Such a random, tortuous structure is clearly not optimized for its desired purpose, and clearly not what one would design with better control over the fabrication process. 6 Figure 1.4. SEM images of typical fuel cell structures. Left shows the membrane-electrode assembly in cross section and right shows the porous surface of the anode in plan view. Pictures and description from Shao et al. [53]. 5.2 5 4 3.4 3 2 0 0 1 2 3 x (µm) 4 5 6 z (µ 1 1.7 m) 2 1 y (µm) y (µm) 3 0 0 0 2.0 x (µm) Figure 1.5. A view of a Figure three-dimensional reconstructed SOFC anode, showing 2 3D anode reconstruction. A view of the 3D reconstruction showing theNi Ni (green), YSZ (transluFigure 3 3D map of the three-phase boundari cent/gray), and pore (blue) phases. Picture and and description (green), YSZ (translucent/grey), pore (blue) phases. from Wilson et al. [65]. represents a set of contiguous TPBs. The majori connected (coloured white/grey). The remaining disconnected TPB segments Using advanced materials and advanced nanoengineering methods, combined with a better under-(having colours oth these intersect the sample boundaries, and henc 6 −3 4.28 × 10 m cm . Although TPB lengths have previously been segments existing outside the sample volume. H 10 Sr MnO –YSZ cathodes , this estimated from 2D images of La 0.2 transport 3 standing of the interplay between microstructure0.8and (as well as electrochemistry), it is of the TPBs contact neither the highly inteconnec is the first direct measurement known to the authors. boundaries, that is, they are actual short segmen We are currently developing detailed electrochemical models possible to completelyofchange SOFC electrodes are made,that resulting potentially much higher poroushow anode and cathode performance connect in TPB length and other microstructural information quantitatively to performance, and greatly reducing In cost kilowatt. performance. theper meantime, a crude assessment can be done 18% of the TPBs are undetermined to determine if the measured TPB density reported above is sample boundaries. A related question is tortuosity o consistent with the measured SOFC performance. We can estimate the expected anode area-specific resistance by multiplying the gas-phase transport pathways within above TPB density by previously measured values of the length- example, a high gas-phase tortuosity 1.2 Organization of Thesis specific resistance (LSR) of the TPB, and an approximate anode exchange of H2 O and H2 ; this can active depth (≈10 µm (refs 17,21)). On the basis of various studies is most readily observed as a limiting including point-contact18,22 or patterned23,24 electrodes, reported current densities (for example, the rap −2 The goal of this thesis LSR is tovalues provide picture of the a traditional seenanode, in Supplementary I to 5 × 106 !ofcm at ≈4 A cmSOFC varya in a range frominternal 5 × 104 framework ◦ 2 700 C, yielding anode area-specific resistance of 0.1–10 ! cm . first attempt at evaluating gas-phase transport within the pore subdomain broad, this suggestions range encompasses typical values SOFC for thisperformance. and based on the data Although presented, to offer that may improve type of SOFC anode25 and is consistent with our total SOFC equation. Recent flux-based Monte C area-specific resistance at 700 ◦ C (≈0.7 ! cm2 , from the data in media suggest that such an approac the Supplementary Information, Fig. S1, which also included for both molecular and Knudsen d electrolyte and cathode resistances). The large uncertainties in porosity is higher than ∼10% (ref. this comparison reinforce the need to separate microstructure illustrates the calculation procedure. from other complicating variables, including utilization gradients, reconstructed volume of the gas po impurities, polarization, and polarization history, all of which have finite-element mesh. In practice, the at the pore–solid interface was replac been shown to be important in overall anode performance18 . 7 Chapter 2 is aimed to provide the reader with an introduction to the kinetic theory of gases as it relates to the work that follows. The mathematical and physical framework for many concepts and quantities used throughout this thesis are established in this chapter. Chapter 3 contains the background and description of the simulation method used to model mass transport in porous media. The chapter includes discussion on how model porous materials are created, ranging from simple straight capillaries, to packed spheres, to 3D reconstruction techniques used to represent real SOFC anodes. Mass transport is simulated using a Monte-Carlo-based tracer-diffusion scheme that is not novel in its implementation. The term tracer-diffusion refers to the technique of tracking the passage of molecules, or tracers, through the system. The significant contribution of this thesis comes from taking advantage of the tracer-diffusion scheme to probe the inner networks of the porous media itself. Most notably, simulations of tracer molecules following ballistic trajectories from one wall collision to the next are used to estimate the pore-size distribution from the distribution of chord lengths between each gassolid collision. These findings are presented in chapter 4. Chapter 5 provides some quantitative data describing the complex networks that make up the SOFC anode microstructure. Until recently, little has been known about the physical makeup and structure of the separate phases that make up the anode. Data pertaining to TPB density and activity can be found in this chapter, as well as volume, surface area, and connectivity information for each of the phases. Finally, chapter 6 is the concluding chapter building from the results of the previous two. Here, an alternative design for a SOFC-MEA is discussed and analyzed using the techniques developed in earlier chapters. Although the proposed design follows from earlier concepts, this thesis presents new and unique results that should provide the basis for future studies in this field. 8 9 Chapter 2 Kinetic Theory and Gas Diffusion “Here we have the kinetic theory picture of a gas—a crowd of molecules, each moving on its own independent path, entirely uncontrolled by forces from the other molecules, although its path may be abruptly altered as regards both speed and direction, whenever it collides with another molecule or strikes the boundary of the containing vessel. The molecules move so swiftly that even gravity has practically no controlling effect on their motions...that we may, without appreciable error, think of the molecules as moving in straight lines at uniform speeds...” [27] 2.1 Introduction The purpose of this chapter is to provide the reader with an introduction to the kinetic theory of gases. The mathematical and physical framework and background for the quantities I discuss in this thesis are established, all with primary emphasis on gas diffusion. For a more complete discussion on kinetic theory, the reader is referred to Kennard [32], Jeans [27], and Present [49]. Kinetic theory originated to explain and correlate the physical properties of gases and gaseous phenomena on the bases of the molecular hypothesis: that matter is composed of small discrete units known as molecules, that the molecule is the smallest quantity of a substance that retains the chemical properties of that matter (whereas an atom is the smallest portion of matter which has the property of remaining essentially intact in every chemical reaction) [32], that all molecules of a given substance are alike, and that the three states of matter differ essentially in the arrangement and state of motion of the molecules [49]. Apart from the molecular hypothesis, kinetic theory is further developed based 10 on conservation of momentum and energy, as well as the use of statistical methods. In regards to the need for statistical descriptions, Present writes: Because of the tremendous number of molecules involved in any practical sample of a gas, the actual behavior of the gas is expected to be the same as would be predicted from the average behavior of the molecules. In order to calculate this average behavior, we must assume that, after taking into account the dynamical laws and the conditions of the problem, the motions of the molecules are governed by pure chance... It is necessary to stress that kinetic theory is a statistical theory. In dealing with systems that are much too complicated to permit explicit predictions of their behavior, a representative collection of similar systems are used to develop averages and most probable behavior. 2.2 Statistical Methods and Distribution Functions Statistical information about a gas or any physical system is embodied in a so-called distribution function. This section presents several useful and very important distribution functions, including the velocity and free-path distributions. 2.2.1 The Maxwell-Boltzmann Law and the Velocity Distribution Function The integro-differential Boltzmann equation is used in this section to present a rigorous derivation of the Maxwell-Boltzmann distribution law for an equilibrium gas. The Boltzmann equation for a pure gas is simply stated, rather than explicitly derived. For a detailed and rigorous derivation and discussion of the Boltzmann equation, the reader is referred specifically to Present [49] and Bird [3], as well as any generally regarded text on kinetic theory. Consider a sample of gas containing N identical molecules that is homogeneous in physical space. Each molecule has a velocity v with components vx , vy , and vz aligned with the Cartesian axes specified by subscript. Just as physical space is defined in the Cartesian sense with position coordinates x, y, and z, velocity space is defined by the three velocity components, as shown in figure 2.1. The velocity 11 distribution function f (v) is then defined by dN = N f (v)dvx dvy dvz = N f (v)dv, (2.1) where dN is the number of molecules in the sample of gas (in all physical space) with velocity components in the range [vx , vx + dvx ], [vy , vy + dvy ], and [vz , vz + dvz ], i.e., the number of molecules contained in the velocity-space control volume in figure 2.1. Since dN and N refer to the molecules in the same volume of physical space, the above expression may be written as dn = f (v)dv, n (2.2) where the number density n has been used. The velocity distribution function is normalized Z∞ f (v)dv = −∞ Z∞ Z∞ Z∞ f (v)dvx dvy dvz = 1, (2.3) −∞ −∞ −∞ and further, f (v) can never be negative, and must have finite bounds in velocity space or tend to zero as c (molecular speed—magnitude of v—chosen to avoid possible confusion between symbols) tends to infinity. vz z dy dz dvy dvz dx dvx dr dv y r physical space dr = dxdydz x vy v velocity space dv = dvx dvydvz vx Figure 2.1. Representative control-volume elements in physical and velocity space. 12 Macroscopic properties of the gas can be obtained from various moments of the velocity distribution. For example, the average value of any molecular quantity Q can be found from the average over the molecules in an element of physical space, 1 Q= N Z QdN, such that substituting for dN from equation (2.1) gives Q= Z Qf (v)dv. (2.4) For the purposes of kinetic theory, the instantaneous state of a gas is completely specified if the distribution function for molecular velocities and positions is known throughout the gas (just as a gas flow in the classical sense would be completely described by the position, velocity, and internal state of every molecule at a particular instant—but a real gas contains such a large number of molecules that we must resort to statistical descriptions in terms of probability distribution functions). Using the velocity distribution (herein denoted simply as f for f (v)), the Boltzmann equation for a simple dilute gas is [3] ∂ ∂ ∂ (nf ) + v · (nf ) + F · (nf ) = ∂t ∂r ∂v Z∞ Z4π −∞ 0 n2 (f ? f1? − f f1 )cr σdΩdv1 . (2.5) If each term is considered to act on a phase-space volume element (i.e., dvdr), each has a physical meaning. The left-most term refers to the rate of change of the number of molecules in the element. The processes that contribute to this rate of change of molecules are represented by the additional terms in the equation. The middle term on the left-hand side describes the convection of molecules across the face of dr by the molecular velocity v. The final term on the left is also a convective term. Just as the middle term represents molecules moving into physical space carried by a velocity, the final term represents molecules moving into velocity space carried by an acceleration produced by an external force F. The velocity, or rate of change of position, is responsible for moving molecules through physical space. Similarly, the acceleration, or rate of change of velocity, is responsible for moving molecules 13 through velocity space. The delineation of physical and velocity space emphasizes the fact that v and r are treated as independent variables. The sole remaining term on the right-hand side of equation (2.5) embodies the scattering of molecules into and out of the phase-space element dvdr as a result of intermolecular collisions. The gas is considered to be dilute, in that the spacing between molecules is large compared to the molecular diameter. This assumption has two major implications: (1) a collision is consequently considered to be an instantaneous event at a fixed position in physical space, but causes a molecule to jump from one point to another in velocity space; and (2) all collision events are assumed to be binary. Since we seek the scattering of molecules into and out of the phase element, we must consider precollision and postcollision velocities of the molecules participating in collisions. In particular, we are concerned with the collision of a molecule of class v with one of class v1 , such that their respective postcollision velocities are v? and v1? . This class of collisions is responsible for scattering molecules out of the phase element, as v, v1 → v? , v1? . The inverse is true for collisions that scatter molecules into the phase element, moving them from v? , v1? → v, v1 . In the way we have defined f to denote the value of the velocity distribution at v (i.e., f = f (v)), the variables f1 , f ? , and f1? are used to denote the values of f at v1 , v? , and v1? , respectively. A molecule of class v may be chosen as a test particle moving with relative speed cr (magnitude of v − v1 ) as viewed in the frame of reference of class v1 molecules. The volume swept out in physical space by the cross section for this class of collision is given by cr σdΩ, where σdΩ is the differential cross section and dΩ is the unit solid angle centered about the post-collision relative velocity. Of the term in parenthesis on the right-hand side of equation (2.5), f ? f1? is indicative of molecules scattered into the phase space and f f1 those scattered out of phase space via binary collisions. Integration of the whole expression over the entire collision cross section, followed by integration of v1 over all velocity space, yields the total number of molecules of class v scattered into the phase element. The full Boltzmann equation given in equation (2.5) is the starting point for the derivation of the Maxwell-Boltzmann distribution law for a gas in equilibrium. The equilibrium state is defined as one in which the distribution of molecular velocities is independent of position and time. For such a state, in the absence of an external force field, the entire left-hand side of equation (2.5) must vanish. The first 14 term vanishes for any steady state, as well as the equilibrium state. The middle term involves gradients of the number density n = n(r, t) and velocity distribution function f . In equilibrium, the density must be uniform in the absence of an external field. Additionally, the equilibrium state is defined such that the velocity distribution is independent of position. Hence, the middle term vanishes. The third term on the left vanishes because there is no external force field. Under these assumptions, the Boltzmann equation of equation (2.5) reduces to 0= Z∞ Z4π −∞ 0 n2 (f ? f1? − f f1 )cr σdΩdv1 . (2.6) One possible solution of equation (2.6) requires the vanishing of the term in parenthesis for all values of v and v1 , to give the condition f ? f1? = f f1 , or equally f (v? )f (v1? ) = f (v)f (v1 ). (2.7) However, since the integrand can have either sign (i.e., positive indicates molecules are scattered into the phase element due to collisions, and negative gives a net number of molecules scattered out of the element), the vanishing of the integral does not require the condition of equation (2.7). Instead, all that is required is a mutual cancellation of the positive and negative contributions to the value of the integral from different parts of the region of integration. Therefore, the condition in equation (2.7) is sufficient for equilibrium, but the analysis to this point cannot conclude that it is necessary. Before proof of necessity via the Boltzmann H-theorem is established, it is helpful to present a moment of the Boltzmann equation by multiplying equation (2.5) by the quantity Q, and then integrating over all velocity space, just as moments of the velocity distribution are found through equation (2.4). The quantity Q relates to a single molecule and is either a constant or function of the molecular velocity v. Omitting details of the derivation, the moment equation produced from multiplication of equation (2.5) by Q and subsequent integration over velocity space is [3] ∂ ∂Q (nQ) + ∇ · (nvQ) − nF · = ∆[Q], ∂t ∂v (2.8) 15 where ∆[Q] = Z∞ Z∞ Z4π −∞ −∞ 0 n2 Q(f ? f1? − f f1 )cr σdΩdv1 dv. (2.9) Though not shown here, the moment equations can be used to obtain the monatomic gas version of the conservation equations of a continuum gas [3]. For the case of a spatially homogeneous pure gas in the absence of a field, the Boltzmann equation reduces to equation (2.6), except the left-hand side now contains n(∂f /∂t). Under the assumptions above, the number density n is constant, all spatial gradients (i.e., ∂/∂r) are zero, and the external forces F are zero. Therefore, the equation of interest simplifies to ∂f =n ∂t Z∞ Z4π (f ? f1? − f f1 )cr σdΩdv1 . (2.10) −∞ 0 The so-called Boltzmann H-theorem makes use of the Boltzmann H-function, defined as H = ln(nf ) = Z∞ f ln(nf )dv, −∞ where I have substituted Q = ln(nf ) according to equation (2.4). Using this value of Q in the simplified version of equation (2.8) corresponding to equation (2.10), gives [3] n ∂H =− ∂t 4 Z∞ Z∞ Z4π −∞ −∞ 0 ln(f ? f1? /f f1 )(f ? f1? − f f1 )cr σdΩdv1 dv. (2.11) The integrand of equation (2.11) is nonnegative for all values of the variables of integration (specifically, if f ? f1? > f f1 , then ln(f ? f1? /f f1 ) > 0 and (f ? f1? − f f1 ) > 0, and alternatively, if f ? f1? < f f1 , then ln(f ? f1? /f f1 ) < 0 and (f ? f1? − f f1 ) < 0). Hence, the quantity ∂H/∂t ≤ 0 where the equality is only true if the integrand vanishes everywhere, for all values of v and v1 . This is true only if f ? f1? = f f1 , and therefore, equation (2.7) is both necessary and sufficient for equilibrium. This equilibrium condition can be written as ln f ? + ln f1? = ln f + ln f1 . This equation shows that ln f is a collisional invariant at the equilibrium state, i.e., it is a conserved 16 quantity for binary molecular collisions. The only other collisional invariants for a simple monatomic gas are the mass m, momentum mv, and kinetic energy mc2 /2, set forth by conservation laws. The quantity ln f must therefore be a superposition of these possibilities, ln f = Amc2 /2 + B · mv + C, (2.12) where A, B, and C are constants. The velocity can be written in terms of the mean or free stream or mass velocity v0 and the thermal or fluctuating or random velocity v0 , as v = v0 + v0 . The square of the speed becomes c2 = v · v = c20 + 2v0 · v0 + c02 (note again that I use the convention of v to represent a velocity vector, and the scalar speed c is the norm of that vector). This result and the expanded relation for v can be substituted into equation (2.12) to give ln f = Amc02 /2 + m(Av0 + B) · v0 + Amc20 /2 + B · mv0 + C. By definition, the random or fluctuating velocities of an equilibrium gas must have no preferred direction so that the distribution is isotropic. This requires the coefficient of the second term to be zero, such that B = −Av0 , and ln f = Amc02 /2 − Amc20 /2 + C, which leads to an expression for the velocity distribution f , f = exp(Amc02 /2 − Amc20 /2 + C) = exp(−Amc20 /2 + C) exp(Amc02 /2). As defined after equation (2.3), the distribution f must be bounded, requiring Am/2 < 0 (note that the first term is simply a finite constant, and the term we must bound is the second term involving c02 ). For convenience, let Am/2 = −β 2 so that f = exp(β 2 c20 + C) exp(−β 2 c02 ). 17 The constant exp(β 2 c20 + C) can be determined from the normalization condition, Z∞ f dv = exp(β 2 c20 + C) Z∞ exp(−β 2 c02 )dv0 = 1, −∞ −∞ where Z∞ Z∞ Z∞ Z∞ exp(−β c )dv = 2 02 0 exp[−β 2 (vx02 + vy02 + vz02 )]dvx0 dvy0 dvz0 = π 3/2 /β 3 , −∞ −∞ −∞ −∞ so exp(β 2 c20 + C) = β 3 /π 3/2 , and f (v) = fM = (β 3 /π 3/2 ) exp(−β 2 c02 ), (2.13) is the equilibrium or Maxwellian distribution of thermal velocities (where c02 = v0 · v0 ). The constant β can be related to the temperature of the gas through the average kinetic energy 1 3 kB T = mc02 2 2 Z∞ 1 = m c02 fM dv0 2 −∞ β3 1 = m 3/2 2 π Z∞ c02 exp(−β 2 c02 )dv0 , −∞ where Z∞ c02 exp(−β 2 c02 )dv0 −∞ = Z∞ Z∞ Z∞ −∞ −∞ −∞ Z∞ =3 exp(−β 2 vy02 )dvy0 −∞ =3   (vx02 + vy02 + vz02 ) exp −β 2 (vx02 + vy02 + vz02 ) dvx0 dvy0 dvz0  π 1/2 β  π 1/2 β  Z∞ −∞ 1/2 π 2β 3 exp(−β 2 vz02 )dvz0  Z∞ −∞ . Therefore, 3 1 3 kB T = m 2 , 2 4 β vx02 exp(−β 2 vx02 )dvx0 18 giving β2 = m . 2kB T (2.14) Since I have assumed an equilibrium gas, the translational and thermodynamic temperatures are equal, so there is no need to distinguish between the two. The distribution of molecular velocities can be used to derive the distribution of molecular speeds. It follows from equation (2.13) and equation (2.14) that the equilibrium distribution of molecular speeds is given by f (c0 ) = (4/π 1/2 )β 3 c02 exp(−β 2 c02 ). (2.15) The most probable thermal speed c0m and mean thermal speed c0 can be found from the above distribution as df (c0 ) = 0 → c0m = 1/β, dc0 c0 =c0 (2.16) c0 f (c0 )dc0 = 2/(π 1/2 β) = (2/π 1/2 )c0m . (2.17) m c0 = Z∞ 0 PDF for molecular speeds, f Hc'L The distributions and other quantities derived above will be used extensively in this work. 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.5 1.0 1.5 Βc' 2.0 2.5 3.0 Figure 2.2. Plot of the probability distribution function for molecular speeds, f (c0 ), from equation (2.15). The maximum or most probable speed occurs at βc0 = 1 as given by equation (2.16). 2.2.1.1 Sampling Molecular Speeds from the Equilibrium Distribution Statistical simulation generally requires sampling values from a prescribed distribution through the use of random numbers. The acceptance-rejection method can be used to sample speeds from f (c0 ). 19 While the method involves a repetitive process which may require the evaluation of a large number of functions and many random fractions, it is quite efficient for this case. In order to make direct use of a random fraction rf distributed uniformly over the range [0, 1], the distribution function is normalized, F= f (c0 ) , f (c0m ) where f (c0m ) corresponds to the maximum value of the distribution function on the interval [0, ∞). This normalization requires F to be on [0, 1]. A value of c0 is chosen at random, on the basis of c0 being uniformly distributed between its limits. However, since the upper limit is ∞, an arbitrary cutoff value of 3/β is used. The fraction of values of c0 lying outside the range [0, 3/β] is 0.00043985. Hence, c0 is sampled from c0 = 3rf /β, which leads to F= f (c0 ) c02 exp(−β 2 c02 ) = 02 = 9rf2 exp(1 − 9rf2 ), 0 f (cm ) cm exp(−β 2 c02 ) m where the two instances of rf2 indicate the same random fraction (as opposed to generating two separate random fractions, one for each instance). A second random fraction rf 2 is then generated and, if F > rf 2 , the value of c0 is accepted. If, on the other hand, F < rf 2 , the value of c0 is rejected and the process is repeated (i.e., generate new rf for c0 , compute F, then generate another rf 2 and compare to F) until a value of c0 is accepted. 2.2.2 Scattering of Molecules from a Surface At equilibrium, the flux of molecules to a surface (incident molecules) must be equal to the flux of molecules leaving the surface (emitted molecules). Hence, there is no net mass flow, just as there can be no net energy or momentum transfer. Because the rates of adsorption and desorption are equal, the coverage is constant. Equilibrium also requires the gas and surface temperatures to be equal. It follows that the angular distributions of the incident and emitted molecules must be identical (and hence, no flow directionality). The conservation laws do not mean that when a molecule hits a surface, resides there for some time, and subsequently desorbs it enters the gas phase with exactly the 20 z rsinθ dφ rsinθ rdθ dθ θ r φ dφ y x Figure 2.3. Spherical polar coordinates in physical space. In velocity space, in a frame of reference moving with the stream velocity, r → c0 and the x, y, and z axes become vx0 , vy0 , and vz0 , respectively. same energy it had before it became adsorbed [35]. The beauty of thermodynamics is that, even though no individual molecule returns to the gas phase with exactly the same energy and direction of travel, when averaged over all molecules in the system, the energy and propagation directions do not vary in time for a system at equilibrium when averaged over a sufficiently long time. A good discussion on the arguments and issues related to ad/desorption phenomena and angular distributions can be found in Kolasinski [35] and Greenwood [21]. Essentially, at equilibrium, since the angular distribution of molecules incident on the surface follows a cosine distribution, based on simple geometric arguments, the angular distribution of the desorbing molecules must also follow a cosine distribution. This is often referred to as the Knudsen cosine law, though Kolasinski notes “that Knudsen’s proposition should attain the status of a ‘law’ is remarkable because (1) it was incorrectly derived and (2) experiments of the type he [Knudsen] performed to prove it should exhibit deviations from a [cosine] distribution.” Nevertheless, subsequent work (see references in Kolasinski, §3.8) established that the second law of thermodynamics demands the validity of the cosine law for the angular distribution of molecules that leave a surface. For a system of spherical-polar coordinates, as shown in figure 2.3, the volume of the velocity space 21 n dθ solid angle dΩ = sinθ dθ dφ θ dA t2 φ dφ t1 Figure 2.4. Molecules crossing a surface dA in direction dΩ with speed c0 must lie in the volume element c0 cos θdAdt. The solid angle dΩ makes an angle θ with the surface normal n. The two vectors tangent to the surface are t1 and t2 . element pictured in figure 2.1 is dv = c02 sin θdθdφdc0 . (2.18) Combining equation (2.2), equation (2.13), and equation (2.18), the fraction of molecules with speed between c0 and c0 +dc0 , which make an angle between θ and θ +dθ with the polar direction, and between φ and φ+dφ with the azimuth (i.e., traveling in the direction defined by the solid angle dΩ = sin θdθdφ) can be written as dn/n = fM c02 sin θdθdφdc0 = (β 3 /π 3/2 )c02 exp(−β 2 c02 ) sin θdθdφdc0 . I can write this more simply as dn/n = f (c0 )dc0 dΩ , 4π where f (c0 ) is given in equation (2.15). Consider now the picture in figure 2.4. The flux of molecules dJ = dn/(dAdt) (molecules/m2 ·s) crossing surface element dA in time dt is dJ = nc0 cos θf (c0 )dc0 dΩ . 4π 22 Figure 2.5. Three-dimensional polar plot of the angular distribution functions. The right image shows a cutaway through the vertical midplane of the left image. The rotational symmetry is due to the uniform distribution of the azimuthal angle over 2π as given by g(φ) = 1/2π. The lopsided shape evident in the cutaway image is the shape of the polar angle distribution function h(θ) = 2 cos θ sin θ. This expression can be used to obtain several important quantities and distributions. The first is the total flux of molecules crossing the surface from one side (relative to the orientation of the surface normal), J ± =n ± Z∞ 0 c f (c )dc 0 0 0 Zπ/2 1 cos θ sin θdθ 2 0 Z2π n± c0 dφ = , 2π 4 (2.19) 0 where the + and − superscripts serve only to signify the flux from each side of the surface and the dependence on the concentration n± on that side. For zero net flux across a surface, J + = J − and n = n+ = n− . 2.2.2.1 Sampling from the Angular Probability Distributions The equation for J + and J − reveals the distribution functions for the polar angle θ and azimuthal angle φ. These distributions give the angular probability of a molecule incident on or leaving from a surface. The normalized distribution functions are g(φ) = 1/(2π), h(θ) = 2 cos θ sin θ. These distributions are illustrated in the plots of figure 2.5 and figure 2.6. For a set of random fractions rf , uniformly distributed on [0, 1], the inverse-cumulative method PDF for the polar angle at a surface, hHΘL 23 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.5 1.0 1.5 Θ HradL Figure 2.6. Plot of the probability distribution function for the polar angle, h(θ) = 2 cos θ sin θ, for molecules incident on or leaving a surface. The most probable angle is π/4. gives a simple way to sample each of the angles, rf = Zφ 1 g(φ)dφ = 2π 0 Zφ dφ = 0 φ → φ = 2πrf , 2π and rf = Zθ h(θ)dθ = 2 0 Zθ cos θ sin θdθ = 2 0 sin Z θ 0 1/2 sin θd(sin θ) = sin2 θ → sin θ = rf . Likewise, since sin2 θ + cos2 θ = 1 and because 1 − rf is no different from rf , the variable θ can be 1/2 sampled through cos θ = rf . Successive values of rf are used to obtain values for each of the angles. Sampling from the angular distributions is tremendously important each time a molecule collides with a solid surface. A new speed and new direction must be randomly generated for each molecule at every collision. 2.2.2.2 Sampling from the Speed Distribution at a Surface The distribution of molecular speeds at a surface is not prescribed by the Maxwell distribution f (c0 ). Instead, equation (2.19) reveals that the distribution of speeds is actually of the form c0 f (c0 ) such that the normalized distribution is f s (c0 ) = c0 f (c0 ) = 2β 4 c03 exp(−β 2 c02 ), c0 (2.20) PDFs for molecular speeds, f Hc'L or f s Hc'L 24 1.0 0.8 f s Hc'L 0.6 0.4 f Hc'L 0.2 0.0 0.0 0.5 1.0 1.5 Βc' 2.0 2.5 3.0 Figure 2.7. Comparison of the probability distribution functions for a Maxwell-Boltzmann speed distribution, f (c0 ) from equation (2.15), and the distribution of speeds across a surface, f s (c0 ) from equation (2.20). where the superscript s indicates surface, f (c0 ) is given by equation (2.15), and the mean speed c0 follows from equation (2.17). The most probable and mean speed for molecules striking the surface are, respectively, 1/2 c0s /β = (3/2)1/2 c0m , m = (3/2) c0s = Z∞ c0 f s (c0 )dc0 = 3π 1/2 /(4β) = (3π/8)c0 = (3π/8)1/2 c0s m. 0 Figure 2.7 gives a comparison of f (c0 ) and f s (c0 ). The overall shape of the distribution remains fairly consistent, but the distribution of speeds at the surface is shifted toward higher speeds. As is the case for the Maxwell speed distribution, sampling of speeds at a surface can be done relatively efficiently using the acceptance-rejection method described previously. 2.2.3 Distribution of Free Paths One of the most useful concepts of kinetic theory is that of the mean free path [49] and the probability of a free path of given length. As with the other distribution functions, knowledge of the distribution of free paths will be particularly useful in modeling gas transport within a porous substrate. Let f (`) represent the probability that a molecule has a free-path length at least equal to `, i.e., that no collision occurs in the distance `. After the molecule has traveled a distance `, the chance of collision within a further distance d` is proportional to d` and can be written simply as αd` where α is 25 a proportionality constant, to be determined later. In general, α will depend on the number density of gas molecules, the molecular size, the speed, and the pressure. Since αd` is the probability of collision between ` and ` + d`, then the chance of no collision occurring within the same distance is simply 1 − αd`. The total probability that no collision occurs in the distance ` + d` is given by the product f (`)(1 − αd`). By definition, this must be equal to f (` + d`) which can be written in differential form as f (`) + df (`) d`, d` so that f (`)(1 − αd`) = f (` + d`) = f (`) + df (`) d`, d` to give the differential equation, df (`) = −αf (`). d` As defined, the parameter α > 0 and obviously f (0) = 1. Imposing these conditions leads to a solution for the distribution of free paths, f (`) = e−α` . The probability of a free path terminating between ` and ` + d` is given by f (`)d` = f (`) − f (` + d`) = − df (`) d` = αf (`)d` = αe−α` d`. d` The mean free path λ can be calculated as λ=`= Z∞ 0 `f (`)d` = α Z∞ 0 `e−α` d` = 1/α. 26 Therefore, the distribution of free paths is given by the exponential function, f (`) = e−`/λ . A plot of f (`) is provided in figure 2.8. PDF for free paths, f H{L 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 10 {Λ Figure 2.8. Plot of the probability distribution function for free paths, f (`) = e−`/λ . The exponential shows that free paths that are many times greater than the mean free path λ are rare. 2.2.3.1 Sampling from the Free-Path Distribution Again, as in the sampling of molecular speeds, I assume the availability of a successive set of random fractions rf that are uniformly distributed on the interval [0, 1]. Free paths can be sampled without the use of the acceptance-rejection method. Instead, the inverse-cumulative method allows for quick sampling since the cumulative distribution can be inverted to solve for the free path in terms of rf (this is not the case for the distribution of speeds in the previous section). For the free-path distribution, a free-path length may be selected by allowing rf = Z` 0 1 f (`)d` = λ Z` 0 e−`/λ d` = 1 − e−`/λ . Noting that 1 − rf is the same as rf , a free path ` can be randomly generated as ` = −λ ln rf . 27 It is very important to note that the chance of a collision is independent of the dynamical history of the molecule. 2.3 Diffusion The diffusion coefficient is the proportionality constant relating the gradient in species concentration to mass flux. It is one of the three main transport coefficients that appear as parameters in the macroscopic conservation equations for mass, momentum, and energy [30]. The rigorous mathematical derivation of transport properties is very complex [4, 9], and it is not my intention to derive them. I am simply presenting the important results used extensively in the following chapters of this thesis. 2.3.1 Normal Diffusion in Capillaries It is convenient to first consider the normal one-dimensional diffusion of gases through capillary tubes. For binary diffusion, the usual equation for the molecular flux of each species can be written as [49, 4] Ji = −D12 ∂ni + xi J, ∂z i = 1 or 2, (2.21) where n = n1 + n2 is the total molecular concentration, xi = ni /n is the mole fraction of species i, and J = J1 + J2 is the net flux of molecules. It is assumed that the capillary axis is aligned with the coordinate z. The parameter D12 is the normal diffusion coefficient and is the focus of this section. Deriving an expression for D12 requires a rigorous, complex, and specialized theory in statistical mechanics and molecular interactions [30, 9]. The detailed derivation begins with the Boltzmann equation of equation (2.5) and is solved using a method of successive approximation known as ChapmanEnskog theory [9]. The solution method is not presented here, though its elegance ought to be noted; I am concerned more with the results. For a binary gas mixture, Chapman-Enskog theory provides p p 3 T 3 /m 3 2πkB T /m12 3 2πkB 12 = , D12 = 16 pπσ 2 Ω(1,1)? 16 nπσ 2 Ω(1,1)? 12 12 12 12 where m12 = m1 m2 /(m1 + m2 ) is the reduced or inertial molecular mass, σ12 is a length scale corre- 28 sponding to the interaction of a molecule of type 1 and one of type 2 (often called the collision diameter), (1,1)? and Ω12 is the so-called collision integral, which depends on temperature and the interaction poten- tial between molecules. The pressure has been related to the molecular concentration through the ideal gas law, p = nkB T. For hard spheres, when 1 and 2 are indistinguishable, the collision integral is defined to equal unity (1,1)? and σ12 = d, where d is the sphere diameter. Deviations of Ω12 from unity are a measure of the deviations of the intermolecular forces from forces between rigid spheres [18]. For indistinguishable rigid spheres, the Chapman-Enskog result reduces to 3 Dii = √ 16 2 p 3 T 3 /πm 8kB 3 = √ pd2 16 2 p 8kB T /πm , nd2 i = 1 or 2, (2.22) where m1 = m2 = m and Dii is called the self-diffusion coefficient. Let me now reiterate some useful results derived in the previous sections. I assumed that molecules in the gas consist of rigid, hard spheres of mass m and diameter d that are constantly in motion. Collisions between molecules are binary and instantaneous, and the molecules travel in ballistic trajectories from one collision to the next. The mean speed of a molecule from equation (2.17) together with equation (2.16), and equation (2.14) is c0 =  8kB T πm 1/2 . This can be substituted into equation (2.22) to give Dii = 3 c0 √ , 16 2 nd2 i = 1 or 2, where the hard-sphere diameter can be related to the gas molecule mean free path by kB T 1 λ= √ =√ . 2πpd2 2πnd2 29 Using this result, another form of the self-diffusion coefficient is Dii = 2.3.2 3 πλc0 , 16 i = 1 or 2. Diffusion as a Random Walk The general diffusion equation can be derived from the Boltzmann moment equation of equation (2.8) and equation (2.9) letting Q = m. Noting that m is a collision invariant (i.e., mass is always conserved in collisions), the collision integral of equation (2.9) vanishes. The appropriate averaged quantities are easily calculated from equation (2.4) to give ∂ (nm) + ∇ · (nmv) − 0 = 0. ∂t Since the mass m is constant, it can be divided out from both remaining terms. In addition, note that the molecular flux J = nv, so that ∂n + ∇ · J = 0. ∂t Generalizing equation (2.21) to three dimensions, and assuming no net flux with uniform diffusivity over all phase space leads to the diffusion equation: ∂n − D∇ · ∇n = 0, ∂t or rearranging, ∂n = D∇2 n. ∂t (2.23) The theory discussed above is based on the concept of a local concentration of gas molecules, requiring a large number of molecules to average out randomness in the instantaneous positions. However, in efforts to relate molecular-level phenomena to observed rates of transport, it is more natural to focus on the motion of a single molecule, or Brownian particle. The position of such a molecule can only be described in terms of probabilities, so that the calculation becomes a stochastic problem. The derivation that follows concerns the motion of a single molecule executing a random walk. 30 The molecule is assumed to move independently of other molecules, and chemical reactions are not considered. Denote the probability density of finding the molecule at position r at time t as P(r, t). The molecule moves in discrete steps of length ` in random directions. The time required for such a step is τ . All directions are equally available, so the probability that a given step will end somewhere on a spherical surface of radius ` centered about the current position is unity. Let R be the vector from the current position to the next position a distance ` away in a random direction. The probability density describing R is S(R) = δ(|R| − `) , 4π`2 where δ(|R| − `) is the Dirac delta function indicating that jumps are only allowed to be of length `. This function is normalized so that integration over all possible R gives unity: Z SdR = Z2πZπ Z∞ 0 0 δ(r − `) 2 1 r sin θdrdθdφ = 2 4π`2 ` 0 Z∞ 0 δ(r − `)r2 dr = 1. (2.24) The probability of finding the molecule at position r at time t + τ , i.e., P(r, t + τ ), is the product of P(r − R, t) and the probability describing the the jump S(R), integrated over all possible R. Mathematically, this translates to P(r, t + τ ) = Z P(r − R, t)S(R)dR. (2.25) After a large number of jumps, and assuming ` and τ are small compared to the length and timescales of interest, P(r, t) can be approximated as a continuous function so that the probabilities in equation (2.25) can be written as Taylor expansions about position r and time t. The expansion in time of the left-hand side gives P(r, t + τ ) = P(r, t) + ∂P τ + O(τ 2 ), ∂t and the expansion for the term in the integral gives 1 P(r − R, t) = P(r, t) − R · ∇P + RR : ∇∇P + O(R3 ), 2 31 where all derivatives are evaluated at r and t. Substituting the leading terms of these expansions into equation (2.25) produces ∂P P(r, t) + τ = P(r, t) ∂t Z SdR − ∇P · Z 1 RSdR + ∇∇P : 2 Z RRSdR. There are three integral terms on the right-hand side of the equation. The first R (2.26) SdR = 1 by equa- tion (2.24). The second represents the average value of the displacement vector R. Since the vector has no directional bias, this average is expected to vanish: Z Z2πZπ Z∞ δ(r − `) 2 RSdR = {r sin θ cos φ, r sin θ sin φ, r cos θ} r sin θdrdθdφ 4π`2 0 0 0 Z2πZπ Z∞ 1 = {sin2 θ cos φ, sin2 θ sin φ, sin θ cos θ}δ(r − `)r3 drdθdφ 4π`2 0 0 0 Zπ Z∞ 1 = {0, 0, 2π sin θ cos θ}δ(r − `)r3 drdθ 4π`2 0 0 Z∞ 1 {0, 0, 0}δ(r − `)r3 dr = 0. = 4π`2 0 The third integral gives the average value for the dyadic product RR: Z RRSdR   sin θ cos φ Z2πZπ Z∞    sin3 θ sin φ cos φ =   0 0 0  sin2 θ cos θ cos φ  3  π sin θ Zπ Z∞   1  = 0  2 4π`  0 0  0  3 1 = 4π`2 sin θ sin φ cos φ 3 2 0  4π/3 Z∞    0 4π/3   0  0 0 sin3 θ sin2 φ sin2 θ cos θ sin φ 0 0 π sin3 θ 0 0 2π sin θ cos2 θ  0 0 4π/3  sin θ cos θ cos φ    δ(r − `) 4 sin2 θ cos θ sin φ   4π`2 r drdθdφ   sin θ cos2 θ  2     δ(r − `)r4 drdθ      2   δ(r − `)r4 dr = ` δ,  3   32 where δ is the identity matrix. Substituting the integral results into equation (2.26) yields P(r, t) + ∂P `2 τ = P(r, t) − 0 + ∇∇P : δ. ∂t 6 Simplifying the above expression and using the identity ∇∇P : δ = ∇2 P gives a diffusion-type equation, ∂P `2 2 = ∇ P. ∂t 6τ This result is analogous to equation (2.23). While the equation above was derived for a single molecule, it applies to all molecules individually. Summing each side over all molecules changes the probability density P to a number density or concentration. Consequently, this is exactly what appears in equation (2.23). Therefore, the self-diffusivity D is related to the jump distance and time as D = `2 /6τ so that ∂P = D∇2 P. ∂t (2.27) In the absence of any directional bias, the probability density P is expected to be spherically symmetric and therefore, P = P(r, t) only. The solution to equation (2.27) must satisfy the normalization condition that the probability of finding the molecule somewhere in space is unity. Also, one expects P → 0 as r → ∞. The PDE is readily solved under the stated conditions to give P(r, t) = 2 1 e−r /4Dt . 8(πDt)3/2 This solution enables us to relate D to measured molecule positions (e.g., as obtained in experiments or simulations). By observing and recording the results for a large number of random-walking molecules, certain average properties of the position can be evaluated. The mean displacement of a molecule follows from hri = Z rPdV = 0, as a result of spherical symmetry. This result contains no information about D. The mean-square 33 displacement is subsequently hr i = 2 Z Z∞ 4π r PdV = 8(πDt)3/2 2 2 r4 e−r /4Dt dr = 6Dt. 0 By measuring the displacements of a large number of molecules over a fixed time interval t, the diffusivity D can be calculated from the average of the squares of these displacements hr2 i as D= hr2 i h|r(t) − r(0)|2 i = . 6t 6t (2.28) This is an incredibly useful result for the evaluation of diffusion coefficients from experimental observations or molecular simulations! This result is heavily relied on in this thesis. 2.3.3 Diffusion in Porous Media and Free-Molecule Diffusion in Capillaries Normal diffusion in porous media is analogous to diffusion in tubes, and can be described by the same equations if suitable parameters characterizing the porous medium are inserted [18]. Typically, an effective diffusion coefficient is substituted in place of D12 to give eff Ji = −D12 ∂ni + xi J, ∂z i = 1 or 2, (2.29) eff where Ji and J now refer to a flux based on unit area of porous medium and D12 is the effective diffusivity, proportional to D12 . The constant of proportionality accounts for the internal structure of the porous medium, and includes the porosity  and tortuosity τ as [18, 4] eff Dij = Dij /τ, i, j = 1 or 2. Like the porosity, the tortuosity is considered to be a macroscopic parameter, characteristic of the porous material itself. It is an all-inclusive measure of the complexity of the pore network, accounting for the twisting and turning of pores, necking, and so on. The benchmark case is a straight, cylindrical pore of constant cross section, which has τ = 1. All other systems have a tortuosity greater than this. Free-molecule or Knudsen diffusion in porous media is drastically different from normal diffusion. 34 The transport of any individual species is unaffected by the presence of others. Therefore, gas-gas collisions are negligible, and the flow is dominated by collisions between gas molecules and the pore surfaces. The importance of the last statement cannot be overstated. The only influence on the flow is the structure of the porous medium. Understanding the internal framework of the porous medium is necessary if Knudsen flow is to be optimized. Since gas-gas interactions can be ignored, the Knudsen diffusivity DiK becomes identical with the permeability coefficient Ki , which is defined by the expression [18] Ji = − Ki ∂pi ∂ni = −Ki , kB T ∂z ∂z where pi is the partial pressure of species i. In place of equation (2.29), Knudsen flow gives Ji = −DiK ∂ni , ∂z (2.30) which has no convective term. The Knudsen diffusivity is generally given by the expression [49, 66, 4] DiK = 2 R 3  8kB T πm 1/2 = 2 0 Rc , 3 (2.31) eff where R has dimensions of length. Just as  and τ appearing in Dij are characteristic of the porous medium, R is a constant of the medium. Knudsen diffusion in capillaries is well covered in the literature [34, 11, 32, 49, 15, 3]. For the simple case of diffusion in a tube, the parameter R is the tube radius. Because of this, many people have extended the assumption that R is equal to the mean pore radius for a wide range of porous materials. While this seems perfectly reasonable for the ideal case of tubelike pores, most engineering applications involve highly complex materials with both macropores and micropores of highly unstructured shapes. It seems insufficient to describe such a material using only a mean pore radius. Understanding more about the pore-size distribution is a necessary step to understanding the diffusive transport within the porous medium. In fact, all of chapter 4 of this thesis is devoted specifically to investigating the pore-size distribution function. The length scale R is further discussed in that chapter. It is interesting to note the pressure dependence of the normal self-diffusion coefficient Dii and the 35 Knudsen diffusivity DiK . The former is inversely proportional to pressure but the latter is independent of pressure. In the transitional flow limit, or intermediate pressure region, the diffusion coefficient is generally calculated using the formula of Bosanquet [6] DiB =  1 1 + Dii DiK −1 . (2.32) This formula can be used for all pressures since the appropriate diffusivity term will dominate—Dii for high pressures, and DiK for low pressures. For all pressures, the effective diffusivity for single species flow in porous media follows as eff DiB = DiB /τ. Note that the self-diffusion coefficient of equation (2.28) was derived by assuming a molecule had all physical space available to it through which to move. However, only a fraction of this space is available in a porous sample, and the effective diffusivity is D. This can be related to DiB using the equation above, D = DiB /τ, so that D/DiB = /τ. (2.33) The above ratio between the effective diffusion coefficient and the coefficient of diffusion in the free gas is a dimensionless constant smaller than unity, which is characteristic of the porous material—denote this constant as the diffusibility. This ratio is used in the following chapters. 36 37 Chapter 3 Mass Transport in Porous Structures “The influence on [diffusion] of not only porosity, but also particle shape, was clearly shown and there can therefore be no unique relationship for all materials, as often supposed, between diffusion rates and porosity.” [14] 3.1 Introduction Gas transport in porous solids has received a great deal of attention because of its importance to a wide range of engineering applications. Processes that intimately depend on diffusion in porous media span from catalytic reactors, to nutrient delivery in biological systems, to methane recovery in coal beds. Because of the complexity of the diffusion processes that take place in porous media, there is an enormous body of literature dealing with experimental determination and the theoretical prediction of effective diffusivities. These complexities involve bulk, Knudsen, and surface diffusion and in realistic systems, takes place in the presence of adsorption, reaction, temperature gradients, and convective transport [33]. In porous catalysts, reactions at pore walls are often transport limited by the delivery of reactants to (and removal of products from) active surfaces. In general, reactant delivery and product removal from reactive sites are frequently hindered by transport limitations due to the complex three-dimensionality of the pore network. These limitations control reaction rates as well as surface coverages of adsorbing species. Despite the importance this has on a broad spectrum of technologies, it is of particular importance to solid-oxide fuel cells (SOFCs). Approaches to understand transport in porous media are as diverse as the problems they attempt 38 to solve, from early theory and experiments by Derjaguin [16], Hoogschagen [23], Millington [45], Currie [13, 14], and Weissberg [62], to the “Dusty-Gas” model of Evans III et al. [18, 19], Mason et al. [41], and Mason and Malinauskas [42], Monte-Carlo type simulations apparently initiated by Evans et al. [17] and Abbasi et al. [1] on up to the more recent Zalc et al. [67, 68], and the emergence of the lattice Boltzmann method used by Chiu et al. [10], Joshi et al. [29, 28], and Asinari et al. [2]. There are several mechanisms governing the complex nature of mass transport in porous solids, including molecular diffusion, Knudsen diffusion, viscous and slip flow, as well as surface diffusion along pore walls [66]. An accurate understanding of mass transport in a SOFC involves several challenges [10]: • mass transfer by diffusion occurs at high temperature (600–1000 ◦ C) and through submicron pore sizes with Knudsen number (Kn = λ/d, where λ is the gas mean free path and d is the characteristic dimension of a pore) of order unity or greater—therefore, continuum theory is no longer valid; • depending on the fuel, a large number of gaseous species may be present and Fick’s law (only valid for binary diffusion) does not apply; • gas molecules can adsorb/desorb on the solid pore surface and undergo electrochemical and reformation reactions; • the complex pore geometry is difficult to characterize and realistically represent in predictive models. A comprehensive model for SOFC mass transfer should include all these effects, and using such a model as a pore-level design tool can have a significant impact on optimizing the SOFC microstructure. Until recently, pore-scale models of SOFC mass transfer were not very practical, due in part to a lack of quantitative data describing the SOFC electrode microstructure in detail. A majority of past investigations consider a simplified one-dimensional (1D) model of the SOFC electrode, and condense the description of the electrode microstructure into a few parameters, most notably the porosity, mean pore size, and tortuosity. The mean pore size may not always be a good descriptor of the pore network and corresponding pore-size distribution (PSD). This is discussed in much more detail in chapter 4. Further, the tortuosity is difficult to measure, and presumably serves as an all-inclusive benchmark of 39 the complexity (i.e., the twisting and turning, necking, connectivity, etc.) of the pore network. It is often used as a fitting parameter to correlate model predictions with experimental data, without regard to details of the microstructure (mostly because these details were not known). The Dusty-Gas model, for example, is commonly used to simulate 1D multicomponent mass transfer outside of the continuum flow regime. Once calibrated, these simplified models have proved to be very useful in predicting SOFC performance and studying the impact of a few input parameters. Unfortunately, these models lack the capability to analyze and optimize SOFC electrode microstructures and are limited to simple geometries. Recent advancements in high-resolution imaging techniques (see Wilson et al. [64] and the references therein) have paved the way for pore-scale mass-transport models of SOFC electrodes. Not only have these recent developments opened the door for studying mass transport in complex porous SOFC microstructures, but details of the network morphology can be exposed (as we shall see in chapter 5). In particular, with regard to modeling mass transport, Monte-Carlo-based tracer methods are quite suitable for the analysis of multicomponent, noncontinuum mass transfer within a complex geometry like porous SOFC electrodes. This chapter reviews some of the numerical methods used to model mass transport, focuses on recent work related to Monte-Carlo-based methods, describes the simulation procedure used in this work, and presents results for a SOFC anode microstructure analyzed in two very different ways: (1) as a particle-based structure represented by a collection of various size spheres; and (2) as a voxel-based mesh of a digitally reconstructed anode. 3.2 Background Traditional computational fluid dynamics (CFD) models for SOFCs are suitable for the macroscopic analysis of flow or heat transfer in gas channels or SOFC stacks (e.g., see figure 3.1). However, the application of these models to particle scales (submicron) is problematic due to noncontinuum effects, complex geometries with various length scales, and in general, a large number of diffusing species. There are several modeling approaches for multicomponent mass diffusion under these conditions. One of the most fundamental techniques is molecular dynamics (MD), but it is also the most computationally expensive. Monte-Carlo (MC) methods are less expensive and work best for free-molecular Knowledge and understanding at certain scales informs design and optimization at other scales Naval Research Laboratory Ansys, Inc. 40 California Institute of Technology Colorado School of Mines Figure 3.1. Knowledge and understanding at certain scales informs design and optimization at other scales. Image courtesy of R. J. Kee. RTC_intro_4/06.p3 or Knudsen flow. In this limit, gas-gas collisions can be neglected, and gas behavior is only influenced by gas-surface interactions. The Dusty-Gas model works very well for 1D-multicomponent diffusion, but extension to complex geometries is a major drawback. Lattice Boltzmann methods (LBM) are This RTC has significantly advanced the state-ofthe-art inand modeling electrochemical processes growing in popularity, are well-suited to this type of problem (for recent work and review, see Chiu Naval Research Laboratory Ansys, Inc. California Institute of Technology Colorado School of Mines et al. [10]). Monte-Carlo tracer algorithms provide an excellent middle ground between MD (and even Major accomplishments and new modeling tools more complex direct-simulation MCfield approaches) and electronic-structure LBM. They are not nearly as computationally • Incorporate electric effects into computation • Build reaction mechanisms based on atomic-level understanding demanding •!Predict as MD, yet they retain more of the underpinning physics han LBM methods. Another homogeneous chemistry of molecular-weight growth • Closely couple elementary thermal and electrochemical reactions advantage is intrinsic elementary to the method itself, in that molecule paths can transport provide insight into the • Couple chemistry withtracing complex porous-media • Distributed charge-transfer within porous electrode structures influence of•!Represent the internal pore architecture. cermet electrode structures at the particle level • Develop major new general-purpose CANTERA software Monte-Carlo methods for predicting gas transport in porous mediapowerful are essentially the same as those • Couple CANTERA and FLUENT as an extraordinarily capability Demonstrate new capabilities bytubes collaboration with fuel-cell developers used to calculate molecular flow rates through of various shapes [15]. These geometries provide an • Assist the development of new cell architectures •!Assist theand development of fullbecause systems excellent starting point model validation they have been the subject of various theoretical analyses [34, 11, 48]. Extension of the MC method to more complex geometries of modeled porous media was initiated by Evans et al. [17] and Abbasi et al. [1]. In the following years, the same MC techniques were adopted by Smith and Huizenga [55] and Huizenga and Smith [24] for systems of packed RTC_intro_4/06.p4 spheres, and favorably compared to experimental results obtained by the same authors. There is a large body of work between these pioneering efforts and recent years, all utilizing a MC approach. For some of the most notable, the reader is referred to Burganos and Sotirchos [8], Reyes and Iglesia [51, 50], 41 Tomadakis and Sotirchos [58], Meyerhoff et al. [43], Riley et al. [52], Transvalidou and Sotirchos [59], Trinh et al. [60], and Zalc et al. [67, 68]. 3.3 Simulation Procedure This section first describes the Monte-Carlo tracer algorithm for diffusing molecules in a cylindrical tube, then extends the method to more complex packed-sphere models and three-dimensional (3D) voxel-based mesh geometries representing real SOFC electrodes. The fundamentals of the MC method are not new, but modifying it to be used in voxel-based geometries is original to this work. 3.3.1 Tracer Diffusion in Tubes Consider an infinite cylinder of radius a centered about the z-axis, so that any point on the surface of the cylinder is given by x = a cos θ, y = a sin θ, z. The polar angle θ is measured counterclockwise from the positive x-axis. The surface of the cylinder is assumed to be a purely diffuse reflector of molecules incident on it. It is a solid surface representing the pore wall (i.e., in this case the pore happens to be a simple cylinder). In the free-molecule limit, gas molecules do not interact with one another, so each may be treated independently. Under this assumption, a molecule, or tracer, is placed at a random position inside the cylinder (i.e., inside the pore). The random position is chosen uniformly from all those available so that doing this for a large number of tracers essentially “fills” the tube with gas molecules. The probability of finding a molecule at a specific radius is proportional to that radius (i.e., more molecules will be located at larger radii). Random starting positions are selected from a random radius and polar angle given by 1/2 r0 = arf 1 , θ0 = 2πrf 2 , 42 where x0 = r0 cos θ0 , y0 = r0 sin θ0 , z0 = Lrf 3 , and rf 1 , rf 2 , and rf 3 represent subsequent random fractions on [0, 1]. In order to simulate an infinite tube, I set a finite range 0 ≤ z ≤ L and apply periodic boundary conditions at z = 0 and z = L. Therefore, the initial z position z0 is just a uniform random position on [0, L]. Let r(t) denote the vector position of the tracer at time t, so that r(0) = {x0 , y0 , z0 } is the starting position at t = 0. The molecule is then assigned a random direction, selected from those available to it (e.g., the starting position has all directions available, so the direction is chosen from the unit sphere), cos θ = −1 + 2rf 1 , sin θ = (1 − cos2 θ)1/2 , φ = 2πrf 2 , (3.1) where θ is the polar angle between 0 and π, φ is the azimuthal angle, and the subsequent random fractions rf 1 and rf 2 are different from those used to define the initial position. The direction vector follows as dr = {dx, dy, dz} = {sin θ cos φ, sin θ sin φ, cos θ}. (3.2) The final step is choosing a molecular speed c0 according to the method described in §2.2.1.1. The molecule now has a starting position and velocity vector from which to start its motion. The next step is to calculate the first collision location. A free-path length ` is generated according to the method described in §2.2.3.1. In addition, the collision location on the cylinder surface is calculated using vector geometry and ray-tracing algorithms. The cylinder is defined by the equation x2 + y 2 = a2 on 0 ≤ z ≤ L, and the caps by x2 + y 2 ≤ a2 at z = 0 and z = L. The ray or vector of the molecule 43 along its direction of travel is given by R(s) = r(0) + sdr, where the desired value for s is the path length to collision with the tube surface. Therefore, I am interested only in the smallest positive value of s. The components of the ray equation are substituted into the equation for the cylinder to give (x0 + sdx)2 + (y0 + sdy)2 − a2 = 0, which produces a quadratic equation in the variable s, As + Bs + C = 0 → s = 2 −B ± √ B 2 − 4AC , 2A with A = dx2 + dy 2 , B = 2(x0 dx + y0 dy), C = x20 + y02 − a2 . Two solutions for s are obtained: one positive (along the direction dr) and one negative (in the direction opposite to dr)—technically, no solution is a possibility since the molecule could be moving exactly parallel to the z axis (though highly unlikely, if this were the case, the molecule would never intersect the tube and it is of little interest). Only positive solutions for s are of interest since these represent the molecule moving along the desired path to collision. If s < `, the molecule collides with the tube before moving a free-path length. In this case, the new vector position is given by r(t) = {x, y, z} = r(0) + sdr, where t = s/c0 is the time required to move a distance s at speed c0 . The only exception arises if the z 44 coordinate of the new position exceeds the computational domain defined on 0 ≤ z ≤ L. In this case, appropriate periodic boundary conditions are applied to translate z to the range of interest. The z coordinate in the computational domain is referred to as the local position, and the actual z position before applying the boundary conditions is called the global position. The local position is used in the simulation so that z remains on the order of x and y, and the global position is used to calculate overall molecule displacement from its initial position. Once the new position r(t) is set, a new set of directions and a new speed are chosen according to §2.2.2.1 and §2.2.2.2, respectively. For the new directions, 1/2 cos θ = rf is aligned with the normal direction (i.e., the radial direction perpendicular to the cylinder axis). In Cartesian coordinates, the direction vector for a molecule leaving the cylinder surface is dr = {x/a, y/a, 0} cos θ + {−y/a, x/a, 0} sin θ cos φ + {0, 0, 1} sin θ sin φ = nx a cos θ − o y y x sin θ cos φ, cos θ + sin θ cos φ, sin θ sin φ . a a a On the other hand, if ` < s, the molecule will move a free-path length ` before colliding with the cylinder. This situation arises when shifting attention from the free-molecule limit so that gas interactions start to play a role on mass transport. Without directly simulating gas-gas collisions, this method allows their inclusion while being less computationally expensive relative to true directsimulation MC or MD algorithms. The new position is given by the same equation for r(t) above, but t = `/c0 . A new direction and speed follow from equation (3.1), equation (3.2), and §2.2.1.1. The simulation procedure above repeats from the updated position with the new direction and speed, calculating the next collision location and so on. This continues for a preset amount of time or until the molecule has experienced a preset number of collisions. Throughout the simulation, various quantities can be monitored and saved, including the tracer positions, the path lengths between each collision (called chord lengths), the distribution of speeds given to the tracer, the tracer displacement from its initial position over the specified time interval, and so on. These data are gathered for a large number of tracers so that a wealth of statistical information is available from a single simulation. The single tube model is of particular importance because it provides a base case, or perfect-pore model of a single, straight, cylindrical pore. Further, results for free-molecule flow in the tube are easily 45 Transmission probability, Q 1.0ç ç 0.8 0.6 0.4 0.2 0.0 0.01 ç ç ç ç çç çç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç çç çç çç ççç ç 0.1 1 ççç çççç ç 10 100 ç 1000 La Figure 3.2. Transmission probability of molecules through a diffuse-reflecting cylindrical tube of length L and radius a for Knudsen flow. The solid dark line is the analytic result from Clausing [11] reproduced in equation (3.3) and the points are from Monte-Carlo (MC) tracer simulations. The MC results are the average of three cases at each L/a using 106 molecules per case. In each case, the molecules have different random starting positions at the tube entrance. validated to analytic results from Clausing [11]. When a negative pressure difference ∆p is maintained across a cylindrical tube of length L and radius a, the net molecular flux J across the tube is given by J =−   1 ∆p Qc0 L , 4 kB T L where T is the gas temperature, and Q is the probability that a gas molecule entering the tube will exit a distance L away. This equation is strikingly similar to equation (2.30) and is written such that the term in parentheses is the Knudsen diffusivity DiK . The remaining terms express the concentration gradient over the length of the tube. The total number of molecules incident on the tube entrance is given by equation (2.19) with n = ∆p/(kB T ) from the ideal gas law. The quantity nQ is the total number of molecules incident on the exit plane at z = L (assuming the axis of the tube is parallel to the z-axis and the entrance is located at z = 0). Clausing’s result expresses Q in terms of the tube dimensions, Q= i i 1 − 2α h 3 1−α h 2 2 2 2 2 1/2 3 2 2 1/2 2 4a + (L − 2a )(L + 4a ) − L + α + L − L(L + 4a ) + 2a , (3.3) 3a2 L 2a2 with α=f ! √ a 2a 7 √ , , L 3L + 2a 7 46 Transmission probability 1.0 B A 0.8 0.6 Aa=1 ç 0.4 2ç 3ç 4ç 5ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç ç 1 2 3 4 5 0.2 0.0 0 6 7 Ba Figure 3.3. Transmission probability of molecules through a 90◦ cylindrical elbow calculated using a Monte-Carlo (MC) algorithm for Knudsen flow. The tubes composing the cylindrical elbow have constant radius a, a length A before the elbow (upstream), and length B after (downstream). The solid line is a fit to results for MC calculations by Davis [15] and the points are MC data from the present work. As B/a → 0, each of the lines approaches the value from figure 3.2 for A/a = L/a. and   a z a L−z , =f , f = L L L L n o n  o  1/2 1/2 (L − z) (L − z)2 + 4a2 − (L − z)2 − z z 2 + 4a2 − z2 2 2 +4a2 ] Lz 2 −[2z−L][z 2 +4a2 ] − L(L−z) −[2(L−z)−L][(L−z) [z 2 +4a2 ]1/2 [(L−z)2 +4a2 ]1/2 , where z is the distance along the tube. For the validation case, the tracer simulation proceeds as described previously, with the exception that the initial z position of the molecules is z0 = 0 (i.e., all molecules begin their travel on one end of the tube). The molecules move between collisions as before, and their maximum penetration depth into the tube is stored. Each molecule is allowed to move until it returns to the plane z = 0 where it started, or reaches some maximum penetration limit. Molecules penetrating a depth L into the tube are said to have been transmitted that distance. Figure 3.2 shows the results of the MC simulation compared to Clausing’s analytic formula. The results are in excellent agreement and it is concluded that the simulation method is valid for Knudsen flow. The code was also validated against other MC results for Knudsen flow through a cylindrical elbow and in a cylindrical annulus [15]. These data are shown in figure 3.3 and figure 3.4, respectively, and are also in excellent agreement. 47 Transmission probability 1.0 L 0.8 ç çç ç 0.6 ç 0.4 çç ç ç ç 0.2 ç 0.0 0 ç ç ç ç ç ç 0.5 ç ç 0.75 ç 0.9 5 10 La ç ri a=0.25 15 20 Figure 3.4. Transmission probability of molecules through a cylindrical annulus calculated using a Monte-Carlo (MC) algorithm for Knudsen flow. The annulus has length L, outer radius a, and inner radius ri . The solid line is a fit to results for MC calculations by Davis [15] and the points are MC data from the present work. The case ri /a = 0 (not shown) reproduces the cylindrical tube results of figure 3.2. 3.3.2 Tracer Diffusion Between Aggregates of Sphere Particles The simulation procedure for molecules moving around and colliding with spherical particles is not inherently different from the method described in the previous section. However, there are a few outstanding differences that will be mentioned here. While the tube model provides a simple, singlepore representation, it is not a good model for most real porous materials. There has been some work extending the single tube idea to create parallel-pore models and other capillary networks [8, 22]. The general idea is that on some size scale, a pore can essentially be described by a cylinder which eventually intersects other pores forming nodes. It seems the most popular models for porous media are packed-sphere models. Here, the sphere particles represent the solid network, and the space between them the pores. By specifying a particle size distribution and certain packing rules, various forms of porous media can be created. These models seem to more closely represent porous materials found in engineering applications because the pore shapes are more irregular and, often times, highly unstructured. There has been much work focused on gas diffusion through porous media modeled as packed spheres. Some of the work has focused on random packed-sphere models [17, 55, 24, 44, 51, 67, 68], while others have concentrated their efforts on ordered particle arrangements [46, 60]. In addition, the same ideas have been extended to packed arrangements of other geometrical bodies [33, 60] and fractal-based geometries [39, 12, 40]. 48 (a) (b) (c) Figure 3.5. Classic examples of ordered, isotropic, monodisperse packed spheres and their unit-cell representations. The three figures show (a) simple cubic (SC) packing at maximum packing ratio (MPR), i.e., 47.6% porosity; (b) body-centered cubic (BCC) at MPR, i.e., 32% porosity; (c) facecentered cubic (FCC) at MPR, i.e., 26% porosity. 3.3.2.1 Generation of Porous Solids Using Packed Spheres Consider the three classic packings of uniform spheres shown with their unit cells in figure 3.5. For modeling purposes, an infinite isotropic porous sample can be created by periodically repeating the unit cell in all directions, so that sphere particles comprise the solid phase and the space between them the void phase. A range of porosities can be achieved for each configuration by allowing the particles to shrink or grow, while maintaining their central coordinate. The maximum packing ratio (MPR) is the state at which the distance between nearest neighbors equals one particle diameter, i.e., nearest neighbors are just touching. Figure 3.6 illustrates two examples of porous solids created from disordered or random sphere particles. The left image started as a BCC array, but particles were allowed to move a random distance up to one radius length in a random direction. There is no requirement for a particle to maintain contact with another particle. This would be necessary for the structural stability of a real structure, but the simulation case includes “floating”/unconnected particles. The image on the right is an aggregation of various-sized spheres. It represents a particle-based reconstruction of the 3D-SOFC anode previously pictured in figure 1.5. The method for creating this reconstruction is discussed in §3.3.3.1. 49 Figure 3.6. Randomized array of uniform-size sphere particles at 39.6% porosity (left), and a packing of various-sized spheres at 17.5% porosity (right). 3.3.2.2 Monte-Carlo Simulation of Tracer Molecules in Particle-Based Geometries As in the tube case, a tracer simulation begins by placing molecules at random positions within the void network. This is done by choosing a random point in the computation domain and checking to see if it lies within one of the particles representing the solid phase. If it does not, this is a valid point. If it does, a new point is chosen and the process repeats. The ratio of valid points to total number of iterations (i.e., the total number of points tried) gives a good estimate of the porosity  of the modeled porous medium. The initial position of the valid point is given by r(0) = {x0 , y0 , z0 }, just as before. The tracer is assigned a random direction dr according to equation (3.1) and equation (3.2), and a molecular speed c0 according to the method described in §2.2.1.1. The tracer then moves a distance equal to the smaller of: a free-path length ` generated according to the method described in §2.2.3.1, or the distance to the nearest surface in the direction of dr. As before, the collision locations are calculated using a ray-tracing method. The ray defining the molecule path is R(s) = r(0) + sdr. Each sphere particle composing the solid-phase network is defined by the equation (x−xci )2 +(y −yci )2 +(z −zci )2 = a2i or in vector notation, |R − Rci |2 − a2i = 0, where Rci and ai represent the center and radius of sphere particle i, respectively. Following the same method described for the tube case, the ray equation is substituted into the sphere equation to solve 50 for the path length to intersection s. This produces |r(0) + sdr − Rci |2 − a2i = 0, to give A = dr · dr = 1, B = 2dr · [r(0) − Rci ], C = [r(0) − Rci ] · [r(0) − Rci ] − a2i , which leads to s= −B ± √ B 2 − 4AC . 2A Only the positive real solutions to this equation are considered, and more specifically, the smallest positive nonzero solution. Possible solutions occur only when the discriminant is positive (a negative discriminant means there is no intersection between the ray and sphere). If the discriminant is positive, the smallest positive nonzero solution is chosen for s. These equations are solved for all possible intersections between the ray and each sphere particle i that falls within its path. The smallest value of s is stored, as well as the index i of the sphere with which the molecule intersected. This represents diffuse reflection from the surface of the sphere, aligned with the local unit normal. Figure 3.7. Pictorial illustration of gas behavior under different flow conditions. 51 If, instead, ` < s, the tracer moves a distance `, the position is updated, and a new set of directions and speed are assigned to the tracer from equation (3.1), equation (3.2), and §2.2.1.1. If s < `, the molecule collides with sphere i and the position is updated. Once the new position r(t) is set, a new set of directions and a new speed are chosen according to §2.2.2.1 and §2.2.2.2, 1/2 respectively. For the new directions, cos θ = rf is aligned with the normal direction (i.e., the radial direction pointing outward from the sphere particle). In Cartesian coordinates, the direction vector for a molecule leaving the sphere surface is dr = {cos Φ sin Θ, sin Φ sin Θ, cos Θ} cos θ +{cos Φ cos Θ, sin Φ cos Θ, − sin Θ} sin θ cos φ +{− sin Φ, cos Φ, 0} sin θ sin φ = {cos Φ sin Θ cos θ + cos Φ cos Θ sin θ cos φ − sin Φ sin θ sin φ, sin Φ sin Θ cos θ + sin Φ cos Θ sin θ cos φ + cos Φ sin θ sin φ, cos Θ cos θ − sin Θ sin θ cos φ}, where cos Θ = (z − zci )/ai , sin Θ = (1 − cos2 Θ)1/2 , cos Φ = (x − xci )/(ai sin Θ), sin Φ = (y − yci )/(ai sin Θ). This process repeats continuously until a certain diffusion time has been reached. The entire procedure repeats for a large number of tracers to provide complete sampling of the void network. Collision data, free-path (or chord) lengths between collisions, displacements, and so on are collected for each tracer. 52 3.3.3 Tracer Diffusion in Three-Dimensional Voxel-Based Geometries The MC methods used to simulate gas diffusion are not limited to solid-phase geometries constructed from sphere particles. In fact, any geometric object could be used. One of the goals of this work is to compare how estimated transport parameters in simple particle-based models relate to known, measured, or other predicted values in actual SOFC structures. If a set of rules can be developed so that particle-based models produce relatively accurate results compared to actual structures, these models will facilitate the ability to better study and understand SOFCs. Consider the actual SOFC anode reconstruction pictured in figure 1.5, created from a sequence of images obtained through focused ionbeam milling and scanning electron microscopy (FIB-SEM) [65]. In order to directly compare results between particle-based models and more realistic representations of the geometry, the first step is to create a geometric model for each case. There has been little work using MC to simulate gas diffusion in binary, pixelated, or voxel-based models [43, 7, 47]. 3.3.3.1 Creating Three-Dimensional Cermet Data from Two-Dimensional FIB-SEM Images The following description on how to construct a modeled cermet from two-dimensional (2D) images is transcribed from personal communication with another member of my research group, Vaughan Thomas. The first step in rendering is to put the data in the form of a locus of tagged points. Each point represents a location within the cermet and the type of matter at that point in the cermet, in this case, nickel, YSZ, or pore. This locus can be readily reconstructed from graphics files of the cross sections of the cermet. In the case of figure 1.5, the source images were provided from Wilson et al. [65] in the form of 2D “.tif” files. These images will be referred to as the Barnett data. The “.tif” data from 82 separate sections were reinterpreted into one data object, where each pixel from the “.tif” files was used to set the location and type of a volume element, a voxel, within a 3D structure. The voxels are then used to determine the type, size, and placement of spherical particles in a particle-based cermet model. The process is detailed below. The Barnett data is organized into 82 “.tif” files. Each “.tif” file represents one 2D image of the 53 Figure 3.8. FIB-SEM data as a “.tif” image (nickel: white, YSZ: gray, pore: black). Figure 3.9. Image comparison with original image on the left and reduced image on the right. cermet and can be read in MATLAB as an ordered 2D array of pixel color data. The Barnett data contains only three colors: white for nickel, 50/50 gray for YSZ, and black for void. The overall dimensions of the cermet in the x, y, and z directions are, respectively, 6 × 5.2 × 3.4 µm. Each “.tif” file represents a slice at constant z, showing the x-y plane. In the x and y directions each file is 430 × 370 pixels, respectively. Based on the aforementioned dimensions, the implied spacing for the pixels is δx = δy = 14 nm, δz = 42 nm, so that the spacing in the z direction is three times the spacing the x and y directions. Using this method, 82 new “.tif” images are generated at the lower resolution, 143 × 123 for x and y. At this point, the implied distance between pixels in the x, y, and z directions are equal, that is, δx = δy = δz = 42 nm. In addition, the pixels can be readily interpreted as cubic voxels in a 3D structure. The next step is to translate these voxels into a series of spheres. Note that the volume can equivalently be viewed as a lattice of equidistant cubes with an edge length equal to δx, or as a lattice 54 of equidistant overlapping spheres with radius √ 2δx/2. Choosing the radii of the replacement spheres has a strong impact on the final geometry and its fidelity to the original cubic voxel reconstruction. If the sphere radius is r = √ 2δx/2, then, at the boundaries between the solid phases and the void phase, a cap volume will project into the void space, distorting the geometry and reducing the volume of void space. At the same time, the apparent surface area of the solid phases is increased and the length of the triple-phase boundary (TPB) will be increased as well. At each interface of a solid-phase voxel in √  2π−3 contact with a void-phase voxel, the void volume is reduced by δVv = δx3 if the solid-phase 18 voxel is represented as a sphere with radius r = √ 2δx/2. Alternatively, this means the volume of the solid-phase voxel has increased by an amount 6δVv . In this case, the void fraction of the sphere-based reconstruction will be smaller than the void fraction of the voxel-based reconstruction. When a resolution of one sphere per voxel is used in the reconstruction, this does not make a material difference, but it will be shown that in other cases this distortion can become significant. An alternative is to make the volume of the spheres match the volume of the voxels. In this case, the radii of the spheres are r = (3/4π)1/3 δx. This approach still changes the geometry of the boundaries between the void- and solid-phase voxels, but preserves the void fraction and the relative solid fractions. However, there is a distortion to the geometry and the TPB. The formula of using one sphere per voxel results in a system with 1,424,709 spheres. Calculations with this many particles would be very computationally intensive and time-consuming. In order to make useful calculations, it is necessary to reduce the complexity of the system without losing important information. A simple solution is to group homogenous voxels and represent them with a single sphere. This is accomplished by evaluating, for each voxel, the matter type of the neighboring voxels. In this way, we determine the largest cube of side length n · δx centered on the voxel being evaluated, where n is a positive integer. Once this neighborhood is determined, a sphere of that matter type is placed, its center being the center of the voxel under examination, and the radius of the sphere being equivalent to (n + 1) · δx. An example of a sphere-based structure based on this technique for a single “.tif” file is shown in figure 3.10. The full 3D reconstructed cermet structure is pictured in figure 3.11. Table 3.1 gives a breakdown of the size and number of particles used in the reconstruction. In order to directly compare results from the particle-based reconstruction shown on the bottom 55 Figure 3.10. Image file and resulting particle array (nickel: green, YSZ: blue, pore: black).LETTERS 5.2 5 4 3.4 3.4 0 0 1 2 3 x (µm) 4 5 6 2.3 1.7 1.1 z (µ 1 m) 2 1 0 0 Figure 2 3D anode reconstruction. A view of the 3D reconstruction showing the Ni (green), YSZ (translucent/grey), and pore (blue) phases. 4.28 × 106 m cm−3 . Although TPB lengths have previously been estimated from 2D images of La0.8 Sr0.2 MnO3 –YSZ cathodes10 , this is the first direct measurement known to the authors. We are currently developing detailed electrochemical models of porous anode and cathode performance that connect TPB length and other microstructural information quantitatively to performance. In the meantime, a crude assessment can be done to determine if the measured TPB density reported above is consistent with the measured SOFC performance. We can estimate the expected anode area-specific resistance by multiplying the above TPB density by previously measured values of the lengthspecific resistance (LSR) of the TPB, and an approximate anode active depth (≈10 µm (refs 17,21)). On the basis of various studies including point-contact18,22 or patterned23,24 electrodes, reported LSR values vary in a range from 5 × 104 to 5 × 106 ! cm at 700 ◦ C, yielding anode area-specific resistance of 0.1–10 ! cm2 . Although broad, this range encompasses typical values for this type of SOFC anode25 and is consistent with our total SOFC area-specific resistance at 700 ◦ C (≈0.7 ! cm2 , from the data in the Supplementary Information, Fig. S1, which also included electrolyte and cathode resistances). The large uncertainties in this comparison reinforce the need to separate microstructure from other complicating variables, including utilization gradients, impurities, polarization, and polarization history, all of which have been shown to be important in overall anode performance18 . For a TPB to contribute to anode electrochemistry (equation (1)), the gaseous, ionic, and electronic phases adjoining the TPB must have a contiguous connection to the rest of the microstructure. That is, the pore must be connected through the surrounding pore network to the fuel stream, the Ni phase to the external electrical circuit, and the YSZ phase to the bulk YSZ electrolyte. As an initial step in evaluating phase connectivity, we have developed an algorithm to identify each of the contiguous TPB regions in the sample. (Note that contiguous TPBs imply contiguous adjoining phases.) In Fig. 3, different colours are assigned to these individual TPB segments. Analysis of the data indicates that 63% of TPBs are interconnected and 19% are short unconnected segments where the three phases may have poor connectivity with the rest of the microstructure; the remaining 0 4.0 2.0 z (µm) 3 2 y (µm) y (µm) 3 0 6.0 x (µm) Figure 3 3D map of the three-phase boundaries in the anode. Each colour represents a set of contiguous TPBs. The majority of the TPB length (63%) is connected (coloured white/grey). The remaining length consists of shorter, disconnected TPB segments (having colours other than white/grey). A fraction of these intersect the sample boundaries, and hence may be connected to larger segments existing outside the sample volume. However, a substantial fraction (19%) of the TPBs contact neither the highly inteconnected white/grey TPBs nor the sample boundaries, that is, they are actual short segments. 18% of the TPBs are undetermined because they contact the sample boundaries. A related question is tortuosity of the ionic, electronic, and gas-phase transport pathways within the electrode network. For example, a high gas-phase tortuosity at the anode inhibits the exchange of H2 O and H2 ; this can limit cell performance and is most readily observed as a limiting current behaviour at high current densities (for example, the rapid drop in cell voltage above ≈4 A cm−2 seen in Supplementary Information, Fig. S1). As a first attempt at evaluating gas-phase tortuosity, we assumed that transport within the pore subdomains is described by Laplace’s equation. Recent flux-based Monte Carlo simulations in porous media suggest that such an approach should generally be valid for both molecular and Knudsen diffusion, provided that the porosity is higher than ∼10% (ref. 26). Figure 4 schematically illustrates the calculation procedure. We first converted the 3D reconstructed volume of the gas pores shown in Fig. 2 into a finite-element mesh. In practice, the no-flux boundary condition at the pore–solid interface was replaced by setting the diffusivity in the solid to be much smaller than that in the pore and solving ∇ · (D(x)∇ψ) = 0 instead, where D(x) is the position-dependent diffusivity and −∇ψ is the flux. The test results using artificial microstructures showed that using a factor of 100 between the two diffusivities ensured that the results were not altered by this approximation. We used FEMLAB, commercial software from Comsol, to obtain ψ. The tortuosity was then recovered by equating the macrohomogeneously defined flux to the volume-average flux within the sample volume: !!! 1 ε ∇ψdV , − ∇ψMF = − τ V V Figure 3.11. Full reconstruction of the Barnett data pictured on the top. Bottom left: voxel-based reconstruction of the 3D-SOFC anode. Bottom right: particle-based reconstruction at 17.5% porosity (nickel: green, YSZ: blue, pore: translucent/gray). right in figure 3.11 to predictions from the actual structure pictured on the top, a voxel-based representation of the SOFC anode can be used. The voxelated porous solid is shown on the bottom left in figure 3.11. Developing voxel-based models of porous structures allows for more complicated geome- tries to be studied, and are seemingly only limited by the resolution of the images from which they where V is the sample volume, ε is the porosity, and τ is the tortuosity. This procedure was carried out for the three primary are constructed. The voxel size must be of the order of the size of the smallest structural features one 543 nature materials VOL 5 JULY 2006 www.nature.com/naturematerials hopes to model. ©2006 Nature Publishing Group !"#$#%&'()***) +,-.-/.0*++1/21)3*45 56 Table 3.1. Breakdown of the size and number of particles used in the particle reconstruction of a SOFC anode Radius (µm) 0.0846 0.1268 0.1691 0.2114 0.2537 0.2959 0.3382 0.3805 0.4228 0.4650 Average 0.1059 3.3.3.2 Number 17,945 4,969 1,871 761 300 127 32 10 6 1 Total 26,022 Surface Normals The main requirement of tracer MC methods used in this work is knowledge of surface normals at poresurface interfaces. Surface normals and tangent vectors are necessary to define the local orientation of a surface and the resulting influence this has on gas-surface interactions. For simple particle-based models of the porous solid, defining the surface normal at a point (i.e., a gas-surface collision location) is straightforward because all particle centers and radii are known. The real difficulty in defining surface orientations is inherent to voxel-based porous solids. The method used here is explained in the paragraphs below for a 2D (pixelated) system. Extension to 3D (voxels) is as simple as applying the method to the extra dimension. To construct the surface normal of the pixel at which a gas-surface collision occurred, the pixel hit by the molecule is not considered isolated, but together with its eight nearest neighbors (twenty-six for voxels). As an example, consider the pixel system shown in figure 3.12, where dark pixels indicate the solid phase and white pixels (including those labeled with ∗) indicate the void phase through which molecules move. The pixel hit by a molecule (in this case, the pixel labeled 0), is always taken as the center with its nearest neighbors surrounding it. This is the pixel for which the surface normal is desired. Assume the molecule originated from one of the void pixels to the right of or above pixel 0, i.e., those labeled with ∗. In this case, the local solid surface is only formed by the pixels labeled 0, 1, and 3. The normals for each of these pixels at their shared interface with a ∗ pixel is determined. For pixel 1, this gives {1, 0, 0} in the local coordinate system (meaning a unit vector pointing to the right). 57 1 2 * * 0 * 3 * Figure 3.12. Illustration of how the surface normal is determined for a two-dimensional pixel system. Dark pixels represent the solid phase and white pixels the void phase. The large arrow indicates the local surface normal at pixel 0 in the direction of the void pixels labeled with ∗. Note the surface normal in the direction of the blank void pixel would be a vector pointing left and down at 45◦ . Pixel 0 is bordered by a ∗ pixel to the top and right, so there are two outward normals giving {1, 1, 0}. For pixel 3 we are only interested in the normal in the direction of the ∗ pixels representing the region from which the molecule was located before collision. Hence, pixel 3 provides {1, 0, 0}. All of these vectors are illustrated as small arrows in figure 3.12. The surface normal at pixel 0 is then taken as √ the normalized sum of these smaller unit vectors, or {3, 1, 0}/ 10 (making an angle of 18.43◦ with the horizontal). This vector is represented as the larger arrow in figure 3.12. Likewise, the same analysis would lead to a surface normal in the direction of the blank pixel pointing left and down at 45◦ , or √ {−1, −1, 0}/ 2. 3.3.3.3 Monte-Carlo Simulation of Tracer Molecules in Voxel-Based Geometries Following a similar procedure to that described in §3.3.1 and §3.3.2.2, begin by choosing a random void voxel, and then a random starting position within that voxel given by r(0) = {x0 , y0 , z0 } = {xci + δx(rf 1 − 0.5), yci + δx(rf 2 − 0.5), zci + δx(rf 3 − 0.5)}, where Rci = {xci , yci , zci } is the center coordinate of void voxel i, δx is the voxel dimension (the voxels are cubic, so δx = δy = δz), and rf 1 , rf 2 and rf 3 are three subsequent random fractions on [0, 1]. A random direction dr is chosen, based on equation (3.1) and equation (3.2). Likewise, a molecular speed c0 is assigned according to the method described in §2.2.1.1. The simulation procedure then deviates somewhat from that described previously. Instead of using ray-tracing algorithms to calculate 58 the next collision location, the molecule takes a small step in its direction of travel. A good step size was found to be δs = 2.5 × 10−3 rf δx, where rf is a random fraction on [0, 1]. From this equation, the step length is uniformly distributed from 0 to 2.5 × 10−3 δx. Choosing δs in this manner provides good results without requiring a large amount of computational time. With this step length, the molecule position is updated according to r(t) = r(0) + δsdr, where t = δs/c0 is the time required to take the step. Then a check is performed to determine which voxel center Rci is closest to the current tracer position. If the molecule is still in a void voxel, it continues moving along its current path taking steps of random length δs. When a step results in a position inside a solid-phase voxel (of type nickel or YSZ), a simulated collision occurs. A new set of directions and a new speed are chosen according to §2.2.2.1 and §2.2.2.2, respectively. For the new 1/2 directions, cos θ = rf is aligned with the unit normal direction n = {nx , ny , nz } as determined by the method described in §3.3.3.2. The unit tangent vectors at the central voxel are calculated as t1 = and  q    {n , −n , 0}/ n2 + n2 y x x y q    {0, nz , −ny }/ n2y + n2z : nx 6= 0, : otherwise, t2 = n × t1 . The new direction vector is therefore dr = n cos θ + t1 sin θ cos φ + t2 sin θ sin φ, which represents diffuse reflection at the voxel described by surface normal n. The tracer is moved back to the center of the last void voxel it was in before collision, and the program continues with the new velocity vector. The program continues for a set diffusion time (constantly updated as t = t+δs/c0 ) and 59 for a specified number of tracer molecules. The various data are recorded throughout the simulation. 3.3.4 General Issues With Tracer Simulations in Systems With Complex Geometry Molecule intersections with the boundaries of the control volume (i.e., the boundaries of the volume representing the porous solid) generally result in elastic reflection back into the domain (in some cases, periodic boundary conditions are more suitable). Tracer displacements must be monitored within the local control volume, as well as in a global coordinate system determined by the boundary conditions on the control volume. The global position is required for displacement calculations. One major drawback of the tracer method is that initial placement of a molecule randomly in the pore space may result in an initial position inside an isolated pore. In reality, gas molecules would never have access to such regions. Since there is no way to a priori determine whether or not an initial position lies within an isolated pore, the statistics from these tracers (distribution of path lengths, displacements, etc.) must be removed from the data. This is accomplished by removing tracers that never undergo a boundary collision with the control volume (since isolated pores never intersect the boundaries). This requires the diffusion time be chosen large enough so that tracers starting deep within the porous solid can diffuse to the boundary. 3.4 Results and Discussion Diffusivity results for MC-tracer simulations in various models of porous media are presented here. Additional results and their implications will be presented and analyzed in the following chapters. 3.4.1 Diffusivity Results for Cylindrical Pores and Other Capillary Models The main purpose for calculating diffusion in cylindrical pores was to use it as a base case and validation of the simulation method. This is because analytic results are more easily obtained for Knudsen flow in a cylinder as compared to more complex geometries. That said, the fundamental results used for validation of the methodology, ray-tracing algorithms, random number generation, and so on were previously presented in figure 3.2, figure 3.3, and figure 3.4. Because of the level of agreement between 60 the results, the aforementioned highlights are assumed to be accurate. No other tests were performed to confirm additional validation (e.g., the random number generator was not tested separately). 3.4.2 Diffusivity Results for Packed-Sphere Models Recall the diffusibility ratio introduced at the end of chapter 2. The ratio is a dimensionless constant between zero and unity, which is characteristic of the porous material. Diffusibility results for the four packed sphere arrangements are shown in figure 3.13. The points are the results of MC simulations of Knudsen flow, and the line is a least-squares best-fit fourth-order polynomial used for clarity. For 1.0 1.0 0.8 SC pack Diffusibility, ΕDDiB Diffusibility, ΕDDiB 0.8 0.6 0.4 çø æ ç æ ç æ 0.2 BCC pack 0.6 0.4 çø æ ç æ 0.2 ç æ ç æ ç æ ç æ ç æ 0.0 çæ 0.0 ç æ ç æ 0.2 0.4 0.6 Porosity, Ε 0.8 0.0 çæ 0.0 1.0 1.0 0.2 0.4 0.6 Porosity, Ε 0.8 1.0 0.8 0.6 0.4 0.2 RAND pack 0.6 ç æ ç æ 0.4 ç æ ç æ 0.2 çø æ ç æ ç æ ç æ ç æ ç æ ç æ 0.0 çæ çæ 0.0 0.2 1.0 ç æ 0.8 FCC pack Diffusibility, ΕDDiB Diffusibility, ΕDDiB ç æ ç æ 0.4 0.6 Porosity, Ε 0.8 1.0 0.0 çæ 0.0 0.2 0.4 0.6 Porosity, Ε 0.8 1.0 Figure 3.13. Diffusibility versus porosity for various computational porous media composed of uniform sphere particles. The solid points are MC data from the present work for Knudsen flow. The line is a least-squares best-fit fourth-order polynomial used for clarity. Data are obtained for geometries in which the sphere particles overlap or are just touching (known as the maximum packing ratio)— this point is indicated with a star (?). Larger porosities are not investigated since this results in “floating”/unconnected particles. 61 Diffusibility, ΕDDiB or DD0 1.0 0.8 0.6 ó — SC — BCC — FCC — RAND ç Currie; spheres - experiment H1960L á Hoogschagen; spheres - experiment H1955L ó ñ Kim; spheres - experiment H1987L + õ Huizenga; spheres - experiment H1986L ó Trinh; spheres - 3D computation H2000L + Ryan; squares - 2D computation H1983L ó + ó + ó + 0.4 ó 0.2 óç ó ç ç ç õ õ õ ñññ ç õ õç á á õ + á + ç ç+ó 0.0 ó 0.0 ó 0.2 0.4 0.6 Porosity, Ε 0.8 1.0 Diffusibility, ΕDDiB or DD0 0.5 0.4 ó SC + á 0.3 0.2 ó ççç ççññ õõ çõ ñ áõõ á+õ õ óç BCC RAND ç ç+ó 0.1 FCC ó 0.0 ó 0.0 0.1 0.2 0.3 Porosity, Ε 0.4 0.5 Figure 3.14. Diffusibility versus porosity for various experimental and computational porous media composed of sphere particles. There is a noticeable difference between random and ordered structures using sphere particles. The bottom plot is a zoom of the highlighted window in the top plot. Diffusion data for SC, BCC, FCC, and RAND structures are the result of Knudsen flow simulations presented in figure 3.13. SOFCs, mass transfer occurs in submicron pores to give Knudsen numbers of order unity or greater [10]. Under this assumption (which will be upheld in the next chapter), only Knudsen diffusion simulations are considered. Data are obtained for the three regular classical structures (SC, BCC, and FCC) at various porosities up to the maximum packing ratio (MPR) for each structure. The porosity is varied by allowing particles to grow and overlap, while maintaining the location of their center. Beyond the 62 MPR (i.e., at higher porosities), particles are unconnected and floating in space. The porosity range for the RAND structure was allowed to extend over the entire range. The best-fit curves are plotted in figure 3.14 and shown with various experimental and computational results. For uniform particles, there is a very clear difference between random and regular arrangements. All three regular structures give higher diffusibilities at any given porosity than the random one. The BCC (red line) arrangement appears to yield the highest diffusibility. The curve for FCC (green) arrangements almost parallels the SC (blue) curve, giving slightly lower diffusibilities over the range of applicable porosities. The fourth-order polynomial fit for the FCC structure results in a negative value for the predicted diffusibility between 0% and 5% porosity. A negative number is unrealistic since the diffusibility can only vary between zero and unity; this is the direct result of the polynomial fit and not from tracer-diffusion simulations. The RAND (purple) curve shows good agreement with experimental results and the results obtained by Trinh et al. [60], who used a different computational approach than the one adopted here (a jump-diffusion scheme). The experimental data of Huizenga and Smith [24] was obtained specifically for Knudsen diffusion in random assemblages of uniform spheres. The RAND curve clearly passes through their data and the mass of other experimental data focused between porosities of 35% to 40%. 3.4.3 Diffusivity Results for Particle- and Voxel-Based SOFC Models Table 3.2 and table 3.3 present results from diffusion simulations based on mean-square displacements for particle- and voxel-based SOFC models, respectively. Perhaps the most important parameter listed in the tables is the tortuosity factor. This is essentially a measure of how much the pore network affects gas diffusion as compared to a diffusing gas under the same conditions in the absence of the porous medium. Previously reported values for the pore network of the SOFC electrode studied here are (τx , τy , τz )Pore = (2.1, 2.2, 1.9) [65]. It is believed that the anisotropy is due to the limited sample volume, and is a first indication of the accuracy that can be achieved for a given sample volume [65]. Similar values were recently reported by Iwai et al. [25] for a SOFC anode having volume fractions close to 26% Ni, 25% YSZ, and 49% pore space (a significant difference from the present case of 26% Ni, 54% YSZ, and 20% pore). Despite the fact that the pore space volume fraction differs by a factor of 63 Table 3.2. Summary of Knudsen transport data for the three-dimensional SOFC particle-reconstructed anode Ni 26.4 YSZ 56.1 Pore 17.5 Mean chord length,a h`i (µm) Mean-squared chord length, h`2 i (µm2 ) Ratio, h`2 i/(2h`i2 ) — — — — — — 0.192 0.092 1.251 β of equation (4.4)b — — 0.442 Diffusibility,c D/DiB Tortuosity, τ τx τy τz — — — — — — — — — — 0.027 6.42 10.14 8.08 4.08 Phase volume (%) The mean chord length gives an estimate of the mean pore radius, R, as appearing in equation (2.31) from the relation h`i = 2R. b The parameter β provides a measure of the nature of redirecting collisions at gas-solid interfaces. It is discussed in chapter 4. c The parameter D is defined by equation (2.28), and DiB in Der from equation (4.3). equation (2.32) using DiK = DiK a almost two and a half, Iwai et al. obtain virtually the same tortuosity factors for a significantly larger sample volume. These authors report values of (τx , τy , τz )Ni = (22.10, 29.46, 6.91), (τx , τy , τz )YSZ = (27.89, 14.95, 9.86), and (τx , τy , τz )Pore = (2.05, 1.99, 1.78). Iwai et al. [25] first reported tortuosity values for all three phases of an SOFC anode; the present work is the second. Though gases do not diffuse in the two solid phases, the tracer-diffusion scheme was used to obtain data for these phases in addition to the pore phase; this is done to provide some information on the complexity (or lack thereof) of the two solid phases and the transport pathways within the electrode network. The large discrepancy in tortuosity values (varying by a factor of ten or more between the pore and solid phases) reported by Iwai et al. is somewhat disturbing. While it is hard to make comparisons between the two different anodes (especially since the YSZ and pore volume fractions are opposite), Iwai et al. obtain tortuosity factors an order of magnitude higher for each of the solid phases. Each of these phases has a slightly higher volume fraction than the pore phase of the anode studied here and by Wilson et al. (approximately 25% versus 19%). It is presumed that Iwai et al. would obtain a pore-phase tortuosity factor for the present anode on the order of or greater than the 64 Table 3.3. Summary of Knudsen transport data for the three-dimensional SOFC voxel -reconstructed anode Ni 25.85 YSZ 54.29 Pore 19.86 Mean chord length,a h`i (µm) Mean-squared chord length, h`2 i (µm2 ) Ratio, h`2 i/(2h`i2 ) 0.525 0.535 0.971 0.557 0.673 1.084 0.261 0.137 1.000 β of equation (4.4)b 0.516 0.419 0.499 Diffusibility,c D/DiB Tortuosity, τ τx τy τz 0.046 5.56 10.50 5.90 3.64 0.335 1.62 1.99 1.57 1.41 0.039 5.05 8.39 5.93 3.26 Phase volume (%) The mean chord length gives an estimate of the mean pore radius, R, as appearing in equation (2.31) from the relation h`i = 2R. For Ni and YSZ, this could represent the mean radius of a solid particle. b The parameter β provides a measure of the nature of redirecting collisions at gas-solid interfaces. It is discussed in chapter 4. c The parameter D is defined by equation (2.28), and DiB in Der from equation (4.3). equation (2.32) using DiK = DiK a tortuosity factors obtained for the Ni and YSZ components—i.e., because of the disparity in pore-phase volume fractions between the two anodes, and since the anode used here has a pore volume fraction more in line with the two solid-phase volume fractions of the other anode, it is reasonable to assume Iwai et al. would calculate a pore-phase tortuosity factor of 20 or greater. However, it is still interesting that Wilson et al. and Iwai et al. obtain tortuosity values around 2 for their respective SOFC anodes, seemingly showing an independence from the volume fraction of pore space in the anode. There is discrepancy in the tortuosity value of 5 reported here, compared to an approximate value of 2 reported by Wilson et al. In addition, present results compared to Iwai et al. only show good agreement for the phases making up the majority of the anode—here the YSZ phase makes up 54% of the anode and tortuosity for that phase agrees with the reported value from Iwai et al. for the pore phase, which makes up 49% of their anode. The two Ni phases should be in better agreement since they each make up the same fraction of the anode. The tortuosity parameter is almost “mystical” in its definition, and various authors calculate it differently. Wilson et al. [65] assume that transport within the pores is described by Laplace’s equation, 65 Diffusibility, ΕDDiB or DD0 1.0 0.8 0.6 ó — SC — BCC — FCC — RAND ç Currie; spheres - experiment H1960L á Hoogschagen; spheres - experiment H1955L ó ñ Kim; spheres - experiment H1987L + õ Huizenga; spheres - experiment H1986L ó Trinh; spheres - 3D computation H2000L + Ryan; squares - 2D computation H1983L ó + æ Particle-based SOFC à Voxel-based SOFC ó ó + + 0.4 ó 0.2 óç ó ç ç ç õ õ õ ñññ ç õ õç á á õ + á + ç ç+ó 0.0 ó 0.0 ó æà 0.2 0.4 0.6 Porosity, Ε 0.8 1.0 Figure 3.15. The data from figure 3.14 with included results from the SOFC models. and recover the tortuosity by equating the macrohomogeneously defined flux to the volume-averaged flux within the sample volume. Recent flux-based Monte-Carlo simulations suggest that such an approach should generally be valid for both molecular and Knudsen diffusion [68]. Iwai et al. [25] use both a random-walk calculation and a lattice Boltzmann method to calculate the tortuosity; both approaches give similar results. The random-walk calculation assumes each walker jumps from a voxel to one of its neighbors. If the neighbor is of the same voxel type, the computation repeats, but if the selected voxel is of different type, the walker stays in the original voxel and waits for the next time step. Mean-squared displacements of the walkers are computed, and the tortuosity follows from the ratio of calculated diffusivity in the absence of the porous medium, to the calculated diffusivity in the presence of the porous medium (using equation (2.28) to calculate diffusivities from mean-squared displacements and time). Values reported here are based on the same principle, but the scheme to model molecular movement is different. Tortuosities are calculated from equation (2.33) as, τ = DiB /D, Der where D is defined by equation (2.28), and DiB in equation (2.32) using DiK = DiK from equa- tion (4.3). The curves in the plots showing D/DiB versus  have slopes equal to 1/τ . When the curves 66 become nearly linear, the tortuosity approaches a constant value with respect to porosity. The diffusibility values for all phases are plotted on figure 3.14 to produce figure 3.15. It is immediately evident that current SOFC electrodes are very far from being optimized for gas transport (at least compared to what can be achieved with packed spheres). Gas diffusion is severely limited within the porous anode, which directly impacts the delivery of fresh reactants to electrochemically active sites near the anode/electrolyte interface. 3.5 Conclusions Comparison between experiments and simulations clearly indicate that the coefficient of diffusion through a porous medium is a function of both internal geometry and porosity. Even when considering only transport through structures of sphere particles, there is a significant difference between regular and random arrangements, and still a clear difference between each of the regular geometries. Extending the computation scheme to the more complex geometries of reconstructed SOFC anodes shows that gas diffusion in traditional electrode designs is greatly hindered by the complexity of the porous network. The current method for producing such electrodes (high-temperature sintering of ceramic powders) is not optimal for gas diffusion; an optimal electrode will combine high porosity and connectivity for enhanced gas transport. In fact, one could achieve higher diffusivities if uniform packed spheres were used as the framework for the SOFC anode. For the geometries studied and a given porosity around 20%, the diffusibility ranges from approximately 0.4 to 0.18 (for the SOFC anode and BCC, respectively)—a factor of four and a half! In order of most optimal for gas transport, the structures are BCC, SC, FCC, RAND, and the SOFC anode. Random structures like RAND and the SOFC anode will contain isolated pores which do not contribute to available space for gas transport. This is perhaps one reason why these structures result in lower gas diffusibilities. 67 Chapter 4 Turning Porous Media Inside Out: What the Pore-Size Distribution Suggests About the Internal Structure 4.1 Introduction Typical porous solids are inherently complex three-dimensional (3D) structures. Apart from macroscopic properties, e.g., porosity and tortuosity, the distribution of pore sizes may arguably have the most signicant impact on the resistance of the structure to a flow moving within the void space. Pore shape, necking, and connectivity are certainly important as well, but these parameters are fundamentally more difficult to measure and characterize independently (at least the first two are, and connectivity can be scrutinized to some extent as we shall see later). Understanding the distribution of pore sizes in a variety of porous media can lead to signicant advancements in a wide range of phenomenon that rely on flow in these structures. The pore space in a typical porous solid consists of a highly convoluted network of interconnected paths. Details of the morphology control the rate and mechanism by which molecules of a gas move within the porous network. Transport mechanisms depend on pore-size distributions (PSDs), which cause local variations in diffusion rates [68]. Diffusion can occur via bulk, Knudsen, or surface processes, controlled by molecule-molecule interactions, molecule-wall collisions, or adsorption/desorbtion phenomena, respectively. For diffusion within the voids (as opposed to surface diffusion), the relationship between gas mean-free path and relevant pore dimensions determine the nature of the transport 68 mechanism. For example, flow within large pores (relative to the gas mean-free path) can exhibit transitional or continuum behavior, while flow within small pores must be analyzed under the Knudsen-flow assumption. Though these phenomena have been well understood for some time, the difficulty arises in the correlation between molecular mechanisms and influence of the physical structure of the pore network. In fact, the network structure itself is not well understood, so the relationship between structure and diffusion fluxes is not easily delineated. Apart from the porosity, most theories and models used to predict diffusion in porous media only consider a mean pore size and the tortuosity parameter [24, 46, 8, 50, 58]. The idea of adopting a mean pore size is born out of well-developed theories for diffusion in capillaries [34, 11, 48, 15]. In the case of a smooth pore-size distribution, we can readily define a mean pore size and the type of distribution (which turns out to be exponential in form for random packed spheres [38]), but some structures have several peaks in their distribution functions, and this is an interesting result because the distribution is not well described by its mean. Using a simple mean (i.e., the first moment of the pore-size distribution) is sometimes not enough to define or accurately represent the distribution. For multimodal distributions, it is necessary to account for transport at each mode. Further, specific engineering applications may benefit from smooth distributions, while others may rely on multimodal distributions targeted at optimizing some specific process, or even a host of processes at the same time. 4.2 Background There is no known way to a priori estimate the mean pore size of a porous solid, and, in any case, a simple one-parameter representation of the PSD can neglect potentially useful and significant information about the pore structure. Little work has been concerned with the shape of the PSD and its effects on mass transport in porous media. It is useful to first address how a PSD is obtained computationally for a particular structure. The PSD is assumed to relate to the distribution of chord lengths in the void spaces of a porous material. The distribution of chord lengths is obtained from tracer-diffusion simulations. Suppose we have a model pore structure, and gas molecules moving within the pore space. In the low-pressure limit 69 (rarefied or Knudsen-flow assumption), collisions between molecules can be neglected, so that the only influence on molecule behavior is the interaction between a molecule and the walls of the structure. Hence, molecules move from one wall collision to the next and their paths can be traced. A molecule, or tracer, is allowed to move for a given time, and the distribution of path lengths between collisions can be stored and analyzed for a large number of tracers dispersed throughout the porous network. This distribution of path lengths, or chord-length distribution, is used to represent the PSD since the chords provide a set of lengths which span the pores. Gas-surface interactions are prescribed by the cosine law [21] as a consequence of the second law of thermodynamics (see, for example Steckelmacher and Lucas [57] and Kolasinski [35], §3.8). The cosine law ensures that incident and reflection angles, and therefore, successive path lengths between collisions, are uncorrelated. Lu and Torquato [38] showed that ballistic tracers moving under the Knudsen-flow assumption in voids between random packed spheres follow free-path distributions equivalent to the PSD of the network. Further, this distribution was found to be exponential of the form, f (`) =   1 ` exp − , h`i h`i (4.1) where ` is the distance between surface collisions under the Knudsen-diffusion assumption (molecules travel in straight lines from wall collision to wall collision), and h`i is the first moment of the distribution, or mean-free path to collision. The function f (`) represents the probability distribution function for finding a pore of size `. Derjaguin [16] was perhaps the first to study the nature of length distributions between wall collisions in the Knudsen regime. Studying the behavior of the gas in the Knudsen limit removes any affect the molecules have on one another, so their paths are simply a measure of the distances between wall collisions. The PSD is purely a structural feature, and studying transport in the Knudsen limit means the only influence on the molecules is interactions with the structure. Derjaguin’s work also included the role of redirecting collisions and was later utilized by Levitz [36, 37] to characterize the void space within a dense random packing of hard spheres. Structural properties in packed beds have also been examined by Reyes and Iglesia [51]. Results of the work discussed above are very interesting but seem to garner little attention, recently appearing with good discussion in Zalc et al. [68]. The collective 70 analysis of these studies reveal that the ubiquitous use of a mean pore size as the relevant length scale is nonrigorous for some void structures. Generally, the length scale R appearing in equation (2.31) is approximated as R = h`i/2, and equation (2.31) becomes DiK = 1 h`ic0 . 3 (4.2) The proposed correction factor of Derjaguin accounts for the second moment of the pore-size distribution as well as the nature of the redirecting collisions to give 1 Der DiK = h`i 3   h`2 i − β c0 , 2h`i2 (4.3) where the nature of redirecting collisions is captured in β=− m=∞ X m=1 hcos γm i, (4.4) and hcos γm i is the average cosine of the angles between trajectory segments separated by m wall collisions [68]. Derjaguin showed that the Knudsen cosine law leads to hcos γm i = (−4/9)m for a stream of molecules striking randomly placed spheres. This results in β = 0.3077. The value of β is different depending on the geometry, reflection law, and diffusion process. Though highly useful, these works do not particularly focus on how the PSD changes for various structures, and consequently, what impact this has on Derjaguin’s originally proposed correction factor (or on mass transport in general). To the author’s knowledge, no work has rigorously examined the differences in PSDs in random and ordered materials. 4.3 Methods Procedures for generating model porous solids and simulating tracer diffusion are described in chapter 3. The Monte-Carlo (MC) computational scheme used herein and discussed previously is similar to that well documented in the literature [15, 17, 55, 67, 68]. Under the Knudsen-flow assumption, tracers follow ballistic trajectories that are independent of each other and may be generated as such. Data 71 is gathered for a large number of tracer molecules moving within the pore network. Among the data stored during this process are the many path lengths between surface collisions. This distribution of path lengths is used as a representation of the pore-size distribution. It is essentially the same as computing the distribution of distances from a surface to all other surfaces in all directions. However, the MC computation allows this to be computed in situ with tracer-diffusion calculations from the previous chapter. 4.4 Results and Discussion The chord-length distribution for each geometry is obtained from Knudsen-diffusion simulations of at least 10,000 molecules placed randomly in the voids of the structure. Molecules move from wall collision to wall collision until experiencing a minimum of 1,000 collisions, and collisions are assumed to be diffuse (post-collision directions are not dependent on incident angles). Each path length, or chord, between collisions is stored. Using a larger number of molecules or allowing them to undergo more collisions does not have any appreciable impact on the distribution. The figures presented in this section will show typical chord distributions for several geometries. These distributions are plotted as normalized probability distribution functions with area under each curve equal to unity. The thick lines in the plots are from tracer data, and thin lines give the exponential distribution (when applicable) of equation (4.1) using the mean from the tracer data. For each distribution, data is presented on two side-by-side figures using different horizontal scales, linear (left) and log (right). The log scale is used to reveal more details at smaller chord or pore sizes. These figures are meant to illustrate the vast differences in estimated pore-size distribution for different structures, even for those structures that appear quite similar. The shape of the distribution is useful in identifying pore types (e.g., cylindrical or tube like pores versus spherical pores versus any other shape), and can prove useful in design applications where one pore type will suit the design guidelines and transport requirements better than another. The distributions are also useful in trying to pinpoint why two structures with the same porosity differ in diffusibility, i.e., why one structure results in more facile gas transport compared to another. Cumulative probability distributions are also presented. 72 1.5 1.5 1.0 0.5 Probability 1.5 1.0 0.5 1 2 3 Normalized chord length, {R 1.5 1.5 Probability 2.0 0.5 1 2 3 Normalized chord length, {R 2.0 1.5 1.5 0.5 1 2 3 Normalized chord length, {R 2.0 1.5 1.5 0.5 0.0 0 1 2 3 Normalized chord length, {R 4 10 0.01 0.1 1 Normalized chord length, {R 10 0.01 0.1 1 Normalized chord length, {R 10 0.5 2.0 1.0 0.01 0.1 1 Normalized chord length, {R 1.0 0.0 0.001 4 10 0.5 2.0 1.0 0.01 0.1 1 Normalized chord length, {R 1.0 0.0 0.001 4 10 0.5 2.0 1.0 0.01 0.1 1 Normalized chord length, {R 1.0 0.0 0.001 4 Probability Probability Probability Probability 0.0 0.001 4 1.5 0.0 0 L/R → ∞ 0.5 2.0 0.0 0 L/R = 1000 1 2 3 Normalized chord length, {R 1.0 2.0 0.0 0 L/R = 0.5 (blue) = 5 (magenta) = 50 (yellow) = 500 (green) Probability 2.0 Probability L/R = 0.2 (blue) = 2 (magenta) = 20 (yellow) = 200 (green) Inside a cylinder (log) 2.0 0.0 0 Probability L/R = 0.1 (blue) = 1 (magenta) = 10 (yellow) = 100 (green) Probability Inside a cylinder (linear) 1.0 0.5 0.0 0.001 Figure 4.1. Probability distributions of chord lengths for molecules “bouncing” around inside cylinders of various aspect ratios L/R. The area under each curve is one. The molecule reflects diffusely from the inner surface of the cylinder after each collision. Curves are produced from data obtained in tracer simulations. 73 L/R = 1000 L/R → ∞ Cumulative probability Cumulative probability L/R = 0.5 (blue) = 5 (magenta) = 50 (yellow) = 500 (green) Cumulative probability L/R = 0.2 (blue) = 2 (magenta) = 20 (yellow) = 200 (green) Cumulative probability L/R = 0.1 (blue) = 1 (magenta) = 10 (yellow) = 100 (green) Cumulative probability Inside a cylinder 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 1.0 0.8 0.6 0.4 0.2 0.0 0 1.0 0.8 0.6 0.4 0.2 0.0 0 1.0 0.8 0.6 0.4 0.2 0.0 0 1.0 0.8 0.6 0.4 0.2 0.0 0 Figure 4.2. Cumulative probability distributions of chord lengths for molecules “bouncing” around inside cylinders of various aspect ratios L/R. 4.4.1 Short and Long Cylinders The chord distributions from the inside of short and long cylinders are shown in figure 4.1. The distributions for short cylinders lack the tail at larger chord lengths that appear for longer cylinders. 74 This is because low-angle trajectories (i.e., those that are close to parallel with the cylinder axis) result in molecules exiting the tube (for the short-tube case) rather than suffering a collision after a long path length (like the long-tube case). All distributions peak at or very near `/R = 2, meaning the most probable length is the cylinder diameter and the most probable direction for a molecule leaving the surface is radially inward. For cases when `/R ≥ 100, the distributions are more or less symmetric, so the most probable length is also the mean. Cumulative distributions are shown in figure 4.2. The effect of cylinder length-to-radius ratio, L/R, on mean chord length is shown more clearly in figure 4.3. The data confirm that when `/R ≥ 100, the mean chord length is equal to the diameter of the cylinder. Interestingly, most pores in random porous media (especially those considered for SOFC electrodes) do not resemble long cylinders with large aspect ratios. Instead, the aspect ratios are probably closer to unity as a result of the isotropic nature of the pore network. 1.0 X{\H2RL or X{2 \H2X{\2 L ç æ ç æ ç æ ç æ ç æ ç æ ç æ ç æ ç æ ç æ ç æ ç æ 0.9 ç æ ç æ 0.8 X{\H2RL ç æ 0.7 ç æ ç æ ç æ ç 0.6 æ ç æ 0.1 X{2 \H2X{\2 L ç æ ç æ ç æ ç æ ç æ ç æ 1 10 100 Length to radius ratio, LR 1000 Figure 4.3. Effect of cylinder aspect ratio, L/R, on mean chord length. 4.4.2 Inside a Single Closed Sphere The chord distribution from the inside of a single diffuse-reflecting sphere appears in figure 4.4. Here, we see a probability that increases with increasing chord length, up to the sphere diameter. The distribution is linear, peaking at the sphere diameter. This means the most probable direction for a molecule leaving the surface is radially inward, and is consistent with what is seen in the cylinder case. The mean of the distribution occurs at h`i/R = 4/3. The cumulative distribution is shown in figure 4.5. 75 Inside a sphere (log) 2.0 2.0 1.5 1.5 Probability Probability Inside a sphere (linear) 1.0 0.5 0.0 0 1 2 3 Normalized chord length, {R 4 1.0 0.5 0.0 0.001 0.01 0.1 1 Normalized chord length, {R 10 Figure 4.4. Probability distribution of chord lengths for molecules “bouncing” around inside a sphere. The area under the curve is one. The molecule reflects diffusely from the inner surface of the sphere after each collision. The curve is produced from data obtained in a tracer simulation. The distribution has a mean of h`i/R = 4/3, and most probable chord length at the sphere diameter, i.e., `/R = 2. Cumulative probability Inside a sphere 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 Normalized chord length, {R 10 Figure 4.5. Cumulative probability distribution of chord lengths for molecules “bouncing” around inside a sphere. 4.4.3 Monodisperse Simple Cubic Packing of Spheres Consider an infinite array of SC-packed spheres according to the unit cell shown in figure 3.5(a). The geometric chord-length distribution (or geometric PSD) obtained for such a structure (at various porosities) is shown in figure 4.6. Each of the thick colored lines represents the PSD obtained from tracer data. The thin lines are used only as reference, and they illustrate an exponential distribution with mean equal to the corresponding mean of the thick-line data. Clearly, the PSDs are multimodal, and assuming a single mean pore size as the relevant length scale is nonrigorous. To explain the multimodal nature, focus on the PSD for the case at R = RMPR . The first mode, or peak, falls in the range [0, 2) for `/R (if you imagine the tail of this mode to continue out past `/R ≈ 0.8). Referencing the unit cell in figure 3.5(a) and using geometric arguments, this peak is caused by tracers moving from the surface of a particle to another particle in a neighboring cell that shares a face with the unit cell (six total). The smallest possible path length in this case is identically 76 1.5 1.5 1.0 0.5 Probability 1 2 3 Normalized chord length, {R 1.0 0.5 0.0 0.1 4 2.0 1.5 1.5 Probability 2.0 1.0 0.5 0.0 0 R = 1.2RMPR  = 0.20 Probability 2.0 1 2 3 Normalized chord length, {R 0.5 2.0 2.0 1.5 1.5 1.0 0.5 0.0 0 1 2 3 Normalized chord length, {R 4 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R 1.0 0.0 0.1 4 Probability R = 1.1RMPR  = 0.33 SC pack (log) 2.0 0.0 0 Probability R = RMPR  = 0.48 Probability SC pack (linear) 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R 1.0 0.5 0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R Figure 4.6. Probability distributions for chord lengths in the void space between simple cubic (SC) arrangements of spheres. Each row represents a different particle radius (and therefore, a different porosity). The radius RMPR is the radius that gives maximum packing ratio, or the point at which nearest-neighbor particles are just touching. In terms of the unit cell pictured in figure 3.5(a), RMPR = L/2, where L is the length of any one side of the cubic unit cell. Thick lines are from tracer simulations and the thin line is an exponential distribution having the same mean. The area under each curve is one. zero (because the particles are touching), and the maximum possible path length that intersects both particles is exactly twice the particle radius. From the same geometrical arguments, the second mode is the result of tracers moving to a neighboring cell that shares an edge with the unit cell (twelve total). Geometry will predict the location of √ √ this mode to be 2 2 − 2 ≤ `/R ≤ 2 2 (or 0.828 . `/R . 2.828); these are obtained by analyzing the closest and furthest approach for a molecule traveling from one particle to another across a shared edge. This is exactly what we see on the plot. Likewise, the third peak results from movement across a shared √ √ vertex (eight total), located near 2 3 − 2 ≤ `/R ≤ 2 2 (or 1.464 . `/R . 2.828). There is clearly a lot 77 R = 1.2RMPR  = 0.20 Cumulative probability R = 1.1RMPR  = 0.33 Cumulative probability R = RMPR  = 0.48 Cumulative probability SC pack 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 1.0 0.8 0.6 0.4 0.2 0.0 0 1.0 0.8 0.6 0.4 0.2 0.0 0 Figure 4.7. Cumulative probability distributions for chord lengths in the void space between simple cubic (SC) arrangements of spheres. of overlap between the second and third modes, and they each decay to the same furthest approach, or maximum chord length. A fourth mode is also evident in the PSD in the approximate range (2.8, 3.8) for `/R. It is caused by tracers traveling more than one unit cell before collision. Trajectories of this nature have a higher relative probability at higher porosities (because there is more open space), and this trend is confirmed since this peak decreases as the particle size increases (decreasing the porosity). The shape of these distributions will result in a cumulative volume distribution that has a stepfunction behavior, meaning large chunks of volume are contained in only a few different pore sizes. This is shown in figure 4.7. 78 3.0 2.5 2.5 Probability 3.0 2.0 1.5 1.0 0.0 0 2.0 1.5 1.0 0.5 1 2 3 Normalized chord length, {R 0.0 0.1 4 3.0 3.0 2.5 2.5 Probability R = 1.1RMPR  = 0.14 BCC pack (log) 0.5 Probability R = RMPR  = 0.32 Probability BCC pack (linear) 2.0 1.5 1.0 2.0 1.5 1.0 0.5 0.5 0.0 0 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R 1 2 3 Normalized chord length, {R 4 0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R Figure 4.8. Probability distributions for chord lengths in the void space between body-centered cubic (BCC) arrangements of spheres. Each row represents a different particle radius (and therefore, a different porosity). The radius RMPR is the radius that gives maximum packing ratio, or the point at which nearest-neighbor particles are just touching. In terms of the unit cell pictured in figure 3.5(b), √ RMPR = 3L/4, where L is the length of any one side of the cubic unit cell. Thick lines are from tracer simulations and the thin line is an exponential distribution having the same mean. The area under each curve is one. R = 1.1RMPR  = 0.14 Cumulative probability R = RMPR  = 0.32 Cumulative probability BCC pack 1.0 0.8 0.6 0.4 0.2 0.0 0 1 2 3 Normalized chord length, {R 4 1 2 3 Normalized chord length, {R 4 1.0 0.8 0.6 0.4 0.2 0.0 0 Figure 4.9. Cumulative probability distributions for chord lengths in the void space between bodycentered cubic (BCC) arrangements of spheres. 79 4.4.4 Monodisperse Body-Centered Cubic Packing of Spheres Refer to figure 4.8 for PSDs of BCC arrangements of spheres. The BCC structure is more dense than a SC arrangement, and therefore the distributions shift slightly toward smaller pore sizes. The same kind of geometrical arguments can be applied to explain peak locations for this set of PSDs, and consequently, the cumulative volume distribution is characterized by a slight stepping behavior as shown in figure 4.9. With respect to the unit cell in figure 3.5(b), the first peak is caused by movement between the center particle and one of the vertex particles, or first nearest-neighbor interactions; the second is characterized by movement between vertex particles sharing an edge of the unit cell, or equivalently, between center particles of neighboring unit cells, known as second nearest-neighbor interactions. Notice how the first and second peaks are very much overlapped because the spacing in these situations does not differ significantly. In fact, there is little difference between the first two modes at low porosities. The third mode has a much lower probability than the first two, and is the result of movement from one vertex particle to another diagonally across a face of the unit cell, or third nearest-neighbor interactions. This is the equivalent of a molecule traveling from center particle to center particle in unit cells sharing an edge. Additional peaks are more prominent at higher porosities, as was true for the SC results These modes arise from molecules traveling greater distances in the increased void space before suffering a collision. 4.4.5 Monodisperse Face-Centered Cubic Packing of Spheres As with the previous two cases, PSDs for FCC structures have multiple modes and are illustrated in figure 4.10. The FCC packing is the most dense arrangement of the three already discussed. The PSDs have relatively smaller secondary (and so on) peaks, and it is worthwhile to note that as the porosity decreases in general over the range of ordered sphere packings, the lesser peaks become increasingly less prominent. That is, the PSDs tend closer to the exponential form shown as the thin line. It becomes increasingly difficult for a tracer molecule to move very far without reaching a surface. However, when the porosity gets too low, the pores become isolated and the peaks of the distribution will be sharper. 80 3 3 Probability 4 2 1 1 2 3 Normalized chord length, {R 2 1 0 0.1 4 4 4 3 3 Probability R = 1.1RMPR  = 0.085 FCC pack (log) 4 0 0 Probability R = RMPR  = 0.26 Probability FCC pack (linear) 2 1 0 0 1 2 3 Normalized chord length, {R 4 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R 2 1 0 0.1 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R Figure 4.10. Probability distributions for chord lengths in the void space between face-centered cubic (FCC) arrangements of spheres. Each row represents a different particle radius (and therefore, a different porosity). The radius RMPR is the radius that gives maximum packing ratio, or the point at which nearest-neighbor particles are just touching. In terms of the unit cell pictured in figure 3.5(c), √ RMPR = 2L/4, where L is the length of any one side of the cubic unit cell. Thick lines are from tracer simulations and the thin line is an exponential distribution having the same mean. The area under each curve is one. R = 1.1RMPR  = 0.085 Cumulative probability R = RMPR  = 0.26 Cumulative probability FCC pack 1.0 0.8 0.6 0.4 0.2 0.0 0 1 2 3 Normalized chord length, {R 4 1 2 3 Normalized chord length, {R 4 1.0 0.8 0.6 0.4 0.2 0.0 0 Figure 4.11. Cumulative probability distributions for chord lengths in the void space between facecentered cubic (FCC) arrangements of spheres. 81 1.5 1.5 1.0 0.5 1 2 3 Normalized chord length, {R 1.0 0.5 0.0 0.1 4 2.0 1.5 1.5 Probability 2.0 1.0 0.5 Probability 0.0 0 R = 1.2RMPR  = 0.15 Probability 2.0 1 2 3 Normalized chord length, {R 0.5 2.0 2.0 1.5 1.5 1.0 0.5 0.0 0 1 2 3 Normalized chord length, {R 4 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R 1.0 0.0 0.1 4 Probability R = RMPR  = 0.40 RAND pack (log) 2.0 0.0 0 Probability R = 0.8RMPR  = 0.66 Probability RAND pack (linear) 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R 1.0 0.5 0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0 Normalized chord length, {R Figure 4.12. Probability distributions for chord lengths in the void space between a random (RAND) arrangement of spheres. Each row represents a different particle radius (and therefore, a different porosity). Since this arrangement started as a BCC array as described in §3.3.2.1, the radius RMPR is the radius that gives maximum packing ratio for the BCC structure. In terms of the unit cell pictured √ in figure 3.5(b), RMPR = 3L/4, where L is the length of any one side of the cubic unit cell. Thick lines are from tracer simulations and the thin line is an exponential distribution having the same mean. The area under each curve is one. 4.4.6 Random Monodisperse Sphere Particles Pore-size distributions for a single random array of sphere particles is presented in figure 4.12. Various porosities were achieved by allowing the particles to grow or shrink uniformly, while holding their position constant. This is the same method used to vary the porosity in SC, BCC, and FCC arrangements. Though results are only presented for one random structure, the reader should be informed that all simulated random structures using uniform size particles produced similar PSDs, regardless of the amount of allowable particle overlap, and regardless of the algorithm used to create the particle packing. 82 R = 1.2RMPR  = 0.15 Cumulative probability R = RMPR  = 0.40 Cumulative probability R = 0.8RMPR  = 0.66 Cumulative probability RAND pack 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 2 4 6 8 Normalized chord length, {R 10 1.0 0.8 0.6 0.4 0.2 0.0 0 1.0 0.8 0.6 0.4 0.2 0.0 0 Figure 4.13. Cumulative probability distributions for chord lengths in the void space between a random (RAND) arrangement of spheres. Randomized structures do not contain repeated patterns of the microstructure, common to the ordered counterparts SC, BCC, and FCC. Hence, the PSD does not have multiple peaks that are evident in the latter cases and the cumulative volume distribution is very smooth (see figure 4.13). Instead, the PSD of random particle-based structures can be approximated quite well with the exponential distribution of equation (4.1). For porosities below  = 0.4 the PSD is essentially exponential in nature. 4.4.7 Real Structures: Particle- and Voxel-Based Reconstruction of a SolidOxide Fuel Cell Anode New imaging techniques have made possible the 3D reconstruction of actual SOFC electrodes (see, for example Wilson et al. [65] and Gostovic et al. [20] for focused ion beam scanning electron microscopy (FIB-SEM), and Izzo Jr. et al. [26] for high-resolution X-ray tomography (XCT)). More details and 83 Figure 4.14. Particle-based model of the reconstructed SOFC anode. Moving clockwise, the top left image shows the entire anode, the top right shows a cutout down the midplanes of the three major axes for particles in that region, the bottom right illustrates the oxide (YSZ) phase (translucent/gray), and the bottom left the metal (Ni) phase (green). Note that particles of different phases, shown separately in the two bottom images, do overlap. Particle 3D-SOFC (log) 2.0 2.0 1.5 1.5 Probability Probability Particle 3D-SOFC (linear) 1.0 0.5 0.0 0 1 2 3 Normalized chord length, {R 4 1.0 0.5 0.0 0.001 0.01 0.1 1 Normalized chord length, {R 10 Figure 4.15. Probability distribution for chord lengths in the pore phase of a particle-based threedimensional SOFC reconstructed anode. Particles are allowed to overlap; chord lengths are normalized by the average particle radius, R = 0.1059 µm; porosity  = 17.5%. A total of 26,022 particles are in a 6.04553 × 5.2 × 3.42439 µm box. The mean chord length is h`i = 0.192 µm. The thick curve is from a tracer simulation and the thin line is an exponential distribution having the same mean. The area under the curve is one. background on these techniques can be found in the next chapter. Using image data from Wilson et al. [65], the 3D-SOFC anode has been modeled using spherical particle-based reconstruction techniques, and a voxel-based mesh approach (discussion of the reconstruction technique can be found in §3.3.3.1, while discussions about the two model microstructures are 84 addressed in the next chapter). The particle-based model used here contains 26,022 spheres of various sizes, with pore-space volume fraction  = 17.5% (compared to 19.5% calculated by Wilson et al. from their reconstruction). Images of the particle-based structure are shown in figure 4.14. The PSD for the particle-based anode is presented in figure 4.15. The shape of the curve is quite peculiar; there are a lot of chords with length `/R . 0.4, while there is a dip in the number between 0.4 < `/R . 2.6, as compared to the exponential distribution (thin line). This feature is presumed to be an artifact of the particle model. Additional results have shown that changing the distribution of particle sizes can greatly reduce the deviation of the PSD from an exponential distribution, and we saw previously that the void space in a random packing of uniform-size spheres is characterized by the exponential distribution. Most of the area under the curve is found at chord lengths less than a micron, meaning it is likely that most pores are on the order of a micron or smaller. The PSD for the voxel-based anode is presented in figure 4.17. The shape of the curve is significantly different from the particle-based model except at the tail of the distribution. The voxel-based model results show a dramatic sawtooth behavior. This is solely an artifact of the resolution of the voxel model. The “teeth” occur at multiples of the voxel dimension, δx = δy = δz = 42.28 nm because molecules must travel through some integer number of voxels between every collision. Consequently, this shows up as a jump in the distribution at integer multiples of the voxel dimension followed by a decay in the probability to the next integer multiple. Several test cases using different voxel dimensions confirms this trend. Note that the smallest chord length that can be resolved is set by the resolution of the model, or the voxel dimension (this appears as the magenta vertical line). Practically all of the area under the curve is found at chord lengths less than a micron, agreeing with observations from the particle model, and likely indicating that most pores are on the order of a micron or smaller. The cumulative probability distribution for chord lengths in the pore phase of the particle- and voxel-based anodes are shown in figure 4.18 and figure 4.19, respectively. Note that the horizontal axis is the normalized chord length, `/R. For the particle model, the normalizing radius is the average particle size, 0.1059 µm. However, the normalizing radius in the voxel model is 1 µm since there existed no particle radius to use for normalization. As such, the cumulative distributions are actual quite similar if considered as a function of the chord length only (as opposed to the normalized chord 85 Figure 4.16. Voxel-based model of the reconstructed SOFC anode. Moving clockwise, the top left image shows the entire anode, the top right shows the void phase (blue), the bottom right illustrates the oxide (YSZ) phase (translucent/gray), and the bottom left the metal (Ni) phase (green). Note, due to graphic-plotting limitations, the structure shown here is for a voxel dimension of 84.56 nm. The actual structure analyzed has a voxel dimension of 42.28 nm, giving a single-voxel volume eight times smaller than those pictured. Voxel 3D-SOFC (log) 4 4 3 3 Probability Probability Voxel 3D-SOFC (linear) 2 1 0 0.0 0.2 0.4 0.6 0.8 Normalized chord length, {R 1.0 2 1 0 0.001 0.01 0.1 1 Normalized chord length, {R 10 Figure 4.17. Probability distribution for chord lengths in the pore phase of a voxel-based threedimensional SOFC reconstructed anode. The chord lengths are normalized using R = 1 µm; porosity  = 19.8%. A total of 1,424,709 cubic voxels are in a 143 × 123 × 81 array with a single voxel having dimension 42.28 nm on a side. The two vertical lines indicate the voxel dimension (magenta) and exactly twice this (orange). The mean chord length is h`i = 0.261 µm. The thick curve is from a tracer simulation and the thin line is an exponential distribution having the same mean. The area under the curve is one. length). Each is smooth, and reaches full probability around a chord length of ` = 1 µm. Again, this confirms that essentially all pore volume is contained in pores smaller than one micron. 86 Cumulative probability Particle 3D-SOFC 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 Normalized chord length, {R 10 Figure 4.18. Cumulative probability distribution for chord lengths in the pore phase of a particle-based three-dimensional SOFC reconstructed anode. Refer to the caption of figure 4.15 for more information regarding the model. Cumulative probability Voxel 3D-SOFC 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 Normalized chord length, {R 10 Figure 4.19. Cumulative probability distribution for chord lengths in the pore phase of a voxel-based three-dimensional SOFC reconstructed anode. Refer to the caption of figure 4.17 for more information regarding the model. 4.4.8 Moments of the Chord Distribution From each of the PSDs in figure 4.6, figure 4.8, figure 4.10, figure 4.12, and figure 4.15, the first and second moments of the distribution can be determined and plotted versus the porosity. The first moment is nothing more than the average chord length, which many have made ubiquitous use of as a descriptor of the PSD, without regard for the shape of the distribution. The first moment of the chord distributions for SC, BCC, FCC, RAND, and the particle-based 3D-SOFC are plotted in figure 4.20. Also shown on the plot is an exponential fit for the four arrangements composed of uniformly sized particles. Over the entire range of porosity, this fit can be used to estimate the mean chord length—or mean pore size—for packed-sphere models using a single sphere size. The predicted mean chord length is obtained from h`i = R exp(−1.50 + 3.83), where R is the particle size. When the particle size is greater than that pertaining to the maximum packing ratio, the particles overlap and there is some deviation from the prediction. Nonetheless, the general result can be substituted into equation (4.2) 87 Mean chord length, X{\R 10.0 5.0 ç æ X{\R = expH-1.50 + 3.83ΕL 2.0 æ Particle 3D-SOFC 1.0 ç æ 0.5 0.2 0.0 ç æ ç æ ç æ ç æ ç æ çø æ ç æ çø æ çø æ ç æ X{\R = 4ŐH3H1 - ΕLL 0.2 ç SC æ ç BCC æ ç FCC æ ç RAND æ 0.4 0.6 Porosity, Ε 0.8 1.0 Figure 4.20. Mean chord lengths of the pore space in sphere-based structures. All points are for uniform packed spheres except for the single black point from the particle-based 3D-SOFC model. Sphere particles are allowed to overlap or just touch (i.e., there are no floating particles, except perhaps in the RAND structure), and the star (?) indicates the maximum packing ratio for the structure (i.e., when nearest particles are just touching). The mean chord lengths are normalized using the particle radius (in the case of the 3D-SOFC model, the average particle radius is used). The solid line is an exponential fit through the data for uniform spheres only (SC, BCC, FCC, and RAND). The dashed line is generally regarded as a fit for nonoverlapping particles. so the Knudsen diffusivity can be written in terms of known, easily measured, or easily calculable quantities, such as the particle size, the porosity, and the mean thermal speed. This substitution gives DiK = 1 R exp(−1.50 + 3.83)c0 ≈ 0.074Rc0 e3.83 . 3 This result allows one to estimate the Knudsen diffusivity through a porous medium composed of uniform packed spheres without a priori information regarding the internal structure. The dashed curve in figure 4.20 is the expression h`i/R = 4/(3(1 − )), and is derived from a ratio of available volume to surface area for nonoverlapping, uniform packed spheres [55]. For particle-based models of porous media using spheres, this expression is generally used to estimate the mean pore diameter using the mean chord length [55, 24, 50]. The estimate is in good agreement with the three points denoted with a star (?), which represent the cases when particles are just touching for the three ordered arrangements. However, this relationship is only valid for nonoverlapping spheres of uniform size and clearly does not agree with the cases that involve overlapping particles. The result presented here is more general, and applies to a greater range of porosities and sphere geometries. The second moment of the distribution is the mean-squared chord length, and these are plotted in figure 4.21 with a power-law fit, h`2 i/R2 = exp(−2.39 + 7.73). The fits for the first and second Mean-squared chord length, X{2 \R2 88 100.0 50.0 10.0 5.0 ç æ Particle 3D-SOFC æ ç æ ç æ ç æ 1.0 0.5 ç ç æ æ ç æ 0.0 ç æ 0.2 çø æ ç æ çø æ çø æ X{2 \R2 = expH-2.39 + 7.73ΕL ç SC æ ç BCC æ ç FCC æ ç RAND æ 0.4 0.6 Porosity, Ε 0.8 1.0 Figure 4.21. Mean-squared chord lengths of the pore space in sphere-based structures. All points are for uniform packed spheres except for the single black point from the particle-based 3D-SOFC model. Sphere particles are allowed to overlap or just touch (i.e., there are no floating particles, except perhaps in the RAND structure), and the star (?) indicates the maximum packing ratio for the structure (i.e., when nearest particles are just touching). The mean-squared chord lengths are normalized using the square of the particle radius (in the case of the 3D-SOFC model, the average particle radius is used). The solid line is an exponential fit through the data for uniform spheres only (SC, BCC, FCC, and RAND). moments can be used in the corrected version of the Knudsen diffusivity appearing in equation (4.3). This substitution produces

Der DiK ≈ 0.074Rc0 e3.83 0.92e0.07 − β = DiK 0.92e0.07 − β . It is important to note that these predictions are obtained from architectures of uniform spheres. The results for the particle 3D-SOFC, which is composed of spheres of various sizes, deviates significantly from these estimates. This is perhaps due to normalization using the mean particle size. It might be more appropriate in such a case to normalize using the particle size in which most of the solid-phase volume is contained. The mean particle size for the particle-reconstructed anode is 0.1059 µm. The size in which the majority of the volume is contained is difficult to determine. Table 3.1 gives a breakdown of the size and number of particles making up the reconstructed anode. However, due to varying amounts of overlap between particles, and the fact that the same amount of overlap between a small and large particle eclipses a higher percentage of the volume of the small particle, estimating the particle size in which most of the volume is found is a challenging effort indeed. For an exponential distribution like that of equation (4.1), the mean is h`i and the variance is Var(`) = h`i2 . Since the variance can be written for any distribution as Var(`) = h`2 i−h`i2 , substitution X{2 \H2X{\2 L ratio of chord distribution 89 1.2 1.0 0.8 æ Particle 3D-SOFC ç æ ç æ ç æ çø æ ç æ æ ç çæ æ çø For exponential, this ratio is unity ç æ çø æ ç æ ç æ Based on function fits for uniform packed spheres 0.6 ç SC æ ç BCC æ ç FCC æ ç RAND æ 0.4 0.2 0.0 0.0 0.2 0.4 0.6 Porosity, Ε 0.8 1.0 Figure 4.22. A measure of the deviation of a chord-length distribution from an exponential distribution, based on the first and second moments. All points are for uniform packed spheres except for the single black point from the particle-based 3D-SOFC model. Sphere particles are allowed to overlap or just touch (i.e., there are no floating particles, except perhaps in the RAND structure), and the star (?) indicates the maximum packing ratio for the structure (i.e., when nearest particles are just touching). The solid line is a fit through the data for uniform spheres only (SC, BCC, FCC, and RAND). The dashed line at a value of h`2 i/(2h`i2 ) = 1 is indicative of the exponential distribution of equation (4.1). requires h`2 i/(2h`i2 ) = 1 for the exponential distribution. This ratio appears in the Derjaguin [16] correction factor found in equation (4.3). That said, we can plot this ratio for the particle models, producing the lines in figure 4.22. In a way, this plot gives a measure of the deviation of a particular distribution from exponential. The solid line is the ratio based on the functional fits and shows a slight upward trend. This ratio is close to unity for the cases studied. Even though many of the distributions had several peaks, those for SC, BCC, FCC, and RAND networks at least follow the general trend of the exponential distribution. In these cases, it is a good predictor of the average behavior of the distributions. 4.5 Conclusions Results clearly show that the pore-size distribution (PSD) of the space between overlapping and nonoverlapping sphere particles is highly dependent on the structure architecture. For any process that relies on transport in porous materials, analysis of the PSD can potentially be used to optimize the process or change the design to achieve better results. A multimodal PSD can be used to predict the various types of transport phenomena occurring within the material. In addition, while most models and predictors use a mean pore size as the only relevant length scale, we see that is nonrigorous for some void 90 structures. Structures characterized by a multipeaked PSD are more properly described by multiple length scales, one corresponding to each peak in the PSD. The pore network can then be characterized as a system of macropores (perhaps those in which the pore size is greater than the mean-free path of the gas moving within), mesopores, and micropores (a small pore size relative to the gas mean-free path). Such materials can be referred to as multiscale porous materials. The distribution of pore sizes in the SOFC models shows that the majority of pores are smaller than a micron. This perhaps is one of the causes of hindered gas transport, and upholds the assumption that Knudsen diffusion will dominate because typical pressures are not enough to force continuum dynamics. “Real” microstructures such as the SOFC anode model are likely characterized by spheres with a size distribution, rather than a single size, or by random packing. This means a multipeaked PSD is atypical, and the exponential distribution might be a good approximation. 91 Chapter 5 Piecing Together the Puzzle: Analysis of Solid-Oxide Fuel Cell Anode Microstructure 5.1 Introduction Current solid-oxide fuel cell (SOFC) electrodes are typically produced from high-temperature sintering of powders embedded with pore formers. The resulting structures are highly tortuous, complex threedimensional (3D) networks, optimal neither for gas transport nor electrochemistry. Not until the recent work of Wilson et al. [65] have we been able to visualize and analyze the microstructure of a SOFC cermet anode. The results of a 3D reconstruction of a cermet anode using focused ion-beam milling and scanning electron microscopy (FIB-SEM) are shown in figure 1.5. The microstructure is highly convoluted and only partially connected—Wilson et al. report that only 63% of the triple-phase boundary (TPB) regions are well connected to the metal, oxide, and pore networks. A large proportion (19%) of the TPB is formed by short segments disconnected from one or more networks, and therefore inactive for electrochemistry. Such a random, tortuous structure is clearly not optimized for its desired purpose, and clearly not what one would design with better control over the fabrication process. There is a lack of quantitative data describing the complex networks that make up the SOFC electrode microstructure. A SOFC demands high connectivity of all three phases, so that each phase can supply (or remove) reactants (or products) to (from) the TPB. That is to say, the metal (Ni) phase must be well connected from the electrolyte to the interconnect (see figure 5.1) for optimal electron transport. Likewise, the pore phase must be well connected for gas transport and the oxide (YSZ) phase 92 e– channel through open pores toward the dense-electrolyte interface, whereupon it may participate in charge-transfer reactions. In a system like Ni-YSZ, charge-transfer reac– tions occur in the vicinity of a three-phase boundary (TPB), which is formed at the intersection of the gas, electrolyte Anode (YSZ), and an electronic conductor (Ni). The products of the charge-transfer process (i.e., gas-phase species and elecO2– e trons) must be transported away from the three-phase rehod Cat gion. The gaseous product species are transported (diffusively and convectively) through the porous structure back – Air toward the flow channel where they mix with the unspent fuel and flow toward the exhaust. For a hydrocarbon fuel the net gas flow is from the electrolyte interface towards the fuel channel (e.g., for each mole of C4 H10 oxidized, five – – moles of H2 O and four moles of CO2 are produced). The electrons are carried in the anode metal toward the intercon– nect. The metallic phase in the anode must be contiguous and offer relatively little electrical resistance. An ode As illustrated in Fig. 3 the ceramic in the anode is an ion conductor, often the same material as the electrolyte. 2– 2– Electrolyte O 2– O Since ions can be conducted up into the porous anode, O the effective length of the three-phase region is extended Cathode greatly. This, in turn, improves the effectiveness of the – charge-transfer process. Because the electrolyte layer must be thin, one of the – electrodes usually serves as the structural support for the MEA. Either electrode can provide this function. The supFig. 3: Expanded view of the three-phase region in a membraneporting electrode is typically around 1 millimeter thick, electrode assembly. Primary particles of electrode and electolyte with the other electrode being only a few tens of microns Figure 5.1. Expanded view of the triple-phase assembly. material are typically on the boundary order of a micron. region in a membrane-electrode thick. The thick electrode should have a thermal-expansion coefficient matchedPicture to that of the electrolyte. Primary particles of electrode and electrolyte material are typically on the orderthat ofisaclosely micron. As the system is cycled between room temperature and opand description from Kee et al. [31]. The electrolyte must have sufficient conductivity to superating temperature, it is important that the thin ceramic port the ion current without substantial ohmic losses. Unelectrolyte not be over-stressed and crack. Any gas leaks fortunately, even “good” ion conductors offer considerable between the fuel and air channels lead to thermal oxidation resistance to ion transport. Resistance can be reduced by in(i.e., combustion), which produces no electricity and can 2− creasing temperature or making electrolyte layers and damage the system. must provide a contiguous path for O operating -ion transport. Additionally, a highoverheat density of electrochemically thinner. Since reducing operating temperature is desirable With hydrocarbon fuels, there may be a preference for for other reasons, electrolytes are usually made to be very anode support. In this case the electron-carrying metal (e.g., thin.regions Achievingfor electrolyte thickness of around 10 micromNi) could serve as a reforming or catalytic partial oxidation active TPB will provide abundant power-producing reactions. eters is practical, even for large-area MEAs. (CPOX) catalyst. On the other hand, a thick anode may also Ceria doped with either gadolinium or samarium has lead to long gas-phase residence times that could promote higher ion conductivity than YSZ, but also suffers some molecular-weight growth and the deleterious formation of electronic conductivity at high temperatures [6, 7]. Ceria carbon deposits. based systems usually operate around 600 ◦ C or below, because electronic conductivity increases to unacceptable lev2.2.2. Cathode els at higher temperatures. Load e Products Fuel e i e e e e e 5.2 Background The cathode’s job is to reduce oxygen using electrons and provide oxygen ions to the electrolyte. In the highly oxidizing high-temperature region where the cathode must Electrodes are designed to accomplish imporPerhaps one of the first studies to make use of digital imageseveral to represent a practical porous operate, metals like Ni are very support rapidly oxidized to nickel tant functions. These include structural support, gas-phase oxides. Therefore, cathodes are usually porous ceramic transport, electron transport, catalytic chemistry, and eleccomposites of an electrolyte (e.g., YSZ) and a mixed ionictrocatalysis. In SOFCs the electrodes are fabricated as was published by Meyerhoff et al. [43]. The authors present (2D) method forcasedeelectronic conductor (MIEC). In this there is no metal porous ceramics or ceramic-metallic composites two-dimensional (i.e., certo serve as the electrocatalyst and current collector. Rather mets), which are porous structures that provide interpenthe MIEC serves this function. For example a perovskite etrating, continuous, three-dimensional metal, oxide, and like strontium-doped lanthanum manganite (LSM) termining effective Knudsen-diffusion coefficients formetal a gas moving within material the pores, and provide some gas networks. The porosity, fraction, and ceramic is an MIEC material that shows good performance [1]. As fraction are each typically around 33%. with the anode the structure must be porous to enable transport of oxygen to the three-phase region. characteristics of the pore texture, 2.2.1. such Anode as the connectivity. More recent 2D analyses applicable to 2.2. Electrodes Within the porous anode fuel is transported from the flow 2.3. Interconnect SOFCs can be found in Wilson and Barnett [63] and discussed more in Wilson et al. [64]. The 2D 3 methods give reasonably accurate measures of volume fractions, surface areas, and TPB densities but are unable to quantify inherently 3D properties like tortuosity and phase connectivity. The development of high-resolution tomography has made possible the ability to reconstruct actual 3D-SOFC electrodes [65, 20, 26, 25], without relying on particle-based sphere-packing models (e.g., see Zalc et al. [68]). However, the question of whether or not a particle-based model is sufficient to represent an actual 3D-SOFC is one of the concerns addressed in this work. Comparison between particle-based models and voxel-based meshes of real anodes will be used to inform better designs and study the effects of 93 microstructure on gas transport. The method described here is based on that first presented by Wilson et al. [64]. The labels “dead-end” and “across” used by Wilson and colleagues are well-intentioned but give classification to sometimes very different or very similar networks for which a more appropriate label should be “unknown.” For example, “dead-end” is defined as a network containing exactly one boundary label. This means it enters the volume analyzed and dead-ends within. This in itself is not misleading. However, the network may connect externally to one of the other “dead-end” regions or even to one of the “across” regions. In fact, these networks could even represent isolated regions if they are unconnected to 3D analyzed, Analysis of SOFC Ni-YSZis Anode Interconnectivity anything outside of the representative volume. Albeit, within the volume which presumed 75 Figure 4. Three-dimensional visualization of the networks constituting the Ni phase. ~a! Shows the entire Ni phase; ~b! highlights the isolated Ni networks in red, and the dead-end networks in blue. The axis units are microns. Figure 5.2. Three-dimensional visualization of the networks constituting the Ni phase. (a) Shows the entire Ni phase; (b) highlights the “isolated” Ni networks in red, and the “dead-end” networks in blue. The axis units are microns. Picture andecreasing description from Wilson et treatment al. [64].of these effects may require finite element analyconnectivity while not that of YSZ significantly. In reality, however, the very low ionic conductivity of YSZ compared to gas diffusivity in pores and electronic conductivity in Ni must be considered. Thus, a TPB with a connection to the electrolyte via a very narrow YSZ pathway may be effectively equivalent to a disconnected TPB, which is not the case for Ni or pores. The higher YSZ content presumably yields similar effective degrees of connectivity of the three phases, such that transport through one individual phase does not become a dominant rate-limiting step. A full sis, as in previous pore tortuosity gas diffusion calculations ~Wilson et al., 2006! C ONCLUSIONS A technique was developed for using 3D image data to identify the separate contiguous networks constituting each 94 to be large enough to be representative of the anode, these are in fact dead-end networks. However, some of the networks labeled as “across” are also dead-end networks. The “across” label is used to define networks intersecting more than one boundary, and are considered contiguous throughout the electrode. There exists a gray area somewhere between what is a “dead-end” and what is a contiguous network spanning “across” the electrode (or at least the volume). Refer to the color visualization of the networks shown in figure 5.2, taken from Wilson et al. [64], and the two regions located approximately at x ∈ (2, 5), y ∈ (4, 5), and z ∈ (3, 3.5). In images (a) and (b) one region lies more-or-less due north and the other just to the right of the first. They both intersect the top (z ≈ 3.5) and back (y ≈ 6) boundaries (top and back respective to the viewpoint of the images), in which case they have been labeled “across.” Inspection of the image resolves more networks sharing a similar difficulty. Clearly these networks are no different than those labeled “dead-end” in that they enter the volume and deadend within. Nothing can be said for the structure outside of the volume, and they too could be isolated regions, true dead ends, or connected to a contiguous network. Based on the volume analyzed, it is possible to have a transport flux through the “across” networks from one boundary to another given appropriate boundary conditions. However, it is impossible for this to occur for “dead-end” networks (as conservation requires no net flux across the single boundary since these networks are not connected within the volume). So, perhaps that is the defining difference. However, it would be more appropriate to assign all such regions an “unknown” label. The catch-22 is that even the most well connected region across the volume could be assigned the “unknown” label. The only certainty is that “isolated” networks will indeed remain isolated for a larger volume. With this noted, the definitions of Wilson et al. to label networks as “isolated,” “dead-end,” or “across” are used in this research. 5.3 The Method Details of the voxel reconstruction of the 3D-SOFC anode from digital images is described in §3.3.3.1. The focus in this chapter is on the methods used to analyze the voxel representation of the SOFC anode. The voxel-based model used here contains 1,424,709 cubic voxels (143 × 123 × 81 in the x × y × z directions, respectively) of side length 42.28 nm, giving a sample size of 6.05 × 5.20 × 3.42 µm. A single voxel has volume 75,560.5 nm3 . Volume fractions of the three phases are 19.86% pore space, 25.85% Ni, 95 and 54.29% YSZ. This is consistent with the original reported values of 19.5%, 25.9%, and 54.6% [65]. Images of the voxelated structure are shown in figure 4.16. Determining the different phase networks starts with the complete set of voxels making up that particular phase. For the purposes of this explanation, this set will be called the original list. The search process begins with the first voxel and studies the six nearest neighbors to that voxel (i.e., the neighboring voxels that share a face with the original). Of course, only interior voxels have six neighbors, and the algorithm is amended as needed for boundary voxels. If any of the nearest neighbors are of the same phase type as the original voxel, they are appended to a new list representing one network of the phase. The same nearest-neighbor search continues with all the voxels added to the new list until the end of the list is reached. At this point, the list represents a set of voxels making up one particular network of the phase. A new voxel is chosen from the original list of all those making up the phase, ensuring that it has not already been added to any of the existing network lists. This voxel is the starting point of a new network, and the process continues as above until all voxels of the original list have been added to a network. The separate networks may contain only one voxel (in which case none of the nearest neighbors were of the same phase type), or there could feasibly be only one network containing all the voxels (which means all the voxels were connected to each other). This method also enables calculation of phase-boundary areas as well as the average connectivity of each voxel to another for any phase-type pair. Extending the scheme to consider more than the six nearest neighbors (for example, by allowing voxels that share an edge or vertex to connect) is not expected to change the network data. The method for calculating the TPB segments is similar to that used by Wilson et al. [64]. Starting with the voxel at coordinate location (or corresponding index) (0, 0, 0), consider the three edges in the direction of the remaining volume to analyze (because edges laying within a boundary plane cannot possibly be part of the TPB as a result of the nature of the geometric problem—this also prevents double counting TPB segments as the edges of each subsequent voxel are studied). At each of the three edges, the four voxels sharing that edge must contain at least one pore, one metal, and one oxide voxel as shown in figure 5.3. When all three are found to share an edge, that edge is added to the TPB. Upon finding a connection, or if all three phases do not share an edge, the search starts over at the next 96 Figure 5.3. Triple-phase boundary (TPB) curves of the voxel anode are made up of line segments. Each segment contributes 42.28 nm to the total TPB length, and is an edge shared by a metal (Ni), oxide (YSZ), and pore voxel. Edges that share endpoints are connected to form curves. edge. When the three edges have been searched, the algorithm moves to the next voxel and the process repeats. After all edge segments (each having length 42.28 nm) have been added to the TPB, these can be grouped into different curves based on whether or not segments are connected. A connection between segments is formed if they share a common endpoint, and TPB curves are constructed by joining these segments until no more connections exist. A different methodology is implemented to determine the phase networks and TPB density for particle-based anodes like the one shown in figure 4.14. Since the center coordinates and radius of all particles are known, a list is created that contains all pairs of intersecting particles, i.e., the distance between their centers is smaller than the sum of their radii. Pairs of particles of the same phase type are linked to other pairs to build up the different phase networks. Pairs of different phase types are used to determine shared surface areas and the TPB, because a TPB requires the intersection of metal and oxide in the presence of the gas (void) phase. Noting that an intersection between particles i and j is the same intersection between j and i, the list is further reduced to avoid this double-counting. At this point, the list only contains pairs of particle indices i and j for which |Rci − Rcj |2 < (ai + aj )2 and Type(i) 6= Type(j), where Rci and ai represent the center and radius of particle i, respectively. The intersection region 97 Pore space Figure 5.4. Triple-phase boundary (TPB) curves of the particle anode are made up of discretized circular arcs. Circles are formed at the intersections of metal (Ni) and oxide (YSZ) particles, and the radius depends on the amount of overlap of the particles. The circles are discretized using 100 points, so that each point represents 1/100 of the length of the original circle. Points contained within another particle are removed and do not contribute to the total TPB length because they are not exposed to the pore space. Arcs are connected to form curves if their closest approach is less than the sum of the representative distances of the points that represent that arc. between any two sphere particles is a circle that can be parametrized as follows: 1. Let n be the unit normal pointing from the center of particle j to the center of particle i, n= Rci − Rcj . |Rci − Rcj | 2. Let R?cij denote the circle center and a?ij the radius. The law of cosines can be used to first ? determine θij , the angle between the vector n and a radial vector pointing from the center of particle i to any point on the circle. Note that this vector has magnitude ai . The complete triangle is formed by this vector, its counterpart pointing from the center of particle j to the same point on the circle (having magnitude aj ), and the vector between the centers of i and j. Using ? this triangle, the angle θij is opposite the side with length aj , such that ? θij = cos−1 a2i + |Rci − Rcj |2 − a2j 2ai |Rci − Rcj | ! . 98 The circle radius and center follow as ? a?ij = ai sin θij , ? R?cij = Rci − nai cos θij . 3. Let u be a unit vector perpendicular to n = {nx , ny , nz } so that u · n = 0 and u follows simply as u=  q    {n , −n , 0}/ n2 + n2 y x x y q    {0, nz , −ny }/ n2y + n2z : nx 6= 0, : otherwise. 4. A third unit vector is found from v = n × u. 5. A point on the circle is given by p = R?cij + a?ij (u cos φ + v sin φ), where φ is the parameter. To discretize the circle into N equally spaced points, the angle φ varies from 0 to 2π − 2π/N in steps of ∆φ = 2π/N . After all the intersection regions have been parametrized by a set of N points per circle, it is necessary to determine which of these points actually represent a segment of the TPB (see figure 5.4 where N = 100). A check is performed to determine if each of the points lies within another particle (i.e., the distance between the point and that particle’s center is less than the particle radius). This test excludes the two particles for which the point lies on their intersection. Points that are contained within another particle cannot possibly contribute to the TPB because they are not exposed to the gas phase, and consequently, they are disregarded. Only points that are not contained contribute to the total TPB, where each point represents a length segment of 2πa?ij /N . Connections between segments are made if the points are separated by a distance that is less than the sum of their representative lengths, and TPB curves are formed by building these connections until no more are found. 99 5.4 Results and Discussion Table 5.1 gives a summary of the volume, surface area, and connectivity information for the networks in each phase in the 3D-SOFC voxel-reconstructed anode. The Ni phase was found to be 90.2% connected between two boundaries of the anode, 8.5% made up of regions that contacted a single boundary and dead-ended within the volume, and 1.3% isolated or stand-alone networks. The respective numbers reported by Wilson et al. [64] are 86%, 12.7%, and 1.3% so there is a discrepancy between the “across” and “dead-end” volume fractions. It is interesting to note that 66.2% of the Ni phase is composed of a single connected network (this agrees with the value of 68% reported in Wilson et al.). Likewise, the YSZ phase is essentially 100% connected and the pore phase is 96.3% connected with the largest network of the pore phase constituting 90.6% of its volume. The numbers for the YSZ and pore phases agree with those reported by Wilson et al. The high connectivity of the YSZ phase is expected since it composes 54.3% of the anode by volume. However, it is interesting that the pore phase has a smaller volume fraction than the Ni phase (19.9% compared to 25.9%) yet it is better connected. Also reported in table 5.1 are the average connectivities of an interior voxel to a voxel of each phase type. The calculated total phase-boundary areas per unit volume are found to be 2.30 µm−1 for Ni (including interfaces with YSZ and pores), 4.21 µm−1 for YSZ (including interfaces with Ni and pores), and 3.41 µm−1 for pores (including interfaces with Ni and YSZ). The breakdown of these volumespecific phase-boundary areas can be found in table 5.1. Previously reported values for the same anode are 1.6 µm−1 for Ni, 3.0 µm−1 for YSZ, and 2.4 µm−1 for pores [65]. Though the numbers are not in excellent agreement by value, the ratio between the three is very close—1:1.835:1.485 (Ni:YSZ:pore, present results), and 1:1.875:1.5 (Ni:YSZ:pore, Wilson et al.). Similarly, table 5.2 gives a summary of the volume, surface area, and connectivity information for the networks in each phase in the 3D-SOFC particle-reconstructed anode. The Ni phase was found to be 87.2% connected between two boundaries of the anode, 11.5% made up of regions that contacted a single boundary and dead-ended within the volume, and 1.3% isolated or stand-alone networks. These numbers are in excellent agreement with the results reported by Wilson et al. [64] and also in line with the values in table 5.1 for the voxel anode mesh. The remaining data listed in each of the tables is very similar, confirming that such a simple particle model can be used in lieu of a detailed 3D data set to 100 Table 5.1. Volume, surface area, and connectivity information for the networks of each phase in the three-dimensional SOFC voxel -reconstructed anode having total volume of 107.65 µm3 Ni 53 18 20 15 YSZ 11 4 3 4 Pore 143 71 50 22 Phase volume (µm3 ) (%) 27.83 25.85 58.45 54.29 21.38 19.86 Volume % of “isolated” networks Volume % of “dead-end” networks Volume % of “across” networks 1.31 8.49 90.19 0.020 0.013 99.97 1.35 2.33 96.32 Number of voxels in largest network Volume % of largest network 243,937 66.23 772,847 99.92 256,191 90.56 Estimated specific phase-boundary area (µm2 /µm3 ) Ni-toYSZ-toPore-toTotal estimated specific phase-boundary area — 1.551 0.746 2.297 1.551 — 2.665 4.216 0.746 2.665 — 3.411 Average connectivity of an interior voxel (max is six neighbors) Ni-toYSZ-toPore-to- 5.623 0.122 0.1596 0.254 5.667 0.5699 0.123 0.211 5.2705 Total number of separate networks Number of “isolated” networks Number of “dead-end” networks Number of “across” networks Table 5.2. Volume, surface area, and connectivity information for the networks of each phase in the three-dimensional SOFC particle-reconstructed anode having total volume of 107.65 µm3 Ni 48 13 21 14 YSZ 11 3 3 5 Pore — — — — Phase volume (µm3 ) (%) 28.42 26.4 60.39 56.1 18.84 17.5 Volume % of “isolated” networks Volume % of “dead-end” networks Volume % of “across” networks 1.32 11.51 87.17 0.024 0.022 99.95 ∼1.5 — — Number of particles in largest network Volume % of largest network 5,247 64.44 17,310 99.90 — — Estimated specific phase-boundary area (µm2 /µm3 ) Ni-toYSZ-toPore-to- — 1.880 — 1.880 — — — — — Total number of separate networks Number of “isolated” networks Number of “dead-end” networks Number of “across” networks 101 represent a SOFC anode, at least in terms of phase volumes and connectivities. Figure 5.5, figure 5.6, and figure 5.7 give a graphical representation of the network phase volumes and network types for the Ni, YSZ, and pore phases, respectively. The top plot in each shows the relationship between a network and those equal or smaller in size, whereas the bottom plot is the cumulative volume distribution of the networks. Steep gradients on the top plots are indicative of many networks of similar size. The bottom plots are useful to see how much of the volume is contained in each of the separate networks making up that phase, as well as what the network classification is. As expected, most of the smaller networks are isolated within the volume, and most of the larger networks are connected across the sample volume. Lines are used to connect the points for clarity. The network type and size distributions compare very well between the voxel and particle models. For all phases, it is optimal to have all networks span across the anode to facilitate transport in that phase. The number of different networks is less important and there probably exists a balance between number of networks and available surface area for reactions. For the pore phase, almost 20% of networks by number (not volume) are composed of a single voxel or two voxels, as indicated by the stack of points at the vertical lines on the top plot of figure 5.7. The bottom plot reveals that these do not contribute significantly to the volume of pore space. Networks composed of only one or two voxels are most likely due to inherent errors in converting the FIB-SEM images to a usable model. The scheme of building networks based only on voxels that share a face also contributes, since expanding this to allow connections between shared edges or vertices will reduce the number of single-voxel networks and cause them to be absorbed into larger networks. Additional computations show that the overall volume percentages of the different connected networks remains unchanged under these circumstances, so the data presented here is not dependent on the methodology used to construct the different phase networks. Figure 5.8, figure 5.10, and figure 5.12 provide a 3D visualization of the “isolated” and “deadend” networks in each of the phases for the voxel anode (and analogously, figure 5.9, figure 5.11, and figure 5.13 for the particle anode). “Isolated” networks are colored red and “dead-end” regions blue. Recall that a “dead-end” network is defined as one that intersects only one boundary of the sample volume, so nothing can be transported across the volume through these networks. However, their true classification remains unknown, since we do not know if these form isolated regions in the entire anode, 102 or are parts of larger connected networks spanning the anode. In fact, two “dead-end” networks of different phase types may still contribute to “active” TPB regions. Assume that a Ni and YSZ network enter the sample volume and dead end within, but still intersect each other. While it is true that ionic and electronic transport paths will ultimately dead end because of this, the TPB created at their intersection (assuming the pore phase is also present) can still be active due to transport of ions and electrons to/from the reactive site. The plots and 3D visualizations presented here were created to allow direct comparison to those presented in Wilson et al. [64], one of which is reproduced in figure 5.2. Wilson et al. only provide visualization data for the Ni phase, and results for all three phases for each of the model reconstructions are presented here. Phase network connectivity has a direct impact on electrochemical kinetics and subsequently, the TPB regions of electrochemical activity. A TPB represents a reactive site where all three phases come together so that all constituents for the electrochemical reaction are present (electrons in the Ni/metal phase, O2− -ions in the YSZ/oxide phase, and gas or fuel in the pore phase). In order for a TPB to be active, it must lie on connected networks so that reactants and products can be transported to and from the site. Therefore, the distinction must be made that active TPB sites are found only on connected or “across” networks, and inactive sites on “isolated” networks. The activity of TPB edges falling on “dead-end” networks is unknown. More TPB length means more active sites, and thus, a higher probability that a reaction will occur. Optimization of the TPB can be achieved through processes that affect the microstructure or by using mixed ionic-electronic conducting materials. Table 5.3 summarizes the calculated TPB densities and activities for the voxel-based anode. It was found that 13.3% of the total TPB length falls on “isolated” regions (and is therefore inactive), 16.8% on “dead-end” networks (unknown activity), and 69.9% on connected networks (electrochemically active). On the scale of an entire anode, it is reasonable to assume that TPB curves of unknown classification will mimic the connectivity of those contained within the sample volume. If it is assumed that 13.3% of the 16.8% of unknown TPB is actually inactive—and the remainder is active—the amended TPB length is 15.5% inactive and 84.5% active. Of the 5.8 µm−2 total TPB density, this results in 0.9 µm−2 inactive and 4.9 µm−2 active. Wilson et al. [64] report an amended value of 3.7 µm−2 active TPB, 103 Table 5.3. Summary of the triple-phase boundary (TPB) density, activity, and curve lengths in the three-dimensional SOFC voxel -reconstructed anode having total volume of 107.65 µm3 —each edge has a length 42.28 nm “isolated” Inactive 1,964 73 “dead-end” Unknown 2,485 82 “across” Active 10,333 147 Total 14,782 302 Length of TPB (µm) (µm/µm3 ) (%) Amended TPB (µm/µm3 ) (%) 83.03 0.77 13.29 0.90 15.52 105.06 0.98 16.81 — — 436.84 4.06 69.90 4.91 84.48 624.93 5.81 100 5.81 100 Number of edges in longest curve Length of longest curve (µm) (µm/µm3 ) As % of curve type (%) 122 5.16 0.048 6.21 536 22.66 0.21 21.57 1,247 52.72 0.49 12.07 — — — — Total number of separate edges Total number of separate curves Table 5.4. Summary of the triple-phase boundary (TPB) density, activity, and curve lengths in the three-dimensional SOFC particle-reconstructed anode having total volume of 107.65 µm3 —each point represents a different distance based on the radius of circle that parametrizes each particle intersection “isolated” Inactivea — 25,311 47 “dead-end” Unknowna — 102,404 123 “across” Activea — 448,616 471 Total 24,434 576,331 641 Length of TPB (µm) (µm/µm3 ) (%) Amended TPB (µm/µm3 ) (%) 73.36 0.68 4.32 0.80 5.085 302.32 2.81 17.79 — — 1,323.73 12.30 77.89 14.98 94.915 1,699.41 15.79 100 15.79 100 Length of longest curve (µm) (µm/µm3 ) As % of curve type (%) 15.66 0.15 21.35 91.41 0.85 30.24 715.85 6.65 54.08 — — — Total number of Ni-YSZ particle intersections Total number of TPB points (using 100/circle) Total number of separate connected curves a The TPB activity classification is based only on the Ni and YSZ networks—if the pore network were included, some of the curves currently labeled as active will change to unknown or inactive. 104 from a total density of 4.2 µm−2 . Differences in reported values despite similar calculation methods may be sourced to the original conversion process of the 2D image data to a usable 3D model. Figure 5.14 gives a graphical representation of the TPB curve lengths and activity types. The top plot shows the relationship between each curve and those equal or smaller in length, whereas the bottom plot is the cumulative length distribution of all TPB curves. Steep gradients in the top plots indicate a high concentration of curves at that length. For the voxel anode (left plots), most curves have lengths between 0.8 and 3 µm, and the bottom plot reveals that these constitute approximately 40% to the overall TPB length. Most of the inactive TPB curves are shorter than 3 µm, and active curves are distributed throughout the range, from 0.1 to 50 µm. The total TPB length is 625 µm, as reported in table 5.3. Figure 5.15 provides a visualization of the entire set of TPB curves in the voxel-based anode and highlights the inactive curves in red, unknown in blue, and active in white/light gray. This figure was created to allow direct visual comparison to Figure 5 in Wilson et al. [64], not reproduced here because the two images are very similar. It also facilitates comparison to the 3D visualization of TPB curves from the particle reconstruction in figure 5.16. The particle-based TPB density and activity data is presented in table 5.4. The particle model phase volume fractions and connectivities are roughly the same as those in the voxel reconstruction. However, the resulting TPB density is much larger, calculated as 15.8 µm−2 . This is approximately 2.7 times greater than the total voxel TPB density (5.8 µm−2 ) which is almost 1.4 times larger than the number reported by Wilson et al. [64] (4.2 µm−2 ). Connectivity and phase-volume information for the pore phase of the particle-based anode were not analyzed in detail, and therefore, the activity classification of the TPB curves is based only on the network information from the Ni and YSZ phases. It is clear, however, that the particle representation grossly exaggerates the amount of TPB in the anode sample. The set of plots on the right of figure 5.14 shows the length distributions of the different TPB curves in the particle reconstruction. Because the resolution of the original data was reduced to 42.28 nm, features smaller than this should be disregarded. Vertical lines on the plots indicate this resolution. The slope of the curve in the top plot is related to the number of TPB curves at that length. The steepest slope indicates the highest density of TPB curves. Most curves have a length less than the resolution of the data, and this is an artifact of the particle model and the method for extracting 105 TPB curves from particle intersections (e.g., one of the causes originates from the fact that particles that barely overlap result in very short TPB curves). There is a lack of TPB curves having length between 0.002 and 0.1 µm, but only points to the right of the vertical lines are considered to be valid TPB curves of the model. Even though a large number of curves have lengths less than 42.28 nm, they contribute little to the overall TPB length. This is evident in the bottom plot which does not start increasing significantly until lengths greater than 0.1 µm. Classifying the pore phase in the particle-based anode is much more difficult than its voxel counterpart. Being able to count voxels, determine their connectivity, and quickly highlight regions of interest is not a possibility in the particle anode. The only thing representing the pore space is the fact that it is the void between the particles composing the Ni and YSZ phases. In an attempt to highlight the “isolated” regions of the pore-space in the particle anode, consider the escape probability of a molecule placed within the anode sample. Each molecule moves only under the influence of Knudsen diffusion, so that after a long time, some molecules will diffuse to one of the six boundary planes of the anode (i.e., they escape), while others never make it to a boundary. For a sufficiently long time, the molecules that never escape are assumed to be trapped in an isolated pore. The difficulty in this method is determining what a sufficiently long diffusion time is, so that molecules starting deep within the anode can diffuse to the boundary. In addition, a large number of simulated molecules must be used to provide good sampling of the pore space. Figure 5.17 highlights results from the method described in the paragraph above. Each red point represents the starting position of a molecule that never intersects a boundary of the sample anode. Hence, these positions represent “isolated” regions in the anode, and a high localized density sheds some light on the size of these regions, at least visually. The ratio of trapped molecules to total number of simulated molecules provides an estimate of the predicted volume fraction of “isolated” networks. Moving from left to right in figure 5.17 shows the effect of using more simulated molecules to fill the pore space. Moving from top to bottom shows the effect of increasing the allowed diffusion time, or analogously, the maximum allowed number of surface collisions per molecule. For example, moving down the first column (so that 105 molecules were used in each simulation), we see the number of red points decrease since molecules in small pores or deep within the anode eventually diffuse to a boundary. 106 max max There is a significant drop from Ncoll = 1000 to Ncoll = 2000, and a smaller, but noticeable drop from max max = 3000. The difference between the bottom three images is minimal, so that Ncoll = 2000 to Ncoll max qualitatively, it appears as though Ncoll = 3000 provides a sufficiently long time for molecules that are not trapped to escape the anode. This finding is confirmed in figure 5.18 which provides a more quantitative representation of the data in figure 5.17. The results show an insensitivity to number of molecules used, though it should be noted that more molecules gives better sampling of the entire pore space. The percentage of trapped max > 3000 and the curve asymptotes to around 1.5%. This means molecules changes little for Ncoll that an estimated 1.5% of the pore space in the particle-based reconstructed anode is “isolated” which compares favorably to 1.35% in the voxel anode. A 3D visualization of the comparison between the red “isolated” networks of the pore space can be made between figure 5.12 (bottom) and figure 5.13. While the result from the particle model is littered with low-density scattering of red points—due in part to the method used—it is remarkable how well the larger significant volumes compare to the voxel model. 5.5 Conclusions A method similar to that used by Wilson et al. [64] was developed and extended to analyze the microstructure of a 3D-SOFC voxelated anode reconstruction. The technique allows determination of volume, surface area, and connectivity information for the networks composing each phase (Ni/metal, YSZ/oxide, pore/gas) in the anode sample. Network size distributions and their level of connectivity are broken down to provide a description of the inner architecture of the anode framework. In addition, TPB edges are identified, connected to form curves, and classified based on expected activity. The results from these calculations are compared to the same type of data from a particle-based reconstruction of the same anode. Since one cannot simply count voxels and edges for the particle model, a new method was develop to analyze the volume fractions, connectivity, and TPB information. Most of the volume and connectivity data between the two models are in agreement, but the major difference is the particle model exaggerates the TPB density by a factor of two to three. The parameters presented in this chapter are necessary to develop reliable SOFC models for numerical simulations and provide insight into previously unknown quantities. It is impractical to conclude that those in the fuel-cell community 107 and beyond need to use FIB-SEM (or other imaging techniques) and detailed 3D reconstructed models for every structure to be examined. Particle-based models provide a simple alternative; they are easy to generate, easy to use for computational purposes, and are a good representation of the actual structure, especially if using them as transport models (for all phases). The sensitivity of the particle model on particle size has a strong influence on TPB density, so these types of models are probably not as useful for reactive simulations unless corrections are made to account for the exaggerated TPB. It would be useful in the future to study more particle-based models to help reduce this sensitivity, or perhaps use an entirely different method to better represent the anode structure—level set methods, for example. Though results are presented and compared for a single anode, the methods can easily be extended to other datasets. At the time, the data provided by Wilson et al. [65] was the first of its kind, but more are expected to become available. In fact, Iwai et al. [25] just published results for a much larger set of data, and it would be interesting to use the methods presented herein to provide a more detailed picture of the 3D structure appearing in that article. 1E-5.1E-4.1E-3. 0.01 100 80 60 40 20 ´ò ´´ ò´ ´ ´ ´ 0.1 1. òç çç 10. 100. ç ç çç òç ç ´ç ç ò ò òç ò ´´ ç ò ´ ç òò Voxel anode HNiL ò´ ´´ ç ´ isolated ò ò ò ´´ ò dead-end òò ò ´ ´ ç across òò 0 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. % Networks smaller in volume % Networks smaller in volume 108 3 ´ isolated ò dead-end 60 ç across 40 20 ç çç òç ç òçòç´ç ´´´ò çç òòòò ç´´´ ´ò´´òò òòò ò ´ò ´ ´´´ 0 1E-5.1E-4.1E-3. 0.01 0.1 ç ò ç ç 1. 3 Network volume HΜm L 10. 100. 60 40 20 ò ò ´ ´ ´ ò ò ç ç ò ò´ òò ´ ´ò ò ç ´´ ò ò ò òç ´ ò ç´ ´ ò ò´ ç ò´ 0 1E-5.1E-4.1E-3. 0.01 0.1 1. ç ò 10. 100. ç òç òç çç çç Particle anode HNiL ´ isolated ò dead-end ç across 1. 10. 100. Network volume HΜm L 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. ç 100 Cumulative volume % Cumulative volume % Voxel anode HNiL 80 0.1 3 Network volume HΜm L 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. ç 100 80 1E-5.1E-4.1E-3. 0.01 100 Particle anode HNiL 80 ´ isolated ò dead-end 60 ç across 40 20 ç ò çç ç çç çò ç ò ´ ò ò ò ò ´´ ò ç ´´ ò ò òòç ç ´´ò´ ç´ ò òò ´ ò ´ 0 1E-5.1E-4.1E-3. 0.01 0.1 ò ç 1. 10. 100. 3 Network volume HΜm L Figure 5.5. Distribution of network phase volumes for Ni in the three-dimensional SOFC reconstructions. The top plots show the relationship between the separate networks (53 in the voxel anode versus 48 in the particle anode) based on the number of networks that are equal or smaller in size. The bottom plots are the cumulative volume percent. For reference, the two vertical lines on the voxel plot are drawn at the single-voxel volume of 75,560.5 nm3 (magenta) and exactly twice this (orange). 1E-5.1E-4.1E-3. 0.01 100 0.1 1. 10. 100. ç ç ´ ç 80 ò 60 ò Voxel anode HYSZL ´ ç ò 40 20 ´ isolated ò dead-end ç across ´ ´ 0 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. % Networks smaller in volume % Networks smaller in volume 109 3 100. ç Voxel anode HYSZL 80 ´ isolated ò dead-end 60 ç across 40 20 ´ç òç´ òò ç ´ ´ 0 1E-5.1E-4.1E-3. 0.01 0.1 1. 3 Network volume HΜm L 10. 0.1 1. 10. 100. ç ç ç ´ ò ´ 80 60 20 Particle anode HYSZL ò ò 40 ´ isolated ò dead-end ç ç across ç ´ 0 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. Network volume HΜm L 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100 100. 3 Cumulative volume % Cumulative volume % Network volume HΜm L 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100 1E-5.1E-4.1E-3. 0.01 100 100. ç Particle anode HYSZL 80 ´ isolated ò dead-end 60 ç across 40 20 òç ç ç ´ çò ´´ 0 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. 3 Network volume HΜm L Figure 5.6. Distribution of network phase volumes for YSZ in the three-dimensional SOFC reconstructions. The top plots show the relationship between the separate networks (11 in both cases) based on the number of networks that are equal or smaller in size. The bottom plots are the cumulative volume percent. For reference, the two vertical lines on the voxel plot are drawn at the single-voxel volume of 75,560.5 nm3 (magenta) and exactly twice this (orange). % Networks smaller in volume 110 1E-5.1E-4.1E-3. 0.01 100 80 60 40 20 0.1 çç ç ´ò ç ç ´ò òòòò ç ò ç ò ´ ç ò ´ ò ç ´ ç ´´ò ò ´ ç ´ ç ò´ ç ò ò ò ´ ´ ò ò ´ ò ò´ ò ´ò´ ´ ò ò ´´´ ç ´ ´ ´ ò ´ ò òç ´ ç ´ ç ´ ò´´ ´ ò ò ´ ò ò´ ´ òç ´ ò ´ ò´ ´ ç ò ´ ò ´ ´ ò ò ò ´ ´ ´ ´ ´ ´ ç ´ ò 0 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. ç ç Voxel anode HPoreL ´ isolated ò dead-end ç across 1. 10. 100. 3 Cumulative volume % Network volume HΜm L 1E-5.1E-4.1E-3. 0.01 0.1 1. 10. 100. ç 100 Voxel anode HPoreL 80 ´ isolated ò dead-end 60 ç across 40 20 çç ç ´ò ç ç ´ò òòòò ç ò ç ò ç´ ´ ò´ò çç ´´´ ò ´ çç ò ò ò ´´´ òò ò ò´ ò ´ò´ ´ ò ò ´´´´´ ´ ò òçç ´ ç ´ ç ´ òò´´ ò´ ò´ ´ ´ç òò ´ òç ç ´ ò ´ 0 1E-5.1E-4.1E-3. 0.01 0.1 ç 1. 10. 100. 3 Network volume HΜm L Figure 5.7. Distribution of network phase volumes for pores in the three-dimensional SOFC voxel reconstruction. The top plot shows the relationship between the 143 separate networks based on the number of networks that are equal or smaller in size. The bottom plot is the cumulative volume percent. For reference, the two vertical lines are drawn at the single-voxel volume of 75,560.5 nm3 (magenta) and exactly twice this (orange). 111 Figure 5.8. The entire Ni phase (top) and highlights of the “isolated” (red) and “dead-end” (blue) networks for the phase (bottom) in the voxel anode. 112 Figure 5.9. The entire Ni phase (top) and highlights of the “isolated” (red) and “dead-end” (blue) networks for the phase (bottom) in the particle anode. 113 Figure 5.10. The entire YSZ phase (top) and highlights of the “isolated” (red) and “dead-end” (blue) networks for the phase (bottom) in the voxel anode. 114 Figure 5.11. The entire YSZ phase (top) and highlights of the “isolated” (red) and “dead-end” (blue) networks for the phase (bottom) in the particle anode. 115 Figure 5.12. The entire pore phase (top) and highlights of the “isolated” (red) and “dead-end” (blue) networks for the phase (bottom) in the voxel anode. 116 Figure 5.13. Highlights of the “isolated” networks for the pore phase in the particle anode. The image shows starting positions of molecules that did not escape the anode after a maximum of 3,000 allowable surface collisions. There are 17,804 points from a total of 106 simulated molecules. The volume percent of “isolated” networks was calculated as 1.35% for the voxel-based anode (see table 5.1) and estimated to be around 1.5% for the particle-based anode (see figure 5.18). 0.1 80 60 40 20 0 0.01 ò ç ´ ò ç 10. 100. 100. ç Voxel anode 80 ç ç ò ç ç ò unknowndead-end çç ççç ç ç activeacross çç ç ò ç ´ ç ò ç çç ò ç ´ ç ´ ç ò ç ç ç ç ç ò ´ ç ò ç ´ ç ç ò ´ ç ò ´ ç ò ´ ´ òç ç ò ç ´ ç ò ò´ ç ´ ç ò ç ´ ò ´ ç ò ç ò ´ ò ç ´ ç ´ ò ´ ç ´ ç ò ´ ç ç ´ ç ´ ç ò ò ´ ç ò ´ ç ç ò ç ´ ´ ò ç ò ç ò ´ ´ ç ò ò ´ ç ò ç ´ ò òç ´ òç òç ´ ´ òç òç ´ òç òç ç ´ inactiveisolated 60 40 20 0 0.01 100. ççççç çòç çç ç çç ç ò ç ç´ ò ç ç ç ò ç ´ ç ´ ò ç ç ç ç ç ò ´ ç ò ç ´ ç ´ ç ò ç ç ò ´ ç ò ´ ´ òç ç ò ç ´ ç ò ò´ ç ´ ç ç ò ´ ç ò ç ´ ç ò ç ò ´ ò ç ò ´ ç ò ç ´ ò ´ ç ´ ç ò ´ ç ´ ò ´ Voxel anode ç ´ ç ´ ç ´ ç ò ´ ç ò ´ ç ò ´ ç ç ò ´ ´ inactiveisolated ´ òç ç ò òç ´ ç ò ò ´ ò òç ´ ç ò unknowndead-end ò ç ´ ò´ ò ç ç ò ç ò ´ ç ò ç ´ ç activeacross ´ ò ò ´ ç ´ ´ ò ´ ò òç ò ç òç ç ´ ç ò ò ç ç ò 0.1 1. 10. Length of curve HΜmL 0.1 1. 10. 0.01 100 Cumulative length % 1. % Curves smaller in length 0.01 100 Cumulative length % % Curves smaller in length 117 ´ ò ç ò ç 0.1 1. 10. Length of curve HΜmL 100. 1E-5.1E-4.1E-3.0.01 0.1 100 80 60 40 20 1. 10. 100. 1000. ç çç ç òò çòç ç ´ ç ò çç ´ òò ò ç ç ç ç ç ç ò ´ ç ò òç ç ´ ç ç ç ´ ç ´ ç ò ç ´ ò ç ç ç ç ò ç ´ ò ò ç ç ´ ò ç ò ç ´ ç ò ç ´ ç ò ç ´ ç ò ç ò ç ç ò ç ò ç ç ç ç òò ç ò ç ´ ç ´ ò ç ç ò ç ´ ç ò ç Particle anode ç òò ç ´ ç ò ç ò ç ò ç òç ç ò ç ´ inactiveisolated ç ´ ò ò ç ´ ò ç ç òç òç´ òç ç ç ç ç òç ´ ò unknowndead-end ç ò ç ò ò ç ç ç ò ç ò ç ´ ç ç activeacross ´ ç ò ç ç ç ò ´ ç ò ç ò ç ò ç ç ç ç ò ´ ò òç ò ´ òçòòç çç ç 0 1E-5.1E-4.1E-3.0.01 0.1 1. 10. 100. 1000. Length of curve HΜmL 1E-5.1E-4.1E-3.0.01 0.1 1. 10. 100. 1000. ç 100 Particle anode 80 ´ inactiveisolated ò unknowndead-end 60 40 20 ç activeacross ç ç çò ç ´ ç ç ò ç ´ òò ç çò ç ç ç ç ò ´ ç ò òç ç ´ ç ç ç ´ ç ´ ç ò ç ´ ò ç ç ç ç ò ç ´ ò ò ç ò ç ´ ò ç ò ç ´ ç ò ç ´ ç ò ç ´ ç ò ç ò ç ò ç ò ç ò ç ç ò ç ò ç ò ç ò ç ´ ò ç ´ ò ç ò ´ ò ç ò òç ´ ç òç òçç ç òç ç ç ç ç ´ ç ò ò ç ò ç ç ò ç ´ ò ç ç ç ò ´ ç ò ç ò ç ç ò ´ç ò ´ òç ç çç òçòòç ç ò ç ò 0 1E-5.1E-4.1E-3.0.01 0.1 1. 10. 100. 1000. Length of curve HΜmL Figure 5.14. Distribution of triple-phase boundary (TPB) length for the three-dimensional SOFC reconstructions. The top plots show the relationship between the separate curves (302 in the voxel anode versus 641 in the particle anode) based on the type and number of curves that are equal or smaller in size. The bottom plots are the cumulative percent of the overall TPB length. The total volume of the reconstructed sample is 107.65 µm3 . For reference, the two vertical lines on the voxel plot are drawn at the single-voxel dimension of 42.28 µm (magenta) and exactly twice this (orange). 118 Figure 5.15. Highlights of the “isolated”/inactive (red), “dead-end”/unknown (blue), and “across”/active (white/light gray) triple-phase boundary curves for the three-dimensional SOFC voxel anode. Figure 5.16. Highlights of the “isolated”/inactive (red), “dead-end”/unknown (blue), and “across”/active (white/light gray) triple-phase boundary (TPB) curves for the three-dimensional SOFC particle anode. The TPB classification is based only on the Ni and YSZ network connectivity, without considering the pore space. Including the pore space information would cause a shift toward an increased number of inactive curves. 119 N = 1 × 105 N = 2 × 105 N = 5 × 105 N = 1 × 106 max Ncoll = 1000 max Ncoll = 2000 max Ncoll = 3000 max Ncoll = 4000 max Ncoll = 5000 Figure 5.17. Three-dimensional visualization showing the dependence of escape probability on number max of simulated molecules N and number of maximum allowed collisions Ncoll . Each red-colored point is the starting position of a molecule that did not escape the anode sample after the prescribed number of allowable collisions. These are presumed to represent “isolated” pore regions in the particle-based anode. By choosing a sufficient number of molecules (to sample the entire phase space) and a sufficient time to move within the anode, those that never intersect a boundary plane are deemed to be trapped within an isolated pore. 120 % Molecules that are trapped 5 4 ç á ç N=1´105 á N=2´105 ó N=5´105 ´ N=1´106 3 2 ç á ´ ç ´ ó ç ç 1 0 0 1000 2000 3000 4000 5000 6000 Max allowable wall collisions Hµ residence time in anodeL Figure 5.18. Sensitivity analysis to determine a sufficient simulation time for molecules to escape the particle-based three-dimensional SOFC anode based on the results in figure 5.17. The maximum allowed collisions is proportional to the residence time of a molecule in the anode and therefore is proportional to the computational time. The results show an insensitivity to the number of simulated molecules (at least in the range tested), but more molecules gives better sampling of the entire phase space. The percentage of trapped molecules is representative of the volume percent of “isolated” networks in the pore phase. The curve asymptotes to a value near 1.5%. 121 Chapter 6 An Alternative Vision: Architecturally Engineered Multiscale Fuel Cell Structures 6.1 Summary and Recommendations Much of the following was extracted from a joint proposal submitted to the Global Climate and Energy Program (GCEP) at Stanford University by my previous advisor, David G. Goodwin (principal investigator), and another Caltech professor and member of my committee, Sossina Haile (coprincipal investigator). The major focus of the proposal was on fabrication itself. The subsequent data presented here is all original and supports the vision described below. The goal of this thesis was to provide a picture of the internal framework of a traditional solidoxide fuel cell (SOFC) anode, and based on the data presented, to offer suggestions to improve SOFC performance, or at least improve SOFC modeling capabilities. Chapter 3 described a method to create model porous media and predict Knudsen diffusivities of a gas moving within the porous network. Among the porous samples studied were two reconstructions of an actual SOFC porous anode. Both models of the modern SOFC anode reveal a highly tortuous structure that severely hinders gas transport. Based on how sensitive the particle-based reconstruction was to particle sizes and resolution, it would be interesting to extend this research by constructing a third model, perhaps using a level set method to map the surfaces of the anode. Chapter 4 and chapter 5 are focused on probing the inner networks of the porous sample. Many of the structures studied elucidate a unique distribution of pore sizes, characterized by chord distributions of molecules moving within the pores. Some distributions were 122 smooth and exponential in behavior, while others had many peaks, perhaps signifying discrete pore sizes. No conclusions were made regarding which type of distribution had the biggest influence on gas diffusion. However, all of the materials having multimodal distributions resulted in larger diffusibility predictions as compared to those with smoother distributions at similar void fractions. In addition to estimating pore-size distributions, each of the phase networks of the three-dimensional (3D) SOFC models was analyzed in detail. Data pertaining to triple-phase boundary (TPB) density and activity were presented, as well as volume, surface area, and connectivity information for each of the phases. Aside from a few recent publications which are not as comprehensive, this is one of the first known works to present such detailed data regarding traditional SOFC anodes. Many of the numbers leave a lot to be desired—as far as optimizing anode performance—so the potential for increasing SOFC performance by redesigning the anode is huge. The requirements that must be met to produce high power-density membrane-electrode assemblies (MEAs) are far beyond what can be achieved with conventional powder-based fabrication methods. If one tries to produce a high-porosity cermet with traditional manufacturing techniques, the result is low connectivity between the solid phases and therefore poor electrochemical activity; if one tries to create more surface area for electrochemistry by creating smaller pores, the result is hindered gas transport. Instead of random cermet architectures resulting from current manufacturing processes, achieving truly high-performance SOFCs will require engineered architectures that span a range of length scales—from the micron scale of the framework through which gas flows, to the nanometer scale of the catalysts. This in turn will prohibit the use of the very high processing temperatures required for sintering ceramic particles; therefore, powder-based fabrication methods must give way to ones in which structures are grown from vapor-phase or solution-phase precursors. A conceptual design of a high-current-density cell is shown in figure 6.1. A thin, dense electrolyte layer separates two 3D engineered ion-conducting oxide lattice frameworks, shown here in cross section. The lattices have “beam” thicknesses of order 1 µm to facilitate ion conduction, separated by distances of the same order, for facile gas diffusion throughout the lattice. The lattice has a regular structure with high porosity and open area, low tortuosity, and full connectivity. The oxide lattices may be tens of microns thick, or even thicker on the anode side if additional surface area for catalytic (nonelectro- H2 and CO), or for structural support. The surfaces of the oxide framework are covered with a thick carpet of nanowires 20 to 50 nm in diameter, which provides a high surface area for nanostructured catalysts. The nanowires may be short, as shown here, or might extend further throughout the 123 network, depending on how much surface area is required. They touch one another, creating multiple redundant current paths, and could even be multiply Figure 6.1. Conceptual design of an engineered, high-performance membrane-electrode assembly. Image courtesy of D. G. Goodwin, Caltech. Figure 4: Conceptual design of an engineered, high-performance MEA. chemical) chemistry is needed to preprocess the fuel (e.g., to reform hydrocarbons to H2 and CO), or for structural support. 5 The surfaces of the oxide framework are covered with a thick carpet of nanowires 20 to 50 nm in diameter, which provides a high surface area for nanostructured catalysts. The nanowires may be short, as shown here, or might extend further throughout the network, depending on how much surface area is required. They touch one another, creating multiple redundant current paths, and could even be multiply branched. The nanowires have an electronically conducting core, and are coated with an appropriate thin redox-active oxide catalyst. The coating could be continuous, as shown here, or in the form of discrete nanoparticles. Such a structure can combine high porosity for gas diffusion, high connectivity for ion conduction, and high surface area for chemistry and electrochemistry. In recent years, synthetic chemists, materials scientists, and engineers have focused considerable attention on developing methods to deliberately engineer the microstructure of materials, creating ordered structures on nanometer to micron length scales. One such structure of particular interest to this work are inverse opals (see figure 6.2). The most impressive results obtained using oxide anodes have been those of Goodenoug (Huang, 2006) as briefly 124described above. Given the very early nature of those result about the electrochemical oxidation pathways and consequently it is not possible to imm the most optimal architecture. We plan to explore using SMM as the catalyst, in additio case, we propose to use well-defined anode structures of Sr2MgMnO6-! to first addres reaction mechanisms, and with the insight gained design the ideal structure. Irrespecti details of the reaction pathway, our ability to vastly enhance the concentration of ac expected to immediately and dramatically enhance anode kinetics. 4.2 Architectures We plan to use three main architectural elements to construct high-performance, e MEAs. Oxide frameworks, composed of inverse opals and similar structures with 1 micron feature sizes, will provide high-porosity fully-connected electrode structu resistances to flow of both gas and oxygen ions, while providing structural support arrays of monodisperse nanoparticles, will be used as catalysts and seeds for CVD gro multishell nanowires, with 10-50 nm diameters, and lengths ranging from ~200 nm to as catalysts and of provide electronic conductivity. From these basic “building blo Figure 6.2. Example of inverse opals. supports SEM image Ce0.5 Zr 0.5 O2 macroporous inverse opals showing varietyand of MEA structures can be built.and description from Umeda et al. [61]. the continuous oxide framework pore network (inset). Image The Oxide Framework. The struct the MEA is the oxide “skeleton” or consists of the electrolyte layer, and th extending from the electrolyte on prototypical framework is shown in F graded inverse opal structure is used f with smaller pores near the electrolyt further away. This is done to provide flow to or from the electrochemically a the electrolyte, and to allow rapid gas the outer layers of the electrode. In Fi layers are shown in each inverse op structures, there will likely be tens or layers. Especially on the anode side, area may be needed to support nanost for fuel reforming or other no catalysis. A sufficient thickness is a Figure 6.3. TheFigure oxide framework. Image courtesy Caltech. 8: The oxide framework . of David G. Goodwin, structural stability. Variants of this basic structure are also of interest. For example, vertical vias (holes) e Inverse opals are 3D periodic honeycomb-like structures, usefresh as photonic the electrode to the electrolyte may originally be useful developed to rapidly for bring reactant gas (fuel o electrochemically active region, where it may then diffuse laterally through the sma crystals. They may be fabricated using The templates of polystyrene spheres are later burned inverse opal. ideal here would be to createthat a multiscale networkout of and pores to efficien and/or air to nanostructured catalyst surfaces. would serve very well as the oxide lattice framework for an engineered, high-performance SOFC-MEA. On the anode side, the framework will generally be fabricated of the same oxide Unlike random-particle arrays, these(e.g. structures a fully connected network at high porosities. electrolyte SDC), provide and the fabrication of the anode + electrolyte will be done in one framework will be added in a second step, and may or may not be composed of the elect A prototypical frameworkmay making use of inverse opal structures is framework shown in figure size-graded be possible to fabricate the cathode from6.3. the A catalyst BSCF, or from a and SDC. inverse opal structure is used for each electrode, with smaller pores near the electrolyte and larger ones further away. This is done to provide facile ion current flow to or from the electrochemically active region near the electrolyte, and to allow rapid gas transport through the outer layers of the electrode. Only four layers are shown, but the actual structure will likely consist of tens or even hundreds of layers. Variants of this basic structure are also of interest. For example, vertical vias (holes) extending 125 1.0 Diffusibility, ΕDDiB or DD0 Diffusibility, ΕDDiB 0.8 1.0 ç æ ç æ Inverse opals ç æ ç æ ç æ 0.6 ç æ ç æ 0.4 ç æ 0.2 0.8 0.6 ç æ ó — SC — BCC — FCC ó ç æ + — RAND ç Currie; spheres - experiment H1960L á Hoogschagen; spheres - experiment H1955L ç æó ñ Kim; spheres - experiment H1987L + õ Huizenga; spheres - experiment H1986L ç æ ó Trinh; spheres - 3D computation H2000L ç æ + Ryan; squares - 2D computation H1983L ó + æ Particle-based SOFC ç æ à Voxel-based SOFC ó —— ç Inverse opals æ ç æ+ 0.4 ó 0.2 óç á ó ñç çç ç õ õ õ ç õ õññ á á +õ + ç æ ç ç+ó 0.0 0.0 ç æ 0.2 0.4 0.6 Porosity, Ε 0.8 1.0 0.0 ó 0.0 ó æà 0.2 ç æ 0.4 0.6 Porosity, Ε 0.8 1.0 Figure 6.4. Diffusibility results for a uniform periodic simple cubic array of inverse opals (left), shown together with data from chapter 3 (right). through the electrode to the electrolyte may be useful to rapidly bring fresh reactant gas down to the electrochemically active region, where it may then diffuse laterally through the smaller pores of the inverse opal. The ideal here would be to create a multiscale network of pores to efficiently distribute fuel and/or air to nanostructured catalyst surfaces. 6.2 A First Glance Inside a New Design Consider the simple case of uniform inverse opals stacked in all directions and varying the amount of overlap to achieve different porosities. The greater the allowed overlap, the higher the porosity (for optimized gas diffusion), but less lattice material is then available for mechanical support or to grow catalyst nanowires (necessary for power-producing reactions). Diffusibility results are shown in figure 6.4. Comparison to current technology (shown as the particle- and voxel-based SOFC) shows significant improvement in predicted diffusibility and dramatically lower tortuosity, as expected. This is because inverse opal structures allow much higher porosity without losing material connectivity and structural support. These results were obtained for a simple cubic (SC) packing of inverse opals. Generally, face-centered cubic (FCC) packing is obtained experimentally—FCC has twelve nearest neighbors rather than six for SC. Using FCC would potentially result in even higher predicted diffusibilities and lower tortuosities since there are more opportunities to escape a given pore. Table 6.1 presents these results for four inverse opal geometries (achieved by varying only the 126 Table 6.1. Comparison of Knudsen transport data for inverse opals and the three-dimensional SOFC reconstructed anodes Particle 17.5 Voxel 19.86 Opal 1 55.6 Opal 2 60.4 Opal 3 63.7 Opal 4 68.7 Mean chord length,a h`i (µm) Mean-squared chord length, h`2 i (µm2 ) Ratio, h`2 i/(2h`i2 ) 0.192 0.092 1.251 0.261 0.137 1.000 1.415 2.386 0.596 1.552 3.119 0.648 1.651 3.744 0.687 1.822 4.772 0.719 β of equation (4.4)b 0.442 0.499 0.474 0.448 0.427 0.398 Diffusibility,c D/DiB Tortuosity, τ 0.027 6.42 0.039 5.05 0.306 1.82 0.469 1.29 0.508 1.25 0.594 1.16 Porosity (%) The mean chord length gives an estimate of the mean pore radius, R, as appearing in equation (2.31) from the relation h`i = 2R. The opal radius, R, is 1 µm and h`i/R for a closed sphere is 4/3 (see §4.4.2). b The parameter β provides a measure of the nature of redirecting collisions at gas-solid interfaces. It is discussed in chapter 4. c Der The parameter D is defined by equation (2.28), and DiB in equation (2.32) using DiK = DiK from equation (4.3). a amount of allowable overlap). The predicted diffusibility for inverse opal architectures is an order of magnitude higher than conventional anodes produced from powder-based fabrication methods. Further, the tortuosity is reduced by a factor of three to five (which, of course, contributes to the increased diffusibility). All of this is achieved by allowing less than 10% overlap (based on the opal radius) between neighboring opals, so there remains a great deal of lattice structure for necessary support and on which to grow carpets of nanowire catalysts necessary for power-producing reactions. As more opal overlap is allowed, the network moves from one of connected spherical pores, to one that resembles a 3D network of perpendicular intersecting cylinders. Increasing overlap shifts the pore-size distribution from that obtained for a single closed sphere (see figure 4.4) to one with multiple peaks. Each peak symbolizes a new transport pathway for a diffusing molecule as the increased opening between inverse opals allows it to move more freely from one inverse opal to the next. Small necking areas between larger pores provide connections for a molecule to move from one pore to the next. 6.3 Concluding Remarks Gas transport analysis revealed that the coefficient of diffusion through a porous medium is a function of both internal geometry and porosity. Gas diffusion in traditional SOFC electrode designs is greatly 127 hindered by the complexity of the pore network. Current methods for producing such electrodes (hightemperature sintering of ceramic powders) result in random convoluted networks in the SOFC anode. These networks are made up of submicron pores which limit the ability of the gas to move freely. It is impractical to conclude that those in the fuel-cell community and beyond need to use FIBSEM (or other imaging techniques) and detailed 3D reconstructed models for every structure to be examined. Particle-based models provide a simple alternative; they are easy to generate, easy to use for computational purposes, and are a good representation of the actual structure, especially if using them as transport models (for all phases). The sensitivity of the particle model on particle size has a strong influence on TPB density, so these types of models are probably not as useful for reactive simulations unless corrections are made to account for the exaggerated TPB. It would be useful in the future to study more particle-based models to help reduce this sensitivity, or perhaps use an entirely different method to better represent the anode structure—level set methods, for example. In any case, fuel cells of the future should take advantage of the opportunity to precisely engineer the components of the MEA. Perhaps the only limit will be our own imaginations and ability to innovate any structure we choose. As micro- and nanofabrication processes improve and expand, so too will the possibilities for engineered structures. This opens the door to an expanse of exciting new research in fuel cell technology and energy in general. 128 Bibliography [1] M. Abbasi, J. Evans, and I. Abramson. Diffusion of gases in porous solids: Monte Carlo simulations in the Knudsen and ordinary diffusion regimes. AIChE J., 29(4):617–624, July 1983. [2] P. Asinari, M. Quaglia, M. von Spakovsky, and B. Kasula. Direct numerical calculation of the kinematic tortuosity of reactive mixture flow in the anode layer of solid oxide fuel cells by the lattice Boltzmann method. J. Power Sources, 170:359–375, 2007. [3] G. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Clarendon Press, 2nd edition, 1994. [4] R. Bird, W. Stewart, and E. Lightfoot. Transport Phenomena. John Wiley & Sons, 2nd edition, 2002. [5] S. Bogart, K. Schulz, L. Brown, and B. Russ. Production of liquid synthetic hydrocarbon fuels from air, water, and nuclear power on ships and at shore bases for military use. Reno, NV, 2006. International Congress on Advances in Nuclear Power Plants. [6] C. Bosanquet. Technical report, British TA Rept. BR-507, 1944. [7] V. Burganos. Gas diffusion in random binary media. J. Chem. Phys., 109(16):6772–6779, October 1998. [8] V. Burganos and S. Sotirchos. Knudsen diffusion in parallel, multidimensional or randomly oriented capillary structures. Chem. Eng. Sci., 44(11):2451–2462, 1989. [9] S. Chapman and T. Cowling. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, 3rd edition, 1970. 129 [10] W. Chiu, A. Joshi, and K. Grew. Lattice Boltzmann model for multi-component mass transfer in a solid oxide fuel cell anode with heterogeneous internal reformation and electrochemistry. Eur. Phys. J. Special Topics, 171:159–165, 2009. [11] P. Clausing. The flow of highly rarefied gases through tubes of arbitrary length. J. Vac. Sci. Tech., 8(5):636–646, 1971. [12] M. Coppens and K. Malek. Dynamic Monte-Carlo simulations of diffusion limited reactions in rough nanopores. Chem. Eng. Sci., 58:4787–4795, 2003. [13] J. Currie. Gaseous diffusion in porous media: Part 1.–A non-steady state method. Brit. J. Appl. Phys., 11:314–317, August 1960. [14] J. Currie. Gaseous diffusion in porous media: Part 2.–Dry granular materials. Brit. J. Appl. Phys., 11:318–324, August 1960. [15] D. Davis. Monte Carlo calculation of molecular flow rates through a cylindrical elbow and pipes of other shapes. J. Appl. Phys., 31(7):1169–1176, July 1960. [16] B. Derjaguin. Measurement of the specific surface of porous and disperse bodies by their resistance to the flow of rarefied gases. Prog. Surf. Sci., 45:337–340, 1946. [17] J. Evans, M. Abbasi, and A. Sarin. A Monte Carlo simulation of the diffusion of gases in porous solids. J. Chem. Phys., 72(5):2967–2973, March 1980. [18] R. Evans III, G. Watson, and E. Mason. Gaseous diffusion in porous media at uniform pressure. J. Chem. Phys., 35(6):2076–2083, December 1961. [19] R. Evans III, G. Watson, and E. Mason. Gaseous diffusion in porous media. II. Effect of pressure gradients. J. Chem. Phys., 36(7):1894–1902, April 1962. [20] D. Gostovic, J. Smith, D. Kundinger, K. Jones, and E. Wachsman. Three-dimensional reconstruction of porous LSCF cathodes. Electrochem. Solid-State Letters, 10(12):B214–B217, 2007. [21] J. Greenwood. The correct and incorrect generation of a cosine distribution of scattered particles for Monte-Carlo modelling of vacuum systems. Vacuum, 67:217–222, 2002. 130 [22] M. Hollewand and L. Gladden. Modelling of diffusion and reaction in porous catalysts using a random three-dimensional network model. Chem. Eng. Sci., 47(7):1761–1770, 1992. [23] J. Hoogschagen. Diffusion in porous catalysts and adsorbents. Ind. Eng. Chem., 47(5):906–912, May 1955. [24] D. Huizenga and D. Smith. Knudsen diffusion in random assemblages of uniform spheres. AIChE J., 32(1):1–6, January 1986. [25] H. Iwai, N. Shikazono, T. Matsui, H. Teshima, M. Kishimoto, R. Kishida, D. Hayashi, K. Matsuzaki, D. Kanno, M. Saito, H. Muroyama, K. Eguchi, N. Kasagi, and H. Yoshida. Quantification of SOFC anode microstructure based on dual beam FIB-SEM technique. J. Power Sources, 195:955–961, 2010. [26] J. Izzo Jr., A. Joshi, K. Grew, W. Chiu, A. Tkachuk, S. Wang, and W. Yun. Nondestructive reconstruction and analysis of SOFC anodes using X-ray computed tomography at sub-50 nm resolution. J. Electrochem. Soc., 155(5):B504–B508, 2008. [27] S. J. Jeans. Kinetic Theory of Gases. Cambridge University Press, 1st edition, 1940. [28] A. Joshi, A. Peracchio, K. Grew, and W. Chiu. Lattice Boltzmann method for continuum, multicomponent mass diffusion in complex 2D geometries. J. Phys. D: Appl. Phys., 40:2961–2971, 2007. [29] A. Joshi, A. Peracchio, K. Grew, and W. Chiu. Lattice Boltzmann method for multi-component, non-continuum mass diffusion. J. Phys. D: Appl. Phys., 40:7593–7600, 2007. [30] R. Kee, M. Coltrin, and P. Glarborg. Chemically Reacting Flow: Theory and Practice. John Wiley & Sons, 2003. [31] R. Kee, H. Zhu, and D. Goodwin. Modeling electrochemistry and solid-oxide fuel cells: I. Basic principles. J. Combust. Soc. Jpn., 47(141):192–204, 2005. [32] E. Kennard. Kinetic Theory of Gases. McGraw-Hill Book Company, 1st edition, 1938. 131 [33] J. Kim, J. Ochoa, and S. Whitaker. Diffusion in anisotropic porous media. Trans. Por. Med., 2:327–356, 1987. [34] M. Knudsen. Die gesetze der molekularströmung und der inneren reibungsströmung der gase durch röhren. Ann. Physik, 28:75–130, 1909. [35] K. Kolasinski. Surface Science. John Wiley & Sons, 1st edition, 2002. [36] P. Levitz. Knudsen diffusion and excitation transfer in random porous media. J. Phys. Chem., 97:3813–3818, 1993. [37] P. Levitz. Off-lattice reconstruction of porous media: critical evaluation, geometrical confinement and molecular transport. Adv. Coll. Int. Sci., 76-77:71–106, 1998. [38] B. Lu and S. Torquato. Chord-length and free-path distribution functions for many-body systems. J. Chem. Phys., 98(8):6472–6482, April 1993. [39] K. Malek and M. Coppens. Effects of surface roughness on self- and transport diffusion in porous media in the knudsen regime. Phys. Rev. Letters, 87(12):125505:1–4, September 2001. [40] K. Malek and M. Coppens. Knudsen self- and Fickian diffusion in rough nanoporous media. J. Chem. Phys., 119(5):2801–2811, August 2003. [41] E. Mason, R. Evans III, and G. Watson. Gaseous diffusion in porous media. III. Thermal transpiration. J. Chem. Phys., 38(8):1808–1826, April 1963. [42] E. Mason and A. Malinauskas. Gas transport in Porous Media: The Dusty-Gas Model, volume 17 of Chem. Eng. Monographs. Elsevier, 1983. [43] K. Meyerhoff, B. Schnorr, and D. Hesse. Characterization of pore texture and determination of effective Knudsen diffusion coefficients by digital image processing. Hung. J. Ind. Chem., 22:227– 232, 1994. [44] C. Miller and S. Torquato. Diffusion-controlled reactions among spherical traps: Effect of polydispersity in trap size. Phys. Rev. B, 40(10):7101–7108, October 1989. [45] R. Millington. Gas diffusion in porous media. Science, 130:100–102, 1959. 132 [46] N. Olague, D. Smith, and M. Ciftcioglu. Knudsen diffusion in ordered sphere packings. AIChE J., 34(11):1907–1909, November 1988. [47] J. Petrasch, P. Wyss, and A. Steinfeld. Tomography-based Monte Carlo determination of radiative properties of reticulate porous ceramics. J. Quantitative Spectroscopy & Radiative Transfer, 105:180–197, 2007. [48] W. Pollard and R. Present. On gaseous self-diffusion in long capillary tubes. Phys. Rev., 73(7):762– 774, April 1948. [49] R. Present. Kinetic Theory of Gases. McGraw-Hill Book Company, 1st edition, 1958. [50] S. Reyes and E. Iglesia. Effective diffusivities in catalyst pellets: New model porous structures and transport simulation techniques. J. Catal., 129:457–472, 1991. [51] S. Reyes and E. Iglesia. Monte Carlo simulations of structural properties of packed beds. Chem. Eng. Sci., 46(4):1089–1099, 1991. [52] M. Riley, F. Muzzio, H. Buettner, and S. Reyes. Diffusion in heterogeneous media: Application to immobilized cell systems. AIChE J., 41:691–700, 1995. [53] Z. Shao, J. Mederos, W. Chueh, and S. Haile. High power-density single-chamber fuel cells operated on methane. J. Power Sources, 162:589–596, 2006. [54] S. Singhal and K. Kendall. High Temperature Solid Oxide Fuel Cells. Elsevier, 2003. [55] D. Smith and D. Huizenga. Simulation of Knudsen diffusion in porous media. 10th IASTED Symp. Appl. Mod. Sim., 1984. [56] H. Spacil. US Patent 3,558,360; filed October 30, 1964, modified November 2, 1967, granted March 31, 1970., 1964. [57] W. Steckelmacher and M. Lucas. Gas flow through a cylindrical tube under free molecular conditions. J. Phys. D: Appl. Phys., 16:1453–1460, 1983. [58] M. Tomadakis and S. Sotirchos. Effective Knudsen diffusivities in structures of randomly overlapping fibers. AIChE J., 37(1):74–86, January 1991. 133 [59] F. Transvalidou and S. Sotirchos. Effective diffusion coefficients in square arrays of filament bundles. AIChE J., 42(9):2426–2438, September 1996. [60] S. Trinh, P. Arce, and B. Locke. Effective diffusivities of point-like molecules in isotropic porous media by Monte Carlo simulation. Trans. Por. Med., 38:241–259, 2000. [61] G. Umeda, W. Chueh, L. Noailles, S. Haile, and B. Dunn. Inverse opal ceria-zirconia: architectural engineering for heterogeneous catalysis. Energy Environ. Sci., 1:484–486, 2008. [62] H. Weissberg. Effective diffusion coefficient in porous media. J. Appl. Phys., 34(9):2636–2639, September 1963. [63] J. Wilson and S. Barnett. Solid oxide fuel cell Ni-YSZ anodes: Effect of composition on microstructure and performance. Electrochem. Solid-State Letters, 11(10):B181–B185, 2008. [64] J. Wilson, M. Gameiro, K. Mischaikow, W. Kalies, P. Voorhees, and S. Barnett. Three-dimensional analysis of solid oxide fuel cell Ni-YSZ anode interconnectivity. Microsc. Microanal., 15:71–77, 2009. [65] J. Wilson, W. Kobsiriphat, R. Mendoza, H. Chen, J. Hiller, D. Miller, K. Thornton, P. Voorhees, S. Adler, and S. Barnett. Three-dimensional reconstruction of a solid-oxide fuel-cell anode. Nature Materials, 5:541–544, July 2006. [66] G. Youngquist. Symposium on flow through porous media: Diffusion and flow of gases in porous solids. Ind. Eng. Chem., 62(8):52–63, 1970. [67] J. Zalc, S. Reyes, and E. Iglesia. Monte-Carlo simulation of surface and gas phase diffusion in complex porous structures. Chem. Eng. Sci., 58:4605–4617, 2003. [68] J. Zalc, S. Reyes, and E. Iglesia. The effects of diffusion mechanism and void structure on transport rates and tortuosity factors in complex porous structures. Chem. Eng. Sci., 59:2947–2960, 2004.