Guidelines for the Rational Design and Engineering of 3D Manufactured Solid Oxide Fuel Cell Composite Electrodes

The growth of 3D printing has opened the scope for designing microstructures for solid oxide fuel cells (SOFCs) with improved power density and lifetime. This technique can introduce structural modiﬁcations at a scale larger than particle size but smaller than cell size, such as by inserting electrolyte pillars of ∼ 5–100 μ m. This study sets the minimum requirements for the rational design of 3D printed electrodes based on an electrochemical model and analytical solutions for functional layers with negligible electronic resistance and no mixed conduction. Results show that this structural modiﬁcation enhances the power density when the ratio k eff between effective conductivity and bulk conductivity of the ionic phase is smaller than 0.5. The maximum performance improvement is predicted as a function of k eff . A design study on a wide range of pillar shapes indicates that improvements are achieved by any structural modiﬁcation which provides ionic conduction up to a characteristic thickness ∼ 10–40 μ m without removing active volume at the electrolyte interface. The best performance is reached for thin ( < ∼ 2 μ m) and long ( > ∼ 80 μ m) pillars when the composite electrode is optimised for maximum three-phase boundary density, pointing toward the design of scaffolds with well-deﬁned geometry and fractal structures. what a benchmark ﬂat electrode can provide. In addition, an analytical approximation was proposed to predict the upper bound in performance improvement. All the results provide key insight for the rational design of electrodes through to a series of characteristic parameters, such as the characteristic active thickness of the coupled reaction/conduction process t ∗ ed , the effective conductivity factor of the ion-conducting phase in the composite elec- trode domain k eff , the ratio between pillar half-width and half-distance among pillar walls ω corresponding to the pillar volume fraction, and dimensionless factors (cid:4) W and (cid:4) h comparing pillar width and pillar thickness with the characteristic thickness.

Many studies in the literature have recognised the importance of the electrode microstructure in affecting the power density and the lifetime of solid oxide fuel cells (SOFCs). [1][2][3][4][5][6][7] Significant enhancements in power density have been achieved through the statistical control of the microstructural properties by tuning porosity, 8,9 volume fractions, 2,8,10 particle size [11][12][13] and even tortuosity of the phases. 14 However, nowadays the design of the electrode microstructure can be optimised with unprecedented precision through the advent of 3D printing and additive manufacturing techniques. [15][16][17] Additive manufacturing allows for the development of electrodes with 3D architectures, thus going beyond what conventional fabrication techniques, such as tape casting and screen printing, can produce. The exploitation of 3D additive manufacturing has been successfully applied by the battery community [18][19][20] and the implementation of this approach to SOFCs is currently under investigation. 15,21 One of the most promising applications consists of the insertion of ionconducting pillars connected to the electrolyte 15,[22][23][24][25] as reported in Figure 1. Such a modification can be regarded as the mesoscale structural control of electrode microstructure, where the term mesoscale here refers to a dimension larger than the typical constituent particle size (<1 μm) but smaller than the characteristic length of the cell, that is, falling in the order of 5-100 μm. 26 Ideally, 3D printing offers the opportunity to controllably manufacture pillars of any shape, with the desired geometric and spacing requirements, in order to provide a preferential pathway for ionic conduction. 22,25 It is important to note that this mesoscale modification differs from the corrugation of the electrolyte, 16,27 such as in the mono-block-layer built (MOLB) design, 26,28 where any enhancement in power density is due to the increase in the geometric electrode-electrolyte interfacial area per unit of planar projected area.
In the past, several experimental approaches have been adopted to produce a mesoscale structural modification as those reported in Figure 1b. Dai et al. 29 reported a simple method which enables the production of dimples by spray depositing a slurry consisting of electrolyte particles and 20 μm microspherical pore formers. The method is simple and inexpensive, but it allows only for a coarse statistical control of the structural modification by varying the pore former content in the slurry. A more accurate design of the electrode-electrolyte interface was reported by Konno et al., 24 who realized deep grooves (about 200 μm * 200 μm) via a blasting method into a thick electrolyte. This technique was improved by Nagato et al., 30,31 who used an excimer laser to produce trenches and pillars 5-40 μm wide. However, nowadays much wider possibilities in producing complex architectures have been opened by 3D printing as recently reported by Farandos et al. 15 who used inkjet 3D printing to grow cylindrical pillars with minimum feature resolution of 35 μm in plane and 1.2 μm in height.
These experimental activities have been followed by modelling studies. Konno et al. 24 reported a 2D numerical model, solving the mass and charge conservation equations in idealised structures with rectangular grooves. Their simulation results were in qualitative agreement with the experimental electrochemical data. A similar model was adopted by Delloro and Viviani 22 and by Geagea et al., 23 who focused on pillar shape and on the operational regimes for which the mesoscale structural control has beneficial effects on the electrode performance. More recently Shimura et al. 25 performed a three-dimensional numerical simulation to assess the effect of a pillar, virtually inserted into tomographic data, on the local nickel connectivity and three-phase boundary (TPB) distribution.
Despite these seminal studies, there is still a lack of fundamental understanding in this area. For example, a proper performance indicator to compare flat and structurally-modified electrodes has not yet been defined. Moreover, the fundamental reasons why the mesoscale structural modification i) enhances the power density, ii) for what microstructural architectures it is beneficial, and iii) how the shape and geometric characteristics of the pillar can be rationally designed have not been comprehensively addressed yet. This is particularly important because modelling is expected to guide the electrode design and to set the minimum requirements for the next generation of 3D printing and other manufacturing techniques. 15,16,21,27 In this paper, a mathematical model is proposed in the Modelling section to provide guidelines for the rational design of the mesoscale structural modification of SOFC electrodes. The model is kept as simple as possible to get analytical solutions whenever possible, which are helpful to predict the upper bounds of performance and relate them to the microstructural properties of the composite electrode. In addition, a comprehensive numerical study on the pillar shape and geometry provides indications on how to intelligently design the pillars (Results and discussion section). The conclusions of the study, along with the connection with practical applications, are summarised in the Conclusion section.

