PDF
Fuel Cell Stack Cooling
Introduction
This tutorial models the thermal management of a polymer electrolyte membrane (PEM) fuel cell stack. Operating the stack with similar temperature profile for all cells is important since an uneven temperature distribution may otherwise result in nonuniform water vapor condensation, and a large cell-to-cell variation in performance.
The stack consists of five cells, interlayered with bipolar plates that carry the liquid cooling fluid. The model solves for the temperature, the electrode and electrolyte phase potentials, the mass transport of the reacting species in each separate gas compartment, and the fluid pressures and corresponding velocities in the gas and liquid flow compartments.
The flow channels of the bipolar plates are not resolved explicitly in the geometry. Instead anisotropic permeabilities are used define the fluid flow patterns inside the stack.
Model Definition
Figure 1 depicts the repetitive unit cell which is used to construct the model geometry. The bipolar plate (BPP), constructed by patterned steel sheets, separates the hydrogen and oxygen (air) gas compartments, and also contains a separate flow compartment for the cooling liquid. The BPP also contains inlet and outlet manifolds to conduct the fluids to larger additional flow channels located outside the model geometry. The stack in this tutorial uses a “Z” pattern for the gas flow channels, with the direction of the air channels depicted in Figure 1. For the hydrogen side, the same pattern is used, mirrored in the yz-plane. For the cooling flow, a uniform flow direction along the y-axis is assumed. Between each BPP, membrane-electrode assemblies (MEAs) are placed, consisting of one gas diffusion layer (GDL) and one gas diffusion electrode (GDE) per gas compartment, and a polymer electrolyte membrane in between.
It should be noted that the cooling flow compartments are formed by the union of the hydrogen and the oxygen sides of the BPPs, when stacking the unit cells on top of each other. Generally in stack design, as a result of only having to provide flow of one gas, the first and end plates have to be designed differently compared to the inner plates. In this tutorial the first and end plates are constructed by adding additional “half” BPPs, using the same design as in Figure 1 but with the inactive gas compartment and manifold removed. As a result of this, the cooing flow to the first and last cells in the stack will become larger than for the inner cells.
.
Figure 1: Repetitive unit cell. The larger arrows indicate the orientation of the flow channels in the oxygen (air) compartment of the bipolar plate. The figure is scaled ten times in the z direction.
The individual channels in the patterned BPPs are not resolved in the model geometry. Instead the BPPs are divided into separate domains corresponding to the main gas flow direction. Homogenized flow equations based on Darcy’s Law are then used to define the flow, using anisotropic permeabilities corresponding to the direction of the channels.
The final geometry, as shown in Figure 2, is constructed by stacking five unit cells on top of each other, with the added front and end cooling “half” BPPs, as discussed above. Additional metal end blocks, used for conducting the current and for providing structural rigidity, are placed outside the stack.
The geometry is finalized as an assembly, with the assembly boundary located in the middle of the membranes. This allows since for nonmatching meshes on each side of the mid-membrane boundaries, and hence for sweeping the mesh in the through-plane direction.
.
Figure 2: Model geometry consisting of five unit cells and two end blocks.
The model is defined using one Hydrogen Fuel Cell, one Heat Transfer and one Darcy’s Law interface.
The Hydrogen Fuel Cell defines the electrochemical reactions, the electrode and electrolyte charge transport, gas phase diffusion and convective flow, as well as membrane water transport. For an introduction to the hydrogen fuel cell interface, see for instance the Mass Transport and Electrochemical Reaction in a Fuel Cell Cathode and the Low-Temperature PEM Fuel Cell with Serpentine Flow Field examples.
The Heat Transfer interface solves for the temperature in the stack, applying the heat sources computed by the Hydrogen Fuel Cell interface through the Electrochemical Heating multiphysics node. See also the Nonisothermal PEM Fuel Cell example. Heat is only assumed to exit the stack through the cooling manifold outlets.
The additional Darcy’s Law interface is used to compute the flow profile of the cooling liquid, which is used in the Heat Transfer interface.
The model is solved for a range of average cell voltages using a stationary solver.
Results and Discussion
Figure 3 shows the electrode phase potential in the stack for an average cell voltage of 0.55 V.
Figure 3: Electrode phase potential.
Figure 4 shows the hydrogen molar fraction. As a result of the consumption of hydrogen, the molar fraction decreases toward the outlet.
Figure 4: Hydrogen molar fraction.
Similarly, Figure 5 shows the oxygen molar fraction in the air gas compartment. Depletion effects are somewhat more severe than for the hydrogen side.
Figure 5: Oxygen molar fraction.
Figure 6 shows the cooling flow streamlines and the corresponding temperature. The temperature increases toward the outlet.
Figure 6: Temperature (color legend) and cooling flow stream lines.
Figure 7 shows the temperature in the MEAs. The temperature is higher toward the membrane.
Figure 7: Temperature in the MEAs.
Figure 8 elucidates the temperature profile in the through-plane direction of the stack further. Here the temperature is plotted along a line in the z direction, placed in the middle of the stack in the x direction, and approximately 1 cm from the outlet cooling manifold in the y direction.
Figure 8: Temperature along a z-directed line through the stack, close to the outlet at half of the width of the cells, for different average cell voltages.
As can be seen, polarizing the stack results in significant temperature gradients in the MEA. We also see lower temperatures in the first and the last cells, as an effect of these cells receiving, on average, more cooling flow compared to the inner cells.
Finally, Figure 9 shows the water activity (relative humidity) in the oxygen gas diffusion layer.
Figure 9: Water activity in the oxygen gas diffusion electrodes.
The relative humidity is related both to water production and the temperature. As a result of the water production in the cell, one would generally expect the relative humidity to be higher toward the outlet. However, at the same time the higher temperature toward the outlet counteracts this, resulting in a more complex behavior, depending both on the straight cooling fluid and the gas “Z”-flow field design.
Also note that the humidity levels are generally higher in the first and last cell of the stack. This is due to the generally lower temperatures in these two cells. This indicates that these two cells will be more susceptible to water flooding. To improve the stack design one could consider modifying the cooling flow patters to reduce the cooling of the first and the end cells.
Application Library path: Fuel_Cell_and_Electrolyzer_Module/Fuel_Cells/stack_cooling
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  3D.
2
In the Select Physics tree, select Electrochemistry > Hydrogen Fuel Cells > Proton Exchange Membrane (fc).
3
Click Add.
4
In the Select Physics tree, select Fluid Flow > Porous Media and Subsurface Flow > Darcy’s Law (dl).
5
Click Add.
6
In the Select Physics tree, select Heat Transfer > Porous Media > Heat Transfer in Porous Media (ht).
7
Click Add.
8
Click  Study.
9
In the Select Study tree, select Preset Studies for Selected Physics Interfaces > Hydrogen Fuel Cell > Stationary with Initialization.
10
Click  Done.
Geometry 1
Insert the geometry sequence from a file.
1
In the Geometry toolbar, click Insert Sequence and choose Insert Sequence.
2
Browse to the model’s Application Libraries folder and double-click the file stack_cooling_geom_sequence.mph.
3
In the Geometry toolbar, click  Build All.
4
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
Disable the analysis of the geometry as the remaining small geometric details are needed.
5
In the Model Builder window, click Geometry 1.
6
In the Settings window for Geometry, locate the Cleanup section.
7
Clear the Automatic detection of small details checkbox.
8
In the Geometry toolbar, click  Build All.
9
In the Model Builder window, collapse the Geometry 1 node.
Global Definitions
Geometry Parameters
Some parameters were imported when you imported the geometry sequence. Load some additional parameters required for setting up the physics from a text file.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, type Geometry Parameters in the Label text field.
Physics Parameters
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Physics Parameters in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file stack_cooling_physics_parameters.txt.
Add Material
1
In the Materials toolbar, click  Add Material to open the Add Material window.
Now add some materials from the Material Libraries.
2
Go to the Add Material window.
3
In the tree, select Built-in > Steel AISI 4340.
4
Right-click and choose Add to Component 1 (comp1).
5
In the tree, select Built-in > Water, liquid.
6
Right-click and choose Add to Component 1 (comp1).
7
In the tree, select Fuel Cell and Electrolyzer > Polymer Electrolytes > Nafion®, EW 1100, Vapor Equilibrated, Protonated.
8
Right-click and choose Add to Component 1 (comp1).
9
In the Materials toolbar, click  Add Material to close the Add Material window.
Materials
Steel AISI 4340 (mat1)
Assign the materials you added to different parts of the geometry.
1
In the Settings window for Material, locate the Geometric Entity Selection section.
2
From the Selection list, choose Current Collector and Feeder Plates.
Water, liquid (mat2)
1
In the Model Builder window, click Water, liquid (mat2).
2
In the Settings window for Material, locate the Geometric Entity Selection section.
3
From the Selection list, choose Cooling Channels and Manifolds.
Nafion®, EW 1100, Vapor Equilibrated, Protonated (mat3)
1
In the Model Builder window, click Nafion®, EW 1100, Vapor Equilibrated, Protonated (mat3).
2
In the Settings window for Material, locate the Geometric Entity Selection section.
3
From the Selection list, choose Membranes.
Hydrogen Fuel Cell (fc)
Now start defining the physics. Begin when the fuel cell interface.
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.
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 Membranes.
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_in.
Water Absorption-Desorption, H2 Side 1
1
In the Model Builder window, click Water Absorption-Desorption, H2 Side 1.
2
In the Settings window for Water Absorption-Desorption, H2 Side, locate the Absorption–Desorption Condition section.
3
From the Electrolyte material list, choose Nafion®, EW 1100, Vapor Equilibrated, Protonated (mat3).
Water Absorption-Desorption, O2 Side 1
1
In the Model Builder window, click Water Absorption-Desorption, O2 Side 1.
2
In the Settings window for Water Absorption-Desorption, O2 Side, locate the Absorption–Desorption Condition section.
3
From the Electrolyte material list, choose Nafion®, EW 1100, Vapor Equilibrated, Protonated (mat3).
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 Hydrogen GDLs.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigmas_GDL.
5
Locate the Gas Transport section. In the εg text field, type epsg_GDL.
6
In the κg text field, type kappag_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 Oxygen GDLs.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigmas_GDL.
5
Locate the Gas Transport section. In the εg text field, type epsg_GDL.
6
In the κg text field, type kappag_GDL.
H2 Gas Diffusion Layer 2 (manifolds)
The gas channels and manifolds are modeled as homogenized domains, that is, the individual channels are not resolved in the geometry. Instead the permeability is used to control the flow velocity and the pressure drop.
1
In the Physics toolbar, click  Domains and choose H2 Gas Diffusion Layer.
2
In the Settings window for H2 Gas Diffusion Layer, type H2 Gas Diffusion Layer 2 (manifolds) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Hydrogen Gas Manifolds.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
5
Locate the Gas Transport section. In the εg text field, type eps_gas_BPP.
6
In the κg text field, type perm_gas_BPP.
O2 Gas Diffusion Layer 2 (manifolds)
1
In the Physics toolbar, click  Domains and choose O2 Gas Diffusion Layer.
2
In the Settings window for O2 Gas Diffusion Layer, type O2 Gas Diffusion Layer 2 (manifolds) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Oxygen Gas Manifolds.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
5
Locate the Gas Transport section. In the εg text field, type eps_gas_BPP.
6
In the κg text field, type perm_gas_BPP.
Current Collector 1 (end blocks)
1
In the Physics toolbar, click  Domains and choose Current Collector.
2
In the Settings window for Current Collector, type Current Collector 1 (end blocks) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Current Collector and Feeder Plates.
4
Locate the Electrode Charge Transport section. From the σs list, choose From material.
Current Collector (with cooling flow)
1
In the Physics toolbar, click  Domains and choose Current Collector.
2
In the Settings window for Current Collector, type Current Collector (with cooling flow) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Non-Gas Cooling Channels and Manifolds.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
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 Hydrogen Gas Diffusion Electrodes.
4
Locate the Electrode Thickness section. In the dgde text field, type L_CL.
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 Electrode Kinetics section.
3
In the i0,ref(T) text field, type i0_H2_ref.
4
Locate the Active Specific Surface Area section. In the av text field, type a_CL.
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 Oxygen Gas Diffusion Electrodes.
4
Locate the Electrode Thickness section. In the dgde text field, type L_CL.
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 Electrode Kinetics section.
3
In the i0,ref(T) text field, type i0_O2_ref.
4
In the αa text field, type alphaa_O2.
5
Locate the Active Specific Surface Area section. In the av text field, type a_CL.
H2 Gas Diffusion Layer 3 (x-directed channels)
As for the manifolds, the gas flow channels in the bipolar plates are also modeled as homogenized porous domains. By using anisotropic permeabilities, the flow direction can be controlled.
1
In the Physics toolbar, click  Domains and choose H2 Gas Diffusion Layer.
2
In the Settings window for H2 Gas Diffusion Layer, type H2 Gas Diffusion Layer 3 (x-directed channels) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Hydrogen Gas Channels, x-directed.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
5
Locate the Gas Transport section. In the εg text field, type eps_gas_BPP.
6
From the list, choose Diagonal.
7
Specify the κg matrix as
perm_gas_BPP
0
0
0
0
0
0
0
perm_gas_BPP
H2 Gas Diffusion Layer 4 (y-directed channels)
1
In the Physics toolbar, click  Domains and choose H2 Gas Diffusion Layer.
2
In the Settings window for H2 Gas Diffusion Layer, type H2 Gas Diffusion Layer 4 (y-directed channels) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Hydrogen Gas Channels, y-directed.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
5
Locate the Gas Transport section. In the εg text field, type eps_gas_BPP.
6
From the list, choose Diagonal.
7
Specify the κg matrix as
0
0
0
0
perm_gas_BPP
0
0
0
perm_gas_BPP
O2 Gas Diffusion Layer 3 (x-directed channels)
1
In the Physics toolbar, click  Domains and choose O2 Gas Diffusion Layer.
2
In the Settings window for O2 Gas Diffusion Layer, type O2 Gas Diffusion Layer 3 (x-directed channels) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Oxygen Gas Channels, x-directed.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
5
Locate the Gas Transport section. In the εg text field, type eps_gas_BPP.
6
From the list, choose Diagonal.
7
Specify the κg matrix as
perm_gas_BPP
0
0
0
0
0
0
0
perm_gas_BPP
O2 Gas Diffusion Layer 4 (y-directed channels)
1
In the Physics toolbar, click  Domains and choose O2 Gas Diffusion Layer.
2
In the Settings window for O2 Gas Diffusion Layer, type O2 Gas Diffusion Layer 4 (y-directed channels) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Oxygen Gas Channels, y-directed.
4
Locate the Electrode Charge Transport section. In the σs text field, type sigma_BPP_eff.
5
Locate the Gas Transport section. In the εg text field, type eps_gas_BPP.
6
From the list, choose Diagonal.
7
Specify the κg matrix as
0
0
0
0
perm_gas_BPP
0
0
0
perm_gas_BPP
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
In the Settings window for Electric Ground, locate the Boundary Selection section.
3
From the Selection list, choose Current Feeder Tab.
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
In the Settings window for Electric Potential, locate the Boundary Selection section.
3
From the Selection list, choose Current Collector Tab.
4
Locate the Electric Potential section. In the ϕs,bnd text field, type E_stack.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) > H2 Gas Phase 1 click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Composition section.
3
From the Mixture specification list, choose Humidified mixture.
4
In the Thum text field, type T_in.
5
In the pA,hum text field, type 1[atm]+p_in_an.
6
Locate the Initial Pressure section. In the p0 text field, type p_in_an.
H2 Gas Phase 1
In the Model Builder window, click H2 Gas Phase 1.
H2 Inlet 1
1
In the Physics toolbar, click  Attributes and choose H2 Inlet.
2
In the Settings window for H2 Inlet, locate the Boundary Selection section.
3
From the Selection list, choose Hydrogen Inlets.
4
Locate the Inlet Flow Type section. From the Inlet flow type list, choose Mixture composition constraint.
H2 Gas Phase 1
In the Model Builder window, click H2 Gas Phase 1.
H2 Outlet 1
1
In the Physics toolbar, click  Attributes and choose H2 Outlet.
2
In the Settings window for H2 Outlet, locate the Boundary Selection section.
3
From the Selection list, choose Hydrogen Outlets.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Hydrogen Fuel Cell (fc) > O2 Gas Phase 1 click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Composition section.
3
From the Mixture specification list, choose Humidified air.
4
In the Thum text field, type T_in.
5
In the pA,hum text field, type 1[atm]+p_in_cath.
6
Locate the Initial Pressure section. In the p0 text field, type p_in_cath.
O2 Gas Phase 1
In the Model Builder window, click O2 Gas Phase 1.
O2 Inlet 1
1
In the Physics toolbar, click  Attributes and choose O2 Inlet.
2
In the Settings window for O2 Inlet, locate the Boundary Selection section.
3
From the Selection list, choose Oxygen Inlets.
4
Locate the Inlet Flow Type section. From the Inlet flow type list, choose Mixture composition constraint.
O2 Gas Phase 1
In the Model Builder window, click O2 Gas Phase 1.
O2 Outlet 1
1
In the Physics toolbar, click  Attributes and choose O2 Outlet.
2
In the Settings window for O2 Outlet, locate the Boundary Selection section.
3
From the Selection list, choose Oxygen Outlets.
To improve convergence, enable stabilization to all gas flow domains. By default, stabilization is turned off in GDLs and GDEs.
4
Click the  Show More Options button in the Model Builder toolbar.
5
In the Show More Options dialog, in the tree, select the checkbox for the node Physics > Stabilization.
6
Click OK.
7
In the Model Builder window, click Hydrogen Fuel Cell (fc).
8
In the Settings window for Hydrogen Fuel Cell, click to expand the Consistent Stabilization section.
9
Find the O2 gas mixture subsection. Clear the Add stabilization only to gas flow channels checkbox.
10
Find the H2 gas mixture subsection. Clear the Add stabilization only to gas flow channels checkbox.
Hydrogen Fuel Cell (fc)
In the Model Builder window, collapse the Component 1 (comp1) > Hydrogen Fuel Cell (fc) node.
Darcy’s Law (dl)
Now set up the Darcy’s law interface, which defines the cooling flow velocity and pressure.
1
In the Model Builder window, under Component 1 (comp1) click Darcy’s Law (dl).
2
In the Settings window for Darcy’s Law, locate the Domain Selection section.
3
From the Selection list, choose Cooling Channels and Manifolds.
Use linear elements to match the element order of the other interfaces, and to save some memory and computational time.
4
Click to expand the Discretization section. From the Pressure list, choose Linear.
Porous Matrix 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) > Porous Medium 1 click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type eps_l_BPP.
4
From the κ list, choose User defined. In the associated text field, type perm_cool_BPP.
Inlet 1
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
In the Settings window for Inlet, locate the Boundary Selection section.
3
From the Selection list, choose Cooling Inlets.
4
Locate the Velocity section. In the U0 text field, type v_cool_in.
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
In the Settings window for Outlet, locate the Boundary Selection section.
3
From the Selection list, choose Cooling Outlets.
4
Locate the Boundary Condition section. From the Boundary condition list, choose Pressure.
5
In the Model Builder window, collapse the Darcy’s Law (dl) node.
Heat Transfer in Porous Media (ht)
In the Model Builder window, under Component 1 (comp1) click Heat Transfer in Porous Media (ht).
Solid 1 - Membranes
1
In the Physics toolbar, click  Domains and choose Solid.
2
In the Settings window for Solid, type Solid 1 - Membranes in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Membranes.
Since the model is stationary, the density and heat capacity need not to be defined.
4
Locate the Thermodynamics, Solid section. From the ρ list, choose User defined. From the Cp list, choose User defined.
Solid 2 - GDLs
1
In the Physics toolbar, click  Domains and choose Solid.
2
In the Settings window for Solid, type Solid 2 - GDLs in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose GDLs.
4
Locate the Heat Conduction, Solid section. From the k list, choose User defined. The thermal conductivity of the gas diffusion layers is anisotropic, featuring a higher conductivity in the in-plane direction.
5
From the list, choose Diagonal.
6
Specify the k matrix as
kappa_GDL_IP
0
0
0
kappa_GDL_IP
0
0
0
kappa_GDL_TP
7
Locate the Thermodynamics, Solid section. From the ρ list, choose User defined. From the Cp list, choose User defined.
Solid 3 - End Blocks
1
In the Physics toolbar, click  Domains and choose Solid.
2
In the Settings window for Solid, type Solid 3 - End Blocks in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Current Collector and Feeder Plates.
Solid 4 - Gas Manifolds
1
In the Physics toolbar, click  Domains and choose Solid.
2
In the Settings window for Solid, type Solid 4 - Gas Manifolds in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Gas Manifolds.
4
Locate the Heat Conduction, Solid section. From the k list, choose User defined. In the associated text field, type kappa_BPP_eff.
5
Locate the Thermodynamics, Solid section. From the ρ list, choose User defined. From the Cp list, choose User defined.
Porous Medium 1
Couple the heat transfer to the velocity in the cooling channels as follows:
Fluid 1
1
In the Model Builder window, click Fluid 1.
2
In the Settings window for Fluid, locate the Heat Convection section.
3
From the u list, choose Total Darcy velocity field (dl/porous1).
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type eps_gas_BPP.
4
Locate the Thermodynamics, Porous Matrix section. From the ρb list, choose User defined. Locate the Heat Conduction, Porous Matrix section. From the kb list, choose User defined. In the associated text field, type kappa_BPP_eff.
5
Locate the Thermodynamics, Porous Matrix section. From the Cp,b list, choose User defined.
Inflow 1
1
In the Physics toolbar, click  Boundaries and choose Inflow.
2
In the Settings window for Inflow, locate the Boundary Selection section.
3
From the Selection list, choose Cooling Inlets.
4
Locate the Upstream Properties section. In the Tustr text field, type T_in.
Outflow 1
1
In the Physics toolbar, click  Boundaries and choose Outflow.
2
In the Settings window for Outflow, locate the Boundary Selection section.
3
From the Selection list, choose Cooling Outlets.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T text field, type T_in.
4
In the Model Builder window, collapse the Heat Transfer in Porous Media (ht) node.
Multiphysics
Couple the heat sources in the fuel cell interface with heat transfer as follows. This will also couple the temperature defined by heat transfer to the fuel cell interface.
Electrochemical Heating 1 (ech1)
In the Physics toolbar, click  Multiphysics Couplings and choose Domain > Electrochemical Heating.
Mesh 1
Use a user-defined mesh for this model. The mesh sequence is based on sweeping the mesh in the through-plane direction of the stack. Since the model is finalized as an assembly in the geometry sequence, with assembly pairs located in the middle of the membranes, the mesh nodes do not need match along these boundaries.
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
From the Selection list, choose Membrane Boundaries.
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-3.
Swept 1
1
In the Mesh toolbar, click  Swept.
2
In the Settings window for Swept, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Membranes.
5
Click to expand the Source Faces section. From the Selection list, choose Hydrogen Gas Diffusion Electrodes.
6
Click to expand the Destination Faces section. From the Selection list, choose Oxygen Gas Diffusion Electrodes.
7
Locate the Mesh Generation section. From the Elements list, choose Prisms.
Distribution 1
1
Right-click Swept 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Distribution section.
3
In the Number of elements text field, type 2.
Size
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Extremely fine.
4
Click  Build All.
Definitions
To facilitate viewing when setting up the sequence, add a View with scaling in the z direction.
View 8
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose View.
Camera
1
In the Model Builder window, expand the View 8 node, then click Camera.
2
In the Settings window for Camera, locate the Camera section.
3
From the View scale list, choose Manual.
4
In the z scale text field, type 50.
5
Click  Update.
Hide for Physics 1
1
In the Model Builder window, right-click View 8 and choose Hide for Physics.
2
In the Settings window for Hide for Physics, locate the Geometric Entity Selection section.
3
From the Selection list, choose Current Collector and Feeder Plates.
Mesh 1
Click the  Zoom Extents button in the Graphics toolbar.
Swept 2
1
In the Mesh toolbar, click  Swept.
2
In the Settings window for Swept, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Oxygen GDLs.
Distribution 1
1
Right-click Swept 2 and choose Distribution.
2
Right-click Distribution 1 and choose Build Selected.
Swept 3
1
In the Mesh toolbar, click  Swept.
2
In the Settings window for Swept, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Hydrogen GDLs.
Distribution 1
1
Right-click Swept 3 and choose Distribution.
2
Right-click Distribution 1 and choose Build Selected.
Swept 4
1
In the Mesh toolbar, click  Swept.
2
In the Settings window for Swept, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Cooling Channels and Manifolds.
Distribution 1
1
Right-click Swept 4 and choose Distribution.
2
In the Settings window for Distribution, locate the Distribution section.
3
In the Number of elements text field, type 2.
4
Click  Build Selected.
Swept 5
In the Mesh toolbar, click  Swept.
Distribution 1
1
Right-click Swept 5 and choose Distribution.
2
In the Settings window for Distribution, locate the Distribution section.
3
In the Number of elements text field, type 2.
4
Click  Build Selected.
5
In the Graphics window toolbar, clicknext to  Go to Default View, then choose Go to View 1.
6
Click the  Zoom Extents button in the Graphics toolbar.
Study 1
Now start setting up the study sequence.
The first initialization step solves for a primary current distribution, that is, a simplified set of the fuel cell equations, not including activation or concentration overpotentials.
Add a second initialization step solves that solves for a secondary current distribution, which includes activation but not concentration overpotentials. The second step will use the results of the first step as initial values for the dependent variables.
Step 3: Current Distribution Initialization 2
1
In the Study toolbar, click  More Study Steps and choose Other > Current Distribution Initialization.
2
Right-click Step 3: Current Distribution Initialization 2 and choose Move Up.
3
In the Settings window for Current Distribution Initialization, locate the Study Settings section.
4
From the Current distribution type list, choose Secondary.
Step 3: Stationary
The third step we will use to only solve for the pressures in the fuel cell interface. Specifying that we will only solve for the pressures will be done later.
1
In the Model Builder window, click Step 3: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the Solve for column of the table, under Component 1 (comp1), select the checkbox for Hydrogen Fuel Cell (fc).
4
In the Solve for column of the table, under Component 1 (comp1), clear the checkboxes for Darcy’s Law (dl) and Heat Transfer in Porous Media (ht).
Step 4: Stationary 2
Add a fourth step to only solve for the convective cooling flow (Darcy’s Law).
1
In the Study toolbar, click  Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the Solve for column of the table, under Component 1 (comp1), clear the checkbox for Hydrogen Fuel Cell (fc).
4
In the Solve for column of the table, under Component 1 (comp1), select the checkbox for Darcy’s Law (dl).
5
In the Solve for column of the table, under Component 1 (comp1), clear the checkbox for Heat Transfer in Porous Media (ht).
Step 5: Stationary 3
Add a fifth step which we will use to compute the final solution, for a range of potentials. We will assume the convective cooling flow (Darcy’s Law) not to be affected by the temperature, and hence we need not to solve for it in this step.
1
In the Study toolbar, click  Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the Solve for column of the table, under Component 1 (comp1), select the checkbox for Hydrogen Fuel Cell (fc).
4
In the Solve for column of the table, under Component 1 (comp1), clear the checkbox for Darcy’s Law (dl).
5
In the Solve for column of the table, under Component 1 (comp1), select the checkbox for Heat Transfer in Porous Media (ht).
6
Click to expand the Study Extensions section. Select the Auxiliary sweep checkbox.
7
Click  Add.
8
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
E_cell_avg (Average cell voltage)
 
