PDF
Species Transport in the Gas Diffusion Layers of a PEM
Introduction
This example focuses on the species transport within the gas diffusion layers (GDLs) of a proton exchange membrane (PEM) fuel cell. Specifically, the influence of anisotropic tortuosity on gas diffusivity in the GDLs is investigated. The geometry models a cell with two adjacent flow channels of different pressures, a situation that may occur in a cell with serpentine flow channels, or in a cell using an interdigitated flow field design. The model uses current balances, mass transport equations (Maxwell–Stefan diffusion for reactants, water and nitrogen gas), and momentum transport (Darcy’s law for the gas flows) to simulate a PEM fuel cell’s behavior.
Model Definition
The modeled section of the fuel cell consists of three domains: an anode (Ωa), a proton exchange membrane (Ωm), and a cathode (Ωc) as indicated in Figure 1.
Figure 1: Model geometry with domain and boundary labels.
Each of the electrodes (gas diffusion layers) is in contact with a gas distributing flow field, featuring channels (∂Ωa/c,inlet) and current collector ribs (∂Ωa/c,cc).
Humidified hydrogen and air are supplied to the inlet channels of the anode and cathode, respectively. Hydrogen reacts and is consumed at the anodic active layer to form protons that carry the ionic current to the cathode.
(1)
At the cathode, oxygen reacts together with the protons to form water at the active layer according to:
(2)
Both feed gases (humidified hydrogen and humidified air) are treated as ideal and are transported through diffusion and convection. The electrodes are treated as homogeneous porous media with uniform morphological properties such as porosity and permeability. The gas within each of the electrodes exists as a continuous phase so Darcy’s law applies.
An agglomerate model describes the electrode reactions in the active catalyst layers. The agglomerates consist of catalyst and carbon particles embedded in polymer electrolyte. The equations for the agglomerate model originate from the analytical solution of a diffusion-reaction problem in a spherical porous particle (Ref. 5 and Ref. 6). At the anodic active catalyst layer, hydrogen is the diffusing and reacting species in the agglomerates, while oxygen is the diffusing and reacting species in the agglomerates at the cathode. An agglomerate model of the cathodic active catalyst layer of a PEM fuel cell has been presented by Broka and others (Ref. 7 and Ref. 8).
Charge Balances
The Hydrogen Fuel Cell interface is used for modeling the potential distributions in the three domains, with the charge-transfer current density of the active layers described as interior boundary conditions, using Thin H2/O2 Gas Diffusion Electrode nodes. The local current density expressions for the anode and cathode, ia and ic, are specified according to the equation
where the index e stands for “a” (anode) or “c” (cathode), Ragg is the agglomerate radius (SI unit: m) and jagg,a and jagg,c (SI unit: A/m3) are the current densities given by the agglomerate model. Further, Lact is the active layer (thin gas diffusion electrode) thickness (m) and εmac its porosity (the macroscopic porosity). The specific surface area of the thin gas diffusion electrode (SI unit: 1/m) is specified according to the following equation as
Agglomerate Model for Anode and Cathode
The agglomerate model describes the current density in an active layer consisting of agglomerates of ionic conductor material and electrically conducting particles covered partially with catalyst. The local current density can be expressed analytically by solving a combination of the diffusion equation and the Butler–Volmer electrode kinetic equation for an agglomerate with constant electric and ionic potentials. The resulting equations for the current density in the anode and cathode are (Ref. 7)
where, again, the index e stands for “a” (anode) or “c” (cathode), and
In these equations, Dagg is the agglomerate gas diffusivity (SI unit: m2/s), ne is a “charge transfer” number (1 for the anode and –2 for the cathode), S is the specific area of the catalyst inside the agglomerate (SI unit: 1/m), and F is Faraday’s constant (SI unit: C/mol). Furthermore, ci,ref are the reference concentrations of the species (SI unit: mol/m3), ci,agg are the corresponding concentrations in the agglomerate surface (SI unit: mol/m3), i0a and i0c are the exchange current densities (SI unit: A/m2), R is the gas constant, T is the temperature (SI unit: K), and the overvoltages at the anode and the cathode are given by
where Eeq (SI unit: V) denotes the equilibrium voltage.
You set the anodic and cathodic reference states equal to the molar fractions at the inlet channels of the anode and cathode, respectively, at 1 atm. The dissolved hydrogen and oxygen concentrations at the surface of the agglomerates are related to the molar fractions of the respective species in the gas phase through Henry’s law
where K is Henry’s constant (SI unit: Pa·m3/mol).
Charge Balances, Continued
The potential difference between the cathode and anode current collectors corresponds to the total cell voltage. Choose the potential at the anode current collector as the reference level by setting it to zero. Then the total cell voltage serves as the boundary condition at the cathode current collector:
For the other boundaries you have electric insulation boundary conditions.
Porous Media Fluid Flow
To model the gas flows in the gas backings, this example uses the Darcy’s law in the Hydrogen Fuel Cell interface. The gas velocity is given by the continuity equation according to
where ρ is the mixture density of the gas phase (SI unit: kg/m3) and u denotes the gas velocity (SI unit: m/s). Darcy’s law for porous media states that the gradient of pressure, the viscosity of the fluid, and the structure of the porous media determine the velocity:
Here kp denotes the electrode’s permeability (SI unit: m2), η represents the gas viscosity (SI unit: Pa·s), and p is the pressure (SI unit: Pa). The ideal gas law gives the gas phase’s mixture density, ρ:
(3)
In this equation, R denotes the gas constant (SI unit: J/(mol·K)), T is the temperature (SI unit: K), M is the molar mass (SI unit: kg/mol), and x is the mole fraction.
At the inlets you specify the pressure as follows:
At the thin gas diffusion electrode boundary for the anode and cathode, the gas velocity is calculated automatically by the Hydrogen Fuel Cell interface, from the total mass flow given by the electrochemical reaction rate and the stoichiometric coefficients of Equation 1 and Equation 2, using Faraday’s law.
Maxwell-Stefan Mass Transport
The model takes into account two species in the anode — H2 and H2O — and three at the cathode — O2, H2O, and N2. The Hydrogen Fuel Cell interface uses Maxwell-Stefan multicomponent diffusion governed by the equations
to solve for the mass fractions, wi. This particular PEM fuel cell model assumes that the temperature-driven diffusion is insignificant and sets the source term, R, to zero. The Maxwell-Stefan binary diffusion coefficients, Dij (SI unit: m2/s), are calculated automatically by the interface. In the porous GDLs, effective binary diffusion coefficients, , need to be considered using a porous media diffusivity correction. This correction is based on the porosity, εg, and tortuosity, τg, giving the following binary diffusion coefficient definition
where the tortuosity is defined as a tensor in the form of a diagonal 2 x 2 matrix
The feed-gas mole fractions are specified at the inlets. At the thin gas diffusion electrode boundaries, the mass fluxes of the species are automatically determined from the electrochemical reaction rate and stoichiometric coefficients, using Faraday’s law. The Stefan velocity contributions are also automatically calculated by the interface.
The Membrane Transport features of the Hydrogen Fuel Cell interface are used to model the transport of water in the ionomer phase in the membrane domain. The molecular flux of water depends both on chemical potential gradient driven permeation and electroosmotic drag, using experimentally estimated parameters available in the Nafion material in the Fuel Cell and Electrolyzer Material Library.
Results and Discussion
Figure 2 shows polarization curves for the PEM fuel cell. Four cases are displayed to visualize the influence of tortuosity in the GDLs on the performance. In the anisotropic cases, the tortuosity settings correspond to more diffusivity limitations either along the height (y) or the thickness (x) of each GDL.
Results indicate that the cell performance is particularly sensitive to an increased tortuosity in the x direction and is due to the fact that reactants need to pass through the layer thickness to reach the thin electrode layer. However, the increased tortuosity in the y direction represents the most realistic scenario since the in-plane GDL structure, parallel to the y direction, is usually more tortuous.
Figure 2: Polarization curves for the PEM fuel cell for different tortuosity conditions in the GDLs.
For the more realistic anisotropic tortuosity case (τxx = 1.1and τyy = 2.8), a small but evident difference increasing with current density is shown compared to the corresponding isotropic tortuosity case (τ = 1.1).
The current density at the active layer at the anode is displayed as a function of cell height in Figure 3. An uneven current distribution can be seen for both plotted cases, with the situation being more pronounced with anisotropic tortuosity.
Figure 3: Current-density distribution at the anodic active layer for the PEM fuel cell operating at 0.4 V for two GDL tortuosity cases.
The current density is the highest in the cell’s upper region but with even less activity in the center of the layer if anisotropic tortuosity is considered. For both cases, the oxygen-reduction reaction rate in the cathode determines the current-density distribution. The maximum current density arises close to the air inlet.
Figure 4 and Figure 5 depict the reactant (oxygen and hydrogen) weight fractions in the cathode and anode gases for the isotropic and anisotropic cases, respectively. The reactant minimum is decreasing substantially in the latter case and is consistent with the minimum in the current-density distribution.
Figure 4: Reactant mass fractions on the anode side (left) and cathode side (right) at 0.4 V with isotropic GDL tortuosity. The reactant in the anode is hydrogen and that in the cathode is oxygen.
Figure 5: Reactant mass fractions on the anode side (left) and cathode side (right) at 0.4 V with anisotropic GDL tortuosity. The reactant in the anode is hydrogen and that in the cathode is oxygen.
Figure 6 and Figure 7 show the current distribution over the fuel cell for the isotropic and anisotropic cases, respectively. The minimum near the center of the cell is slightly lower in the latter case, but the difference is small and local. This displays that over the cell level, the system is not especially affected by higher in-plane tortuosity for species diffusion in the GDLs. A behavior also manifested in the small loss in performance (Figure 2).
Figure 6: Current density (surface plot) and current vector field (arrow plot) in fuel cell operating at 0.4 V with isotropic GDL tortuosity. The anode is on the left and the cathode is on the right.
.
Figure 7: Current density (surface plot) and current vector field (arrow plot) in fuel cell operating at 0.4 V with anisotropic GDL tortuosity. The anode is on the left and the cathode is on the right.
Other processes accounted for in the model, mainly convective mass-transport flux, dictate the fuel cell performance. In Figure 8, the velocity field for the PEM is shown. The field is close to identical for all investigated tortuosities. It should be noted that the gas permeability is set be tortuosity independent in this model.
Figure 8: Gas velocity field in the anode and cathode compartments at 0.4 V.
References
1. W. He, J.S. Yi, and T.V. Nguyen, “Two-Phase Flow Model of the Cathode of PEM Fuel Cells Using Interdigitated Flow Fields,” AIChe J., vol. 46, pp. 2053–2063, 2000.
2. C. Marr and X. Li, “Composition and Performance Modelling of Catalyst Layer in a Proton Exchange Membrane Fuel Cell,” J. Power Sources, vol. 77, pp. 17–27, 1999.
3. P. Futerko and I.-M. Hsing, “Two-Dimensional Finite Element Method Study of the Resistance of Membranes in Polymer Electrolyte Fuel Cells,” Electrochimica Acta, vol. 45, pp. 1741–1751, 2000.
4. D.M. Bernardi and M.W. Verbrugge, “Mathematical Model of a Gas Diffusion Electrode Bonded to a Polymer Electrolyte,” AIChe J., vol. 37, pp. 1151–1163, 1991.
5. H. Scott Fogler, Elements of Chemical Reaction Engineering, 3rd ed., Prentice Hall, 1999.
6. R.B. Bird, W.E. Stewart, and E.N. Lightfoot, Transport Phenomena, John Wiley & Sons, 1960.
7. K. Broka and P. Ekdunge, “Modelling the PEM fuel cell cathode,” J. Appl. Electrochem., vol. 27, pp. 281–289, 1997.
8. K. Dannenberg, P. Ekdunge, and G. Lindbergh, “Mathematical model of the PEMFC,” J. Appl. Electrochem., vol. 30, pp. 1377–1387, 2000.
Application Library path: Fuel_Cell_and_Electrolyzer_Module/Fuel_Cells/pem_gdl_species_transport_2d
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Electrochemistry > Hydrogen Fuel Cells > Proton Exchange Membrane (fc).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select Preset Studies for Selected Physics Interfaces > Stationary with Initialization.
6
Click  Done.
Geometry 1
Create the geometry using rectangles.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 2.5e-4.
4
In the Height text field, type 2e-3.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 1e-4.
4
In the Height text field, type 2e-3.
5
Locate the Position section. In the x text field, type 2.5e-4.
Rectangle 3 (r3)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 2.5e-4.
4
In the Height text field, type 2e-3.
5
Locate the Position section. In the x text field, type 3.5e-4.
Rectangle 4 (r4)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 1e-4.
4
In the Height text field, type 5e-4.
5
Locate the Position section. In the x text field, type -1e-4.
Array 1 (arr1)
1
In the Geometry toolbar, click  Transforms and choose Array.
2
Select the object r4 only.
3
In the Settings window for Array, locate the Size section.
4
In the x size text field, type 2.
5
In the y size text field, type 2.
6
Locate the Displacement section. In the x text field, type 7e-4.
7
In the y text field, type 1.5e-3.
Union 1 (uni1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Union.
2
Select the objects arr1(1,1), arr1(1,2), and r1 only.
3
In the Settings window for Union, locate the Union section.
4
Clear the Keep interior boundaries checkbox.
Union 2 (uni2)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Union.
2
Select the objects arr1(2,1), arr1(2,2), and r3 only.
3
In the Settings window for Union, locate the Union section.
4
Clear the Keep interior boundaries checkbox.
Fillet 1 (fil1)
1
In the Geometry toolbar, click  Fillet.
2
In the Settings window for Fillet, locate the Points section.
3
Click to clear the  Activate Selection toggle button for Vertices to fillet.
4
Click the  Zoom Extents button in the Graphics toolbar.
5
On the object uni1, select Points 6 and 7 only.
6
On the object uni2, select Points 4 and 5 only.
7
Click to select the  Activate Selection toggle button for Vertices to fillet.
8
On the object uni1, select Points 6 and 7 only.
9
On the object uni2, select Points 4 and 5 only.
10
Locate the Radius section. In the Radius text field, type 1e-5.
11
In the Geometry toolbar, click  Build All.
Global Definitions
Proceed to load a set of global model parameters from a text file provided with the Application Library.
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file pem_gdl_species_transport_2d_parameters.txt.
Definitions
Set up a number of selections on the geometry. These will be used later when setting up the physics.
Anode GDL
1
In the Definitions toolbar, click  Explicit.
2
Select Domain 1 only.
3
In the Settings window for Explicit, type Anode GDL in the Label text field.
Membrane
1
In the Definitions toolbar, click  Explicit.
2
Select Domain 2 only.
3
In the Settings window for Explicit, type Membrane in the Label text field.
Cathode GDL
1
In the Definitions toolbar, click  Explicit.
2
Select Domain 3 only.
3
In the Settings window for Explicit, type Cathode GDL in the Label text field.
Anode GDE
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, locate the Input Entities section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 10 only.
5
In the Label text field, type Anode GDE.
Cathode GDE
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, locate the Input Entities section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 13 only.
5
In the Label text field, type Cathode GDE.
Anode GDL Variables
Define a number of domain- and boundary-specific variables. These will be used both for setting up the physics and in postprocessing.
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, type Anode GDL Variables in the Label text field.
3
Locate the Geometric Entity Selection section. From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Anode GDL.
5
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
wReact
fc.wH2
 
Mass fraction, reactant
Cathode GDL Variables
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, type Cathode GDL Variables in the Label text field.
3
Locate the Geometric Entity Selection section. From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Cathode GDL.
5
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
wReact
fc.wO2
 
Mass fraction, reactant
Anode GDE Variables
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, type Anode GDE Variables in the Label text field.
3
Locate the Geometric Entity Selection section. From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose Anode GDE.
5
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
cH2_agg
fc.p*fc.xH2/KH2
 
Henry's law hydrogen agglomerate concentration
eta_a
fc.phis-fc.phil-E_eq_a
V
Anodic overpotential
beta_a
cH2_agg-cH2_ref*exp(-2*F_const*eta_a/(R_const*T))
 
 
lda_a
sqrt(i0_a*S*R_agg^2/(2*F_const*cH2_ref*D_agg))
 
Anodic current density subexpression
i_a
K*(1-lda_a*coth(lda_a))*beta_a*(R_agg/3)
 
Anode current density
Cathode GDE Variables
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, type Cathode GDE Variables in the Label text field.
3
Locate the Geometric Entity Selection section. From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose Cathode GDE.
5
Locate the Variables section. In the table, enter the following settings:
Name
Expression
Unit
Description
cO2_agg
fc.p*fc.xO2/KO2
 
Henry's law oxygen agglomerate concentration
eta_c
fc.phis-fc.phil-E_eq_c
V
Cathodic overvoltage
lda_c
sqrt(i0_c*S*R_agg^2*exp(-F_const*eta_c/(2*R_const*T))/(4*F_const*cO2_ref*D_agg))
 
Cathodic current density subexpression
i_c
-2*K*(1-lda_c*coth(lda_c))*cO2_agg*(R_agg/3)
 
Cathode current density
Average collector
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Average.
2
In the Settings window for Average, type Average collector in the Label text field.
3
In the Operator name text field, type aveop_ca.
4
Locate the Source Selection section. From the Geometric entity level list, choose Boundary.
5
Select Boundaries 8, 23, and 24 only.
Materials
This model uses polymer electrolyte material (Nafion 1100, vapor equilibrated) that is available in the material library. Add the material twice, and assign it first to the membrane domain, and next to the membrane boundaries adjacent to the gas phase nodes.
Add Material
1
In the Materials toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Fuel Cell and Electrolyzer > Polymer Electrolytes > Nafion®, EW 1100, Vapor Equilibrated, Protonated.
4
Right-click and choose Add to Component 1 (comp1).
Materials
Nafion®, EW 1100, Vapor Equilibrated, Protonated (mat1)
1
In the Model Builder window, under Component 1 (comp1) > Materials click Nafion®, EW 1100, Vapor Equilibrated, Protonated (mat1).
2
In the Settings window for Material, locate the Geometric Entity Selection section.
3
Click  Clear Selection.
4
Select Domain 2 only.
Add Material
1
Go to the Add Material window.
2
In the tree, select Fuel Cell and Electrolyzer > Polymer Electrolytes > Nafion®, EW 1100, Vapor Equilibrated, Protonated.
3
Right-click and choose Add to Component 1 (comp1).
4
In the Materials toolbar, click  Add Material to close the Add Material window.
Materials
Nafion®, EW 1100, Vapor Equilibrated, Protonated 1 (mat2)
1
In the Settings window for Material, locate the Geometric Entity Selection section.
2
From the Geometric entity level list, choose Boundary.
3
Select Boundaries 10 and 13 only.
Hydrogen Fuel Cell (fc)
Set up the current distribution and transport model. Include mass transport using Maxwell-Stefan diffusion and momentum transport using Darcy’s Law in both the anode and cathode gas mixtures, along with electroosmotic water drag in the membrane. Also, define the reference pressure level in the interface properties. Note that the default gas species are hydrogen and water on the anode side, and oxygen, nitrogen and water on the cathode side. Start with adding the relevant domain nodes.
1
In the Model Builder window, under Component 1 (comp1) click Hydrogen Fuel Cell (fc).
2
In the Settings window for Hydrogen Fuel Cell, locate the H2 Gas Mixture section.
3
Find the Transport mechanisms subsection. Select the Use Darcy’s Law for momentum transport checkbox.
4
Locate the O2 Gas Mixture section. Select the Use Darcy’s Law for momentum transport checkbox.
5
Click to expand the Electrolyte and Membrane Transport section. Select the Electroosmotic water drag checkbox.
6
Click to expand the Reference Pressure Level section. In the pref text field, type 0.
Membrane 1
1
In the Physics toolbar, click  Domains and choose Membrane.
2
In the Settings window for Membrane, locate the Domain Selection section.
3
From the Selection list, choose Membrane.
H2 Gas Diffusion Layer 1
1
In the Physics toolbar, click  Domains and choose H2 Gas Diffusion Layer.
2
In the Settings window for H2 Gas Diffusion Layer, locate the Domain Selection section.
3
From the Selection list, choose Anode GDL.
O2 Gas Diffusion Layer 1
1
In the Physics toolbar, click  Domains and choose O2 Gas Diffusion Layer.
2
In the Settings window for O2 Gas Diffusion Layer, locate the Domain Selection section.
3
From the Selection list, choose Cathode GDL.
Thin H2 Gas Diffusion Electrode 1
1
In the Physics toolbar, click  Boundaries and choose Thin H2 Gas Diffusion Electrode.
2
In the Settings window for Thin H2 Gas Diffusion Electrode, locate the Boundary Selection section.
3
From the Selection list, choose Anode GDE.
Thin O2 Gas Diffusion Electrode 1
1
In the Physics toolbar, click  Boundaries and choose Thin O2 Gas Diffusion Electrode.
2
In the Settings window for Thin O2 Gas Diffusion Electrode, locate the Boundary Selection section.
3
From the Selection list, choose Cathode GDE.
In the Electrolyte Phase node, the electrolyte conductivity is set to be taken from the Materials node. Inspect the settings in the H2 Gas Phase and O2 Gas Phase nodes. Note that the density and viscosity of the gas mixture, and the binary diffusion coefficients are calculated automatically when the respective default settings are used.
The properties for electroosmotic water drag in the Membrane node and in the child nodes that added by default are automatically taken from the Materials node.
Initial Values 1
1
In the Model Builder window, expand the Membrane 1 node, then click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T0 text field, type T.
H2 Gas Diffusion Layer 1
Set up the properties of the H2 Gas Diffusion Layer and O2 Gas Diffusion Layer nodes.
1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) click H2 Gas Diffusion Layer 1.
2
In the Settings window for H2 Gas Diffusion Layer, locate the Electrode Charge Transport section.
3
In the σs text field, type kappa_s.
4
Locate the Gas Transport section. From the Effective diffusivity correction list, choose Tortuosity.
5
In the εg text field, type eps_mac.
6
From the list, choose Diagonal.
7
Specify the τg matrix as
tortx
0
0
torty
8
In the κg text field, type kappa_p.
O2 Gas Diffusion Layer 1
1
In the Model Builder window, click O2 Gas Diffusion Layer 1.
2
In the Settings window for O2 Gas Diffusion Layer, locate the Electrode Charge Transport section.
3
In the σs text field, type kappa_s.
4
Locate the Gas Transport section. From the Effective diffusivity correction list, choose Tortuosity.
5
In the εg text field, type eps_mac.
6
From the list, choose Diagonal.
7
Specify the τg matrix as
tortx
0
0
torty
8
In the κg text field, type kappa_p.
Thin H2 Gas Diffusion Electrode 1
Set up the properties of the Thin H2 Gas Diffusion Electrode node. The details of electrode kinetics are set in the child node.
1
In the Model Builder window, click Thin H2 Gas Diffusion Electrode 1.
2
In the Settings window for Thin H2 Gas Diffusion Electrode, locate the Electrode Thickness section.
3
In the dgde text field, type l_act.
Thin H2 Gas Diffusion Electrode Reaction 1
1
In the Model Builder window, click Thin H2 Gas Diffusion Electrode Reaction 1.
2
In the Settings window for Thin H2 Gas Diffusion Electrode Reaction, locate the Equilibrium Potential section.
3
From the Eeq list, choose User defined. Due to the user-defined kinetics expression, the equilibrium potential parameter is not in use in this model. Therefore keep the default value of 0 V.
4
Locate the Electrode Kinetics section. From the iloc,expr list, choose User defined. In the associated text field, type i_a.
5
Locate the Active Specific Surface Area section. In the av text field, type Av.
Thin O2 Gas Diffusion Electrode 1
Similarly, set up the properties of the Thin O2 Gas Diffusion Electrode node. The details of electrode kinetics are set in the child node.
1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) click Thin O2 Gas Diffusion Electrode 1.
2
In the Settings window for Thin O2 Gas Diffusion Electrode, locate the Electrode Thickness section.
3
In the dgde text field, type l_act.
Thin O2 Gas Diffusion Electrode Reaction 1
1
In the Model Builder window, click Thin O2 Gas Diffusion Electrode Reaction 1.
2
In the Settings window for Thin O2 Gas Diffusion Electrode Reaction, locate the Equilibrium Potential section.
3
From the Eeq list, choose User defined. Also for the oxygen reduction reaction, the equilibrium potential parameter is not in use in this model. Therefore keep the default value of 0 V.
4
Locate the Electrode Kinetics section. From the iloc,expr list, choose User defined. In the associated text field, type i_c.
5
Locate the Active Specific Surface Area section. In the av text field, type Av.
Finally, set up the exterior boundary conditions and the initial values.
Electronic Conducting Phase 1
Set the initial value for electric potential in the cathode electrode to the cell potential. (The default zero initial values are used for both potentials in the rest of the geometry.)
1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) click Electronic Conducting Phase 1.
Initial Values, O2 Domains 1
1
In the Physics toolbar, click  Attributes and choose Initial Values, O2 Domains.
2
Select Domain 3 only.
3
In the Settings window for Initial Values, O2 Domains, locate the Initial Values section.
4
In the ϕs text field, type V_cell.
Electronic Conducting Phase 1
In the Model Builder window, click Electronic Conducting Phase 1.
Electric Ground 1
1
In the Physics toolbar, click  Attributes and choose Electric Ground.
2
Select Boundaries 8, 23, and 24 only.
Electronic Conducting Phase 1
In the Model Builder window, click Electronic Conducting Phase 1.
Electric Potential 1
1
In the Physics toolbar, click  Attributes and choose Electric Potential.
2
Select Boundary 17 only.
3
In the Settings window for Electric Potential, locate the Electric Potential section.
4
In the ϕs,bnd text field, type V_cell.
H2 Gas Phase 1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) click H2 Gas Phase 1.
H2 Inlet 1
1
In the Physics toolbar, click  Attributes and choose H2 Inlet.
2
Select Boundary 1 only.
3
In the Settings window for H2 Inlet, locate the Inlet Flow Type section.
4
From the Inlet flow type list, choose Mixture composition constraint.
5
Locate the Flow Boundary Condition section. From the list, choose Pressure.
6
In the p0 text field, type p_a_in.
H2 Gas Phase 1
In the Model Builder window, click H2 Gas Phase 1.
H2 Inlet 2
1
In the Physics toolbar, click  Attributes and choose H2 Inlet.
2
Select Boundary 4 only.
3
In the Settings window for H2 Inlet, locate the Inlet Flow Type section.
4
From the Inlet flow type list, choose Mixture composition constraint.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Composition section.
3
In the x0,H2O text field, type xH2Oa_in.
4
Locate the Initial Pressure section. In the p0 text field, type p_ref.
O2 Gas Phase 1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) click O2 Gas Phase 1.
O2 Inlet 1
1
In the Physics toolbar, click  Attributes and choose O2 Inlet.
2
Select Boundary 22 only.
3
In the Settings window for O2 Inlet, locate the Inlet Flow Type section.
4
From the Inlet flow type list, choose Mixture composition constraint.
5
Locate the Flow Boundary Condition section. From the list, choose Pressure.
6
In the p0 text field, type p_c_in.
O2 Gas Phase 1
In the Model Builder window, click O2 Gas Phase 1.
O2 Inlet 2
1
In the Physics toolbar, click  Attributes and choose O2 Inlet.
2
Select Boundary 21 only.
3
In the Settings window for O2 Inlet, locate the Inlet Flow Type section.
4
From the Inlet flow type list, choose Mixture composition constraint.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Composition section.
3
In the x0,H2O text field, type xH2Oc_in.
4
In the x0,N2 text field, type xN2_in.
5
Locate the Initial Pressure section. In the p0 text field, type p_ref.
Global Definitions
Default Model Inputs
Since we will be using the same temperature everywhere in the model, we will define the temperature only once in the Default Model Inputs node. This node may be accessed by multiple physics nodes.
1
In the Model Builder window, under Global Definitions click Default Model Inputs.
2
In the Settings window for Default Model Inputs, locate the Browse Model Inputs section.
3
In the tree, select General > Temperature (K) - minput.T.
4
Find the Expression for remaining selection subsection. In the Temperature text field, type T.
Mesh 1
Steep gradients are expected close to the electrodes and at the corner points where the ribs of the flow plates are suppressed into the GDL. Create a mesh with finer resolution at these parts of the geometry.
Size 1
1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundaries 10 and 13 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type 2.5e-5.
Size 2
1
In the Model Builder window, right-click Mesh 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Point.
4
Select Points 5, 6, 16, and 17 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type 1e-5.
Free Triangular 1
1
In the Mesh toolbar, click  Free Triangular.
2
Right-click Mesh 1 and choose Build All.
The finalized mesh should now look as follows:
Study 1
The problem is now ready for solving. In the first step, solve for secondary current distribution initialization.
Step 1: Current Distribution Initialization
1
In the Model Builder window, under Study 1 click Step 1: Current Distribution Initialization.
2
In the Settings window for Current Distribution Initialization, locate the Study Settings section.
3
From the Current distribution type list, choose Secondary.
Step 2: Stationary
1
In the Model Builder window, click Step 2: Stationary.
2
In the Settings window for Stationary, click to expand the Study Extensions section.
3
Select the Auxiliary sweep checkbox.
4
Click  Add.
5
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
V_cell (Cell voltage)
range(1,-0.1,0.4)
V
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
tortx (Tortuosity, x direction)
1.1 2.8
 