Modelling
The model takes into account only the most significant electrochemical phenomena occurring within the functional layer of a SOFC porous composite electrode. Consider a classic porous anode made of a homogeneous mixture of electron-conducting particles (e.g., Ni) and ion-conducting particles (e.g., yttria-stabilized or scandia-stabilized zirconia, YSZ and ScSZ, respectively) placed on the top of an electrolyte as in Figure 1a. Both the fuel and electronic current enter the anode from the top surface, moving through the pores and the connected electronic paths, respectively (note that the sign of the current is conventionally taken as opposite to the electron and ion flow). The electronic current is converted into ionic current at the three-phase boundaries (TPBs) spread within the electrode. 32,33 The ionic current is carried by the ion-conducting particles dispersed within the electrode and reaches the electrolyte, where current is then transferred to the cathode. No electrochemical reaction occurs within the dense electrolyte phase.
The electrochemical model, developed in the next section, mathematically describes the transport and conversion of current within the electrode and the electrolyte according to the following assumptions based upon well-established behaviors from the literature: r steady-state conditions; r uniform temperature; 22,26 r no mixed ionic-electronic conduction in each solid phase. In addition, the electronic conductivity of the electron-conducting phase is neglected, being several orders of magnitude larger than the ionic conductivity of the electrolyte phase (i.e., σ e >> σ O ); 22,34,35 r the charge-transfer reaction occurs at the percolating TPBs distributed within the composite electrode volume. The kinetic expression relating the current density to the activation overpotential is linearised 35,36 in order to obtain an analytical solution of the model for flat electrodes (see Analytical solution for a flat electrode section); r the gas transport in the stagnation layer and within the pores is neglected since the contribution of concentration losses is generally minor within the functional layer of electrodes. [35][36][37][38][39] For example, for the anodes considered in this study, this assumption holds for pores larger than 0.1 μm and porosity larger than 30%: in the worst case scenario, the current would be ∼2% lower when considering gas transport for a 20 μm-thick flat functional layer; r the microstructure of the composite electrode domain is not explicitly resolved and is regarded as a continuum, 22,34,40 characterized by its effective microstructural properties (i.e., the TPB density and the effective ionic conductivity factor) which are determined with appropriate microstructural models 12,41 or tomographic reconstruction; [42][43][44] r only a 2D cross section is simulated to reduce the computational requirements of the numerical model 22,24,26,35 although the mesoscale structural modification may have a three-dimensional geometry. This choice does not affect the general trends predicted by the model, so that the dimensionless output is readily transferable to 3D structures of the same nature (see Rectangular pillars section for the verification of this assumption).
Incidentally, these assumptions may make the model look similar to the fin model developed by Tanner et al. 35 However, in Tanner et al. model the charge-transfer reaction is assumed to take place at the electrolyte interface without any contribution from the composite electrode domain, which makes their fin model more appropriate to simulate infiltrated electrodes 45 rather than 3D manufactured electrodes.
Although model assumptions seem restrictive, in practice most of the electrodes commonly used in SOFC technology are consistent with these assumptions, such as Ni/YSZ and Ni/ScSZ anodes as well as LSM/YSZ cathodes (LSM stands for strontium-doped lanthanum manganite). A successful application of a similar modelling approach based on the same assumptions is reported by Bertei et al. 39 for SOFC composite anodes. Therefore, with these assumptions the model is representative of the main electrochemical phenomena that determine the electrode performance in functional layers while being simultaneously simple enough to allow for a fundamental insight and, whenever possible, useful analytical solutions. The ionic current density i O is calculated according to the Ohm's law: 33,34 where V O and σ O are the potential and ionic conductivity of the ionconducting phase, respectively. k eff represents the effective (or relative) conductivity factor of the ion-conducting phase within the composite domain. k eff is defined as the ratio between the effective conductivity of the ionic phase and its bulk conductivity, 9,46 that is, Sometimes k eff is defined as the ratio between the volume fraction and the tortuosity factor of the ion-conducting phase in the composite domain, 47,48 which is a definition equivalent to that adopted in this study. This dimensionless factor is always smaller than 1 and depends on the microstructural properties of the composite electrode, in particular on the volume fraction, connectivity, tortuosity and neck size of ion-conducting particles, 49 while being independent of the bulk conductivity σ O . Eq. 1 considers another microstructural parameter of the composite electrode domain, which is the TPB length density L TPB . The TPB density depends on porosity, composition and particle size of the composite domain. 9,12,41 The current density per unit of TPB length i TPB in Eq. 1 is calculated according to a linear electrochemical kinetic expression: 50 where the term in brackets represents the activation overpotential while F, R and T are the Faraday constant, the gas constant and the absolute temperature, respectively. i 0 T P B is the exchange current density per unit of TPB length and represents the kinetic constant of the charge-transfer reaction. More generally, Eq. 3 represents a linear local current-voltage relationship, which is often more appropriate for SOFCs than a non-linear Butler-Volmer kinetic expression. 51 The mathematical problem is closed by boundary conditions. Symmetric boundary conditions (i.e., no flux) are applied at the lateral boundaries. A no flux condition is applied at the top of the composite domain, indicating that the current enters the electrode in electronic form only. A Dirichlet boundary condition V O = -η tot is applied at the bottom of the computational domain, 22 representing the interface of the electrolyte with the cathode. Such a boundary condition sets the overpotential η tot applied to the system composed by electrode and electrolyte. This set of equations and boundary conditions therefore implicitly ensures that the electronic current is completely converted into ionic current within the composite electrode domain.
A numerical solution of Eqs. 1-3 is required for complex 2D geometries as those represented in Figure 1b. The total current density i tot converted within the electrode is numerically calculated by integrating for the ionic current density leaving from the bottom of the computational domain. 35 Note that since model equations are linear in V O , i tot scales linearly with η tot . Figure  1a), the ionic potential V O can be assumed uniform in each planar cross section 48,52,53 . In such a case, the mathematical problem represented by Eqs. 1-3 becomes one-dimensional, allowing for an analytical solution of the ionic potential along the thickness of the electrode as reported in other studies, for example, by Costamagna et al. 36,54 The total current density i f lat tot is determined by considering the electrode resistance R ed and the electrolyte resistance t ey /σ O in series as follows:

