Theoretical analysis of the effect of isotropy on the effective diffusion coefficient in the porous and agglomerated phase of the electrodes of a PEMFC

In polymer membrane fuel cells (PEMFC), the pore microstructure and the effective diffusion coefficient (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff) of the catalytic layer have a significant impact on the overall performance of the fuel cell. In this work, numerical methods to simulate PEMFC catalytic layers were used to study the effect of isotropy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{xy}$$\end{document}Ixy) on the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff. The proposed methodology studies reconstructed systems by Simulated Annealing imaging with different surface fractions of microstructures composed by two diffusive phases: agglomerates and pores. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff is determined numerically by the Finite Volume Method solved for Fick's First Law of Diffusion. The results show that the proposed methodology can effectively quantify the effect of isotropy on the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff for both diffusion phases. Two trends were obtained in the magnitude of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff concerning the change in isotropy: (1) an analytical equation is proposed in this article for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff} \ge 5\% D_{0}$$\end{document}Deff≥5%D0 and (2) numerical solutions are determined for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff} < 5\% D_{0} .$$\end{document}Deff<5%D0. In our analytical equation are both a lineal and a logarithmic sweep. When the surface fraction is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\emptyset =$$\end{document}∅= 50%, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff decreases more linearly than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\emptyset = 10\%$$\end{document}∅=10% at the beginning of the isotropy change, which indicates that small changes in isotropy in the particulate material modify it drastically; under these conditions the diffusion coefficient in the pore is predominant. (3) When the surface fraction is less than 50%, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{eff}$$\end{document}Deff decreases more exponentially at the beginning and more linearly at the end of the isotropy change, which shows that small isotropy changes in the bar-aligned material drastically alter it. In this trend, diffusion in the agglomerate is less affected by isotropy. The proposed methodology can be used as a design tool to improve the mass transport in porous PEMFC electrodes.


Methods
The scope of this work is constrained to isotropic biphasic agglomerates, specifically those consisting of a randomly distributed phase (Pt/C catalyst) within a matrix phase (ionomer).The digital system is spatially defined by a mesh of control volumes (CV).These CV have recognizable phases and diffusivity coefficient values defined at each node.The SA method 35,37,38 is used to stir the initial system into a reference system sequentially.The materials are statistically characterized by intrinsic correlation functions in the SA reconstruction methodology.Because of this, it is possible to determine the isotropy ( I xy ) at a specific SA stirring interval, defined by the iteration number of the SA method ( I SA ).The D eff is established by the FVM and is relates to I xy of the current system.A block diagram of the methodology proposed in this research is presented in Fig. 1.
The source system (initial microstructure) is conditioned to random bars ordered from south to north in the sense of the dominant direction of the mass diffusion.The target system (reference or final microstructure) is conditioned to a dispersed system randomly distributed that is a homogeneous, periodic, and isotropic CL.Our initial and reference microstructures consider only one filter: the one-point correlation function (pore fraction).

