Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Three phase boundary length and effective diffusivity in modeled sintered composite solid oxide fuel… Metcalfe, Thomas Craig 2008

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2009_spring_metcalfe_thomas.pdf [ 3.24MB ]
Metadata
JSON: 24-1.0066851.json
JSON-LD: 24-1.0066851-ld.json
RDF/XML (Pretty): 24-1.0066851-rdf.xml
RDF/JSON: 24-1.0066851-rdf.json
Turtle: 24-1.0066851-turtle.txt
N-Triples: 24-1.0066851-rdf-ntriples.txt
Original Record: 24-1.0066851-source.json
Full Text
24-1.0066851-fulltext.txt
Citation
24-1.0066851.ris

Full Text

 THREE PHASE BOUNDARY LENGTH AND EFFECTIVE DIFFUSIVITY IN MODELED SINTERED COMPOSITE SOLID OXIDE FUEL CELL ELECTRODES  by THOMAS CRAIG METCALFE BSc. Mechanical Engineering, University of Alberta, 2003  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2008 © Thomas Craig Metcalfe, 2008   ii Abstract Solid oxide fuel cells with graded electrodes consisting of multiple composite layers yield generally lower polarization resistances than single layer composite electrodes. Optimization of the performance of solid oxide fuel cells with graded electrode composition and/or microstructure requires an evaluation of both the three phase boundary length per unit volume and the effective diffusion coefficient in order to provide insight into how these properties vary over the design space. A numerical methodology for studying the three phase boundary length and effective diffusivity in composite electrode layers with controlled properties is developed. A three dimensional solid model of a sintered composite electrode is generated for which the mean particle diameter, composition, and total porosity may be specified as independent variables. The total three phase boundary length for the modeled electrode is calculated and tomographic methods are used to estimate the fraction of this length over which the electrochemical reactions can theoretically occur. Furthermore, the open porosity of the modeled electrode is identified and the effective diffusion coefficient is extracted from the solution of the concentration of the diffusing species within the open porosity. Selected example electrode models are used to illustrate the application of the methods developed, and the resulting connected three phase boundary length and diffusion coefficients are compared. A significant result is the need for thickness-specific effective diffusivity to be determined, rather than the general volume averaged property, for electrodes with porosity between the upper and lower percolation thresholds. As the demand for current increases, more of the connected three phase boundaries become active, and therefore a greater fraction of the electrode layer is utilized for a given geometry, resulting in a higher apparent effective diffusivity compared to the same electrode geometry operating at a lower current. The methods developed in this work may be used within a macroscopic electrode performance model to investigate optimal designs for solid oxide fuel cell electrodes with stepwise graded composition and/or microstructure.    iii Table of Contents Abstract ...............................................................................................................................ii Table of Contents ...............................................................................................................iii List of Figures ....................................................................................................................vi List of Tables ......................................................................................................................xi Nomenclature ....................................................................................................................xii Acknowledgements ..........................................................................................................xvi Chapter 1 Introduction ........................................................................................................1 1.1 Solid Oxide Fuel Cells ..................................................................................................1 1.1.1 Role in Society ......................................................................................................1 1.1.2 Operating Principle ................................................................................................2 1.1.3 Outlook ..................................................................................................................3 1.2 Solid Oxide Fuel Cells with Functionally Graded Electrodes ........................................4 1.3 Electrode Performance, Composition and Structure.....................................................7 1.3.1 Charge and Mass Transport in Solid Oxide Fuel Cell Electrodes ...........................8 1.3.2 Electrode Performance ........................................................................................10 1.3.3 Three Phase Boundary Length ............................................................................12 1.3.4 Diffusion of Gases in Porous Electrodes .............................................................18 1.4 Scope of Thesis .........................................................................................................20 Chapter 2 Geometric Modeling of Sintered Composite Electrodes ...............................21 2.1 Introduction ................................................................................................................21 2.2 Definition of Design Space .........................................................................................21 2.3 Description of Face Centered Cubic Lattice Parameters ............................................23 2.4 Description of Sintered Particle Intersection ...............................................................28   iv 2.5 Dispersion Study for Variation in Three Phase Boundary Length by Varying Number of Lattice Sites .....................................................................................................................31 2.6 Parametric Solid Model of Sintered Composite Electrodes ........................................34 2.7 Calculation of Total Porosity ......................................................................................37 2.8 Verification of Total Porosity Calculation ....................................................................39 Chapter 3 Three Phase Boundary Length in Sintered Composite Electrodes ..............41 3.1 Introduction ................................................................................................................41 3.2 Total Three Phase Boundary Length Calculation Method ..........................................41 3.3 Verification of Total Geometric Three Phase Boundary Length ..................................48 3.4 Identification of Connected Three Phase Boundary Voxels ........................................49 3.5 Estimation of Connected Three Phase Boundary Length ...........................................55 3.6 Required Binary Cross Section Image Resolution ......................................................56 Chapter 4 Modeling Effective Diffusivity in Sintered Composite Electrodes ................61 4.1 Introduction ................................................................................................................61 4.2 Image Generation and Calculation of Open Porosity ..................................................64 4.3 Finite Element Mesh Generation ................................................................................66 4.4 Boundary Conditions..................................................................................................69 4.4.1 Distribution of Faradaic Current ...........................................................................71 4.5 Evaluation of Effective Diffusion Coefficient ...............................................................74 4.6 Finite Element Mesh Refinement ...............................................................................78 4.7 Verification .................................................................................................................80 Chapter 5 On the Connected Three Phase Boundary Length and Effective Diffusion Coefficient in Sintered Composite Electrodes ................................................................84 5.1 Introduction ................................................................................................................84 5.2 Connected Three Phase Boundary Length ................................................................84   v 5.2.1 Example Modeled Electrodes ..............................................................................84 5.2.2 Comparison to Predictions using Models from the Literature ...............................92 5.3 Effective Diffusivity in Sintered Composite Electrodes ...............................................94 5.3.1 Example Electrode Models ..................................................................................94 5.3.2 Comparison of Calculated Tortuosity Factor to Measured Values Available in the Literature. ................................................................................................................... 102 Chapter 6 Conclusions and Future Work ...................................................................... 105 6.1 Summary ................................................................................................................. 105 6.2 Recommendations for Future Work ......................................................................... 107 References ....................................................................................................................... 109 Appendix A Computer Code for Electrode Geometry Modeling .................................. 112 Appendix B Computer Code for Calculation of Total Three Phase Boundary Length and Total Porosity ........................................................................................................... 135 Appendix C Computer Code for Calculation of Effective Diffusion Coefficient .......... 149    vi List of Figures Figure 1-1: Exploded view of a single planar solid oxide fuel cell showing anode, electrolyte, cathode and interconnect. Multiple repeating units are combined to form a stack. ........................ 2 Figure 1-2: Schematic of ion conducting particles (white) and electron conducting particles (gray) in a composite cathode near the electrolyte surface. Electrons from the external circuit are conducted through the electron conducting particles to the reaction site, where oxygen ions result from the oxygen reduction reaction. ............................................................... 5 Figure 1-3: Contact between electron and ion conducting particles. Schematic is redrawn based on [29]. ............................................................................................................................. 13 Figure 1-4: Connected three phase boundary length predicted using equation (11), where the average coordination number is a function of porosity according to work by Nakagaki and Sunada [34]. The result predicts increasing three phase boundary length as porosity decreases to zero. ............................................................................................................................. 17 Figure 2-1: Face centered cubic packing of hard spheres illustrating the relationship between edge length (lattice constant) and particle diameter. .............................................................. 23 Figure 2-2: Face centered cubic lattice based on reference particle diameter and reference contact angle, θ⁰ = 12pi radians. ............................................................................................. 25 Figure 2-3: Cross section of face centered cubic packing of monosized spheres with minimum diameter. ...................................................................................................................... 26 Figure 2-4: Face centered cubic packing with maximum particle diameters, corresponding to fully dense packing ............................................................................................................... 27 Figure 2-5: Intersection of two arbitrary particles, i and j, showing the shape defined by N,M,Q, and R to be revolved about line AN to create sinter neck on particle i. ................................... 29 Figure 2-6: Geometry of intersection between particles i and j. ........................................................ 30 Figure 2-7: Dispersion of three phase boundary length for different face centered cubic lattice sizes. Mean particle diameter for Case 1 and Case 3 is 2µm. Mean particle diameter for Case 2 is 2.9µm. ...................................................................................................................... 32 Figure 2-8: Dispersion in porosity for different face centered cubic lattice sizes. Mean particle diameter is 2µm in both cases. ...................................................................................... 34 Figure 2-9: Correlation of fraction of particles to remove from face centered cubic lattice to total porosity for three mean particle diameters. There are a total of 7813 lattice sites for each data point. .................................................................................................................... 36 Figure 2-10: Example solid model of sintered, composite, solid oxide fuel cell electrode. The total porosity is 0.55. The mean particle diameter is 7.43 µm. The electronic conductor (black particles) solid volume fraction is 0.27. ......................................................................... 37 Figure 2-11: Sinter neck arc is revolved about x axis to form volume of revolution for purpose of porosity calculation. ..................................................................................................... 38 Figure 2-12: Spherical cap for particle i. .......................................................................................... 39 Figure 2-13: Geometry used to verify total porosity calculation. The envelope volume is 371µm3 .... 40   vii Figure 3-1: Definition of sinter neck locations. 1,QXYI  represents the sinter neck located in the XY plane, quadrant one, which will only intersect with location 3,QXYI  on the neighbour particle if the latter is present on the lattice. Therefore, each sinter neck location on particle i has a unique corresponding sinter neck location for each of the possible twelve neighbour particles. ...................................................................................................... 42 Figure 3-2: Arbitrary particle with four neighbours showing interference between two of the sinter necks. ........................................................................................................................... 43 Figure 3-3: Interfering sinter necks represented by parametric curves in three dimensional space. The curve with origin at point O has four interfering sinter necks labeled a, b, c, and d. Curve a may also interfere with curve c; likewise, curve b may also interfere with curve d. Curve a cannot interfere with curve b and curve c cannot interfere with curve d. Thick lines are three phase boundary lines. ............................................................................. 44 Figure 3-4: Calculation of angle of overlap between sinter neck with origin at O and interfering sinter neck curves b and d. Curve b has radius bsr , and likewise, curve d has radius dsr , , representing the radii of the sinter necks at these locations. ........................................... 45 Figure 3-5: Example of particles occupying edge, face, and corner lattice sites. ............................... 47 Figure 3-6: Three phase boundary test geometries. (a) Single FCC unit cell. (b) Eight FCC unit cells. ..................................................................................................................................... 48 Figure 3-7: First cross section of geometry in Figure 3-6 (b). (a) Actual cross section image as generated using Unigraphics NX5. Electronic conductor is black. Ionic conductor is gray. Pore is white. (b) Binary image representing electronic conducting phase as white. (c) Binary image representing ionic conducting phase as white. ......................................... 50 Figure 3-8: Three phase boundary voxels for geometry in Figure 3-6 (a). The noisy voxels that do not fall along the circumference of intersection of opposite phase particles require filtering. 51 Figure 3-9: Three phase boundary voxels with noise removed. (a) Three phase boundary voxels to be removed, blue cubes, do not fall on the circumference of intersection between opposite phase particles, indicated by the red circles. (b) Three phase boundary voxels remaining after filtering. Note that the images in (a) and (b) are flipped vertically compared to Figure 3-8 and Figure 3-10. .......................................................................................... 52 Figure 3-10: Connected three phase boundary voxels for geometry shown in Figure 3-6 (a). ............ 53 Figure 3-11: (a) Three phase boundary voxels superimposed on a close up image of the geometry in (b). Both images are in the same orientation as Figure 3-9. ........................................... 53 Figure 3-12: Example scenario in which two three-phase boundary voxels are created where only one physically exists. In (a), E represents the electron conductor phase, I represents the ion conductor phase, and P represents the pore phase. (b) through (d) are the binary representations of (a) for the electron conductor phase, ion conductor phase, and pore phase, respectively. (e) and (f) are the expanded matrices in which ones are placed around the electronic conducting phase in (e) and ionic conducting phase in (f). (g) is the expansion matrix for the pore phase; however, note that ones are not placed around the pore phase. (h) is generated by multiplying matrices (e), (f) and (g) term by term.......... 54 Figure 3-13: A cross section image of the geometry used in Chapter Four, shown here at two different image resolutions. In (a), the resolution is 100 by 100 pixels. The two ionic conducting   viii particles (gray) enclosed in the dashed circle are separated by only one pixel. In (b), the image resolution is 300 by 300 pixels. The same two ionic conducting particles are enclosed in a dashed circle and are separated by six pixels. These images section the particles at their major diameter. ................................................................................... 57 Figure 4-1: Example modeled sintered composite cathode having properties summarized in Table 4-1. Black particles represent electronic conductor and gray particles represent ionic conductor. .................................................................................................................... 61 Figure 4-2: Example cross section image of geometry in Figure 4-1, used for determining open porosity. (a) Modeled geometry cross section where black particles represent electron conducting solid phase, gray particles represent ion conducting solid phase and white represents pore. (b) Electron conducting solid phase represented by white pixels. (c) Ion conducting solid phase represented by white pixels. ...................................................... 65 Figure 4-3: Discretized face centered cubic unit cell showing index of arbitrary pore voxel shown in red wireframe having index [5, 1, 1]. Note that pore voxels in this image are not shown and are instead represented by the absence of the either of the two solid phases. Black voxels are the electron conducting phase. Grey voxels are the ion conducting phase. .... 66 Figure 4-4: Node configuration for single open pore voxel (finite element). ..................................... 67 Figure 4-5: (a) First 25000 mesh elements showing disconnected elements at location of gas channel. The full mesh contains 348101 elements. (b) Same 25000 mesh elements as in (a) but with an additional layer of 10000 elements inserted to simulate the gas channel and to connect all pore networks. The direction of increasing z is the direction of the electrolyte. The direction of decreasing z is the direction of the current collector and gas channel. .. 68 Figure 4-6: (a) Entire mesh consisting of 556375 nodes (for the initially linear elements), shown as blue dots, representing the open porosity of the modeled cathode in Figure 4-1. The Dirichlet boundary condition is applied to all nodes in the xy plane at z equal to zero. The locations of Neumann boundary conditions, shown as red dots, are applied to the boundary surfaces containing a connected three phase boundary voxel. All remaining boundary surfaces are insulated. ................................................................................... 70 Figure 4-7: Faradaic current distribution within LSM / YSZ composite cathode depicted in Figure 4-1. Standard electrolyte and anode properties are used as depicted in Table 4-1. At z/L equal to one is the cathode – electrolyte interface. At z/L equal to zero is the cathode – gas channel interface. L is the cathode thickness, equal to 95µm. .................................. 73 Figure 4-8: Oxygen concentration for modeled cathode in Figure 4-1. A position of zero represents the location of the gas channel. A position of 95x10-6m represents the interface with the electrolyte.  Inset: zoom in on concentration profile showing a difference of  3.6x10- 6mol/m3 in the concentration between the gas channel and the cathode-electrolyte interface. ...................................................................................................................... 75 Figure 4-9: Variation in superficial diffusive flux through the thickness of the modeled electrode. ... 76 Figure 4-10: Curves N1 and N2 fitted to the superficial diffusive flux. Curve N1 is constant on the interval [z0, z1] and equal to the mean of the superficial diffusive flux on this interval. Curve N2 is a fourth order polynomial fitted to the superficial diffusive flux on the interval [z1, z2]. ............................................................................................................. 77 Figure 4-11: Image resolution of (a) 100 by 100 pixels, (b) 76 by 76 pixels, and (c) 50 by 50 pixels. 79   ix Figure 4-12: (a) Solid model of verification geometry. The height in the z direction is 14mm, and the width in the x and y directions are also 14µm. (b) Mesh of pore space for verification geometry, including layer of elements representing gas channel. Mesh contains 60480 quadratic brick elements including the 3136 elements representing the gas channel....... 81 Figure 4-13: Numerical solution for oxygen concentration in verification geometry. The concentration at z equal to 14µm is 2.2735mol/m3. ............................................................................. 83 Figure 5-1: Example modeled sintered composite electrodes. The properties for each of the five example geometries are summarized in Table 5-1. ........................................................ 85 Figure 5-2: (a) Total number of three phase boundary voxels for example electrode I. (b) Connected three phase boundary voxels for example electrode I. An example region of disconnected three phase boundary voxels is shown in the dashed circle. The gas channel / current collector is the xy plane at z equal to zero. .................................................................... 86 Figure 5-3: Example electrode I with the ionic conducting particles suppressed to reveal the electronic conducting particles. (a) All electronic conducting particles present. (b) Only connected electronic conducting particles present. The gas channel / current collector is the xy plane at z equal to zero. ......................................................................................................... 87 Figure 5-4: Connected three phase boundary voxels of example electrode II. The gas channel / current collector is the xy plane at z equal to zero. .................................................................... 89 Figure 5-5: Connected three phase boundary voxels for (a) example electrode III and (b) example electrode IV. The electrolyte is the xy plane at z equal to zero. ..................................... 90 Figure 5-6: Open porosity over the range from 0.05 total porosity to 0.55 total porosity. The lower percolation threshold appears at approximately 10% total porosity and the upper percolation threshold appears at approximately 20% total porosity. Results obtained by creating fifteen solid electrode models, with total porosity as shown, and sectioning with 100 images, each image having a resolution of 100 by 100 pixels. The corresponding data points for example electrodes I through V are also shown for reference......................... 91 Figure 5-7: Finite element mesh nodes, shown as blue dots, and locations of connected three phase boundaries, shown as red dots, for example electrode II. The gas channel is the xy plane at z equal to zero. The electrolyte interface is the xy plane at the maximum z coordinate. ..................................................................................................................................... 96 Figure 5-8: (a) Superficial oxygen diffusive flux for example electrode II. The value of the superficial flux decreases linearly from the gas channel interface to the electrolyte interface. (b) Cross section area of the pore space in the xy plane through the thickness of the electrode, z. .................................................................................................................. 97 Figure 5-9: Nodes, shown as blue dots, and location of Neumann boundary conditions, shown as red dots, for example electrode II. The gas channel is the xy plane at z equal to zero. The electrolyte interface is the xy plane at the maximum z coordinate. ................................ 98 Figure 5-10: Mesh nodes, shown as blue dots, and location of Neumann boundary conditions, shown as red dots, for example electrode V. (b) Corresponding finite element mesh. The element edges are not shown...................................................................................................... 99 Figure 5-11: Mean superficial diffusive flux for example electrode V. ........................................... 100   x Figure 5-12: Average concentration in xy plane through thickness of example electrode V.  Inset: zoom in on concentration profile showing difference of 1.6x10-3 mol/m3 in the concentration between the gas channel and the cathode-electrolyte interface. .............. 100 Figure 5-13: Mesh nodes, shown as blue dots, and locations of Neumann boundary conditions, shown as red dots, for alternate analysis of example electrode V. The Neumann boundary conditions are only applied in the xy plane at the maximum z coordinate and in the z direction. .................................................................................................................... 101 Figure 5-14: Comparison of tortuosity factor calculated in this work to the tortuosity factor measured by Williford et al. using three different measurement techniques: Mercury porosimetry and the Wike-Kallenbach method using CO2 in N2 and also H2 in N2 [20] ................... 103    xi List of Tables Table 2-1: Composite Electrode Design Space................................................................................. 21 Table 2-2: Summary of results for verification of total porosity calculation...................................... 40 Table 3-1: Summary of Total Geometric Three Phase Boundary Length Verification....................... 49 Table 3-2: Estimated Connected Three Phase Boundary Length using Connectivity Coefficient. ..... 56 Table 3-3: Estimated Connected Three Phase Boundary Length using Connectivity Coefficient. ..... 59 Table 4-1: Electrolyte, anode and modeled cathode properties used in development of the method for calculating the effective diffusion coefficient for modeled sintered composite electrodes. ..................................................................................................................................... 62 Table 4-2: Verification of total porosity for three different electrode geometries discretized using 100 cross section images, each image having 10000 pixels. ................................................. 64 Table 4-3: Summary of effective diffusion coefficient calculated using three different mesh densities based on three different image resolutions as shown in Figure 4-11............................... 80 Table 5-1: Summary of properties defining five examples of modeled electrodes. Also summarized are porosity, three phase boundary length, and the corresponding connectivity. ............. 85 Table 5-2: Comparison of connected three phase boundary length calculated in this work to connected three phase boundary length using combined model based on available literature. ......... 93 Table 5-3: Comparison of effective diffusion coefficient for three example electrodes. The bulk binary diffusion coefficient is 2.01x10-4m2/s for the O2-N2 system at 1123K and 1atm. Fuel cell operating point is 0.7V. .................................................................................. 95    xii Nomenclature a  Face centered cubic lattice constant m projA  Projected area of electrode m2 zporeA ,  Cross section area of pore space at position z m2 ic  Concentration of species i mol m-3 o ic  Reference concentration of species i mol m-3 max,,2 zaveO c  Average oxygen concentration at maximum z coordinate mol m-3 d  Diameter of solid particle m 0d  Reference particle diameter on face centered cubic lattice m meand  Mean particle diameter on face centered cubic lattice m d* Approximate molecular diameter m mind  Minimum particle diameter on face centered cubic lattice m maxd  Maximum particle diameter on face centered cubic lattice m eff jiD −  Effective binary diffusion coefficient for species i and j m2 s-1 jiD −  Bulk binary diffusion coefficient for species i and j m2 s-1 D Distance between sinter neck origin and point on interfering sinter neck perimeter m F Faraday’s constant C mol-1 if  Solid volume fraction of species i h Height of spherical cap m Fi  Faradaic current density A m-2 1,QXYI  Sinter neck located in XY plane, quadrant one zFi ,  Faradaic current density at position z A m-2 zFsi ,  Scaled Faradaic current density at position z A m-2 effk  Effective electrical conductivity S m-1 0 ik  Bulk electrical conductivity of species i S m-1   xiii ijL  Distance between centers of particles i and j.  m tpbl  Contact length between two particles of opposite phase m Mq Molecular weight of species q g mol-1 Mq Three dimensional matrix containing locations of phase q pN  Number of lattice sites on face centered cubic lattice fccn Number of face centered cubic unit cells per side of modeled geometry 2O N  Molar oxygen flux mol m-2 s-1 NA Avogadro’s constant νn  Number density of gas molecules [number] m-3 n Number of particles within electrode qn  Number of q type particles in electrode nel Number of moles of electrons per mole of consumed species voxcN ,  Number of connected three phase boundary voxels voxtN ,  Total number of three phase boundary voxels zvoxcN ,,  Number of connected three phase boundary voxels within element layer z zfacen , Number of boundary faces having centroid within element layer z zavgON ,,2 Average superficial oxygen diffusive flux in xy plane at position z mol m -2  s-1 avgON ,2  Arithmetic mean superficial oxygen diffusive flux mol m -2  s-1 weightedON ,2  Weighted mean superficial oxygen diffusive flux mol m -2  s-1 P Pressure Pa pi Partial pressure of species i Pa Pi Fraction of i type particles part of a percolated cluster P(s,t) Point on parametric circle with parameters s and t xQ  Generic source term for species x [species] m-3 s-1 sr  Sinter neck radius m   xiv nr  Radius of sinter neck curvature m pR  Polarization resistance Ω m2 R Gas constant J mol-1 K-1 S Surface area available for electrochemical reaction m2 m-3 Si Area of boundary surface i m2 ijS Distance from center of particle i to plane of intersection with particle j m T Temperature K sV  Solid volume of particles in modeled electrode geometry m3 EV  Envelope volume of modeled geometry m3 nV  Volume of sinter neck  m3 cV  Volume of spherical cap overlapping sinter neck m3 v Position vector vi Unit vector within plane of intersection between two particles. Subscript i ranges from 1 to 4 il ,ν  Length of voxel in i direction. i is either x, y or z. Vq Fuller diffusion volume of species q m3 tpbw  Width of three phase boundary line m xi x coordinate of center of particle i m Dx  Ratio of open porosity to tortuosity factor yi y coordinate of center of particle i m zi  z coordinate of center of particle i m Z  Average coordination number qpZ − Average number of p type particle contacts on particle type q qZ  Average number of contacts on particle type q pα  Ratio of particle radii pqα  Charge transfer coefficient, electrode p, direction q Tε  Total porosity   xv ε  Open porosity mfpλ  Mean free path length m cλ  Connected three phase boundary length m m-3 totalλ  Total three phase boundary length m m-3 η  local overpotential V φ  Electric potential V ijφ  Angle between position vectors vi and vj radian 2 dσ  Variance in particle diameter on face centered cubic lattice m2 0θ  Reference overlap angle between contacting particles radian θ  Particle contact angle radian iτ  Tortuosity of solid phase i τ  Tortuosity factor     xvi Acknowledgements I am grateful to all of my colleagues, friends, and family, who have supported me during this project. I would like to thank Justin Gammage, whose encouragement was pivotal to my decision to pursue graduate studies and who provided the occasional boating escape fueled by C2H5OH. I also want to thank Kylee Metcalfe, my cousin and house-mate in Vancouver, whose friendship and support over the past two years I could always count on. I am most thankful to my research supervisor, Dr. Olivera Kesler, for her guidance, trust and support, both professionally and personally, without which this work would not have been possible. I would also like to thank Tony Rivard, François Gitzhofer, and Nicolas Abatzoglou for collaborating on this project. I would like to thank my parents, brother and family-in-law for their patience, support and understanding, and of course for helping me to move (twice) between Toronto and Vancouver. Lastly, I do not know the words to fully express on paper my thanks to my wife, Jennifer, for her love and support, even over the long distance that was between us. Thank you.    xvii                        This work is dedicated to the continued study of clean energy.   1 Chapter 1 Introduction 1.1 Solid Oxide Fuel Cells 1.1.1 Role in Society The increasing availability of electricity over the past century has brought about significant increases in living standards and facilitated unprecedented economic development. The demand for electricity and the global population will continue to rise into the future. By 2050 it is expected that the global population will reach 12 billion and the demand for global energy services will increase by as much as an order of magnitude [1]. Conventional power generation methods are generally based upon vapor power (Rankine) or gas turbine (Brayton) thermodynamic cycles. The maximum efficiency is therefore limited to that of the reversible Carnot cycle. Combustion processes used to generate heat required for these cycles produce emissions including oxides of nitrogen, sulfur dioxide, and carbon dioxide, which contribute to smog, acid rain, and global warming. It is unlikely that the demand for electricity will subside as the global population rises and living standards increase. It is more likely that the demand for electricity will grow, and therefore, sustainable methods to produce electricity must be considered. Energy conversion devices capable of producing electricity in a clean and efficient manner within the earth’s limits will be required as society evolves. Part of the increasing energy demand could potentially be met with solid oxide fuel cells (SOFCs), which are an energy conversion device not limited by the Carnot cycle, as they convert chemical energy directly into electrical energy. Therefore, energy conversion in fuel cells is generally more efficient compared to energy conversion in combustion-based systems. Practical applications of solid oxide fuel cells include stationary power markets where the high quality by-product heat can be used in combined heat and power facilities. The ability of SOFCs to utilize hydrocarbon fuels allows use of existing fuel delivery infrastructure, enabling a transition from non-renewable energy supplies to renewable energy supplies. Such   2 systems have been developed for demonstration purposes, consisting of an SOFC generator and micro-gas turbine [2], and show promise for future development of solid oxide fuel cell technology. 1.1.2 Operating Principle The single solid oxide fuel cell consists of two porous electrodes separated by a dense electrolyte. Figure 1-1 depicts a single solid oxide fuel cell. Fuel is oxidized at the anode, releasing electrons to an external circuit. Oxygen in air is reduced at the cathode, where electrons are accepted from the external circuit. Oxygen ions, produced during the oxygen reduction reaction, are conducted through the electrolyte to complete the circuit.   In contrast to polymer electrolyte fuel cells currently under extensive development, solid oxide fuel cells are not limited to highly pure hydrogen or methanol as the only fuels. In fact, any combustible compound may theoretically be oxidized at the anode of a solid oxide fuel cell, provided that the anode materials catalyze the reactions, which can include Figure 1-1: Exploded view of a single planar solid oxide fuel cell showing anode, electrolyte, cathode and interconnect. Multiple repeating units are combined to form a stack. Repeating unit Interconnect (current collector) Cathode Electrolyte Anode Fuel Air   3 electrochemical oxidation, internal steam reforming, and water-gas shift reactions. Since most fuels contain at least some hydrogen, the electrochemical oxidation of hydrogen is almost always present, and hydrogen is the simplest fuel that can be used in the fuel cell. Therefore, this work focuses on the oxidation of hydrogen at the anode and reduction of oxygen at the cathode according to (1) and (2), respectively. −− +→+ eOHOH 22 2 2  (1) −− →+ 22 22 1 OeO  (2) 1.1.3 Outlook A significant challenge to the solid oxide fuel cell research community is to increase cell performance in the intermediate operating temperature range of 600℃ to 800℃. In this range, it is expected that cell durability can be increased, due in part to a reduction in thermal stresses compared to operation at the more traditional high temperature range of 800℃ to 1000℃. Operation below approximately 800℃ also mitigates oxidation of stainless steel; therefore, cell material costs can be reduced if this material can be used for the interconnect rather than the more expensive lanthanum chromites that are traditionally used. The primary disadvantage of reducing cell operating temperature is the corresponding reduction in output current resulting from slower electrode reaction kinetics and also from a decreased ionic conductivity - a material property with an exponential dependence on temperature. For example, a common electrolyte material is yttria stabilized zirconia, which has a conductivity approximately two orders of magnitude higher at 1000℃ compared to its conductivity at 600℃  [3]. Therefore, either new materials are required for lower temperature operation, or the cell design must be optimized. This work supports the latter objective, to provide greater insight into the relationship between electrode structure, composition, and performance.   4 1.2 Solid Oxide Fuel Cells with Functionally Graded Electrodes A common electrode configuration for solid oxide fuel cells is a porous composite of two solid phases; one phase is both an electronic conductor and catalyst, while the second phase is an oxygen ion conductor. In common cathodes and anodes, the electrochemical reactions may occur only where three phases meet: the electron conductor, ion conductor, and reactant within a pore, as illustrated in Figure 1-2. All locations where these three phases meet are termed three phase boundaries, and in order for the electrochemical half cell reaction to occur, each phase must be connected from the reaction site to the source of the transported species. Specifically, there must be an electron conduction path to the current collector (interconnect), an oxygen ion conduction path to the electrolyte, and a gas diffusion path to the air channel in a cathode or fuel channel in an anode. The connected three phase boundaries are thus a subset of the total three phase boundaries within a composite electrode. Extending the three phase boundaries to the bulk of the electrode is the primary benefit of composite electrodes.   5  A further distinction can be made by considering active three phase boundaries as a subset of connected three phase boundaries. The mobility of oxygen ions is lower than that of electrons and is directly proportional to electrical conductivity. For example, a typical composite cathode consists of strontium doped lanthanum manganite (LSM), an electronic conductor, and yttria stabilized zirconia (YSZ), an oxygen ion conductor, which have conductivities of 2.66·104 S/m and 4 S/m at 1123K, respectively [4]. Therefore, the distribution of ionic current is not uniform throughout the electrode, but rather, has increasing magnitude toward the plane of intersection with the electrolyte [5]. Therefore, connected three phase boundaries nearer to the plane of intersection with the electrolyte are generally more active than connected three phase boundaries nearer to the plane of intersection with the current collector, due to the lower mobility of oxygen ions compared to that for electrons. The uneven distribution of ionic current within composite electrodes suggests that near the cathode and anode interface with the electrolyte there must be high catalytic activity, and a Figure 1-2: Schematic of ion conducting particles (white) and electron conducting particles (gray) in a composite cathode near the electrolyte surface. Electrons from the external circuit are conducted through the electron conducting particles to the reaction site, where oxygen ions result from the oxygen reduction reaction.   6 high three phase boundary length, to facilitate the reduction of oxygen and oxidation of fuel, respectively. Near the cathode and anode interface with the current collector it is desired to have a material with high electronic conductivity to minimize ohmic losses within the electrode. Furthermore, a microstructure having larger pore and particle sizes near the electrode-current collector interface compared to near the electrode-electrolyte interface would be expected to have at least two benefits compared to an electrode having constant, uniformly distributed, particle and pore sizes. First, the larger pores will improve the diffusive transport of gas species into the bulk of the electrode. Second, the smaller pore and particle sizes near the electrode-electrolyte interface will increase the three phase boundary length in this region of higher catalytic activity. The contrasting composition and microstructure requirements at opposite electrode interfaces suggest that introducing stepwise or continuous composition and microstructure gradients would improve electrode performance. An idealized grading scheme has been proposed by Liu et al. [6], in which a tri-layered cathode was fabricated with the two layers nearest the electrolyte having small particle sizes and pore sizes, while the third layer had larger particle sizes and pore sizes. The layer nearest the electrolyte had a composition of 40% electron conducting lanthanum strontium manganite (LSM) and 60% ion conducting gadolinia doped ceria (GDC), with the LSM providing the additional benefit of a closer thermal expansion coefficient (TEC) match to the yttria stabilized zirconia (YSZ) electrolyte compared to the TEC of the material added in subsequent layers. This is in addition to the LSM being active for the oxygen reduction reaction. The next layer was a composite of 30% LSM, 30% mixed electron and ion conducting lanthanum strontium cobaltite (LSC), and 40% GDC to increase the electronic conductivity of this layer. The final layer further increased the electronic conductivity with a composition of 60% LSC and 40% GDC, while having larger particle and pore sizes to facilitate oxygen gas transport to the first two layers, which have a higher electrochemical activity. There have been similar studies on cathodes with functionally graded composition using LSM and YSZ [7], LSM, LSC, and YSZ [8-12] and by Xu et al.[13] using LSM with samaria doped ceria (SDC) as the ion conducting phase. In each case the polarization resistance of   7 cells with graded cathodes was lower than that of cells with single layer, uniform composite cathodes, in the intermediate and high temperature range. The kinetics of the oxygen reduction reaction are generally slower than those of the hydrogen oxidation reaction at the anode, and therefore, the use of anodes with graded properties has largely focused on microstructural gradients to improve mass transport. Anodes in planar anode-supported solid oxide fuel cells are generally fabricated with two layers. The first is a support layer on the order of 1mm thick and the second is a functional layer on the order of tens of micrometers thick. The support layer provides structure to the cell; however, the relatively large thickness may reduce the effective diffusivity. Therefore, porosity gradients have been introduced [14- 16], characterized by larger pore sizes at the gas channel to smaller pore sizes toward the functional layer. The experimental work described above shows promise toward fabricating solid oxide fuel cells with improved performance; however, the work has not explored a sufficiently wide range of electrode design parameters to fully characterize the effect of graded properties on cell performance. Furthermore, from a modeling perspective, the geometry of a functionally graded electrode is complex and difficult to describe mathematically. Due to the large computational resources required to directly model charge, mass, and heat transport though real solid oxide fuel cell geometries, macroscopic models utilizing effective coefficients have generally been used. However, the addition of graded properties complicates this approach, as the porous medium can no longer be considered uniform. This work seeks to define a method by which the geometry and composition dependent variables, such as three phase boundary length and effective diffusivity, may be evaluated over a sufficiently large design space. This goal is intended to allow macroscopic performance models to be applied to functionally graded electrodes consisting of discrete layers, each layer having uniform, homogeneous properties. 1.3 Electrode Performance, Composition and Structure In this work, electrode performance is evaluated by determining the polarization resistance of the electrode. At a given operating point, the polarization resistance is the slope of a graph depicting cell potential on the ordinate and Faradaic current on the abscissa. The Faradaic   8 current depends in large part on the composition of materials used and the electrode microstructure. The following review outlines relevant literature in an effort to better understand the relationship between composition, microstructure, and electrode performance. 1.3.1 Charge and Mass Transport in Solid Oxide Fuel Cell Electrodes The steady state transport of charge and mass within solid oxide fuel cell electrodes may be described using the general transport equation, for species x in equation (3), where the divergence of the flux of species x, xN , is equal to the rate of production or consumption of species x, xQ . xx QN =⋅∇  (3) Applying porous electrode theory, solid oxide fuel cell electrodes may be modeled as continuum materials characterized by effective properties that are defined at every point, as described by Newman and Tobias [5]. Equation (3) may be written for both charge and mass transport as given in equations (4) and (5), respectively, based upon work by Gazzarri and Kesler [17] as part of an alternating current impedance model for identifying degradation modes in solid oxide fuel cells. For equation (4), effk  is the effective electrical conductivity, φ  is the electric potential, S  is the surface area available for electrochemical reaction, and ( )ηFi  is the local Faradaic current density, which depends on the local overpotential, η . Equation (4) may be applied to describe transport of either electrons or oxygen ions within the cathode and anode. For equation (5), effjiD −  is the effective binary diffusion coefficient for gas species i and j, R is the universal gas constant, T is the temperature in Kelvin, F is Faraday’s constant, ip  is the partial pressure of gas species i, and eln  is the number of moles of electrons transferred per mole of gas species i. ( )ηφ Feff iSk =∇− 2  (4) ( ) Fn iSp RT D el F i eff ji η =∇− − 2 (5)   9 The effective conductivity in equation (4) is obtained through application of a correction to the conductivity of the bulk material. The specific correction to apply is not generally agreed upon in the literature. For example, Deseure et al. [18] correct the bulk conductivity according to equation (6), where Tε  is the total porosity, if is the solid volume fraction of solid phase i, 0ik  is the bulk electrical conductivity of solid phase i, and iτ is the tortuosity of solid phase i. The tortuosity depends upon the morphology of the electrode and is generally described as the ratio of the average transport path length to the length of the porous medium along the major transport direction. Tortuosity values of solid phases have not been reported in the literature, to the knowledge of the author, and are generally approximated by the gas phase tortuosity. The gas phase tortuosity also depends upon electrode morphology and has experimentally measured values for porous sintered ceramics that typically range from two to six [19, 20], while the value of three is suggested for anodes with porosities greater than 30% [20]. In addition to there being less certainty regarding the gas phase tortuosity for porosities below approximately 30%, the assumption of substituting the solid phase tortuosity for the gas phase tortuosity becomes less plausible as the solid volume and pore volume fractions diverge. 01 i i T i eff kfk τ ε− = (6) A second example of a correction for bulk conductivity is used by Kenney et al. [21], according to equation (7), where iP  is the fraction of i phase particles belonging to a percolated cluster. The solid phase tortuosity is not considered, but rather, the fraction of i phase particles belonging to a percolated cluster is approximated by the solid volume fraction of phase i, for volume fractions below the upper percolation threshold and above the lower percolation threshold [22]. Above the upper percolation threshold, the fraction of i phase particles belonging to a percolated cluster is one. The form of the conductivity correction remains an open problem. ( ) 01 iiiTeff kPfk ε−=  (7)   10 In equation (5), Fick’s law is used to describe the molar flux of species i in a binary system of two gases i and j. For gas systems with more than two components, for example, anodes with methane as the fuel, Fick’s law is replaced by the Maxwell-Stefan [23] equations or the Dusty Gas model [24]. For this work, the binary gas systems H2-H2O at the anode and O2-N2 at the cathode are considered, and Fick’s law is used to describe the diffusion of hydrogen within the anode and oxygen within the cathode. The correction to the bulk binary diffusion coefficient, jiD − , is the ratio of electrode open porosity, ε , to gas phase tortuosity factor, τ , as shown in equation (8). ji eff ji DD −− = τ ε  (8) The porous medium is generally modeled based on a bundle of sinuous, but parallel, capillaries or pores. The tortuosity factor is the square of the tortuosity; which is defined as the ratio of the average pore length to the length of the porous medium along the major diffusion axis [25]. The tortuosity factor accounts for both the additional path length and the increase in molar diffusive flux through the capillaries compared to the interstitial axial molar diffusive flux. As noted in [25], both effects arise from the deviation, on average, of the pore direction compared to the major diffusion direction. With regard to correcting the binary diffusivity in porous solid oxide fuel cell electrodes, this definition does not consider the chemical and/or electrochemical reactions occurring within the bulk of the porous electrode. Therefore, the physical meaning of the tortuosity factor is not well defined in the context of fuel cell electrodes, and the quantity might better be considered as a general diffusivity correction factor. 1.3.2 Electrode Performance In both the charge and mass transport equations, (4) and (5) respectively, the source term contains the Faradaic current density. This quantity is commonly modeled in the solid oxide fuel cell literature using the Bulter-Volmer equation, (9), written here for the anode and modified to include the concentration of the species that is oxidized in the reaction, 2H c , with respect to a reference concentration, 20 Hc , as well as the species that is reduced in the reverse   11 reaction, OHc 2 , with respect to a reference concentration, OHc 2 0 . The Faradaic current density describes the rate of charge transfer, ( )ηFi , as a function of the overpotential, η , which is defined as the difference in potential between the electron conducting phase and ion conducting phase with respect to the potential of a reference electrode. In equation (9), aaα and acα  are kinetic parameters which describe the relative proportion of the electrochemical half cell reaction, (1), in the anodic and cathodic directions, respectively. The exchange current density, oi , is a third kinetic parameter describing the current density corresponding to the rate of reaction in the anodic and cathodic directions of equation (1) at equilibrium. ( )             −−      = ηαηαη RT F c c RT F c c ii ac OH OHaa H H oF expexp 2 2 2 2 00  (9) The overpotential, η, represents a departure from the equilibrium potential of the solid oxide fuel cell. The magnitude of the overpotential is affected by three primary mechanisms: ohmic resistance of the electrode materials, kinetics of the electrochemical reaction, and mass transfer within the electrode. Using an electrical analogy, a series resistance describes the ohmic losses and the concept of a polarization resistance, PR , is used to describe the kinetic and mass transport losses within the electrode. The polarization resistance may be described mathematically as the reciprocal of the derivative of the Faradaic current density with respect to the overpotential as in equation (10). A larger magnitude of PR indicates that there are higher kinetic and mass transport losses occurring within the electrode. For performance optimization of solid oxide fuel cell electrodes, a more complete understanding of the effect of electrode composition and microstructure on polarization resistance is required. ( ) η η ∂ ∂ = − F p iR 1  (10) The polarization resistance has been shown to vary inversely with the length of the three phase boundary line on patterned electrodes involving an interface between a purely electronic conductor (electrode) and purely ionic conductor (electrolyte) [26]. However, reducing the particle size to increase the number of three phase boundaries for a given volume may not necessarily decrease the polarization resistance. This is because smaller   12 particle sizes tend to lead to lower effective diffusion coefficients. The electrochemically active surface area (product of three phase boundary length and three phase boundary width) multiplies the Faradaic current density in equations (4) and (5). The width of the three phase boundary line depends on the specific reaction mechanism limiting the rate of the half-cell reaction, which itself is dependent on the cell operating point. For this work, a value of 0.05µm is used for the width of the three phase boundary line for cathodes comprised of lanthanum strontium manganite and yttria stabilized zirconia [27]. This three phase boundary width is specific to the electrode materials used; however, as a first approximation in the absence of data available in the literature, a width of 0.05µm is assumed to also apply to anodes consisting of nickel and yttria stabilized zirconia. Both variables, S and ( )ηFi , depend on electrode composition and microstructure and therefore a greater understanding of their effect on electrode performance is required to model electrodes with graded properties. 1.3.3 Three Phase Boundary Length The models in the literature to calculate the three phase boundary length in composite electrodes are largely based on the random packing of spherical particles that overlap nearest neighbor particles. The classic models developed by Sunde [28] and later by Costamagna et al. [29] model electronic and ionic conducting particles as a bimodal distribution of spheres contacting at an angle θ, as shown in Figure 1-3. In the work by Costamagna et al., the authors do not calculate the three phase boundary length, but rather, they make the assumption of a bulk diffusion limited reaction mechanism and take the intersection area between the contacting particles to be the electrochemically active surface. However, the method used to estimate the fraction of connected area and fraction of connected length are common to both models.   13  Using the model by Sunde as an example, the connected three phase boundary length, cλ , is calculated according to equation (11), written for the case where the electronic and ionic conducting particles are of equal diameter for clarity of description. tpbl  is the contact circumference between two particles of opposite phase, n  is the total number of particles in the electrode, and ef  and if  are the solid volume fractions of the electronic and ionic conducting phases, respectively. Z  is the average coordination number of the electrode and eP  and iP  are the fraction of electronic conducting particles and ionic conducting particles that are part of a connected cluster of the respective solid phases. ( )( ) ieietpbc PPfZnfl=λ  (11) Equation (11) reads that the connected three phase boundary length is equal to the product of the length of intersection per particle contact of opposite phase, tpbl , the number of electronic conducting particles within the electrode, ( )enf , the number of ionic conducting Figure 1-3: Contact between electron and ion conducting particles. Schematic is redrawn based on [29]. θ Ion conducting particle Electron conducting particle   14 particles in contact with an electronic conducting particle, ( )ifZ , and the fractions of both electronic and ionic conducting particles that form electrically connected clusters, iePP . The coordination numbers as outlined by Bouvard and Lange [30] play a central role to the development of equation (11) and are summarized in equations (12) through (15), where the subscript e refers to the electron conducting phase and the subscript i refers to the ion conducting phase. Therefore, eeZ −  is the average number of electron conducting particles in contact with an electron conducting particle, ieZ −  is the average number of ion conducting particles in contact with an electron conducting particle, eiZ −  is the average number of electron conducting particles in contact with an ion conducting particle and iiZ −  is the average number of ion conducting particles in contact with an ion conducting particle. ieeee ZZZ −− +=  (12) iieii ZZZ −− +=  (13) iiee ZfZfZ +=  (14) 1=+ ie ff  (15) eZ , therefore, represents the average number of particles contacting an arbitrary electronic conducting particle and iZ  represents the average number of particles contacting an arbitrary ionic conducting particle. Z  represents the overall average number of contacts on a particle. In order to express the unknown coordination numbers in equations (12) and (13) as a function of the overall average coordination number, the solid volume fraction, and the ratio of particle radii, the following assumption is made as shown in equation (16). This assumption reads that the fraction of ion conducting particles in contact with an electron conducting particle may be represented by either the left hand side or the right hand side of equation (16). Z Zn Z Z ii e ie = −  (16)   15 Z Zn Z Zn iiee +=1  (17) e ie e ee Z Z Z Z −− +=1  (18) If equations (12) and (14) are rearranged as shown in equations (17) and (18), it can be seen that the assumption holds if the following is true. Z Zn Z Z Z Zn Z Z ii e e-iee e ee == − and  (19) Bouvard and Lange found that the fraction of particles in connected clusters is a function of eeZ −  or iiZ − for the electronic and ionic conducting particles, respectively. For example, the fraction of electronic conducting particles that are part of a percolated cluster is given by equation (20), for volume fractions of particles between the upper and lower percolation thresholds. eeZ −  is given by equation (21) and is analogous to iiZ − where the subscript e is replaced with the subscript i [31]. αp is the ratio of particle diameters of the two phases. 4.05.2 2 4 1               − −= −ee e Z P  (20) ( )[ ]21 pee e ee nn nZ Z α−+ = −  (21) At the lower percolation threshold eeZ − is two. Below this threshold there is no connectivity of the electronic conducting particles, whereas, for eeZ −  greater than four, all electronic conducting particles belong to a connected network. Kuo and Gupta, [32], found that equation (20) better fitted experimental data by setting eeZ − equal to 1.764 at the lower percolation threshold, for ratios of ionic to electronic conducting particle radii between 0.154 and 6.464, the region for which the composite powders do not segregate. Setting eeZ −  to 1.764 at the lower percolation threshold, as opposed to setting eeZ − equal to two at the lower percolation threshold, appears to be more common in the literature [31, 33]. However, it   16 should be noted that equation (20) was developed using numerically generated random packings of spheres with a bimodal size distribution, in which the packing density was between 0.57 and 0.61. Furthermore, Kuo and Gupta only considered numerically generated random packings in which the packing density was equal to 0.637. Therefore, the application of equation (20) to solid oxide fuel cell electrodes with porosity outside the approximate range 0.36<ε <0.43, requires that the variation of coordination number with porosity be known. There have been attempts to describe the average coordination number in random spherical packings for monodisperse particles [34, 35]. However, these models, and those above, all contain the implicit assumption that the porosity of the electrode is completely open, allowing gas to diffuse to all reaction sites. Above the upper percolation threshold this assumption may be valid; however, without consideration of the open porosity, the predicted three phase boundary length continues to increase as porosity decreases, as illustrated in Figure 1-4. Here, the average coordination number varies according to equation (22), for porosity less than 0.82 [34]. ( ) 82.061.1 48.1 ≤= − εεZ  (22) Physically, the trend in Figure 1-4 indicates that the three phase boundary length increases as porosity decreases. However, for packing densities above the maximum theoretical value of 0.74, corresponding to porosities below 0.26, the coordination number is 12 for spherical particles and cannot increase above this value. Therefore, a linear trend in three phase boundary length with porosity occurs. This implies that all porosity is open, allowing gas diffusion to all reaction sites; however, the open porosity also has a percolation threshold, below which percolated pore networks do not exist, which contrasts the trend seen in Figure 1-4.   17  A recent Monte Carlo percolation model has been developed by Martinez and Brouwer [36] with consideration to percolation of the pore space. A two dimensional randomized structure was developed for which the electrode is discretized on a square grid, each square being assigned one of three phases, and the connectivity of each phase is tracked to directly determine the number of connected three phase boundaries. The electrode geometry was not modeled by Martinez and Brouwer; however, their work is unique in that it attempts to model the variation in connected three phase boundaries with consideration to the connectivity of the pore phase. A three dimensional reconstruction of a solid oxide fuel cell anode was accomplished by Wilson et al. [37] using focused ion beam milling of a real anode and imaging each section with a scanning electron microscope. This method represents a significant improvement to 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 1014 Porosity l tp b (m /m 3 ) Figure 1-4: Connected three phase boundary length predicted using equation (11), where the average coordination number is a function of porosity according to work by Nakagaki and Sunada [34]. The result predicts increasing three phase boundary length as porosity decreases to zero. Th re e Ph a se  Bo u n da ry  L en gt h (m /m 3 )   18 the modeling of solid oxide fuel cell geometry. With a digital reconstruction of the geometry, the three phase boundary line may be directly modeled, subject to the difficulties of discerning each phase from the grayscale SEM images, and the connectivity of each identified three dimensional pixel (voxel) may be evaluated. Transport modeling of electrons, oxygen ions, and gaseous diffusing species becomes conceivable using the reconstructed geometric domain. To date, this technique has not been applied to a wide range of electrode structures, and therefore, a characterization of connected three phase boundaries over a range of electrode compositions and microstructures has not been realized. Furthermore, the resolution obtained using the focused ion beam milling technique, corresponding to the spacing between successive milling operations, is currently limited to approximately 25nm. It is not known if this resolution is sufficient to reconstruct the three phase boundary length such that further increases in resolution will yield the same three phase boundary length. The use of deterministic electrode models, rather than real electrodes, in which the model may be sectioned with a resolution less than 25nm, may prove useful in at least three ways. First, it may be possible to elucidate the required FIB milling resolution to ensure three phase boundary length convergence. Second, the use of deterministic geometry models allows verification of the algorithms used to evaluate connectivity between the two solid phases, and pore phase, in real electrodes for which the connectivity of each phase is unknown. Finally, it allows the evaluation of effective properties for a broad range of electrode microstructures, thus enabling rational design of graded electrodes at a lower experimental cost. 1.3.4 Diffusion of Gases in Porous Electrodes Diffusion of gases through porous electrodes generally falls into one, or a combination of, the following three categories: ordinary, Knudsen, and surface diffusion. Ordinary diffusion is characterized by the diffusing molecules colliding most frequently with other diffusing molecules, whereas Knudsen diffusion is characterized by diffusing molecules colliding most frequently with the pore walls. Surface diffusion involves adsorbed gas species diffusing on the surface of the pore wall; however, the surface diffusion mechanism is not considered in this work, as it is not expected to significantly contribute to the overall gas transport [38].   19 The mean free path of a diffusing molecule, mfpλ , is estimated using equation (23), which approximates the diffusing molecule as a hard sphere with diameter d* [39]. nv is the number of gas molecules per unit volume. If the ideal gas law is used to estimate the number of gas molecules per unit volume, the mean free path is given by equation (24), where NA is Avogadro’s number and P is the pressure. v mfp nd *2 1 pi λ =  (23) PNd RT A mfp *2pi λ =  (24) At an operating temperature of 1123K and pressure of 1atm, the mean free path of oxygen and nitrogen, assuming both molecular diameters are approximately 0.3nm, is 0.4µm. This agrees with the mean free path for these molecules used by Kenney and Karan [40]. By comparison, hydrogen, at the same temperature and pressure, has a mean free path of 0.9µm if a molecular diameter of 0.2nm is assumed. Therefore, pore sizes less than 1µm are likely to require consideration of Knudsen diffusion. The Knudsen number is useful in determining the pore size for which the Knudsen diffusion mechanism becomes significant and is defined as the ratio of the mean free path to a representative length scale. In a porous electrode, the representative length scale would be some measure of the average pore size. When the Knudsen number is much greater than one, Knudsen diffusion dominates, whereas, when the Knudsen number is much less than one, ordinary diffusion dominates. When the Knudsen number is approximately equal to one, then both diffusion mechanisms should be considered.  It is common in the solid oxide fuel cell literature to consider Knudesn diffusion by defining an average pore diameter by assuming cylindrical pores [41-44]. The modeled geometry used in this work, described in Chapter Two, approximates electrode structures that are typical of tape cast or screen printed electrodes. The porosity consists of missing particles, representative of using pore forming material during the fabrication process, in addition to the void space between the solid particles. Therefore, it becomes less clear what should be the average pore size when calculating the Knudsen diffusion coefficient. As a first step in development of a model to describe effective diffusivity in solid oxide fuel cell   20 electrodes, consideration of Fickian diffusion is sufficient to be applicable to a significant portion of electrode microstructures reported in the literature, where the pore sizes are larger than the mean free path of the diffusing molecules. 1.4 Scope of Thesis This work seeks to improve understanding of the relationship between electrode performance, composition and microstructure through the development of a method to determine the connected three phase boundary length and effective diffusion coefficient in deterministic geometries. The effective electrical conductivity is left for later study. Three dimensional solid models of sintered composite electrodes with varied composition, mean particle size, and total porosity form the basis for evaluating the corresponding change in connected three phase boundary length and effective diffusion coefficient. Tomographic methods developed for three dimensional reconstruction of real electrodes are applied to well defined, deterministic electrode geometries. The specific method of evaluating the connectivity of each phase in the reconstructed model, namely, the ion conductor, electron conductor, and pore, is obtained through collaboration with researchers Tony Rivard, François Gitzhofer, and Nicolas Abatzoglou at the University of Sherbrooke. Chapter Two outlines the method used to create the three dimensional solid model of sintered composite electrodes for a desired material composition, total porosity, and mean particle diameter. The algorithm developed to calculate the total geometric three phase boundary length for the modeled geometry is outlined in Chapter Three. In addition, Chapter Three describes the tomographic method used to reconstruct the three dimensional geometry and subsequently, to evaluate the connectivity of each phase: ion conductor, electron conductor, and pore. Chapter Four describes the algorithm developed to evaluate the effective diffusion coefficient of the modeled geometry. The results for connected three phase boundary length and effective diffusion coefficient for selected cases of the electrode design space defined in Chapter Two are discussed in Chapter Five. Chapter Six summarizes the conclusions derived from this work and the recommendations and near future work planned based on that presented herein.   21 Chapter 2 Geometric Modeling of Sintered Composite Electrodes 2.1 Introduction The electrode geometry model described in what follows is intended to be representative of sintered, composite solid oxide fuel cell electrodes fabricated using traditional wet ceramic techniques, namely, tape casting and screen printing. The resulting geometry model is not used directly, but rather, cross section images are created, as described in Chapter Three, for use in digitally reconstructing the geometry in a form amenable to analysis of the connectivity of each phase. 2.2 Definition of Design Space The electrode design space summarized in Table 2-1 describes the domain for the three independent variables chosen to characterize the electrode geometry. These are the total porosity, electronic conductor solid volume fraction and the mean particle diameter. Independent Variable Minimum Value Maximum Value Total Porosity 0.05 0.55 Electronic conductor solid volume fraction 0.1 0.9 Mean particle diameter 1.0 µm 10 µm At low porosities, there are few pore paths connecting the fuel channel, for anodes, or air channel, for cathodes, to the reaction sites. As the total porosity increases, a larger fraction of the pores becomes connected and fuel, or air, may diffuse to a greater number of reaction sites. The porosity values chosen are expected to span the range where there are very few, if any, connected pore paths near 0.05 to a case where all pore paths are fully connected near 0.55. The electronic conductor solid volume fraction is arbitrarily selected as an independent variable. Since there are only two solid phases comprising the composite electrode layers Table 2-1: Composite Electrode Design Space   22 considered in this work, the ionic conductor solid volume fraction is the complement of the electronic conductor solid volume fraction. Both solid phases must be present in sufficient quantity to form an electrical connection between a given electronic conducting particle and the current collector, and likewise, between a given ionic conducting particle and the electrolyte. The lower site percolation threshold for a face centered cubic lattice is 0.12 [45], indicating the fraction of lattice sites that must be occupied to create the first continuous path from the current collector to the electrolyte surface. Therefore, a value of 0.1 electronic conductor solid volume fraction, corresponding to a site fraction of 0.1 if there are no missing particles from the lattice, is just below the lower site percolation threshold of 0.12. Likewise, the upper limit of 0.9 electronic conductor solid volume fraction is just below the lower site percolation threshold for the ionic conductor. The face centered cubic lattice is a cubic close packed structure and is used as a first approximation to modeling the packing of spherical particles in a solid oxide fuel cell electrode. This is based on the observation that hard spheres placed within a container will generally assume a close packed structure. The specific structure of face centered cubic close packing was chosen over hexagonal close packing, and over randomly stacked close packed planes for computational simplicity, while the randomness of particle and pore placements was introduced numerically. The particle diameter for both solid phases was chosen to vary from 1µm to 10 µm, based on observations of commonly-observed feature sizes in electrode micrographs from the solid oxide fuel cell literature. The three phase boundary length is calculated on a volumetric basis [m/m3], and as such a factor of two decrease in particle diameter results in a factor of four increase in three phase boundary length. This relationship suggests that selecting particles with as small a diameter as physically possible is desirable; however, there is a lower limit. For pore sizes on the order of the mean free path of the diffusing species, ranging from approximately 0.4µm to approximately 0.9µm, as discussed in Chapter One, the diffusing molecules begin to collide more frequently with the pore walls than with other diffusing molecules. As described in Section 2.6, porosity is introduced by removing particles from the face centered cubic lattice, and therefore, pores introduced by removing particles from the lattice cannot be smaller than 1µm.  However, the interstitial pores may be smaller than 1µm, and therefore, there is uncertainty in the diffusion coefficient for these cases if only ordinary   23 diffusion is considered. Therefore, the effective diffusion coefficients calculated for this thesis will be for electrodes having interstitial pores greater than approximately 1µm. 2.3 Description of Face Centered Cubic Lattice Parameters To simulate sintered, composite electrodes across the porosity range described in Section 2.2, packing densities as high as 0.95 are required. The maximum theoretical packing density of monosized hard spheres is 0.74, corresponding to close packing with either the face centered cubic lattice or hexagonal close packed lattice. To reduce model complexity, the face centered cubic lattice is chosen to describe the particle positions, as illustrated in Figure 2-1. An alternative is to determine a stable equilibrium position for each individual particle as particles are consecutively added to a container [46]; however, the random population of a face centered cubic lattice has been chosen, similar to work by Ioselevich et al. [22].  Figure 2-1: Face centered cubic packing of hard spheres illustrating the relationship between edge length (lattice constant) and particle diameter. a d2   24 The lattice constant, a , is dictated by the size of the particles in the packing. If the particles are monosized, as in Figure 2-1, then the lattice constant is given by equation (25), where d is the diameter of each particle. da 2=  (25) Physically, ceramic particles used to fabricate real electrodes are not monosized; rather, they have a distribution about a mean diameter. In this work, the particles have been modeled with a distribution in diameter.  However, for the lattice to maintain a regular and repeating structure, two reference parameters describing the degree to which the electrode is sintered must be introduced. Specifically, these reference parameters are the reference particle contact angle, 0θ , and reference particle diameter, 0d , as illustrated in Figure 2-2. For this work a value of 12pi radians for 0θ  is adopted in accordance with Costamagna et al. [29]. It will be shown below that the reference particle diameter can be derived by specifying both the reference contact angle and the desired mean particle diameter. Using the reference contact angle and reference particle diameter, the lattice constant is defined as shown by equation (26).   25  00 cos2 θda =  (26) The reference particle diameter as illustrated in Figure 2-2 is not useful from an electrode design perspective. Generally, ceramic powders are chosen based on their mean diameter. Specific to this work, the modeled geometry represents the electrode structure post sintering. During sintering, particle diameters decrease as the material diffuses toward the regions of particle contact. The scope of this work does not include the physics of the sintering process; instead, the degree of sintering is defined by the reference particle contact angle. The mean particle diameter that exists after the sintering process is a user input. Constraints are imposed on the minimum and maximum diameter for a particle on the face centered cubic lattice. The minimum diameter is constrained to be the diameter for which a monosized packing of diameter dmin, and reference contact angle θ⁰ , has all particle contacts occurring as points, as illustrated in Figure 2-3. Without this constraint, it would be Figure 2-2: Face centered cubic lattice based on reference particle diameter and reference contact angle, θ⁰ = 12pi radians. 2 0d  0θ a   26 numerically feasible for particles to exist in the packing completely unsupported by their nearest neighbours, an unphysical situation.  The minimum diameter is therefore described in terms of the reference particle diameter and reference contact angle. 00 min cosθdd =  (27) The maximum diameter is constrained to be the diameter for which a monosized packing of diameter dmax and reference contact angle θ⁰ would be fully dense. This scenario is illustrated in Figure 2-4.   Figure 2-3: Cross section of face centered cubic packing of monosized spheres with minimum diameter. 2 mind  mind 2 mind  a   27  The maximum diameter is therefore described in terms of the reference particle diameter and contact angle by equation (28). 00 max cos2 θdd =  (28) The mean particle diameter, dmean, is defined as the arithmetic mean between the minimum and maximum permitted particle diameters. The resulting expression for the reference particle diameter is given by equation (29). ( ) 00 cos21 2 θ+ = meandd  (29) Combining equations (26) and (29), the lattice constant for a modeled electrode with variable particle diameters can be described in terms of the desired mean particle diameter according to equation (30). Figure 2-4: Face centered cubic packing with maximum particle diameters, corresponding to fully dense packing maxda =   28 ( ) meanda 21 22+=  (30) 2.4 Description of Sintered Particle Intersection The geometric details of particle intersection are required to construct a three dimensional solid model of sintered composite electrodes. Once the lattice has been defined, as in Section 2.2, the sinter neck radius, sr , defines the circumference of intersection between two particles and depends upon the intersecting particle diameters as well as the radius of neck curvature, nr , as illustrated in Figure 2-5.  The radius of neck curvature has been defined to be equal to one tenth of the radius of the reference particle, with diameter given by equation (29). All intersecting particles form sinter necks with the same radius of curvature in the model. This definition is intended to be general, in that the geometrical model is to be independent of the specific electronic and ionic conducting materials used. The circumference of the circle of intersection between two particles, having radius sr ,  is the three phase boundary length for that particular particle intersection. Therefore, the expressions developed in this chapter are also fundamental to the calculation of the geometric three phase boundary length in Chapter Three. Figure 2-5 illustrates the intersection of two particles, i and j. Half of the sinter neck is formed by revolving the shape defined by points N,M,Q, and R about line AN  by 2pi radians. The other half is formed by revolving the respective polygon for particle j. It is therefore required to find the coordinates for points N,M,Q, and R.    29  The intersection between the same two particles is illustrated in greater detail in Figure 2-6 for the purpose of defining the required modeling expressions. Figure 2-5: Intersection of two arbitrary particles, i and j, showing the shape defined by N,M,Q, and R to be revolved about line AN to create sinter neck on particle i. R Q A N M Particle j Particle i nr  RNrs =   30  Point N lies in the plane of intersection between particles i and j and does not occur at the midpoint of AB unless the two particles have equal diameter. The points A and B are the center points of particles i and j , respectively. The length of line AB  is the distance between the center coordinates [ ix , iy , iz ] and [ jx , jy , jz ] of both particles. The value of subscript i, corresponding to the particle number, ranges from one to the total number of particles in the composite electrode. The subscript j, corresponding to the neighbouring particle number with respect to particle i, ranges from one to twelve, the maximum number of nearest neighbour particles on the face centered cubic lattice. The distance, ijL , between the centers of neighbouring particles i and j is given by equation (31). Figure 2-6: Geometry of intersection between particles i and j. C A B M N Q R 4 pi  Particle j Particle i ANS ij = ABLij = RNrs =   31 ( ) ( ) ( )222 ijijijij zzyyxxL −+−+−=  (31) The distance from the center of particle i to the point on AB  in the plane of intersection, point  N, is denoted ijS , and is given by equation (32). ijS  is derived according to the geometry of intersection of two particles, i and j, as shown in Figure 2-6. Note that in the case where particles i and j have equal diameters, the intersection plane lies at the midpoint of AB . ij ijij ij L ddL S 8 4 222 +− =  (32) Once the location of the plane of intersection is known, the Pythagorean theorem is used to determine the sinter neck radius, resulting in equation (33). The sinter neck radius therefore depends on the diameter of both particle i and j. nijniins rSrddrr −−++= 222 4442 1  (33) Equations (31), (32) and (33) are sufficient to describe the geometry of particle intersection for each particle on the lattice once the mean particle diameters have been assigned. The variance in the particle diameter is determined by first establishing how large the packing must be to obtain statistically significant results, as described in the following section. 2.5 Dispersion Study for Variation in Three Phase Boundary Length by Varying Number of Lattice Sites The use of random properties destroys the deterministic nature of the modeled geometry in the sense that two geometries with the same input parameters will not be strictly identical. The solid volume fractions, total porosity, and particle size distribution can only be controlled within a certain tolerance, determined by the number of lattice sites in the model. The modeled geometries will be used to evaluate the three phase boundary length and effective diffusion coefficient in sintered composite electrodes. Therefore, the dispersion in both the calculated total three phase boundary length and total porosity has been studied as   32 the number of lattice sites increases. Figure 2-7 shows the total three phase boundary length as a function of the particle packing size for three different sets of input parameters. For each case in Figure 2-7, ten electrode models were created for each packing size on the abscissa and the total three phase boundary length was calculated according to the method described in Chapter Three.  Case 3 has the greatest amount of relative dispersion at the largest packing sizes because this case has the fewest intersections between ionic and electronic conducting particles. This is because it has the highest porosity and smallest solid volume fraction of ionic conductor. Note that equal results for Case 3 would have been obtained for the total three phase boundary length if it was the electronic conducting phase that had a solid volume fraction of 0.1. Case 1 has the largest number of three phase boundaries per unit volume because it has Figure 2-7: Dispersion of three phase boundary length for different face centered cubic lattice sizes. Mean particle diameter for Case 1 and Case 3 is 2µm. Mean particle diameter for Case 2 is 2.9µm. ● Case 1:   50.0= −e f ; 05.0=Tε ▲ Case 2:   39.0= −e f ; 37.0=Tε  ♦ Case 3:   90.0= −e f ; 55.0=Tε Chosen packing size, pN = 7813   33 the largest number of particles on the lattice and equal solid volume fractions, increasing the frequency of intersections of opposite phase particles. Case 2 represents a typical nickel- yttria stabilized zirconia (Ni-YSZ) cermet anode in a solid oxide fuel cell. The porosity is 37% and the nickel solid volume fraction is 0.39. This case has an intermediate number of three phase boundaries compared to cases 1 and 3. A convenient way of describing the packing size is to consider the number of face centered cubic unit cells that comprise an edge of the overall geometry. If the number of unit cells on an edge is ,fccn then the number of lattice sites in the packing is given by equation (34). ( ) ( )131 23 +++= fccfccfccp nnnN  (34) Based on the ten data points at a packing size with twelve unit cell lengths per edge of the overall cubic geometry envelope (7813 total lattice sites, indicated with an arrow in Figure 2-7), the scatter in three phase boundary length is within ±1.5% of the mean value of 6.95x1011 m/m3 for Case 1 and within ±3.3% of the mean value of 1.20x1011 m/m3 for Case 3. Similarly, the porosity was calculated according to the method described in section 2.7, for the minimum and maximum porosity values of 0.05 and 0.55. The porosity of the modeled electrode results from removing particles from the lattice and calculating the volume of removed particles plus the volume of the void space between particles occupying lattice sites. The porosity was calculated for a range of packing sizes, as shown in Figure 2-8. The packing size of 7813 lattice sites (or 12 unit cell lengths per side) determined in Figure 2-7 shows scatter of less than ±5% for both porosity values.   34  2.6 Parametric Solid Model of Sintered Composite Electrodes The first step in creating the modeled geometry is to define the diameter and positions of the particle centers according to the face centered cubic lattice. The particles are centered at the lattice sites and the diameters are normally distributed about the mean particle diameter with variance, 2dσ , described by equation (35). The variance is defined such that a particle having diameter less than mind or greater than maxd is not likely to occur. It was shown in Section 2.5 that the required number of lattice sites is 7813; therefore, there must be four standard deviations between the mean diameter and both the minimum and maximum diameter to statistically ensure that the particle size distribution remains between mind < meand <  maxd . 101 102 103 104 0 0.05 0.15 0.25 0.35 0.45 0.55 0.65 Number of Face Centered Cubic Lattice Sites To ta l P o ro si ty Figure 2-8: Dispersion in porosity for different face centered cubic lattice sizes. Mean particle diameter is 2µm in both cases. ●   Case 1: 55.0=Tε ▲ Case 2: 05.0=Tε   35 ( )2min2 4 1 ddmeand −=σ  (35) Each particle on the lattice must be assigned a phase, either an electronic conducting particle or ionic conducting particle. The phase assignment is random; with each particle having equal probability to be assigned an electronic conductor or ionic conductor, independent of the particle position.  If composition gradients were to be introduced then there could be deliberate bias introduced in phase assignment as a function of position. However, for this work, uniform compositions are desired to enable characterization of specific electrode compositions and microstructures. Porosity is introduced to the packing by removing particles at random. Therefore, the porosity is not strictly an input parameter to the geometry model. The fraction of unoccupied lattice sites is correlated to the desired porosity as shown in Figure 2-9. It can be seen that the mean particle diameter does not influence the correlation, and it is therefore valid over the entire mean particle diameter range under consideration from 1µm to 10µm. A cubic spline interpolation is used to evaluate the required fraction of particles to remove from the lattice to obtain the desired total porosity. The code written to create the geometry is described in detail in Appendix A. The Matlab programming language is used to generate the particle positions, diameters, and phase. Using these data, the position, diameter, and phase of the nearest neighbours for every particle in the lattice is calculated. The creation of the geometry is accomplished using the Unigraphics NXOpen Application Programming Interface, allowing the specific geometry creation commands to be written using Visual Basic.NET, in the Microsoft Visual Studio programming environment. The Visual Basic.NET code is then read into Unigraphics NX5 via the integrated Journaling feature. Because the number of particles modeled in each electrode geometry is large, 7813 particles if zero particles are removed, the file size can become very large and unmanageable. Therefore, each electrode geometry is divided into 27 sub cubes and each sub cube is modeled individually. A sample geometry is shown in Figure 2-10.    36  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Fraction of Particles Removed from FCC Lattice To ta l P o ro si ty   d m  = 1 micrometer d m  = 5 micrometer d m  = 10 micrometer Figure 2-9: Correlation of fraction of particles to remove from face centered cubic lattice to total porosity for three mean particle diameters. There are a total of 7813 lattice sites for each data point. ▲ µ1=meand ●   µm5=meand ■   µm10=meand   37  2.7 Calculation of Total Porosity The porosity of the electrode geometry is calculated by computing the cumulative sum of the solid volume for each particle and then applying equation (36), where the solid volume is Vs and the envelope volume for the geometry is VE. E s T V V −=1ε  (36) Figure 2-10: Example solid model of sintered, composite, solid oxide fuel cell electrode. The total porosity is 0.55. The mean particle diameter is 7.43 µm. The electronic conductor (black particles) solid volume fraction is 0.27.   38 The solid volume is the sum of the volumes of each of the composite geometries consisting of each spherical particle and its sinter necks. The volume of the sinter neck is obtained using a volume of revolution as shown in Figure 2-11.  The equation of the circle containing the sinter neck curvature is described by equation (37), and the location of point y for 00 yyry n ≤≤− at location x on interval [a, b] is given by equation (38). ( ) ( ) 22020 nryyxx =−+−  (37) ( ) 0202 yxxry n +−−−=  (38) The volume of revolution is therefore given by equation (39) and is numerically integrated in Matlab using the built in function quad.m, which approximates the sinter neck curvature using quadratic polynomials according to Simpson’s rule. ∫= b an dxyV 2pi  (39) The volume of the spherical cap for particle i that is subtracted from the sum of the spherical particle volume plus the sinter neck volume, to avoid double counting the cap  Figure 2-11: Sinter neck arc is revolved about x axis to form volume of revolution for purpose of porosity calculation. x y   39 volume, is illustrated in Figure 2-12 and has height .Hrh i −=  The distance AMH =  may be found using similar triangles ACN and AQM with the result in equation (40). The volume of the spherical cap, cV , for particle i is given by equation (41). ni i ji rr rSH + = ,  (40) ( )hrhV ic −= 33 1 2pi  (41)  The code written to calculate the total porosity is described in detail in Appendix B. 2.8 Verification of Total Porosity Calculation The geometry used to verify the total porosity calculation is shown in Figure 2-13. The porosity was calculated both according to the algorithm described in the previous section and Figure 2-12: Spherical cap for particle i. 2 i i d r = Spherical cap for particle i Particle i AMH = C A Q B N M ANS ij = nr   40 using the measuring tools in Unigraphics NX5. The results of the comparison are summarized in Table 2-2.  Note that particles which have center coordinates located on the envelope: either corner, edge or face, must have the solid volume of the sphere reduced to one eighth, one quarter, or one half, respectively. Furthermore, the volume of revolution of the sinter necks is calculated only over the angle for which the sinter neck is exposed to the pore space. The method used to calculate this angle is described in Chapter Three, Section 3.2. Calculation Method Solid Volume sV Total Porosity Calculated According to Equation (36) Algorithm described in Section 2.7 301 µm3 0.194 Measuring tools in Unigraphics NX5 299 µm3 0.189  The geometry model developed in this Chapter is used to describe sintered composite electrodes with variable composition, porosity, and particle size. These geometry models form the basis from which the connected three phase boundary length and effective diffusion coefficient are estimated, as described in Chapters Three and Four respectively. Figure 2-13: Geometry used to verify total porosity calculation. The envelope volume is 371µm3 Table 2-2: Summary of results for verification of total porosity calculation   41 Chapter 3 Three Phase Boundary Length in Sintered Composite Electrodes 3.1 Introduction The total geometric three phase boundary length for a geometry created using the method of Chapter Two may be calculated based on the diameter, phase, and position of each particle in the modeled electrode, and the properties of its nearest neighbours. The connected three phase boundary length is calculated using the tomographic method described in Section 3.4, using the algorithm developed in large part by the project collaborators at the University of Sherbrooke. There are two purposes to calculate the total three phase boundary length based on the continuous geometry. First, this value serves as an upper limit to the connected three phase boundary length. Second, in the digital reconstruction of the geometry from cross section images, the fraction of all identified three dimensional pixels, termed voxels, representing the three phase boundaries that are connected, is the same fraction of the total length that is connected. This provides a convenient way of estimating the connected three phase boundary length in modeled geometries where the total three phase boundary length is known. 3.2 Total Three Phase Boundary Length Calculation Method A consequence of constraining the particle centers to the face centered cubic lattice is that the particles may only intersect, and thus possess a sinter neck, at specific locations. An arbitrary particle i on the lattice has at most twelve nearest neighbours, and the sinter necks are centered on a line connecting the center of particle i to the center of neighbour particle j. Figure 3-1 illustrates the use of local Cartesian coordinate planes positioned at the center of particle i to describe the locations at which sinter necks will occur if the corresponding nearest neighbour particle is present on the lattice. For example, 1,QXYI , on a given particle i, in its first quadrant (Q1), will intersect only with 3,QXYI  on its neighbour particle,  j, if particle j is present on the lattice.   42  For the purpose of calculating the total three phase boundary length, the method developed in this work is to determine, for each particle i on the lattice, the diameter and position of each nearest neighbour particle of opposite phase that is present on the lattice. It is only the opposite phase intersections that are of importance to calculating the three phase boundary length. In the case of Figure 3-1, the diameters of the particle shown, and of all neighbour particles, have been deliberately set to mind to ensure that the sinter necks do not overlap. As described in Chapter Two, the sinter neck radius depends upon the diameter of both particle i and the neighbour particle, j, and it is possible that the sinter necks formed with other neighbour particles overlap each other. An example is shown in Figure 3-2, where 2,QXYI and Figure 3-1: Definition of sinter neck locations. 1,QXYI  represents the sinter neck located in the XY plane, quadrant one, which will only intersect with location 3,QXYI  on the neighbour particle if the latter is present on the lattice. Therefore, each sinter neck location on particle i has a unique corresponding sinter neck location for each of the possible twelve neighbour particles. 1,QXYI  2,QXYI 3,QXYI  4,QXYI 1,QYZI 2,QYZI 1,QXZI  2,QXZI X Y Z   43 1,QYZI interfere. The consequence of this is that the three phase boundary lengths of these two sinter necks cannot be calculated simply by multiplying the respective sinter neck radii by pi2  radians.   Each sinter neck may have from zero up to four interfering sinter necks. The latter case is shown in Figure 3-3, where the curve with origin at point O = [Ox, Oy, Oz] represents the sinter neck formed by particle i intersecting with an arbitrary neighbour particle of opposite phase. The four remaining curves, labeled a, b, c, and d, are sinter necks formed by the same particle i with four other neighbour particles. This scenario becomes common at lower porosities, where a large fraction of lattice sites are occupied with particles. For the case shown in Figure 3-3, the three phase boundary length is a small fraction of the overall circumference, and is illustrated with thick lines. Figure 3-2: Arbitrary particle with four neighbours showing interference between two of the sinter necks. 1,QYZI  2,QXY I 2,QYZI 4,QXZI  X Y Z   44  As discussed in Chapter Two, the particle diameters are restricted to the interval [ ]maxmin , dd , and this imposes  limitations to how the sinter necks may interfere with each other. For example, curve a and curve c may each individually interfere with the curve having origin O. Further, curve a and curve c may interfere with each other in addition to the curve with origin O. The same possibilities exist for curves b and d. It should be noted that curves a and b or curves c and d cannot interfere with each other due to the maximum diameter constraint. For example, curve c and d would only interfere if the maximum diameter on the lattice were greater than maxd , defined in Chapter Two to be equal to the lattice constant. If the diameters of particle i, with sinter neck centered at O, and the particles 1 1.5 2 2.5 0.5 1 1.5 2 2.5 x Figure 3-3: Interfering sinter necks represented by parametric curves in three dimensional space. The curve with origin at point O has four interfering sinter necks labeled a, b, c, and d. Curve a may also interfere with curve c; likewise, curve b may also interfere with curve d. Curve a cannot interfere with curve b and curve c cannot interfere with curve d. Thick lines are three phase boundary lines. Three phase boundary Three phase boundary a b c d O   45 with sinter neck curves c and d were all of maxd , then curves c and d would touch at a point, but would not overlap. Using the observations stated in the previous paragraphs, the arc length of the three phase boundary for the intersection of particle i with its neighbours may be calculated. The method to calculate the total geometric three phase boundary length is illustrated for the case of a sinter neck having two interfering sinter necks as shown in Figure 3-4.  -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Figure 3-4: Calculation of angle of overlap between sinter neck with origin at O and interfering sinter neck curves b and d. Curve b has radius bsr , and likewise, curve d has radius dsr , , representing the radii of the sinter necks at these locations. Three phase boundary v3 O v1 v2 v4 d b sr   46 In order to calculate the three phase boundary length illustrated in Figure 3-4, the angle between vectors 1v  and 2v  and between 3v and 4v  must be determined. Note that vectors 1v  and 2v are the position vectors from point O to the respective points at which the curve having origin at point O intersects with curve b.  Likewise, 3v and 4v  are the position vectors from point O to the respective points at which the curve having origin at point O intersects with curve d. The intersection points of the curve having origin O with curve b are calculated by finding the zeros of equation (42), for all points Px,k, Py,k, and Pz,k  comprising curve b, where the subscript k represents the k’th point describing the curve. ( ) ( ) ( ) szzyyxx rOPOPOPD −−+−+−= 222  (42) Equation (42) describes the distance, D, from a point on the interfering sinter neck curve, for example, curve b, to the point O minus the sinter neck radius sr . Therefore, where equation (42) is zero, the point P lies on the curve centered at O. All points P comprising curve b are calculated by representing curve b as a parametric circle in three dimensions according to equation (43), where B = [Bx, By, Bz] is the origin of curve b, the vectors av and bv are two orthogonal unit vectors spanning the plane on which curve b lies and bsr , is the radius of curve b. ( ) ba vv tsBtsP ++=,  piα α α 20for cos sin , , <≤ = = bs bs rt rs  (43) The angle between two vectors iv and jv  is denoted ji,φ  and is obtained using equation (44).         ⋅ = − ji ji vv vv1 , cosjiφ  (44) For the case in Figure 3-4, the angles 2,1φ , 4,3φ , 3,1φ and 4,2φ are calculated. If the sum of angles 2,1φ  and 4,3φ  is greater than angle 3,1φ , the two sinter necks represented by curves b   47 and d must overlap each other. Therefore, the angle over which the three phase boundary in Figure 3-4 occurs is either ( )4,32,12 φφpi +−  for ( ) 3,14,32,1 φφφ <+  or ( ) 4,24,32,12 φφφpi ++−  for ( ) 3,14,32,1 φφφ >+ . The method described above is conducted for every sinter neck between opposite phased particles in the modeled geometry and cumulatively summed to obtain the overall geometric three phase boundary length. Note that the interfering sinter necks need not be of opposite phase, as the interference occurs regardless of composition. One final check must be performed to ensure the correct total three phase boundary length. For particles occupying lattice sites on an edge, face, or corner of the outer envelope, the three phase boundary length for these particles may extend outside the envelope region. For example, Figure 3-5 highlights particles occupying edge, face, and corner lattice sites. Depending on the specific edge or face, there are specific sinter neck locations which are contained completely within the envelope volume, and specific sinter neck locations which are truncated by the envelope volume by a value of pi radians. Similarly, the corner lattice sites have all sinter neck locations within the cube envelope truncated by pi  radians.  Figure 3-5: Example of particles occupying edge, face, and corner lattice sites. Corner lattice site Edge lattice site Face lattice site   48 The method to calculate the total geometric three phase boundary length as described in this section has been programmed in Matlab such that the total three phase boundary length is calculated at the time of creating the particle positions, diameters, and phase, which are then input to the geometry creation algorithm described in Chapter Two. The code is described in detail in Appendix B. 3.3 Verification of Total Geometric Three Phase Boundary Length The algorithm to calculate the three phase boundary length described in Section 3.2 has been verified using two different test geometries. The geometries were created in Unigraphics NX5 as described in Chapter Two, and the total three phase boundary length was calculated according to the algorithm described above. The total three phase boundary length may also be measured using the Unigraphics NX5 graphical user interface by measuring the arc length of the edges of the opposite phase particles exposed to the pore region. Figure 3-6 illustrates the two test geometries and Table 3-1 summarizes the results.    (a)       (b) Figure 3-6: Three phase boundary test geometries. (a) Single FCC unit cell. (b) Eight FCC unit cells.   49 Geometry Calculation Method Total Three Phase Boundary Length Figure 3-6 (a) Algorithm described in section 3.2 13.2 µm Unigraphics NX5 graphical user interface 13.1 µm Figure 3-6 (b) Algorithm described in section 3.2 57.5 µm Unigraphics NX5 graphical user interface 57.1 µm  It is observed from the results in Table 3-1 that the algorithm described in Section 3.2 agrees with the value obtained using the Unigraphics NX5 measuring tools. The specific method used by Unigraphics to calculate arc lengths is not described in the documentation; however, a test case conducted by modeling an arc with known radius and central angle yields the correct result when measured with the same length measurement tool used to produce the results in Table 3-1. 3.4 Identification of Connected Three Phase Boundary Voxels The connected three phase boundary length is calculated using an algorithm developed by the project collaborators at the University of Sherbrooke. Their algorithm is based on work by Hoshen and Kopelman [47] and is capable of evaluating the percolation, or connectivity, of each of the three phases comprising the composite electrode. Relevant to this work, the method is executed by first converting cross section images of the modeled electrode into a three dimensional matrix describing the particle phases and positions. The cross section images consist of two parallel stacks, each containing binary images of one of the solid phases. For example, one of the stacks represents the electronic conductor as ones (white) and the other two phases are represented by zeros (black). The second stack represents the ionic conductor as ones and the other two phases as zeros. Figure 3-7 shows the first cross section image of the geometry in Figure 3-6 (b) and the resulting two binary images. Table 3-1: Summary of Total Geometric Three Phase Boundary Length Verification   50  The binary image stacks are each read into an m by n by p matrix, where m is the number of cross section images, n is the image width in pixels, and p is the image height in pixels. A third binary matrix representing the pore phase is obtained by performing the appropriate logical operation on the two solid phase matrices. The two solid phase matrices are then expanded such that a layer of ones is added to the surface of the solid phase. Then, each of the three matrices are multiplied term by term, and the resulting ones in the product matrix represent the locations at which the three phases meet; thus, these are three phase boundary voxels. As an example, Figure 3-8 shows the three phase boundary voxels present after analyzing the geometry in Figure 3-6 (a). (a)      (b)    (c) Figure 3-7: First cross section of geometry in Figure 3-6 (b). (a) Actual cross section image as generated using Unigraphics NX5. Electronic conductor is black. Ionic conductor is gray. Pore is white. (b) Binary image representing electronic conducting phase as white. (c) Binary image representing ionic conducting phase as white.   51  The noisy voxels in Figure 3-8 arise during the image capture process in Unigraphics. The ionic conducting particles are colored gray and there are black pixels randomly placed along the perimeter of the cross section of the gray particles. Although they are not visible in the Figure 3-6 (b) they are present and visible on the monitor while creating the cross section images. It is possible to filter out these noisy voxels by considering that three phase boundaries may only occur at the intersection between opposite phase particles. Figure 3-9 (a) shows that the three phase boundary voxels lie on the circumference of the intersections between opposite phase particles. Voxels which do not lie on the circumference are removed.  Figure 3-8: Three phase boundary voxels for geometry in Figure 3-6 (a). The noisy voxels that do not fall along the circumference of intersection of opposite phase particles require filtering. noise   52  Each three phase boundary voxel identified is then checked to determine whether each phase is connected to the source of the transported species. Specifically, there must be a path for electron conduction to the current collector, ion conduction to the electrolyte, and gas species diffusion to the gas channel. If these three conditions are satisfied, the three phase boundary voxel is considered connected. This identification of connected three phase boundary voxels is performed using the code developed by the project collaborators, with the result shown in Figure 3-10. (a)        (b) Figure 3-9: Three phase boundary voxels with noise removed. (a) Three phase boundary voxels to be removed, blue cubes, do not fall on the circumference of intersection between opposite phase particles, indicated by the red circles. (b) Three phase boundary voxels remaining after filtering. Note that the images in (a) and (b) are flipped vertically compared to Figure 3-8 and Figure 3-10.   53  The method described above to identify the three phase boundary voxels will find a larger number than exist physically, in addition to the noisy voxels identified in Figure 3-8. Consider an example case in Figure 3-11, where the two solid phases are adjacent and both solid phases are exposed to the pore phase. In Figure 3-12, the matrix operations corresponding to one plane in Figure 3-11, that intersects the three phase boundary line, are shown. In this case, two three-phase boundary voxels are found where only one actually exists.  Figure 3-10: Connected three phase boundary voxels for geometry shown in Figure 3-6 (a).  Figure 3-11: (a) Three phase boundary voxels superimposed on a close up image of the geometry in (b). Both images are in the same orientation as Figure 3-9. (a) (b)   54  Figure 3-12: Example scenario in which two three-phase boundary voxels are created where only one physically exists. In (a), E represents the electron conductor phase, I represents the ion conductor phase, and P represents the pore phase. (b) through (d) are the binary representations of (a) for the electron conductor phase, ion conductor phase, and pore phase, respectively. (e) and (f) are the expanded matrices in which ones are placed around the electronic conducting phase in (e) and ionic conducting phase in (f). (g) is the expansion matrix for the pore phase; however, note that ones are not placed around the pore phase. (h) is generated by multiplying matrices (e), (f) and (g) term by term.   55 In three dimensions, each solid phase is expanded in three dimensions and thus there are also cases where a three phase boundary voxel may be found although the two solid phases are not physically in contact; but rather, separated by one voxel width. This leads to the restriction that an identified three phase boundary voxel is not considered a three phase boundary voxel unless the following is true. A three phase boundary voxel must share a face with one of the solid phases and either an edge on that face with the other solid phase, or a different face (having an orthogonal normal vector to that of the first face) with the other solid phase. Filters to account for these two situations are applied to the three phase boundary voxels identified prior to identifying the connected voxels. 3.5 Estimation of Connected Three Phase Boundary Length Estimating the connected three phase boundary length as a summation of the contribution from each individual voxel is required for real electrodes where the total three phase boundary length is not known. However, for deterministic electrodes, the total three phase boundary length is known a priori. Therefore, the total three phase boundary length, totalλ , calculated according to Section 3.2, corresponds to the total number of three phase boundary voxels obtained from the input images. The fraction of the three phase boundary length that is connected,  cλ , is then the same fraction of three phase boundary voxels that are connected as shown in equation (45), where voxcN ,  is the number of connected three phase boundary voxels and voxtN ,  is the total number of three phase boundary voxels. Table 3-2 summarizes the results for the two test geometries of Figure 3-6. total voxt voxc c N N λλ , , =  (45)     56  Geometry from Figure 3-6(a) Geometry from Figure 3-6(b) Total No. Three Phase Boundary Voxels including Noise 3384 24723 Total No. Three Phase Boundary Voxels excluding Noise, voxtN ,  1880 7614 No. Connected Three Phase Boundary Voxels, voxcN ,  720 783 Connectivity Coefficient ( voxtvoxc NN ,, ) 0.383 0.103 Actual, Total Three Phase Boundary Length, totalλ † 13.2µm 57.5µm Actual, Connected Three Phase Boundary Length†† 5.00µm 6.02µm Estimated Connected Three Phase Boundary Length, cλ ††† 5.06µm 5.91µm † Calculated using method described in Section 3.2. †† Measured with Unigraphics NX5 length measurement tool. ††† Calculated according to equation (45). 3.6 Required Binary Cross Section Image Resolution The evaluation of connectivity for the two solid phases comprising the electrode cannot be accomplished for image resolutions below that for which two particles are separated by less than three voxels. In Figure 3-13(a), the cross section image has a resolution of 100 by 100 pixels. Figure 3-13(b) depicts the same cross section image with a resolution of 300 by 300 pixels. It can be seen that the two gray coloured, ionic conducting particles enclosed in the dashed circle in Figure 3-13(a) are separated by a single pixel. The code written by the project collaborators to evaluate connectivity of the solid phases uses three dimensional expansion matrices (similar in concept to the two dimensional expansion matrices shown in Figure 3-12(e) and Figure 3-12(f)) for each solid phase. Therefore, if two particles that do not physically touch in the modeled geometry are separated by fewer than three voxels, then they will be considered to be connected. The robust solution is to modify the code so that the connectivity is evaluated based upon the binary representations of the solid phases obtained directly from the images (similar in concept to the two dimensional matrices in Figure Table 3-2: Estimated Connected Three Phase Boundary Length using Connectivity Coefficient.   57 3-12(b) and Figure 3-12(c)). However, this modification has not yet been implemented. Therefore, the solution implemented in this work is to increase the image resolution such that it is unlikely for two particles, which do not physically touch in the modeled geometry, to be separated by fewer than three voxels. Figure 3-13(b) shows in two dimensions that the same two ionic conducting particles identified in Figure 3-13(a) are now separated by six pixels. Therefore, for this work, an image resolution of 300 by 300 pixels has been used for calculating the connected three phase boundary length. The number of cross section images is 300; therefore, the individual voxels are cubic.  The above argument justifies the image size of 300 by 300 pixels to resolve individual particles that do not physically touch. It does not, however, justify this image resolution as fine enough to converge to the actual three phase boundary length. The cross section images of the modeled electrode geometries have twelve face centered cubic unit cells (thirteen particle center coordinates) per edge of the geometry envelope. The image resolution relative to the size of the modeled geometry can be defined as the number of pixels that discretize the  Figure 3-13: A cross section image of the geometry used in Chapter Four, shown here at two different image resolutions. In (a), the resolution is 100 by 100 pixels. The two ionic conducting particles (gray) enclosed in the dashed circle are separated by only one pixel. In (b), the image resolution is 300 by 300 pixels. The same two ionic conducting particles are enclosed in a dashed circle and are separated by six pixels. These images section the particles at their major diameter. (a) (b)   58 edge length of the modeled geometry. Specifically, for a geometry cross section with twelve unit cells on an edge imaged at a resolution of 300 by 300 pixels, the relative image resolution is 300 / 12 = 25 pixels per lattice constant. The image resolution for the single face centered cubic unit cell test geometry in Figure 3-6(a) is 400 by 400 pixels and the image resolution for the eight face centered cubic unit cell test geometry (having two unit cells per edge of the geometry envelope) in Figure 3-6(b) is 800 by 800 pixels. Therefore, the cross section images for each of the test geometries discretize the lattice constant into 400 divisions (i.e. pixels). This resolution was chosen because it corresponds to the highest feasible image resolution, for the modeled electrode geometries having twelve unit cells per edge of the geometry envelope, of 4800 by 4800 pixels. This is the feasible upper limit due to both reasonable file storage limitations (≈ one terabyte) and the calculation limitations of the Parallel-Mammouth supercomputing cluster configured with 872 Dell PowerEdge 750 servers, each with Intel Pentium 4, 3.2GHz processors at the University of Sherbrooke. As of 2005, this was Canada’s fastest computing cluster, capable of a sustained rate of 6.888 trillion mathematical operations per second. This computing cluster is required to analyze the connectivity of each phase due to the large number of voxels present. For 4800 images with resolution of 4800 by 4800 pixels, the matrices to be analyzed have 1011 elements. Table 3-3 summarizes the results for the connected three phase boundary length for both test geometries with the equivalent resolution of Figure 3-13(b). Specifically, the image resolution used for the test geometries of Figure 3-6, in Table 3-3, is 25 pixels per lattice constant.    59  Geometry from Figure 3-6(a) Geometry from Figure 3-6(b) Total No. Three Phase Boundary Voxels without Noise, voxtN , †  115 486 No. Connected Three Phase Boundary Voxels, voxcN ,  42 54 Connectivity Coefficient 0.365 0.111 Actual, Total Three Phase Boundary Length, totalλ †† 13.2µm 57.5µm Actual, Connected Three Phase Boundary Length††† 5.00µm 6.02µm Estimated Connected Three Phase Boundary Length, cλ †††† 4.82µm 6.39µm † At the image resolution of 25 pixels per lattice constant the noisy voxels as defined in Figure 3-8 do not exist. ††   Calculated using method described in section 3.2. ††† Measured with Unigraphics NX5 measuring tool. †††† Calculated according to equation (45). The percent difference between the estimated connected three phase boundary length using equation (45) compared to the actual connected three phase boundary length as measured using the Unigraphics NX5 measuring tool is -3.6% and 6.1% for the geometries of Figure 3-6(a) and Figure 3-6(b), respectively. Figure 3-13 indicates that the image resolution of the modeled electrode geometries of 25 pixels per lattice constant is sufficient to resolve two neighbouring particles that do not physically contact each other. However, the results in Table 3-2 and Table 3-3 indicate that this resolution of 25 pixels per lattice constant does not provide the same estimated connected three phase boundary length as the maximum feasible image resolution of 400 pixels per lattice constant. The use of the lower resolution of 25 pixels per lattice constant in this work is based on the following two reasons. First, the maximum percent difference in the estimated to actual connected three phase boundary length for the two test geometries in Figure 3-6 is 6.1%, considered satisfactory for the purpose of comparing the connected three phase boundary length between modeled electrodes having different porosity, composition, Table 3-3: Estimated Connected Three Phase Boundary Length using Connectivity Coefficient.   60 and particle size. Second, the code written by the project collaborators to evaluate the connectivity of the three electrode phases uses the Matlab programming language, which requires a license for every node used on the computing cluster. Therefore, the code is to be translated to the non-proprietary C programming language to make full use of the parallel computing resources at the University of Sherbrooke. The translation from Matlab to C has not been completed at the time of writing. Therefore, for the purpose of developing the method for estimating connected three phase boundary length in modeled composite electrodes, a resolution of 25 pixels per lattice constant is used. The workstation used in this work has two Intel quad core 3.0GHz E5450 Xeon processors and 32GB of memory. Solving the connectivity of each of the three phases, using the code developed by the project collaborators, at an image resolution of 300 by 300 pixels takes approximately three days to complete.    61 Chapter 4 Modeling Effective Diffusivity in Sintered Composite Electrodes 4.1 Introduction A method to estimate the effective diffusion coefficient for modeled electrode geometries is developed in this chapter. In development of this method, the modeled electrode geometry is taken to be the cathode in a single solid oxide fuel cell and is illustrated in Figure 4-1. For the analysis, an electrolyte and anode with typical properties are used, as summarized in Table 4-1, along with the properties of the modeled cathode. The code written to calculate the effective diffusion coefficient is detailed in Appendix C.   Figure 4-1: Example modeled sintered composite cathode having properties summarized in Table 4-1. Black particles represent electronic conductor and gray particles represent ionic conductor. Black particles, electron conductor Gray particles, ion conductor   62 Component Open Porosity Composition Thickness Connected Three Phase Boundary Length† Electrolyte 0 100% Yttria stabilized zirconia (YSZ) 20µm N/A Anode 0.37 39% Ni / 61% YSZ 40µm 14x1010m/m3†† Cathode 0.35 90% LSM ††† / 10% YSZ 95µm 0.6x10 9m/m3†† † Calculated according to method developed in Chapter Three. †† Connected three phase boundary length does not include the three phase boundaries at the interface between the modeled electrode and electrolyte, only three phase boundaries within the electrode are modeled. ††† Lanthanum strontium manganite (LSM). Modeling the current distribution within a solid oxide fuel cell electrode allows calculation of the molar flux of the electrochemically reacting species, which is required in the calculation of the effective diffusion coefficient. For both the anode and cathode, the three phase boundaries that are formed by contact of the electrode to the electrolyte are not included in the model. Choosing the modeled geometry to be the cathode is arbitrary, and the method described below is equally applicable if the modeled electrode geometry were taken to be the anode. The effective diffusion coefficient, once determined, may subsequently be used in a macroscopic performance model, where the effective property accounts for the electrode geometry. At steady state, the concentration of oxygen in a binary (O2-N2) gaseous system within the open pore space of the cathode is described by Laplace’s equation, (46), where 22 NOD −  is the binary diffusion coefficient and 2O c is the local concentration of oxygen. The boundary conditions required to solve this partial differential equation are described in Section 4.4. ( ) 0 222 =∇−⋅∇ − ONO cD  (46) Table 4-1: Electrolyte, anode and modeled cathode properties used in development of the method for calculating the effective diffusion coefficient for modeled sintered composite electrodes.   63 Equation (46), as written, is only valid for electrode microstructures where the diffusive transport is dominated by the ordinary diffusion mechanism. Therefore, application of the method developed in this Chapter is limited to diffusion layers or support layers in solid oxide fuel cell electrodes for which the ordinary diffusion mechanism dominates gas transport. The objective of Chapter Four is to determine a correction factor for the bulk binary diffusion coefficient that accounts for the geometry of the sintered composite electrode. The effective binary diffusion coefficient in porous media, eff NOD 22 − , is traditionally calculated according to equation (47) [24], where ε  is the open porosity and τ  is the gas phase tortuosity factor. However, as described in Chapter One, Section 1.3.1, the gas phase tortuosity factor is not generally known for solid oxide fuel cell electrodes. 2222 NO eff NO DD −− = τ ε  (47) The bulk binary diffusion coefficient, 22 NOD − , with units of m 2 ·s-1, is calculated according to equation (48) [48]. T is the temperature in Kelvin, p is the pressure in bar, 2O V  and 2N V  are the Fuller diffusion volumes for oxygen and nitrogen, and 22 NOM − is given by equation (49), where 2O M and 2N M are the molecular weights of oxygen and nitrogen, respectively. At 1123K and 1.01bar, the bulk binary diffusion coefficient for the O2-N2 system is 2.01x10-4m2·s-1. 2 3 1 3 1 2 1 75.1 2 2222 22 00143.0 100 1     + = − − NONO NO VVpM TD  (48) 1 22 22 112 − −         += NO NO MM M  (49)   64 4.2 Image Generation and Calculation of Open Porosity The first step toward calculating the effective diffusion coefficient is to identify the open porosity within the modeled electrode geometry. The open porosity is the domain over which equation (46) is solved. The modeled electrode geometry is discretized according to the same procedure as described in Chapter Three, section 3.4; however, the image resolution is coarsened. The pore volume in an idealized electrode having spherical pores varies according to the cube of the pore radius, whereas, the three phase boundary length varies directly with the pore radius. Therefore, fewer voxels are required to model the pore space than are required to model the three phase boundary length. Table 4-2 verifies this assumption by comparing the calculated total porosity of the discretized electrode compared to the calculated porosity using the method of Section 2.7. The total porosity from the discretized geometry is obtained by dividing the total number of pore voxels identified in the binary images by the total number of voxels in the discretized geometry. It can be seen from Table 4-2 that the percent difference between these two methods increases as the total porosity of the geometry decreases. This is explained on the basis that at low porosities, a greater fraction of the porosity is in contact with the solid phase. Therefore, the discretization error introduced by assigning a voxel to either the pore phase or solid phase is more significant at lower porosity. The percent difference may be reduced by using a greater number of cross section images together with a finer image resolution; however, the computational cost also increases as the number of pore voxels increases.  Total Porosity using Discretized Geometry Total Porosity using Method from Section 2.7 Percent Difference Example Geometry 1 0.128 0.121 5.8% Example Geometry 2 0.238 0.229 3.9% Example Geometry 3 0.557 0.550 1.3%  Table 4-2: Verification of total porosity for three different electrode geometries discretized using 100 cross section images, each image having 10000 pixels.   65 Figure 4-2 illustrates the image resolution of 100 by 100 pixels, for a total of 10000 pixels per image, used as a starting point for evaluating the effective diffusivity of the modeled electrodes. There are 100 cross section images having the same resolution. The binary image stacks are each read into a three dimensional matrix, where each matrix element contains a one or zero depending on if the pixel in the binary image is a one (white) or zero (black). For example, if the three dimensional matrix for the electron conducting solid phase has a one in a particular matrix element, then the electron conducting solid phase is present at that location. The three dimensional matrix representing the total porosity, MP, is obtained using the logical equation (50) where Me and Mi are the three dimensional matrices for the electron conducting solid phase and ion conducting solid phase, respectively.  ) OR NOT( ieP MMM =  (50) The open porosity is identified using the algorithm developed by the project collaborators and the three dimensional element indices for the open pore voxels are returned. As an example, Figure 4-3 illustrates a pore voxel having index [5, 1, 1] corresponding to the x, y, and z indices, respectively. If the plane of the electrode surface is the xy plane at z equal to one, where one represents the first layer of voxels parallel to the xy plane, then the pore voxel with index [5, 1, 1] is open as it is in direct contact with the gas channel. All pore  Figure 4-2: Example cross section image of geometry in Figure 4-1, used for determining open porosity. (a) Modeled geometry cross section where black particles represent electron conducting solid phase, gray particles represent ion conducting solid phase and white represents pore. (b) Electron conducting solid phase represented by white pixels. (c) Ion conducting solid phase represented by white pixels. (a) (b) (c)   66 voxels connected to the xy plane at z equal to one are open pore voxels and the matrix index defines the finite element positions, upon which the finite element mesh is generated.  4.3 Finite Element Mesh Generation The finite element mesh is obtained directly from the indices of the open pore voxels. The nodes are the Cartesian coordinates of the eight vertices of each cubic voxel, and are numbered according to the convention defined in Figure 4-4. The coordinates for node one are obtained by subtracting one from the x, y, and z element indices and multiplying by the physical length of the voxel edge. For example, if the edge length of the face centered cubic unit cell in Figure 4-3 is 2.39µm and there are ten voxels per edge, then the length of a single voxel edge is 0.239µm. Therefore, the physical coordinates of node number one for the voxel with index [5, 1, 1] in Figure 4-3 is the point [0.96µm, 0µm, 0µm]. Once the first node coordinates are obtained for each open pore voxel, the remaining seven node coordinates are  Figure 4-3: Discretized face centered cubic unit cell showing index of arbitrary pore voxel shown in red wireframe having index [5, 1, 1]. Note that pore voxels in this image are not shown and are instead represented by the absence of the either of the two solid phases. Black voxels are the electron conducting phase. Grey voxels are the ion conducting phase. Location of pore voxel with index [5, 1, 1] x y z First layer of voxels in xy plane is adjacent to gas channel   67 obtained by knowing the voxel edge lengths, vlx, vly, and vlz, which result from the chosen image resolution and spacing between images. The matrix containing the node coordinates has size nn by three, where nn is the number of nodes in the mesh and the three columns contain the x, y, and z node coordinates, respectively.   Following the calculation of the node coordinates, the element coordination matrix must be determined. The coordination matrix has size np by eight, where np is the number of finite elements (open pore voxels). The columns contain the row numbers of the node matrix for the eight nodes defining the element vertices. The coordination matrix therefore specifies which nodes define each element. This initial mesh, consisting of linear brick elements, is then converted to quadratic brick elements within Comsol Multiphysics. An example of a partial mesh containing the first 25000 elements of the modeled cathode depicted in Figure 4-1 is shown in Figure 4-5(a). It can be seen in Figure 4-5(a) that the first layer of elements, representing the open porosity adjacent to the gas channel, consists of disconnected pore networks. Because this location represents the electrode surface at the gas channel, all pores are physically connected through the gas channel. Therefore, a layer of elements is inserted to connect all pore networks at this plane as shown in Figure 4-5(b). Node Three Node Four Node Five Node Six Node Seven Node Eight Node One Node Two x y z Figure 4-4: Node configuration for single open pore voxel (finite element). vlx vly vlz   68   Figure 4-5: (a) First 25000 mesh elements showing disconnected elements at location of gas channel. The full mesh contains 348101 elements. (b) Same 25000 mesh elements as in (a) but with an additional layer of 10000 elements inserted to simulate the gas channel and to connect all pore networks. The direction of increasing z is the direction of the electrolyte. The direction of decreasing z is the direction of the current collector and gas channel. Direction of electrolyte Direction of current collector / gas channel First Layer of elements in the xy plane represents the gas channel (a) (b) z z z -z x x y y   69 4.4 Boundary Conditions The mesh nodes for the open porosity of the entire cathode depicted in Figure 4-1, for which the partial mesh is shown in Figure 4-5, are shown in Figure 4-6. The Dirichlet boundary condition fixes the concentration of oxygen at the bounding surface representing the gas channel in the xy plane at z equal to zero. The oxygen concentration at the gas channel is assumed constant, representing the cathode of a typical research button cell exposed to excess air, and fixed at the atmospheric partial pressure of 0.21atm, which is equivalent to an oxygen concentration, 0 2Oc , of 2.279 mol·m -3 at 1123K. The Dirichlet boundary condition is therefore given by equation (51). ( ) 0 22 0,, OO cyxc =  (51) Oxygen molecules diffuse from the gas channel through the open porosity to a connected three phase boundary, where they are reduced to oxygen ions. The connected three phase boundary locations are also shown in Figure 4-6, corresponding to the locations of the non- zero Neumann boundary conditions, which specify the concentration gradient normal to the boundary surfaces. The boundary surfaces are comprised of quadrilateral face elements which are the boundary faces of the brick elements representing an open pore voxel. Therefore, to identify the boundary surfaces that are to have a non-zero Neumann boundary condition, corresponding to the consumption of oxygen, the faces of the connected three phase boundary voxels that have the same four nodes as a quadrilateral face element must be identified. This is accomplished by searching for the boundary surface numbers assigned by Comsol Multiphysics that have face elements with the same node coordinates as the connected three phase boundary voxel faces. The exception is that the faces of the connected three phase boundary voxels in the zx plane and zy plane at the extreme y and x coordinates, respectively, are considered insulated. As described in Chapter Three, the connected three phase boundaries are identified by discretizing the modeled geometry using 300 cross section images, each image having a resolution of 300 by 300 pixels. This resolution is sufficient to resolve two particles that do not physically touch in the modeled geometry. Because the modeled geometry is discretized   70 using 100 cross section images, each image having 100 by 100 pixels, for the calculation of the effective diffusion coefficient, there may be more connected three phase boundary voxels identified than exist physically. Therefore, the positions of the connected three phase boundary voxels are identified using 300 cross section images (with image resolution equal to 300 by 300 pixels) and are then converted to their equivalent positions in a three dimensional matrix having 1003 elements.   The molar flux of oxygen molecules, 2O N , within the open porosity is described using Fick’s law as shown in equation (52). 2222 ONOO cDN ∇−= −  (52)  Figure 4-6: (a) Entire mesh consisting of 556375 nodes (for the initially linear elements), shown as blue dots, representing the open porosity of the modeled cathode in Figure 4-1. The Dirichlet boundary condition is applied to all nodes in the xy plane at z equal to zero. The locations of Neumann boundary conditions, shown as red dots, are applied to the boundary surfaces containing a connected three phase boundary voxel. All remaining boundary surfaces are insulated. x z z = 0 Direction of electrolyte Direction of current collector / gas channel z -z Neumann boundary conditions are applied at the locations of the connected three phase boundary voxels Dirichlet boundary condition at z equal to zero   71 Therefore, the form of the Neumann boundary condition is given by equation (53). 22 2 2 NO O O D N c − −=∇  (53) The molar flux of oxygen is directly linked to the production of oxygen ions according to Faraday’s law, equation (54), where ( )ηFi  is the Faradaic current density as a function of the local overpotential, η , nel is the number of moles of electrons transferred per mole of oxygen (i.e. nel = 4) and F is Faraday’s constant, equal to 96485 C·molel-1. ( ) Fn iN el F O η = 2  (54) 4.4.1 Distribution of Faradaic Current Within the cathode, oxygen ions are produced according to equation (2). Therefore, ionic current is consumed as a result of the convention that current is positive for positively charged species. The Faradaic current within the cathode is therefore negative and the reaction is proceeding in the forward direction as written in equation (2). The distribution of Faradaic current, for the cathode represented in Figure 4-1, is shown in Figure 4-7, and is calculated using a one-dimensional version of the model developed by Gazzarri and Kesler [17] according to equation (55). oi  is the exchange current density in A·m-2, caα  is the charge transfer coefficient of the cathode in the anodic direction (opposite direction to that for which equation (2) is written), ccα  is the charge transfer coefficient of the cathode in the cathodic direction (same direction as that for which equation (2) is written), and R is the universal gas constant, equal to 8.314J·mol-1·K-1. −− →+ 22 22 1 OeO  (2) ( ) ( ) ( )               −      −= ηαηαη RT F c c RT Fii cc O O caoF expexp 0 2 2  (55)   72 It can be seen in Figure 4-7 that the magnitude of the Faradaic current density, and therefore of the molar oxygen flux, increases as the z coordinate nears the cathode – electrolyte interface. The ionic current has greater magnitude near the cathode – electrolyte interface because of the lower ionic conductivity of the yttria stabilized zirconia compared to the electronic conductivity of the lanthanum strontium manganite, as described in Chapter One. The Faradaic current density ranges from -220A·m-2 to -222A·m-2, as shown in Figure 4-7, and is low for a solid oxide fuel cell operating at 0.7V. In general, the Faradaic current density for an electrode having properties listed in Table 4-1 would be expected to be larger by one or two orders of magnitude. The Faradaic current density is lower here because the three phase boundaries that are formed through contact between the electrode and electrolyte are not modeled. The connected three phase boundary voxels at this location would account for a significant fraction of the total number of connected three phase boundaries because the ion conducting phase is not percolated throughout the electrode. The molar oxygen flux at connected three phase boundaries is therefore lower than if the three phase boundaries formed at the electrode-electrolyte interface were considered. However, for modeled geometries where the open porosity is fully percolated, it is not expected that the effective diffusion coefficient is dependent on the value of the molar oxygen flux at the connected three phase boundaries. This is because the concentration gradient and the overall average diffusive flux are expected to decrease in the same proportion, and the correction factor is a geometric property when the porosity is fully percolated.   73  The geometry for which the effective diffusion coefficient is desired has three spatial dimensions; however, the Faradaic current density is known in one spatial dimension, that being the direction normal to the plane of the electrode. The electrode thickness in a solid oxide fuel cell is generally on the order of 20-100µm, whereas the planar width, or diameter in the case of a typical research button cell, is on the order of 25mm or greater. Therefore, the in-plane conduction of oxygen ions becomes negligible compared to conduction of oxygen ions normal to the plane of the electrode [49]. For this reason, the magnitude of the molar oxygen flux at a given point z within the electrode is dictated by the Faradaic current density at point z in the electrode. This molar flux is applied normal to each boundary surface containing a connected three phase boundary voxel at location z, where the normal to the boundary surface may be in the ±x, ±y or ±z directions. Figure 4-7: Faradaic current distribution within LSM / YSZ composite cathode depicted in Figure 4-1. Standard electrolyte and anode properties are used as depicted in Table 4-1. At z/L equal to one is the cathode – electrolyte interface. At z/L equal to zero is the cathode – gas channel interface. L is the cathode thickness, equal to 95µm. Cathode – electrolyte interface Cathode – gas channel interface   74 As described in Chapter One, the surface area available for electrochemical reaction is equal to the three phase boundary length multiplied by the three phase boundary width, taken to be 0.05µm. However, the total surface area of the bounding surfaces, containing connected three phase boundary voxels, over which the Neumann boundary conditions are applied in the modeled geometry, is not equal to the total surface area available for electrochemical reaction. Therefore, the Faradaic current at each bounding surface containing a connected three phase boundary must be appropriately scaled. The scaling of the Faradaic current is based on the equality in equation (56). This equality is applied for each layer of voxels parallel to the xy plane. zFi ,  is the Faradaic current density at the midpoint of layer z as determined from equation (55). The midpoint of layer z corresponds to the z coordinate of the centroid of the open pore voxel. zvoxcN ,,  is the number of connected three phase boundary voxels in layer z and voxtN ,  is the total number of three phase boundary voxels in the modeled electrode. Therefore, the ratio of zvoxcN ,,  to voxtN ,  represents the ratio of the total three phase boundary length that is connected and located within layer z. The connected surface area within layer z is obtained by multiplying the connected three phase boundary length (with unit m/m3) by the width of the three phase boundary, tpbw ,  and the modeled electrode envelope volume, VE. Equation (56) is used to determine the equivalent, or scaled, Faradaic current density, zFsi , , that is to be applied to each of the nface,z boundary surfaces having centroid located within layer z, where each surface i has area Si. ( )( )( )        =        ∑ = zfacen i izFsEtpbtpb voxt zvoxc zF SiVwN N i , 1 , , ,, , λ  (56) The Neumann boundary condition, equation (53), is therefore specified for all bounding surfaces containing a connected three phase boundary voxel by applying Faraday’s law, equation (54), to the scaled Faradaic current density, zFsi , . 4.5 Evaluation of Effective Diffusion Coefficient Solving equation (46), subject to the boundary conditions given by equations (51) and (53), provides the concentration of oxygen over the open porosity. The average concentration   75 in the xy plane for each voxel layer, z, is plotted for the cathode in Figure 4-1, and shown in Figure 4-8. The difference in concentration between the gas channel and the electrolyte is 3.6x10-6mol/m3 and is only visible within the inset of Figure 4-8.  The concentration gradient in the z direction is approximated using equation (57), where zmaxave,,2O c is the average concentration at the maximum z coordinate, maxz , and 0 2Oc  is the reference oxygen concentration at the gas channel. The concentration gradient in Figure 4-8 is -0.0377mol/m4, where maxz is 95µm and minz  is zero. Figure 4-8: Oxygen concentration for modeled cathode in Figure 4-1. A position of zero represents the location of the gas channel. A position of 95x10-6m represents the interface with the electrolyte.  Inset: zoom in on concentration profile showing a difference of 3.6x10-6mol/m3 in the concentration between the gas channel and the cathode-electrolyte interface. 3.6x10-6mol/m3 0.50 1.00 1.50 2.00 2.50 3.00 O xy ge n  Co n ce n tra tio n  [m o l/m 3 ] 95µm   76 minmax 0 zmaxave,, zz 222 − − ≈ ∂ ∂ OOO cc z c  (57) A one dimensional effective diffusion coefficient is desired in the z direction, normal to the plane of the electrode, containing the geometric details of the porous electrode. The superficial diffusive flux at location z, zavgON ,,2 , is the z component of the diffusive flux integrated over the cross sectional area of the open pore space at location z, zporeA , , and normalized by the projected area of the electrode, projA , as shown in equation (58). Figure 4-9 illustrates the superficial diffusive flux at each layer of finite elements parallel to the xy plane. Note that equation (57) refers to the estimated concentration gradient over the entire electrode thickness, whereas, equation (58) uses the local oxygen concentration gradient in the z direction. ( ) zporeA ONO proj zavgO dA z c D A N zpore ,,, , 2 222 1 ∫ ∂ ∂ −= −  (58)   Figure 4-9: Variation in superficial diffusive flux through the thickness of the modeled electrode. Mean superficial diffusive flux, avgON ,2 , is 0.97x10 -6mol/m2/s   77 Note that the superficial diffusive flux decreases near the interface with the electrolyte, in the vicinity of the connected three phase boundaries. The connected three phase boundaries extend approximately 20µm into the electrode from the electrolyte for the geometry of Figure 4-1. In this region, oxygen is being consumed, and therefore, the diffusive flux is not constant. However, a one dimensional, constant, property is desired and therefore the arithmetic mean of the superficial diffusive flux, avgON ,2 , is calculated. It is shown in Figure 4-9 that the mean superficial diffusive flux is 0.97x10-6mol/m2/s. Alternatively, a composite curve, consisting of curves N1 and N2, was fitted to the data in Figure 4-9, as shown in Figure 4-10. This curve was then integrated over the length from the gas channel, at z0, to the electrolyte, at z2, and divided by the thickness of the electrode, L, as shown in equation (59). This weighted mean superficial diffusive flux, weightedON ,2 , is 0.98x10-6mol/m2/s, which does not differ significantly from the arithmetic mean. Therefore, the arithmetic mean is used for simplicity.  Figure 4-10: Curves N1 and N2 fitted to the superficial diffusive flux. Curve N1 is constant on the interval [z0, z1] and equal to the mean of the superficial diffusive flux on this interval. Curve N2 is a fourth order polynomial fitted to the superficial diffusive flux on the interval [z1, z2]. 10 6 1 1000.1 zzzN ≤≤×= −  3 5 9 21 23 2 1060.1 2.52 1086.5 1020.2 −×= −= ×= ×−= ≤≤+++= d c b a zzzdczbzazN  2N 1N 1z  2z  0z   78     += ∫ ∫ 1 0 2 1 2 21, 1 z z z z weightedO dzNdzNL N  (59) Fick’s law is used to relate the mean superficial diffusive flux to the concentration gradient obtained using equation (57). The effective diffusion coefficient is therefore calculated according to equation (60), resulting in a value of 2.61x10-5m2/s. 1 minmax 0 zmaxave,, , zz 22 222 − −         − − = OO avgO eff NO cc ND  (60) The correction factor in equation (47) is the ratio of open porosity to gas phase tortuosity factor, denoted here by Dx , in equation (61). The correction factor applied to the bulk diffusion coefficient therefore has a value of 0.13 for the geometry in Figure 4-1. The open porosity is 0.34 and the gas phase tortuosity factor therefore has a value of 2.7. 22 22 NO eff NO D D D x − − =  (61) 4.6 Finite Element Mesh Refinement The effective diffusion coefficient was calculated to be 2.61x10-5m2/s in Section 4.5. The number of degrees of freedom in the finite element problem was 3640942, equal to the number of nodes for the 358101 quadratic brick elements. Rather than increase the number of elements to determine if the mesh size is large enough, the modeled geometry of Figure 4-1 was discretized first using 76 images, each image having a resolution of 76 by 76 pixels, and was then discretized again using 50 images, each image having a resolution of 50 by 50 pixels. This procedure is not strictly a mesh refinement for the reason that the geometry of the open pore space changes slightly as the image resolution is coarsened. However, the finite element mesh is directly obtained from the cross section images, and in all cases, the finite element mesh discretizes the same modeled geometry. Figure 4-11 shows the same cross section for each of the three image resolutions of 100 by 100 pixels, 76 by 76 pixels, and 50 by 50 pixels.   79  Table 4-3 compares the calculated effective diffusion coefficient for each image resolution in Figure 4-11. Compared to mesh case A, mesh case B results in a mesh with approximately half the number of elements; however, the effective diffusion coefficient only differs by 0.7%. Compared to mesh case A, mesh case C has approximately one seventh the number of elements and the effective diffusion coefficient differs by 12%. The open porosity summarized in Table 4-3 results from the number of open pore voxels, equal to the number of open pore finite elements less the number of gas channel finite elements, divided by the total number of voxels contained within the geometry envelope. The open porosity for mesh case B agrees with that for mesh case A to two significant figures. However, the open porosity for mesh case C differs from mesh case A by 2%, indicating that mesh case C results in a more significant change in the finite element geometry than does mesh case B. For  this work, the effective diffusion coefficient is therefore calculated by discretizing the modeled electrode geometries into 76 images, each image having a resolution of 76 by 76 pixels.   Figure 4-11: Image resolution of (a) 100 by 100 pixels, (b) 76 by 76 pixels, and (c) 50 by 50 pixels. (a) (b) (c)   80  Mesh Case Image Resolution No. Elements Open Porosity Effective Diffusion Coefficient Percent Difference in Effective Diffusion Coefficient† A 100 by 100 pixels 358101 0.35 2.61x10 -5m2/s N/A B 76 by 76 pixels 159942 0.35 2.59x10 -5m2/s -0.7% C 50 by 50 pixels 48335 0.37 2.94x10 -5m2/s 12% †Percent difference compared to effective diffusion coefficient calculated using mesh based on image resolution of 10000 pixels. 4.7 Verification Verification for the calculation of effective diffusivity was accomplished by computing the effective diffusion coefficient for an idealized electrode structure consisting of four straight rectangular pores as depicted in Figure 4-12(a). This geometry was sectioned into 56 cross section images, each image having a resolution of 56 by 56 pixels. This number of images, and pixel resolution, was chosen such that each rectangular pore is evenly discretized into 56 voxels in the z direction, 16 voxels in the x direction, and 16 voxels in the y direction. The mesh created according to Section 4.3 for the open porosity is shown in Figure 4-12(b). Note that a layer of elements has been inserted to connect the open pores at the location of the gas channel. Table 4-3: Summary of effective diffusion coefficient calculated using three different mesh densities based on three different image resolutions as shown in Figure 4-11.   81   The Dirichlet boundary condition for the concentration of oxygen, 0 2O c , is specified over the xy plane at z equal to zero, and has the value 2.279mol/m3. The Neumann boundary conditions are specified at the xy plane at z equal to 14µm. The molar oxygen flux, 2O N , at the location of the Neumann boundary conditions corresponds to a typical superficial current density in a solid oxide fuel cell of 10000A/m2. The superficial current density is the total output current normalized by the projected area of the electrode. The verification geometry of Figure 4-12(a) has a projected area of 1.96x10-10m2. Therefore, the current density normalized to the total area, 10000A/m2, corresponds to a local current density at the ends of the four pores, at z=14µm, of 30625A/m2.  Faraday’s law is used to relate this current density to a molar oxygen flux,  2O N , with magnitude of 0.0794mol/m2/s. Note that 2O N  is the one dimensional analog to that given by equation (58) because the diffusive flux is constant along  Figure 4-12: (a) Solid model of verification geometry. The height in the z direction is 14mm, and the width in the x and y directions are also 14µm. (b) Mesh of pore space for verification geometry, including layer of elements representing gas channel. Mesh contains 60480 quadratic brick elements including the 3136 elements representing the gas channel. (a) (b) z x Z=14µm Z=0 4µm 4µm Gas channel is xy plane at z equal to zero   82 the length of the pore. The gaseous species within the pore space is oxygen and nitrogen with a bulk diffusion coefficient, 22 NO D − , at 1123K of 2.011x10-4m2/s. The analytical solution is obtained by considering one of the rectangular pores as a one dimensional line extending from z equal to zero to z equal to 14µm. Laplace’s equation, (62), is then solved subject to the boundary conditions given by equations (63) and (64), with the solution given by equation (65). The concentration of oxygen at z equal to 14µm is 2.2735mol/m3. Note that the concentration is the same at the end of each of the four pores. The concentration gradient has a value of -395mol/m4. 02 2 2 = dz cd O  (62) ( ) 0 22 0 OO cc =  (63) 22 22 µm14 NO O z O D N dz dc − = − =  (64) ( ) 0 2 22 2 2 O NO O O czD N zc + − = −  (65) The superficial diffusive flux,  avgON ,2 , is obtained using equation (66), which is analogous to equation (58). Equation (66) normalizes the molar oxygen flux at z equal to 14µm by the ratio of the actual pore cross section area to the projected electrode area. 12 4 1 ,,, smmol 0260.01 22 −− = ⋅⋅== ∑ i iporezO proj Oavg ANA N  (66) The effective diffusion coefficient is obtained from Fick’s law using the concentration gradient of -395mol/m4 and superficial diffusive flux of 0.0260mol/m2/s, resulting in an effective diffusion coefficient of 6.61x10-5m2/s. Using equation (47), where the open porosity is 0.327, the tortuosity is calculated to be unity. This is the value expected for straight parallel pores because the average diffusion length is the same as the electrode thickness.   83 Therefore, for this case, the correction factor applied to the bulk diffusion coefficient to obtain the effective diffusion coefficient is simply the open porosity of the electrode. The numerical solution for the oxygen concentration within the pore space, subject to the boundary conditions given by equations (63) and (64), is shown in Figure 4-13.   The concentration at z equal to 14µm agrees with the analytical solution of 2.2735mol/m3. The concentration gradient is -391mol/m4 and the resulting effective diffusion coefficient is 6.63x10-5m2/s. This value differs from the effective diffusion coefficient of 6.61x10-5m2/s calculated analytically by 0.3%. Figure 4-13: Numerical solution for oxygen concentration in verification geometry. The concentration at z equal to 14µm is 2.2735mol/m3. 2O c  [mol/m3] z x 2.2735 2.274 2.275 2.276 2.277 2.278 2.279   84 Chapter 5 On the Connected Three Phase Boundary Length and Effective Diffusion Coefficient in Sintered Composite Electrodes 5.1 Introduction The methods developed in Chapter Three to determine the connected three phase boundary length and in Chapter Four to determine the effective diffusion coefficient are applied to specific examples of modeled sintered composite electrodes. In addition to illustrating the application of these methods, the objective is to gain insight into the resulting connected three phase boundary length and effective diffusivity for different electrode designs. The current methods available in the literature do not allow determination of these quantities for a wide range of electrode properties, in part because composite electrode layers that are feasible within a graded electrode are not necessarily feasible as the only layer comprising an electrode. There has been little modeling work on solid oxide fuel cells with functionally graded electrodes. Therefore, example electrodes that have a diverse range of properties are considered in this Chapter, and the merit of each as a layer in a graded solid oxide fuel cell is considered. 5.2 Connected Three Phase Boundary Length 5.2.1 Example Modeled Electrodes Figure 5-1 illustrates the five example electrodes used to compare the connected three phase boundary length between modeled composite electrodes with different properties. The total porosity, electronic conductor solid volume fraction, and mean particle diameter are the defining electrode properties and are summarized in Table 5-1, along with the connected three phase boundary length and the connectivity coefficient.   85    Example No. Total Porosity Solid Volume Fraction of Electronic Conductor Mean Particle Diameter Open Porosity Connectivity Coefficient Connected Three Phase Boundary Length I 37% 0.39 2.93µm 37% 0.52 13.9x1010m/m3 II 19% 0.44 10.0µm 19% 0.79 2.46x1010m/m3 III 41% 0.10 4.21µm 41% 0.04 0.18x1010m/m3 IV 35% 0.90 6.79µm 35% 0.03 0.06x1010m/m3 V 12% 0.61 8.07µm 7% 0.40 1.80x1010m/m3 Table 5-1: Summary of properties defining five examples of modeled electrodes. Also summarized are porosity, three phase boundary length, and the corresponding connectivity. I II III IV V Figure 5-1: Example modeled sintered composite electrodes. The properties for each of the five example geometries are summarized in Table 5-1.   86 Example I has the properties of a typical nickel-yttria stabilized zirconia (Ni-YSZ) anode. The porosity of 37% is fully open, the solid volume fraction of the electron conducting phase is 39%, and the mean particle diameter is 2.93µm. The fraction of the total three phase boundary length that is connected is 52%, for a connected three phase boundary length of 1.39x1011m/m3. The three phase boundary voxels that are present in this geometry are shown in Figure 5-2(a). The connected three phase boundary voxels are shown in Figure 5-2(b), where it can be seen that there are regions of disconnected three phase boundary voxels. If all three phases are fully percolated, then the fraction of three phase boundaries that are connected is one. In the case of example I, the porosity is fully percolated and the ionic conducting phase is almost certainly percolated, because this phase is present in greater proportion than the electron conducting phase. Therefore, the electron conducting phase is likely not fully connected.     Figure 5-2: (a) Total number of three phase boundary voxels for example electrode I. (b) Connected three phase boundary voxels for example electrode I. An example region of disconnected three phase boundary voxels is shown in the dashed circle. The gas channel / current collector is the xy plane at z equal to zero. Example region of disconnected three phase boundary voxels z x x z (b) (a)   87 The electron conducting phase is confirmed to be the phase that is disconnected in Figure 5-3, where the ion conducting phase has been suppressed to reveal the electron conducting phase. The electronic conducting particles that do not form a connected path have been manually removed. It can be seen that the region of disconnected electron conducting particles corresponds to the region of disconnected three phase boundaries that is highlighted in Figure 5-2(b).     Examples I and II are similar in that the porosity if fully percolated and the   The solid volume fraction of electronic conducting particles of 0.39 is close to the measured value for the upper percolation threshold for nickel in a Ni-YSZ anode of approximately 38% [50]. The lower percolation threshold for sites on a face centered cubic lattice is 0.12 [45], and the fraction of lattice sites in example electrode I occupied by the electronic conductor phase is 0.25. Therefore, the electronic conducting phase is expected to (a) (b) Ionic conducting particles suppressed from lattice in this region Ionic conducting particles and disconnected electronic conducting particles suppressed from lattice in this region Figure 5-3: Example electrode I wit  the ionic conducting particles suppressed to reveal the electronic conducting particles. (a) All electronic conducting particles present. (b) Only connected electronic conducting particles present. The gas channel / current collector is the xy plane at z equal to zero. z x x z   88 be above the lower percolation threshold and either at or near the upper percolation threshold. Two possibilities for the observed disconnected regions in example electrode I are offered. First, a physical description is possible whereby the electron conducting phase is constricted by the oxygen ion conducting phase. Assignment of the phase of each particle in the model is random, but the distribution of phases is not necessarily uniform in a finite volume. Therefore, in regions of agglomerated ion conducting particles, there will not only be a reduction in the number of three phase boundaries due to the presence of the agglomerate, but it may be possible for the agglomerate to constrict the second solid phase and disconnect larger regions of three phase boundaries. A second possibility for the observed disconnected three phase boundary regions is that the lattice size of the model does not have the same percolation threshold as an infinite face centered cubic lattice of 0.19. The lattice size was chosen, as a starting point, based on the dispersion in total three phase boundary length being less than ±5% of the mean value obtained from ten numerically generated geometries at each lattice size, over a range of lattice sizes as described in Chapter Two. Therefore, a second dispersion study using the connected three phase boundary length is required to assess the repeatability in connected three phase boundary length for different models with the same input parameters. Example II is similar to example I in that the porosity of 19% is fully percolated and the two solid phases are above the theoretical upper site percolation threshold for a face centered cubic lattice of 0.12. The fraction of three phase boundaries that are connected is 79%, which is not as high as expected for three fully percolated phases; however, it is greater than the connectivity coefficient for example electrode I. The connected three phase boundary voxels for electrode II are shown in Figure 5-4. There do not appear to be any large disconnected regions; however, it is possible that near the faces of the geometry envelope, there are disconnected three phase boundary voxels. This reinforces the requirement to conduct a dispersion study using the connected three phase boundary length to determine the minimum lattice size, and determine the required number of lattice sites to obtain a connectivity coefficient of approximately one.   89  Although the fraction of three phase boundaries that are connected in example I is lower than in example II, the connected three phase boundary length of example I (13.9x1010m/m3) is larger than that of example II (2.46x1010m/m3). The reason is two-fold. First, the mean particle diameter of 2.93µm for example I is smaller by approximately a factor of three compared to the mean particle diameter of 10µm for example II. On a volumetric basis, if the particle diameter is decreased by a factor of three, then the three phase boundary length per unit volume increases by a factor of nine. Second, the lower porosity of example II reduces the total number of three phase boundaries present in the modeled electrode. Examples III and IV are similar in that one of the solid phases is present at a composition of 0.10, resulting in the electron conducting phase for example III and ion conducting phase for example IV not being percolated. Therefore, the connected three phase boundaries are confined to the interface with the gas channel for example III and the electrolyte for example IV, as shown in Figure 5-5. Figure 5-4: Connected three phase boundary voxels of example electrode II. The gas channel / current collector is the xy plane at z equal to zero. z x   90  Although both electrode III and IV have low fractions of connected three phase boundaries, electrode III has an even lower electrochemical activity than electrode IV. This is because the three phase boundaries are confined to the location of the gas channel, resulting in large ohmic losses due to the conduction path for the oxygen ions being approximately equal to the electrode thickness. In contrast, example IV has connected three phase boundaries confined to the interface with the electrolyte, similarly to traditional solid oxide fuel cell cathodes consisting of lanthanum strontium manganite (LSM) as the only solid phase. The connected three phase boundary length of example IV (0.06x1010m/m3) is less than the connected three phase boundary length of example III (0.18x1010m/m3) because example III has a smaller particle diameter. Example V is the only case presented that does not have fully percolated porosity. The open porosity is 7%, compared to the total porosity of 12%. However, the fraction of connected three phase boundaries is relatively large, 40%. This result occurs because the porosity is comprised both of the particles missing from the lattice and the void space  Figure 5-5: Connected three phase boundary voxels for (a) example electrode III and (b) example electrode IV. The electrolyte is the xy plane at z equal to zero. (a) (b) x x z z   91 between the solid particles present on the lattice. The missing particles on the lattice are representative of the voids in sintered electrodes due to burn out of pore former particles. The void space between the solid particles represents incomplete sintering, or also the reduction of NiO to Ni in the case of the anode. The significance is that a relatively high connectivity coefficient may result for electrodes having porosities below the upper percolation threshold, compared to electrodes where the solid volume fraction of one of the solid phases is below its upper percolation threshold. The upper percolation threshold for porosity is approximately 20% total porosity, and the lower percolation threshold is approximately 10% total porosity, as shown in Figure 5-6, for modeled electrodes within the design space described in Chapter Two. Therefore, it may be feasible to have a thin layer next to the electrolyte with porosity below the upper percolation threshold. The three phase boundary length could further be increased if the particle diameters are reduced.  0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.00 0.10 0.20 0.30 0.40 0.50 0.60 O p e n  P o ro si ty Total Porosity Figure 5-6: Open porosity over the range from 0.05 total porosity to 0.55 total porosity. The lower percolation threshold appears at approximately 10% total porosity and the upper percolation threshold appears at approximately 20% total porosity. Results obtained by creating fifteen solid electrode models, with total porosity as shown, and sectioning with 100 images, each image having a resolution of 100 by 100 pixels. The corresponding data points for example electrodes I through V are also shown for reference. Electrode I Electrode III Electrode II Electrode IV Electrode V   92 5.2.2 Comparison to Predictions using Models from the Literature A three phase boundary length model based upon available literature [28, 34, 35] is used as a comparison to the connected three phase boundary length calculated using methods presented in this work. Equation (11) is used to calculate the connected three phase boundary length. Equation (22) is used to estimate the average coordination number, and the percolation probabilities, Pe and Pi, are calculated according to equation (20). It should be noted that for monosized particles, eeZ −  is equal to Zne . Equations (11), (20) and (22) are rewritten below for clarity of description. The number of particles on the lattice, n, is taken to be the total number of lattice sites less the number of particles removed from the lattice to introduce porosity. All particles are taken to have the mean particle diameter because equation (22) is only valid for monosized spheres. The contact angle between two opposite phased particles is taken to be the reference contact angle used in this work of pi/12 radians. Table 5-2  compares the connected three phase boundary length calculated using equation (11) and the connected three phase boundary length calculated in this work. ( )( ) ieietpbc PPfZnfl=λ  (11) 4.05.2 2 41               − −= −ee e Z P  (20) ( ) 82.061.1 48.1 ≤= − εεZ  (22)    93  Example No. Predicted Connected Three Phase Boundary Length - Eq.(1-11) Connected Three Phase Boundary Length (this work) Percent Difference Compared to This Work I 15x1010m/m3 14x1010m/m3 9.3% II 4.7x1010m/m3 2.5x1010m/m3 91% III 0 0.18x1010m/m3 -100% IV 0 0.06x1010m/m3 -100% V 5.6x1010m/m3 1.8x1010m/m3 210%  For examples I and II, the connected three phase boundary length calculated using equation (11)  differs from the connected three phase boundary length calculated in this work by approximately 10% and a factor of two, respectively. The porosity of both electrodes is fully open, and therefore the difference is likely due to the more realistic treatment of the geometry in this work compared to the monosized random packing of spherical particles assumed in equation (11). Equation (11) does not consider interference between particle contacts, which causes the three phase boundary length to be overestimated. The 9.3% difference between the predicted three phase boundary lengths from equation (11) and this work is not likely to be indicative of the predictive capability of equation (11). This is because the connectivity coefficient of example electrode I is 0.52, which is lower than expected for three percolated phases within the electrode. Conclusive comparison between the two methods for example electrode I is therefore not possible until a dispersion study to determine the required lattice size based on the connected three phase boundary length is conducted. A further shortcoming of equation (11) is that percolation of the porosity phase is not accounted for, but rather, the porosity is implied to be fully open. The percolation probabilities have only been developed for porosities between the approximate range Table 5-2: Comparison of connected three phase boundary length calculated in this work to connected three phase boundary length using combined model based on available literature.   94 between 0.36 and 0.43. Therefore, examples II and V fall outside the range of applicability of the models available in the literature. However, these two electrodes may have merit as a feasible composite layer within a stepwise graded electrode, near the electrolyte, provided that the effective diffusion coefficient is sufficiently large so that mass transfer does not significantly limit the output current at the desired operating point. Examples III and IV are predicted to have zero connected three phase boundary length because the coordination numbers, eeZ −  and iiZ −  of the electronic and ionic conducting phases, respectively, are below two. Therefore, the percolation probabilities in an infinite lattice are zero. Electrode III is not expected to be feasible for any layer within a solid oxide fuel cell electrode. This is because of the absence of reaction sites (connected three phase boundaries) near the electrolyte. Example electrode IV, however, is commonly used as a single layer electrode, and reasonable performances can be achieved at high temperature, despite the reactions being confined only to the interface with the electrolyte. 5.3 Effective Diffusivity in Sintered Composite Electrodes 5.3.1 Example Electrode Models The method for calculation of the effective diffusion coefficient was described in Chapter Four using example electrode IV for illustration purposes. The open porosity of example electrode IV is 35%, which is consistent with most single layer composite solid oxide fuel cell electrodes presented in the literature, which range from approximately 30% to 40% porosity. In the literature, data on the tortuosity factor in solid oxide fuel cell electrodes for porosities less than 25% is non-existent, and therefore, example electrodes II and V, from Figure 5-1, were evaluated because the open porosities are 0.19 and 0.07, respectively. Table 5-3 summarizes the effective diffusion coefficients for these three example electrodes.    95  Example No. Open Porosity Edge Length of Finite Element Effective Diffusion Coefficient Tortuosity Factor Bulk Diffusion Coefficient Correction Factor†† IV 0.35 1.85µm 2.59x10-5m2/s 2.7 0.13 II 0.19 1.25µm 4.21x10-6m2/s 9.1 0.02 V† 0.07 1.11µm 1.45x10-7m2/s 100 7x10-4 †Mesh created based on image resolution of 100 by 100 pixels, rather than 76 by 76 pixels as used for example electrode geometries II and IV. †† Ratio of porosity to tortuosity factor. The edge length of the finite elements in each of the three electrode geometries is greater than 1µm. The discretization of the geometry using images with a resolution of 76 by 76 pixels does not precisely reconstruct the modeled geometry.  A pore that is smaller than 1µm may not be assigned to the pore phase in the image if the width of the pore is less than the width of a pixel. The contribution of this discretization error to the effective diffusion coefficient is assumed to be negligible based on the mesh refinement study in Chapter Four. The mean free path of oxygen and nitrogen was estimated to be approximately 0.4µm in Chapter One, and therefore, the ordinary diffusion mechanism is considered to dominate the diffusive transport of oxygen within the pore space. Example electrode IV was calculated to have a tortuosity factor of 2.7, and the effective diffusion coefficient of 2.59x10-5m2/s is 13% of the bulk diffusion coefficient. The results in Table 5-3 indicate that the tortuosity factor increases, and thus the effective diffusion coefficient decreases, as the open porosity decreases. Figure 5-7 shows the mesh nodes and location of connected three phase boundaries for example electrode II. The distribution of connected three phase boundaries appears to be uniform in that there do not appear to be large regions of pore space without a connected three phase boundary. This is due to the Table 5-3: Comparison of effective diffusion coefficient for three example electrodes. The bulk binary diffusion coefficient is 2.01x10-4m2/s for the O2-N2 system at 1123K and 1atm. Fuel cell operating point is 0.7V.   96 roughly equal solid volume fractions of the electronic and ionic conducting particles, equal to 0.44 and 0.56, respectively.  The superficial diffusive flux decreases linearly from a maximum value at the gas channel to approximately zero at the electrolyte interface, as shown in Figure 5-8(a). The cross section area of the pore space remains approximately constant through the entire thickness of the electrode, as shown in Figure 5-8(b), and so the decrease in superficial diffusive flux is due to the consumption of oxygen within the bulk of the electrode. Figure 5-7: Finite element mesh nodes, shown as blue dots, and locations of connected three phase boundaries, shown as red dots, for example electrode II. The gas channel is the xy plane at z equal to zero. The electrolyte interface is the xy plane at the maximum z coordinate. z x   97  To investigate the impact of fuel cell operating point (i.e. rate of electrochemical reaction) on the apparent effective diffusivity of composite electrodes, example electrode II was subsequently modeled by applying the Neumann boundary conditions not at the locations of the connected three phase boundary voxels, but rather, only in the xy plane at the maximum z coordinate, as shown in Figure 5-9. The superficial diffusive flux in this case is constant throughout the electrode thickness and equal to 1.8x10-6mol/m2/s, compared to 1.5x10- 6 mol/m2/s for the Neumann boundary conditions applied at the locations shown in Figure 5-7. The tortuosity factor for the case where the boundary conditions were applied as shown in Figure 5-9 was 8.6, compared to 9.1 for the application of Neumann boundary conditions according to Figure 5-7. This result suggests that when the porosity is fully percolated, reactions occurring within the bulk of the electrode, rather than only at the interface with the electrolyte, do not significantly affect the effective diffusion coefficient.   Figure 5-8: (a) Superficial oxygen diffusive flux for example electrode II. The value of the superficial flux decreases linearly from the gas channel interface to the electrolyte interface. (b) Cross section area of the pore space in the xy plane through the thickness of the electrode, z. Position [10-4 m] (a) (b) Su pe rfic ial Dif fus ive  Fl ux [10 -5  m ol/m 2 /s ] Position [10-4 m] Po re Cro ss Se ctio n A rea  [1 03 µm 2 ] 2 5 6 4 3 1 1.0 1.5 2.0 2.5 3.0 3.5 0.50 0.75 1.00 1.25 0.25 0.20 0.40 0.60 0.80 1.00 1.20 1.40 0.5 0 0   98  Example electrode V does not have fully percolated porosity. The nodes and locations of the Neumann boundary conditions are shown in Figure 5-10(a), and the corresponding finite element mesh is shown in Figure 5-10(b). The effective diffusion coefficient is calculated to be 1.45x10-7m2/s, which results in a tortuosity factor of 100. This is much greater than both example electrode II, at the upper porosity percolation threshold, and example electrode IV, in the common porosity range for single layer composite electrodes. Figure 5-9: Nodes, shown as blue dots, and location of Neumann boundary conditions, shown as red dots, for example electrode II. The gas channel is the xy plane at z equal to zero. The electrolyte interface is the xy plane at the maximum z coordinate. z x   99  The superficial diffusive flux decreases to approximately zero at approximately 85µm into the electrode from the gas channel, as shown in Figure 5-11. Pores whose axes deviate greatly from the z direction do not contribute greatly to the superficial diffusive flux in the z direction. It can be seen in Figure 5-12 that the concentration does not decrease evenly throughout the electrode, indicating that the direction of the interstitial flux (local flux within the open porosity) deviates significantly from the major diffusion direction, with the average deviation increasing toward increasing z as a result of pores that terminate within the bulk of the electrode. Figure 5-10: Mesh nodes, shown as blue dots, and location of Neumann boundary conditions, shown as red dots, for example electrode V. (b) Corresponding finite element mesh. The element edges are not shown. (a) (b) z x   100    Figure 5-11: Mean superficial diffusive flux for example electrode V. Figure 5-12: Average concentration in xy plane through thickness of example electrode V. Inset: zoom in on concentration profile showing difference of 1.6x10-3 mol/m3 in the concentration between the gas channel and the cathode-electrolyte interface. 1.6x10-3 mol/m3 113µm   101 For comparison, example electrode V was also modeled by specifying Neumann boundary conditions only in the xy plane at the maximum z coordinate, as shown in Figure 5-13. In this case, the effective diffusion coefficient was lower, having a value of 7.58x10-8m2/s and a tortuosity factor of 186. This indicates that the consumption of oxygen within the electrode does affect the apparent effective diffusion coefficient when the porosity is not fully percolated. The reason for this is that there are pores that do not extend the full length from the gas channel to the electrolyte. Therefore, the open porosity is effectively limited to the value of 7% calculated from the full geometry, whereas when the boundary conditions are applied at all reaction sites, the pores that terminate before the electrolyte can still contribute to the current density. At steady state, if there is no connected three phase boundary at a pore that terminates within the bulk of the electrode, there is no concentration gradient and thus no net diffusive flux there.  To further illustrate the effect of dead ended pores, example electrode V was solved with the Neumann boundary conditions applied at the locations of the connected three phase boundary voxels but with an operating point of 0.3V, compared to the 0.7V used for all other analyses described above. The tortuosity factor was calculated to be 2, compared to a value Figure 5-13: Mesh nodes, shown as blue dots, and locations of Neumann boundary conditions, shown as red dots, for alternate analysis of example electrode V. The Neumann boundary conditions are only applied in the xy plane at the maximum z coordinate and in the z direction. z x   102 of 100 at an operating point of 0.7V. Operating at 0.3V increases the overpotential, η, defined as the difference between the operating potential and the equilibrium potential, and therefore the Faradaic current density increases, as implied by equation (55) rewritten below. The effective diffusion coefficient is still low, equal to 7.06x10-6m2/s, due to the low porosity of 0.07. However, it is clear that there is a dependence of the apparent effective diffusion coefficient on operating point for composite electrodes with porosity that is not fully percolated. The dead end pores are more active at higher current operating points because the current is more equally distributed throughout the electrode. Thus, the electrode becomes more fully utilized at higher operating currents, allowing dead-ended pores to contribute to the net diffusive flux at high current densitites, but not at lower current densities. ( ) ( ) ( )               −      −= ηαηαη RT F c c RT Fii cc O O caoF expexp 0 2 2  (55)  The tortuosity factor for example electrode II was calculated at an operating point of 0.3V to be 8.9, which only differs by 2% compared to the value of 9.1, calculated at an operating point of 0.7V. This result further supports the observation that the apparent effective diffusion coefficient changes at different fuel cell operating points in electrodes with porosity between the upper and lower percolation thresholds. The effective diffusion coefficient does not change not significantly in electrodes with porosity above the upper percolation threshold. 5.3.2 Comparison of Calculated Tortuosity Factor to Measured Values Available in the Literature. The tortuosity factor calculated for example electrodes II, IV, and V at an operating point of 0.7V are plotted in Figure 5-14 together with measured data from the literature [20]. This comparison is intended to show that the tortuosity factor for example electrode IV, which has porosity in the range of the measured values, is comparable to that measured in a real solid oxide fuel cell anode. The comparison is not intended to provide validation of the model for two reasons. First, the method used to extract the tortuosity from the measured effective diffusion coefficients and the definition of tortuosity [factor] used by the authors was not   103 presented in [20].  Second, the microstructure of the anode used in the measurements is not described in detail, so it is unclear to what extent the Knudsen diffusion and ordinary diffusion mechanisms contribute to the effective diffusion coefficient. The authors only mention that Knudsen diffusion was considered. The comparison in Figure 5-14 also illustrates the increase in tortuosity factor with decreasing porosity, which is especially apparent at porosities below the upper percolation threshold of approximately 20%.  The method developed to estimate the connected three phase boundary length is applicable to all particle size scales that may be modeled using the methods of Chapter Two, from the nano to the micro scale in solid oxide fuel cells. It is apparent from example electrode I that a study of the dispersion in connected three phase boundary length is required to ensure repeatable results. The method developed to estimate the effective diffusion coefficient is Figure 5-14: Comparison of tortuosity factor calculated in this work to the tortuosity factor measured by Williford et al. using three different measurement techniques: Mercury porosimetry and the Wike-Kallenbach method using CO2 in N2 and also H2 in N2 [20]   104 limited to electrode geometries for which the ordinary diffusion mechanism dominates, such as the diffusion layers or support layers in a graded solid oxide fuel cell electrode. Diffusion modeling of functional layers, in which the electrochemical reactions occur, requires consideration of the Knudsen diffusion mechanism if the pore sizes result in a Knudsen number near unity. However, despite this limitation, it has been shown that estimation of the effective diffusion coefficient in composite electrodes with porosity below the upper percolation threshold, a condition which may have merit in a layer immediately adjacent to the electrolyte, requires consideration of the location of the reaction sites and thickness of the layer. Otherwise, the species flux within pores that terminate before the interface of the electrode with the electrolyte will not be considered, and the incorrect concentration gradient and superficial diffusive flux would be obtained.    105 Chapter 6 Conclusions and Future Work 6.1 Summary A method for which composition and microstructure dependent variables may be evaluated over a feasible electrode design space has been developed. Specifically, the connected three phase boundary length and effective diffusivity may be calculated for modeled sintered composite geometries comprising a single layer, or for stepwise graded solid oxide fuel cell electrodes when combined with an electrochemical performance model. Evaluation of the three phase boundary length and effective diffusivity required detailed electrode geometry modeling. Both the three phase boundary length and the effective diffusivity are dependent on the porosity and composition of the electrode as well as the diameter of the sintered particles within the electrode. The requirement to model a wide range of electrode designs for characterization of the effect of composition and microstructure on their performance, to enable rational design of graded electrodes, precludes the use of any model currently available in the literature. This is the first work to provide a detailed treatment of the sinter necks between neighbouring sintered particles. The random nature of a porous composite solid oxide fuel cell electrode has been approximated by assigning the diameter and phase of each particle at random, as well as randomly removing particles from a close packed lattice to introduce the desired porosity. The use of a face centered cubic lattice to specify the particle positions restricts the range between the minimum and maximum particle diameter that may be used for a given lattice; however, this geometry creation method serves as a starting point to develop the methods required to calculate the three phase boundary length and effective diffusivity of sintered composite electrodes. The connected three phase boundary length was calculated based on the fraction of the total three phase boundary length in the continuous geometry that is connected. The fraction of connected three phase boundary length was estimated by discretizing the modeled geometry using tomographic methods and determining the fraction of three phase boundary   106 voxels that are connected. This method is suitable for modeled deterministic electrodes where the total three phase boundary length is known a priori; however, this method of calculating the connected three phase boundary length is not applicable to the digital reconstruction of real solid oxide fuel cell electrodes. Despite this limitation, the method for calculating the connected three phase boundary length presented in this work is sufficient to characterize a broad range of electrode designs and develop insight into how the three phase boundary length changes with composition and microstructure. The effective ordinary diffusion coefficient was calculated by solving Laplace’s equation for the local concentration of the diffusing species within the open porosity of the modeled electrode. Using Fick’s law, the effective diffusion coefficient is the ratio of the average superficial diffusive flux to the concentration gradient, in the direction normal to the plane of the electrode. The geometric details of the modeled electrode are therefore contained within the effective coefficient, which may then be used in a macroscopic electrode performance model. The Knudsen diffusion mechanism was not included, and therefore, a constraint is imposed on the size of the particles within the electrode such that the interstitial pore sizes must be greater than the mean free path of the diffusing molecules, generally on the order of 1µm in solid oxide fuel cells. Relevant to modeling graded electrodes, the method developed in this work is suitable for diffusion layers where the ordinary diffusion mechanism dominates gas transport. Despite the limited number of electrode designs for which the three phase boundary length and effective diffusion coefficient have been calculated in this work, some insight has been obtained into the relationship between electrode composition and microstructure and electrode performance. Most significantly, electrode layers having porosity between the lower and upper percolation threshold have different apparent effective diffusion coefficients at different fuel cell operating points. As the demand for current increases, more of the connected three phase boundary sites become active, and therefore a greater fraction of the electrode layer is utilized, resulting in a higher apparent effective diffusivity compared to the same electrode geometry operating at a lower current (rate of reaction).   107 6.2 Recommendations for Future Work The most important studies to be completed in the near future are the dispersion studies for the connected three phase boundary length and the open porosity. Specifically, the size of the lattice in the modeled electrode used in this work was chosen, as a starting point, based on convergence of the total three phase boundary length to within ±5% of the mean value over ten different geometries with the same input parameters. The fraction of connected three phase boundaries may have a large degree of scatter between different modeled electrodes created using the same input parameters if the number of lattice sites is not large enough. Similarly, a modeled electrode with porosity between the upper and lower percolation thresholds may not yield the same open porosity as a different modeled electrode with the same input parameters. This work is currently in progress. The geometry modeling can be improved by removing the constraint imposed by the face centered cubic lattice on the particle positions. In this way, particles of predetermined size may assume any stable position available within a container of specified shape. Once the particle positions are determined and the phase is assigned, the sinter neck geometry can be modeled and resulting total three phase boundary length can be calculated. The tomographic method to calculate the connected three phase boundary length and effective diffusion coefficient can be applied in this situation without modification. The exclusion of the Knudsen diffusion mechanism limits the applicability of the method for calculating the effective diffusion coefficient presented in this work. However, cell performance can still be modeled at fuel cell operating points in the region where mass transport losses are not significant. In order to model the performance of graded electrode layers having a Knudsen number less than or approximately equal to one, the Knudsen diffusion mechanism or both Knudsen and Ordinary diffusion mechanisms, respectively, must be considered. The application of the tomographic methods to calculate the three phase boundary length of real solid oxide fuel cell electrodes requires the ability to estimate the length of the connected three phase boundaries from the individual three phase boundary voxels. Use of a characteristic length for each voxel, depending on the local voxel arrangement, would allow   108 the connected three phase boundary length to be estimated in any electrode that can be fabricated. Furthermore, the modeled geometries created using methods developed in this work could be used to determine the minimum spacing required between successive focused ion beam milling operations that is required to converge to the actual three phase boundary length. The ultimate objective of this work is to apply the methods developed for determining connected three phase boundary length and effective diffusivity to optimize performance of solid oxide fuel cells having functionally graded electrode properties. Once individual composite electrode layers are characterized using the methods presented in this work, the detailed electrode structure is contained in the calculated three phase boundary lengths and the effective diffusion coefficients. Therefore, the individual layers may be systematically combined to design stepwise graded electrode structures using a computer experimental design. The performance of the graded electrode can then be modeled in a macroscopic performance model, and the response variable, for example, polarization resistance, can then be determined for each experiment. If a reliable meta-model for the response variable can be obtained, then the resulting function can be minimized to determine the optimum electrode design scenarios.    109 References [1] Stambouli AB, Traversa E. Solid oxide fuel cells (SOFCs): a review of an environmentally clean and efficient source of energy. Renewable & Sustainable Energy Reviews 2002;6:433. [2] Hassmann K. SOFC Power Plants, the Siemens-Westinghouse Approach. Fuel Cells 2001;1:78. [3] Singhal SC, Kendall K. Electrolytes. High Temperature Solid Oxide Fuel Cells: Fundamentals, Design and applications. Elsevier, 2003. [4] Kenney B, Karan K. Mathematical micromodel of an SOFC composite cathode. Hydrogen and Fuel Cells 2004, 2004. p.1. [5] Newman JS, Tobias CW. Theoretical Analysis of Current Distribution in Porous Electrodes. Journal of the Electrochemical Society 1962;109:1183. [6] Liu Y, Compson C, Liu ML. Nanostructured and functionally graded cathodes for intermediate temperature solid oxide fuel cells. Journal of Power Sources 2004;138:194. [7] Barthel K, Rambert S, Siegmann S. Microstructure and polarization resistance of thermally sprayed composite cathodes for solid oxide fuel cell use. Journal of Thermal Spray Technology 2000;9:343. [8] Cassidy M, Bagger C, Brandon N, Day M. Improved cathode performance using graded structures. Fourth European Solid Oxide Fuel Cell Forum. Proceedings. Lucern, Switzerland: European Fuel Cell Forum, 2000. [9] Holtappels P, Bagger C. Fabrication and performance of advanced multi-layer SOFC cathodes. Journal of the European Ceramic Society 2002;22:41. [10] Hart NT, Brandon NP, Day MJ, Shemilt JE. Functionally graded cathodes for solid oxide fuel cells. Journal of Materials Science 2001;36:1077. [11] Xia CR, Rauch W, Wellborn W, Liu ML. Functionally graded cathodes for honeycomb solid oxide fuel cells. Electrochemical and Solid State Letters 2002;5:A217. [12] Zha SW, Zhang YL, Liu ML. Functionally graded cathodes fabricated by sol-gel/slurry coating for honeycomb SOFCs. Solid State Ionics 2005;176:25. [13] Xu XY, Xia CR, Mao GL, Peng D. Fabrication and perfonnance of functionally graded cathodes for IT-SOFCs based on doped ceria electrolytes. Solid State Ionics 2005;176:1513. [14] Holtappels P, Sorof C, Verbraeken MC, Rambert S, Vogt U. Preparation of porosity-graded SOFC anode substrates. Fuel Cells 2006;6:113. [15] Hart NT, Brandon NP, Day MJ, Lapena-Rey N. Functionally graded composite cathodes for solid oxide fuel cells. Journal of Power Sources 2002;106:42. [16] Gerk C, Willert-Porada M. Development of Graded Composite Electrodes for the SOFC. Materials science forum 1999;308-311:806. [17] Gazzarri JI, Kesler O. Electrochemical AC impedance model of a solid oxide fuel cell and its application to diagnosis of multiple degradation modes. Journal of Power Sources 2007;167:100. [18] Deseure J, Bultel Y, Dessemond L, Siebert E. Theoretical optimisation of a SOFC composite cathode. Electrochimica Acta 2005;50:2037.   110 [19] Cadle PJ, Satterfi.Cn. Uniformity of Diffusivity in a Nickel Base Steam-Hydrocarbon Reforming Catalyst. Industrial & Engineering Chemistry Fundamentals 1968;7:189. [20] Williford RE. Diffusion limitations in the porous anodes of SOFCs. Journal of the Electrochemical Society 2003;150:A1067. [21] Kenney B, Karan K. Impact of nonuniform potential in SOFC composite cathodes on the determination of electrochemical kinetic parameters - A numerical analysis. Journal of the Electrochemical Society 2006;153:A1172. [22] Ioselevich A, Kornyshev AA, Lehnert W. Statistical geometry of reaction space in porous cermet anodes based on ion-conducting electrolytes - Patterns of degradation. Solid State Ionics 1999;124:221. [23] Krishna R, Wesselingh JA. Review article number 50 - The Maxwell-Stefan approach to mass transfer. Chemical Engineering Science 1997;52:861. [24] Mason EA, Malinauskas AP. Gas Transport in Porous Media: The Dusty-Gas Model. Amsterdam: Elesvier Science Publishers B.V., 1983. [25] Epstein N. On Tortuosity and the Tortuosity Factor in Flow and Diffusion Through Porous- Media. Chemical Engineering Science 1989;44:777. [26] Radhakrishnan R, Virkar AV, Singhal SC. Estimation of charge-transfer resistivity of La0.8Sr0.2MnO3 cathode on Y0.16Zr0.84O2 electrolyte using patterned electrodes. Journal of the Electrochemical Society 2005;152:A210. [27] vanHeuveln FH. Electrode properties of Sr-doped LaMnO3 on yttria-stabilized zirconia .1. Three-phase boundary area. Journal of the Electrochemical Society 1997;144:126. [28] Sunde S. Monte Carlo simulations of polarization resistance of composite electrodes for solid oxide fuel cells. Journal of the Electrochemical Society 1996;143:1930. [29] Costamagna P, Costa P, Antonucci V. Micro-modelling of solid oxide fuel cell electrodes. Electrochimica Acta 1998;43:375. [30] Bouvard D, Lange FF. Relation between percolation and particle coordination in binary powder mixtures. Acta Metallurgica Et Materialia 1991;39:3083. [31] Chan SH, Xia ZT. Anode micro model of solid oxide fuel cell. Journal of the Electrochemical Society 2001;148:A388. [32] Kuo CH, Gupta PK. Rigidity and Conductivity Percolation Thresholds in Particulate Composites. Acta Metallurgica Et Materialia 1995;43:397. [33] Mori H, Nonaka N, Mizuno M, Abe H, Naito M. Development of a Geometrical Model for Optimizing Porous Anode Microstructure of Solid Oxide Fuel Cells. Journal of Chemical Engineering of Japan 2008;41:246. [34] Nakagaki M, Sunada H. Theoretical studies on structures of the sedimentation bed of spherical particles. Yakugaku Zasshi 1968;88:651. [35] Suzuki M, Makino K, Yamada M, Iinoya K. A study of coordination-number in a random packed system of monosized sphere particles. Kagaku Kogaku Ronbunshu 1980;6:59. [36] Martinez AS, Brouwer J. Percolation modeling investigation of TPB formation in a solid oxide fuel cell electrode-electrolyte interface. Electrochimica Acta 2008;53:3597.   111 [37] Wilson JR, Kobsiriphat W, Mendoza R, Chen HY, Hiller JM, Miller DJ, Thornton K, Voorhees PW, Adler SB, Barnett SA. Three-dimensional reconstruction of a solid-oxide fuel-cell anode. Nature Materials 2006;5:541. [38] Hines AL, Maddox RN. Mass Transfer: Fundamentals and Applications: Prentice Hall, 1985. [39] Bird GA. Definition of Mean Free-Path for Real Gases. Physics of Fluids 1983;26:3222. [40] Kenney B, Karan K. Engineering of microstructure and design of a planar porous composite SOFC cathode: A numerical analysis. Solid State Ionics 2007;178:297. [41] Chan SH, Khor KA, Xia ZT. A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness. Journal of Power Sources 2001;93:130. [42] Zhu HY, Kee RJ. A general mathematical model for analyzing the performance of fuel-cell membrane-electrode assemblies. Journal of Power Sources 2003;117:61. [43] Yakabe H, Hishinuma M, Uratani M, Matsuzaki Y, Yasuda I. Evaluation and modeling of performance of anode-supported solid oxide fuel cell. Journal of Power Sources 2000;86:423. [44] Bove R, Ubertini S. Modeling solid oxide fuel cell operation: Approaches, techniques and results. Journal of Power Sources 2006;159:543. [45] Stauffer D. Introduction to Percolation Theory. Philadelphia: Taylor & Francis Inc., 1985. [46] Nandakumar K, Shu YQ, Chuang KT. Predicting geometrical properties of random packed beds from computer simulation. Aiche Journal 1999;45:2286. [47] Hoshen J, Kopelman R. Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Physical Review B 1976;14:3438. [48] Poling BE, Prausnitz JP, O'Connell JP. The Properties of Liquids and Gases. New York: McGraw-Hill, 2000. [49] Gazzarri JI, Kesler O. Non-destructive delamination detection in solid oxide fuel cells. Journal of Power Sources 2007;167:430. [50] Dees DW, Claar TD, Easler TE, Fee DC, Mrazek FC. Conductivity of porous Ni/ZrO2-Y2O3 cermets. Journal of the Electrochemical Society 1987;134:2141. [51] Gazzarri J. Impedance Model of a Solid Oxide Fuel Cell for Degradation Diagnosis. Mechanical Engineering, Doctor of Philosophy. Vancouver: University of British Columbia, 2007.    112 Appendix A Computer Code for Electrode Geometry Modeling Matlab Code for Numerical Geometry Creation The first step in creating the electrode solid model is to specify the position, size, and phase of each particle on the lattice. The lattice is defined by the reference contact angle Theta0, between particles and the reference particle diameter, D0. The radius of curvature of the sinter neck, rs, is the same for each particle intersection. The porosity is obtained by specifying the fraction of particles to remove from the lattice. The solid volume fraction of the electronic conductor, Ecomp, is a user input, as are the fraction of particles to remove, fracremove, and the reference particle diameter and contact angle. Theta0 = pi/12;   %reference contact angle (set to pi/12) D0 = <input>;   %Reference particle diameter (user input) rs = 0.1*(D0/2)*ones(1,12); %Sinter neck curvature radius fracRemove = <input>;  %Fraction of particles to remove (user input) Ecomp = <input>;   %Solid vol. fraction of e-conductor (user input) aConst = 12;   %Number of fcc unit cells on geometry edge  The coordinates of the particle centers are determined by the function setCoordinates.m and is called using the following command, where Pos is a matrix having three columns (representing the x, y and z coordinates), and 7813 rows (the lattice size). Subt is an array containing the row numbers in Pos to remove to introduce porosity. a is the lattice constant, and Numsites is the number of lattice sites, equal to 7813. [Pos subt a NumSites] = setCoordinates(aConst, Theta0, D0, fracRemove);  The function setCoordinates.m is as follows: %************************************************************************** function [Pos,subt,a,NumSites]=setCoordinates(aConst,Theta0,D0,fracremove) a = sqrt(2)*D0*cos(Theta0);     % Lattice constant  % Center position of particles on fcc lattice, Type 1 n = 0;     %counter for z = 0:a:aConst*a     for y = 0:a:aConst*a         for x = 0:a:aConst*a            n = n + 1;            Pos(n,:) = [x y z];         end     end end % N1full = n;   %number of type one particles   113  % Center position of particles on fcc lattice, Type 2 for z = 0:a:aConst*a     for y = a/2:a:(2*aConst-1)*a/2         for x = a/2:a:(2*aConst-1)*a/2            n = n + 1;            Pos(n,:) = [x y z];         end     end end % N2full = n - N1full;  %number of type two particles  % Center position of particles on fcc lattice, Type 3 for z = a/2:a:(2*aConst-1)*a/2     for y = a/2:a:(2*aConst-1)*a/2         for x = 0:a:aConst*a            n = n + 1;            Pos(n,:) = [x y z];         end     end end % N3full = n - (N2full + N1full); %number of type three particles  % Center position of particles, on fcc lattice Type 4 for z = a/2:a:(2*aConst-1)*a/2     for y = 0:a:aConst*a         for x = a/2:a:(2*aConst-1)*a/2            n = n + 1;            Pos(n,:) = [x y z];         end     end end % N4full = n - (N3full + N2full + N1full); %number of type four particles  NumSites = n;      % number of lattice sites removeparticles = round(fracremove*NumSites); %integer # of particles to remove subtRandom = randperm(NumSites);  %Random permutation of 1:7813 subt = subtRandom(1:removeparticles);  %Particles to be removed %**************************************************************************  Next, the particles identified to be removed are removed: Pos(subt',:) = [];    %Remove particles with index same as array subt N = length(Pos);  %Number of particles remaining on the lattice  The particle diameters are set using setDiameters.m, where D is an array containing the diameter of each particle with center coordinate listed in matrix Pos. The array index in D corresponds to the row number in Pos. The particle dimeters have four standard deviations between the minimum diameter, Dmin, and mean diameter, Dmean, and between the mean diameter and the maximum diameter, Dmax. [D stdeviation] = setDiameters(Theta0, NumSites, D0);  The function setDiameters.m is as follows:   114 %**************************************************************************  function [D stdeviation] = setDiameters(Theta0, NumSites, D0) Dmin = D0*cos(Theta0);         % Min diameter Dmax = sqrt(2)*Dmin;           % Max diameter Dmean = 0.5*(Dmin+Dmax);       % Mean diameter between Dmin and Dmax D = Dmean + sqrt((1/4*(Dmean - Dmin))^2)*randn(1,NumSites); %normal distribution stdeviation = sqrt((1/4*(Dmean - Dmin))^2); %four std deviations: Dmin->Dmean %**************************************************************************  The indices in D corresponding to the particle diameters for particles that have been removed from the lattice must be removed: D(subt) = [];       % also remove the diameter index  The solid phase of each particle is assigned using the function setPhase.m as follows. sPhase = setPhase(N, Ecomp);  The function setPhase.m is as follows. An array with a random permutation of the integers from 1 to the number of particles on the lattice, N, is created. Then the electronic conductor solid volume fraction is used to determine how many particles (elements of the array sPhase) contain a 1, corresponding to the electronic conducting phase. %************************************************************************** function sPhase = setPhase(N, Ecomp) indexRandom = randperm(N);  %random permutation of the integers 1:N sPhase = zeros(1,N);   %zero represents ionic conductor. %reassign the phase of N*Ecomp particles to be the electronic conducting %phase. This phase is marked by a value of 1. for ii = 1:round(N*Ecomp)     sPhase(indexRandom(ii)) = 1; end %**************************************************************************  A new matrix, matrx, is the transpose of the particle position matrix Pos. matrx = [Pos(:,1)'; Pos(:,2)'; Pos(:,3)'];  Now the properties of the neighbour particles must be determined using the function Neighbors.m as follows. PN is a matrix with 36 columns corresponding to the x, y, and z coordinates for all 12 possible neighbors for a given particle on the lattice. The first three columns are for neighbour particle position 1, columns 4 through 6 are for neighbor particle position 2 etc., the number of rows is equal to the number of particles on the lattice. A -1 in the column represents that the particle is not present on the lattice. DN is a matrix with 12 columns, containing the diameter of all twelve neighbor particles, where a -1 in a particular column indicates that no neighbor particle is present at that location. The number of rows is   115 also equal to the number of particles on the lattice. PhN is a matrix with 12 columns, each column identifying the phase of each neighbor particle. The number of rows is equal to the number of particles on the lattice. PM = matrx'; DM = D'; [PN, DN, PhN, tol] = Neighbors(PM, DM, N, D0, a, sPhase);  The function Neighbors.m is as follows: %************************************************************************** function [PositionNeighbors, DiameterNeighbors, PhaseNeighbors, tol]...     = Neighbors(PM, DM, N, D0, a, Phase)  % Want index of all neighbor particles. There are a max of 12. The % neighbors are searched in the following order. XY Plane Quandrant 1 then % 2 then 3 then 4, XZ Plane Quadrant 1 then 2 then 3 then 4, YZ Plane % Quadrant 1 then 2 then 3 then 4. The index k ranges from 1 to 12, where the XY plane, quadrant 1 is index 1, XY plane quadrant 2 is index 2 etc…  k = 1;  xK = 3*k - 2;    zK = 3*k;  tol = D0*1e-2;  for i = 1:N     %(+x,+y,0)     if k == 1     Ind = find(PM(:,1)>(PM(i,1) + a/2 - tol) & PM(:,1)<(PM(i,1) + a/2 + tol)&...                PM(:,2)>(PM(i,2) + a/2 - tol) & PM(:,2)<(PM(i,2) + a/2 + tol)&...                PM(:,3)>(PM(i,3) + 0   - tol) & PM(:,3)<(PM(i,3) + 0   + tol));         if isempty(Ind) == 0             PositionNeighbors(i,xK:zK) = PM(Ind,:);             DiameterNeighbors(i,k) = DM(Ind);             PhaseNeighbors(i,k) = Phase(Ind);         elseif isempty(Ind) == 1             PositionNeighbors(i,xK:zK) = -1;             DiameterNeighbors(i,k) = -1;             PhaseNeighbors(i,k) = -1;         end     end      k = 2; xK = 3*k - 2; zK = 3*k;     %(-x,+y,0)     if k == 2     Ind = find(PM(:,1)>(PM(i,1) - a/2 - tol) & PM(:,1)<(PM(i,1) - a/2 + tol)&...                PM(:,2)>(PM(i,2) + a/2 - tol) & PM(:,2)<(PM(i,2) + a/2 + tol)&...                PM(:,3)>(PM(i,3) + 0   - tol) & PM(:,3)<(PM(i,3) + 0   + tol));         if isempty(Ind) == 0             PositionNeighbors(i,xK:zK) = PM(Ind,:);             DiameterNeighbors(i,k) = DM(Ind);             PhaseNeighbors(i,k) = Phase(Ind);         elseif isempty(Ind) == 1             PositionNeighbors(i,xK:zK) = -1;             DiameterNeighbors(i,k) = -1;             PhaseNeighbors(i,k) = -1;         end     end      k = 3; xK = 3*k - 2; zK = 3*k;     %(-x,-y,0)   116     if k == 3     Ind = find(PM(:,1)>(PM(i,1) - a/2 - tol) & PM(:,1)<(PM(i,1) - a/2 + tol)&...                PM(:,2)>(PM(i,2) - a/2 - tol) & PM(:,2)<(PM(i,2) - a/2 + tol)&...                 PM(:,3)>(PM(i,3) + 0   - tol) & PM(:,3)<(PM(i,3) + 0   + tol));         if isempty(Ind) == 0             PositionNeighbors(i,xK:zK) = PM(Ind,:);             DiameterNeighbors(i,k) = DM(Ind);             PhaseNeighbors(i,k) = Phase(Ind);         elseif isempty(Ind) == 1             PositionNeighbors(i,xK:zK) = -1;             DiameterNeighbors(i,k) = -1;             PhaseNeighbors(i,k) = -1;         end     end      k = 4; xK = 3*k - 2; zK = 3*k;     %(+x,-y,0)     if k == 4     Ind = find(PM(:,1)>(PM(i,1) + a/2 - tol) & PM(:,1)<(PM(i,1) + a/2 + tol)&...                PM(:,2)>(PM(i,2) - a/2 - tol) & PM(:,2)<(PM(i,2) - a/2 + tol)&...                PM(:,3)>(PM(i,3) + 0   - tol) & PM(:,3)<(PM(i,3) + 0   + tol));         if isempty(Ind) == 0         PositionNeighbors(i,xK:zK) = PM(Ind,:);         DiameterNeighbors(i,k) = DM(Ind);         PhaseNeighbors(i,k) = Phase(Ind);     elseif isempty(Ind) == 1         PositionNeighbors(i,xK:zK) = -1;         DiameterNeighbors(i,k) = -1;         PhaseNeighbors(i,k) = -1;         end     end  <<< The code for the two remaining planes (YZ and ZX) is NOT included here as it is analogous to the code above for the XY plane >>>> %**************************************************************************  The magnitude of the largest particle coordinate is maxP. maxP = max(max(PM));  Now that the positions, diameters, and phases for all particles and their nearest neighbours are known, the text files containing this information need to be written. An electrode geometry with 7813 lattice sites leads the solid model part file in Unigraphics to be very large and unmanageable. Therefore, the cubic geometry envelope is divided into 27 sub cubes, and data files containing the particle position, size, and phase, as well as the neighbor particle position, size, and phase for each sub cube are written. This code to accomplish this is not included here. Visual Basic.NET code for Solid Model Creation   117 Unigraphics NX5 has an integrated Journaling feature that is capable of reading the geometry creation commands from a program written in the Visual Basic.NET programming language. A program has been written to create a solid model using the output data files from the numerical geometry described above, which consists of the particle positions, sizes, and phases as well as the neighbor particle positions, sizes, and phases. The entire program is not included here as it is over 126 pages (8-½” by 11”) long, single spaced with size 10 Courier New font. The full program is available from the author upon request (tcm@mie.utoronto.ca). The program is explained for creation of a single sub cube, but has been implemented such that 13 different computers simultaneously create two sub cubes each, and a 14th computer creates the 27th sub cube. This is because the geometry creation program for one sub cube takes up to 8 hours to complete. The computers that are part of the PACE program at the University of British Columbia were used to create the geometry described in this Appendix. The following Imports allow access to the NXOpen Application Programming Interface. Imports NXOpen Imports NXOpen.Features Imports NXOpen.Preferences Imports NXOpen.UF Imports NXOpen.UI Imports NXOpen.Utilities  The module is called SinteredSpherePacking and the Unigraphics sessions: Session and UFSession, must be declared. Module SinteredSpherePacking     Dim theSession As Session     Dim theUfSession As UFSession  The number of particles in the sub cubes is entered into the variable NumSpheres. This value is obtained from the output of the numerical geometry creation described above.     Dim NumSpheres() As Integer = {<user input>}     Dim maxCol As Integer = max(NumSpheres())  The solid model part file name is contained in partFileName. Data file names are stored in the variables as follows. The data files contain the particle sizes, positions, and phases for each particle and also for the neighbor particles of each particle on the lattice.     Dim partFileName() As String = {"SubCube1"}    ‘part file name     Dim PNfileName() As String = {"PN1csvlist.dat"} ‘neighbour positions     Dim DNfileName() As String = {"DN1csvlist.dat"} ‘neighbour diameters   118     Dim posFileName() As String = {"M1csvlist.dat"} ‘particle positions     Dim diamFileName As String = "Diamcsvlist.dat" ‘particle diameters     Dim sPhaseFileName As String = "sPhase1csvlist.dat" ‘particle phases  The particle diameters are read, using the following code, into the variable diam.       Dim diam(26, maxCol)       Dim FileReaderDiam As StreamReader       Dim LineInDiam As String       Dim FieldArrayDiam       Dim jDiam As Integer = 0       Dim countDiam As Integer        Dim path2Diam As String = diamFileName       Dim DiamName As String = Path.Combine(path1, path2Diam)       FileReaderDiam = File.OpenText(DiamName)       LineInDiam = FileReaderDiam.ReadLine()              While LineInDiam <> ""                 FieldArrayDiam = Split(LineInDiam, ",")                 If LineInDiam <> "" Then                     For countDiam = 0 To maxCol                         diam(jDiam, CountDiam) = FieldArrayDiam(CountDiam)                     Next                 End If                 jDiam = jDiam + 1                 LineInDiam = FileReaderDiam.ReadLine()             End While             FileReaderDiam.Close()  The particle phases are read using analogous code into the variable sPhase. It should be noted that VisualBasic.NET identifies the first position within an array as position zero; whereas, Matlab identifies the first position within an array (or vector) as one. The neighbor particle diameters are read into variable ND as follows.     Dim ND(NumSpheres, 11)     Dim FileReaderD As StreamReader     Dim LineInD As String     Dim FieldArrayD     Dim jD As Integer = 0     Dim countD As Integer     Dim DNfile As String = DNfileName(numLoops)     Dim path2DN As String = DNfile     Dim DNname As String = Path.Combine(path1, path2DN)     FileReaderD = File.OpenText(DNname)     LineInD = FileReaderD.ReadLine()              While LineInD <> ""                    FieldArrayD = Split(LineInD, ",")                    If LineInD <> "" Then                        For countD = 0 To 11                            ND(jD, CountD) = FieldArrayD(CountD)                        Next                    End If                    jD = jD + 1        LineInD = FileReaderD.ReadLine()       End While   119       FileReaderD.Close()  The neighbour particle positions are read into variable NP using analogous code to that for ND. The x, y, and z coordinates of the particle positions are read into arrays PosX, PosY, and PosZ, respectively as follows. Dim PosX(NumSpheres()) As String Dim PosY(NumSpheres()) As String Dim PosZ(NumSpheres(n)) As String Dim FileReader As StreamReader       Dim LineIn As String       Dim jPos As Integer = 0        Dim MFile As String = posFileName(numLoops)       Dim path2MFile As String = MFile       Dim MFileName As String = Path.Combine(path1, path2MFile)        FileReader = File.OpenText(MFileName) LineIn = FileReader.ReadLine()            While LineIn <> ""                    If LineIn <> "" Then                          If jPos = 0 Then                             PosX = Split(LineIn, ",")                          End If                          If jPos = 1 Then                             PosY = Split(LineIn, ",")                          End If If jPos = 2 Then PosZ = Split(LineIn, ",") End If LineIn = FileReader.ReadLine() End If jPos = jPos + 1 End While FileReader.Close() A new solid model part file is created as follows.      Dim UFPart As NXOpen.Tag      Dim units As Integer = 2      Dim name As String      Dim partFile As String = partFileName(numLoops)      Dim path2Part As String = partFile      Dim partName As String = Path.Combine(path1, path2Part)      theUfSession.Part.[New](partName, units, UFPart)      Dim workPart As Part = theSession.Parts.Work   workPart.ModelingViews.WorkView.RenderingStyle= View.RenderingStyleType.Shaded workPart.Preferences.ShadeVisualization.ShadedViewTolerance = Preferences.PartVisualizationShade.ShadedViewToleranceType.Fine Dim object_in_part As NXOpen.Tag = theSession.Parts.Work.Tag  The following loop executes j times, where j is the number of particles on the lattice. The particles are initially modeled as a sphere. If the jth index of sPhase is a 1, then the particle is coloured black to indicate an electronic conducting particle.   120 Dim j As Integer  For j = 0 To NumSpheres() Dim Position()() As Double = New Double(NumSpheres())() {} Position(j) = New Double(2) {PosX(j), PosY(j), PosZ(j)} Dim r As Double = diam(, j) / 2  'Create sphere and assign color theUfSession.Modl.CreateList(face_list) theUfSession.Modl.CreateSphere1(FeatureSigns.Nullsign, Position(j), diam(, j), sphere_obj_id) Dim sphere_body_tag As NXOpen.Tag theUfSession.Modl.AskFeatBody(sphere_obj_id, sphere_body_tag) Dim sphere_body As Body = CType(NXObjectManager.Get(sphere_body_tag), Body) theUfSession.Modl.AskBodyFaces(sphere_body_tag, face_list)                      If sPhase(, j) = 1 Then                         sphere_body.Color = 216 'Black                         sphere_body.RedisplayObject()                     ElseIf sPhase(, j) = 0 Then                         sphere_body.Color = 130 'Dark Gray                         sphere_body.RedisplayObject()                     End If  For each particle j, the sinter necks for all intersections with its nearest neighbours must be created. The following illustrates the creation of only one sinterneck that would occur if particle j has a neighbour in the XY plane, quadrant one. The remaining sinter necks are created in an analogous way. rs_XY_1 is the sinter neck curvature radius equal to one tenth the reference particle diameter. Dim xK, yK, zK As Double Dim rs_XY_1 As Double = <input> 'neighbor 1 sinter radius in xy plane  The following code creates the revolved sinter neck for particle j in the XY plane, quadrant one. LengthPOrigin is the distance between nearest neighbour particle centers. SplaneInt is the distance from center of particle j to plane of intersection between the two intersection particles. rN is the radius of the neighbour particle. r is the radius of particle j. The shape to be revolved to create the sinter neck is described in Chapter Two. For k = 0 To 11 xK = 3 * k + 0 'x coordinate index in neighbour particle center, x yK = 3 * k + 1 'y coordinate index in neighbour particle center, y zK = 3 * k + 2 'z coordinate index in neighbour particle center, z  If ND(j, k)>=0 Then        Dim rN As Double = ND(j, k) / 2 Dim LengthPOrigin As Double = Math.sqrt((NP(j, xK) - PosX(j)) ^ 2 + (NP(j, yK) - PosY(j)) ^ 2 + (NP(j, zK) - PosZ(j)) ^ 2) Dim SplaneInt As Double = Math.abs(1 / (2 * LengthPOrigin) * (LengthPOrigin ^ 2 - rN ^ 2 + r ^ 2))   121             Dim Lm As Double = Math.abs(SplaneInt * r / (r + rs_XY_1))             Dim Hn As Double = Math.sqrt((r + rs_XY_1) ^ 2 - SplaneInt ^ 2)             Dim Hm As Double = r / (r + rs_XY_1) * Hn             Dim horizN As Double = Math.sqrt(2) / 2 * SplaneInt             Dim vertN As Double = Math.sqrt(2) / 2 * SplaneInt             Dim phi As Double = Math.acos(SplaneInt / (r + rs_XY_1))             Dim thetaPhi As Double = 45 / 180 * Math.pi + phi  Dim sParcx_XY_1 As Double = (PosX(j) + Math.sqrt(2) / 2 * (Lm - Hm)) 'arc start x Dim sParcy_XY_1 As Double = (PosY(j) + Math.sqrt(2) / 2 * (Lm + Hm)) 'arc start y             Dim sParcz_XY_1 As Double = PosZ(j) 'arc start z Dim eParcx_XY_1 As Double = PosX(j) + horizN - Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'arc end x Dim eParcy_XY_1 As Double = PosY(j) + vertN + Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'arc end y             Dim eParcz_XY_1 As Double = PosZ(j) 'arc end z Dim rscx_XY_1 As Double = PosX(j) + Math.cos(thetaPhi) * (r + rs_XY_1) 'sinter radius center point x Dim rscy_XY_1 As Double = PosY(j) + Math.sin(thetaPhi) * (r + rs_XY_1) 'sinter radius center point y              Dim rscz_XY_1 As Double = PosZ(j) 'sinter radius center point z              Dim pOarcx_XY_1 As Double = (1 / 2) * (sParcx_XY_1 + eParcx_XY_1) 'point on arc x              Dim pOarcy_XY_1 As Double = Math.sqrt(rs_XY_1 ^ 2 - (pOarcx_XY_1 - rscx_XY_1) ^ 2) + rscy_XY_1 'point on arc y              Dim pOarcz_XY_1 As Double = PosZ(j) 'point on arc z  Dim ePline2x_XY_1 As Double = sParcx_XY_1 + Math.sqrt(2) / 2 * Hm 'end point line 2 x  Dim ePline2y_XY_1 As Double = sParcy_XY_1 - Math.sqrt(2) / 2 * Hm 'end point line 2 y  Dim ePline2z_XY_1 As Double = PosZ(j) 'end point line 2 z  Dim ePline1x_XY_1 As Double = eParcx_XY_1 + Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'end point line 1 x  Dim ePline1y_XY_1 As Double = eParcy_XY_1 - Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'end point line 1 y  Dim ePline1z_XY_1 As Double = PosZ(j) 'end point line 1 z  Dim mLine1_XY_1 As Double = (ePline1y_XY_1 - eParcy_XY_1) / (ePline1x_XY_1 - eParcx_XY_1) ' slope of line 1  Dim bLine1_XY_1 As Double = eParcy_XY_1 - mLine1_XY_1 * eParcx_XY_1 ' intercept of line 1  Dim pOline1x_XY_1 As Double = 0.5 * (eParcx_XY_1 + ePline1x_XY_1) 'point on line 1 x  Dim pOline1y_XY_1 As Double = mLine1_XY_1 * pOline1x_XY_1 + bLine1_XY_1 'point on line 1 y  Dim pOline1z_XY_1 As Double = PosZ(j) 'point on line 1 z  Dim mLine2_XY_1 As Double = (ePline2y_XY_1 - sParcy_XY_1) / (ePline2x_XY_1 - sParcx_XY_1) 'slope line 2  Dim bLine2_XY_1 As Double = (sParcy_XY_1 - mLine2_XY_1 * sParcx_XY_1) 'intercept line 2  Dim pOline2x_XY_1 As Double = 0.5 * (sParcx_XY_1 + ePline2x_XY_1) 'point on line 2 x  Dim pOline2y_XY_1 As Double = mLine2_XY_1 * pOline2x_XY_1 + bLine2_XY_1 'point on line 2 y  Dim pOline2z_XY_1 As Double = PosZ(j) 'point on line 2 z  Dim mLine3_XY_1 As Double = (ePline1y_XY_1 - ePline2y_XY_1) / (ePline1x_XY_1 - ePline2x_XY_1) 'slope line 3  Dim bLine3_XY_1 As Double = ePline2y_XY_1 - mLine3_XY_1 * ePline2x_XY_1 'intercept line 3   122  Dim pOline3x_XY_1 As Double = 0.5 * (ePline1x_XY_1 + ePline2x_XY_1) 'point on line 3 x  Dim pOline3y_XY_1 As Double = mLine3_XY_1 * pOline3x_XY_1 + bLine3_XY_1 'point on line 3 y  Dim pOline3z_XY_1 As Double = PosZ(j) 'point on line 3 z   Dim sParcx_XY_1_Tag, sParcy_XY_1_Tag, sParcz_XY_1_Tag As NXOpen.Tag 'tags for start point arc coordinates  Dim eParcx_XY_1_Tag, eParcy_XY_1_Tag, eParcz_XY_1_Tag As NXOpen.Tag 'tags for end point arc coordinates  Dim rscax_XY_1_Tag, rscay_XY_1_Tag, rscaz_XY_1_Tag As NXOpen.Tag 'tags for center coordinates of arc  Dim sParc_XY_1_Tag, eParc_XY_1_Tag, cParc_XY_1_Tag As NXOpen.Tag 'tags for arc start, end, center  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcx_XY_1, sParcx_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcy_XY_1, sParcy_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcz_XY_1, sParcz_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, rscx_XY_1, rscax_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, rscy_XY_1, rscay_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, rscz_XY_1, rscaz_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcx_XY_1, eParcx_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcy_XY_1, eParcy_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcz_XY_1, eParcz_XY_1_Tag)  Dim sParc_XY_1_TagArray() As NXOpen.Tag = {sParcx_XY_1_Tag, sParcy_XY_1_Tag, sParcz_XY_1_Tag}  Dim eParc_XY_1_TagArray() As NXOpen.Tag = {eParcx_XY_1_Tag, eParcy_XY_1_Tag, eParcz_XY_1_Tag}  Dim cParc_XY_1_TagArray() As NXOpen.Tag = {rscax_XY_1_Tag, rscay_XY_1_Tag, rscaz_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sParc_XY_1_TagArray, sParc_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, eParc_XY_1_TagArray, eParc_XY_1_Tag)   123  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, cParc_XY_1_TagArray, cParc_XY_1_Tag)  Dim arcPnts_XY_1_TagArray() As NXOpen.Tag = {cParc_XY_1_Tag, sParc_XY_1_Tag, eParc_XY_1_Tag}              Dim arc_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateArcCenter2Pnts(object_in_part, SmartObject.UpdateOption.WithinModeling, arcPnts_XY_1_TagArray, arc_XY_1_Tag)   Dim sPline2x_XY_1_Tag, sPline2y_XY_1_Tag, sPline2z_XY_1_Tag As NXOpen.Tag 'tags for line 2 start point coordinates  Dim ePline2x_XY_1_Tag, ePline2y_XY_1_Tag, ePline2z_XY_1_Tag As NXOpen.Tag 'tags for line 2 end point coordinates  Dim sPline2_XY_1_Tag, ePline2_XY_1_Tag As NXOpen.Tag 'tags for line 2 start and end points  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcx_XY_1, sPline2x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcy_XY_1, sPline2y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcz_XY_1, sPline2z_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2x_XY_1, ePline2x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2y_XY_1, ePline2y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2z_XY_1, ePline2z_XY_1_Tag)  Dim sPline2_XY_1_TagArray() As NXOpen.Tag = {sPline2x_XY_1_Tag, sPline2y_XY_1_Tag, sPline2z_XY_1_Tag}  Dim ePline2_XY_1_TagArray() As NXOpen.Tag = {ePline2x_XY_1_Tag, ePline2y_XY_1_Tag, ePline2z_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sPline2_XY_1_TagArray, sPline2_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2_XY_1_TagArray, ePline2_XY_1_Tag)  Dim line2Pnts_XY_1_TagArray() As NXOpen.Tag = {sPline2_XY_1_Tag, ePline2_XY_1_Tag}  Dim line2_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateLineTwoPoints(object_in_part, SmartObject.UpdateOption.WithinModeling, line2Pnts_XY_1_TagArray, line2_XY_1_Tag)   124   Dim sPline1x_XY_1_Tag, sPline1y_XY_1_Tag, sPline1z_XY_1_Tag As NXOpen.Tag 'tags for line 1 start coordinates  Dim ePline1x_XY_1_Tag, ePline1y_XY_1_Tag, ePline1z_XY_1_Tag As NXOpen.Tag 'tags for line 2 end coordinates  Dim sPline1_XY_1_Tag, ePline1_XY_1_Tag As NXOpen.Tag 'tags for line 1 start and end points  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcx_XY_1, sPline1x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcy_XY_1, sPline1y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcz_XY_1, sPline1z_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1x_XY_1, ePline1x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1y_XY_1, ePline1y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1z_XY_1, ePline1z_XY_1_Tag)  Dim sPline1_XY_1_TagArray() As NXOpen.Tag = {sPline1x_XY_1_Tag, sPline1y_XY_1_Tag, sPline1z_XY_1_Tag}  Dim ePline1_XY_1_TagArray() As NXOpen.Tag = {ePline1x_XY_1_Tag, ePline1y_XY_1_Tag, ePline1z_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sPline1_XY_1_TagArray, sPline1_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1_XY_1_TagArray, ePline1_XY_1_Tag)  Dim line1Pnts_XY_1_TagArray() As NXOpen.Tag = {sPline1_XY_1_Tag, ePline1_XY_1_Tag}              Dim line1_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateLineTwoPoints(object_in_part, SmartObject.UpdateOption.WithinModeling, line1Pnts_XY_1_TagArray, line1_XY_1_Tag)   Dim sPline3x_XY_1_Tag, sPline3y_XY_1_Tag, sPline3z_XY_1_Tag As NXOpen.Tag 'tags for line 3 start coordinates  Dim ePline3x_XY_1_Tag, ePline3y_XY_1_Tag, ePline3z_XY_1_Tag As NXOpen.Tag 'tags for line 3 end coordinates  Dim sPline3_XY_1_Tag, ePline3_XY_1_Tag As NXOpen.Tag 'tags for line 3 start and end points  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1x_XY_1, sPline3x_XY_1_Tag)   125  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1y_XY_1, sPline3y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1z_XY_1, sPline3z_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2x_XY_1, ePline3x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2y_XY_1, ePline3y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2z_XY_1, ePline3z_XY_1_Tag)  Dim sPline3_XY_1_TagArray() As NXOpen.Tag = {sPline3x_XY_1_Tag, sPline3y_XY_1_Tag, sPline3z_XY_1_Tag}  Dim ePline3_XY_1_TagArray() As NXOpen.Tag = {ePline3x_XY_1_Tag, ePline3y_XY_1_Tag, ePline3z_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sPline3_XY_1_TagArray, sPline3_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline3_XY_1_TagArray, ePline3_XY_1_Tag)  Dim line3Pnts_XY_1_TagArray() As NXOpen.Tag = {sPline3_XY_1_Tag, ePline3_XY_1_Tag}              Dim line3_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateLineTwoPoints(object_in_part, SmartObject.UpdateOption.WithinModeling, line3Pnts_XY_1_TagArray, line3_XY_1_Tag)   Dim section1_XY_1 As Section              section1_XY_1 = workPart.Sections.CreateSection()              Dim nullFeatures_Feature_XY_1 As Features.Feature = Nothing              Dim revolveBuilder1_XY_1 As Features.RevolveBuilder  revolveBuilder1_XY_1 = workPart.Features.CreateRevolveBuilder(nullFeatures_Feature_X Y_1)              revolveBuilder1_XY_1.Section = section1_XY_1              revolveBuilder1_XY_1.Tolerance = 0.001   Dim booleanOperation1_XY_1 As GeometricUtilities.BooleanOperation              booleanOperation1_XY_1 = revolveBuilder1_XY_1.BooleanOperation  Dim limits1_XY_1 As GeometricUtilities.Limits              limits1_XY_1 = revolveBuilder1_XY_1.Limits  Dim angularLimits1_XY_1 As GeometricUtilities.AngularLimits = CType(limits1_XY_1, GeometricUtilities.AngularLimits)              Dim extend1_XY_1 As GeometricUtilities.Extend              extend1_XY_1 = angularLimits1_XY_1.StartExtend              extend1_XY_1.SetValue("0")              Dim extend2_XY_1 As GeometricUtilities.Extend              extend2_XY_1 = angularLimits1_XY_1.EndExtend   126              extend2_XY_1.SetValue("360")              Dim featureOffset1_XY_1 As GeometricUtilities.FeatureOffset              featureOffset1_XY_1 = revolveBuilder1_XY_1.Offset              Dim section2_XY_1 As Section              section2_XY_1 = revolveBuilder1_XY_1.Section               Dim axis1_XY_1 As Axis              axis1_XY_1 = revolveBuilder1_XY_1.Axis              Dim featureOptions1_XY_1 As GeometricUtilities.FeatureOptions              featureOptions1_XY_1 = revolveBuilder1_XY_1.FeatureOptions               Dim curves1_XY_1(0) As Curve              Dim arc_XY_1 As Arc = CType(NXObjectManager.Get(arc_XY_1_Tag), Arc)              curves1_XY_1(0) = arc_XY_1              Dim curveDumbRule1_XY_1 As CurveDumbRule  curveDumbRule1_XY_1 = workPart.ScRuleFactory.CreateRuleCurveDumb(curves1_XY_1)              section2_XY_1.AllowSelfIntersection(False)              Dim rules1_XY_1(0) As SelectionIntentRule              rules1_XY_1(0) = curveDumbRule1_XY_1              Dim nullNXObject_XY_1 As NXObject = Nothing  Dim helpPoint1_XY_1 As Point3d = New Point3d(pOarcx_XY_1, pOarcy_XY_1, pOarcz_XY_1)  section2_XY_1.AddToSection(rules1_XY_1, arc_XY_1, nullNXObject_XY_1, nullNXObject_XY_1, helpPoint1_XY_1, Section.Mode.Create)              revolveBuilder1_XY_1.Section = section2_XY_1               Dim curves2_XY_1(0) As Curve  Dim line1_XY_1 As Line = CType(NXObjectManager.Get(line1_XY_1_Tag), Line)              curves2_XY_1(0) = line1_XY_1              Dim curveDumbRule2_XY_1 As CurveDumbRule  curveDumbRule2_XY_1 = workPart.ScRuleFactory.CreateRuleCurveDumb(curves2_XY_1)              section2_XY_1.AllowSelfIntersection(False)              Dim rules2_XY_1(0) As SelectionIntentRule              rules2_XY_1(0) = curveDumbRule2_XY_1  Dim helpPoint2_XY_1 As Point3d = New Point3d(pOline1x_XY_1, pOline1y_XY_1, pOline1z_XY_1)  section2_XY_1.AddToSection(rules2_XY_1, line1_XY_1, nullNXObject_XY_1, nullNXObject_XY_1, helpPoint2_XY_1, Section.Mode.Create)              revolveBuilder1_XY_1.Section = section2_XY_1               Dim curves3_XY_1(0) As Curve  Dim line2_XY_1 As Line = CType(NXObjectManager.Get(line2_XY_1_Tag), Line)              curves3_XY_1(0) = line2_XY_1  Dim curveDumbRule3_XY_1 As CurveDumbRule  curveDumbRule3_XY_1 = workPart.ScRuleFactory.CreateRuleCurveDumb(curves3_XY_1)              section2_XY_1.AllowSelfIntersection(False)              Dim rules3_XY_1(0) As SelectionIntentRule              rules3_XY_1(0) = curveDumbRule3_XY_1  Dim helpPoint3_XY_1 As Point3d = New Point3d(pOline2x_XY_1, pOline2y_XY_1, pOline2z_XY_1)  section2_XY_1.AddToSection(rules3_XY_1, line2_XY_1,   127 nullNXObject_XY_1, nullNXObject_XY_1, helpPoint3_XY_1, Section.Mode.Create)              revolveBuilder1_XY_1.Section = section2_XY_1               Dim curves4_XY_1(0) As Curve  Dim line3_XY_1 As Line = CType(NXObjectManager.Get(line3_XY_1_Tag), Line)              curves4_XY_1(0) = line3_XY_1              Dim curveDumbRule4_XY_1 As CurveDumbRule  curveDumbRule4_XY_1 = workPart.ScRuleFactory.CreateRuleCurveDumb(curves4_XY_1)              section2_XY_1.AllowSelfIntersection(False)              Dim rules4_XY_1(0) As SelectionIntentRule              rules4_XY_1(0) = curveDumbRule4_XY_1  Dim helpPoint4_XY_1 As Point3d = New Point3d(pOline3x_XY_1, pOline3y_XY_1, pOline3z_XY_1)  section2_XY_1.AddToSection(rules4_XY_1, line3_XY_1, nullNXObject_XY_1, nullNXObject_XY_1, helpPoint4_XY_1, Section.Mode.Create)              revolveBuilder1_XY_1.Section = section2_XY_1              Dim direction_XY_1 As Direction  direction_XY_1 = workPart.Directions.CreateDirection(line3_XY_1, Sense.Forward, SmartObject.UpdateOption.WithinModeling)              Dim nullPoint_XY_1 As Point = Nothing               Dim axis2_XY_1 As Axis  axis2_XY_1 = workPart.Axes.CreateAxis(nullPoint_XY_1, direction_XY_1, SmartObject.UpdateOption.WithinModeling)              revolveBuilder1_XY_1.Axis = axis2_XY_1              Dim axis3_XY_1 As Axis  axis3_XY_1 = workPart.Axes.CreateAxis(nullPoint_XY_1, direction_XY_1, SmartObject.UpdateOption.WithinModeling)  booleanOperation1_XY_1.Type = GeometricUtilities.BooleanOperation.BooleanType.Unite              Dim targetBodies_XY_1(0) As Body  Dim body_XY_1 As Body = CType(NXObjectManager.Get(sphere_body_tag), Body)              targetBodies_XY_1(0) = body_XY_1  booleanOperation1_XY_1.SetTargetBodies(targetBodies_XY_1)              Dim feature_XY_1 As Features.Feature              feature_XY_1 = revolveBuilder1_XY_1.CommitFeature()              revolveBuilder1_XY_1.Destroy()       End If  <<<< The remaining 11 revolve features are excluded but are created in an analogous manner>>>>  Next k  Once all 12 sinter neck locations have been searched, and if a neighbour is present then the revolved feature representing the sinter neck is created, the part of the spherical particle that extends past the sinter neck must be trimmed at the plane of intersection. Otherwise, there would be interference in the modeled solid representing the intersecting particles. This would cause cross section images with overlapping solid phases.   128 The trimming of the solid at the plane of intersection is programmed in the following way. Again, only one trim feature is shown here for example, but there are up to 12 trim features required for each particle. For k = 0 To 11 xK = 3 * k + 0 'x coordinate index in neighbour particle center, x yK = 3 * k + 1 'y coordinate index in neighbour particle center, y zK = 3 * k + 2 'z coordinate index in neighbour particle center, z  If ND(j, k)>=0 Then        Dim rN As Double = ND(j, k) / 2 Dim LengthPOrigin As Double = Math.sqrt((NP(j, xK) - PosX(j)) ^ 2 + (NP(j, yK) - PosY(j)) ^ 2 + (NP(j, zK) - PosZ(j)) ^ 2) Dim SplaneInt As Double = Math.abs(1 / (2 * LengthPOrigin) * (LengthPOrigin ^ 2 - rN ^ 2 + r ^ 2))             Dim Lm As Double = Math.abs(SplaneInt * r / (r + rs_XY_1))             Dim Hn As Double = Math.sqrt((r + rs_XY_1) ^ 2 - SplaneInt ^ 2)             Dim Hm As Double = r / (r + rs_XY_1) * Hn             Dim horizN As Double = Math.sqrt(2) / 2 * SplaneInt             Dim vertN As Double = Math.sqrt(2) / 2 * SplaneInt             Dim phi As Double = Math.acos(SplaneInt / (r + rs_XY_1))             Dim thetaPhi As Double = 45 / 180 * Math.pi + phi                             'Create Datum Plane for Trim Reference             If k = 0 Then    Dim sParcx_XY_1 As Double = (PosX(j) + Math.sqrt(2) / 2 * (Lm - Hm)) 'arc start x    Dim sParcy_XY_1 As Double = (PosY(j) + Math.sqrt(2) / 2 * (Lm + Hm)) 'arc start y                Dim sParcz_XY_1 As Double = PosZ(j) 'arc start z    Dim eParcx_XY_1 As Double = PosX(j) + horizN - Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'arc end x    Dim eParcy_XY_1 As Double = PosY(j) + vertN + Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'arc end y                Dim eParcz_XY_1 As Double = PosZ(j) 'arc end z    Dim rscx_XY_1 As Double = PosX(j) + Math.cos(thetaPhi) * (r + rs_XY_1) 'sinter radius center point x    Dim rscy_XY_1 As Double = PosY(j) + Math.sin(thetaPhi) * (r + rs_XY_1) 'sinter radius center point y                Dim rscz_XY_1 As Double = PosZ(j) 'sinter radius center point z    Dim pOarcx_XY_1 As Double = (1 / 2) * (sParcx_XY_1 + eParcx_XY_1) 'point on arc x    Dim pOarcy_XY_1 As Double = Math.sqrt(rs_XY_1 ^ 2 - (pOarcx_XY_1 - rscx_XY_1) ^ 2) + rscy_XY_1 'point on arc y                Dim pOarcz_XY_1 As Double = PosZ(j) 'point on arc z    Dim ePline2x_XY_1 As Double = sParcx_XY_1 + Math.sqrt(2) / 2 * Hm 'end point line 2 x    Dim ePline2y_XY_1 As Double = sParcy_XY_1 - Math.sqrt(2) / 2 * Hm 'end point line 2 y                Dim ePline2z_XY_1 As Double = PosZ(j) 'end point line 2 z    Dim ePline1x_XY_1 As Double = eParcx_XY_1 + Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'end point line 1 x    Dim ePline1y_XY_1 As Double = eParcy_XY_1 - Math.sqrt(2) / 2 * (Hn - rs_XY_1) 'end point line 1 y                Dim ePline1z_XY_1 As Double = PosZ(j) 'end point line 1 z    Dim mLine1_XY_1 As Double = (ePline1y_XY_1 - eParcy_XY_1) / (ePline1x_XY_1 - eParcx_XY_1) ' slope of line 1    Dim bLine1_XY_1 As Double = eParcy_XY_1 - mLine1_XY_1 * eParcx_XY_1 ' intercept of line 1   129    Dim pOline1x_XY_1 As Double = 0.5 * (eParcx_XY_1 + ePline1x_XY_1) 'point on line 1 x    Dim pOline1y_XY_1 As Double = mLine1_XY_1 * pOline1x_XY_1 + bLine1_XY_1 'point on line 1 y    Dim pOline1z_XY_1 As Double = PosZ(j) 'point on line 1 z    Dim mLine2_XY_1 As Double = (ePline2y_XY_1 - sParcy_XY_1) / (ePline2x_XY_1 - sParcx_XY_1) 'slope line 2    Dim bLine2_XY_1 As Double = (sParcy_XY_1 - mLine2_XY_1 * sParcx_XY_1) 'intercept line 2    Dim pOline2x_XY_1 As Double = 0.5 * (sParcx_XY_1 + ePline2x_XY_1) 'point on line 2 x    Dim pOline2y_XY_1 As Double = mLine2_XY_1 * pOline2x_XY_1 + bLine2_XY_1 'point on line 2 y    Dim pOline2z_XY_1 As Double = PosZ(j) 'point on line 2 z    Dim mLine3_XY_1 As Double = (ePline1y_XY_1 - ePline2y_XY_1) / (ePline1x_XY_1 - ePline2x_XY_1) 'slope line 3    Dim bLine3_XY_1 As Double = ePline2y_XY_1 - mLine3_XY_1 * ePline2x_XY_1 'intercept line 3    Dim pOline3x_XY_1 As Double = 0.5 * (ePline1x_XY_1 + ePline2x_XY_1) 'point on line 3 x    Dim pOline3y_XY_1 As Double = mLine3_XY_1 * pOline3x_XY_1 + bLine3_XY_1 'point on line 3 y    Dim pOline3z_XY_1 As Double = PosZ(j) 'point on line 3 z     Dim sParcx_XY_1_Tag, sParcy_XY_1_Tag, sParcz_XY_1_Tag As NXOpen.Tag 'tags for start point arc coordinates    Dim eParcx_XY_1_Tag, eParcy_XY_1_Tag, eParcz_XY_1_Tag As NXOpen.Tag 'tags for end point arc coordinates    Dim rscax_XY_1_Tag, rscay_XY_1_Tag, rscaz_XY_1_Tag As NXOpen.Tag 'tags for center coordinates of arc    Dim sParc_XY_1_Tag, eParc_XY_1_Tag, cParc_XY_1_Tag As NXOpen.Tag 'tags for arc start, end, center  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcx_XY_1, sParcx_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcy_XY_1, sParcy_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcz_XY_1, sParcz_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, rscx_XY_1, rscax_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, rscy_XY_1, rscay_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, rscz_XY_1, rscaz_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcx_XY_1, eParcx_XY_1_Tag)   130  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcy_XY_1, eParcy_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcz_XY_1, eParcz_XY_1_Tag)    Dim sParc_XY_1_TagArray() As NXOpen.Tag = {sParcx_XY_1_Tag, sParcy_XY_1_Tag, sParcz_XY_1_Tag}    Dim eParc_XY_1_TagArray() As NXOpen.Tag = {eParcx_XY_1_Tag, eParcy_XY_1_Tag, eParcz_XY_1_Tag}    Dim cParc_XY_1_TagArray() As NXOpen.Tag = {rscax_XY_1_Tag, rscay_XY_1_Tag, rscaz_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sParc_XY_1_TagArray, sParc_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, eParc_XY_1_TagArray, eParc_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, cParc_XY_1_TagArray, cParc_XY_1_Tag)    Dim arcPnts_XY_1_TagArray() As NXOpen.Tag = {cParc_XY_1_Tag, sParc_XY_1_Tag, eParc_XY_1_Tag}                Dim arc_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateArcCenter2Pnts(object_in_part, SmartObject.UpdateOption.WithinModeling, arcPnts_XY_1_TagArray, arc_XY_1_Tag)     Dim sPline2x_XY_1_Tag, sPline2y_XY_1_Tag, sPline2z_XY_1_Tag As NXOpen.Tag 'tags for line 2 start point coordinates    Dim ePline2x_XY_1_Tag, ePline2y_XY_1_Tag, ePline2z_XY_1_Tag As NXOpen.Tag 'tags for line 2 end point coordinates    Dim sPline2_XY_1_Tag, ePline2_XY_1_Tag As NXOpen.Tag 'tags for line 2 start and end points  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcx_XY_1, sPline2x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcy_XY_1, sPline2y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, sParcz_XY_1, sPline2z_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2x_XY_1, ePline2x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2y_XY_1, ePline2y_XY_1_Tag)   131  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2z_XY_1, ePline2z_XY_1_Tag)    Dim sPline2_XY_1_TagArray() As NXOpen.Tag = {sPline2x_XY_1_Tag, sPline2y_XY_1_Tag, sPline2z_XY_1_Tag}    Dim ePline2_XY_1_TagArray() As NXOpen.Tag = {ePline2x_XY_1_Tag, ePline2y_XY_1_Tag, ePline2z_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sPline2_XY_1_TagArray, sPline2_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2_XY_1_TagArray, ePline2_XY_1_Tag)    Dim line2Pnts_XY_1_TagArray() As NXOpen.Tag = {sPline2_XY_1_Tag, ePline2_XY_1_Tag}                Dim line2_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateLineTwoPoints(object_in_part, SmartObject.UpdateOption.WithinModeling, line2Pnts_XY_1_TagArray, line2_XY_1_Tag)     Dim sPline1x_XY_1_Tag, sPline1y_XY_1_Tag, sPline1z_XY_1_Tag As NXOpen.Tag 'tags for line 1 start coordinates    Dim ePline1x_XY_1_Tag, ePline1y_XY_1_Tag, ePline1z_XY_1_Tag As NXOpen.Tag 'tags for line 2 end coordinates    Dim sPline1_XY_1_Tag, ePline1_XY_1_Tag As NXOpen.Tag 'tags for line 1 start and end points  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcx_XY_1, sPline1x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcy_XY_1, sPline1y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, eParcz_XY_1, sPline1z_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1x_XY_1, ePline1x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1y_XY_1, ePline1y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1z_XY_1, ePline1z_XY_1_Tag)    Dim sPline1_XY_1_TagArray() As NXOpen.Tag = {sPline1x_XY_1_Tag, sPline1y_XY_1_Tag, sPline1z_XY_1_Tag}    Dim ePline1_XY_1_TagArray() As NXOpen.Tag = {ePline1x_XY_1_Tag, ePline1y_XY_1_Tag, ePline1z_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part,   132 SmartObject.UpdateOption.WithinModeling, sPline1_XY_1_TagArray, sPline1_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1_XY_1_TagArray, ePline1_XY_1_Tag)    Dim line1Pnts_XY_1_TagArray() As NXOpen.Tag = {sPline1_XY_1_Tag, ePline1_XY_1_Tag}                Dim line1_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateLineTwoPoints(object_in_part, SmartObject.UpdateOption.WithinModeling, line1Pnts_XY_1_TagArray, line1_XY_1_Tag)     Dim sPline3x_XY_1_Tag, sPline3y_XY_1_Tag, sPline3z_XY_1_Tag As NXOpen.Tag 'tags for line 3 start coordinates    Dim ePline3x_XY_1_Tag, ePline3y_XY_1_Tag, ePline3z_XY_1_Tag As NXOpen.Tag 'tags for line 3 end coordinates    Dim sPline3_XY_1_Tag, ePline3_XY_1_Tag As NXOpen.Tag 'tags for line 3 start and end points  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1x_XY_1, sPline3x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1y_XY_1, sPline3y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline1z_XY_1, sPline3z_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2x_XY_1, ePline3x_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2y_XY_1, ePline3y_XY_1_Tag)  theUFSession.So.CreateScalarDouble(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline2z_XY_1, ePline3z_XY_1_Tag)     Dim sPline3_XY_1_TagArray() As NXOpen.Tag = {sPline3x_XY_1_Tag, sPline3y_XY_1_Tag, sPline3z_XY_1_Tag}     Dim ePline3_XY_1_TagArray() As NXOpen.Tag = {ePline3x_XY_1_Tag, ePline3y_XY_1_Tag, ePline3z_XY_1_Tag}  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, sPline3_XY_1_TagArray, sPline3_XY_1_Tag)  theUFSession.So.CreatePoint3Scalars(object_in_part, SmartObject.UpdateOption.WithinModeling, ePline3_XY_1_TagArray, ePline3_XY_1_Tag)    Dim line3Pnts_XY_1_TagArray() As NXOpen.Tag = {sPline3_XY_1_Tag, ePline3_XY_1_Tag}                Dim line3_XY_1_Tag As NXOpen.Tag  theUFSession.So.CreateLineTwoPoints(object_in_part,   133 SmartObject.UpdateOption.WithinModeling, line3Pnts_XY_1_TagArray, line3_XY_1_Tag)     Dim nullFeatures_Feature_XY_1 As Features.Feature = Nothing    Dim arc_XY_1 As Arc = CType(NXObjectManager.Get(arc_XY_1_Tag), Arc)                Dim direction_XY_1 As Direction    Dim line3_XY_1 As Line = CType(NXObjectManager.Get(line3_XY_1_Tag), Line)    direction_XY_1 = workPart.Directions.CreateDirection(line3_XY_1, Sense.Forward, SmartObject.UpdateOption.WithinModeling)    Dim body_XY_1 As Body = CType(NXObjectManager.Get(sphere_body_tag), Body)                 Dim datumPlaneBuilder1_XY_1 As Features.DatumPlaneBuilder    datumPlaneBuilder1_XY_1 = workPart.Features.CreateDatumPlaneBuilder(nullFeatures_Featur e_XY_1)                Dim plane1_XY_1 As Plane                plane1_XY_1 = datumPlaneBuilder1_XY_1.GetPlane()                Dim scalar1_XY_1 As Scalar    scalar1_XY_1 = workPart.Scalars.CreateScalar(1.0, Scalar.DimensionalityType.None, SmartObject.UpdateOption.WithinModeling)                Dim pointD1_XY_1 As Point    pointD1_XY_1 = workPart.Points.CreatePoint(arc_XY_1, scalar1_XY_1, SmartObject.UpdateOption.WithinModeling)    Dim directionD1_XY_1 As Direction    directionD1_XY_1 = workPart.Directions.CreateDirection(arc_XY_1, scalar1_XY_1, Direction.OnCurveOption.Tangent, Sense.Forward, SmartObject.UpdateOption.WithinModeling)  plane1_XY_1.SetMethod(PlaneTypes.MethodType.PointDir)                Dim geom1_XY_1(1) As NXObject                geom1_XY_1(0) = pointD1_XY_1                geom1_XY_1(1) = directionD1_XY_1                plane1_XY_1.SetGeometry(geom1_XY_1)  plane1_XY_1.SetAlternate(PlaneTypes.AlternateType.One)                plane1_XY_1.Evaluate()  plane1_XY_1.SetMethod(PlaneTypes.MethodType.PointDir)                Dim geom2_XY_1(1) As NXObject                geom2_XY_1(0) = pointD1_XY_1                geom2_XY_1(1) = direction_XY_1                plane1_XY_1.SetGeometry(geom2_XY_1)  plane1_XY_1.SetAlternate(PlaneTypes.AlternateType.One) plane1_XY_1.Evaluate()                Dim success1_XY_1 As Boolean                success1_XY_1 = direction_XY_1.ReverseDirection()  plane1_XY_1.SetMethod(PlaneTypes.MethodType.PointDir)                Dim geom3_XY_1(1) As NXObject                geom3_XY_1(0) = pointD1_XY_1                geom3_XY_1(1) = direction_XY_1                plane1_XY_1.SetGeometry(geom3_XY_1)  plane1_XY_1.SetAlternate(PlaneTypes.AlternateType.One)                plane1_XY_1.Evaluate()                Dim flip1_XY_1 As Boolean   134                flip1_XY_1 = plane1_XY_1.Flip                Dim featureD1_XY_1 As Features.Feature                featureD1_XY_1 = datumPlaneBuilder1_XY_1.CommitFeature()    Dim datumPlaneFeature1_XY_1 As Features.DatumPlaneFeature = CType(featureD1_XY_1, Features.DatumPlaneFeature)                Dim datumPlane1_XY_1 As DatumPlane                datumPlane1_XY_1 = datumPlaneFeature1_XY_1.DatumPlane                datumPlane1_XY_1.SetReverseSection(False)                datumPlaneBuilder1_XY_1.Destroy()                 Dim trimBodyBuilder1_XY_1 As Features.TrimBodyBuilder    trimBodyBuilder1_XY_1 = workPart.Features.CreateTrimBodyBuilder(nullFeatures_Feature_ XY_1)                Dim target1_XY_1() As Body                target1_XY_1 = trimBodyBuilder1_XY_1.GetTargets()                Dim nXObject1_XY_1 As NXObject                nXObject1_XY_1 = trimBodyBuilder1_XY_1.Tool                Dim direction3T_XY_1 As Features.TrimBodyBuilder.DirectionType                direction3T_XY_1 = trimBodyBuilder1_XY_1.TrimDirection                trimBodyBuilder1_XY_1.AddTarget(body_XY_1)                Dim scCollector1_XY_1 As ScCollector                scCollector1_XY_1 = workPart.ScCollectors.CreateCollector()                Dim faces1_XY_1(0) As DatumPlane                faces1_XY_1(0) = datumPlane1_XY_1                Dim faceDumbRule1_XY_1 As FaceDumbRule    faceDumbRule1_XY_1 = workPart.ScRuleFactory.CreateRuleFaceDatum(faces1_XY_1)                Dim rules1T_XY_1(0) As SelectionIntentRule                rules1T_XY_1(0) = faceDumbRule1_XY_1                scCollector1_XY_1.ReplaceRules(rules1T_XY_1, False)                trimBodyBuilder1_XY_1.Tool = scCollector1_XY_1                Dim feature2T_XY_1 As Features.Feature                feature2T_XY_1 = trimBodyBuilder1_XY_1.CommitFeature()                trimBodyBuilder1_XY_1.Destroy() End If  <<<< The remaining 11 trim features for the remaining sinter neck locations are excluded but are created in an analogous manner>>>>  Next k Next j End Module The revolve and trim features have now been created for the first particle j. The for loop for j now increments to the next j, representing the next particle on the lattice.   135 Appendix B Computer Code for Calculation of Total Three Phase Boundary Length and Total Porosity Matlab Code for Calculation of Total Three Phase Boundary Length The total three phase boundary length is calculated based on the position, size, and phase of each particle on the lattice together with the positions, sizes, and phases of the nearest neighbor particles for each particle on the lattice. The function calcTPBlength.m is used to calculate the total three phase boundary length. As described in Appendix One, PM is the matrix containing the x, y, and z coordinates of the particle centers in each column. The number of rows in PM is equal to the number of particles on the lattice. DM is an array containing the diameters of each particle on the lattice. sPhase is an array containing a 1 or a 0 in each element, corresponding to whether the particle is an electronic conductor or ionic conductor, respectively. PN is a matrix describing the position of each of the nearest neighbour particles for each particle on the lattice. There are 36 columns, corresponding to the x, y, and z coordinates for all 12 possible nearest neighbours. If a nearest neighbour position is unoccupied, a value of -1 is used for the coordinate as an indicator that the particle is not present on the lattice. DN is a matrix containing the diameters of the neighbour particles for each particle on the lattice and contains a -1 where no neighbor is present. rs is the radius of curvature for the sinter necks. N is the number of particles on the lattice. a is the lattice constant and maxP is the maximum coordinate magnitude, which corresponds to the length of one side of the overall modeled geometry. function TPB = calcTPBlength(PM, DM, sPhase, PN, DN, rs, N, a, maxP)  The first task in calculating the total three phase boundary length is to identify all the opposite phase intersections between neighbour particles. This is accomplished using calcIntersections.m as shown below. Intersect is a matrix with six columns. Each row represents an intersection between opposite phase particles and the six columns contain the following information: the particle number on the lattice (corresponding to the row number in PM containing its center coordinates) pNumC, the location of intersection on the particle   136 (plane and quadrant) lNumC, the neighbor particle number (corresponding to the row in PM containing its center coordinates) pNumN, the location of intersection on the neighbour particle (plane and quadrant) lNumN, the phase of the particle sPhase(pNumC) and the phase of the neighbour particle sPhase(pNumN). IntF is a matrix which maps the locations on a particle to the locations of intersection with its nearest neighbour particles. A description of these locations is presented in Chapter Three. %******************************************************************* function [Intersect IntF] = calcIntersections(PM, sPhase, PN, N);  % Each intersection interferes with a max of four other intersections: %XY Plane: 1 = Q1, 2 = Q2, 3 = Q3, 4 = Q4 IntF(1,:) = [   9   5   12  8]; IntF(2,:) = [   9   6   12  7]; IntF(3,:) = [   10  6   11  7]; IntF(4,:) = [   8   11  5   10];  %XZ Plane: 5 = Q1, 6 = Q2, 7 = Q3, 8 = Q4 IntF(5,:) = [   1   9   4   10]; IntF(6,:) = [   9   2   10  3]; IntF(7,:) = [   2   12  11  3]; IntF(8,:) = [   1   12  11  4];  %YZ Plane: 9 = Q1, 10 = Q2, 11 = Q3, 12 = Q4 IntF(9,:) = [   5   1   6   2]; IntF(10,:)= [   5   4   3   6]; IntF(11,:)= [   7   3   8   4]; IntF(12,:)= [   1   8   7   2];  % MPL is an array where the index ranges from 1 to 12. For example, a 3 in index 1 represents that location 1 on a particle intersects with location 3 on the neighbour particle. MPL = [3; 4; 1; 2; 7; 8; 5; 6; 11; 12; 9; 10]; %map particle location  The following code identifies all opposite phase intersections in the lattice. IntersectCount = 1; for i = 1:N      for k = 1:12         kx = 3*k - 2;   % index of x coordinate in PM of neighbor particle         ky = 3*k - 1;   % index of y coordinate in PM of neighbor particle         kz = 3*k - 0;   % index of z coordinate in PM of neighbor particle          if PN(i,kx) < 0             % No neighbor here so advance to next neighbor position             continue         else             % Record the particle number of neighbor.             pNumN = find(PM(:,1)==PN(i,kx)&PM(:,2)==PN(i,ky)&PM(:,3)==PN(i,kz));             % Record particle number of current particle             pNumC = i;             % Record location of intersection on current particle             lNumC = k;             % Record location of intersection on neighbor particle   137             lNumN = MPL(k); Intersection(IntersectCount,:) = [pNumC lNumC pNumN lNumN ... sPhase(pNumC) sPhase(pNumN)];         end         IntersectCount = IntersectCount + 1;     end end  sizeIntersect = size(Intersection); if  sizeIntersect(1) == 0    Intersection = [];    IntF = [];    return end  %add a row for stop criterion for the following while loop; sizeIntersect = size(Intersection); Intersection(sizeIntersect(1)+1,:)=[max(Intersection(:,1))+1, 0 , 0 , 0, 0, 0]; sizeIntersect = size(Intersection);  % Remove rows in the Intersection matrix that correspond to the same particle % pair intersection. i = 1; while Intersection(i,1) <= N     tRN = Intersection(i,1);     tRL = Intersection(i,2);     Ind = find((Intersection(:,3) == tRN) & (Intersection(:,4) == tRL));     if isempty(Ind) == 0         Intersection(Ind,:) = [];     end     i = i + 1;     testWhile = size(Intersection);     if i == testWhile(1)         break     end end  sizeIntersect = size(Intersection); Intersection(sizeIntersect(1),:) = []; %Remove flag row  % Remove rows in Intersection Matrix for which the intersecting particles are % of the same phase. count = 1; sizeIntersect = size(Intersection); for i = 1:sizeIntersect(1)     if Intersection(i,5) == Intersection(i,6)         RemoveRow(count) = i;         count = count + 1;     end end  %The matrix Intersection now contains only intersections of opposite phase %particles with no duplicates. Intersection(RemoveRow,:) = []; %*******************************************************************  Calculating the parameterized circles representing the sinter neck perimeters as described in Chapter Three, section 3.2, requires definition of the min and max angle corresponding to   138 0 and 2pi for a full circle, as well as the resolution (corresponding to the number of points used to represent the circle). %for circle parameterization, s = R*sin(theta), t =R*cos(theta) %[Point on Circle] = [Point on plane] + s*v1 + t*v2 where v1 and v2 are %orthogonal unit vectors spanning the plane (here, the 'plane' refers to the %plane of intersection of the sinter neck) and R is the radius of this %intersection. tol = 1e-3; increment = 0.001; mintheta = increment; maxtheta = 2*pi;  Local Cartesian planes are defined at the center of each particle to define the permissible particle intersection locations. Unit normal vectors to the planes of intersection as well as unit vectors parallel to the plane of intersection are required to identify the locations at which sinter necks may interfere with each other, less any interference by other particles or the geometry envelope. %unit normal vectors, n, and unit vectors, u, for all 12 planes of %intersection. n is normal to plane of intersection. u is parallel to plane %of intersection. %XY Plane: 1 = Q1, 2 = Q2, 3 = Q3, 4 = Q4 n(1,:) = [ sqrt(2)/2, sqrt(2)/2, 0];    u(1,:) = [-sqrt(2)/2, sqrt(2)/2, 0]; n(2,:) = [-sqrt(2)/2, sqrt(2)/2, 0];    u(2,:) = [+sqrt(2)/2, sqrt(2)/2, 0]; n(3,:) = [-sqrt(2)/2,-sqrt(2)/2, 0];    u(3,:) = [+sqrt(2)/2,-sqrt(2)/2, 0]; n(4,:) = [ sqrt(2)/2,-sqrt(2)/2, 0];    u(4,:) = [-sqrt(2)/2,-sqrt(2)/2, 0];  %XZ Plane: 5 = Q1, 6 = Q2, 7 = Q3, 8 = Q4 n(5,:) = [ sqrt(2)/2, 0, sqrt(2)/2];    u(5,:) = [-sqrt(2)/2, 0, sqrt(2)/2]; n(6,:) = [-sqrt(2)/2, 0, sqrt(2)/2];    u(6,:) = [+sqrt(2)/2, 0, sqrt(2)/2]; n(7,:) = [-sqrt(2)/2, 0,-sqrt(2)/2];    u(7,:) = [+sqrt(2)/2, 0,-sqrt(2)/2]; n(8,:) = [ sqrt(2)/2, 0,-sqrt(2)/2];    u(8,:) = [-sqrt(2)/2, 0,-sqrt(2)/2];  %YZ Plane: 9 = Q1, 10 = Q2, 11 = Q3, 12 = Q4 n(9,:) = [ 0, sqrt(2)/2, sqrt(2)/2];    u(9,:) = [ 0,-sqrt(2)/2, sqrt(2)/2]; n(10,:)= [ 0,-sqrt(2)/2, sqrt(2)/2];    u(10,:)= [ 0,+sqrt(2)/2, sqrt(2)/2]; n(11,:)= [ 0,-sqrt(2)/2,-sqrt(2)/2];    u(11,:)= [ 0,+sqrt(2)/2,-sqrt(2)/2]; n(12,:)= [ 0, sqrt(2)/2,-sqrt(2)/2];    u(12,:)= [ 0,-sqrt(2)/2,-sqrt(2)/2];  For every opposite phase intersection, the radius of the sinter neck must be calculated. The three phase boundary length is the circumference of intersection between the two opposite phase particles. for j = 1:size(Intersect,1)      i = Intersect(j,1);    %Particle number => row number in PM     r = DM(i)/2;    %Radius of particle     kInt = Intersect(j,2);  %location of intersection on particle i      kxInt = 3*kInt - 2;   % index of x coordinate in PM of neighbor particle     kyInt = 3*kInt - 1;   % index of y coordinate in PM of neighbor particle   139     kzInt = 3*kInt - 0;   % index of z coordinate in PM of neighbor particle      % Distance between center coordinates of particles i and k     L_Int = sqrt((PN(i,kxInt) - PM(i,1))^2 + ...                  (PN(i,kyInt) - PM(i,2))^2 + ...                  (PN(i,kzInt) - PM(i,3))^2);      rN = DN(i,kInt)/2;       % Radius of neighbor particle      % Distance from particle center i to plane of intersection of i     % with kInt is SplaneInt     SplaneInt = abs(1/(2*L_Int) * (L_Int ^ 2 - rN ^ 2 + r ^ 2));  In vector form, the position vector, vLmInt, from the center of particle i to the origin of the circle describing the sinter neck (Point O in Figure 3-4, and point N in Figure 2-6) is represented by the distance SplaneInt in the direction of the normal vector to the plane of intersection.     vLmInt = SplaneInt*n(kInt,:);  The corresponding coordinates of the center of the sinter neck (point O in Figure 3-4 and point N in Figure 2-6) are therefore:     pLmInt = PM(i,:) + vLmInt; The radius of the sinter neck, HmInt, formed between particle i and neighbour with location kInt is:     HmInt = sqrt((r + rs(kInt)) ^ 2 - SplaneInt ^ 2) - rs(kInt);  The three phase boundary length would be equal to HmInt multiplied by 2pi radians if there were no interference by other sinter necks on particle i. Therefore, the interference between any sinter necks must be determined.     k = IntF(kInt,:);     kx = 3.*k - 2;       % index of x coordinate in PM of neighbor particle     ky = 3.*k - 1;  % index of y coordinate in PM of neighbor particle     kz = 3.*k - 0;  % index of z coordinate in PM of neighbor particle  There can be up to four interfering sinter necks for any given sinter neck. for w = 1:4         % Check if there is a neighbor particle at location IntF(k,1).         % If true then this sinterneck interferes with sinterneck of         % particle i, location k. Therefore, if true, need to calculate the         % equation of the circle of intersection of sinter neck on         % particle i, location IntF(k,1).         if PN(i,kx(w)) >= 0             kCheck1(w) = k(w);             % Distance between center coordinates of particles i and k             L_Int1 = sqrt((PN(i,kx(w)) - PM(i,1))^2 + ...                           (PN(i,ky(w)) - PM(i,2))^2 + ...   140                           (PN(i,kz(w)) - PM(i,3))^2);             rN_Int1 = DN(i,k(w))/2;   % Radius of neighbor particle k             Splane_Int1 = abs(1/(2*L_Int1)*(L_Int1 ^ 2 - rN_Int1 ^ 2 + r ^ 2));             vLm_Int1 = Splane_Int1*n(k(w),:);             pLm_Int1 = PM(i,:) + vLm_Int1;             Hm_Int1 = sqrt((r + rs(k(w))) ^ 2 - Splane_Int1 ^ 2) - rs(k(w));             vHm_Int1 = Hm_Int1*u(k(w),:); Two orthogonal unit vectors spanning the plane of intersection between the two neighbour particles are v1u_Int1 and v2u_Int1.             v1u_Int1 = vHm_Int1/norm(vHm_Int1);             v2u_Int1 = cross(u(k(w),:),n(k(w),:));  The points describing the interfering sinter neck perimeter are:             s = Hm_Int1*sin(theta);  t = Hm_Int1*cos(theta);             Pntx_Int1 = pLm_Int1(1) + s.*v1u_Int1(1) + t.*v2u_Int1(1);             Pnty_Int1 = pLm_Int1(2) + s.*v1u_Int1(2) + t.*v2u_Int1(2);             Pntz_Int1 = pLm_Int1(3) + s.*v1u_Int1(3) + t.*v2u_Int1(3);  The location of intersection with the interfering sinter neck is calculated using the function sinterIntersect.m as follows: theta_Int1 = interIntersect(pLmInt,Pntx_Int1,Pnty_Int1,Pntz_Int1,...                  HmInt, theta, mintheta, maxtheta, v1u_Int1, v2u_Int1, ... vHm_Int1, Hm_Int1);  %*************************************************************************** function theta_Int1 = sinterIntersect(pLmInt,Pntx,Pnty,Pntz, HmInt, theta, mintheta, maxtheta, v1u_Int1, v2u_Int1, vHm_Int1, Hm_Int1)  The distance, fcn, between the center of sinter neck and the points on the circle describing the interfering sinter neck is: for i = 1:length(Pntx)     fcn(i) = sqrt((Pntx(i)-pLmInt(1))^2 + (Pnty(i)-pLmInt(2))^2 + ...         (Pntz(i)-pLmInt(3))^2) - HmInt; end  %fit cubic spline cs = csapi(theta,fcn);  If the global minimum of the function fcn is greater than zero, then there is no interfering sinter neck. If this minimum is less than zero or equal to zero then the interfering sinter neck intersects, or touches, the sinter neck, respectively. The angle, theta_Int1, for which the distance, fcn, is zero, is used to calculate the coordinates for which the interfering sinter necks overlap. check = fnmin(cs,[mintheta maxtheta]); if check > 0     %no intersection   141     theta_Int1 = []; else     z = fnzeros(cs,[mintheta maxtheta]);     theta_Int1 = z(1,:);     if numel(theta_Int1) == 1         %intersect at point         theta_Int1 = [];     end end %***************************************************************************  The coordinates of the intersection points between two interfering sinter necks are pInt1 and pInt2. Then, vectors int1_v1 and int1_v1 from the center of the sinter neck to the intersection points are calculated and the angle between them, ang_Int1, is determined. if isempty(theta_Int1) == 0                 p1x_Int1 = pLm_Int1(1)+Hm_Int1*sin(theta_Int1(1))*v1u_Int1(1)...                     + Hm_Int1*cos(theta_Int1(1))*v2u_Int1(1);                 p1y_Int1 = pLm_Int1(2)+Hm_Int1*sin(theta_Int1(1))*v1u_Int1(2)...                     + Hm_Int1*cos(theta_Int1(1))*v2u_Int1(2);                 p1z_Int1 = pLm_Int1(3)+Hm_Int1*sin(theta_Int1(1))*v1u_Int1(3)...                     + Hm_Int1*cos(theta_Int1(1))*v2u_Int1(3);                 pInt1(w,:) = [p1x_Int1 p1y_Int1 p1z_Int1];                  p2x_Int1 = pLm_Int1(1)+Hm_Int1*sin(theta_Int1(2))*v1u_Int1(1)...                     + Hm_Int1*cos(theta_Int1(2))*v2u_Int1(1);                 p2y_Int1 = pLm_Int1(2)+Hm_Int1*sin(theta_Int1(2))*v1u_Int1(2)...                     + Hm_Int1*cos(theta_Int1(2))*v2u_Int1(2);                 p2z_Int1 = pLm_Int1(3)+Hm_Int1*sin(theta_Int1(2))*v1u_Int1(3)...                     + Hm_Int1*cos(theta_Int1(2))*v2u_Int1(3);                 pInt2(w,:) = [p2x_Int1 p2y_Int1 p2z_Int1];                 % Create vectors from center to intersection points 1 and 2                 % and calculate the angle between them                 int1_v1(w,:) = pInt1(w,:) - pLmInt;                 int1_v2(w,:) = pInt2(w,:) - pLmInt; ang_Int1(w) = atan2(norm(cross(int1_v1(w,:),int1_v2(w,:))),...                     dot(int1_v1(w,:),int1_v2(w,:)));             else                 kCheck1(w) = 0;                 ang_Int1(w) = 0;                 int1_v1(w,:) = zeros(1,3);                 int1_v2(w,:) = zeros(1,3);                 pInt2(w,:) = zeros(1,3);                 pInt1(w,:) = zeros(1,3);             end         else             kCheck1(w) = 0;             ang_Int1(w) = 0;             int1_v1(w,:) = zeros(1,3);             int1_v2(w,:) = zeros(1,3);             pInt2(w,:) = zeros(1,3);             pInt1(w,:) = zeros(1,3);         end  Here, the w for loop iterates through the remaining three possible interfering sinter neck locations and calculates the intersection points and angle as described above.   142     end  As described in Chapter Three, Figure 3-3, it is possible for the interfering sinter necks to interfere with each other. Therefore, the net angles of interference, subt_ang1 and subt_ang2, must be determined. In Figure 3-3, kCheck1(1) and kCheck1(2) represent the angle between the intersection points for curves a and c or for curves b and d. and kCheck1(3) and kCheck1(4) represent the angles between the intersection points for curves b and d or for curves a and c.      if kCheck1(1) > 0 && kCheck1(2) > 0          kCheck_ang1a = acos(dot(int1_v1(1,:),int1_v1(2,:))/...              (norm(int1_v1(1,:))*norm(int1_v1(2,:))));          kCheck_ang2a = acos(dot(int1_v1(1,:),int1_v2(2,:))/...              (norm(int1_v1(1,:))*norm(int1_v2(2,:))));          kCheck_ang3a = acos(dot(int1_v2(1,:),int1_v1(2,:))/...              (norm(int1_v2(1,:))*norm(int1_v1(2,:))));          kCheck_ang4a = acos(dot(int1_v2(1,:),int1_v2(2,:))/...              (norm(int1_v2(1,:))*norm(int1_v2(2,:))));           %check if the interfering sinternecks overlap each other          kCheck_a = sort([kCheck_ang1a kCheck_ang2a kCheck_ang3a kCheck_ang4a]);          if kCheck_a(end) > sum(ang_Int1(1:2))%(kCheck_a(2) + kCheck_a(3))              subt_ang1 = 0;          else              subt_ang1 = kCheck_a(1);          end       else          subt_ang1 = 0;      end      if kCheck1(3) > 0 && kCheck1(4) > 0          kCheck_ang1b = acos(dot(int1_v1(3,:),int1_v1(4,:))/...              (norm(int1_v1(3,:))*norm(int1_v1(4,:))));          kCheck_ang2b = acos(dot(int1_v1(3,:),int1_v2(4,:))/...              (norm(int1_v1(3,:))*norm(int1_v2(4,:))));          kCheck_ang3b = acos(dot(int1_v2(3,:),int1_v1(4,:))/...              (norm(int1_v2(3,:))*norm(int1_v1(4,:))));          kCheck_ang4b = acos(dot(int1_v2(3,:),int1_v2(4,:))/...              (norm(int1_v2(3,:))*norm(int1_v2(4,:))));          kCheck_b = sort([kCheck_ang1b kCheck_ang2b kCheck_ang3b kCheck_ang4b]);           %check if the interfering sinternecks overlap each other          if kCheck_b(end) > sum(ang_Int1(3:4))%(kCheck_b(2) + kCheck_b(3))              subt_ang2 = 0;          else              subt_ang2 = kCheck_b(1);           end      else          subt_ang2 = 0;      end  The final task for calculating the total three phase boundary length is to determine if the sinter neck is on a face, edge, or corner of the overall geometry envelope. The three phase   143 boundary length for the sinter neck is then reduced by pi radians if the sinter neck location also intersects the geometry envelope. % Check for lattice site to be internal to packing:     if ((PM(i,1) > (0 + tol)) ...                               % x is not zero      && ~(PM(i,1) > (maxP + tol) && PM(i,1) < (maxP - tol))) ... % x is not maxP      &&((PM(i,2) > (0 + tol)) ...                                % y is not zero      &&~(PM(i,2) > (maxP + tol) && PM(i,2) < (maxP - tol)))... % y is not maxP      &&((PM(i,3) > (0 + tol) ...                                 % z is not zero      &&~(PM(i,3) > (maxP + tol) && PM(i,3) < (maxP - tol))))     % z is not maxP         %This is an internal lattice site (not on any face, corner, or edge)         % Calculate the angle exposed to pore space         thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         %Calculate TPB length         countInternal = countInternal + 1;     end     % Check for lattice site be located parallel YZ plane at x = 0 or x     % = maxP. Further check if the intersection location lies in the     % same plane.     if ((PM(i,1) < (0 + tol)) ...                                   % x is zero      ||((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol)))) ... % x is maxP      &&((PM(i,2) > (0 + tol)) ...                                % y is not zero     &&~((PM(i,2) > (maxP + tol)) && (PM(i,2) < (maxP - tol))))...% y is not maxP      &&((PM(i,3) > (0 + tol)...                                  % z is not zero      &&~(PM(i,3) > (maxP + tol) && PM(i,3) < (maxP - tol))))     % z is not maxP         if any(ismember(9:12,kInt)) == 1             % Calculate the angle exposed to pore space             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);             countInternalPlane = countInternalPlane + 1;             %disp('YZ internal plane and neighbor int is in YZ plane')         else             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);             countInternalPlane = countInternalPlane + 1;             %disp('YZ internal plane and neighbor int is NOT in YZ plane')         end     end     % Check for lattice site be located parallel XZ plane at y = 0 or y     % = maxP. Further check if the intersection location lies in the     % same plane.     if ((PM(i,2) < (0 + tol)) ...                                 % y is zero      ||((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol)))) ... % y is maxP      &&((PM(i,1) > (0 + tol)) ...                                % x is not zero     &&~((PM(i,1) > (maxP + tol)) && (PM(i,1) < (maxP - tol))))...% x is not maxP      &&((PM(i,3) > (0 + tol))...                                 % z is not zero     &&~((PM(i,3) > (maxP + tol)) && (PM(i,3) < (maxP - tol))))  % z is not maxP         if any(ismember(5:8,kInt)) == 1             % Calculate the angle exposed to pore space             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);             countInternalPlane = countInternalPlane + 1;             %disp('XZ internal plane and neighbor int is in XZ plane')         else             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);             countInternalPlane = countInternalPlane + 1;             %disp('XZ internal plane and neighbor int is NOT in XZ plane')         end     end     % Check for lattice site be located parallel XY plane at z = 0 or z     % = maxP. Further check if the intersection location lies in the     % same plane.   144     if ((PM(i,3) < (0 + tol)) ...                              % z is zero      ||((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol)))) ...  % z is maxP      &&((PM(i,2) > (0 + tol)) ...                                % y is not zero     &&~((PM(i,2) > (maxP + tol)) && (PM(i,2) < (maxP - tol))))...% y is not maxP      &&((PM(i,1) > (0 + tol)) ...                               % x is not zero     &&~((PM(i,1) > (maxP + tol)) && (PM(i,1) < (maxP - tol))))   % x is not maxP     % This line used to read:     % &&~((PM(i,1) > (maxP + tol)) && (PM(i,3) < (maxP - tol)))) % x is not maxP         if any(ismember(1:4,kInt)) == 1             % Calculate the angle exposed to pore space             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);             countInternalPlane = countInternalPlane + 1;             %disp('XY internal plane and neighbor int is in XY plane')         else             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);             countInternalPlane = countInternalPlane + 1;             %disp('XY internal plane and neighbor int is NOT in XY plane')         end     end     %/check if lattice site is on an edge parallel to x axis at y=0, z=0     if ((PM(i,1) > (0 + tol)) ...           % (x is not zero AND x is not maxP)    && ~((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol))))...      && (PM(i,2) < (0 + tol)) ...            % AND y is zero      && (PM(i,3) < (0 + tol))                % AND z is zero         if kInt == 9             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is on an edge parallel to x axis at y=0,z=maxP     if ((PM(i,1) > (0 + tol)) ...       % (x is not zero AND x is not maxP)    && ~((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol))))...      && (PM(i,2) < (0 + tol)) ...       % AND y is zero      &&((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol))) % AND z is maxP         if kInt == 12             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is on an edge parallel to x axis at y=maxP,z=maxp     if ((PM(i,1) > (0 + tol)) ...        % (x is not zero AND x is not maxP)    && ~((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol))))...     && ((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol)))... % AND y is maxP     && ((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol)))  % AND z is maxP         if kInt == 11             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is on an edge parallel to x axis at y=maxP,z=0     if ((PM(i,1) > (0 + tol)) ...         % (x is not zero AND x is not maxP)    && ~((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol))))...     && ((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol)))... % AND y is maxP     &&  (PM(i,3) < (0 + tol))              % AND z is zero         if kInt == 10             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);   145         end     end     %/check if lattice site is parallel to y axis at x = 0, z=0     if ((PM(i,2) > (0 + tol)) ...        % (y is not zero AND y is not maxP)    && ~((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol))))...      && (PM(i,1) < (0 + tol)) ...       % AND x is zero      && (PM(i,3) < (0 + tol))          % AND z is zero         if kInt == 5             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is parallel to y axis, x = maxP, z = 0     if ((PM(i,2) > (0 + tol)) ...       % (y is not zero AND y is not maxP)    && ~((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol))))...     && ((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol)))... % AND x is maxP      && (PM(i,3) < (0 + tol))                                % AND z is zero         if kInt == 6             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is parallel to y axis at x = maxP, z=maxP     %///// checked if statement     if ((PM(i,2) > (0 + tol)) ...         % (y is not zero AND y is not maxP)    && ~((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol))))...     && ((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol))) ... % x is maxP     && ((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol)))  % z is maxP         if kInt == 7             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is parallel to y axis at x=0,z=maxP     if ((PM(i,2) > (0 + tol)) ...       % (y is not zero AND y is not maxP)    && ~((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol))))...     &&  (PM(i,1) < (0 + tol)) ...          % x is zero     && ((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol))) % z is maxP         if kInt == 8             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is parallel to z axis at x = 0, y = 0     %///// checked if statement     if ((PM(i,3) > (0 + tol)) ...        % (z is not zero AND z is not maxP)    && ~((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol))))...      && (PM(i,1) < (0 + tol)) ...        % x is zero      && (PM(i,2) < (0 + tol))            % y is zero         if kInt == 1             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check is lattice site is parallel to z axis at y=0,x=maxP     if ((PM(i,3) > (0 + tol)) ...       % (z is not zero AND z is not maxP)   146    && ~((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol))))...     && ((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol))) ... % x is maxP      && (PM(i,2) < (0 + tol))                   % y is zero         if kInt == 2             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is parallel to z axis at x = maxP,y=maxP     if ((PM(i,3) > (0 + tol)) ...         % (z is not zero AND z is not maxP)    && ~((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol))))...     && ((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol)))...% AND x is maxP     && ((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol)))  % y is maxP         if kInt == 3             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     %/check if lattice site is parallel to z axis at x = 0, y = maxP     if ((PM(i,3) > (0 + tol)) ...          % (z is not zero AND z is not maxP)    && ~((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol))))...     &&  (PM(i,1) < (0 + tol))  ...         % AND x is zero     && ((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol))) % y is maxP         if kInt == 4             thetaExposed = 2*pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         else             thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         end     end     % Determine if lattice site is on a corner     if ((PM(i,1) < (0 + tol))...                 % (x is zero OR x is maxP)      ||((PM(i,1) < (maxP + tol)) && (PM(i,1) > (maxP - tol)))) ...      &&((PM(i,2) < (0 + tol)) ...                % AND (y is zero OR y is maxP)      ||((PM(i,2) < (maxP + tol)) && (PM(i,2) > (maxP - tol)))) ...      &&((PM(i,3) < (0 + tol)) ...               % AND (z is zero OR z is maxP)      ||((PM(i,3) < (maxP + tol)) && (PM(i,3) > (maxP - tol)))) ...         %This is a corner particle on zy plane         % Calculate the angle exposed to pore space         thetaExposed = pi - sum(ang_Int1) + (subt_ang1 + subt_ang2);         countCorner = countCorner + 1;     end     TPBL(j) = thetaExposed * HmInt; end  TPB = (sum(TPBL)*1e-6)/((maxP*1e-6)^3); % Three phase boundary length on a %volumetric basis Matlab Code for Calculating Total Porosity The total porosity is calculated by first determining the solid volume of the particles (including the sinter necks). The pore volume is the complement of the solid volume. The function used to calculate the solid volume is SolidVolume.m. VS is an array containing the solid volume of a sphere with the same diameter as the particle corresponding to the array   147 index. VSrev is the volume of the revolved sinter neck. VScap is the volume of the spherical cap that would otherwise be double counted since it intersects both the sinter neck and the neighbour particle, as described in Chapter Two, Figure 2-12. The function SolidVolume.m is not fully reproduced here, only the important differences compared to the three phase boundary length, described above. %************************************************************************ function SV = SolidVolume(PM, DM, DN, PN, rs, maxP, tol)  VS = zeros(1,length(PM)); VSrev = zeros(length(PM),12); VScap = zeros(length(PM),12); Vparticle = zeros(1,length(PM));  for i = 1:length(PM)  Not shown here are the checks to determine which particles are internal to the electrode and which particles are located on the face, edge or corner of the lattice because they are the same as those described above. The solid volume of particle i internal to the electrode is given by: VS(i) = 4/3*pi*(DM(i)/2)^3;  The solid volume of particle i on a face of the geometry envelope is given by: VS(i) = (1/2)*4/3*pi*(DM(i)/2)^3;  The solid volume of particle i on an edge of the geometry envelope is given by: VS(i) = (1/4)*4/3*pi*(DM(i)/2)^3;  The solid volume of particle i at a corner of the electrode is given by: VS(i) = (1/8)*4/3*pi*(DM(i)/2)^3;  The volume of the sinter neck is calculated as follows:     %Volume of sinter neck calculated from volume of revolution. This     %is done for each neighbor particle (each sinter neck)     for k = 1:12         kx = 3*k - 2;   % index of x coordinate in PM of neighbor particle         ky = 3*k - 1;   % index of y coordinate in PM of neighbor particle         kz = 3*k - 0;   % index of z coordinate in PM of neighbor particle          if PN(i,kx) < 0             % no neighbor here so advance to next neighbor position             continue         else          % Distance between centers of the particle and its neighbor: L   148          L = sqrt((PN(i,kx) - PM(i,1))^2 + ...                   (PN(i,ky) - PM(i,2))^2 + ...                   (PN(i,kz) - PM(i,3))^2);          rN = DN(i,k)/2;        % radius of neighbor particle          r = DM(i)/2;                          % radius of particle          % distance from particle center to plane of intersection with          % neighbor: SplaneInt          SplaneInt = abs(1/(2*L) * (L ^ 2 - rN ^ 2 + r ^ 2));          % Angle from line joining two centers counterclockwise to line          % joining center of particle to center of sinterneck radius          phi = acos(SplaneInt / (r + rs(k)));          % y coordinate of center of sinter neck radius of curvature: y0          y0 = tan(phi) * SplaneInt;          % x coordinate of center of sinter neck radius of curvature: x0          x0 = SplaneInt;          % distance to left boundary of volume of revolution: a          a = cos(phi) * r;          % distance to right boundary of volume of revolution: b          b = SplaneInt;          % equation of sinter neck radius: y          y = @(x)pi.*(-(rs(k).^2-x.^2+2*x.*x0-x0.^2).^(1/2)+y0).^2;          % Perform integration for volume of revolution          VSrev(i,k) = quad(y,a,b);%quad(y,a,b,1.0e-6);   The volume of the spherical cap is:          SphereCapHeight = r - a;          VScap(i,k) = 1/3*pi*SphereCapHeight^2*(3*r - SphereCapHeight);  Not shown here are the checks to determine if the spherical cap and sinter neck volume lie on an edge, face, or corner.         end     end      % Volume of particle, add the sinter neck volume and subtract the     % overlapping volume     Vparticle(i) = VS(i) + sum(VSrev(i,:)) - sum(VScap(i,:)); end  SV = sum(Vparticle); %Total solid volume in electrode model. %************************************************************************    149 Appendix C Computer Code for Calculation of Effective Diffusion Coefficient The following code has not been included in this appendix: • The code written by the author to create the cross section images of the modeled geometry and that to subsequently create the binary representations of these cross section images. • The code written by the author to remove the noisy voxels erroneously generated by Unigraphics. • The code written to convert the three phase boundary coordinates obtained using high resolution images (for example, 300 by 300 pixels) to their equivalent positions in a smaller matrix resulting from lower resolution images (for example, 76 by 76 pixels). • The code written by the author to filter out non-physical three phase boundary voxels that was added to the code written by the project collaborators. • The code written by the project collaborators to evaluate the connectivity of each of the three phases: electron conductor, ion conductor, and pore.  The calculation of the effective diffusion coefficient is described from the point where the open porosity matrix elements have been written to a data file porosity_open.out and the matrix elements for the connected three phase boundary voxels have been converted to their equivalent indices within a matrix of the same size as the open porosity, contained in the file tpbl_connected.out. The size of the open porosity matrix in this work is 76 by 76 by 76, as determined using the mesh refinement study described in Chapter Four. Matlab Code for Finite Element Mesh Generation The open porosity matrix indices are loaded from the file porosity_open.out, which contains four columns. The first is simply a counter used in the program code written by the project collaborators. Therefore, the open porosity elements are taken from the last three columns, containing the three indices cooresponding to the x, y, and z directions,   150 respectively. The coordinates of the open pore voxels start at 2 and range up to 77 as a result of the method of coding by the project collaborators. For example, a coordinate of 2 is in matrix element [1, 1, 1]. The physical length, char_length of the edge of one open pore voxel (finite element) is a user input. load 'porosity_open.out' porosity_open = porosity_open(:,2:end);  char_length = <input>;  The gas channel is represented by a layer of elements. p_open_add = porosity_open; X = 2:76+1; Y = 2:76+1; Nodesz = ones(1,(numel(X))^2); [XX YY] = meshgrid(X,Y); Nodesx = reshape(XX,1,numel(XX)); Nodesy = reshape(YY,1,numel(YY)); Nodes_add = [Nodesx' Nodesy' Nodesz'];  The node coordinates for node number one, using the convention in Figure 4-4, are generated as follows. All nodes, Nodes, in the initial finite element mesh are created as follows: p_open_add = [p_open_add' Nodes_add']'; Nodes1(:,1) = p_open_add(:,1) - 2; %#ok<COLND> Nodes1(:,2) = p_open_add(:,2) - 2; %#ok<COLND> Nodes1(:,3) = p_open_add(:,3) - 1; %#ok<COLND>  Nodes2 = [Nodes1(:,1)+1 Nodes1(:,2)     Nodes1(:,3)  ]; Nodes3 = [Nodes1(:,1)   Nodes1(:,2)+1   Nodes1(:,3)  ]; Nodes4 = [Nodes1(:,1)+1 Nodes1(:,2)+1   Nodes1(:,3)  ]; Nodes5 = [Nodes1(:,1)   Nodes1(:,2)     Nodes1(:,3)+1]; Nodes6 = [Nodes1(:,1)+1 Nodes1(:,2)     Nodes1(:,3)+1]; Nodes7 = [Nodes1(:,1)   Nodes1(:,2)+1   Nodes1(:,3)+1]; Nodes8 = [Nodes1(:,1)+1 Nodes1(:,2)+1   Nodes1(:,3)+1];  Nodes = [Nodes1' Nodes2' Nodes3' Nodes4' Nodes5' Nodes6' Nodes7' Nodes8'];  The duplicate nodes are removed: Nd = Nodes'; clear Nodes Nodes = unique(Nd,'rows');  The coordination matrix, vert_Ind,  is calculated as follows. This loop runs in parallel on four processors to speed up computation time. The coordination matrix has size np by eight, where np is the number of finite elements (open pore voxels). The columns contain the row   151 numbers of the node matrix for the eight nodes defining the element vertices. The coordination matrix therefore specifies which nodes define each element. elemCount = 1; %counter vert_Ind_Initial = zeros(numel(Nodes(:,1)),8); %preallocate matrix matlabpool open parfor numLoops = 1:numel(Nodes(:,1))     test_node = Nodes(numLoops,:);     Ind1 = find(Nodes(:,1) == (test_node(1) + 1) & ...                 Nodes(:,2) == (test_node(2) + 0) & ...                 Nodes(:,3) == (test_node(3) + 0));     Ind2 = find(Nodes(:,1) == (test_node(1) + 0) & ...                 Nodes(:,2) == (test_node(2) + 1) & ...                 Nodes(:,3) == (test_node(3) + 0));     Ind3 = find(Nodes(:,1) == (test_node(1) + 1) & ...                 Nodes(:,2) == (test_node(2) + 1) & ...                 Nodes(:,3) == (test_node(3) + 0));     Ind4 = find(Nodes(:,1) == (test_node(1) + 0) & ...                 Nodes(:,2) == (test_node(2) + 0) & ...                 Nodes(:,3) == (test_node(3) + 1));     Ind5 = find(Nodes(:,1) == (test_node(1) + 1) & ...                 Nodes(:,2) == (test_node(2) + 0) & ...                 Nodes(:,3) == (test_node(3) + 1));     Ind6 = find(Nodes(:,1) == (test_node(1) + 0) & ...                 Nodes(:,2) == (test_node(2) + 1) & ...                 Nodes(:,3) == (test_node(3) + 1));     Ind7 = find(Nodes(:,1) == (test_node(1) + 1) & ...                 Nodes(:,2) == (test_node(2) + 1) & ...                 Nodes(:,3) == (test_node(3) + 1));      % if all seven neighbor vertices comprising the hexahedral element     % exist then add to the vert_Ind_Initial matrix     if size(Ind1) ~= zeros(1,2) & size(Ind2) ~= zeros(1,2) & ...        size(Ind3) ~= zeros(1,2) & size(Ind4) ~= zeros(1,2) & ...        size(Ind5) ~= zeros(1,2) & size(Ind6) ~= zeros(1,2) & ...        size(Ind7) ~= zeros(1,2)        Ind = [Ind1 Ind2 Ind3 Ind4 Ind5 Ind6 Ind7];        vert_Ind_Initial(numLoops,:) = [numLoops Ind];        elemCount = elemCount + 1;     else         vert_Ind_Initial(numLoops,:) = zeros(1,8);     end end  matlabpool close  The nodes are converted to their physical coordinates, Nodes_um,  by multiplying the matrix index by the physical edge length of a voxel (finite element). Nodes_um = Nodes*char_length;  The preallocation of vert_Ind_Initial included more rows than there are elements, so the rows containing only zeros are removed. vert_Ind_unique = unique(vert_Ind_Initial,'rows'); vert_Ind_unique_sorted = sort(vert_Ind_unique,1,'ascend');   152 if vert_Ind_unique_sorted(1,:) == zeros(1,8)     vert_Ind_unique_sorted(1,:) = []; end clear vert_Ind_Initial; vert_Ind = vert_Ind_unique_sorted;  The node matrix, Nodes_um and coordination matrix, vert_Ind, are saved and later recalled when formulating the finite element problem. COMSOL / Matlab Code to Calculate Effective Diffusion Coefficient The effective diffusion coefficient is determined in an iterative manner. The initial value for the effective diffusion coefficient is estimated using the ratio of open porosity to tortuosity factor, using a value of three for the tortuosity factor. This intial estimate is used to estimate the current distribution through the electrode using the performance model developed by Gazzarri and Kesler, the code for which can be found elsewhere [51]. This performance model also contains the code to calculate the bulk diffusion coefficient, and is called using the function curr_dist.m. % Initialize effective diffusion coefficient D_eff_i = 1;  % Set D_O2N2_flag (1 for first iteration and 0 for remaining iterations) D_O2N2_flag = 1;  %/ Inputs load tpbl_total.out %tpbl_total.out contains matrix indices of all tpb voxels adresse_tpbl = tpbl_total(2:end); %coordinates are in the last three columns tpb = size(tpbl_total,1); %number of tpb voxels.  char_length = <input>; %voxel (finite element) edge length cell_length = 76*char_length; %Images have resolution of 76 by 76 pixels.  activite = input('input number of conntected tpb voxels: '); tpbl = input('input total three phase boundary length: '); porc = input('input open porosity: '); atpbl =(activite/tpb)*tpbl; %this is the connected tpb length.  tpbw = (0.03+0.07)/2*1e-6; %width of three phase boundary line  Dbin = curr_dist(D_O2N2_flag, D_eff_i, cell_length, atpbl, porc, tpbw);  At this stage, the bulk diffusion coefficient, Dbin, is known, and the distribution in Faradaic current has been written to a .mat file. The new effective diffusion coefficient, D_eff_new, is now calculated using the method described in Chapter Four by calling the function Deff.m. D_eff_new =Deff(Dbin,char_length,atpbl,tpbl,tpbw,tpb,D_O2N2_flag,cell_length);    153 The function Deff.m calculates the new effective diffusion coefficient as follows. *------------------------------------------------------------------------------- function D_eff_new = Deff(Dbin,char_length,atpbl,tpbl,tpbw,tpb,D_O2N2_flag,cell_length) atpba = tpbw*atpbl*(cell_length)^3; %connected three phase boundary area Upon first iteration of Deff.m, the finite element problem is constructed. Subsequent iterations only update the flux boundary conditions according to the updated distribution in Faradaic current obtained using the new diffusion coefficient in the performance model. if D_O2N2_flag     Nodes = load('UD2NodesCoord.mat','Nodes_um'); %load nodes     Nodes = Nodes.Nodes_um;     vert_Ind = load('UD2NodesCoord.mat','vert_Ind'); %load coordination matrix     vert_Ind = vert_Ind.vert_Ind;  COMSOL uses a structure, named femD here, to store the information relevant for the finite element problem. The mesh is stored in femD.mesh, which is initialized using linear hexahedral elements, ‘hex’, according to the nomenclature used by COMSOL. The built in COMSOL function meshenrich is required by COMSOL and essentially ‘cleans’ the input mesh. Specifically, it adds domain information relevant to geometry creation and fixes any large aspect ratio elements.      el1 = cell(1,0);     el1{1} = struct('type','hex','elem',vert_Ind');     m = femmesh(Nodes',el1);     %femD.shape = 1;     m = meshenrich(m);     femD.mesh = m;  femD.geom is a geometry object created by COMSOL from the input mesh in femD.mesh.     [femD,OP] = mesh2geom(femD, 'srcdata','mesh', 'destfield',{'draw'}, ...                      'srcfem',1, 'destfem',1,'drawtag','g2');     femD.geom = OP;  The node information p is the node coordinates and the element information el contains all information on the face elements, edge elements, and vertex elements, in addition to the hexahedral (brick) elements. el stores this information in a way that the geometry boundary numbers assigned by COMSOL can be obtained.     el = get(femD.mesh,'el');     p = get(femD.mesh,'p');    154 The boundary (or face) numbers assigned by COMSOL that are in a location containing a connected three phase boundary voxel are identified using the function TPB_Faces.m. Ind1 is node number one for the connected three phase boundary voxel, similar to that described in Chatper Four, Figure 4-4. fxy1 is an array containing the boundary numbers assigned by COMSOL that are parallel to the xy plane (equivalent to the face in Figure 4-4 containing nodes 1 through 4). The remaining faces, fxy2, fyz1, fyz2, fzx1, and fzx2 are the remaining five faces of the connected three phase boundary voxel. These six arrays contain the geometry boundary numbers corresponding to the faces on the connected tpb voxels.  [Ind1 fxy1 fyz1 fzx1 fxy2 fyz2 fzx2] = TPB_Faces(el, p, char_length, atpba, cell_length);  The function TPB_Faces.m is as follows: %***************************************************************************** function [Ind1 fxy1 fyz1 fzx1 fxy2 fyz2 fzx2] = TPB_Faces(el, p, char_length, atpba, cell_length);  p = p'; %transpose of element node coordinates  tol=char_length*1e-3; %tolerance used in searching through node coords.  load tpbl_connected.out %data file containing indices of connected tpb’s. Elem_Ind = tpbl_connected(:,2:end); size_elem = size(Elem_Ind);  The physical coordinates of the connected three phase boundary voxel vertices, Ind1 through Ind8, are calculated as follows: Ind1(:,1) = (Elem_Ind(:,1) – 1)*char_length; Ind1(:,2) = (Elem_Ind(:,2) – 1)*char_length; Ind1(:,3) = (Elem_Ind(:,3) – 0)*char_length;  Ind2 = [Ind1(:,1)+1 Ind1(:,2)     Ind1(:,3)  ]*char_length; Ind3 = [Ind1(:,1)   Ind1(:,2)+1   Ind1(:,3)  ]*char_length; Ind4 = [Ind1(:,1)+1 Ind1(:,2)+1   Ind1(:,3)  ]*char_length; Ind5 = [Ind1(:,1)   Ind1(:,2)     Ind1(:,3)+1]*char_length; Ind6 = [Ind1(:,1)+1 Ind1(:,2)     Ind1(:,3)+1]*char_length; Ind7 = [Ind1(:,1)   Ind1(:,2)+1   Ind1(:,3)+1]*char_length; Ind8 = [Ind1(:,1)+1 Ind1(:,2)+1   Ind1(:,3)+1]*char_length;  The following loop runs in parallel to speed up computation time. The xy1, xy2, yz1, yz2, zx1 and zx2 faces of each connected three phase boundary voxel that lie in the plane of a boundary face are identified and the boundary face number is returned. If the face is on the   155 outer geometry envelope parallel to the YZ and ZX planes, the face number is not returned as these faces are insulated. The function findface.m is described below. fxy1 = []; fxy2 = []; fzx1 = []; fzx2 = []; fyz1 = []; fyz2 = [];  matlabpool open %run loop in parallel on four processors (default number)  parfor ii = 1:size_elem(1)      % Determine if face xy1 exists      fc1 = findface(Ind1(ii,:), Ind2(ii,:), Ind3(ii,:), Ind4(ii,:), el, p);      if ~isempty(fc1) && ~(Ind1(ii,1) < char_length+tol)          fxy1 = [fxy1 fc1];      end      % Determine if face zx1 exists and is not on the outer envelope      fc2 = findface(Ind1(ii,:), Ind2(ii,:), Ind5(ii,:), Ind6(ii,:), el, p);      if ~isempty(fc2) && ~(Ind1(ii,2) < (0 + tol))          fzx1 = [fzx1 fc2];      end      % Determine if face yz1 exists      fc3 = findface(Ind1(ii,:), Ind3(ii,:), Ind5(ii,:), Ind7(ii,:),el, p);      if ~isempty(fc3) && ~(Ind1(ii,1) < (0 + tol))          fyz1 = [fyz1 fc3];      end      % Determine if face xy2 exists      fc4 = findface(Ind5(ii,:), Ind6(ii,:), Ind7(ii,:), Ind8(ii,:),el, p);      if ~isempty(fc4)          fxy2 = [fxy2 fc4];      end      % Determine if face zx2 exists      fc5 = findface(Ind3(ii,:), Ind4(ii,:), Ind7(ii,:), Ind8(ii,:),el, p);      if ~isempty(fc5) && ~(Ind3(ii,2) < (cell_length + tol) && ...                            Ind3(ii,2) > (cell_length - tol))          fzx2 = [fzx2 fc5];      end      % Determine if face yz2 exists      fc6 = findface(Ind2(ii,:), Ind4(ii,:), Ind6(ii,:), Ind8(ii,:),el, p);      if ~isempty(fc6) && ~(Ind2(ii,1) < (cell_length + tol) && ...                            Ind2(ii,1) > (cell_length - tol))          fyz2 = [fyz2 fc6];      end end matlabpool close %*****************************************************************************  The function findface.m returns the geometry face number that contains the same four nodes as the particular face on the connected three phase boundary voxel. %***************************************************************************** function fc = findface(Ind1, Ind2, Ind3, Ind4, el, p)  p = p'; %transpose of node coordinate matrix. %must identify which field of the structure ‘el’ contains the quadrilateral face %element information. for ii = 1:numel(el)     % 'quad' are the quadrilateral face elements     if strcmp(el{ii}.type,'quad') == 1   156         qel = ii;         q = el{qel}.elem;     end end  The variable q is a matrix with four rows (each row contains the row number of the matrix p corresponding to the x, y, and z coordinates of the face vertex). There are four vertices per face element so there are four rows. The number of columns of p corresponds to the number of boundary faces in the geometry. Therefore, the boundary number containing the connected three phase boundary voxel face is equal to the column number of q. xG = p(1,:); yG = p(2,:); zG = p(3,:); tol = 1e-9;   v =  (xG(q(1,:)) > Ind1(1) - tol & xG(q(1,:)) < Ind1(1) + tol ...      |  xG(q(2,:)) > Ind1(1) - tol & xG(q(2,:)) < Ind1(1) + tol ...      |  xG(q(3,:)) > Ind1(1) - tol & xG(q(3,:)) < Ind1(1) + tol ...      |  xG(q(4,:)) > Ind1(1) - tol & xG(q(4,:)) < Ind1(1) + tol) ...      & (yG(q(1,:)) > Ind1(2) - tol & yG(q(1,:)) < Ind1(2) + tol ...      |  yG(q(2,:)) > Ind1(2) - tol & yG(q(2,:)) < Ind1(2) + tol ...      |  yG(q(3,:)) > Ind1(2) - tol & yG(q(3,:)) < Ind1(2) + tol ...      |  yG(q(4,:)) > Ind1(2) - tol & yG(q(4,:)) < Ind1(2) + tol)...      & (xG(q(1,:)) > Ind2(1) - tol & xG(q(1,:)) < Ind2(1) + tol ...      |  xG(q(2,:)) > Ind2(1) - tol & xG(q(2,:)) < Ind2(1) + tol ...      |  xG(q(3,:)) > Ind2(1) - tol & xG(q(3,:)) < Ind2(1) + tol ...      |  xG(q(4,:)) > Ind2(1) - tol & xG(q(4,:)) < Ind2(1) + tol) ...      & (yG(q(1,:)) > Ind2(2) - tol & yG(q(1,:)) < Ind2(2) + tol ...      |  yG(q(2,:)) > Ind2(2) - tol & yG(q(2,:)) < Ind2(2) + tol ...      |  yG(q(3,:)) > Ind2(2) - tol & yG(q(3,:)) < Ind2(2) + tol ...      |  yG(q(4,:)) > Ind2(2) - tol & yG(q(4,:)) < Ind2(2) + tol) ...      & (xG(q(1,:)) > Ind3(1) - tol & xG(q(1,:)) < Ind3(1) + tol ...      |  xG(q(2,:)) > Ind3(1) - tol & xG(q(2,:)) < Ind3(1) + tol ...      |  xG(q(3,:)) > Ind3(1) - tol & xG(q(3,:)) < Ind3(1) + tol ...      |  xG(q(4,:)) > Ind3(1) - tol & xG(q(4,:)) < Ind3(1) + tol) ...      & (yG(q(1,:)) > Ind3(2) - tol & yG(q(1,:)) < Ind3(2) + tol ...      |  yG(q(2,:)) > Ind3(2) - tol & yG(q(2,:)) < Ind3(2) + tol ...      |  yG(q(3,:)) > Ind3(2) - tol & yG(q(3,:)) < Ind3(2) + tol ...      |  yG(q(4,:)) > Ind3(2) - tol & yG(q(4,:)) < Ind3(2) + tol) ...      & (xG(q(1,:)) > Ind4(1) - tol & xG(q(1,:)) < Ind4(1) + tol ...      |  xG(q(2,:)) > Ind4(1) - tol & xG(q(2,:)) < Ind4(1) + tol ...      |  xG(q(3,:)) > Ind4(1) - tol & xG(q(3,:)) < Ind4(1) + tol ...      |  xG(q(4,:)) > Ind4(1) - tol & xG(q(4,:)) < Ind4(1) + tol) ...      & (yG(q(1,:)) > Ind4(2) - tol & yG(q(1,:)) < Ind4(2) + tol ...      |  yG(q(2,:)) > Ind4(2) - tol & yG(q(2,:)) < Ind4(2) + tol ...      |  yG(q(3,:)) > Ind4(2) - tol & yG(q(3,:)) < Ind4(2) + tol ...      |  yG(q(4,:)) > Ind4(2) - tol & yG(q(4,:)) < Ind4(2) + tol) ...      & (zG(q(1,:)) > Ind1(3) - tol & zG(q(1,:)) < Ind1(3) + tol ...      |  zG(q(2,:)) > Ind1(3) - tol & zG(q(2,:)) < Ind1(3) + tol ...      |  zG(q(3,:)) > Ind1(3) - tol & zG(q(3,:)) < Ind1(3) + tol ...      |  zG(q(4,:)) > Ind1(3) - tol & zG(q(4,:)) < Ind1(3) + tol) ...      & (zG(q(1,:)) > Ind2(3) - tol & zG(q(1,:)) < Ind2(3) + tol ...      |  zG(q(2,:)) > Ind2(3) - tol & zG(q(2,:)) < Ind2(3) + tol ...      |  zG(q(3,:)) > Ind2(3) - tol & zG(q(3,:)) < Ind2(3) + tol ...      |  zG(q(4,:)) > Ind2(3) - tol & zG(q(4,:)) < Ind2(3) + tol) ...      & (zG(q(1,:)) > Ind3(3) - tol & zG(q(1,:)) < Ind3(3) + tol ...      |  zG(q(2,:)) > Ind3(3) - tol & zG(q(2,:)) < Ind3(3) + tol ...      |  zG(q(3,:)) > Ind3(3) - tol & zG(q(3,:)) < Ind3(3) + tol ...      |  zG(q(4,:)) > Ind3(3) - tol & zG(q(4,:)) < Ind3(3) + tol) ...   157      & (zG(q(1,:)) > Ind4(3) - tol & zG(q(1,:)) < Ind4(3) + tol ...      |  zG(q(2,:)) > Ind4(3) - tol & zG(q(2,:)) < Ind4(3) + tol ...      |  zG(q(3,:)) > Ind4(3) - tol & zG(q(3,:)) < Ind4(3) + tol ...      |  zG(q(4,:)) > Ind4(3) - tol & zG(q(4,:)) < Ind4(3) + tol);  if isempty(v) == 1      disp('face does not exist')  else     fc = unique(el{qel}.dom(find(v)));  end %*****************************************************************************  There may be more than one face element on a boundary face, so the unique face number is obtained with the command fc = unique(el{qel}.dom(find(v))), as shown above. At this stage, the face numbers that contain a connected three phase boundary have been identified. These are the faces that the Neumann boundary conditions are applied to. The face to apply the Dirichlet boundary condition for the concentration of oxygen at the gas channel is identified using findDirichlet_quad.m. The gas channel is always defined to be the face at z equal to zero. Dirichlet_BC = findDirichlet_quad(el,p);  %***************************************************************************** function Dirichlet_BC = findDirichlet_quad(el,p)  z = p(3,:); tol = 1e-9;  for ii = 1:numel(el)     % 'quad' are the quadrilateral face elements     if strcmp(el{ii}.type,'quad')         qel = ii;         q = el{qel}.elem;     end end  cz = z(q(1,:)) < tol & z(q(2,:)) < tol & z(q(3,:)) < tol & z(q(4,:)) < tol;  Dirichlet_BC = unique(el{qel}.dom(find(cz))); %*****************************************************************************  The distribution in Faradaic current, which has been converted to a flux using Faraday’s law, is imported from the matfile cs.mat. The flux is contained in a cubic spline, cs.     load cs.mat cs  The flux must be scaled because the connected three phase boundary area is not the same as the cumulative area of the faces that contain three phase boundary voxels. The scaling requires that the area of each boundary face in fxy1, fxy2, yyz1, fyz2, fzx1, and fzx2 be calculated. Also required is the z coordinate of the centroid of each boundary face. The   158 centroid and area of each face cannot be calculated until the COMSOL function meshextend is called. The meshextend function cannot be called until the finite element problem is fully defined. Therefore, a dummy flux is used for the interim.     %specify dummy fluxes     N_yz1 = fnval(cs,Ind1(1,3)*(1:numel(fyz1))); numel_yz1 = numel(N_yz1);     N_yz2 = fnval(cs,Ind1(1,3)*(1:numel(fyz2))); numel_yz2 = numel(N_yz2);     N_xy1 = fnval(cs,Ind1(1,3)*(1:numel(fxy1))); numel_xy1 = numel(N_xy1);     N_xy2 = fnval(cs,Ind1(1,3)*(1:numel(fxy2))); numel_xy2 = numel(N_xy2);     N_zx1 = fnval(cs,Ind1(1,3)*(1:numel(fzx1))); numel_zx1 = numel(N_zx1);     N_zx2 = fnval(cs,Ind1(1,3)*(1:numel(fzx2))); numel_zx2 = numel(N_zx2);     num_N_BC = sum([numel_yz1 numel_yz2 numel_xy1 numel_xy2  numel_zx1 ...          numel_zx2]);  The following code describes the finite element problem. % determine the locations that the insulatory ‘NO’, concentration ‘C’ and %flux ‘N” boundary conditions are to be applied. The number of  boundary %conditions is one ‘C’, num_N_BC ‘N’ and the remainder are insulatory.     for ii = 1:num_N_BC+2         if ii == 1             bnd_type{ii} = 'N0';         elseif ii == 2             bnd_type{ii} = 'C';         else             bnd_type{ii} = 'N';         end     end      % Determine face element information     for ii = 1:numel(el)         % 'tri' are the triangular face elements         if strcmp(el{ii}.type,'quad')             qel = ii;             q = el{qel}.elem;         end     end      % Number of boundary faces total     num_faces = numel(unique(el{qel}.dom));      % Values for Dirichlet Boundary conditions     bnd_c0 = zeros(1,num_N_BC+2);     bnd_c0(2) = 2.279; %concentration of oxygen at gas channel     bnd_c0 = num2cell(bnd_c0);      % Initial values for Neumann Boundary Conditions     bnd_N = num2cell([0 0 N_yz1 N_yz2 N_xy1 N_xy2 N_zx1 N_zx2]);      % Indices for all boundary conditions (which face to apply what BC)     bnd_ind = boundaryIndices(numel_yz1, numel_yz2, numel_xy1,          ...         numel_xy2,  numel_zx1, numel_zx2, fyz1, fyz2, ...         fxy1, fxy2, fzx1, fzx2, Dirichlet_BC, num_faces);  The function boundaryIndices.m determines which boundary faces get which boundary condition. The argument for bnd_ind is an array of the face numbers assigned by COMSOL.   159 %***************************************************************************** function bnd_ind = boundaryIndices(n1, n2, n3, n4,  n5, n6, f1, f2, ...     f3, f4, f5, f6, d1, nf)  bnd_ind = ones(1,nf);  bnd_ind(d1) = 2;  bnd_ind(f1) = 3:n1+2; bnd_ind(f2) = n1+2+1:n1+2+n2; bnd_ind(f3) = n1+2+n2+1:n1+2+n2+n3; bnd_ind(f4) = n1+2+n2+n3+1:n1+2+n2+n3+n4; bnd_ind(f5) = n1+2+n2+n3+n4+1:n1+2+n2+n3+n4+n5; bnd_ind(f6) = n1+2+n2+n3+n4+n5+1:n1+2+n2+n3+n4+n5+n6; %*****************************************************************************      % set up steady state diffusion application mode in COMSOL      clear appl     appl.mode.class = 'FlDiffusion';        appl.assignsuffix = '_di';     clear prop;                             prop.analysis='static';     appl.prop = prop;                       clear bnd     % Specify type of BC                    % Specify Dirichlet BC     bnd.type = bnd_type;                    bnd.c0 = bnd_c0;     % Specify Neumann BC                    % Specifiy face to apply BC     bnd.N = bnd_N;                          bnd.ind = bnd_ind;     appl.bnd = bnd;                         clear equ     % O2-N2 binary diffusivity              %On one subdomain (open porosity)     equ.D = Dbin;                           equ.ind = 1;     equ.init = bnd_c0(2);     appl.equ = equ;                         femD.appl{1} = appl;     femD.border = 1;                        femD.outform = 'general';      % Multiphysics     femD=multiphysics(femD);      % Extend mesh     femD.xmesh=meshextend(femD);  %Get initial value init = asseminit(femD);      % Update model     femD.sol = init;  The meshextend function has now been called. One of the tasks of this function is to update the mesh to quadratic elements. Now, the centroid and area of each boundary face on which a Neumann boundary condition is to be applied can be calculated using postprocessing commands.     %preallocated arrays for area (xy1_area) and centroids (z_c_xy1)     xy1_area = zeros(1,numel(fxy1)); z_c_xy1 = xy1_area;     xy2_area = zeros(1,numel(fxy2)); z_c_xy2 = xy2_area;     zx1_area = zeros(1,numel(fzx1)); z_c_zx1 = zx1_area;     zx2_area = zeros(1,numel(fzx2)); z_c_zx2 = zx2_area;     yz1_area = zeros(1,numel(fyz1)); z_c_yz1 = yz1_area;     yz2_area = zeros(1,numel(fyz2)); z_c_yz2 = yz2_area;   160      for ii = 1:size(xy1_area,2)  The area is obtained by integrating the value 1 over the face         xy1_area(ii) = postint(femD,'1','unit','m^2','dl',fxy1(ii),'edim',2);  The z coordinate of centroid is obtained by integrating z over the face and dividing by the area of the face.         z_c_xy1(ii)=postint(femD,'z','unit','m^3','dl',fxy1(ii),'edim',2)/...             xy1_area(ii);     end     for ii = 1:size(xy2_area,2)         xy2_area(ii) = postint(femD,'1','unit','m^2','dl',fxy2(ii),'edim',2);         z_c_xy2(ii)= postint(femD,'z','unit','m^3','dl',fxy2(ii),'edim',2)/...             xy2_area(ii);     end     for ii = 1:size(zx1_area,2)         zx1_area(ii) = postint(femD,'1','unit','m^2','dl',fzx1(ii),'edim',2);         z_c_zx1(ii)= postint(femD,'z','unit','m^3','dl',fzx1(ii),'edim',2)/...             zx1_area(ii);     end     for ii = 1:size(zx2_area,2)         zx2_area(ii) = postint(femD,'1','unit','m^2','dl',fzx2(ii),'edim',2);         z_c_zx2(ii)= postint(femD,'z','unit','m^3','dl',fzx2(ii),'edim',2)/...             zx2_area(ii);     end     for ii = 1:size(yz1_area,2)         yz1_area(ii) = postint(femD,'1','unit','m^2','dl',fyz1(ii),'edim',2);         z_c_yz1(ii)= postint(femD,'z','unit','m^3','dl',fyz1(ii),'edim',2)/...             yz1_area(ii);     end     for ii = 1:size(yz2_area,2)         yz2_area(ii) = postint(femD,'1','unit','m^2','dl',fyz2(ii),'edim',2);         z_c_yz2(ii)= postint(femD,'z','unit','m^3','dl',fyz2(ii),'edim',2)/...             yz2_area(ii);     end  In order to scale the flux at each layer of elements parallel to the xy plane, the number of connected three phase boundary voxels at each layer of elements parallel to the xy plane must be determined. numel_ctpb is an array with 76 elements (there are 76 layers of elements), and each index contains the number of connected three phase boundary voxels in the corresponding element layer parallel to the xy plane.     load tpbl_connected.out %load the coordinates of the tpb voxles     % determine number of connected three phase boundary voxels per layer     numel_ctpb = zeros(1,76);     for ii = 1:76        numel_ctpb(ii) = numel(find(tpbl_connected(:,3) == ii));     end      tol = 1e-8;     % need to group the faces. Each group contains faces with z coordinate   161     % of the centroid located within a particular element layer (1 <=z <= 76). a     % layer includes the top xy face but not the bottom xy face.     %Preallocate:     S_area = zeros(1,76);     S_xy1 = cell(1,76); gxy1 = S_xy1;     S_xy2 = cell(1,76); gxy2 = S_xy2;     S_zx1 = cell(1,76); gzx1 = S_zx1;     S_zx2 = cell(1,76); gzx2 = S_zx2;     S_yz1 = cell(1,76); gyz1 = S_yz1;     S_yz2 = cell(1,76); gyz2 = S_yz2;      for ii = 1:76   gxy1{ii} contains the index number of the centroid array of the boundary faces that have centroid located between (ii+1)*char_length + tol and (ii)*char_length + tol. i.e. the face has a centroid in the ii layer of elements parallel to the xy plane. S_xy1{ii} contains the total area of all xy1 faces having a centroid in the ii layer of elements parallel to the xy plane.         gxy1{ii} = find(z_c_xy1 < ((ii+1)*char_length + tol) & ...             z_c_xy1 > ((ii)*char_length + tol));         S_xy1{ii} = sum(xy1_area(gxy1{ii}));          gxy2{ii} = find(z_c_xy2 < ((ii+1)*char_length + tol) & ...             z_c_xy2 > ((ii)*char_length + tol));         S_xy2{ii} = sum(xy2_area(gxy2{ii}));          gzx1{ii} = find(z_c_zx1 < ((ii+1)*char_length + tol) & ...             z_c_zx1 > ((ii)*char_length + tol));         S_zx1{ii} = sum(zx1_area(gzx1{ii}));          gzx2{ii} = find(z_c_zx2 < ((ii+1)*char_length + tol) & ...             z_c_zx2 > ((ii)*char_length + tol));         S_zx2{ii} = sum(zx2_area(gzx2{ii}));          gyz1{ii} = find(z_c_yz1 < ((ii+1)*char_length + tol) & ...             z_c_yz1 > ((ii)*char_length + tol));         S_yz1{ii} = sum(yz1_area(gyz1{ii}));          gyz2{ii} = find(z_c_yz2 < ((ii+1)*char_length + tol) & ...             z_c_yz2 > ((ii)*char_length + tol));         S_yz2{ii} = sum(yz2_area(gyz2{ii}));   S_area(ii) contains the sum of the area of all faces in the (ii) layer of elements parallel to the xy plane.         S_area(ii) = sum([S_xy1{ii} S_xy2{ii} S_zx1{ii} ...                           S_zx2{ii} S_yz1{ii} S_yz2{ii}]);     end  The value of the Faradaic current density, i_f, at the midpoint of each layer of elements parallel to the xy plane is:   162     i_f = zeros(1,76);     for ii = 1:76;         %faradaic current density at midpoint of layer z         i_f(ii) = fnval(cs,0.5*(2*ii+1)*char_length)*4*96485;     end  The scaled Faradaic current density is therefore:     cnctd_area = numel_ctpb/tpb*(tpbl*tpbw*cell_length^3);     i_fs = i_f.*cnctd_area./S_area;  The scaled Faradaic current density is then converted to the molar oxygen flux using Faraday’s law.     % Faraday's law to obtain scaled flux     N_scale = i_fs/4/96485;  The dummy fluxes are updated:     clear N_yz1 N_yz2 N_zx1 N_zx2 N_xy1 N_xy2      N_yz1 = zeros(1,numel_yz1);     N_yz2 = zeros(1,numel_yz2);     N_xy1 = zeros(1,numel_xy1);     N_xy2 = zeros(1,numel_xy2);     N_zx1 = zeros(1,numel_zx1);     N_zx2 = zeros(1,numel_zx2);      %assign new flux values for each bounding surface     for ii = 1:76         N_yz1(gyz1{ii}) = N_scale(ii);         N_yz2(gyz2{ii}) = N_scale(ii);         N_zx1(gzx1{ii}) = N_scale(ii);         N_zx2(gzx2{ii}) = N_scale(ii);         N_xy1(gxy1{ii}) = N_scale(ii);         N_xy2(gxy2{ii}) = N_scale(ii);     end  Now the flux boundary conditions are updated to the structure femD.     clear bnd_N     bnd_N = num2cell([0 0 N_yz1 N_yz2 N_xy1 N_xy2 N_zx1 N_zx2]);      %/ update femD structure     femD.appl{1}.bnd.N = bnd_N;      % Multiphysics     femD=multiphysics(femD);      % Extend mesh     femD.xmesh=meshextend(femD);  The problem is now ready to be solved. But first the fem structure is written to .mat files so that upon the next iteration of this function, Deff.m, only the updated flux values need be computed.   163     f_appl = femD.appl{1};     f_brdr = femD.border;     f_outform = femD.outform;     fmesh = femD.mesh;     fgeom = femD.geom;     save femDappl.mat f_appl f_brdr f_outform     clear f_appl f_brdr f_outform save numelvars.mat numel_yz1 numel_yz2 numel_xy1 numel_xy2 numel_zx1 ... numel_zx2 numel_ctpb     save areas.mat cnctd_area S_area     save faces.mat gyz1 gyz2 gzx1 gzx2 gxy1 gxy2     save femDmesh.mat fmesh     save femDgeom.mat fgeom     clear fmesh fgeom  The next block of code is only executed on subsequent iterations of this function Deff.m. else     load femDappl.mat     load numelvars.mat     load areas.mat     load faces.mat     load femDmesh.mat     load femDgeom.mat     clear cs N_yz1 N_yz2 N_zx1 N_zx2 N_xy1 N_xy2     load cs.mat cs      femD.appl{1} = f_appl;     femD.border = f_brdr;     femD.outform = f_outform;     femD.mesh = fmesh;     femD.geom = fgeom;      i_f = zeros(1,76);     for ii = 1:76         %faradaic current density at midpoint of layer z         i_f(ii) = fnval(cs,0.5*(2*ii+1)*char_length)*4*96485;     end     % find the faces that have centroid at layer z (1 <= z <= 76)     cnctd_area = numel_ctpb/tpb*tpbl*tpbw*cell_length^3;     % calculate the scaled faradaic current density     i_fs = i_f.*cnctd_area./S_area;     % Faraday's law to obtain scaled flux     N_scale = i_fs/4/96485;      N_yz1 = zeros(1,numel_yz1);     N_yz2 = zeros(1,numel_yz2);     N_xy1 = zeros(1,numel_xy1);     N_xy2 = zeros(1,numel_xy2);     N_zx1 = zeros(1,numel_zx1);     N_zx2 = zeros(1,numel_zx2);      %assign new flux values for each bounding surface     for ii = 1:76         N_yz1(gyz1{ii}) = N_scale(ii);         N_yz2(gyz2{ii}) = N_scale(ii);         N_zx1(gzx1{ii}) = N_scale(ii);         N_zx2(gzx2{ii}) = N_scale(ii);         N_xy1(gxy1{ii}) = N_scale(ii);         N_xy2(gxy2{ii}) = N_scale(ii);   164     end      clear bnd_N     bnd_N = num2cell([0 0 N_yz1 N_yz2 N_xy1 N_xy2 N_zx1 N_zx2]);      %update femD structure     clear femD.appl{1}.bnd.N     femD.appl{1}.bnd.N = bnd_N;      % Multiphysics     femD=multiphysics(femD);      % Extend mesh     femD.xmesh=meshextend(femD);  end  The iterative conjugate gradients solver is run using the following command. femD.sol=femstatic(femD,'solcomp',{'c'},'outcomp',{'c'},'linsolver','cg',... 'prefun’,'ssor');  The postprocessing involves calculating the cross section area of the pore space, plane_area, at each element layer. The concentration is integrated over this area for each layer of elements and stored in c_conc. The z component of the diffusive flux is also integrated over this area and stored in diff_z_flux. z_inc is the z coordinate of the element layer parallel to the xy plane to be evaluated. This plane is specified by the matrix fem_dom, which contains two vectors within the plane, and the offset parallel to the plane. ndiv = 76; % number of locations to evaluate average concentration. % Determine locations along z axis to evaluate average concentration in yz % plane tol = char_length*0.001; minZ = char_length; z_inc = linspace(minZ, cell_length-tol, ndiv);  %initalize vectors c and plane_area c_conc = zeros(1,ndiv); plane_area = zeros(1,ndiv); diff_z_flux = zeros(1,ndiv);  for k = 1:ndiv     % 3 by 3 matrix; each column containing a point in the plane     fem_dom=[cell_length,0,0;0,cell_length,0; z_inc(k),z_inc(k),z_inc(k)];     pd_area = postcrossplot(femD,2,fem_dom,'surfdata','1', ...         'outtype', 'postdata', 'axislabel', {'x','y','z'});     % integrate 'c' over the plane specified by fem_dom     if isempty(pd_area.p)         plane_area(k:end) = 1;         break     end     plane_area(k) = meshintegrate(pd_area.p, pd_area.t, pd_area.d);   165      % interpolate concentration 'c' in specified plane: fem_dom     pd_conc = postcrossplot(femD,2,fem_dom,'surfdata','c', ...         'outtype', 'postdata', 'axislabel', {'x','y','z'});     % integrate 'c' over the plane specified by fem_dom     c_conc(k) = meshintegrate(pd_conc.p, pd_conc.t, pd_conc.d);     pd_diff_z_flux = postcrossplot(femD,2,fem_dom,'surfdata',...         'dflux_c_z_di', 'outtype', 'postdata', 'axislabel', ... {'x','y','z'});     diff_z_flux(k) = meshintegrate(pd_diff_z_flux.p, pd_diff_z_flux.t, ...         pd_diff_z_flux.d);  end  The concentration gradient is calculated according to; c_left = c_conc(1)/plane_area(1); c_right = c_conc(end)/plane_area(end); c_grad = (c_right - c_left)/((z_inc(end)+char_length) - z_inc(1));  The average superficial diffusive flux is: avg_diff_Nz = sum(diff_z_flux./cell_length^2)/numel(diff_z_flux);  The new effective diffusion coefficient is: D_eff_new = -avg_diff_Nz/c_grad; *------------------------------------------------------------------------------- The remainder of the script to carry out the iterations using the functions curr_dist.m and Deff.m is shown below. % Reset flag to continue iteration until converges D_O2N2_flag = 0;  % initialize error and set error tolerance error = abs(D_eff_new - D_eff_i)/D_eff_i; error_tol = 1e-2; disp('error is: '); disp(error) iter = 1; maxiter = 5;  disp('about to iterate until convergence or 5 iterations') while error >= error_tol && iter <= maxiter          D_eff_i = D_eff_new;          %Dbin = main1D(D_O2N2_flag, cell_length, D_eff_i);         Dbin = curr_dist(D_O2N2_flag, D_eff_i, cell_length, atpbl, porc, tpbw);          D_eff_new = Deff(Dbin, char_length, atpbl, tpbl, tpbw, tpb, D_O2N2_flag, cell_length);          error = abs(D_eff_new - D_eff_i)/D_eff_i;          iter = iter + 1 end

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066851/manifest

Comment

Related Items