Analytical solution for a flat electrode.-In a flat electrode (
where t ey is the electrolyte thickness. By solving Eqs. 1-3, the electrode polarisation resistance R ed is equal to: where the minimum electrode resistance R min ed and the dimensionless Thiele modulus are calculated as follows: In Eq. 6b, t ed is the electrode thickness while t * ed is the characteristic thickness of the coupled reaction/conduction process. 55 Both the Thiele modulus and the characteristic thickness are typically defined in coupled reaction/diffusion problems in chemical engineering. 36,56,57 In particular, the characteristic thickness t * ed depends on transport properties (σ O ), kinetic properties (i 0 T P B ) and microstructural characteristics of the composite electrode (L TPB and k eff ), while being independent of the electrode thickness (note that the electrode thickness t ed can arbitrarily be smaller, equal or larger than t * ed ). The characteristic thickness t * ed can be regarded as the decay factor of the electrode resistance R ed as a function of electrode thickness t ed according to Eqs. 5 and 6b. More practically, t * ed is an indicator of the characteristic distance from the electrolyte interface required for electrochemical reaction and ionic conduction to take place in ideal conditions. For an ideal semi-infinite electrode (i.e., t ed → ∞), which provides the upper bound of total current, the current converted within the characteristic thickness t * ed is equal to the 63.2% (i.e., 1-e -1 ) of the total current converted in the electrode.

Average approximation for semi-infinite electrodes with mesoscale modification.-
The analytical solution for a flat electrode reported in Eqs. 4-6 can be applied to estimate the resistance in the limit case of semi-infinite electrode thickness even for electrodes with mesoscale structural modification.
The geometric ratio ω is defined as: where W is the half-width of pillars and w is the half-distance between pillar walls, see Figure 1c. In the limit of t ed → ∞, in 2D ω is equal to the volume fraction of pillars, while in 3D the volume fraction of pillars is equal to ω = W 2 /(W + w) 2 . In this condition the electrode can be regarded as a homogeneous continuum, with average microstructural properties, identified with a bar sign, evaluated as follows: Eq. 8a indicates that the average TPB density L T P B is weighted for the volume fraction of composite electrode domain 1-ω, that is, the volume fraction of electrode wherein the electrochemical reaction can take place. Similarly, the average effective conductivity factor of the ion-conducting phase k ef f is calculated according to the Wiener equation 58,59 , by averaging for the effective conductivity factor of the dense pillars, equal to 1, and the effective conductivity factor of the composite electrode domain, equal to k eff .
These average microstructural properties are used to calculate approximate values of the electrode resistance and current density with mesoscale structural modification in the limit of semi-infinite thickness, as follows: Eqs. 9a and 9b are formally equal to Eqs. 6a and 4, respectively, derived previously for a flat electrode, however with the use of the average microstructural properties L T P B and k ef f as defined in Eq. 8. Table I. Reference conditions considered in all the simulations. These conditions refer to the best Ni/ScSZ composite anode produced by Somalu et al., 10 that is, the anode that showed the lowest polarisation resistance. Microstructural and materialspecific parameters can be found in Bertei et al. 39 Parameter Value The values obtained from the application of Eqs. 7-9 are an approximate solution of the 2D model reported in the Model equations section. However, Eqs. 7-9 allow for a quick estimation of the current density by analytical means, providing also an accurate assessment of the upper bound of performance as discussed later.

Results and Discussion
Flat electrode: the benchmark.-It is worth analysing the electrochemical performance of a flat electrode because any improvement provided by a mesoscale structural modification must be evaluated in comparison with the current density that a benchmark flat electrode would exhibit in the same conditions.
In this study, unless otherwise specified, the reference conditions for all the simulations are those reported in Table I. These conditions refer to the best porous composite Ni/ScSZ anode produced and tested by Somalu et al., 10 which showed the lowest polarisation resistance in a series of flat electrodes fabricated with different Ni volume fractions and porosities. The choice of materials, i.e., Ni and ScSZ, is particularly indicated for operation at intermediate temperature, 60,61 which is set to 973 K. The thickness of the electrolyte, made of ScSZ, is set to 10 μm, which is typical for anode-supported planar SOFCs. 2,32,62,63 Material-specific and microstructural properties were evaluated by Bertei et al. 39 Although the absolute value of the results depends on the choice of this specific set of parameters, the emergent trends and the scaling laws discussed afterwards are general and highly applicable to a wide range of conditions and materials. Figure 2 shows the current density predicted by Eqs. 4-6 as a function of electrode thickness for a flat configuration. The current density monotonously increases as the thickness increases, approach-  Table I. ing an asymptotic value i f lat,asympt tot in the limit of semi-infinite electrode thickness (t ed → ∞). 35 The Figure also reports the characteristic thickness t * ed , equal to 20.4 μm. Note that in a hypothetical electrode showing the same minimum electrode resistance R min ed with halved k eff and doubled L TPB , the characteristic thickness t * ed is smaller by a factor 2 and equal to 10.2 μm, see the dashed line in Figure 2. Indeed, the minimum electrode resistance R min ed scales with 1/ L T P B k ef f while the characteristic thickness t * ed scales with k ef f /L T P B according to Eqs. 6a and 6c, respectively.
It is also interesting to note that, in the asymptotic limit, the electrolyte area specific resistance is t ey /σ O = 2.07 · 10 −6 m 2 while the electrode polarisation resistance is equal to R min ed = 2.35 · 10 −5 m 2 . The ratio of the two resistances, equal to the fuel cell Biot number Bi FC = (t ey /σ O )/R min ed as defined by Konno et al., 24 is equal to 0.088. This means that the dominant resistance is from the electrode, thus the whole cell performance would be significantly enhanced by focussing on the reduction of the electrode polarisation resistance, such as through a mesoscale structural modification.
Any improvements produced by a mesoscale structural modification of the electrode must be compared with a benchmark current density of the flat electrode evaluated in the same conditions. The asymptotic current density i f lat,asympt tot is a well-defined and unambiguous benchmark. Thus, the performance indicator proposed in this study to quantify the improvement due to mesoscale structural modification is the relative increase in current density δi tot , defined as follows: [10] where i tot is the total current density produced in the electrode with mesoscale modification. Note that defining δi tot in this way allows one to quantify any improvements against the maximum current density that a flat electrode can theoretically yield. It is worth mentioning that such an asymptotic value i f lat,asympt tot is approached if electronic ohmic losses and gas transport limitations are neglected, 9,35,36 as assumed in this study (see Modelling section). However, as the thickness increases, these assumptions are no longer met in real electrodes. This is not a problem for the practical application of Eq. 10: instead of the asymptotic current density i f lat,asympt tot , the 90% of its value might have been considered as a benchmark for the definition of the performance indicator δi tot . As shown in Figure  2, the 90% of i f lat,asympt tot is met for relatively thin electrodes, where both electronic ohmic losses and gas transport limitations are typically negligible. 36,64,37,39 Thus, the applicability of Eq. 10 can be extended by considering any desired benchmark value for the flat electrode, provided that such a benchmark is well-defined and kept constant. In this study, the asymptotic current density i f lat,asympt tot is used to quantify the relative improvement δi tot according to Eq. 10 because, in such a case, δi tot shows the improvement over the maximum current density that a flat electrode can ideally exhibit.
Rectangular pillars.-The first mesoscale structural modification analyzed in this study is the insertion of rectangular pillars as in the first sketch in Figure 1b. Rectangular pillars have been tested in a few experimental studies. 15,30 Figure 3a shows the computational domain in greater detail: t h is the pillar thickness, W is the half-width and w is the half-distance between pillar walls. Figure 3b shows the relative improvement in current density δi tot as a function of the geometric ratio ω defined in Eq. 7 for different values of W in the limit of semi-infinite electrode thickness (t ed → ∞) and semi-infinite pillar thickness (t h → ∞) as obtained from the numerical solution of Eqs. 1-3. The Figure shows that δi tot exhibits a parabolic behavior with a defined maximum as a function of ω. This is not surprising because as ω approaches 0 the electrode becomes flat, while as ω approaches 1 the composite electrode volume drops to zero and no electrochemical reaction occurs. The maximum improvement in current density over the flat electrode performance is 27% for ω ≈ 0.4. As W increases, the general improvement given by the mesoscale structural modification drops. This is explained by analysing the inset in Figure 3b, which shows that the ionic current is not evenly distributed within the pillar width for large values of W: since the inner region of the pillar is no longer optimally exploited for ionic conduction, δi tot decreases. Accordingly, the maximum of δi tot shifts to larger values of ω as W increases (see dashed line in Figure 3b) to compensate for the lower utilisation of the pillar width. For W > 60 μm, pillars are detrimental for performance as δi tot ≤ 0 in the whole range of ω, thus such structures should be avoided under the conditions used for this study.
A dimensionless factor W can be defined as: Such a dimensionless factor, which compares the pillar half-width with the characteristic thickness of the composite electrode domain, is useful to identify two regimes of ionic conduction within the pillar width. For W << 1 the whole width of the pillar is utilised for conduction and the ionic current density distribution is uniform, while for W >> 1 the inner region of the pillar is not optimally utilised. Interestingly, this means that the maximum improvement in current density is achieved for W << 1 while for W >> 1 the mesoscale structural modification can in fact be detrimental, in agreement with the numerical results in Figure 3b. This outcome is supported by published studies, for example by Nagato et al., 30 who proved experimentally that pillars with smaller width result in higher power density than wider pillars. In addition, similar predictions were reached by Tanner et al. 35 Interestingly, the dot-dashed line in Figure 3b shows that the average approximation reported in Eqs. 8 and 9 reproduces remarkably well the trend of δi tot as a function of ω for W ≤ 0.1. In addition, results of a 3D simulation in the limit of W → 0, reported with open diamonds, lie on the curve predicted with the analytical average approximation, confirming that 2D simulations are representative of the general trends of a three-dimensional structural modification. Thus, the upper bound of performance and the optimal ω (i.e., the optimal pillar spacing) can be analytically predicted by using Eqs. 8-9 as a function of electrode microstructural and material properties. This is a useful result for a rapid short-cut assessment and design of mesoscale structural modifications. Figure 4 shows the effect of the pillar thickness t h on δi tot in the limit of t ed → ∞ for two representative cases: W = 5 μm (i.e., W < 1, Figure 4a) and W = 30 μm (i.e., W > 1, Figure 4b). For any values of t h , δi tot shows a maximum as a function of ω, as already reported in Figure 3b in the limit of t h → ∞. More importantly, for both W < 1 and W > 1 the improvement over flat electrode δi tot decreases as the pillar thickness t h decreases. Similarly to Eq. 11,  another dimensionless factor can be defined as follows: [12] which compares the pillar thickness with the characteristic thickness of the electrode. Figure 4 shows that, for any value of W and ω, h must be much larger than 1 in order to maximise δi tot . When both W and h are considered together they help determine the maximum performance improvement caused by the mesoscale modification as discussed later. The effect of the electrode thickness t ed on δi tot is reported in Figure  5 for a representative case W = 5 μm, w = 8 μm and t h = 20 μm. Figure  5 shows that δi tot increases as t ed increases and approaches an asymptotic value in the limit t ed → ∞. Therefore, given the pillar geometric parameters, while a small electrode thickness can be even detrimental for performance as δi tot < 0 for t ed < 40 μm, there is no significant improvement in making the electrode thicker beyond t ed ≈ 4 · t h .
In summary, the results reported in Figures 3-5 lead to the following key outcomes informing design guidelines: r the mesoscale structural modification produces the largest improvement in performance in the limit of thin and long pillars, that is, for W = W/t * ed < 0.1 and h = t h /t * ed >> 1. Eventually, this matches the concept of scaffold electrodes, 65-70 that is, electrodes fabricated with a porous ionic backbone which is infiltrated with electronconducting particles or even a mixture of electron-conducting and ionconducting particles. In other words, the optimal 3D printed electrode with rectangular pillars results in the production of a scaffold with well-defined geometry; r the upper bound in performance can be analytically predicted through the approximated model reported in Eqs. 8 and 9. This represents an effective rapid method for the optimisation of the mesoscale structural modification.

Microstructural properties of the composite electrode domain.-
The microstructural properties of the composite electrode domain, k eff and L TPB , play a significant role in determining the electrochemical performance of a flat electrode, as already discussed above. Similarly, the electrode microstructural properties are crucial also when a mesoscale structural modification is applied. As demonstrated in the previous section, the upper bound of improvement is achieved in the limit of t ed = t h → ∞ and is accurately predicted by the average approximation described in Eqs. 8 and 9. Hence, the effect of the microstructure of the composite electrode domain on δi tot is addressed in this section through the average approximation model.
As already discussed in the previous section, the maximum improvement in performance δi tot occurs in the limit t ed = t h → ∞ for a specific value of ω, that is, for an optimal volume fraction of pillars ω opt . Such an optimal condition, which represents the upper bound of performance that the mesoscale structural modification can yield, can where R min ed is the minimum polarisation resistance of the flat electrode (see Eq. 6a) and R min ed,opt is the minimum polarisation resistance of the electrode with mesoscale modification corresponding to the upper bound of improvement in δi tot (i.e., calculated for ω = ω opt ). Note that Eq. 13 is valid only for k eff ≤ 0.5. For k eff > 0.5 R  Figure 6 shows ω opt and the ratio R min ed,opt /R min ed as a function of the effective conductivity factor of the ion-conducting phase in the composite electrode domain k eff . As k eff increases from 0 to 0.5, the optimal volume fraction of pillars ω opt decreases from 0.5 to 0 while the ratio R min ed,opt /R min ed increases from 0 to 1. Note that the lower the ratio R min ed,opt /R min ed , the larger the improvement δi tot over the flat electrode in accordance with Eqs. 9b and 10.
Several considerations can therefore be drawn from the analysis of Eq. 13 and Figure 6: r the optimal volume fraction of pillars ω opt depends only on k eff (see Eq. 13a). Other microstructural or material parameters, such as L TPB , σ O or i 0 T P B , do not play any role in determining the optimal volume fraction of the pillars under the current assumptions: the decision of introducing or not a mesoscale modification depends solely on the effective conductivity factor of the ion-conducting phase in the composite electrode domain; r the mesoscale structural modification is beneficial only for k eff ≤ 0.5, while it is always detrimental for k eff > 0.5, because the minimum resistance of the electrode with pillars R min ed,opt exceeds the resistance of the benchmark flat electrode R min ed (see Figure 6). Thus, for k eff > 0.5 the improvement of electrode performance must be achieved by different means other than mesoscale structural modification; r as k eff decreases, the optimal volume fraction of pillars ω opt increases (see Figure 6), meaning that pillars are required to compensate for the lower effective ionic conductivity of the composite electrode domain; r the mesoscale structural modification is particularly beneficial for electrodes with low effective ionic conductivity factor k eff , since the reduction in electrode resistance R   benchmark flat electrode R min ed , is relatively more significant (see Figure 6). However, assuming that the other parameters remain constant (in particular L TPB ), since the resistance of the flat electrode R min ed scales with 1/ k ef f (see Eq. 6a), a flat electrode with low effective ionic conductivity is not a good benchmark as it would show poor performance. Thus, a useful framework to improve performance is i) the microstructure of the flat electrode should be improved first, and then ii) an optimal mesoscale structural modification can be applied to further enhance the performance; r the improvement of the microstructure of the composite electrode domain envisaged in the previous point must be done rationally. Since R min ed scales with 1/ L T P B k e f f (see Eq. 6a), the microstructure of the composite electrode domain could be modified in order to reduce k eff and proportionally increase L TPB to keep the same R min ed . By doing so, the benefits of a mesoscale modification would be amplified, since R min ed,opt /R min ed decreases as k eff decreases (see Figure 6). This implies that a composite electrode intended for mesoscale structural modification should be optimised for maximum TPB density rather than for maximum k eff , because the effective ionic conduction is provided by the ion-conducting pillars. As an additional consequence, the microstructure of a composite electrode optimised for flat electrode configuration may not be the best one to be used as a basis for a mesoscale structural modification.
The last two points highlight the importance of a proper design of the microstructure of the composite electrode domain even when a mesoscale structural modification is adopted. In addition, the accurate quantification of the effective conductivity factor is of a paramount concern and appropriate methods, 71,72 based on the tomographic reconstruction of the microstructure, must be employed to avoid any incorrect evaluation.
Shape of the mesoscale structural modification.-In the previous sections, rectangular pillars were considered in the analysis. Although this is the pattern currently produced with different techniques, 15,24,30,31 additive manufacturing can potentially print pillars with unique and different shapes, 21 going beyond the capability of existing manufacturing methods such as screen printing and tape casting.
The introduction of a terrace of width W t and thickness t t onto a rectangular pillar, as reported in Figure 7a, is the first shape modification analyzed in this section. Figure 7b shows on the left axis the increase in current density beyond the benchmark flat electrode, as a function of the pillar thickness t h , for a pillar with terrace of dimensions W t = t t = 5 μm (δi terr tot , solid line). For comparison, the result for a rectangular pillar without terrace is reported (δi no terr tot , dashed line). The difference between the relative improvement given by the terrace beyond an equivalent rectangular pillar, i.e., δi terr tot − δi no terr tot , is reported on the right axis. The two insets show the distribution of current density converted per unit of TPB i TPB in the active volume for two cases, t h = 20 μm and t h = 40 μm. The insets in Figure 7b show that the electrochemical reaction takes place close to the electrode/electrolyte interface, being distributed mainly within a distance roughly equal to the characteristic thickness t * ed . This current distribution matches what expected in a flat electrode configuration. 36 The mesoscale structural modification makes the i TPB distribution distorted toward the pillar, allowing for more current to be converted in the proximity of the pillar up to a thickness slightly larger than t * ed . Such an effect, consistent with the work reported by Delloro and Viviani, 22 is due to the enhanced conduction that the pillar provides, thus opening a faster ionic conduction path toward the electrolyte.
When the terrace is placed at t h ≈ t * ed , such an extra lateral fin boosts additional current conversion, leading to an improvement in current density δi terr tot larger than that obtained without terrace δi no terr tot , as shown by the right axis in Figure 7b. On the other hand, placing the terrace too far above from the active region does not make any difference with the case with no terrace: δi terr tot − δi no terr tot decreases and approaches 0 as t h increases beyond t * ed in Figure 7b. Similarly, placing the terrace too close to the electrolyte interface (i.e., t h << t * ed ) produces detrimental effects as δi terr tot − δi no terr tot becomes negative, meaning that more current is produced when there is no terrace. This detrimental effect is due to the fact that the terrace removes electrochemically active volume from the composite electrode domain, thus placing the terrace too close to the electrolyte interface, exactly where more current i TPB is produced, reduces the active region wherein current is converted, leading globally to a lower performance. Figure 7c, which shows the difference δi terr tot − δi no terr tot as a function of the terrace thickness t t , provides an additional insight. As t t increases, the improvement in current density over the pillar without terrace increases and reaches as maximum, beyond which δi terr tot −δi no terr tot decreases and becomes negative. There is a simple reason behind this trend. For small t t values, the terrace is too thin to provide an effective conduction path to the main body of the pillar and to collect ionic current from the reacting zone in the composite electrode domain.
As t t increases, the terrace acts as a horizontal conducting fin, effectively providing ionic conduction. However, as the terrace thickness increases, more active volume is removed in the composite domain, especially closer to the electrolyte interface where electrode volume is required by the electrochemical reaction to take place. Thus, the beneficial effect of a more effective lateral ionic conduction is quickly counterbalanced by the removal of active volume in the composite domain.
Hence, the results of Figure 7 provide a clear insight into the role of the terrace and, more generally, on the pillar shape. The terrace is useful when it provides ionic conduction far from the electrolyte interface, especially at a distance comparable with the characteristic thickness t * ed , in order to collect more ionic current and thus provide a faster conduction pathway toward the electrolyte. Meanwhile, the mesoscale structural modification shall not remove active volume in the composite electrode domain, especially at the electrolyte interface where there is the maximum current conversion rate. In other words, the mesoscale structural modification must provide ionic conduction where needed, that is, farther from the electrolyte interface, without removing active electrode volume where needed, that is, at the electrolyte interface.
The proper balance of these two opposite requirements is a golden rule for the rational design of the pillar shape. In particular, maximising the ionic conduction farther from the electrolyte interface without removing active volume points toward a pillar shape which has to mimic the ohmic overpotential profile along the electrode thickness. In a flat electrode the ohmic overpotential profile is opposite to the activation overpotential distribution, 36 being 0 at the electrolyte interface and increasing exponentially along the electrode thickness. 36,54 Such an exponential profile is mathematically mimicked by the following equation: [14] which is used in Figure 8a to draw the outline of the pillar, resulting in a rounded terrace shape which resembles the ohmic overpotential distribution in a flat configuration. 36,54 In Eq. 14, t r is a roundness characteristic length: the smaller t r , the sharper the outline. Figure 8b shows the results of a design study on parameters W t and t r for W = 5 μm, w = 20 μm and t h = 25 μm. In particular, given a value W t , the right axis reports the value of t r corresponding to the maximum δi tot , which is reported on the left axis. As W t increases, δi tot increases until a broad maximum is reached. Notably, the values of δi tot obtained for an electrode with rounded terrace are larger than the maximum values provided by rectangular pillars with and without terrace for the same geometric parameters, although the improvement in δi tot is only minor, in the order of 1% or less. On the other hand, the roundness characteristic length t r continuously decreases, meaning that the optimal shape becomes sharper as W t increases as depicted in the insets of Figure 8b.
The results of the rounded pillar shape further confirm the role of the terrace: the mesoscale structural modification must provide conduction farther from the electrolyte without removing active volume at the electrolyte interface, leading to a rather sharp optimal shape which enhances the performance obtained by rectangular pillars. In accordance with this, the triangular shape depicted in the fourth sketch in Figure 1b is expected to be the worst design among the series analyzed, since it removes composite electrode volume at the electrolyte interface without providing any improvement through a conducting fin into the active region within a distance t * ed from the electrolyte. Such a rounded terrace shape is an ideal design paradigm, which can be practically mimicked by fabricating dimples according to the method proposed by Dai et al. 29 In practice, in a 3D structure the terraces can even touch each other (i.e., W t = w), provided that the third dimension z (i.e., the through-plane direction in Figure 8a) is exploited to provide gas and electronic current underneath the terrace. In addition, the ability to produce and the mechanical requirements of such a sharp shape are something which needs to be addressed, especially during operation and fabrication if additive manufacturing is to be employed. 15 While clearly of practical importance, these aspects are out of the scope of the current study. The main point of the design analysis on the pillar shape is the underlying rationale that should guide the optimisation of the mesoscale structural modification, which is producing a shape that is capable of providing ionic conduction farther from the electrolyte up to a thickness in the order of t * ed without removing active volume at the base of the pillar. Hierarchical and fractal-like structures resembling trees or human lungs can be foreseen as well, 73 provided the feasibility in producing mechanically stable dense features at the small resolution required by the coupled reaction/conduction process, which has a mesoscale characteristic length equal to t * ed .