Simulated annealing
The SA methodology is an optimization method used to determine local or global optima in complex problems 39,40 .In the iteration of the SA method, a pixel of a different phase is exchanged, and the exchange is evaluated and accepted based on a cooling probability based on a fictitious temperature.This approach transforms the original system into a target system.This technique has been extensively employed in reconstructing stochastic heterogeneous materials (SHM) [32][33][34][41][42][43] , demonstrating its versatility and effectiveness. Howevr, it is essential to define some fundamental bases.
The statistical analysis considers the generation of an ensemble (Ω) composed of W realizations ꞷ; in this work, W = 15.Statistically, the realization ꞷ represents the domain V of the studied system (CL) within the measurable space ℜ d , partitioned by a j phase.In this context, the region in the space V j (ω), the surface fraction φ j (ω) , and the index function I j (ω) can be defined.
(1) The realization ω can be interpreted as a CL of the set Ω of CL´s.All these are microstructures with a finite number of nodes on the x-axis ( Nx ) and the y-axis ( Ny ); the specific phase of each node can be identified.This work studies two-phase isometric materials in a square domain of Nx = Ny pixels.In this way, it is possible to represent the agglomerate matrix phase (white pixels) and the pore dispersion phase (black pixels).When examining the white phase, the index function 1 can be rewritten as: Once the domain has been defined, applying different correlation functions to obtain statistical information to characterize the material is possible.In this work, the two-point correlation functions (S (2) j (r) ) and the linear path correlation function ( P j ) were employed in both phases of SA reconstruction.
Considering a material with domain V ⊆ ℜ d , volume V, and composed of j = 2 disjoint random phases.The S (2) j (r) function is the probability that the initial point and the endpoint of a line of length r fall in the same phase j.Considering an isotropic and homogeneous medium.The S (2) j (r) can be defined as a function of distance r 28 : where 〈〉 refers to the statistical average from evaluating the whole domain ꞷ.I J is the index function of the computational domain (Eqs. 1 or 2).The function S (2) j (r) provides a measure of how the endpoints of a vector r in the studied phase are correlated.On the other hand, the P j function is defined by: Both correlation functions are used for the SA reconstruction; they are considered for both phases ( j 0 and j 1 ).Equation ( 5) defines the error between the source system and the current system ( E SA ): where, F(r) characterizes the source system, and F ′ (r) characterizes the SA current system.The pixel phase shift is accepted if the probability P(�E SA ) is satisfied, In this work, the SA methodology is applied pixel by pixel, simplifying the theoretical study of the evolution of an initially bar-ordered system as its transitions to a random dispersion system.It should be noted that the user defines both the initial system (ordered bars) and the final system (random dispersion).

Isotropy
The isotropy of a material can be expressed as a function of the orientation of the tensor describing the transport phenomenon 44 or as a geometric constraint influencing the covariance matrix in a random field 45 .The relationship between the sum of covariance models in factor spaces and zonal isotropy and the relationship between geometric isotropy and an isotropic covariance function are relevant aspects in this context 46 .It is important to note that geometric and zonal isotropy are interrelated.
The concept of geometric isotropy is a valuable tool for describing the behavior of microstructures based on both accurate models, such as images acquired through scanning electron microscopy and synthetic CLs.This approach for describing microstructures is applicable to CL structures of varying dimensions.However, one of the fundamental challenges is the selection of appropriate scaling factors, which is a significant hurdle to overcome in the field of research.In this study, we will use the S  where, E xy is the mean square error of the correlation functions.Since the user defines the initial and final microstructure of the study, it is feasible to normalize the difference for a generalized comparison.The proposed isotropy (I xy ) is shown in Eq. ( 9).
(2) I j (ω) = 1, if the position pixel x, y is white agglomerate matrix .0, otherwise it is black pore dispertion .
(3) S (2) www.nature.com/scientificreports/E xy,initial corresponds to the source system of ordered bars (higher anisotropy) and E xy,actual is the mean of the difference between squares E xy of the actual SA iteration system.The bar microstructure is represented by I xy = 0 , as SA shakes the system I xy tends to 1 38 .