5
Click  Add.
6
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
torty (Tortuosity, y direction)
1.1 2.8
 
7
From the Sweep type list, choose All combinations.
8
In the Study toolbar, click  Compute.
Results
The following steps reproduce the figures found in the Results and Discussion section.
Polarization curves
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Polarization curves in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
4
Click to expand the Title section. From the Title type list, choose None.
5
Locate the Plot Settings section.
6
Select the x-axis label checkbox. In the associated text field, type Average cell current density (A/cm<sup>2</sup>).
7
Select the y-axis label checkbox. In the associated text field, type Cell voltage (V).
8
Select the Flip the x- and y-axes checkbox.
Global 1
1
Right-click Polarization curves and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
aveop_ca(-fc.nIs)
A/cm^2
Average current density (A/cm>sup>2</sup>)
4
Click to expand the Legends section. From the Legends list, choose Manual.
5
In the table, enter the following settings:
Legends
Isotropic, \tau<sub>xx</sub>=1.1 & \tau<sub>yy</sub>=1.1
Anisotropic, \tau<sub>xx</sub>=1.1 & \tau<sub>yy</sub>=2.8
Anisotropic, \tau<sub>xx</sub>=2.8 & \tau<sub>yy</sub>=1.1
Isotropic, \tau<sub>xx</sub>=2.8 & \tau<sub>yy</sub>=2.8
6
In the Polarization curves toolbar, click  Plot.
7
Click the  Zoom Extents button in the Graphics toolbar.
GDL Current Density Distribution
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type GDL Current Density Distribution in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
4
From the Parameter value (tortx) list, choose 1.1.
5
From the Parameter value (torty) list, choose 1.1.
Surface 1
1
Right-click GDL Current Density Distribution and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Hydrogen Fuel Cell > fc.IsMag - Electrode current density magnitude - A/m².
3
Click to expand the Range section. Select the Manual color range checkbox.
4
In the Maximum text field, type 1e4.
5
Click to expand the Coloring and Style section. From the Color table list, choose Prism.
Arrow Surface 1
1
In the Model Builder window, right-click GDL Current Density Distribution and choose Arrow Surface.
2
In the Settings window for Arrow Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Hydrogen Fuel Cell > fc.Isx,fc.Isy - Electrode current density vector.
3
Click to expand the Coloring and Style section.
4
Select the Scale factor checkbox. In the associated text field, type 5e-9.
5
In the GDL Current Density Distribution toolbar, click  Plot.
6
Click the  Zoom Extents button in the Graphics toolbar.
GDL Current Density Distribution
1
In the Model Builder window, click GDL Current Density Distribution.
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Parameter value (torty) list, choose 2.8.
4
In the GDL Current Density Distribution toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
Anode Reaction Current Density
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Anode Reaction Current Density in the Label text field.
3
Locate the Data section. From the Dataset list, choose None.
4
Click to expand the Title section. From the Title type list, choose None.
5
Locate the Legend section. From the Position list, choose Lower right.
Line Graph 1
1
Right-click Anode Reaction Current Density and choose Line Graph.
2
In the Settings window for Line Graph, locate the Data section.
3
From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
4
From the Parameter selection (tortx) list, choose First.
5
From the Parameter selection (torty) list, choose First.
6
From the Parameter selection (V_cell) list, choose Last.
7
Locate the Selection section. From the Selection list, choose Anode GDE.
8
Click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Definitions > Variables > i_a - Anode current density - A/m².
9
Locate the y-Axis Data section. In the Expression text field, type i_a*Av*l_act.
10
Select the Description checkbox. In the associated text field, type Anode reaction current density.
11
Click Replace Expression in the upper-right corner of the x-Axis Data section. From the menu, choose Component 1 (comp1) > Geometry > Coordinate > y - y-coordinate.
12
Click to expand the Legends section. Select the Show legends checkbox.
13
From the Legends list, choose Manual.
14
In the table, enter the following settings:
Legends
Isotropic, \tau<sub>xx</sub>=1.1 & \tau<sub>yy</sub>=1.1
Line Graph 2
1
Right-click Line Graph 1 and choose Duplicate.
2
In the Settings window for Line Graph, locate the Data section.
3
From the Parameter selection (torty) list, choose Last.
4
Locate the Legends section. In the table, enter the following settings:
Legends
Anisotropic, \tau<sub>xx</sub>=1.1 & \tau<sub>yy</sub>=2.8
5
In the Anode Reaction Current Density toolbar, click  Plot.
6
Click the  Zoom Extents button in the Graphics toolbar.
Velocity Field
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Velocity Field in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
4
From the Parameter value (tortx) list, choose 1.1.
Surface 1
1
Right-click Velocity Field and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Hydrogen Fuel Cell > fc.U - Velocity magnitude - m/s.
3
Locate the Coloring and Style section. From the Color table list, choose Prism.
Arrow Surface 1
1
In the Model Builder window, right-click Velocity Field and choose Arrow Surface.
2
In the Settings window for Arrow Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Hydrogen Fuel Cell > fc.u,fc.v - Velocity field.
3
Locate the Coloring and Style section.
4
Select the Scale factor checkbox. In the associated text field, type 0.0035.
5
In the Velocity Field toolbar, click  Plot.
6
Click the  Zoom Extents button in the Graphics toolbar.
Reactant Mass Fraction
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Reactant Mass Fraction in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol3).
4
From the Parameter value (tortx) list, choose 1.1.
5
From the Parameter value (torty) list, choose 1.1.
Surface 1
1
Right-click Reactant Mass Fraction and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type wReact.
4
Locate the Coloring and Style section. From the Color table list, choose Prism.
5
In the Reactant Mass Fraction toolbar, click  Plot.
6
Click the  Zoom Extents button in the Graphics toolbar.
Reactant Mass Fraction
1
In the Settings window for 2D Plot Group, locate the Data section.
2
From the Parameter value (torty) list, choose 2.8.
3
In the Model Builder window, expand the Reactant Mass Fraction node.
Electrode Potential with Respect to Ground (fc), Electrolyte Potential (fc)
1
In the Model Builder window, under Results, Ctrl-click to select Electrode Potential with Respect to Ground (fc) and Electrolyte Potential (fc).
2
Right-click and choose Delete.