Conclusions
In this study, the rational design of the electrode/electrolyte interface, here referred as mesoscale structural modification, in SOFC porous composite electrode functional layers was addressed through modelling means for the limit case of negligible electronic resistance and no mixed ionic-electronic conduction. A simple 2D numerical model, taking into account ionic conduction and electrochemical conversion of current, was developed and applied to simulate the enhancement in current density beyond what a benchmark flat electrode can provide. In addition, an analytical approximation was proposed to predict the upper bound in performance improvement. All the results provide key insight for the rational design of electrodes through to a series of characteristic parameters, such as the characteristic active thickness of the coupled reaction/conduction process t * ed , the effective conductivity factor of the ion-conducting phase in the composite electrode domain k eff , the ratio between pillar half-width and half-distance among pillar walls ω corresponding to the pillar volume fraction, and dimensionless factors W and h comparing pillar width and pillar thickness with the characteristic thickness.
The main conclusions of the study can be summarised as follows: r the fundamental reason why the mesoscale structural modification enhances the electrochemical performance stems from the enhanced ion-conducting pathway provided by pillars; r the mesoscale structural modification is beneficial only for k eff ≤ 0.5, while it is detrimental for k eff > 0.5 under the assumptions of this study. The optimal volume fraction of pillars ω opt , which depends on k eff only, as well as the upper bound in performance can be predicted analytically according to Eq. 13; r a composite electrode intended for mesoscale structural modification must be optimised for maximum TPB density rather than for maximum effective ionic conductivity; r the mesoscale structural modification produces the largest improvement in performance in the limit of thin and long pillars ( W < 0.1 and h << 1). This means that the optimal mesoscale modification tends toward the concept of scaffold electrodes with well-defined geometry; r the shape of the mesoscale structural modification should be chosen in order to provide ionic conduction farther from the electrolyte up to a distance of about t * ed without removing active volume at the electrolyte interface. This resembles the distribution of ohmic overpotential.
The understanding of the fundamental reason about when and why the mesoscale structural modification is beneficial, the correlation between performance improvement and microstructure of the composite electrode, the definition of characteristic parameters as well as the analytical expressions for the prediction of upper bounds are useful tools to guide the rational design of the mesoscale modification of SOFC electrodes and to set the minimum requirements for any fabrication technique adopted for such a scope, such as additive manufacturing. There are still some points that deserve further investigation in the future, such as the mesoscale structural design of electrodes with finite electronic resistivity, with mixed ionic-electronic conduction or with limiting gas transport losses and, more importantly, the mechanical requirements of the optimal shapes during electrode fabrication and operation.

Acknowledgments
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 654915 and the EPSRC grant EP/M014045/1. Helpful discussions with Dr. Samuel J. Cooper and Nick Farandos (Imperial College London) are gratefully acknowledged.
Disclaimer: the European Commission Research Executive Agency is not responsible for any use that may be made of the information this paper contains. ω ratio between pillar half-width and half-distance between pillar walls (Eq. 7) (-)

Superscripts
variables calculated according to the average approximation (Eqs. 8 and 9) asympt asymptotic value for semi-infinite thickness in flat electrode flat calculated in flat electrode max maximum no terr without terrace (i.e., rectangular pillar) terr with terrace Subscripts opt optimal conditions (i.e., upper bound of performance)