V
9
In the table, click to select the cell at row number 1 and column number 3.
10
Click  Range.
11
In the Range dialog, type E_cell_avg_start in the Start text field.
12
In the Step text field, type -0.1.
13
In the Stop text field, type 0.55.
14
Click Replace.
Solution 1 (sol1)
Set the third step to solve only for the fuel cell pressures as follows:
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Dependent Variables 3.
3
In the Settings window for Dependent Variables, locate the General section.
4
From the Defined by study step list, choose User defined.
5
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 3 node, then click Chemical Potential (comp1.fc.mu0).
6
In the Settings window for Field, locate the General section.
7
Clear the Solve for this field checkbox.
8
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 3 click Electrolyte Potential (comp1.fc.phil).
9
In the Settings window for Field, locate the General section.
10
Clear the Solve for this field checkbox.
11
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 3 click Electric Potential (comp1.fc.phis).
12
In the Settings window for Field, locate the General section.
13
Clear the Solve for this field checkbox.
14
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 3 click Mass Fraction (comp1.fc.wH2O_H2).
15
In the Settings window for Field, locate the General section.
16
Clear the Solve for this field checkbox.
17
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 3 click Mass Fraction (comp1.fc.wH2O_O2).
18
In the Settings window for Field, locate the General section.
19
Clear the Solve for this field checkbox.
20
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 3 click Mass Fraction (comp1.fc.wN2_O2).
21
In the Settings window for Field, locate the General section.
22
Clear the Solve for this field checkbox.
23
In the Study toolbar, click  Compute.
Results
Set default units for result presentation.
Preferred Units 1
1
In the Results toolbar, click  Configurations and choose Preferred Units.
2
In the Settings window for Preferred Units, locate the Units section.
3
Click  Add Physical Quantity.
4
In the Physical Quantity dialog, select General > Temperature (K) in the tree.
5
Click OK.
6
In the Settings window for Preferred Units, locate the Units section.
7
In the table, enter the following settings:
Quantity
Unit
Preferred unit
Temperature
K
°C
8
In the Model Builder window, click Preferred Units 1.
9
Click  Apply.
Electrode Potential with Respect to Ground (fc)
Reproduce the figures from the Results and Discussion section as follows:
1
In the Model Builder window, expand the Electrode Potential with Respect to Ground (fc) node.
Arrow Volume 1, Multislice 1
1
In the Model Builder window, under Results > Electrode Potential with Respect to Ground (fc), Ctrl-click to select Multislice 1 and Arrow Volume 1.
2
Right-click and choose Disable.
Electrode Potential with Respect to Ground (fc)
1
In the Model Builder window, click Electrode Potential with Respect to Ground (fc).
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges checkbox.
Surface 1
1
Right-click Electrode Potential with Respect to Ground (fc) 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.phis - Electric potential - V.
3
Click to expand the Coloring and Style section. From the Color table list, choose Cynanthus.
4
In the Electrode Potential with Respect to Ground (fc) toolbar, click  Plot.
Mole Fraction, H2, Surface (fc)
1
In the Model Builder window, under Results click Mole Fraction, H2, Surface (fc).
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges checkbox.
Surface 1
1
In the Model Builder window, expand the Mole Fraction, H2, Surface (fc) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Cynanthus.
4
In the Mole Fraction, H2, Surface (fc) toolbar, click  Plot.
Mole Fraction, O2, Surface (fc)
1
In the Model Builder window, under Results click Mole Fraction, O2, Surface (fc).
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges checkbox.
Surface 1
1
In the Model Builder window, expand the Mole Fraction, O2, Surface (fc) node, then click Surface 1.
2
In the Settings window for Surface, click to expand the Coloring and Style section.
3
From the Color table list, choose Cynanthus.
4
In the Mole Fraction, O2, Surface (fc) toolbar, click  Plot.
Multislice 1
1
In the Model Builder window, expand the Water Activity (Relative Humidity) (fc) node.
2
Right-click Multislice 1 and choose Disable.
Water Activity in Oxygen GDEs
1
In the Model Builder window, under Results click Water Activity (Relative Humidity) (fc).
2
In the Settings window for 3D Plot Group, type Water Activity in Oxygen GDEs in the Label text field.
3
Click to expand the Title section. From the Title type list, choose Label.
Surface 1
1
Right-click Water Activity in Oxygen GDEs 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.aw - Water activity (relative humidity) - 1.
3
Click to expand the Coloring and Style section. From the Color table list, choose Cynanthus.
Selection 1
1
Right-click Surface 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Oxygen Gas Diffusion Electrodes.
Water Activity in Oxygen GDEs
1
In the Model Builder window, under Results click Water Activity in Oxygen GDEs.
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges checkbox.
4
In the Water Activity in Oxygen GDEs toolbar, click  Plot.
Temperature in MEAs
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Temperature in MEAs in the Label text field.
3
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
Surface 1
1
Right-click Temperature in MEAs 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) > Heat Transfer in Porous Media > Temperature > T - Temperature - K.
3
Locate the Coloring and Style section. From the Color table list, choose HeatCameraLight.
Selection 1
1
Right-click Surface 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose GDLs.
Surface 2
1
In the Model Builder window, under Results > Temperature in MEAs right-click Surface 1 and choose Duplicate.
2
In the Settings window for Surface, click to expand the Title section.
3
From the Title type list, choose None.
4
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
Selection 1
1
In the Model Builder window, expand the Surface 2 node, then click Selection 1.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Membranes.
4
In the Temperature in MEAs toolbar, click  Plot.
Temperature in Cooling Channels
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Temperature in Cooling Channels in the Label text field.
3
Locate the Title section. From the Title type list, choose Label.
4
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
5
Locate the Color Legend section. Select the Show units checkbox.
Streamline Multislice 1
1
In the Temperature in Cooling Channels toolbar, click  More Plots and choose Streamline Multislice.
2
In the Settings window for Streamline Multislice, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Darcy’s Law > Velocity and pressure > dl.u,dl.v,dl.w - Total Darcy velocity field.
3
Locate the Multiplane Data section. Find the X-planes subsection. In the Planes text field, type 0.
4
Find the Y-planes subsection. In the Planes text field, type 0.
5
Find the Z-planes subsection. From the Entry method list, choose Coordinates.
6
In the Coordinates text field, type range(D_cc+D_bpp*0.5,D_cell,D_cc+D_bpp*0.5+N_cells*D_cell).
7
Locate the Streamline Positioning section. From the Positioning list, choose Uniform density.
8
Locate the Coloring and Style section. Find the Point style subsection. From the Type list, choose Arrow.
9
From the Arrow distribution list, choose Equal inverse time.
Color Expression 1
1
Right-click Streamline Multislice 1 and choose Color Expression.
2
In the Settings window for Color Expression, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Heat Transfer in Porous Media > Temperature > T - Temperature - K.
3
Locate the Coloring and Style section. From the Color table list, choose HeatCameraLight.
Cut Line 3D 1
1
In the Results toolbar, click  Cut Line 3D.
2
In the Settings window for Cut Line 3D, locate the Line Data section.
3
In row Point 1, set X to W_plate/2.
4
In row Point 1, set Y to 0.95*L_plate.
5
In row Point 1, set Z to D_cc.
6
In row Point 2, set X to W_plate/2.
7
In row Point 2, set Y to 0.95*L_plate.
8
In row Point 2, set Z to D_cc+D_cell*N_cells+D_bpp.
Mid-Stack Temperature Toward Cooling Outlets
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Mid-Stack Temperature Toward Cooling Outlets in the Label text field.
3
Click to expand the Title section. From the Title type list, choose Label.
Line Graph 1
Right-click Mid-Stack Temperature Toward Cooling Outlets and choose Line Graph.
Mid-Stack Temperature Toward Cooling Outlets
Locate the Data section. From the Dataset list, choose Cut Line 3D 1.
Line Graph 1
1
In the Model Builder window, click Line Graph 1.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Heat Transfer in Porous Media > Temperature > T - Temperature - K.
3
Locate the x-Axis Data section. From the Parameter list, choose Expression.
4
In the Expression text field, type z.
5
Click to expand the Legends section. Select the Show legends checkbox.
6
In the Mid-Stack Temperature Toward Cooling Outlets toolbar, click  Plot.