Effective diffusion coefficient
The solution of the FVM is iterative and requires the implementation of the physical model in the nodal point 46,47 .The molar flux due to diffusion can be calculated from the concentration gradient and a diffusivity coefficient, as are indicated by the first Fick's law described below 48 : where J F is the flux density of the species diffusing through a unit area, D is the diffusion coefficient, C is the differential of the species concentration, and x is the diffusion length.Considering the general equation of transport phenomena described in classical fluid mechanics and that the study is limited to a stationary diffusion state, without a convective and transient phenomenon, as well as without generation or consumption of species, the diffusion in a volume can be expressed as in Eq. ( 11): Using the Gaussian divergence theorem, we have that the integration over the CV is equal to the integration over the surface bounding that CV.For two-dimension mesh, the diffusion can be expressed as follows: where A is the CV cross-sectional area, ∂x is the distance of the volume or control element between boundaries, C is the species concentration, D the phase diffusivity coefficient as a function of the flow direction, the subindices e , w, s, and n are the east, west, south, and north direction, respectively.As well as P is the VC central position.The tri-diagonal matrix algorithm (TDMA) is used to determine the pattern of the complete spatial concentration distribution ( C x,y ).
In this study, we define the effective flux ( J eff ) as the sum of the local fluxes across horizontal mesh layers.We then utilize the numerical J eff value and modify the proportionality coefficient in the fundamental Eq. ( 10), but renaming it to D eff , as it is shown in Eq. ( 13), The material properties of the CL mesh define, in a deterministic way, the ETC ( D eff ). Figure 2 shows FVM boundary conditions for J eff numerical determination.
In this work, we simulate both pore (empty medium) and agglomerate phases (Pt/C + Nafion).The diffusion coefficient of oxygen in an empty medium is D O2,Pore = 3.23 × 10 −5 m 2 s −1 and the diffusion coefficient of (10) oxygen in the agglomerate is the same that in the ionomer D O2,Agg = 8.45 × 10 −9 m 2 s −149 .The exact analytical solution for a parallel diffusion resistance is given by the Eq. ( 14).
where D 0 is the exact value of the diffusion coefficient when the both phases are vertically aligned ( I xy = 0.0).

Results
The process of reconstruction of the catalytic layers starts with the generation of the synthetic microstructures comprised by the pore phase (black color phase) and the agglomerate (white color phase), the microstructures are bars aligned vertically from north to south and have the following microstructural characteristics: surface fraction ( φ ) of 50% pore-50% agglomerate, 40% pore-60% agglomerate, 30% pore-70% agglomerate, 20% pore-80% agglomerate and 10% pore-90% agglomerate.The domain of all samples is 350 × 350 pixels where each pixel is equivalent to 0.1 µm, so the samples are 35 µm × 35 µm.
The SA method is employed to obtain an ensemble composed of W = 15 random repetitions for each studied CL.The results present the CL isotropy and the D eff determined for both pore and agglomerate phases.Figure 3 shows representative images of the studied CL microstructures.The first column shows porosity, and the C1 column shows the aligned bars when I SA = 0 iteration, C2 column is for I SA = 5000 iteration, C3 column is for I SA = 20000 iteration, and the C4 column shows the dispersed random material (reconstructed image).
It is important to highlight that all samples' SA reconstruction temperature was the same ( T = 1x10 −6 ) and the reconstructions were stirring until the convergence of the annealing error E SA = 1x10 −6 .The C1 column indicates the bar-ordered CL system, I xy = 0.00 ± 0.00 , while the C4 column represents the homogeneously agitated CL system I xy = 1.00 ± 0.00 .C2 and C3 columns were selected to analyze their behavior in the results shown in Fig. 6.Each image in Fig. 3 represents a Ω ensemble, but the values are the averages and standard deviation obtained for W = 15.
Figure 4 presents plots based on I SA iterations: (A) for the error variation of the reconstruction process E SA (Eq.5), (B) for the isotropy variation I xy (Eq. 9)and (C) for D eff (Eq.13), all graphs based on the I SA statistical iteration.The values presented in this figure include all the universe of ω realizations of the total samples studied.
Figure 4A shows that the E SA value tends to zero, which indicates that the actual microstructure of the CL is converging to the values of the objective microstructure of the CL, as observed in the representative images of Fig. 3.As the I SA iterations increase, the overall error between the reconstructed microstructure and the target microstructure decreases, the reconstruction can be seen in Fig. 3, where in column 1 we have the initial microstructure arranged in aligned bars and in the following columns 2 to 4 the dispersion of the pore phase and the agglomerated phase is observed.
The I xy values in Fig. 4B tend to unity, indicating that the error between the correlation functions S j,x (r) and S j,y (r) are decreasing, as indicated in Eq. ( 9).It is observed that I xy is more dependent of Φ than E SA , this trend is similar to that reported by 38 .Figure 4C shows the response of D eff coefficient depending on the reconstruction process iteration I SA .It is observed that the D eff variation is higher for the first I SA iterations.
Figure 5 shows the response of D eff compared to the evolution of isotropy I xy .It must be pointed out that the numerical solution considers the black color as the pore phase and the white color as the agglomerate phase, both for a synthesized microstructure of a catalytic layer.The values presented in this figure includes all  www.nature.com/scientificreports/possibilities (points), average values (solid lines) and the curves of D eff after solving the new proposed Eq. ( 15) (black dotted lines).
Figure 5A shows the evolution of D eff as a function of isotropy I xy .The results from all examined samples ( � ∈ W = 15ω ) are presented with markers and solid lines denoting their averages.It can be observed that when the I xy = 0 the microstructure corresponds to the pore phase and the agglomerate are vertically aligned.In this condition, it is observed that D eff values are equal to the analytical solution of a parallel diffusion resist- ances; these values their equations are highlighted in the gray shading of Fig. 5A.The values of D 0 obtained by exact analytical solution of Eq. ( 14) are: D 0 (∅ = 50) = 1.6 × 10 −5 m 2 s −1 @ (50%D O2,Pore -50%D O2,Agg ) and D 0 (∅ = 10) = 3.2 × 10 −6 m 2 s −1 @ (10%D O2,Pore -90%D O2,Agg ).The theoretically based Eq. ( 15) is proposed for the range D eff ≥ 5%D 0 .This range is selected by the authors to reduce the variance between the numerical average values and the analytical solution values.
where D 0 is the exact solution of Eq. ( 14), I xy 0 is the lowest numerical I xy value different to cero (in this work I xy 0 = 5.40 × 10 −3 ) and m 1 y m 2 are analytical slopes obtained from numerical data, related by the porosity ∅ in Eqs. ( 16) and ( 17). ( 15) Table 1 shows the solution of Eqs.(14, 16 and 17).These values are applied in Eq. ( 15) to obtain the black dotted lines in Fig. 5A.
Ceballos et al. 30 reports effective oxygen diffusion values for a catalytic layer of a PEMFC with 30% pore phase surface fraction; the results reported in their work are obtained by solving different effective diffusion models, they report oxygen D eff in a range of 1.87 × 10 −6 m 2 s −1 for Nam and Kabiany model and 8.91 × 10 −6 m 2 s −1 for Tomadakis and Storirchos model; In this work, for the same pore phase surface fraction condition an analytical result of 9.7 × 10 −6 m 2 s −1 is obtained.The variations in the reported values are attributed to the different microstructures of the evaluated electrodes.
As isotropy increases, D eff values exhibit an exponential decrease, the exponential decrease being particu- larly pronounced in the sample φ = 10% (green color) and less pronounced in the sample φ = 50% (red color).Figure 5B shows a zoom section for details of the range D eff ≤ 1.5 × 10 −7 m 2 s −1 and I xy ≥ 0.3 .In this figure,   50 report experimentally obtained magnitudes of the effective diffusivity coefficient for catalytic layers.The reported averages of D eff for catalytic layers range from 4.2 ± 0.9 × 10 −7 m 2 s −1 @ 25 °C to 4.6 ± 0.5 × 10 −7 m 2 s −1 @ 75 °C.Experimental samples exhibited a surface fraction of the pore phase within a range of φ = 30% to φ = 40%.In this study, the reported values for different surface fractions of the pore phase and agglomerates are in the range of magnitudes of the order of 10 −6 when I xy is equal to zero and 10 −8 when I xy is equal to one.Our numerical model takes into account the diffusive coefficient of oxygen within the agglomerate, attributing that as the pore phase decreases the D eff values approach the reported value of D O2,Agg = 8.45 × 10 −9 m 2 s −1 .Similarly, Shen et al. 51 report experimentally obtained magnitudes of D eff for a catalytic layer with different thicknesses; the average value obtained is 1.47 ± 0.05 × 10 −7 m 2 s −1 , a figure that also falls within the range of results obtained in this study.Chen et al. 52 presents results of the effective diffusivity coefficient of oxygen for a sample with a surface fraction of 17% for the pore phase using the lattice Boltzmann method, reporting a D eff of 7.71 × 10 −8 m 2 s −1 .In our work, the result obtained for a similar surface fraction of the pore phase is 1.4 × 10 −8 m 2 s −1 when it is 20% D O2,Pore -80% D O2,Agg .The variation in results is attributed to the difference in the isotropy value between synthetic and experimental materials.Microstructural isotropy directly influences the behavior of the effective diffusivity coefficient.The results show that the exponential tendency is preserved at higher isotropy ranges.The box highlighted in the gray shading of Fig. 5B presents D eff @I xy = 1.
As evidenced by the results, the impact of oxygen diffusion within the agglomerate exhibits insignificance concerning the effective diffusion coefficient, D eff , particularly when the surface fraction of the porous phase is predominant.Conversely, in scenarios where the surface fraction of the agglomerate phase predominates, it becomes apparent that the diffusion within the agglomerate significantly influences the ultimate determination of the diffusion coefficient.This observation is elucidated in the zoom section in Fig. 5.
Figure 6 shows D eff as a function of pore surface fraction φ , for three selected isotropies referred to in columns of Figs. 3 and 5: the solid black line for I xy = 0.0 , the gay line for I xy = 1.0 and dashed black line for I xy = 0.5.
Figure 6 shows the results obtained for D eff as a function of the superficial fraction of the agglomerate phase.The results from all examined samples ( � ∈ W = 15ω ) are presented with markers and solid lines denoting their averages.For I xy = 0.0 and I xy = 0.5 the D eff ranges are on the left side of the graph; for I xy = 1.0 the D eff ranges are on the right side of the graph.
The results are presented within φ = 10% to φ = 50%, with 10% intervals.The solid black line, labeled as I xy = 0.0, corresponds to the average of all samples when the microstructure is vertically aligned.Numerical results indicate that the effective diffusion coefficient is directly related to the porosity of the catalyst layer microstructure.Higher porosity has a positive impact on D eff enhancement, this trend agrees with empirical correlation models of the effective diffusion coefficient in porous materials (Bruggeman 53 , Neale and Nader 54 , Tomadakis and Sotirchos 55 , Mezedur et al. 56 , Zamel et al. 57 , Das et al. 58 ).The results suggest that the isotropy of the microstructure also have a significant effect on the D eff as shown in Fig. 6.An increasing trend is observed when the porosity of the microstructure is more significant.When the porosity is minor, the diffusion of the agglomerate becomes more relevant in the effective diffusion coefficient, because the higher the surface fraction of the agglomerate, the higher the value of the predominant oxygen diffusion is that of the diffusion in the ionomer, which can be seen in Fig. 6 50,59 D eff values range from 3.26 × 10 -6 m 2 s −1 @ φ = 10% to 1.60 × 10 -5 m 2 s −1 @ φ = 50%.This behavior is attributed to the reduced diffusive resistance when the phases are vertically aligned.As isotropy in the microstructure increases, D eff diminishes due to the agitation within the microstructure.When I xy = 0.5, D eff values are within the range of 3.74 × 10 -8 m 2 s −1 @ φ = 10% and 6.53 × 10 -6 m 2 s −1 @ φ = 50%.In the case of a randomly monodisperse microstructure ( I xy = 1.0), the value of D eff varies between 1.08 × 10 -8 m 2 s −1 @ φ = 10% and 6.22 × 10 -8 m 2 s −1 @ φ = 50%.Accordingly, microstructures with lower isotropy (aligned bars) and higher porosity have lower diffusion resistance, which translates as higher D eff .

Conclusions
In this study, a theorical investigation has been carried out to determine the effect of the isotropy on the effective diffusion coefficient of PEMFC catalyst layers that change from aligned to dispersed pattern microstructure.The CL´s studied in this work are composed of two diffusive phases, the porosity, and agglomerated phases.We apply two advanced numerical methodologies: (1) Finite volume method to define a specific oxygen diffusion coefficient into each mesh node: D O2,Pore for pores and D O2,Agg for agglomerates. (2)Simulated Annealing to simulate the evolution of unidirectional aligned bars to random mono-dispersed CL´s.
The analysis confirms that the aligned bar is the better structure for improving the CL mass transport phenomena when the diffusion is made in the parallel direction to the bars because there is less resistance to the transport of the oxygen molecule.Further, this statistical strategy allows us to determine the effect of geometric isotropy ( I xy ) on the effective diffusion coefficient ( D eff ).When the superficial fraction of the pore and the agglomerate are both equal to 50% in a monodisperse random configuration ( I xy = 1.0 + 0.0), D eff exhibits values within a range with a maximum of 6.20 × 10 -8 m 2 s −1 and a minimum close to the oxygen diffusion coefficient in the ionomer.As the superficial fraction of the pore decreases, the value of D eff decreases significantly.For instance, when the superficial fraction of the pore phase is 10% values close to the diffusion coefficient of oxygen in the ionomer are obtained.As the surface fraction of the agglomerate is predominant, the D eff is mostly influenced by the D O2,Agg , which impacts the overall performance of the fuel cell.
In conclusion, it has been observed that in samples with a surface fraction of the porous phase equal or higher than 50% at isotropy = 0, anisotropic material, the effective diffusion coefficient presents a higher value than when the isotropy is equal to 1.This decrease is attributed to the dispersion of the phases by the SA agitation process, the diffusion coefficient decreases exponentially until reaching a maximum of isotropy.At surface fractions greater than or equal to 50% of the pore phase the diffusion coefficient is influenced by the diffusion of oxygen in the pore, since it is the predominant phase.
On the other hand, as isotropic values approach unity, a change in the behavior of the effective diffusion is observed, with oxygen diffusion in the agglomerate being the predominant factor.This change becomes more evident in samples with porous phase surface fractions below 50%, where the effective diffusion coefficient shows an exponential decrease as the isotropy process progresses.
Ultimately, as the isotropic values approach unity, the effective diffusion coefficient falls into ranges close to the diffusion of oxygen in the ionomer, indicating that this phase becomes predominant in the microstructure.These findings are relevant to understanding and controlling diffusion processes in anisotropic materials, providing valuable information for the design and optimization of devices and systems that rely on the diffusion species.
The results obtained in this study provide information on the behavior of the Diffusion coefficients in different porosity conditions and how this directly affects the effective transport coefficient.In turn, it provides a methodology that is aimed at a better design and manufacture of catalytic layers of PEM fuel cells, leading to a better overall performance of the cell.

Figure 1 .
Figure 1.Block diagram of the implemented methodology.

( 2 )
j (r) to generate a vertical axis representation S (2) j,y (r) versus the horizontal deviation S

Figure 2 .
Figure 2. FVM boundary conditions and graphical index of the local fluxes used for the global J Eff .

Figure 3 .
Figure 3. Representative images of the SA process of the studied CL microstructures.

Figure 4 .
Figure 4. E SA error response (A), I xy isotropy response (B), and D eff coefficient (C), the three graphs as a function of I SA iteration.

29 Figure 5 .
Figure 5. D eff coefficient as a function of the isotropy evolution I xy : (A) for all possibilities and analytical solutions magnitudes @ I xy = 0 and (B) for a zoom section and numerical results @ I xy = 1.

Figure 6 .
Figure 6.D eff coefficient as a function of agglomerate surface fraction φ , for different values of isotropy I xy .

Table 1 .
Numerical values used in equation 15.the diffusion coefficient of oxygen in the agglomerate (data input: D O2,Agg = 8.45 × 10 −9 m 2 s −1 ) is noted by the dashed blue line.Zhao et al.